source: trunk/LMDZ.PLUTO/libf/phypluto/ch4surf.F @ 3959

      SUBROUTINE ch4surf(ngrid,nlayer,nq,ptimestep,capcal,tsurf,
     &   pdtsurf,pplev,pdpsurf,pq,pdq,pqsurf,pdqsurf,pdqch4,pdqsch4)
  
      use callkeys_mod, only: dayfrac, thresh_non2, no_ch4frost
      use comcstfi_mod, only: g, r
      use comgeomfi_h 
      use comsaison_h, only: fract
      use planete_mod, only: z0
      use tracer_h, only: igcm_ch4_gas,igcm_ch4_ice,igcm_n2,mmol,lw_ch4
      IMPLICIT NONE

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

#include "dimensions.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 capcal(ngrid)           ! surface heat capacity (J m-2 K-1)
      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(ngrid)
      REAL zpsrf(ngrid)
      REAL zq_ch4(ngrid)
      REAL zq_n2surf(ngrid)
      REAL ztsurf(ngrid)
      REAL gamm(ngrid) ! 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=0.3  ! 6
      v=0.4  ! 3 
      uv=sqrt(u**2+v**2)
       
      pdqsch4(:)=0.      
      pdqch4(:)=0.      

      !! Update CH4, pressure
      DO ig=1,ngrid
          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(ngrid,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,ngrid

          !! 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)

         !! option to not form any CH4 frost
         if (no_ch4frost.and.pqsurf(ig,igcm_ch4_ice).lt.5000.) then

           pdqch4(ig)=0.
           pdqsch4(ig)=0.

         else

           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 (lw_ch4*pdqsch4(ig)*ptimestep/capcal(ig).gt.1.) then
            pdqsch4(ig)=1./(lw_ch4*ptimestep)*capcal(ig)
           endif
           if (lw_ch4*pdqsch4(ig)*ptimestep/capcal(ig).lt.-1.) then
            pdqsch4(ig)=-1./(lw_ch4*ptimestep)*capcal(ig)
           endif
           !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) 

         endif

      ENDDO

      RETURN
      END