!$gpum horizontal klon
MODULE lmdz_reevap

CONTAINS

   SUBROUTINE reevap(klon, klev, temp, qv, ql, qi, &
      &         d_temp_eva, d_qv_eva, d_ql_eva, d_qi_eva)

!=============================================================================
! This routine evaporates (and sublimates) the condensed water at the
! beginning of the physics in order to work with the total water
! specific content of water which is a conserved variable
! (one of the so-called Bett's variable)
!
! reference: Madeleine et al. 2020, doi:10.1029/2020MS002046
! contact: C. Rio, catherine.rio@meteo.fr
!=============================================================================

      USE lmdz_reevap_ini, only: fl_cor_ebil, iflag_ice_thermo
      USE lmdz_reevap_ini, only: RLVTT, RLSTT, RCPD, RVTMP2

      IMPLICIT none

!======================================================================
! Declarations
!======================================================================

! Input variables:
!-----------------
      INTEGER, INTENT(IN) :: klon, klev  ! horizontal and vertical dimensions
      REAL, DIMENSION(klon, klev), INTENT(IN) :: temp ! temperature [K]
      REAL, DIMENSION(klon, klev), INTENT(IN) :: qv ! specific humidity (vapor) [kg/kg]
      REAL, DIMENSION(klon, klev), INTENT(IN) :: ql ! specific liquid water [kg/kg]
      REAL, DIMENSION(klon, klev), INTENT(IN) :: qi ! specific ice water [kg/kg]

! Output variables:
!-----------------
      REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_temp_eva ! temperature increment due to evap. [K]
      REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_qv_eva ! specific humidity increment due to evap [kg/kg]
      REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_ql_eva ! specific liquid water increment due to evap [kg/kg]
      REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_qi_eva ! specific ice water increment due to evap [kg/kg]

! Local variables
!-----------------
      INTEGER :: i, k
      REAL    :: za, zb, zlvdcp, zlsdcp, zdelta

!======================================================================
! Computation of evaporation and sublimation
!======================================================================

      DO k = 1, klev
         DO i = 1, klon
            IF (fl_cor_ebil .GT. 0) THEN
               ! the mass of condensates is taken into account when applying enthalpy conservation
               ! (first principle of thermodynamics)
               zlvdcp = RLVTT/RCPD/(1.0 + RVTMP2*(qv(i, k) + ql(i, k) + qi(i, k)))
               zlsdcp = RLSTT/RCPD/(1.0 + RVTMP2*(qv(i, k) + ql(i, k) + qi(i, k)))
            ELSE
               zlvdcp = RLVTT/RCPD/(1.0 + RVTMP2*qv(i, k))
               zlsdcp = RLSTT/RCPD/(1.0 + RVTMP2*qv(i, k))
            END IF

            IF (iflag_ice_thermo .EQ. 0) THEN
               ! if iflag_ice_thermo = 0, only the liquid phase is considered thermodynamically
               ! active in LMDZ, no only liquid water is evaporated
               zlsdcp=zlvdcp
               zdelta = 0.0
               zb = MAX(0.0, ql(i, k))
               za = - MAX(0.0 , ql(i,k)) &
                  * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
               d_temp_eva(i, k) = za
               d_qv_eva(i, k) = zb
               d_ql_eva(i, k) = -ql(i, k)
               d_qi_eva(i, k) = 0.

            ELSE
               ! both liquid and ice are evaporated / sublimated
               zb = MAX(0.0, ql(i, k) + qi(i, k))
               za = -MAX(0.0, ql(i, k))*zlvdcp &
                    - MAX(0.0, qi(i, k))*zlsdcp
               d_temp_eva(i, k) = za
               d_qv_eva(i, k) = zb
               d_ql_eva(i, k) = -ql(i, k)
               d_qi_eva(i, k) = -qi(i, k)
            END IF

         END DO
      END DO

      RETURN

   END SUBROUTINE reevap

END MODULE lmdz_reevap