subroutine evap(dtime,pt, pq, pdq, pdt,  
     $     dqevap,dtevap, qevap, tevap) 
        
      use watercommon_h 

      implicit none

#include "dimensions.h"
#include "dimphys.h"
#include "tracer.h"
#include "comcstfi.h"

!==================================================================
!     
!     Purpose
!     -------
!     Evaporation of all water in the atmopshere.
!     
!     Authors
!     -------
!     Adapted from the LMDTERRE code by B. Charnay (2010)
!     Original author Z. X. Li (1993)
!     
!==================================================================

! Arguments:
      REAL pt(ngridmx,nlayermx) 
      REAL pq(ngridmx,nlayermx,nqmx) 
      REAL pdt(ngridmx,nlayermx) 
      REAL pdq(ngridmx,nlayermx,nqmx) 
      REAL dqevap(ngridmx,nlayermx)
      REAL dtevap(ngridmx,nlayermx)
      REAL qevap(ngridmx,nlayermx,nqmx)
      REAL dtime
 
! Local:
      REAL tevap(ngridmx,nlayermx) 
      REAL zlvdcp
      REAL zlsdcp
      REAL zdelta
      INTEGER l,ig

!
! Re-evaporer l'eau liquide nuageuse
!

      DO l=1,nlayermx
        DO ig=1,ngridmx
         qevap(ig,l,igcm_h2o_vap)=pq(ig,l,igcm_h2o_vap)	
     s          	          +pdq(ig,l,igcm_h2o_vap)*dtime
         qevap(ig,l,igcm_h2o_ice)=pq(ig,l,igcm_h2o_ice)	
     s                            +pdq(ig,l,igcm_h2o_ice)*dtime
         tevap(ig,l)=pt(ig,l)+pdt(ig,l)*dtime
	ENDDO
      ENDDO

      DO l = 1, nlayermx  
      	DO ig = 1, ngridmx
         zlvdcp=RLVTT/RCPD!/(1.0+RVTMP2*qevap(ig,l,igcm_h2o_vap))
         zlsdcp=RLSTT/RCPD!/(1.0+RVTMP2*qevap(ig,l,igcm_h2o_vap))
         ! ignoring qevap term creates huge difference when qevap large!!!

         zdelta = MAX(0.,SIGN(1.,T_h2O_ice_liq-tevap(ig,l))) ! what is this?
                                                  ! for division between water / liquid
	 dqevap(ig,l) = MAX(0.0,qevap(ig,l,igcm_h2o_ice))/dtime
         dtevap(ig,l) = - dqevap(ig,l)*RLVTT/RCPD ! exactly as in largescale.F
!         dtevap(ig,l) = - dqevap(ig,l)	
!     s                       * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
         qevap(ig,l,igcm_h2o_vap) = qevap(ig,l,igcm_h2o_vap)	
     s     	                    +dqevap(ig,l)*dtime
	 qevap(ig,l,igcm_h2o_ice) = 0.0
         tevap(ig,l) = tevap(ig,l)+dtevap(ig,l)*dtime

      	ENDDO
      ENDDO

      RETURN
      END