module vdifc_pluto_mod

    implicit none

    contains

    SUBROUTINE vdifc_pluto(ngrid,nlay,nq,ppopsk,                   &
                     ptimestep,pcapcal,lecrit,               &
                     pplay,pplev,pzlay,pzlev,pz0,            &
                     pu,pv,ph,pq,pt,ptsrf,pemis,pqsurf,      &
                     pdufi,pdvfi,pdhfi,pdqfi,pdtfi,pfluxsrf, &
                     pdudif,pdvdif,pdhdif,pdtsrf,pq2,        &
                     pdqdif,pdqsdif,qsat_ch4,qsat_ch4_l1)

      use comgeomfi_h
      use callkeys_mod, only: carbox, methane, condcosurf, condensn2, condmetsurf,&
                              kmix_proffix, vertdiff, tracer, kmixmin, no_n2frost
      use datafile_mod, only: datadir
      use surfdat_h, only: phisfi
      use comcstfi_mod, only: g, r, rcp, cpp
      USE tracer_h, only: igcm_ch4_gas, igcm_ch4_ice, igcm_co_gas, igcm_co_ice,&
                          igcm_n2, lw_ch4, lw_co, lw_n2

      implicit none

!=======================================================================
!
!   subject:
!   --------
!   Turbulent diffusion (mixing) for potential T, U, V and tracer
!
!   Shema implicite
!   On commence par rajouter au variables x la tendance physique
!   et on resoult en fait:
!      x(t+1) =  x(t) + dt * (dx/dt)phys(t)  +  dt * (dx/dt)difv(t+1)
!
!   author:
!   ------
!      Hourdin/Forget/Fournier
!=======================================================================

!-----------------------------------------------------------------------
!   declarations:
!   -------------

#include "dimensions.h"

!
!   arguments:
!   ----------

      INTEGER ngrid,nlay
      REAL ptimestep
      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
      REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
      REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
      REAL ptsrf(ngrid),pemis(ngrid)
      REAL pdufi(ngrid,nlay),pdvfi(ngrid,nlay),pdhfi(ngrid,nlay)
      REAL pdtfi(ngrid,nlay)
      REAL pt(ngrid,nlay)
      REAL pfluxsrf(ngrid)
      REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay),pdhdif(ngrid,nlay)
      REAL pdtsrf(ngrid),pcapcal(ngrid)
      REAL pq2(ngrid,nlay+1)
      REAL qsat_ch4(ngrid)
      REAL qsat_co(ngrid)
      REAL qsat_ch4_l1(ngrid)
      REAL zq1temp_ch4(ngrid)

!    Argument added for condensation:
!      REAL n2ice (ngrid)
      REAL ppopsk(ngrid,nlay)
      logical lecrit
      REAL pz0

!    Traceurs :
      integer nq
      REAL pqsurf(ngrid,nq)
      real pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
      real pdqdif(ngrid,nlay,nq)
      real pdqdifeddy(ngrid,nlay,nq)
      real pdqsdif(ngrid,nq),pdqsdif1(ngrid,nq)

!   local:
!   ------

      INTEGER ilev,ig,ilay,nlev

      REAL z4st,zdplanck(ngrid)
      REAL zkv(ngrid,nlay+1),zkh(ngrid,nlay+1)
      REAL zcdv(ngrid),zcdh(ngrid),sat2(ngrid)
      REAL zcdv_true(ngrid),zcdh_true(ngrid)
      REAL zu(ngrid,nlay),zv(ngrid,nlay)
      REAL zh(ngrid,nlay),zt(ngrid,nlay)
      REAL ztsrf2(ngrid),sat(ngrid),sat1(ngrid)
      REAL z1(ngrid),z2(ngrid)
      REAL za(ngrid,nlay),zb(ngrid,nlay)
      REAL zb0(ngrid,nlay)
      REAL zc(ngrid,nlay),zd(ngrid,nlay)
      REAL zcst1
      REAL zu2
      EXTERNAL SSUM,SCOPY
      REAL SSUM
      LOGICAL firstcall
      SAVE firstcall
      real,dimension(:),save,allocatable :: qsat_co_factor ! factor to prevent co frost formation if no n2 frost
!$OMP THREADPRIVATE(qsat_co_factor)

      !!read fixed profile for kmix
      integer Nfine,ifine
      parameter(Nfine=701)
      character(len=100) :: file_path
      real,save :: levdat(Nfine),kmixdat(Nfine)
      real :: kmix_z(nlay)   ! kmix from kmix_proffix

!     variable added for N2 condensation:
!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      REAL hh , zhcond(ngrid,nlay)
!      REAL latcond,tcond1p4Pa
      REAL tcond1p4Pa
      REAL acond,bcond
      SAVE acond,bcond
!      DATA latcond,tcond1p4Pa/2.5e5,38/
      DATA tcond1p4Pa/38/

!    Tracers :
!    ~~~~~~~
      INTEGER iq
      REAL zq(ngrid,nlay,nq)
      REAL zq1temp_co(ngrid)
      REAL rho(ngrid) ! near surface air density
      DATA firstcall/.true./

!    ** un petit test de coherence
!       --------------------------

      IF (firstcall) THEN
         IF(ngrid.NE.ngrid) THEN
            write(*,*) 'STOP dans vdifc'
            write(*,*) 'probleme de dimensions :'
            write(*,*) 'ngrid  =',ngrid
            write(*,*) 'ngrid  =',ngrid
            STOP
         ENDIF
!        To compute: Tcond= 1./(bcond-acond*log(.7143*p)) (p in pascal)
         write(*,*) 'In vdifc: Tcond(P=1.4Pa)=',tcond1p4Pa,' Lcond=',lw_n2
         bcond=1./tcond1p4Pa
         acond=r/lw_n2
         write(*,*) '          acond,bcond',acond,bcond

        firstcall=.false.

        ! If fixed profile of kmix
        IF (kmix_proffix) then
            file_path=trim(datadir)//'/gas_prop/kmix.txt'
            open(114,file=file_path,form='formatted')
            do ifine=1,Nfine
              read(114,*) levdat(ifine), kmixdat(ifine)
            enddo
            close(114)
        ENDIF

        ! If fixed distribution of N2, then no CO frost either
        ALLOCATE(qsat_co_factor(ngrid))
        qsat_co_factor(:)=1.
        IF (no_n2frost) then
           DO ig=1,ngrid
              if (pqsurf(ig,igcm_n2).eq.0.) then
                   qsat_co_factor(ig) = 1.e6
              endif
           ENDDO
        ENDIF

      ENDIF

!-----------------------------------------------------------------------
!    1. initialisation
!    -----------------

      nlev=nlay+1

!    ** calcul de rho*dz et dt*rho/dz=dt*rho**2 g/dp
!       avec rho=p/RT=p/ (R Theta) (p/ps)**kappa
!       ----------------------------------------

      DO ilay=1,nlay
         DO ig=1,ngrid
            za(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/g
         ENDDO
      ENDDO

      zcst1=4.*g*ptimestep/(r*r)
      DO ilev=2,nlev-1
         DO ig=1,ngrid
            zb0(ig,ilev)=pplev(ig,ilev)* &
           (pplev(ig,1)/pplev(ig,ilev))**rcp / &
           (ph(ig,ilev-1)+ph(ig,ilev))
            zb0(ig,ilev)=zcst1*zb0(ig,ilev)*zb0(ig,ilev)/ &
           (pplay(ig,ilev-1)-pplay(ig,ilev))
!            write(300,*)'zb0',zb0(ig,ilev)
         ENDDO
      ENDDO
      DO ig=1,ngrid
	  	 zb0(ig,1)=ptimestep*pplev(ig,1)/(r*ptsrf(ig))
      ENDDO

!    ** diagnostique pour l'initialisation
!       ----------------------------------

!      IF(lecrit) THEN
!         ig=ngrid/2+1
!         write(*,*) 'Pression (mbar) ,altitude (km),u,v,theta, rho dz'
!         DO ilay=1,nlay
!            WRITE(*,'(6f11.5)') &
!           .01*pplay(ig,ilay),.001*pzlay(ig,ilay), &
!           pu(ig,ilay),pv(ig,ilay),ph(ig,ilay),za(ig,ilay)
!         ENDDO
!         write(*,*) 'Pression (mbar) ,altitude (km),zb'
!         DO ilev=1,nlay
!            WRITE(*,'(3f15.7)') &
!           .01*pplev(ig,ilev),.001*pzlev(ig,ilev), &
!           zb0(ig,ilev)
!         ENDDO
!      ENDIF

!     Potential Condensation temperature:
!     -----------------------------------

        if (condensn2) then
          DO ilev=1,nlay
            DO ig=1,ngrid
              zhcond(ig,ilev) = &
       (1./(bcond-acond*log(.7143*pplay(ig,ilev))))/ppopsk(ig,ilev)
            END DO
          END DO
        else
          DO ilev=1,nlay
            DO ig=1,ngrid
              zhcond(ig,ilev) = 0.
            END DO
          END DO
        end if


!-----------------------------------------------------------------------
!   2. ajout des tendances physiques
!      -----------------------------

      DO ilev=1,nlay
         DO ig=1,ngrid
            zt(ig,ilev)=pt(ig,ilev)+pdtfi(ig,ilev)*ptimestep
            zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
            zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
            zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
            zh(ig,ilev)=max(zh(ig,ilev),zhcond(ig,ilev))
         ENDDO
      ENDDO
      if(tracer) then
        DO iq =1, nq
         DO ilev=1,nlay
           DO ig=1,ngrid
              zq(ig,ilev,iq)=pq(ig,ilev,iq)+pdqfi(ig,ilev,iq)*ptimestep
           ENDDO
         ENDDO
        ENDDO
      end if

!-----------------------------------------------------------------------
!   3. schema de turbulence
!      --------------------

!    ** source d'energie cinetique turbulente a la surface
!       (condition aux limites du schema de diffusion turbulente
!       dans la couche limite
!       ---------------------
      CALL vdif_cd( ngrid,nlay,pz0,g,pzlay,pu,pv,ptsrf,ph &
                  ,zcdv_true,zcdh_true)
      DO ig=1,ngrid
        zu2=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1)
        !TB16: GCM wind for flat hemisphere
        IF (phisfi(ig).eq.0.) zu2=max(2.,zu2)

        zcdv(ig)=zcdv_true(ig)*sqrt(zu2)
        zcdh(ig)=zcdh_true(ig)*sqrt(zu2)
      ENDDO

!    ** schema de diffusion turbulente dans la couche limite
!       ----------------------------------------------------

        CALL vdif_kc(ngrid,nlay,ptimestep,g,pzlev,pzlay &
                   ,pu,pv,ph,zcdv_true &
                   ,pq2,zkv,zkh)


!    Adding eddy mixing to mimic 3D general circulation in 1D
!    RW FF 2010
      if ((ngrid.eq.1)) then
        !kmixmin is the minimum eddy mix coeff in 1D

        ! If fixed profile of kmix
        IF (kmix_proffix) then
            !! Interpolate on the model vertical grid
            CALL interp_line(levdat,kmixdat,Nfine,&
                       pzlay(1,:)/1000.,kmix_z(:),nlay)

            do ilev=1,nlay
               zkh(1,ilev) = max(kmix_z(ilev),zkh(1,ilev))
               zkv(1,ilev) = max(kmix_z(ilev),zkv(1,ilev))
               !zkh(1,ilev) = kmixmin
               !zkv(1,ilev) = kmixmin
            end do
        ELSE
            do ilev=1,nlay
               zkh(1,ilev) = max(kmixmin,zkh(1,ilev))
               zkv(1,ilev) = max(kmixmin,zkv(1,ilev))
               !zkh(1,ilev) = kmixmin
               !zkv(1,ilev) = kmixmin
            end do
        ENDIF
      endif ! ngrid.eq.1

!!    Temporary:
!      zkh = zkh*0.1
!      zkv = zkv*0.1

!    ** diagnostique pour le schema de turbulence
!       -----------------------------------------

!      IF(lecrit) THEN
!         write(*,*) !         write(*,*) 'Diagnostic for the vertical turbulent mixing'
!         write(*,*) 'Cd for momentum and potential temperature'

!         write(*,*) zcdv(ngrid/2+1),zcdh(ngrid/2+1)
!         write(*,*) 'Mixing coefficient for momentum and pot.temp.'
!         DO ilev=1,nlay
!            write(*,*) zkv(ngrid/2+1,ilev),zkh(ngrid/2+1,ilev)
!         ENDDO
!      ENDIF

!-----------------------------------------------------------------------
!   4. inversion pour l'implicite sur u
!      --------------------------------

!    ** l'equation est
!       u(t+1) =  u(t) + dt * {(du/dt)phys}(t)  +  dt * {(du/dt)difv}(t+1)
!       avec
!       /zu/ = u(t) + dt * {(du/dt)phys}(t)   (voir paragraphe 2.)
!       et
!       dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1)
!       donc les entrees sont /zcdv/ pour la condition a la limite sol
!       et /zkv/ = Ku

      CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))
      CALL multipl(ngrid,zcdv,zb0,zb)

      DO ig=1,ngrid
         z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
         zc(ig,nlay)=za(ig,nlay)*zu(ig,nlay)*z1(ig)
         zd(ig,nlay)=zb(ig,nlay)*z1(ig)
      ENDDO

      DO ilay=nlay-1,1,-1
         DO ig=1,ngrid
            z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+   &
              zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
            zc(ig,ilay)=(za(ig,ilay)*zu(ig,ilay)+ &
              zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
            zd(ig,ilay)=zb(ig,ilay)*z1(ig)
         ENDDO
      ENDDO

      DO ig=1,ngrid
         zu(ig,1)=zc(ig,1)
      ENDDO
      DO ilay=2,nlay
         DO ig=1,ngrid
            zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)
         ENDDO
      ENDDO

!-----------------------------------------------------------------------
!   5. inversion pour l'implicite sur v
!      --------------------------------

!    ** l'equation est
!       v(t+1) =  v(t) + dt * {(dv/dt)phys}(t)  +  dt * {(dv/dt)difv}(t+1)
!       avec
!       /zv/ = v(t) + dt * {(dv/dt)phys}(t)   (voir paragraphe 2.)
!       et
!       dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1)
!       donc les entrees sont /zcdv/ pour la condition a la limite sol
!       et /zkv/ = Kv

      DO ig=1,ngrid
         z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
         zc(ig,nlay)=za(ig,nlay)*zv(ig,nlay)*z1(ig)
         zd(ig,nlay)=zb(ig,nlay)*z1(ig)
      ENDDO

      DO ilay=nlay-1,1,-1
         DO ig=1,ngrid
            z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ &
              zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
            zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+   &
              zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
            zd(ig,ilay)=zb(ig,ilay)*z1(ig)
         ENDDO
      ENDDO

      DO ig=1,ngrid
         zv(ig,1)=zc(ig,1)
      ENDDO
      DO ilay=2,nlay
         DO ig=1,ngrid
            zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)
         ENDDO
      ENDDO

!-----------------------------------------------------------------------
!   6. inversion pour l'implicite sur h sans oublier le couplage
!      avec le sol (conduction)
!      ------------------------

!    ** l'equation est
!       h(t+1) =  h(t) + dt * {(dh/dt)phys}(t)  +  dt * {(dh/dt)difv}(t+1)
!       avec
!       /zh/ = h(t) + dt * {(dh/dt)phys}(t)   (voir paragraphe 2.)
!       et
!       dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1)
!       donc les entrees sont /zcdh/ pour la condition de raccord au sol
!       et /zkh/ = Kh
!       -------------

      CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
      CALL multipl(ngrid,zcdh,zb0,zb)

      DO ig=1,ngrid
         z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
         zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay)*z1(ig)
         zd(ig,nlay)=zb(ig,nlay)*z1(ig)
      ENDDO

      DO ilay=nlay-1,1,-1
         DO ig=1,ngrid
            z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ &
              zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
            zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+   &
              zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
            zd(ig,ilay)=zb(ig,ilay)*z1(ig)
         ENDDO
      ENDDO

!    ** calcul de (d Planck / dT) a la temperature d'interface
!       ------------------------------------------------------

      z4st=4.*5.67e-8*ptimestep
      DO ig=1,ngrid
         zdplanck(ig)=z4st*pemis(ig)*ptsrf(ig)*ptsrf(ig)*ptsrf(ig)
      ENDDO

!    ** calcul de la temperature_d'interface et de sa tendance.
!       on ecrit que la somme des flux est nulle a l'interface
!       a t + \delta t,
!       !'est a dire le flux radiatif a {t + \delta t}
!       + le flux turbulent a {t + \delta t} &
!            qui 'ecrit K (T1-Tsurf) avec T1 = d1 Tsurf + c1
!            (notation K dt = /cpp*b/)
!       + le flux dans le sol a t
!       + l'evolution du flux dans le sol lorsque la temperature d'interface
!       passe de sa valeur a t a sa valeur a {t + \delta t}.
!       ----------------------------------------------------

      DO ig=1,ngrid
         z1(ig)=pcapcal(ig)*ptsrf(ig)+cpp*zb(ig,1)*zc(ig,1) &
          +zdplanck(ig)*ptsrf(ig)+ pfluxsrf(ig)*ptimestep
         z2(ig)= pcapcal(ig)+cpp*zb(ig,1)*(1.-zd(ig,1))+zdplanck(ig)
         ztsrf2(ig)=z1(ig)/z2(ig)
         pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep

!        Modif speciale N2 condensation:
!        tconds = 1./(bcond-acond*log(.7143*pplev(ig,1)))
!        if ((condensn2).and. &
!          ((n2ice(ig).ne.0).or.(ztsrf2(ig).lt.tconds)))then
!           zh(ig,1)=zc(ig,1)+zd(ig,1)*tconds
!        else
            zh(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig)
!        end if
      ENDDO

!    ** et a partir de la temperature au sol on remonte
!       -----------------------------------------------

      DO ilay=2,nlay
         DO ig=1,ngrid
            hh = max( zh(ig,ilay-1) , zhcond(ig,ilay-1) ) ! modif n2cond
            zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*hh
         ENDDO
      ENDDO


!-----------------------------------------------------------------------
!   TRACERS
!   -------

      if(tracer) then

!     Using the wind modified by friction for lifting and  sublimation
!     ----------------------------------------------------------------
!     This is computed above

!        DO ig=1,ngrid
!          zu2=zu(ig,1)*zu(ig,1)+zv(ig,1)*zv(ig,1)
!          zcdv(ig)=zcdv_true(ig)*sqrt(zu2)
!          zcdh(ig)=zcdh_true(ig)*sqrt(zu2)
!        ENDDO

!       Calcul flux vertical au bas de la premiere couche (cf dust on Mars)
!       -----------------------------------------------------------
        do ig=1,ngrid
          rho(ig) = zb0(ig,1) /ptimestep
!          zb(ig,1) = 0.
        end do

        pdqsdif(:,:) = 0
        pdqdif(:,:,:)=0.


!     TB: Eddy lifting of tracers :
! ****************************************************************
!        DO ig=1,ngrid
!!          option : injection only on an equatorial band
!!          if (abs(lati(ig))*180./pi.le.25.) then
!            pdqsdif(ig,igcm_eddy1e6) =-1.e-12
!           pdqsdif(ig,igcm_eddy1e7) =-1.e-12
!            pdqsdif(ig,igcm_eddy5e7) =-1.e-12
!            pdqsdif(ig,igcm_eddy1e8) =-1.e-12
!            pdqsdif(ig,igcm_eddy5e8) =-1.e-12
!          endif
!        ENDDO

!        pdqdifeddy(:,:,:)=0.
!       injection a 50 km
!        DO ig=1,ngrid
!            pdqdifeddy(ig,17,igcm_eddy1e6)=1e-12
!            pdqdifeddy(ig,17,igcm_eddy1e7)=1e-12
!            pdqdifeddy(ig,17,igcm_eddy5e7)=1e-12
!            pdqdifeddy(ig,17,igcm_eddy1e8)=1e-12
!            pdqdifeddy(ig,17,igcm_eddy5e8)=1e-12
!        ENDDO

!      Inversion pour l'implicite sur q
!       --------------------------------
        do iq=1,nq
          CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))

          if ((methane).and.(iq.eq.igcm_ch4_gas)) then
!            This line is required to account for turbulent transport
!            from surface (e.g. ice) to mid-layer of atmosphere:
             CALL multipl(ngrid,zcdv,zb0,zb(1,1))
          else if ((carbox).and.(iq.eq.igcm_co_gas)) then
             CALL multipl(ngrid,zcdv,zb0,zb(1,1))
          else ! (re)-initialize zb(:,1)
             zb(1:ngrid,1)=0
          end if

          DO ig=1,ngrid
               z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
               zc(ig,nlay)=za(ig,nlay)*zq(ig,nlay,iq)*z1(ig)
               zd(ig,nlay)=zb(ig,nlay)*z1(ig)
          ENDDO

          DO ilay=nlay-1,2,-1
               DO ig=1,ngrid
                z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ &
                zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
                zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+    &
                zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
                zd(ig,ilay)=zb(ig,ilay)*z1(ig)
               ENDDO
          ENDDO

            ! special case for methane and CO ice tracer: do not include
            ! ice tracer from surface (which is set when handling
            ! vapour case (see further down).
          if (methane.and.(iq.eq.igcm_ch4_ice)) then
            DO ig=1,ngrid
                z1(ig)=1./(za(ig,1)+zb(ig,1)+   &
                 zb(ig,2)*(1.-zd(ig,2)))
                zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+ &
                 zb(ig,2)*zc(ig,2)) *z1(ig)
            ENDDO
          else if (carbox.and.(iq.eq.igcm_co_ice)) then
            DO ig=1,ngrid
                z1(ig)=1./(za(ig,1)+zb(ig,1)+   &
                 zb(ig,2)*(1.-zd(ig,2)))
                zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+ &
                 zb(ig,2)*zc(ig,2)) *z1(ig)
            ENDDO

          else ! general case
            DO ig=1,ngrid
               z1(ig)=1./(za(ig,1)+zb(ig,1)+    &
              zb(ig,2)*(1.-zd(ig,2)))
               zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+  &
              zb(ig,2)*zc(ig,2) +   &
             (-pdqsdif(ig,iq)) *ptimestep) *z1(ig)  !tracer flux from surface
            ENDDO
          endif ! of if (methane.and.(iq.eq.igcm_ch4_ice))

!           Calculation for turbulent exchange with the surface (for ice)
          IF (methane.and.(iq.eq.igcm_ch4_gas)) then

            !! calcul de la valeur de q a la surface  :
            call methanesat(ngrid,ptsrf,pplev(1,1),   &
                                        qsat_ch4(:),pqsurf(:,igcm_n2))

            !! For output:
            call methanesat(ngrid,zt(:,1),pplev(1,1), &
                                      qsat_ch4_l1(:),pqsurf(:,igcm_n2))

            !! Prevent CH4 condensation at the surface
            if (.not.condmetsurf) then
               qsat_ch4=qsat_ch4*1.e6
            endif

            DO ig=1,ngrid
              zd(ig,1)=zb(ig,1)*z1(ig)
              zq1temp_ch4(ig)=zc(ig,1)+ zd(ig,1)*qsat_ch4(ig)
              pdqsdif(ig,igcm_ch4_ice)=rho(ig)*zcdv(ig) &
                            *(zq1temp_ch4(ig)-qsat_ch4(ig))
            END DO

            DO ig=1,ngrid
                if ((-pdqsdif(ig,igcm_ch4_ice)*ptimestep) &
                  .gt.(pqsurf(ig,igcm_ch4_ice))) then

                  pdqsdif(ig,igcm_ch4_ice)= &
                              -pqsurf(ig,igcm_ch4_ice)/ptimestep

                  z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2)))

                  zc(ig,1)=(za(ig,1)*zq(ig,1,igcm_ch4_gas)+ &
                 zb(ig,2)*zc(ig,2) +    &
                 (-pdqsdif(ig,igcm_ch4_ice)) *ptimestep) *z1(ig)

                  zq1temp_ch4(ig)=zc(ig,1)
                endif

             zq(ig,1,igcm_ch4_gas)=zq1temp_ch4(ig)

!            TB: MODIF speciale  pour  CH4
             pdtsrf(ig)=pdtsrf(ig)+ &
             lw_ch4*pdqsdif(ig,igcm_ch4_ice)/pcapcal(ig)


            ENDDO ! of DO ig=1,ngrid

          ELSE IF (carbox.and.(iq.eq.igcm_co_gas)) then

            !! Calculating saturation mixing ratio at surface 
            call cosat(ngrid,ptsrf,pplev(1,1),qsat_co, &
                   pqsurf(:,igcm_n2))

            !! Prevent CO condensation at the surface
            if (.not.condcosurf) then
               qsat_co(:)=qsat_co(:)*1.e6
            endif
            if (no_n2frost) then
               qsat_co(:)=qsat_co(:)*qsat_co_factor(:)     
            endif

            DO ig=1,ngrid
              zd(ig,1)=zb(ig,1)*z1(ig)
              zq1temp_co(ig)=zc(ig,1)+ zd(ig,1)*qsat_co(ig)
              pdqsdif(ig,igcm_co_ice)=rho(ig)*zcdv(ig) &
                            *(zq1temp_co(ig)-qsat_co(ig))
            END DO


            DO ig=1,ngrid
!   -------------------------------------------------------------
!           TEMPORAIRE : pour initialiser CO si glacier N2
!               meme si il n'y a pas de CO disponible
!             if (pqsurf(ig,igcm_n2).le.10.) then
!   -------------------------------------------------------------
!
                if ((-pdqsdif(ig,igcm_co_ice)*ptimestep) &
                  .gt.(pqsurf(ig,igcm_co_ice))) then
                  pdqsdif(ig,igcm_co_ice)= &
                              -pqsurf(ig,igcm_co_ice)/ptimestep
                  z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2)))
                  zc(ig,1)=(za(ig,1)*zq(ig,1,igcm_co_gas)+  &
                 zb(ig,2)*zc(ig,2) +    &
                 (-pdqsdif(ig,igcm_co_ice)) *ptimestep) *z1(ig)
                  zq1temp_co(ig)=zc(ig,1)
                endif
!             endif

               zq(ig,1,igcm_co_gas)=zq1temp_co(ig)

!            MODIF speciale  pour  CO / corrected by FF 2016
             pdtsrf(ig)=pdtsrf(ig)+ &
             lw_co*pdqsdif(ig,igcm_co_ice)/pcapcal(ig)

            ENDDO ! of DO ig=1,ngrid

          ELSE ! if (methane)

            DO ig=1,ngrid
               zq(ig,1,iq)=zc(ig,1)
            ENDDO

          END IF ! of IF ((methane).and.(iq.eq.igcm_ch4_gas))

          !! Diffusion verticale : shut down vertical transport if vertdiff = false
          if (vertdiff) then
           DO ilay=2,nlay
             DO ig=1,ngrid
                zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
            ENDDO
           ENDDO
          endif

        enddo ! of do iq=1,nq
      end if ! of if(tracer)

!-----------------------------------------------------------------------
!   8. calcul final des tendances de la diffusion verticale
!      ----------------------------------------------------
      DO ilev = 1, nlay
         DO ig=1,ngrid
            pdudif(ig,ilev)=(    zu(ig,ilev)-   &
           (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep)    )/ptimestep
            pdvdif(ig,ilev)=(    zv(ig,ilev)-   &
           (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep)    )/ptimestep
           hh = max(ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep ,  &
                zhcond(ig,ilev))        ! modif n2cond
            pdhdif(ig,ilev)=( zh(ig,ilev)- hh )/ptimestep
         ENDDO
      ENDDO

      if (tracer) then
        DO iq = 1, nq
          DO ilev = 1, nlay
            DO ig=1,ngrid
              pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)-   &
           (pq(ig,ilev,iq) + pdqfi(ig,ilev,iq)*ptimestep))/ptimestep
!             pdqdif(ig,ilev,iq)=pdqdifeddy(ig,ilev,iq)+(zq(ig,ilev,iq)-    &
!           (pq(ig,ilev,iq) + pdqfi(ig,ilev,iq)*ptimestep))/ptimestep
            ENDDO
         ENDDO
        ENDDO
      end if

!    ** diagnostique final
!       ------------------

      IF(lecrit) THEN
         write(*,*) 'In vdif'
         write(*,*) 'Ts (t) and Ts (t+st)'
         WRITE(*,'(a10,3a15)') &
        'theta(t)','theta(t+dt)','u(t)','u(t+dt)'
         write(*,*) ptsrf(ngrid/2+1),ztsrf2(ngrid/2+1)
         DO ilev=1,nlay
            WRITE(*,'(4f15.7)') &
           ph(ngrid/2+1,ilev),zh(ngrid/2+1,ilev), &
           pu(ngrid/2+1,ilev),zu(ngrid/2+1,ilev)

         ENDDO
      ENDIF

      RETURN
    !   END

    end subroutine vdifc_pluto

end module vdifc_pluto_mod