!
! $Header: /home/cvsroot/LMDZ4/libf/phylmd/clmain.F,v 1.3 2005/02/07 16:41:35 fairhead Exp $
!
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 Tout ce qui a trait au traceurs est dans phytrac maintenant
cAA pour l'instant le calcul de la couche limite pour les traceurs
cAA se fait avec cltrac et ne tient pas compte de la differentiation
cAA des sous-fraction de sol.
cAA REM bis :
cAA----------
cAA Pour pouvoir extraire les coefficient d'echanges et le vent 
cAA dans la premiere couche, 3 champs supplementaires ont ete crees
cAA zcoefh,zu1 et zv1. Pour l'instant nous avons moyenne les valeurs
cAA de ces trois champs sur les 4 subsurfaces du modele. Dans l'avenir 
cAA si les informations des subsurfaces doivent etre prises en compte
cAA il faudra sortir ces memes champs en leur ajoutant une dimension, 
cAA c'est a dire nbsrf (nbre de subsurface).
!      USE ioipsl
!      USE interface_surf
      use dimphy, only: klon, klev
      use mod_grid_phy_lmdz, only: nbp_lev
      use cpdet_phy_mod, only: t2tpot
      use turb_mod, only :yustar
      
#ifdef CPP_XIOS      
      use xios_output_mod, only: send_xios_field
#endif      
      
      IMPLICIT none
c======================================================================
c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
c Objet: interface de "couche limite" (diffusion verticale)
c Arguments:
c dtime----input-R- interval du temps (secondes)
c itap-----input-I- numero du pas de temps
c t--------input-R- temperature (K)
c u--------input-R- vitesse u
c v--------input-R- vitesse v
c ts-------input-R- temperature du sol (en Kelvin)
c paprs----input-R- pression a intercouche (Pa)
c pplay----input-R- pression au milieu de couche (Pa)
c radsol---input-R- flux radiatif net (positif vers le sol) en W/m**2
c rlat-----input-R- latitude en degree
c dx-----input-R- resolution des mailles en x (m)
c dy-----input-R- resolution des mailles en y (m)
c
c q2-----inoutput-R- $q^2$ TKE at inter-layers
c
c d_t------output-R- le changement pour "t"
c d_u------output-R- le changement pour "u"
c d_v------output-R- le changement pour "v"
c d_ts-----output-R- le changement pour "ts"
c flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)
c                    (orientation positive vers le bas)
c flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
c flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal
c dflux_t derive du flux sensible
cAA on rajoute en output yu1 et yv1 qui sont les vents dans 
cAA la premiere couche
c======================================================================
c$$$ PB ajout pour soil
      include "dimsoil.h"
      include "iniprint.h"
      include "clesphys.h"
      include "compbl.h"
c
c Parametres d'entree
      REAL, INTENT(IN) :: dtime
      INTEGER, INTENT(IN) :: itap
      REAL, INTENT(IN) :: t(klon,klev)
      REAL, INTENT(IN) :: u(klon,klev), v(klon,klev)
      REAL, INTENT(IN) :: paprs(klon,klev+1), pplay(klon,klev)
! ADAPTATION GCM POUR CP(T)
      REAL, INTENT(IN) :: ppk(klon,klev)
      REAL, INTENT(IN) :: rlon(klon), rlat(klon), dx(klon), dy(klon)
      REAL, INTENT(IN) :: rmu0(klon)         ! cosine of solar zenith angle
      LOGICAL, INTENT(IN) :: debut, lafin
      REAL, INTENT(IN) :: ts(klon)
      REAL, INTENT(IN) :: sollw(klon), solsw(klon), sollwdown(klon)
      
c Paramètres IN/OUT
      REAL, INTENT(INOUT) :: fder(klon)
      REAL, INTENT(INOUT) :: flux_u(klon,klev), flux_v(klon,klev)
      REAL, INTENT(INOUT) :: radsol(klon)
      REAL, INTENT(INOUT) :: q2(klon,klev+1)
      
c Parametres de sorties
      REAL, INTENT(OUT) :: d_t(klon, klev)
      REAL, INTENT(OUT) :: d_u(klon, klev), d_v(klon, klev)
      REAL, INTENT(OUT) :: flux_t(klon,klev)
      REAL, INTENT(OUT) :: dflux_t(klon)
      REAL, INTENT(OUT) :: cdragh(klon), cdragm(klon)
      REAL, INTENT(OUT) :: d_ts(klon)
      REAL, INTENT(OUT) :: albe(klon)
cAA
      REAL, INTENT(OUT) :: zcoefh(klon,klev)
      REAL, INTENT(OUT) :: zu1(klon)
      REAL, INTENT(OUT) :: zv1(klon)
cAA
c$$$ PB ajout pour soil
      REAL, INTENT(INOUT) :: ftsoil(klon,nsoilmx)

c======================================================================
      EXTERNAL clqh, clvent, coefkz
c======================================================================
c Parametre locaux
      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
      include "YOMCST.h"
      REAL u1lay(klon), v1lay(klon)
      REAL delp(klon,klev)
      INTEGER i, k
      INTEGER ni(klon), knon, j     
c======================================================================
      REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.
c======================================================================
c
      LOGICAL zxli ! utiliser un jeu de fonctions simples
      PARAMETER (zxli=.FALSE.)
c
      REAL zt, zdelta, zcor
C
      character (len = 20) :: modname = 'clmain'
      LOGICAL check
      PARAMETER (check=.false.)
C
      if (check) THEN
          write(*,*) modname,'  klon=',klon
          call flush(6)
      endif
          
      DO k = 1, klev   ! epaisseur de couche
      DO i = 1, klon
         delp(i,k) = paprs(i,k)-paprs(i,k+1)
      ENDDO
      ENDDO
      DO i = 1, klon  ! vent de la premiere couche
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 chercher les indices:
      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$$$ PB ajour pour 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 calculer Cdrag et les coefficients d'echange
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c-------------------------------------------------
c  Calcul anciens du LMD. 
c  dans les 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 POUR CP(T)
           print*," coefkzmin: ADAPTATION NON FAITE..."
cIM cf FH: 201103 BEG

c   Calcul d'une diffusion minimale pour les conditions tres stables.
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 adapte a Mars Richard Fournier et 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 POUR 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


c   Bug introduit volontairement pour converger avec les resultats
c  du papier sur les thermiques.
         if (1.eq.1) then
         y_cd_m(1:knon) = ycoefm(1:knon,1)
         y_cd_h(1:knon) = ycoefh(1:knon,1)
         else
         y_cd_h(1:knon) = ycoefm(1:knon,1)
         y_cd_m(1:knon) = ycoefh(1:knon,1)
         endif

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

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

c   iflag_pbl peut etre utilise comme longueur de melange

         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 calculer la diffusion des vitesses "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 pour le couplage
      ytaux = y_flux_u(:,1)
      ytauy = y_flux_v(:,1)

c FH modif sur le cdrag temperature
c$$$PB : déplace dans clcdrag
c$$$      do i=1,knon
c$$$         ycoefh(i,1)=ycoefm(i,1)*0.8
c$$$      enddo

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c calculer la diffusion de "q" et de "h"
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! ADAPTATION GCM POUR 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: tendance sur Tsurf forcee = 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$$$ PB ajout pour soil
      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

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 mod_grid_phy_lmdz, only: nbp_lon, nbp_lat, nbp_lev
      use cpdet_phy_mod, only: t2tpot,tpot2t,cpdet

      IMPLICIT none
c======================================================================
c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
c Objet: diffusion verticale de "h"
c======================================================================
      include "YOMCST.h"
      include "dimsoil.h"
      include "iniprint.h"

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======================================================================
      include "compbl.h"
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 tsol_rad(klon), tsurf_new(klon), alb_new(klon)

      character (len = 20) :: modname = 'Debut clqh'
      LOGICAL check
      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 tsol_rad, 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
      
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
      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======================================================================
      include "iniprint.h"

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======================================================================
      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
      
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
#ifdef CPP_XIOS      
      use xios_output_mod, only: send_xios_field
#endif

      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======================================================================
      include "YOMCST.h"
      include "iniprint.h"
      include "compbl.h"
      include "clesphys.h"
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
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

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
      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======================================================================
      include "YOMCST.h"
      include "iniprint.h"
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