!$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 &
#ifdef ISO
               ,xtv,xtl,xti,d_xtv_eva,d_xtl_eva,d_xti_eva &
#endif
                    )

!=============================================================================
! 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
#ifdef ISO
      USE infotrac_phy, ONLY: ntiso
#ifdef ISOVERIF
      USE isotopes_verif_mod
!, ONLY: errmax,errmaxrel, iso_verif_o18_aberrant_nostop,deltaD,deltaO
      USE isotopes_mod, ONLY: iso_eau,iso_hdo,iso_o18,ridicule
#ifdef ISOTRAC
      USE isotrac_routines_mod, ONLY: iso_verif_traceur_pbidouille   
#endif
#endif
#endif

      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]
#ifdef ISO
      REAL, DIMENSION(ntiso,klon,klev), INTENT(IN) :: xtv ! isotopic ratios
      REAL, DIMENSION(ntiso,klon,klev), INTENT(IN) :: xtl
      REAL, DIMENSION(ntiso,klon,klev), INTENT(IN) :: xti
#endif

! 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]
#ifdef ISO
      REAL, DIMENSION(ntiso,klon,klev), INTENT(OUT) :: d_xtv_eva ! isotopic increments
      REAL, DIMENSION(ntiso,klon,klev), INTENT(OUT) :: d_xtl_eva
      REAL, DIMENSION(ntiso,klon,klev), INTENT(OUT) :: d_xti_eva
#endif

! Local variables
!-----------------
      INTEGER :: i, k, ixt
      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.

#ifdef ISO
               DO ixt=1,ntiso
                  zb = MAX(0.0,xtl(ixt,i,k))
                  d_xtv_eva(ixt,i,k) = zb
                  d_xtl_eva(ixt,i,k) = -xtl(ixt,i,k)
                  d_xti_eva(ixt,i,k) = 0.
               ENDDO
#ifdef ISOVERIF
               DO ixt=1,ntiso
                  CALL iso_verif_noNaN(xtv(ixt,i,k), &
                 &     'reevap 112: apres evap tot')
               ENDDO
               IF (iso_eau.GT.0) THEN
                  CALL iso_verif_egalite_choix( &
                 &     xtv(iso_eau,i,k),qv(i,k), &
                 &     'reevap 117+, après reevap totale',errmax,errmaxrel)
                  CALL iso_verif_egalite_choix( &
                 &     xtl(iso_eau,i,k),ql(i,k), &
                 &     'reevap 120+, après reevap totale',errmax,errmaxrel)
               ENDIF !if (iso_eau.gt.0) then       
               IF ((iso_HDO.GT.0).AND.(iso_O18.GT.0)) THEN
                  IF (qv(i,k).GT.ridicule) THEN
                     IF (iso_verif_o18_aberrant_nostop( &
                 &       xtv(iso_HDO,i,k)/qv(i,k), &
                 &       xtv(iso_O18,i,k)/qv(i,k), &
                 &       'reevap 127: apres reevap totale').eq.1) then
                        WRITE(*,*) 'i,k,qv(i,k)=',i,k,qv(i,k)
                        WRITE(*,*) 'd_qv_eva(i,k)=',d_qv_eva(i,k)
                        WRITE(*,*) 'deltaD(d_qv_eva(i,k))=',deltaD(d_xtv_eva(iso_HDO,i,k)/d_qv_eva(i,k))
                        WRITE(*,*) 'deltaO18(d_qv_eva(i,k))=',deltaO(d_xtv_eva(iso_O18,i,k)/d_qv_eva(i,k))
                        STOP
                     ENDIF !  if (iso_verif_o18_aberrant_nostop 
                  ENDIF !if (qv(i,k).gt.ridicule) then   
               ENDIF!if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then 
#ifdef ISOTRAC      
               CALL iso_verif_traceur(xtv(1,i,k), &
             &      'reevap 138a')
               CALL iso_verif_traceur_pbidouille(xtv(1,i,k), &
             &      'reevap 138b')
#endif
#endif
#endif

            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)
#ifdef ISO
               DO ixt=1,ntiso
                  zb = MAX(0.0,xtl(ixt,i,k)+xti(ixt,i,k))
                  d_xtv_eva(ixt,i,k) = zb
                  d_xtl_eva(ixt,i,k) = -xtl(ixt,i,k)
                  d_xti_eva(ixt,i,k) = -xti(ixt,i,k)
               ENDDO

#ifdef ISOVERIF
               DO ixt=1,ntiso
                  CALL iso_verif_noNaN(xtv(ixt,i,k), &
                 &     'reevap 165: apres evap tot')
               ENDDO
               IF (iso_eau.GT.0) THEN
                  CALL iso_verif_egalite_choix( &
           &           xtv(iso_eau,i,k),qv(i,k), &
           &          'reevap 170, après réévap totale',errmax,errmaxrel)
                  CALL iso_verif_egalite_choix( &
           &           xtl(iso_eau,i,k),ql(i,k), &
           &          'reevap 173, après réévap totale',errmax,errmaxrel)
                  CALL iso_verif_egalite_choix( &
           &           xti(iso_eau,i,k),qi(i,k), &
           &          'reevap 173b, après réévap totale',errmax,errmaxrel)
               ENDIF !if (iso_eau.gt.0) then
     
               IF ((iso_HDO.GT.0).AND.(iso_O18.GT.0)) THEN  
                  IF (q_seri(i,k).GT.ridicule) THEN  
                     IF (iso_verif_o18_aberrant_nostop( &
                 &       xtv(iso_HDO,i,k)/qv(i,k), &
                 &       xtv(iso_O18,i,k)/qv(i,k), &
                 &      'reevap 184: apres reevap totale').EQ.1) THEN
                        WRITE(*,*) 'i,k,qv(i,k)=',i,k,qv(i,k)                       
                        STOP
                     ENDIF !  if (iso_verif_o18_aberrant_nostop 
                  ENDIF !if (qv(i,k).gt.errmax) then   
               ENDIF !if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then        
#ifdef ISOTRAC      
               CALL iso_verif_traceur(xtv(1,i,k), &
         &          'reevap 192c')
               CALL iso_verif_traceur_pbidouille(xtv(1,i,k), &
         &          'reevap 192d')
#endif                  
#endif                  
#endif

            END IF

         END DO
      END DO

      RETURN

   END SUBROUTINE reevap

END MODULE lmdz_reevap