MODULE clmain_mod
      
      IMPLICIT NONE
      
      integer,save :: iflag_pbl ! PBL scheme number
C$OMP THREADPRIVATE(iflag_pbl)
      
      CONTAINS
c
c
      SUBROUTINE clmain(dtime,itap,
     .                  t,u,v,
     .                  rmu0, 
     .                  ts,
     .                  ftsoil,
     .                  paprs,pplay,ppk,radsol,albe,
     .                  solsw, sollw, sollwdown, fder,
     .                  rlon, rlat, dx, dy, 
     .                  q2,
     .                  debut, lafin, 
     .                  d_t,d_u,d_v,d_ts,
     .                  flux_t,flux_u,flux_v,cdragh,cdragm,
     .                  dflux_t,
     .                  zcoefh,zu1,zv1) 
cAA REM:
cAA-----
cAA Tracers are not handled here but in phytrac.
cAA REM bis :
cAA----------
cAA To extract exchange coefficients and wind in the first layer 
cAA 3 dedicated fields have been added: zcoefh,zu1 et zv1.
      use dimphy, only: klon, klev
      use mod_grid_phy_lmdz, only: nbp_lev
      use cpdet_phy_mod, only: t2tpot
      use turb_mod, only :yustar
      use soil_mod, only: nsoilmx
      use clesphys_mod, only: soil_model, ok_kzmin
      use YOMCST_mod, only: RD, RG
      use print_control_mod, only: lunout, prt_level
      
#ifdef CPP_XIOS      
      use xios_output_mod, only: send_xios_field
#endif      
      
      IMPLICIT none
c======================================================================
c Interface for the "boundary layer" (vertical diffusion)
c======================================================================

c
c Input arguments:
      REAL, INTENT(IN) :: dtime ! physics time step (s)
      INTEGER, INTENT(IN) :: itap ! physics time step counter
      REAL, INTENT(IN) :: t(klon,klev) ! atmospheric temperature (K)
      REAL, INTENT(IN) :: u(klon,klev) ! zonal wind (m/s)
      REAL, INTENT(IN) :: v(klon,klev) ! meridional wind (m/s)
      REAL, INTENT(IN) :: paprs(klon,klev+1) ! pressure at layer boundaries (Pa)
      REAL, INTENT(IN) :: pplay(klon,klev) ! pressure at mid-layer (Pa)
! ADAPTATION GCM FOR CP(T)
      REAL, INTENT(IN) :: ppk(klon,klev)
      REAL, INTENT(IN) :: rlon(klon) ! longitudes (deg)
      REAL, INTENT(IN) :: rlat(klon) ! latitudes (deg)
      REAL, INTENT(IN) :: dx(klon) ! mesh size along x (m)
      REAL, INTENT(IN) :: dy(klon) ! mesh size along y (m)
      REAL, INTENT(IN) :: rmu0(klon)         ! cosine of solar zenith angle
      LOGICAL, INTENT(IN) :: debut ! .true. if first call to physics
      LOGICAL, INTENT(IN) :: lafin ! .true. if last call to physics
      REAL, INTENT(IN) :: ts(klon) ! surface temperature (K)
      REAL, INTENT(IN) :: sollw(klon), solsw(klon), sollwdown(klon)
      
c IN/OUT arguments:
      REAL, INTENT(INOUT) :: fder(klon)
      REAL, INTENT(INOUT) :: flux_u(klon,klev) ! zonal wind stress (kg m/s)/(m**2 s) or Pa
      REAL, INTENT(INOUT) :: flux_v(klon,klev) ! meridional wind stress (kg m/s)/(m**2 s) or Pa
      REAL, INTENT(INOUT) :: radsol(klon) ! net radiative flux (positive towards ground) W/m**2
      REAL, INTENT(INOUT) :: q2(klon,klev+1) ! $q^2$ TKE at inter-layers
      
c output arguments
      REAL, INTENT(OUT) :: d_t(klon, klev) ! temperature increment (K)
      REAL, INTENT(OUT) :: d_u(klon, klev) ! zonal wind increment (m/s)
      REAL, INTENT(OUT) :: d_v(klon, klev) ! meridional wind increment (m/s)
      REAL, INTENT(OUT) :: flux_t(klon,klev) ! latent heat flux (CpT) J/m**2/s (W/m**2)
                                             ! (positive when downwards)
      REAL, INTENT(OUT) :: dflux_t(klon) ! derivative of sensible heat flux
      REAL, INTENT(OUT) :: cdragh(klon)
      REAL, INTENT(OUT) :: cdragm(klon)
      REAL, INTENT(OUT) :: d_ts(klon) ! surface temperature increment (K)
      REAL, INTENT(INOUT) :: albe(klon) ! albedo of the surface
cAA
      REAL, INTENT(OUT) :: zcoefh(klon,klev)
      REAL, INTENT(OUT) :: zu1(klon) ! zonal wind in 1st layer (m/s)
      REAL, INTENT(OUT) :: zv1(klon) ! meridional wind in 1st layer (m/s)
cAA
      REAL, INTENT(INOUT) :: ftsoil(klon,nsoilmx) ! subsurface temperatures (K)

c======================================================================
c Local variables
      REAL ytsoil(klon,nsoilmx)
      REAL yts(klon)
      REAL yalb(klon)
      REAL yu1(klon), yv1(klon)
      real ysollw(klon), ysolsw(klon), ysollwdown(klon)
      real yfder(klon), ytaux(klon), ytauy(klon)
      REAL yrads(klon)
C
      REAL y_d_ts(klon)
      REAL y_d_t(klon, klev)
      REAL y_d_u(klon, klev), y_d_v(klon, klev)
      REAL y_flux_t(klon,klev)
      REAL y_flux_u(klon,klev), y_flux_v(klon,klev)
      REAL y_dflux_t(klon)
      REAL ycoefh(klon,klev), ycoefm(klon,klev)
      REAL yu(klon,klev), yv(klon,klev)
      REAL yt(klon,klev)
      REAL ypaprs(klon,klev+1), ypplay(klon,klev), ydelp(klon,klev)
c
      REAL ycoefm0(klon,klev), ycoefh0(klon,klev)

      real yzlay(klon,klev),yzlev(klon,klev+1)
      real yteta(klon,klev)
      real ykmm(klon,klev+1),ykmn(klon,klev+1)
      real ykmq(klon,klev+1)
      real y_cd_m(klon),y_cd_h(klon)
c
      REAL u1lay(klon), v1lay(klon)
      REAL delp(klon,klev)
      INTEGER i, k
      INTEGER ni(klon), knon, j     
c======================================================================
      REAL zx_alf1, zx_alf2 ! ambient values used for extrapolation
c======================================================================
c
      LOGICAL,PARAMETER :: zxli=.FALSE. ! use simple functions testcase
c
      REAL zt, zdelta, zcor
C
      character (len = 20) :: modname = 'clmain'
      LOGICAL,PARAMETER :: check=.false.
C
      if (check) THEN
          write(*,*) trim(modname),'  klon=',klon
      endif
          
      DO k = 1, klev   ! thickness of atmospheric layers
      DO i = 1, klon
         delp(i,k) = paprs(i,k)-paprs(i,k+1)
      ENDDO
      ENDDO
      DO i = 1, klon  !  wind in the first layer
ccc         zx_alf1 = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
         zx_alf1 = 1.0
         zx_alf2 = 1.0 - zx_alf1
         u1lay(i) = u(i,1)*zx_alf1 + u(i,2)*zx_alf2
         v1lay(i) = v(i,1)*zx_alf1 + v(i,2)*zx_alf2
      ENDDO
c
c initialisation:
c
      DO i = 1, klon
         cdragh(i) = 0.0
         cdragm(i) = 0.0
         dflux_t(i) = 0.0
         zu1(i) = 0.0
         zv1(i) = 0.0
      ENDDO
      yts = 0.0
      yalb = 0.0
      yfder = 0.0
      ytaux = 0.0
      ytauy = 0.0
      ysolsw = 0.0
      ysollw = 0.0
      ysollwdown = 0.0
      yu1 = 0.0
      yv1 = 0.0
      yrads = 0.0
      ypaprs = 0.0
      ypplay = 0.0
      ydelp = 0.0
      yu = 0.0
      yv = 0.0
      yt = 0.0
      y_flux_u = 0.0
      y_flux_v = 0.0
      ycoefh=0.0
      ycoefm=0.0
C$$ PB
      y_dflux_t = 0.0
      ytsoil = 999999.
      DO i = 1, klon
         d_ts(i) = 0.0
      ENDDO
      flux_t = 0.
      flux_u = 0.
      flux_v = 0.
      DO k = 1, klev
      DO i = 1, klon
         d_t(i,k) = 0.0
         d_u(i,k) = 0.0
         d_v(i,k) = 0.0
         zcoefh(i,k) = 0.0
      ENDDO
      ENDDO
c
c identify indexes:
      DO j = 1, klon
         ni(j) = j
      ENDDO
      knon = klon

      DO j = 1, knon
      i = ni(j)
        yts(j) = ts(i)
        yalb(j) = albe(i)
        yfder(j) = fder(i)
        ytaux(j) = flux_u(i,1)
        ytauy(j) = flux_v(i,1)
        ysolsw(j) = solsw(i)
        ysollw(j) = sollw(i)
        ysollwdown(j) = sollwdown(i)
        yu1(j) = u1lay(i)
        yv1(j) = v1lay(i)
        yrads(j) =  ysolsw(j)+ ysollw(j)
        ypaprs(j,klev+1) = paprs(i,klev+1)
      END DO
C
c$$$ for the soil
      DO k = 1, nsoilmx
        DO j = 1, knon
          i = ni(j)
          ytsoil(j,k) = ftsoil(i,k)
        END DO  
      END DO 
      DO k = 1, klev
      DO j = 1, knon
      i = ni(j)
        ypaprs(j,k) = paprs(i,k)
        ypplay(j,k) = pplay(i,k)
        ydelp(j,k) = delp(i,k)
        yu(j,k) = u(i,k)
        yv(j,k) = v(i,k)
        yt(j,k) = t(i,k)
      ENDDO
      ENDDO
c
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c compute Cdrag and exchange coefficients
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c-------------------------------------------------
c  Old LMD computations. 
c  in routines coefkz, coefkz2, coefkzmin
c-------------------------------------------------

        CALL coefkz(knon, ypaprs, ypplay, ppk,
     .            yts, yu, yv, yt,
     .            ycoefm, ycoefh)

c       CALL coefkz2(knon, ypaprs, ypplay,yt,
c    .                  ycoefm0, ycoefh0)
c       DO k = 1, klev
c       DO i = 1, knon
c        ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
c        ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
c       ENDDO
c       ENDDO

c
cIM: 261103
        if (ok_kzmin) THEN
! ADAPTATION GCM FOR CP(T)
           print*," coefkzmin: NOT ADAPTATED..."
cIM cf FH: 201103 BEG

c   Compute minimal diffusion for very stable conditions.
c        call coefkzmin(knon,ypaprs,ypplay,yu,yv,yt,ycoefm
c    .   ,ycoefm0,ycoefh0)
c      call dump2d(nbp_lon,nbp_lat-2,ycoefm(2:klon-1,2), 'KZ         ')
c      call dump2d(nbp_lon,nbp_lat-2,ycoefm0(2:klon-1,2),'KZMIN      ')
 
         ycoefm0 = 1.e-3
         ycoefh0 = 1.e-3

         if ( 1.eq.1 ) then
          DO k = 1, klev
          DO i = 1, knon
           ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
           ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
          ENDDO
          ENDDO
         endif
cIM cf FH: 201103 END
        endif !ok_kzmin
cIM: 261103

      IF (iflag_pbl.ge.3) then
c-------------------------------------------------
c MELLOR ET YAMADA adapted to Mars Richard Fournier and Frederic Hourdin
c-------------------------------------------------

         yzlay(1:knon,1)=
     .   RD*yt(1:knon,1)/(0.5*(ypaprs(1:knon,1)+ypplay(1:knon,1)))
     .   *(ypaprs(1:knon,1)-ypplay(1:knon,1))/RG
         do k=2,klev
            yzlay(1:knon,k)=
     .      yzlay(1:knon,k-1)+RD*0.5*(yt(1:knon,k-1)+yt(1:knon,k))
     .      /ypaprs(1:knon,k)*(ypplay(1:knon,k-1)-ypplay(1:knon,k))/RG
         enddo

! ADAPTATION GCM FOR CP(T)
         call t2tpot(knon*nbp_lev,yt,yteta,ppk)

         yzlev(1:knon,1)=0.
         yzlev(1:knon,klev+1)=2.*yzlay(1:knon,klev)-yzlay(1:knon,klev-1)
         do k=2,klev
            yzlev(1:knon,k)=0.5*(yzlay(1:knon,k)+yzlay(1:knon,k-1))
         enddo


         y_cd_m(1:knon) = ycoefm(1:knon,1)
         y_cd_h(1:knon) = ycoefh(1:knon,1)

         call ustarhb(knon,yu,yv,y_cd_m, yustar)

        if (prt_level > 9) THEN
          WRITE(lunout,*)'USTAR = ',yustar
        ENDIF

c   iflag_pbl can be used as mixing length

         if (iflag_pbl.ge.11) then
            call vdif_kcay(knon,dtime,rg,rd,ypaprs,yt
     s      ,yzlev,yzlay,yu,yv,yteta
     s      ,y_cd_m,ykmm,ykmn,yustar,
     s      iflag_pbl)
         else
            call yamada4(knon,dtime,rg,rd,ypaprs,yt
     s      ,yzlev,yzlay,yu,yv,yteta
     s      ,y_cd_m,q2,ykmm,ykmn,ykmq,yustar,
     s      iflag_pbl)
         endif

         ycoefm(1:knon,1)=y_cd_m(1:knon)
         ycoefh(1:knon,1)=y_cd_h(1:knon)
         ycoefm(1:knon,2:klev)=ykmm(1:knon,2:klev)
         ycoefh(1:knon,2:klev)=ykmn(1:knon,2:klev)

c-------------------------------------------------
      ENDIF

c        print*,"Kzm=",ycoefm(100,:)
c        print*,"Kzh=",ycoefh(100,:)

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c compute diffusion for winds "u" et "v"
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yu,ypaprs,ypplay,ydelp,
     s            y_d_u,y_flux_u)
      CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yv,ypaprs,ypplay,ydelp,
     s            y_d_v,y_flux_v)
 
c for the coupling
      ytaux = y_flux_u(:,1)
      ytauy = y_flux_v(:,1)


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c compute diffusion for "q" and "h"
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! ADAPTATION GCM FOR CP(T)
      CALL clqh(dtime, itap, debut,lafin,
     e          rlon, rlat, dx, dy,
     e          knon, 
     e          soil_model, ytsoil,
     e          rmu0, 
     e          yu1, yv1, ycoefh,
     e          yt,yts,ypaprs,ypplay,ppk,
     e          ydelp,yrads,yalb, 
     e          yfder, ytaux, ytauy,
     e          ysollw, ysollwdown, ysolsw,
     s          y_d_t, y_d_ts,
     s          y_flux_t, y_dflux_t) 

      DO j = 1, knon
         i = ni(j)
         d_ts(i) = y_d_ts(j)
c----------------------------------------
c VENUS TEST: tendancy on Tsurf forced = 0
c        d_ts(i) = 0.
c----------------------------------------
         albe(i) = yalb(j)
         cdragh(i) = cdragh(i) + ycoefh(j,1)
         cdragm(i) = cdragm(i) + ycoefm(j,1)
         dflux_t(i) = dflux_t(i) + y_dflux_t(j)
         zu1(i) = zu1(i) + yu1(j)
         zv1(i) = zv1(i) + yv1(j)
      END DO

c for the subsurface
      DO k = 1, nsoilmx
        DO j = 1, knon
         i = ni(j)
         ftsoil(i, k) = ytsoil(j,k)
        ENDDO
      END DO
      
      DO k = 1, klev
        DO j = 1, knon
         i = ni(j)
         flux_t(i,k) = y_flux_t(j,k)
         flux_u(i,k) = y_flux_u(j,k)
         flux_v(i,k) = y_flux_v(j,k)
         d_t(i,k) = d_t(i,k) + y_d_t(j,k)
         d_u(i,k) = d_u(i,k) + y_d_u(j,k)
         d_v(i,k) = d_v(i,k) + y_d_v(j,k)
         zcoefh(i,k) = zcoefh(i,k) + ycoefh(j,k)
        ENDDO
      ENDDO

      END SUBROUTINE clmain

C=================================================================
C=================================================================
C=================================================================
C=================================================================

      SUBROUTINE clqh(dtime,itime, debut,lafin,
     e                rlon, rlat, dx, dy, 
     e                knon, 
     $                soil_model,tsoil,
     e                rmu0, 
     e                u1lay,v1lay,coef,
     e                t,ts,paprs,pplay,ppk,
     e                delp,radsol,albedo, 
     e                fder, taux, tauy,
     $                sollw, sollwdown, swnet, 
     s                d_t, d_ts, 
     s                flux_t, dflux_s)

      use interface_surf, only: interfsurf_hq
      use dimphy, only: klon, klev
      use soil_mod, only: nsoilmx
      use mod_grid_phy_lmdz, only: nbp_lon, nbp_lat, nbp_lev
      use cpdet_phy_mod, only: t2tpot,tpot2t,cpdet
      use YOMCST_mod, only: RD, RG, RKAPPA, RMD

      IMPLICIT none
c======================================================================
c Compute vertical diffusion for "h"
c======================================================================

c Arguments:
c Parametres d'entree
      INTEGER, INTENT(IN) :: knon
      REAL, INTENT(IN) :: dtime              ! intervalle du temps (s)
      REAL, INTENT(IN) :: u1lay(klon)        ! vitesse u de la 1ere couche (m/s)
      REAL, INTENT(IN) :: v1lay(klon)        ! vitesse v de la 1ere couche (m/s)
      REAL, INTENT(IN) :: coef(klon,klev)    ! le coefficient d'echange (m**2/s)
c                               multiplie par le cisaillement du 
c                               vent (dV/dz); la premiere valeur
c                               indique la valeur de Cdrag (sans unite)
      REAL, INTENT(IN) :: t(klon,klev)       ! temperature (K)
      REAL, INTENT(IN) :: ts(klon)           ! temperature du sol (K)
      REAL, INTENT(IN) :: paprs(klon,klev+1) ! pression a inter-couche (Pa)
      REAL, INTENT(IN) :: pplay(klon,klev)   ! pression au milieu de couche (Pa)
! ADAPTATION GCM POUR CP(T)
      REAL, INTENT(IN) :: ppk(klon,klev)     ! fonction d'Exner milieu de couche
      REAL, INTENT(IN) :: delp(klon,klev)    ! epaisseur de couche en pression (Pa)
      REAL, INTENT(INOUT) :: radsol(klon)       ! ray. net au sol (Solaire+IR) W/m2
      REAL, INTENT(IN) :: rmu0(klon)         ! cosinus de l'angle solaire zenithal
      REAL, INTENT(IN) :: rlon(klon), rlat(klon), dx(klon), dy(klon)
      INTEGER, INTENT(IN) :: itime
      LOGICAL, INTENT(IN) :: debut, lafin
      REAL, INTENT(INOUT) :: fder(klon)
      REAL, INTENT(IN) :: taux(klon), tauy(klon)
      REAL, INTENT(IN) :: sollw(klon), sollwdown(klon)
      REAL, INTENT(IN) :: swnet(klon)
      
      
c$$$C PB ajout pour soil
      LOGICAL, INTENT(IN) :: soil_model
      REAL, INTENT(IN) :: tsoil(klon, nsoilmx)
c
c Parametre de sorties
      REAL, INTENT(OUT) :: albedo(klon)       ! albedo de la surface
      REAL, INTENT(OUT) :: d_t(klon,klev)     ! incrementation de "t"
      REAL, INTENT(OUT) :: d_ts(klon)         ! incrementation de "ts"
      REAL, INTENT(OUT) :: flux_t(klon,klev)  ! (diagnostic) flux de la chaleur
c                               sensible, flux de Cp*T, positif vers
c                               le bas: j/(m**2 s) c.a.d.: W/m2
      REAL, INTENT(OUT) :: dflux_s(klon) ! derivee du flux sensible dF/dTs
c======================================================================
c Variables locales
      INTEGER i, k
      REAL zx_ch(klon,klev)
      REAL zx_dh(klon,klev)
      REAL zx_buf2(klon)
      REAL zx_coef(klon,klev)
      REAL local_h(klon,klev) ! enthalpie potentielle
      REAL local_ts(klon)
      REAL swdown(klon)
c======================================================================
c contre-gradient pour la chaleur sensible: Kelvin/metre
! ADAPTATION GCM POUR CP(T)
      REAL gamt(klon,klev),zt(klon,klev)
      REAL z_gamah(klon,klev)
      REAL zdelz
c======================================================================
c Rajout pour l'interface
      real zlev1(klon)
      real temp_air(klon)
      real epot_air(klon)
      real tq_cdrag(klon), petAcoef(klon)
      real petBcoef(klon)
      real p1lay(klon), pkh1(klon)

! Parametres de sortie
      real fluxsens(klon)
      real tsurf_new(klon), alb_new(klon)

      character (len = 20) :: modname = 'Debut clqh'
      LOGICAL,PARAMETER :: check=.false.

C
c----------------------
c ATMOSPHERE PROFONDE
      real p_lim,p_ref
      real rsmu_mod(klon,klev),zrsmu_mod(klon,klev)
      real ratio_mod(klon,klev) 
! theta_mod = theta * ratio_mod => ratio_mod = mmm0/mmm(p)

      p_lim = 6e6
      p_ref = 8.9e6

      DO k = 1, klev
      DO i = 1, klon
        ratio_mod(i,k) = 1.
        rsmu_mod(i,k) = RD
c ATM PROFONDE DESACTIVEE !
       if (    (1 .EQ. 0)  .AND.(pplay(i,k).gt.p_lim)) then
          ratio_mod(i,k) = RMD /
     &  (RMD+0.56*(log(pplay(i,k))-log(p_lim))/(log(p_ref)-log(p_lim)))
           ratio_mod(i,k) = max(ratio_mod(i,k),RMD/(RMD+0.56))
           rsmu_mod(i,k) = RD*ratio_mod(i,k)
       endif
      ENDDO
      ENDDO
c----------------------

      if (iflag_pbl.eq.1) then
        do k = 3, klev
          do i = 1, knon
            gamt(i,k)=  -1.0e-03
          enddo
        enddo
        do i = 1, knon
          gamt(i,2) = -2.5e-03
! ADAPTATION GCM POUR CP(T)
          gamt(i,1) = 0.0e0
        enddo
      else
        do k = 1, klev
          do i = 1, knon
            gamt(i,k) = 0.0
          enddo
        enddo
      endif

      DO i = 1, knon
         local_ts(i) = ts(i)
      ENDDO
! ADAPTATION GCM POUR CP(T)
      DO k = 2,klev 
      DO i = 1, knon
            zt(i,k)    = (t(i,k)+t(i,k-1)) * 0.5
c----------------------
c ATMOSPHERE PROFONDE
            zrsmu_mod(i,k)    = (rsmu_mod(i,k)+rsmu_mod(i,k-1)) * 0.5
c----------------------
      ENDDO
      ENDDO
                                                   
c contre-gradient en potentiel:
! ADAPTATION GCM POUR CP(T)
c en fait, les valeurs mises pour gamt sont pour la T pot...
c Donc on garde les memes...
      z_gamah = gamt

c passage en enthalpie potentielle
      call t2tpot(knon*nbp_lev,t,local_h,ppk)
c     print*,"tpot en entree de clqh=",local_h(klon/2,:)

      DO k = 1, klev
      DO i = 1, knon
c h = tpot*cp
         local_h(i,k)= local_h(i,k)*cpdet(t(i,k))
      ENDDO
      ENDDO
c     print*,"enthalpie potentielle en entree de clqh=",
c    .        local_h(klon/2,:)
c
c Convertir les coefficients en variables convenables au calcul:
c
c
      DO k = 2, klev
      DO i = 1, knon
         zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
c----------------------
c ATMOSPHERE PROFONDE
     .                  *(paprs(i,k)/zt(i,k)/zrsmu_mod(i,k))**2
c----------------------
         zx_coef(i,k) = zx_coef(i,k) * dtime*RG
      ENDDO
      ENDDO
c
c Preparer les flux lies aux contre-gardients  OBSOLETE
c
      DO k = 2, klev
      DO i = 1, knon
         zdelz = RD * t(i,k) / RG /paprs(i,k)
     .              *(pplay(i,k-1)-pplay(i,k))
! ADAPTATION GCM POUR CP(T)
         z_gamah(i,k) = z_gamah(i,k)*cpdet(zt(i,k))*zdelz
      ENDDO
      ENDDO
c     print*,"contregradient d(enth pot) en entree de clqh=",
c    .        z_gamah(klon/2,:)

      DO i = 1, knon
! ADAPTATION GCM POUR CP(T)
         zx_buf2(i) = delp(i,klev) + zx_coef(i,klev)
         zx_ch(i,klev) = (local_h(i,klev)*delp(i,klev)
     .                   -zx_coef(i,klev)*z_gamah(i,klev))/zx_buf2(i)
         zx_dh(i,klev) = zx_coef(i,klev) / zx_buf2(i)
      ENDDO
      DO k = klev-1, 2 , -1
      DO i = 1, knon
! ADAPTATION GCM POUR CP(T)
         zx_buf2(i) = delp(i,k)+zx_coef(i,k)
     .               +zx_coef(i,k+1)*(1.-zx_dh(i,k+1))
         zx_ch(i,k) = (local_h(i,k)*delp(i,k)
     .                 +zx_coef(i,k+1)*zx_ch(i,k+1)
     .                 +zx_coef(i,k+1)*z_gamah(i,k+1)
     .                 -zx_coef(i,k)*z_gamah(i,k))/zx_buf2(i)
         zx_dh(i,k) = zx_coef(i,k) / zx_buf2(i)
      ENDDO
      ENDDO
C
C nouvelle formulation JL Dufresne
C
C h1 = zx_ch(i,1) + zx_dh(i,1) * Flux_H(i,1) * dt
C
      DO i = 1, knon
! ADAPTATION GCM POUR CP(T)
         zx_buf2(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dh(i,2))
         zx_ch(i,1) = (local_h(i,1)*delp(i,1)
     .                 +zx_coef(i,2)*(z_gamah(i,2)+zx_ch(i,2)))
     .                /zx_buf2(i)
         zx_dh(i,1) = -1. * RG / zx_buf2(i)
      ENDDO

C Appel a interfsurf (appel generique) routine d'interface avec la surface

c initialisation
       petAcoef =0. 
       petBcoef =0.
       p1lay =0.

c      do i = 1, knon
        petAcoef(1:knon) = zx_ch(1:knon,1)
        petBcoef(1:knon) = zx_dh(1:knon,1)
        tq_cdrag(1:knon) =coef(1:knon,1)
        temp_air(1:knon) =t(1:knon,1)
        epot_air(1:knon) =local_h(1:knon,1)
        pkh1(1:knon)  = ppk(1:knon,1)
     .                 *(paprs(1:knon,1)/pplay(1:knon,1))**RKAPPA
        p1lay(1:knon) = pplay(1:knon,1)
        zlev1(1:knon) = delp(1:knon,1)
        swdown(1:knon) = swnet(1:knon)
c      enddo

! ADAPTATION GCM POUR CP(T)
      CALL interfsurf_hq(itime, dtime, rmu0,
     e klon, nbp_lon, nbp_lat-1, knon, 
     e rlon, rlat, dx, dy, 
     e debut, lafin, soil_model, nsoilmx,tsoil, 
     e zlev1,  u1lay, v1lay, temp_air, epot_air,  
     e tq_cdrag, petAcoef, petBcoef,
     e sollw, sollwdown, swnet, swdown,
     e fder, taux, tauy, 
     e albedo, 
     e ts, pkh1, p1lay, radsol,
     s fluxsens, dflux_s,              
     s tsurf_new, alb_new) 


      do i = 1, knon
        flux_t(i,1) = fluxsens(i)
        d_ts(i) = tsurf_new(i) - ts(i)
        albedo(i) = alb_new(i)
      enddo

c==== une fois on a zx_h_ts, on peut faire l'iteration ========
      DO i = 1, knon
         local_h(i,1) = zx_ch(i,1) + zx_dh(i,1)*flux_t(i,1)*dtime
      ENDDO
      DO k = 2, klev
      DO i = 1, knon
        local_h(i,k) = zx_ch(i,k) + zx_dh(i,k)*local_h(i,k-1)
      ENDDO
      ENDDO
c======================================================================
c== flux_t est le flux de cpt (energie sensible): j/(m**2 s) (+ vers bas)
! ADAPTATION GCM POUR CP(T)
      DO k = 2, klev
      DO i = 1, knon
        flux_t(i,k) = (zx_coef(i,k)/RG/dtime)
     .                * (local_h(i,k)-local_h(i,k-1)+z_gamah(i,k))
      ENDDO
      ENDDO
c======================================================================
C Calcul tendances
! ADAPTATION GCM POUR CP(T)
c     print*,"enthalpie potentielle en sortie de clqh=",
c    .        local_h(klon/2,:)
c tpot = h/cp
      DO k = 1, klev
      DO i = 1, knon
         local_h(i,k) = local_h(i,k)/cpdet(t(i,k))
      ENDDO
      ENDDO
      call tpot2t(knon*nbp_lev,local_h,d_t,ppk)

c     print*,"tpot en sortie de clqh=",local_h(klon/2,:)
c     print*,"T en sortie de clqh=",d_t(klon/2,:)
      DO k = 1, klev
      DO i = 1, knon
         d_t(i,k) = d_t(i,k) - t(i,k)
      ENDDO
      ENDDO
c

      END SUBROUTINE clqh
      
c======================================================================
c======================================================================
c======================================================================
c======================================================================
c======================================================================
c======================================================================

      SUBROUTINE clvent(knon,dtime, u1lay,v1lay,coef,t,ven,
     e                  paprs,pplay,delp,
     s                  d_ven,flux_v)

      use dimphy, only: klon, klev
      use YOMCST_mod, only: RD, RG

      IMPLICIT none
c======================================================================
c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
c Objet: diffusion vertical de la vitesse "ven"
c======================================================================
c Arguments:
c dtime----input-R- intervalle du temps (en second)
c u1lay----input-R- vent u de la premiere couche (m/s)
c v1lay----input-R- vent v de la premiere couche (m/s)
c coef-----input-R- le coefficient d'echange (m**2/s) multiplie par
c                   le cisaillement du vent (dV/dz); la premiere
c                   valeur indique la valeur de Cdrag (sans unite)
c t--------input-R- temperature (K)
c ven------input-R- vitesse horizontale (m/s)
c paprs----input-R- pression a inter-couche (Pa)
c pplay----input-R- pression au milieu de couche (Pa)
c delp-----input-R- epaisseur de couche (Pa)
c
c
c d_ven----output-R- le changement de "ven"
c flux_v---output-R- (diagnostic) flux du vent: (kg m/s)/(m**2 s)
c======================================================================

c Parametres d'entree
      INTEGER, INTENT(IN) :: knon
      REAL, INTENT(IN) :: dtime
      REAL, INTENT(IN) :: u1lay(klon) , v1lay(klon)
      REAL, INTENT(IN) :: coef(klon, klev)
      REAL, INTENT(IN) :: t(klon, klev), ven(klon, klev)
      REAL, INTENT(IN) :: paprs(klon, klev+1), pplay(klon, klev)
      REAL, INTENT(IN) :: delp(klon, klev)
      
c Parametres de sorties
      REAL, INTENT(OUT) :: d_ven(klon, klev)
      REAL, INTENT(OUT) :: flux_v(klon, klev)
c======================================================================
c     include "YOMCST.h"
c======================================================================
c Parametres locaux
      INTEGER i, k
      REAL zx_cv(klon,2:klev)
      REAL zx_dv(klon,2:klev)
      REAL zx_buf(klon)
      REAL zx_coef(klon,klev)
      REAL local_ven(klon,klev)
      REAL zx_alf1(klon), zx_alf2(klon)
c======================================================================
      DO k = 1, klev
      DO i = 1, knon
         local_ven(i,k) = ven(i,k)
      ENDDO
      ENDDO
c======================================================================
      DO i = 1, knon
ccc         zx_alf1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
         zx_alf1(i) = 1.0
         zx_alf2(i) = 1.0 - zx_alf1(i)
         zx_coef(i,1) = coef(i,1)
     .                 * SQRT(u1lay(i)**2+v1lay(i)**2)
     .                 * pplay(i,1)/(RD*t(i,1))
         zx_coef(i,1) = zx_coef(i,1) * dtime*RG
      ENDDO
c======================================================================
      DO k = 2, klev
      DO i = 1, knon
         zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
     .                  *(paprs(i,k)*2/(t(i,k)+t(i,k-1))/RD)**2
         zx_coef(i,k) = zx_coef(i,k) * dtime*RG
      ENDDO
      ENDDO
c======================================================================
      DO i = 1, knon
         zx_buf(i) = delp(i,1) + zx_coef(i,1)*zx_alf1(i)+zx_coef(i,2)
         zx_cv(i,2) = local_ven(i,1)*delp(i,1) / zx_buf(i)
         zx_dv(i,2) = (zx_coef(i,2)-zx_alf2(i)*zx_coef(i,1))
     .                /zx_buf(i)
      ENDDO
      DO k = 3, klev
      DO i = 1, knon
         zx_buf(i) = delp(i,k-1) + zx_coef(i,k)
     .                         + zx_coef(i,k-1)*(1.-zx_dv(i,k-1))
         zx_cv(i,k) = (local_ven(i,k-1)*delp(i,k-1)
     .                  +zx_coef(i,k-1)*zx_cv(i,k-1) )/zx_buf(i)
         zx_dv(i,k) = zx_coef(i,k)/zx_buf(i)
      ENDDO
      ENDDO
      DO i = 1, knon
         local_ven(i,klev) = ( local_ven(i,klev)*delp(i,klev)
     .                        +zx_coef(i,klev)*zx_cv(i,klev) )
     .                   / ( delp(i,klev) + zx_coef(i,klev)
     .                       -zx_coef(i,klev)*zx_dv(i,klev) )
      ENDDO
      DO k = klev-1, 1, -1
      DO i = 1, knon
         local_ven(i,k) = zx_cv(i,k+1) + zx_dv(i,k+1)*local_ven(i,k+1)
      ENDDO
      ENDDO
c======================================================================
c== flux_v est le flux de moment angulaire (positif vers bas)
c== dont l'unite est: (kg m/s)/(m**2 s)
      DO i = 1, knon
         flux_v(i,1) = zx_coef(i,1)/(RG*dtime)
     .                 *(local_ven(i,1)*zx_alf1(i)
     .                  +local_ven(i,2)*zx_alf2(i))
      ENDDO
      DO k = 2, klev
      DO i = 1, knon
         flux_v(i,k) = zx_coef(i,k)/(RG*dtime)
     .               * (local_ven(i,k)-local_ven(i,k-1))
      ENDDO
      ENDDO
c
      DO k = 1, klev
      DO i = 1, knon
         d_ven(i,k) = local_ven(i,k) - ven(i,k)
      ENDDO
      ENDDO
c
      END SUBROUTINE clvent
      
c======================================================================
c======================================================================
c======================================================================
c======================================================================
c======================================================================
c======================================================================

      SUBROUTINE coefkz(knon, paprs, pplay, ppk,
     .                  ts,u,v,t,
     .                  pcfm, pcfh)

      use dimphy, only: klon, klev
      use cpdet_phy_mod, only: cpdet,t2tpot
      use clesphys_mod, only: ksta, lmixmin
#ifdef CPP_XIOS      
      use xios_output_mod, only: send_xios_field
#endif
      use YOMCST_mod, only: R, RD, RG, RKAPPA
      use print_control_mod, only: lunout, prt_level
      IMPLICIT none
c======================================================================
c Auteur(s) F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) date: 19930922
c           (une version strictement identique a l'ancien modele)
c Objet: calculer le coefficient du frottement du sol (Cdrag) et les
c        coefficients d'echange turbulent dans l'atmosphere.
c Arguments:
c knon-----input-I- nombre de points a traiter
c paprs----input-R- pression a chaque intercouche (en Pa)
c pplay----input-R- pression au milieu de chaque couche (en Pa)
c ts-------input-R- temperature du sol (en Kelvin)
c u--------input-R- vitesse u
c v--------input-R- vitesse v
c t--------input-R- temperature (K)
c
c itop-----output-I- numero de couche du sommet de la couche limite
c pcfm-----output-R- coefficients a calculer (vitesse)
c pcfh-----output-R- coefficients a calculer (chaleur et humidite)
c======================================================================

c
c Arguments:
c
c Parametres d'entree
      INTEGER, INTENT(IN) :: knon
      REAL, INTENT(IN) :: ts(klon)
      REAL, INTENT(IN) :: pplay(klon, klev)
      REAL, INTENT(IN) ::paprs(klon, klev+1)
! ADAPTATION GCM POUR CP(T)
      REAL, INTENT(IN) :: ppk(klon, klev)
      REAL, INTENT(IN) :: u(klon, klev), v(klon, klev), t(klon, klev)
      
c Parametre de sorties
      REAL, INTENT(OUT) :: pcfm(klon, klev), pcfh(klon, klev)
c
c Quelques constantes et options:
c
      REAL, PARAMETER :: cepdu2=((1.e-5)**2)
c     REAL, PARAMETER :: cepdu2 =(0.1)**2
      REAL, PARAMETER :: ckap=0.4
      REAL, PARAMETER :: cb=5.0
      REAL, PARAMETER :: cc=5.0
      REAL, PARAMETER :: cd=5.0
      REAL, PARAMETER :: clam=160.0
      REAL, PARAMETER :: ric=0.4 ! nombre de Richardson critique
      REAL, PARAMETER :: prandtl=0.4
      INTEGER isommet ! le sommet de la couche limite
      
      LOGICAL, PARAMETER :: tvirtu=.TRUE. ! calculer Ri d'une maniere plus performante
      LOGICAL, PARAMETER :: opt_ec=.FALSE. ! formule du Centre Europeen dans l'atmosphere
      
c      REAL cepdu2, ckap, cb, cc, cd, clam
c TEST VENUS
c     PARAMETER (cepdu2 =(0.1)**2)
c      PARAMETER (cepdu2 =(1.e-5)**2)
c      PARAMETER (CKAP=0.4)
c      PARAMETER (cb=5.0)
c      PARAMETER (cc=5.0)
c      PARAMETER (cd=5.0)
c      PARAMETER (clam=160.0)
c      REAL ric ! nombre de Richardson critique
c      PARAMETER(ric=0.4)
c      REAL prandtl
c      PARAMETER (prandtl=0.4)
c      INTEGER isommet ! le sommet de la couche limite

c      LOGICAL tvirtu ! calculer Ri d'une maniere plus performante
c      PARAMETER (tvirtu=.TRUE.)
c      LOGICAL opt_ec ! formule du Centre Europeen dans l'atmosphere
c      PARAMETER (opt_ec=.FALSE.)

c
c Variables locales:
c
      INTEGER itop(klon)
      INTEGER i, k
      REAL zgeop(klon,klev)
! ADAPTATION GCM POUR CP(T)
      REAL zmgeom(klon,klev),zpk(klon,klev)
      REAL zt(klon,klev),ztvu(klon,klev),ztvd(klon,klev)
      real ztetav(klon,klev),ztetavu(klon,klev),ztetavd(klon,klev)
      REAL zri(klon,klev),zri1(klon),z1(klon)
      REAL pcfm1(klon), pcfh1(klon)
c
      REAL zdphi, zdu2, zcdn, zl2
      REAL zscf
      REAL zdelta, zcvm5, zcor
      REAL z2geomf, zalh2, zalm2, zscfh, zscfm
cIM
      LOGICAL, PARAMETER :: check=.false.
c
c contre-gradient pour la chaleur sensible: Kelvin/metre
      REAL gamt(2:klev)
      REAL gamh(2:klev)
c
      LOGICAL, SAVE :: appel1er = .TRUE.
C$OMP THREADPRIVATE(appel1er)
c
c Fonctions thermodynamiques et fonctions d'instabilite
      REAL fsta, fins, x
      LOGICAL, PARAMETER :: zxli=.FALSE. ! utiliser un jeu de fonctions simples PARAMETER (zxli=.FALSE.)
c
      fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))
      fins(x) = SQRT(1.0-18.0*x)
c
c----------------------
c ATMOSPHERE PROFONDE
      real p_lim,p_ref
      real modgam(klon,klev)

      p_lim = 6e6
      p_ref = 8.9e6

      DO k = 1, klev
      DO i = 1, klon
        modgam(i,k) = 0.
c ATM PROFONDE DESACTIVEE !
       if (    (1 .EQ. 0)  .AND.(pplay(i,k).gt.p_lim)) then
           modgam(i,k) = 0.56E-3/(log(p_ref)-log(p_lim))/R
       endif
      ENDDO
      ENDDO
c----------------------

      isommet=klev

      IF (appel1er) THEN
        if (prt_level > 9) THEN
          WRITE(lunout,*)'coefkz, opt_ec:', opt_ec
          WRITE(lunout,*)'coefkz, isommet:', isommet
          WRITE(lunout,*)'coefkz, tvirtu:', tvirtu
          appel1er = .FALSE.
        endif
      ENDIF
c
c Initialiser les sorties
c
      DO k = 1, klev
      DO i = 1, knon
         pcfm(i,k) = 0.0
         pcfh(i,k) = 0.0
      ENDDO
      ENDDO
      DO i = 1, knon
         itop(i) = 0
      ENDDO
c
c Prescrire la valeur de contre-gradient
c
      if (iflag_pbl.eq.1) then
         DO k = 3, klev
            gamt(k) = -1.0E-03
         ENDDO
         gamt(2) = -2.5E-03
      else
         DO k = 2, klev
            gamt(k) = 0.0
         ENDDO
      ENDIF

c
c Calculer les geopotentiels de chaque couche
c
      DO i = 1, knon
         zgeop(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))
     .                   * (paprs(i,1)-pplay(i,1))
      ENDDO
      DO k = 2, klev
      DO i = 1, knon
         zgeop(i,k) = zgeop(i,k-1)
     .              + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)
     .                   * (pplay(i,k-1)-pplay(i,k))
      ENDDO
      ENDDO
c
c Calculer les coefficients turbulents dans l'atmosphere
! ADAPTATION GCM POUR CP(T)
c tout a ete modifie...
c

      DO k = 2,klev 
      DO i = 1, knon
            zt(i,k)    = (t(i,k)+t(i,k-1)) * 0.5
            zmgeom(i,k)= zgeop(i,k)-zgeop(i,k-1)
            zdphi      = zmgeom(i,k)/2.
c----------------------
c ATMOSPHERE PROFONDE
            ztvd(i,k)  = t(i,k)   + zdphi*
     &                      (1./cpdet(zt(i,k))+modgam(i,k))
            ztvu(i,k)  = t(i,k-1) - zdphi*
     &                      (1./cpdet(zt(i,k))+modgam(i,k))
c----------------------
            zpk(i,k)   = ppk(i,k)*(paprs(i,k)/pplay(i,k))**RKAPPA
      ENDDO
      ENDDO
      DO i = 1, knon
        itop(i) = isommet
        zt(i,1)   = ts(i)
        ztvu(i,1) = ts(i) 
c----------------------
c ATMOSPHERE PROFONDE
        ztvd(i,1) = t(i,1)+zgeop(i,1)*
     &                      (1./cpdet(zt(i,1))+modgam(i,1))
c----------------------
        zpk(i,1)  = ppk(i,1)*(paprs(i,1)/pplay(i,1))**RKAPPA
      ENDDO
      call t2tpot(klon*klev,zt,ztetav,zpk)
      call t2tpot(klon*klev,ztvu,ztetavu,zpk)
      call t2tpot(klon*klev,ztvd,ztetavd,zpk)

      DO k = 2, isommet
        DO i = 1, knon
            zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2
     .                     +(v(i,k)-v(i,k-1))**2)
c contre-gradient en potentiel:
! ADAPTATION GCM POUR CP(T)
c en fait, les valeurs mises pour gamt sont pour la T pot...
c Donc on garde les memes...
            gamh(k) = gamt(k)
c
c           calculer le nombre de Richardson:
c
            IF (tvirtu) THEN
            zri(i,k) =(  zmgeom(i,k)*(ztetavd(i,k)-ztetavu(i,k))
     .            + zmgeom(i,k)*zmgeom(i,k)/RG*gamh(k))   ! contregradient
     .           /(zdu2*ztetav(i,k))
c
            ELSE ! calcul de Richardson compatible LMD5
        print*,"calcul de Richardson sans tvirtu..." 
        print*,"Pas prevu... A revoir"
        stop
            ENDIF
c
c           finalement, les coefficients d'echange sont obtenus:
c
            zcdn=SQRT(zdu2) / zmgeom(i,k) * RG
c
          IF (opt_ec) THEN
            z2geomf=zgeop(i,k-1)+zgeop(i,k)
            zalm2=(0.5*ckap/RG*z2geomf
     .             /(1.+0.5*ckap/rg/clam*z2geomf))**2
            zalh2=(0.5*ckap/rg*z2geomf
     .             /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2
            IF (zri(i,k).LT.0.0) THEN  ! situation instable
               zscf = ((zgeop(i,k)/zgeop(i,k-1))**(1./3.)-1.)**3
     .                / (zmgeom(i,k)/RG)**3 / (zgeop(i,k-1)/RG)
               zscf = SQRT(-zri(i,k)*zscf)
               zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf)
               zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf)
               pcfm(i,k)=zcdn*zalm2*(1.-2.0*cb*zri(i,k)*zscfm)
               pcfh(i,k)=zcdn*zalh2*(1.-3.0*cb*zri(i,k)*zscfh)
            ELSE ! situation stable
               zscf=SQRT(1.+cd*zri(i,k))
               pcfm(i,k)=zcdn*zalm2/(1.+2.0*cb*zri(i,k)/zscf)
               pcfh(i,k)=zcdn*zalh2/(1.+3.0*cb*zri(i,k)*zscf)
            ENDIF
          ELSE
            zl2=(lmixmin*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1))
     .                          /(paprs(i,2)-paprs(i,itop(i)+1)) ))**2
            pcfm(i,k)=sqrt(max(zcdn*zcdn*(ric-zri(i,k))/ric, ksta))
            pcfm(i,k)= zl2* pcfm(i,k)
            pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
          ENDIF
        ENDDO
      ENDDO
c Richardson au sol:
      DO i = 1, knon
            zdu2=MAX(cepdu2,u(i,1)**2+v(i,1)**2)
            zri(i,1) = zgeop(i,1)*(ztetavd(i,1)-ztetavu(i,1))
     .              /(zdu2*ztetav(i,1))
      ENDDO
c
c Calculer le frottement au sol (Cdrag)
! ADAPTATION GCM POUR CP(T)
c
      DO i = 1, knon
       z1(i) = zgeop(i,1)
       zri1(i) = zri(i,1)
      ENDDO
c
      CALL clcdrag(klon, knon, zxli, 
     $             z1, zri1,
     $             pcfm1, pcfh1) 
C
      DO i = 1, knon
       pcfm(i,1)=pcfm1(i)
       pcfh(i,1)=pcfh1(i)
      ENDDO
c
c Au-dela du sommet, pas de diffusion turbulente:
c
      DO i = 1, knon
         IF (itop(i)+1 .LE. klev) THEN
            DO k = itop(i)+1, klev
               pcfh(i,k) = 0.0
               pcfm(i,k) = 0.0
            ENDDO
         ENDIF
      ENDDO
c
c VENUS TEST :
c      pcfm(:,:)= 0.15
c      pcfh(:,:)= 0.15
c
c VENUS TEST : frottement de surface beaucoup plus grand
c      pcfm(:,1)= pcfm(:,1)*20.
c      pcfh(:,1)= pcfh(:,1)*20.

c#ifdef CPP_XIOS      
c     CALL send_xios_field("Ri",zri)
c#endif

      END SUBROUTINE coefkz

C=================================================================
C=================================================================
C=================================================================
C=================================================================

      SUBROUTINE coefkz2(knon, paprs, pplay,t,
     .                  pcfm, pcfh)

      use dimphy, only: klon, klev
      use mod_grid_phy_lmdz, only: nbp_lev
      use cpdet_phy_mod, only: cpdet
      use YOMCST_mod, only: RD
      IMPLICIT none
c======================================================================
c J'introduit un peu de diffusion sauf dans les endroits
c ou une forte inversion est presente
c On peut dire qu'il represente la convection peu profonde
c
c Arguments:
c knon-----input-I- nombre de points a traiter
c paprs----input-R- pression a chaque intercouche (en Pa)
c pplay----input-R- pression au milieu de chaque couche (en Pa)
c t--------input-R- temperature (K)
c
c pcfm-----output-R- coefficients a calculer (vitesse)
c pcfh-----output-R- coefficients a calculer (chaleur et humidite)
c======================================================================
c
c Arguments:
c
      INTEGER,INTENT(IN) :: knon
      REAL,INTENT(IN) :: paprs(klon,klev+1)
      REAL,INTENT(IN) :: pplay(klon,klev)
      REAL,INTENT(IN) :: t(klon,klev)
c
      REAL,INTENT(OUT) :: pcfm(klon,klev), pcfh(klon,klev)
c
c Variables locales:
c
      INTEGER i, k, invb(knon)
      REAL zl2(knon), zt
      REAL zdthmin(knon), zdthdp
c
c Initialiser les sorties
c
      DO k = 1, klev
      DO i = 1, knon
         pcfm(i,k) = 0.0
         pcfh(i,k) = 0.0
      ENDDO
      ENDDO
c
c Chercher la zone d'inversion forte
c
      DO i = 1, knon
         invb(i) = klev
         zdthmin(i)=0.0
      ENDDO
      DO k = 2, klev/2-1
      DO i = 1, knon
! ADAPTATION GCM POUR CP(T)
         zt = 0.5*(t(i,k)+t(i,k+1))
         zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1))
     .          - RD * zt/cpdet(zt)/paprs(i,k+1)
         zdthdp = zdthdp * 100.0
         IF (pplay(i,k).GT.0.8*paprs(i,1) .AND.
     .       zdthdp.LT.zdthmin(i) ) THEN
            zdthmin(i) = zdthdp
            invb(i) = k
         ENDIF
      ENDDO
      ENDDO
c
      END SUBROUTINE coefkz2

      END MODULE clmain_mod