SUBROUTINE cosurf(ngrid,nlayer,nq,ptimestep, 
     &  tsurf,pplev,pdpsurf,pq,pdq,pqsurf,pdqsurf,pdqco,pdqsco)
       
      use callkeys_mod, only: dayfrac
      use comcstfi_mod, only: g, r
      use comgeomfi_h 
      use comsaison_h, only: fract
      use planete_mod, only: z0
      use tracer_h, only: igcm_co_gas,igcm_co_ice,igcm_ch4_ice,
     &                            igcm_n2,mmol
      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 :
!                             pdqco, pdqs_co  

!-----------------------------------------------------------------------
!     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 qsurf_n2(ngrid)

      ! Output
      REAL pdqco(ngrid)
      REAL pdqsco(ngrid)

      ! local
      REAL qsat(ngrid)
      REAL zpsrf(ngrid)
      REAL zq_co(ngrid)
      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=z0 !1.e-2 ! rugosity
      cdrag=(vonk/log(alt/z00))**2
      
      u=6.
      v=3.
      uv=sqrt(u**2+v**2)
       
      DO ig=1,ngrid
         pdqsco(ig)=0.      
         pdqco(ig)=0.      
      ENDDO

      DO ig=1,ngrid
          zpsrf(ig)=pplev(ig,1)
          zq_co(ig)=pq(ig,1,igcm_co_gas)
!    &    +  pdq(ig,1,igcm_co_gas)*ptimestep ! pas utile si fast 
      ENDDO

      call cosat(ngrid,tsurf,zpsrf,qsat,pqsurf(:,igcm_n2),
     &                                       pqsurf(:,igcm_ch4_ice))

      DO ig=1,ngrid
         rho = zpsrf(ig) / (r *  tsurf(ig) )
         
         pdqsco(ig)=(-rho*uv*cdrag*(qsat(ig)-zq_co(ig)))

         if ((-pdqsco(ig)*ptimestep).gt.(pqsurf(ig,igcm_co_ice))) then
                  !write(*,*)'cosurf : on sublime plus que qsurf_co!'
                  pdqsco(ig)=-pqsurf(ig,igcm_co_ice)/ptimestep
         endif
         pdqco(ig)=-pdqsco(ig)*g/zpsrf(ig) 
      ENDDO
      RETURN
      END