SUBROUTINE ch4surf(ngrid,nlayer,nq,ptimestep,tsurf,pdtsurf,
     &           pplev,pdpsurf,pq,pdq,pqsurf,pdqsurf,pdqch4,pdqsch4)
  
      use planet_h         
      use comgeomfi_h 
      IMPLICIT NONE

c----------------
c   declarations:
c   -------------

#include "dimensions.h"
#include "dimphys.h"
#include "comcstfi.h"
#include "comsaison.h"
#include "callkeys.h"
#include "surfdat.h"
#include "tracer.h"
!#include "planete.h"

! Routine for nogcm : sublimation/condensation scheme at the surface
! Output : tendancy for methane mixing ratio and surface reservoir :
!                             pdqch4, pdqs_ch4  

!-----------------------------------------------------------------------
!     Arguments


      INTEGER ngrid,nlayer,nq
      REAL ptimestep  
      INTEGER ig,iq

      ! input :
      REAL tsurf(ngrid)
      REAL pplev(ngrid,nlayer+1)
      REAL pdpsurf(ngrid)
      REAL pq(ngrid,nlayer,nq)
      REAL pdq(ngrid,nlayer,nq)
      REAL pqsurf(ngrid,nq)
      REAL pdqsurf(ngrid,nq)
      REAL pdtsurf(ngrid)

      ! Output
      REAL pdqch4(ngrid)
      REAL pdqsch4(ngrid)

      ! local
      REAL qsat(ngridmx)
      REAL zpsrf(ngridmx)
      REAL zq_ch4(ngridmx)
      REAL zq_n2surf(ngridmx)
      REAL ztsurf(ngridmx)
      REAL gamm(ngridmx) ! activity coefficient
      REAL rho,u,v,uv,z00,cdrag,alt
      REAL vonk       ! Von Karman Constant 
      SAVE vonk        
      DATA vonk/0.4/

      ! Calculation of turbulent flux : F=rho*cdrag*uv*(qsat-zq)
   
      ! Calcul de cdrag
      alt=5.   ! m
      !z00=1.e-2 ! rugosity
      z00=z0
      cdrag=(vonk/log(alt/z00))**2
      
      u=6.  ! 6
      v=3.  ! 3 
      uv=sqrt(u**2+v**2)
       
      pdqsch4(:)=0.      
      pdqch4(:)=0.      

      !! Update CH4, pressure
      DO ig=1,ngridmx
          zpsrf(ig)=pplev(ig,1)
          zq_ch4(ig)=pq(ig,1,igcm_ch4_gas)
     &    +  pdq(ig,1,igcm_ch4_gas)*ptimestep 
          zq_n2surf(ig)=pqsurf(ig,igcm_n2)
     &    +  pdqsurf(ig,igcm_n2)*ptimestep 
          ztsurf(ig)=tsurf(ig) !+pdtsurf(ig)*ptimestep
      ENDDO

      !! Get qsat for CH4
      call methanesat(ngridmx,ztsurf,zpsrf,qsat,zq_n2surf(:))

      !! Dayfrac: Fraction of the daytime where we do not condense CH4 in N2
          ! corresponds to cold layer of N2 pushing CH4 or depleting CH4 near the surface
          ! By default, we do not take this into account : dayfrac=0. We condense all the time
          ! dayfrac = 1 : we do not condense CH4 in N2 during daytime
          ! dayfrac = 2 : all saturated, so we never condense CH4 in N2
          ! dayfrac = 0.5 : we do not condense CH4 during half of the daytime (we condense during all night + half of daytime)

      !! Loop
      DO ig=1,ngridmx

         !! Take into account activity coefficient
         gamm(ig)=499.9-21.8*ztsurf(ig)+0.249*ztsurf(ig)**2
     &             -1.3*(zq_ch4(ig)*
     &             mmol(igcm_n2)/mmol(igcm_ch4_gas)*100.-0.6)/0.3
         gamm(ig)=max(gamm(ig),1.)
         qsat(ig)=qsat(ig)*gamm(ig)

         rho = zpsrf(ig) / (r *  tsurf(ig) )
         !! Condensation Flux
         pdqsch4(ig)=(-rho*uv*cdrag*(qsat(ig)-zq_ch4(ig)))

         if (dayfrac.gt.0.and.zq_n2surf(ig).gt.thresh_non2) then
            if (dayfrac.gt.1.) then
               pdqsch4(ig)=min(pdqsch4(ig),0.)
            else
               if (pdqsch4(ig).gt.0.) then ! condensation
                 pdqsch4(ig)=pdqsch4(ig)*(1.-fract(ig)*dayfrac)
               endif
            endif
         endif
                      
         !! Conserve mass if reservoir depleted                   
         if ((-pdqsch4(ig)*ptimestep).gt.
     &                          (pqsurf(ig,igcm_ch4_ice))) then
            pdqsch4(ig)=-pqsurf(ig,igcm_ch4_ice)/ptimestep
         endif

         if (pdqsch4(ig)*ptimestep.gt.zq_ch4(ig)*zpsrf(ig)/g) then
            pdqsch4(ig)=zq_ch4(ig)/ptimestep*zpsrf(ig)/g
         endif

         !! Security to avoid large changes in temperatures due to
         !latent heat                   
         if (pdqsch4(ig)*ptimestep.gt.0.25) then
            pdqsch4(ig)=0.25/ptimestep
         endif
         if (pdqsch4(ig)*ptimestep.lt.-0.25) then
            pdqsch4(ig)=-0.25/ptimestep
         endif

         !! Atm tendency
         pdqch4(ig)=-pdqsch4(ig)*g/zpsrf(ig) 

      ENDDO

      RETURN
      END