!
! $Header$
!
c
c
      SUBROUTINE clmain(dtime,itap,date0,pctsrf,pctsrf_new,
     .                  t,q,u,v,
     .                  jour, rmu0, co2_ppm,
     .                  ok_veget, ocean, npas, nexca, ts,
     .                  soil_model,cdmmax, cdhmax,
     .                  ksta, ksta_ter, ok_kzmin, ftsoil,qsol,
cIM BAD .                  paprs,pplay,radsol,snow,qsurf,evap,albe,alblw,
     .                  paprs,pplay,snow,qsurf,evap,albe,alblw,
     .                  fluxlat,
     .                  rain_f, snow_f, solsw, sollw, sollwdown, fder,
     .                  rlon, rlat, cufi, cvfi, rugos,
     .                  debut, lafin, agesno,rugoro,
     .                  d_t,d_q,d_u,d_v,d_ts,
     .                  flux_t,flux_q,flux_u,flux_v,cdragh,cdragm,
     .                  q2,
     .                  dflux_t,dflux_q,
     .                  zcoefh,zu1,zv1, t2m, q2m, u10m, v10m,
cIM cf. AM : pbl
     .                  pblh,capCL,oliqCL,cteiCL,pblT,
     .                  therm,trmb1,trmb2,trmb3,plcl,
     .                  fqcalving,ffonte, run_off_lic_0,
cIM "slab" ocean
     .                  flux_o, flux_g, tslab, seaice)
cAA .                  itr, tr, flux_surf, d_tr)
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
      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 date0----input-R- jour initial
c t--------input-R- temperature (K)
c q--------input-R- vapeur d'eau (kg/kg)
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 rugos----input-R- longeur de rugosite (en m)
c cufi-----input-R- resolution des mailles en x (m)
c cvfi-----input-R- resolution des mailles en y (m)
c
c d_t------output-R- le changement pour "t"
c d_q------output-R- le changement pour "q"
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_q---output-R- flux de vapeur d'eau (kg/m**2/s)
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
c dflux_q derive du flux latent
cIM "slab" ocean
c flux_g---output-R-  flux glace (pour OCEAN='slab  ')
c flux_o---output-R-  flux ocean (pour OCEAN='slab  ')
c tslab-in/output-R temperature du slab ocean (en Kelvin) ! uniqmnt pour slab
c seaice---output-R-  glace de mer (kg/m2) (pour OCEAN='slab  ')
ccc
c ffonte----Flux thermique utilise pour fondre la neige
c fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la
c           hauteur de neige, en kg/m2/s
cAA on rajoute en output yu1 et yv1 qui sont les vents dans 
cAA la premiere couche
cAA ces 4 variables sont maintenant traites dans phytrac
c itr--------input-I- nombre de traceurs
c tr---------input-R- q. de traceurs
c flux_surf--input-R- flux de traceurs a la surface
c d_tr-------output-R tendance de traceurs
cIM cf. AM : PBL
c trmb1-------deep_cape
c trmb2--------inhibition 
c trmb3-------Point Omega
c Cape(klon)-------Cape du thermique
c EauLiq(klon)-------Eau liqu integr du thermique
c ctei(klon)-------Critere d'instab d'entrainmt des nuages de CL
c lcl------- Niveau de condensation
c pblh------- HCL
c pblT------- T au nveau HCL
c======================================================================
#include "dimensions.h"
#include "dimphy.h"
#include "indicesol.h"
c$$$ PB ajout pour soil
#include "dimsoil.h"
#include "iniprint.h"
#include "compbl.h"
c
      REAL dtime
      real date0
      integer itap
      REAL t(klon,klev), q(klon,klev)
      REAL u(klon,klev), v(klon,klev)
cIM 230604 BAD  REAL radsol(klon) ??? 
      REAL paprs(klon,klev+1), pplay(klon,klev)
      REAL rlon(klon), rlat(klon), cufi(klon), cvfi(klon)
      REAL d_t(klon, klev), d_q(klon, klev)
      REAL d_u(klon, klev), d_v(klon, klev)
      REAL flux_t(klon,klev, nbsrf), flux_q(klon,klev, nbsrf)
      REAL dflux_t(klon), dflux_q(klon)
cIM "slab" ocean
      REAL flux_o(klon), flux_g(klon)
      REAL y_flux_o(klon), y_flux_g(klon)
      REAL tslab(klon), ytslab(klon)
      REAL seaice(klon), y_seaice(klon)
cIM cf JLD
      REAL y_fqcalving(klon), y_ffonte(klon)
      REAL fqcalving(klon,nbsrf), ffonte(klon,nbsrf)
      REAL run_off_lic_0(klon), y_run_off_lic_0(klon)

      REAL flux_u(klon,klev, nbsrf), flux_v(klon,klev, nbsrf)
      REAL rugmer(klon), agesno(klon,nbsrf),rugoro(klon)
      REAL cdragh(klon), cdragm(klon)
      integer jour            ! jour de l'annee en cours
      real rmu0(klon)         ! cosinus de l'angle solaire zenithal
      REAL co2_ppm            ! taux CO2 atmosphere
      LOGICAL debut, lafin, ok_veget
      character*6 ocean
      integer npas, nexca
cAA      INTEGER itr
cAA      REAL tr(klon,klev,nbtr)
cAA      REAL d_tr(klon,klev,nbtr)
cAA      REAL flux_surf(klon,nbtr)
c
      REAL pctsrf(klon,nbsrf)
      REAL ts(klon,nbsrf)
      REAL d_ts(klon,nbsrf)
      REAL snow(klon,nbsrf)
      REAL qsurf(klon,nbsrf)
      REAL evap(klon,nbsrf)
      REAL albe(klon,nbsrf)
      REAL alblw(klon,nbsrf)
c$$$ PB
      REAL fluxlat(klon,nbsrf)
C
      real rain_f(klon), snow_f(klon)
      REAL fder(klon)
cIM cf. JLD   REAL sollw(klon), solsw(klon), sollwdown(klon)
      REAL sollw(klon,nbsrf), solsw(klon,nbsrf), sollwdown(klon)
      REAL rugos(klon,nbsrf)
C la nouvelle repartition des surfaces sortie de l'interface
      REAL pctsrf_new(klon,nbsrf)
cAA
      REAL zcoefh(klon,klev)
      REAL zu1(klon)
      REAL zv1(klon)
cAA
c$$$ PB ajout pour soil
      LOGICAL soil_model
cIM ajout seuils cdrm, cdrh
      REAL cdmmax, cdhmax
cIM: 261103
      REAL ksta, ksta_ter
      LOGICAL ok_kzmin
cIM: 261103
      REAL ftsoil(klon,nsoilmx,nbsrf)
      REAL ytsoil(klon,nsoilmx)
      REAL qsol(klon)
c======================================================================
      EXTERNAL clqh, clvent, coefkz, calbeta, cltrac
c======================================================================
      REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)
      REAL yalb(klon)
      REAL yalblw(klon)
      REAL yu1(klon), yv1(klon)
      real ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)
      real yrain_f(klon), ysnow_f(klon)
      real ysollw(klon), ysolsw(klon), ysollwdown(klon)
      real yfder(klon), ytaux(klon), ytauy(klon)
      REAL yrugm(klon), yrads(klon),yrugoro(klon)
c$$$ PB
      REAL yfluxlat(klon)
C
      REAL y_d_ts(klon)
      REAL y_d_t(klon, klev), y_d_q(klon, klev)
      REAL y_d_u(klon, klev), y_d_v(klon, klev)
      REAL y_flux_t(klon,klev), y_flux_q(klon,klev)
      REAL y_flux_u(klon,klev), y_flux_v(klon,klev)
      REAL y_dflux_t(klon), y_dflux_q(klon)
      REAL ycoefh(klon,klev), ycoefm(klon,klev)
      REAL yu(klon,klev), yv(klon,klev)
      REAL yt(klon,klev), yq(klon,klev)
      REAL ypaprs(klon,klev+1), ypplay(klon,klev), ydelp(klon,klev)
cAA      REAL ytr(klon,klev,nbtr)
cAA      REAL y_d_tr(klon,klev,nbtr)
cAA      REAL yflxsrf(klon,nbtr)
c
c
      LOGICAL ok_nonloc
      PARAMETER (ok_nonloc=.FALSE.)
      REAL ycoefm0(klon,klev), ycoefh0(klon,klev)

cIM 081204 hcl_Anne ? BEG
      real yzlay(klon,klev),yzlev(klon,klev+1),yteta(klon,klev)
      real ykmm(klon,klev+1),ykmn(klon,klev+1)
      real ykmq(klon,klev+1)
      real yq2(klon,klev+1),q2(klon,klev+1,nbsrf)
      real q2diag(klon,klev+1)
cIM 081204   real yustar(klon),y_cd_m(klon),y_cd_h(klon)
cIM 081204 hcl_Anne ? END
c
#include "YOMCST.h"
#include "YOETHF.h"
#include "FCTTRE.h"
      REAL u1lay(klon), v1lay(klon)
      REAL delp(klon,klev)
      INTEGER i, k, nsrf 
cAA   INTEGER it
      INTEGER ni(klon), knon, j
c Introduction d'une variable "pourcentage potentiel" pour tenir compte
c des eventuelles apparitions et/ou disparitions de la glace de mer
      REAL pctsrf_pot(klon,nbsrf)
      
c======================================================================
      REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.
c======================================================================
c
c maf pour sorties IOISPL en cas de debugagage
c
      CHARACTER*80 cldebug
      SAVE cldebug
      CHARACTER*8 cl_surf(nbsrf)
      SAVE cl_surf
      INTEGER nhoridbg, nidbg
      SAVE nhoridbg, nidbg
      INTEGER ndexbg(iim*(jjm+1))
      REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1), zjulian
      REAL tabindx(klon)
      REAL debugtab(iim,jjm+1)
      LOGICAL first_appel
      SAVE first_appel
      DATA first_appel/.true./
      LOGICAL debugindex
      SAVE debugindex
      DATA debugindex/.false./
      integer idayref
#include "temps.h"
      REAL t2m(klon,nbsrf), q2m(klon,nbsrf)
      REAL u10m(klon,nbsrf), v10m(klon,nbsrf)
c
      REAL yt2m(klon), yq2m(klon), yu10m(klon)
      REAL yustar(klon)
c -- LOOP
       REAL yu10mx(klon)
       REAL yu10my(klon)
       REAL ywindsp(klon)
c -- LOOP
c
      REAL yt10m(klon), yq10m(klon)
cIM cf. AM : pbl, hbtm2 (Comme les autres diagnostics on cumule ds physic ce qui 
c   permet de sortir les grdeurs par sous surface)
      REAL pblh(klon,nbsrf)
      REAL plcl(klon,nbsrf)
      REAL capCL(klon,nbsrf)
      REAL oliqCL(klon,nbsrf)
      REAL cteiCL(klon,nbsrf)
      REAL pblT(klon,nbsrf)
      REAL therm(klon,nbsrf)
      REAL trmb1(klon,nbsrf)
      REAL trmb2(klon,nbsrf)
      REAL trmb3(klon,nbsrf)
      REAL ypblh(klon)
      REAL ylcl(klon)
      REAL ycapCL(klon)
      REAL yoliqCL(klon)
      REAL ycteiCL(klon)
      REAL ypblT(klon)
      REAL ytherm(klon)
      REAL ytrmb1(klon)
      REAL ytrmb2(klon)
      REAL ytrmb3(klon)
      REAL y_cd_h(klon), y_cd_m(klon)
c     REAL ygamt(klon,2:klev) ! contre-gradient pour temperature
c     REAL ygamq(klon,2:klev) ! contre-gradient pour humidite
      REAL uzon(klon), vmer(klon)
      REAL tair1(klon), qair1(klon), tairsol(klon)
      REAL psfce(klon), patm(klon)
c
      REAL qairsol(klon), zgeo1(klon)
      REAL rugo1(klon)
c
      LOGICAL zxli ! utiliser un jeu de fonctions simples
      PARAMETER (zxli=.FALSE.)
c
      REAL zt, zqs, zdelta, zcor
      REAL t_coup
      PARAMETER(t_coup=273.15)
C
      character (len = 20) :: modname = 'clmain'
      LOGICAL check
      PARAMETER (check=.false.)
C
      if (check) THEN
          write(*,*) modname,'  klon=',klon
          call flush(6)
      endif
      IF (debugindex .and. first_appel) THEN
          first_appel=.false.
!
! initialisation sorties netcdf
!
          idayref = day_ini
          CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
          CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,zx_lon)
          DO i = 1, iim
            zx_lon(i,1) = rlon(i+1)
            zx_lon(i,jjm+1) = rlon(i+1)
          ENDDO
          CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,zx_lat)
          cldebug='sous_index'
          CALL histbeg(cldebug, iim,zx_lon(:,1),jjm+1,zx_lat(1,:),
     $        1,iim,1,jjm
     $        +1, itau_phy,zjulian,dtime,nhoridbg,nidbg) 
! no vertical axis
          cl_surf(1)='ter'
          cl_surf(2)='lic'
          cl_surf(3)='oce'
          cl_surf(4)='sic'
          DO nsrf=1,nbsrf
            CALL histdef(nidbg, cl_surf(nsrf),cl_surf(nsrf), "-",iim,
     $          jjm+1,nhoridbg, 1, 1, 1, -99, 32, "inst", dtime,dtime) 
          END DO
          CALL histend(nidbg)
          CALL histsync(nidbg)
      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
         rugmer(i) = 0.0
         cdragh(i) = 0.0
         cdragm(i) = 0.0
         dflux_t(i) = 0.0
         dflux_q(i) = 0.0
         zu1(i) = 0.0
         zv1(i) = 0.0
      ENDDO
      ypct = 0.0
      yts = 0.0
      ysnow = 0.0
      yqsurf = 0.0
      yalb = 0.0
      yalblw = 0.0
      yrain_f = 0.0
      ysnow_f = 0.0
      yfder = 0.0
      ytaux = 0.0
      ytauy = 0.0
      ysolsw = 0.0
      ysollw = 0.0
      ysollwdown = 0.0
      yrugos = 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
      yq = 0.0
      pctsrf_new = 0.0
      y_flux_u = 0.0
      y_flux_v = 0.0
C$$ PB
      y_dflux_t = 0.0
      y_dflux_q = 0.0
      ytsoil = 999999.
      yrugoro = 0.
c -- LOOP
      yu10mx = 0.0
      yu10my = 0.0
      ywindsp = 0.0
c -- LOOP
      DO nsrf = 1, nbsrf
      DO i = 1, klon
         d_ts(i,nsrf) = 0.0
      ENDDO
      END DO
C§§§ PB
      yfluxlat=0.
      flux_t = 0.
      flux_q = 0.
      flux_u = 0.
      flux_v = 0.
      DO k = 1, klev
      DO i = 1, klon
         d_t(i,k) = 0.0
         d_q(i,k) = 0.0
c$$$         flux_t(i,k) = 0.0
c$$$         flux_q(i,k) = 0.0
         d_u(i,k) = 0.0
         d_v(i,k) = 0.0
c$$$         flux_u(i,k) = 0.0
c$$$         flux_v(i,k) = 0.0
         zcoefh(i,k) = 0.0
      ENDDO
      ENDDO
cAA      IF (itr.GE.1) THEN
cAA      DO it = 1, itr
cAA      DO k = 1, klev
cAA      DO i = 1, klon
cAA         d_tr(i,k,it) = 0.0
cAA      ENDDO
cAA      ENDDO
cAA      ENDDO
cAA      ENDIF

c
c Boucler sur toutes les sous-fractions du sol:
c
C Initialisation des "pourcentages potentiels". On considere ici qu'on 
C peut avoir potentiellementdela glace sur tout le domaine oceanique 
C (a affiner)

      pctsrf_pot = pctsrf
      pctsrf_pot(:,is_oce) = 1. - zmasq(:)
      pctsrf_pot(:,is_sic) = 1. - zmasq(:)

      DO 99999 nsrf = 1, nbsrf

c chercher les indices:
      DO j = 1, klon
         ni(j) = 0
      ENDDO
      knon = 0
      DO i = 1, klon

C pour determiner le domaine a traiter on utilise les surfaces "potentielles"
C  
      IF (pctsrf_pot(i,nsrf).GT.epsfra) THEN
         knon = knon + 1
         ni(knon) = i
      ENDIF
      ENDDO
c
      if (check) THEN
          write(*,*)'CLMAIN, nsrf, knon =',nsrf, knon
          call flush(6)
      endif
c
c variables pour avoir une sortie IOIPSL des INDEX
c
      IF (debugindex) THEN 
          tabindx(:)=0.
c          tabindx(1:knon)=(/FLOAT(i),i=1:knon/)
          DO i=1,knon
            tabindx(1:knon)=FLOAT(i)
          END DO 
          debugtab(:,:)=0.
          ndexbg(:)=0
          CALL gath2cpl(tabindx,debugtab,klon,knon,iim,jjm,ni)
          CALL histwrite(nidbg,cl_surf(nsrf),itap,debugtab,iim*(jjm+1)
     $        ,ndexbg)
      ENDIF 
      IF (knon.EQ.0) GOTO 99999
      DO j = 1, knon
      i = ni(j)
        ypct(j) = pctsrf(i,nsrf)
        yts(j) = ts(i,nsrf)
cIM "slab" ocean
        ytslab(i) = tslab(i)
c
        ysnow(j) = snow(i,nsrf)
        yqsurf(j) = qsurf(i,nsrf)
        yalb(j) = albe(i,nsrf)
        yalblw(j) = alblw(i,nsrf)
        yrain_f(j) = rain_f(i)
        ysnow_f(j) = snow_f(i)
        yagesno(j) = agesno(i,nsrf)
        yfder(j) = fder(i)
        ytaux(j) = flux_u(i,1,nsrf)
        ytauy(j) = flux_v(i,1,nsrf)
        ysolsw(j) = solsw(i,nsrf)
        ysollw(j) = sollw(i,nsrf)
        ysollwdown(j) = sollwdown(i)
        yrugos(j) = rugos(i,nsrf)
        yrugoro(j) = rugoro(i)
        yu1(j) = u1lay(i)
        yv1(j) = v1lay(i)
        yrads(j) =  ysolsw(j)+ ysollw(j)
        ypaprs(j,klev+1) = paprs(i,klev+1)
        y_run_off_lic_0(j) = run_off_lic_0(i)
c -- LOOP
       yu10mx(j) = u10m(i,nsrf)
       yu10my(j) = v10m(i,nsrf)
       ywindsp(j) = SQRT(yu10mx(j)*yu10mx(j) + yu10my(j)*yu10my(j) )
c -- LOOP
      END DO
C
C     IF bucket model for continent, copy soil water content
      IF ( nsrf .eq. is_ter .and. .not. ok_veget ) THEN 
          DO j = 1, knon
            i = ni(j)
            yqsol(j) = qsol(i)
          END DO
      ELSE 
          yqsol(:)=0.
      ENDIF 
c$$$ PB ajour pour soil
      DO k = 1, nsoilmx
        DO j = 1, knon
          i = ni(j)
          ytsoil(j,k) = ftsoil(i,k,nsrf)
        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)
        yq(j,k) = q(i,k)
      ENDDO
      ENDDO
c
c
c calculer Cdrag et les coefficients d'echange
      CALL coefkz(nsrf, knon, ypaprs, ypplay,
cIM 261103
     .     ksta, ksta_ter,
cIM 261103
     .            yts, yrugos, yu, yv, yt, yq,
     .            yqsurf, 
     .            ycoefm, ycoefh)
cIM 081204 BEG
cCR test
      if (iflag_pbl.eq.1) then
cIM 081204 END
        CALL coefkz2(nsrf, knon, ypaprs, ypplay,yt,
     .                  ycoefm0, ycoefh0)
        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
c
cIM cf JLD : on seuille ycoefm et ycoefh
      if (nsrf.eq.is_oce) then
         do j=1,knon
c           ycoefm(j,1)=min(ycoefm(j,1),1.1E-3)
            ycoefm(j,1)=min(ycoefm(j,1),cdmmax)
c           ycoefh(j,1)=min(ycoefh(j,1),1.1E-3)
            ycoefh(j,1)=min(ycoefh(j,1),cdhmax)
         enddo
      endif

c
cIM: 261103
      if (ok_kzmin) THEN
cIM cf FH: 201103 BEG
c   Calcul d'une diffusion minimale pour les conditions tres stables.
      call coefkzmin(knon,ypaprs,ypplay,yu,yv,yt,yq,ycoefm
     .   ,ycoefm0,ycoefh0)
c      call dump2d(iim,jjm-1,ycoefm(2:klon-1,2), 'KZ         ')
c      call dump2d(iim,jjm-1,ycoefm0(2:klon-1,2),'KZMIN      ')
 
       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

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c MELLOR ET YAMADA adapte a Mars Richard Fournier et Frederic Hourdin
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

         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
         do k=1,klev
            yteta(1:knon,k)=
     .      yt(1:knon,k)*(ypaprs(1:knon,1)/ypplay(1:knon,k))**rkappa
     .      *(1.+0.61*yq(1:knon,k))
         enddo
         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
         DO k = 1, klev+1
            DO j = 1, knon
               i = ni(j)
               yq2(j,k)=q2(i,k,nsrf)
            enddo
         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 longuer 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,yq2,q2diag,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,yq2,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)


      ENDIF

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
      CALL clqh(dtime, itap, date0,jour, debut,lafin,
     e          rlon, rlat, cufi, cvfi,
     e          knon, nsrf, ni, pctsrf,
     e          soil_model, ytsoil,yqsol,
     e          ok_veget, ocean, npas, nexca,
     e          rmu0, co2_ppm, yrugos, yrugoro,
     e          yu1, yv1, ycoefh,
     e          yt,yq,yts,ypaprs,ypplay,
     e          ydelp,yrads,yalb, yalblw, ysnow, yqsurf, 
     e          yrain_f, ysnow_f, yfder, ytaux, ytauy,
c -- LOOP
     e          ywindsp, 
c -- LOOP
c$$$     e          ysollw, ysolsw,
     e          ysollw, ysollwdown, ysolsw,yfluxlat,
     s          pctsrf_new, yagesno,
     s          y_d_t, y_d_q, y_d_ts, yz0_new,
     s          y_flux_t, y_flux_q, y_dflux_t, y_dflux_q,
     s          y_fqcalving,y_ffonte,y_run_off_lic_0,
cIM "slab" ocean
     s          y_flux_o, y_flux_g, ytslab, y_seaice)
c
c calculer la longueur de rugosite sur ocean
      yrugm=0.
      IF (nsrf.EQ.is_oce) THEN
      DO j = 1, knon
         yrugm(j) = 0.018*ycoefm(j,1) * (yu1(j)**2+yv1(j)**2)/RG 
     $      +  0.11*14e-6 / sqrt(ycoefm(j,1) * (yu1(j)**2+yv1(j)**2))
         yrugm(j) = MAX(1.5e-05,yrugm(j))
      ENDDO
      ENDIF
      DO j = 1, knon
         y_dflux_t(j) = y_dflux_t(j) * ypct(j)
         y_dflux_q(j) = y_dflux_q(j) * ypct(j)
         yu1(j) = yu1(j) *  ypct(j)
         yv1(j) = yv1(j) *  ypct(j)
      ENDDO
c
      DO k = 1, klev
        DO j = 1, knon
          i = ni(j)
          ycoefh(j,k) = ycoefh(j,k) * ypct(j)
          ycoefm(j,k) = ycoefm(j,k) * ypct(j)
          y_d_t(j,k) = y_d_t(j,k) * ypct(j)
          y_d_q(j,k) = y_d_q(j,k) * ypct(j)
C§§§ PB
          flux_t(i,k,nsrf) = y_flux_t(j,k)
          flux_q(i,k,nsrf) = y_flux_q(j,k)
          flux_u(i,k,nsrf) = y_flux_u(j,k)
          flux_v(i,k,nsrf) = y_flux_v(j,k)
c$$$ PB        y_flux_t(j,k) = y_flux_t(j,k) * ypct(j)
c$$$ PB        y_flux_q(j,k) = y_flux_q(j,k) * ypct(j)
          y_d_u(j,k) = y_d_u(j,k) * ypct(j)
          y_d_v(j,k) = y_d_v(j,k) * ypct(j)
c$$$ PB        y_flux_u(j,k) = y_flux_u(j,k) * ypct(j)
c$$$ PB        y_flux_v(j,k) = y_flux_v(j,k) * ypct(j)
        ENDDO
      ENDDO


      evap(:,nsrf) = - flux_q(:,1,nsrf)
c
      albe(:, nsrf) = 0.
      alblw(:, nsrf) = 0.
      snow(:, nsrf) = 0.
      qsurf(:, nsrf) = 0.
      rugos(:, nsrf) = 0.
      fluxlat(:,nsrf) = 0.
      DO j = 1, knon
         i = ni(j)
         d_ts(i,nsrf) = y_d_ts(j)
         albe(i,nsrf) = yalb(j)
         alblw(i,nsrf) = yalblw(j)
         snow(i,nsrf) = ysnow(j)
         qsurf(i,nsrf) = yqsurf(j)
         rugos(i,nsrf) = yz0_new(j)
         fluxlat(i,nsrf) = yfluxlat(j)
c$$$ pb         rugmer(i) = yrugm(j)
         IF (nsrf .EQ. is_oce) then 
           rugmer(i) = yrugm(j)
           rugos(i,nsrf) = yrugm(j)
         endif	
cIM cf JLD ??
         fqcalving(i,nsrf) = y_fqcalving(j)        
         ffonte(i,nsrf) = y_ffonte(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)
         dflux_q(i) = dflux_q(i) + y_dflux_q(j)
         zu1(i) = zu1(i) + yu1(j)
         zv1(i) = zv1(i) + yv1(j)
      END DO
      IF ( nsrf .eq. is_ter ) THEN 
      DO j = 1, knon
         i = ni(j)
         qsol(i) = yqsol(j)
      END DO
      END IF 
      IF ( nsrf .eq. is_lic ) THEN 
        DO j = 1, knon
          i = ni(j)
          run_off_lic_0(i) = y_run_off_lic_0(j)
        END DO
      END IF 
c$$$ PB ajout pour soil
      ftsoil(:,:,nsrf) = 0.
      DO k = 1, nsoilmx
        DO j = 1, knon
          i = ni(j)
          ftsoil(i, k, nsrf) = ytsoil(j,k)
        END DO 
      END DO 
c
#ifdef CRAY
      DO k = 1, klev
      DO j = 1, knon
      i = ni(j)
#else
      DO j = 1, knon
      i = ni(j)
      DO k = 1, klev
#endif
         d_t(i,k) = d_t(i,k) + y_d_t(j,k)
         d_q(i,k) = d_q(i,k) + y_d_q(j,k)
c$$$ PB        flux_t(i,k) = flux_t(i,k) + y_flux_t(j,k)
c$$$         flux_q(i,k) = flux_q(i,k) + y_flux_q(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)
c$$$  PB       flux_u(i,k) = flux_u(i,k) + y_flux_u(j,k)
c$$$         flux_v(i,k) = flux_v(i,k) + y_flux_v(j,k)
         zcoefh(i,k) = zcoefh(i,k) + ycoefh(j,k)
      ENDDO
      ENDDO
c
c
#undef T2m     
#define T2m     
#ifdef T2m
ccc diagnostic t,q a 2m et u, v a 10m
c
      DO j=1, knon
        i = ni(j)
        uzon(j) = yu(j,1) + y_d_u(j,1)
        vmer(j) = yv(j,1) + y_d_v(j,1)
        tair1(j) = yt(j,1) + y_d_t(j,1)
        qair1(j) = yq(j,1) + y_d_q(j,1)
        zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1)))
     &                   * (ypaprs(j,1)-ypplay(j,1))
        tairsol(j) = yts(j) + y_d_ts(j)
        rugo1(j) = yrugos(j)
        IF(nsrf.EQ.is_oce) THEN
         rugo1(j) = rugos(i,nsrf)
        ENDIF
        psfce(j)=ypaprs(j,1)
        patm(j)=ypplay(j,1)
c
        qairsol(j) = yqsurf(j)
      ENDDO
c
      CALL stdlevvar(klon, knon, nsrf, zxli,
     &               uzon, vmer, tair1, qair1, zgeo1,
     &               tairsol, qairsol, rugo1, psfce, patm,
cIM  &               yt2m, yq2m, yu10m)
     &               yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
cIM 081204 END
c
c
      DO j=1, knon
       i = ni(j)
       t2m(i,nsrf)=yt2m(j)

c
       q2m(i,nsrf)=yq2m(j)
c
c u10m, v10m : composantes du vent a 10m sans spirale de Ekman
       u10m(i,nsrf)=(yu10m(j) * uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)
       v10m(i,nsrf)=(yu10m(j) * vmer(j))/sqrt(uzon(j)**2+vmer(j)**2)
c
      ENDDO
c
cIM cf AM : pbl, HBTM
      DO i = 1, knon
         y_cd_h(i) = ycoefh(i,1)
         y_cd_m(i) = ycoefm(i,1)
      ENDDO
c     print*,'appel hbtm2'
      CALL HBTM(knon, ypaprs, ypplay,
     .          yt2m,yt10m,yq2m,yq10m,yustar,
     .          y_flux_t,y_flux_q,yu,yv,yt,yq,
     .          ypblh,ycapCL,yoliqCL,ycteiCL,ypblT,
     .          ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl)
c     print*,'fin hbtm2'
c
      DO j=1, knon
       i = ni(j)
       pblh(i,nsrf)   = ypblh(j)
       plcl(i,nsrf)   = ylcl(j)
       capCL(i,nsrf)  = ycapCL(j)
       oliqCL(i,nsrf) = yoliqCL(j)
       cteiCL(i,nsrf) = ycteiCL(j)
       pblT(i,nsrf)   = ypblT(j)
       therm(i,nsrf)  = ytherm(j)
       trmb1(i,nsrf)  = ytrmb1(j)
       trmb2(i,nsrf)  = ytrmb2(j)
       trmb3(i,nsrf)  = ytrmb3(j)
      ENDDO
c
#else
cIM cf AM : pbl, HBTM
      DO j=1, knon
       i = ni(j)
       pblh(i,nsrf)   = 0.
       plcl(i,nsrf)   = 0.
       capCL(i,nsrf)  = 0.
       oliqCL(i,nsrf) = 0.
       cteiCL(i,nsrf) = 0.
       pblT(i,nsrf)   = 0.
       therm(i,nsrf)  = 0.
       trmb1(i,nsrf)  = 0.
       trmb2(i,nsrf)  = 0.
       trmb3(i,nsrf)  = 0.
      ENDDO 
      DO j=1, knon
       i = ni(j) 
       t2m(i,nsrf)=0.
       q2m(i,nsrf)=0.
       u10m(i,nsrf)=0.
       v10m(i,nsrf)=0.
      ENDDO 
#endif

      do j=1,knon
         do k=1,klev+1
         i=ni(j)
         q2(i,k,nsrf)=yq2(j,k)
         enddo
      enddo
cIM "slab" ocean 
      IF(OCEAN.EQ.'slab  '.OR.OCEAN.EQ.'force ') THEN
       IF (nsrf.EQ.is_oce) THEN
        DO j = 1, knon
c on projette sur la grille globale
         i = ni(j)
         IF(pctsrf_new(i,is_oce).GT.epsfra) THEN
          flux_o(i) = y_flux_o(j)
         ELSE
          flux_o(i) = 0.
         ENDIF
        ENDDO
       ENDIF
c
       IF (nsrf.EQ.is_sic) THEN
        DO j = 1, knon
         i = ni(j)
cIM 230604 on pondere lorsque l'on fait le bilan au sol :  flux_g(i) = y_flux_g(j)*ypct(j)
         IF(pctsrf_new(i,is_sic).GT.epsfra) THEN
          flux_g(i) = y_flux_g(j)
         ELSE
          flux_g(i) = 0.
         ENDIF
        ENDDO
       ENDIF !nsrf.EQ.is_sic
      ENDIF !OCEAN
c
      IF(OCEAN.EQ.'slab  ') THEN
       IF(nsrf.EQ.is_oce) then
        tslab(1:klon) = ytslab(1:klon)
        seaice(1:klon) = y_seaice(1:klon)
       ENDIF !nsrf
      ENDIF !OCEAN
99999 CONTINUE
C
C On utilise les nouvelles surfaces
C A rajouter: conservation de l'albedo
C
      rugos(:,is_oce) = rugmer
      pctsrf = pctsrf_new

      RETURN
      END
      SUBROUTINE clqh(dtime,itime, date0,jour,debut,lafin,
     e                rlon, rlat, cufi, cvfi, 
     e                knon, nisurf, knindex, pctsrf,
     $                soil_model,tsoil,qsol,
     e                ok_veget, ocean, npas, nexca,
     e                rmu0, co2_ppm, rugos, rugoro,
     e                u1lay,v1lay,coef,
     e                t,q,ts,paprs,pplay,
     e                delp,radsol,albedo,alblw,snow,qsurf, 
     e                precip_rain, precip_snow, fder, taux, tauy,
c -- LOOP
     e                ywindsp,
c -- LOOP
     $                sollw, sollwdown, swnet,fluxlat, 
     s                pctsrf_new, agesno,
     s                d_t, d_q, d_ts, z0_new, 
     s                flux_t, flux_q,dflux_s,dflux_l,
     s                fqcalving,ffonte,run_off_lic_0,
cIM "slab" ocean
     s                flux_o,flux_g,tslab,seaice)

      USE interface_surf

      IMPLICIT none
c======================================================================
c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
c Objet: diffusion verticale de "q" et de "h"
c======================================================================
#include "dimensions.h"
#include "dimphy.h"
#include "YOMCST.h"
#include "YOETHF.h"
#include "FCTTRE.h"
#include "indicesol.h"
#include "dimsoil.h"
#include "iniprint.h"

c Arguments:
      INTEGER knon
      REAL dtime              ! intervalle du temps (s)
      real date0
      REAL u1lay(klon)        ! vitesse u de la 1ere couche (m/s)
      REAL v1lay(klon)        ! vitesse v de la 1ere couche (m/s)
      REAL 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 t(klon,klev)       ! temperature (K)
      REAL q(klon,klev)       ! humidite specifique (kg/kg)
      REAL ts(klon)           ! temperature du sol (K)
      REAL evap(klon)         ! evaporation au sol
      REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)
      REAL pplay(klon,klev)   ! pression au milieu de couche (Pa)
      REAL delp(klon,klev)    ! epaisseur de couche en pression (Pa)
      REAL radsol(klon)       ! ray. net au sol (Solaire+IR) W/m2
      REAL albedo(klon)       ! albedo de la surface
      REAL alblw(klon)
      REAL snow(klon)         ! hauteur de neige
      REAL qsurf(klon)         ! humidite de l'air au dessus de la surface
      real precip_rain(klon), precip_snow(klon)
      REAL agesno(klon)
      REAL rugoro(klon)
      REAL run_off_lic_0(klon)! runof glacier au pas de temps precedent
      integer jour            ! jour de l'annee en cours
      real rmu0(klon)         ! cosinus de l'angle solaire zenithal
      real rugos(klon)        ! rugosite
      integer knindex(klon)
      real pctsrf(klon,nbsrf)
      real rlon(klon), rlat(klon), cufi(klon), cvfi(klon)
      logical ok_veget 
      REAL co2_ppm            ! taux CO2 atmosphere
      character*6 ocean
      integer npas, nexca
c -- LOOP
       REAL yu10mx(klon)
       REAL yu10my(klon)
       REAL ywindsp(klon)
c -- LOOP


c
      REAL d_t(klon,klev)     ! incrementation de "t"
      REAL d_q(klon,klev)     ! incrementation de "q"
      REAL d_ts(klon)         ! incrementation de "ts"
      REAL 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 flux_q(klon,klev)  ! flux de la vapeur d'eau:kg/(m**2 s)
      REAL dflux_s(klon) ! derivee du flux sensible dF/dTs
      REAL dflux_l(klon) ! derivee du flux latent dF/dTs
cIM cf JLD
c Flux thermique utiliser pour fondre la neige
      REAL ffonte(klon)
c Flux d'eau "perdue" par la surface et nécessaire pour que limiter la
c hauteur de neige, en kg/m2/s
      REAL fqcalving(klon)
cIM "slab" ocean
      REAL tslab(klon)  !temperature du slab ocean (K) (OCEAN='slab  ')
      REAL seaice(klon) ! glace de mer en kg/m2
      REAL flux_o(klon) ! flux entre l'ocean et l'atmosphere W/m2
      REAL flux_g(klon) ! flux entre l'ocean et la glace de mer W/m2
c
c======================================================================
      REAL t_grnd  ! temperature de rappel pour glace de mer
      PARAMETER (t_grnd=271.35)
      REAL t_coup
      PARAMETER(t_coup=273.15)
c======================================================================
      INTEGER i, k
      REAL zx_cq(klon,klev)
      REAL zx_dq(klon,klev)
      REAL zx_ch(klon,klev)
      REAL zx_dh(klon,klev)
      REAL zx_buf1(klon)
      REAL zx_buf2(klon)
      REAL zx_coef(klon,klev)
      REAL local_h(klon,klev) ! enthalpie potentielle
      REAL local_q(klon,klev)
      REAL local_ts(klon)
      REAL psref(klon) ! pression de reference pour temperature potent.
      REAL zx_pkh(klon,klev), zx_pkf(klon,klev)
c======================================================================
c contre-gradient pour la vapeur d'eau: (kg/kg)/metre
      REAL gamq(klon,2:klev)
c contre-gradient pour la chaleur sensible: Kelvin/metre
      REAL gamt(klon,2:klev)
      REAL z_gamaq(klon,2:klev), z_gamah(klon,2:klev)
      REAL zdelz
c======================================================================
#include "compbl.h"
c======================================================================
c Rajout pour l'interface
      integer itime
      integer nisurf
      logical debut, lafin
      real zlev1(klon)
      real fder(klon), taux(klon), tauy(klon)
      real temp_air(klon), spechum(klon)
      real epot_air(klon), ccanopy(klon)
      real tq_cdrag(klon), petAcoef(klon), peqAcoef(klon)
      real petBcoef(klon), peqBcoef(klon)
      real sollw(klon), sollwdown(klon), swnet(klon), swdown(klon)
      real p1lay(klon)
c$$$C PB ajout pour soil
      LOGICAL soil_model
      REAL tsoil(klon, nsoilmx)
      REAL qsol(klon)

! Parametres de sortie
      real fluxsens(klon), fluxlat(klon)
      real tsol_rad(klon), tsurf_new(klon), alb_new(klon)
      real emis_new(klon), z0_new(klon)
      real pctsrf_new(klon,nbsrf)
c JLD
      real zzpk
C
      character (len = 20) :: modname = 'Debut clqh'
      LOGICAL check
      PARAMETER (check=.false.)
C
      if (check) THEN
          write(*,*) modname,' nisurf=',nisurf
          call flush(6)
      endif
c
      if (check) THEN
       WRITE(*,*)' qsurf (min, max)'
     $     , minval(qsurf(1:knon)), maxval(qsurf(1:knon))
       call flush(6)
      ENDIF
C
C
      if (iflag_pbl.eq.1) then
        do k = 3, klev
          do i = 1, knon
            gamq(i,k)= 0.0
            gamt(i,k)=  -1.0e-03
          enddo
        enddo
        do i = 1, knon
          gamq(i,2) = 0.0
          gamt(i,2) = -2.5e-03
        enddo
      else
        do k = 2, klev
          do i = 1, knon
            gamq(i,k) = 0.0
            gamt(i,k) = 0.0
          enddo
        enddo
      endif

      DO i = 1, knon
         psref(i) = paprs(i,1) !pression de reference est celle au sol
         local_ts(i) = ts(i)
      ENDDO
      DO k = 1, klev
      DO i = 1, knon
         zx_pkh(i,k) = (psref(i)/paprs(i,k))**RKAPPA
         zx_pkf(i,k) = (psref(i)/pplay(i,k))**RKAPPA
         local_h(i,k) = RCPD * t(i,k) * zx_pkf(i,k)
         local_q(i,k) = q(i,k)
      ENDDO
      ENDDO
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))
     .                  *(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
c Preparer les flux lies aux contre-gardients
c
      DO k = 2, klev
      DO i = 1, knon
         zdelz = RD * (t(i,k-1)+t(i,k))/2.0 / RG /paprs(i,k)
     .              *(pplay(i,k-1)-pplay(i,k))
         z_gamaq(i,k) = gamq(i,k) * zdelz
         z_gamah(i,k) = gamt(i,k) * zdelz *RCPD * zx_pkh(i,k)
      ENDDO
      ENDDO
      DO i = 1, knon
         zx_buf1(i) = zx_coef(i,klev) + delp(i,klev)
         zx_cq(i,klev) = (local_q(i,klev)*delp(i,klev)
     .                   -zx_coef(i,klev)*z_gamaq(i,klev))/zx_buf1(i)
         zx_dq(i,klev) = zx_coef(i,klev) / zx_buf1(i)
c
         zzpk=(pplay(i,klev)/psref(i))**RKAPPA
         zx_buf2(i) = zzpk*delp(i,klev) + zx_coef(i,klev)
         zx_ch(i,klev) = (local_h(i,klev)*zzpk*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
         zx_buf1(i) = delp(i,k)+zx_coef(i,k)
     .               +zx_coef(i,k+1)*(1.-zx_dq(i,k+1))
         zx_cq(i,k) = (local_q(i,k)*delp(i,k)
     .                 +zx_coef(i,k+1)*zx_cq(i,k+1)
     .                 +zx_coef(i,k+1)*z_gamaq(i,k+1)
     .                 -zx_coef(i,k)*z_gamaq(i,k))/zx_buf1(i)
         zx_dq(i,k) = zx_coef(i,k) / zx_buf1(i)
c
         zzpk=(pplay(i,k)/psref(i))**RKAPPA
         zx_buf2(i) = zzpk*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)*zzpk*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 q1 = zx_cq(i,1) + zx_dq(i,1) * Flux_Q(i,1) * dt
C h1 = zx_ch(i,1) + zx_dh(i,1) * Flux_H(i,1) * dt
C
      DO i = 1, knon
         zx_buf1(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dq(i,2))
         zx_cq(i,1) = (local_q(i,1)*delp(i,1)
     .                 +zx_coef(i,2)*(z_gamaq(i,2)+zx_cq(i,2)))
     .                /zx_buf1(i)
         zx_dq(i,1) = -1. * RG / zx_buf1(i)
c
         zzpk=(pplay(i,1)/psref(i))**RKAPPA
         zx_buf2(i) = zzpk*delp(i,1) + zx_coef(i,2)*(1.-zx_dh(i,2))
         zx_ch(i,1) = (local_h(i,1)*zzpk*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. 
	peqAcoef = 0.
	petBcoef =0.
	peqBcoef = 0.
	p1lay =0.
	
c      do i = 1, knon
        petAcoef(1:knon) = zx_ch(1:knon,1)
        peqAcoef(1:knon) = zx_cq(1:knon,1)
        petBcoef(1:knon) =  zx_dh(1:knon,1)
        peqBcoef(1:knon) = zx_dq(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)
        spechum(1:knon)=q(1:knon,1)
        p1lay(1:knon) = pplay(1:knon,1)
        zlev1(1:knon) = delp(1:knon,1)
c        swnet = swdown * (1. - albedo)
c
cIM swdown=flux SW incident sur terres
cIM swdown=flux SW net sur les autres surfaces 
cIM     swdown(1:knon) = swnet(1:knon)
        if(nisurf.eq.is_ter) THEN 
         swdown(1:knon) = swnet(1:knon)/(1-albedo(1:knon)) 
        else 
         swdown(1:knon) = swnet(1:knon)
        endif
c      enddo
      ccanopy = co2_ppm

      CALL interfsurf(itime, dtime, date0, jour, rmu0,
     e klon, iim, jjm, nisurf, knon, knindex, pctsrf, 
     e rlon, rlat, cufi, cvfi, 
     e debut, lafin, ok_veget, soil_model, nsoilmx,tsoil, qsol,
     e zlev1,  u1lay, v1lay, temp_air, spechum, epot_air, ccanopy, 
     e tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef,
     e precip_rain, precip_snow, sollw, sollwdown, swnet, swdown,
     e fder, taux, tauy, 
c -- LOOP
     e ywindsp,
c -- LOOP
     e rugos, rugoro,
     e albedo, snow, qsurf,
     e ts, p1lay, psref, radsol,
     e ocean, npas, nexca, zmasq,
     s evap, fluxsens, fluxlat, dflux_l, dflux_s,              
     s tsol_rad, tsurf_new, alb_new, alblw, emis_new, z0_new, 
     s pctsrf_new, agesno,fqcalving,ffonte, run_off_lic_0,
cIM "slab" ocean 
     s flux_o, flux_g, tslab, seaice)


      do i = 1, knon
        flux_t(i,1) = fluxsens(i)
        flux_q(i,1) = - evap(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
         local_q(i,1) = zx_cq(i,1) + zx_dq(i,1)*flux_q(i,1)*dtime
      ENDDO
      DO k = 2, klev
      DO i = 1, knon
        local_q(i,k) = zx_cq(i,k) + zx_dq(i,k)*local_q(i,k-1)
        local_h(i,k) = zx_ch(i,k) + zx_dh(i,k)*local_h(i,k-1)
      ENDDO
      ENDDO
c======================================================================
c== flux_q est le flux de vapeur d'eau: kg/(m**2 s)  positive vers bas
c== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
      DO k = 2, klev
      DO i = 1, knon
        flux_q(i,k) = (zx_coef(i,k)/RG/dtime)
     .                * (local_q(i,k)-local_q(i,k-1)+z_gamaq(i,k))
        flux_t(i,k) = (zx_coef(i,k)/RG/dtime)
     .                * (local_h(i,k)-local_h(i,k-1)+z_gamah(i,k))
     .                / zx_pkh(i,k)
      ENDDO
      ENDDO
c======================================================================
C Calcul tendances
      DO k = 1, klev
      DO i = 1, knon
         d_t(i,k) = local_h(i,k)/zx_pkf(i,k)/RCPD - t(i,k)
         d_q(i,k) = local_q(i,k) - q(i,k)
      ENDDO
      ENDDO
c

      RETURN
      END
      SUBROUTINE clvent(knon,dtime, u1lay,v1lay,coef,t,ven,
     e                  paprs,pplay,delp,
     s                  d_ven,flux_v)
      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 "dimensions.h"
#include "dimphy.h"
#include "iniprint.h"
      INTEGER knon
      REAL dtime
      REAL u1lay(klon), v1lay(klon)
      REAL coef(klon,klev)
      REAL t(klon,klev), ven(klon,klev)
      REAL paprs(klon,klev+1), pplay(klon,klev), delp(klon,klev)
      REAL d_ven(klon,klev)
      REAL flux_v(klon,klev)
c======================================================================
#include "YOMCST.h"
c======================================================================
      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)
     .                 * (1.0+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
      RETURN
      END
      SUBROUTINE coefkz(nsrf, knon, paprs, pplay,
cIM 261103
     .                  ksta, ksta_ter,
cIM 261103
     .                  ts, rugos,
     .                  u,v,t,q,
     .                  qsurf, 
     .                  pcfm, pcfh)
      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 nsrf-----input-I- indicateur de la nature du sol
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 rugos----input-R- longeur de rugosite (en m)
c u--------input-R- vitesse u
c v--------input-R- vitesse v
c t--------input-R- temperature (K)
c q--------input-R- vapeur d'eau (kg/kg)
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 "dimensions.h"
#include "dimphy.h"
#include "YOMCST.h"
#include "indicesol.h"
#include "iniprint.h"
c
c Arguments:
c
      INTEGER knon, nsrf
      REAL ts(klon)
      REAL paprs(klon,klev+1), pplay(klon,klev)
      REAL u(klon,klev), v(klon,klev), t(klon,klev), q(klon,klev)
      REAL rugos(klon)
c
      REAL pcfm(klon,klev), pcfh(klon,klev)
      INTEGER itop(klon)
c
c Quelques constantes et options:
c
      REAL cepdu2, ckap, cb, cc, cd, clam
      PARAMETER (cepdu2 =(0.1)**2)
      PARAMETER (CKAP=0.4)
      PARAMETER (cb=5.0)
      PARAMETER (cc=5.0)
      PARAMETER (cd=5.0)
      PARAMETER (clam=160.0)
      REAL ratqs ! largeur de distribution de vapeur d'eau
      PARAMETER (ratqs=0.05)
      LOGICAL richum ! utilise le nombre de Richardson humide
      PARAMETER (richum=.TRUE.)
      REAL ric ! nombre de Richardson critique
      PARAMETER(ric=0.4)
      REAL prandtl
      PARAMETER (prandtl=0.4)
      REAL kstable ! diffusion minimale (situation stable)
      ! GKtest
      ! PARAMETER (kstable=1.0e-10)
      REAL ksta, ksta_ter
cIM: 261103     REAL kstable_ter, kstable_sinon
cIM: 211003 cf GK   PARAMETER (kstable_ter = 1.0e-6)
cIM: 261103     PARAMETER (kstable_ter = 1.0e-8)
cIM: 261103   PARAMETER (kstable_ter = 1.0e-10)
cIM: 261103   PARAMETER (kstable_sinon = 1.0e-10)
      ! fin GKtest
      REAL mixlen ! constante controlant longueur de melange
      PARAMETER (mixlen=35.0)
      INTEGER isommet ! le sommet de la couche limite
      PARAMETER (isommet=klev)
      LOGICAL tvirtu ! calculer Ri d'une maniere plus performante
      PARAMETER (tvirtu=.TRUE.)
      LOGICAL opt_ec ! formule du Centre Europeen dans l'atmosphere
      PARAMETER (opt_ec=.FALSE.)

#include "compbl.h"
c
c Variables locales:
c
      INTEGER i, k, kk !IM 120704
      REAL zgeop(klon,klev)
      REAL zmgeom(klon)
      REAL zri(klon)
      REAL zl2(klon)

      REAL u1(klon), v1(klon), t1(klon), q1(klon), z1(klon)
      REAL pcfm1(klon), pcfh1(klon)
c
      REAL zdphi, zdu2, ztvd, ztvu, zcdn
      REAL zscf
      REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs
      REAL z2geomf, zalh2, zalm2, zscfh, zscfm
      REAL t_coup
      PARAMETER (t_coup=273.15)
cIM
      LOGICAL check
      PARAMETER (check=.false.)
c
c contre-gradient pour la chaleur sensible: Kelvin/metre
      REAL gamt(2:klev)
      real qsurf(klon) 
c
      LOGICAL appel1er
      SAVE appel1er
c
c Fonctions thermodynamiques et fonctions d'instabilite
      REAL fsta, fins, x
      LOGICAL zxli ! utiliser un jeu de fonctions simples
      PARAMETER (zxli=.FALSE.)
c
#include "YOETHF.h"
#include "FCTTRE.h"
      fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))
      fins(x) = SQRT(1.0-18.0*x)
c
      DATA appel1er /.TRUE./
c
      IF (appel1er) THEN
        if (prt_level > 9) THEN
          WRITE(lunout,*)'coefkz, opt_ec:', opt_ec
          WRITE(lunout,*)'coefkz, richum:', richum
          IF (richum) WRITE(lunout,*)'coefkz, ratqs:', ratqs
          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
cIM cf JLD/ GKtest
      IF ( nsrf .NE. is_oce ) THEN
cIM 261103     kstable = kstable_ter
        kstable = ksta_ter
      ELSE
cIM 261103     kstable = kstable_sinon
        kstable = ksta
      ENDIF
cIM cf JLD/ GKtest fin
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 le frottement au sol (Cdrag)
c
      DO i = 1, knon
       u1(i) = u(i,1)
       v1(i) = v(i,1)
       t1(i) = t(i,1)
       q1(i) = q(i,1)
       z1(i) = zgeop(i,1)
      ENDDO
c
      CALL clcdrag(klon, knon, nsrf, zxli, 
     $             u1, v1, t1, q1, z1,
     $             ts, qsurf, rugos,
     $             pcfm1, pcfh1) 
cIM  $             ts, qsurf, rugos,
C
      DO i = 1, knon
       pcfm(i,1)=pcfm1(i)
       pcfh(i,1)=pcfh1(i)
      ENDDO
c
c Calculer les coefficients turbulents dans l'atmosphere
c
      DO i = 1, knon
         itop(i) = isommet
      ENDDO


      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)
            zmgeom(i)=zgeop(i,k)-zgeop(i,k-1)
            zdphi =zmgeom(i) / 2.0
            zt = (t(i,k)+t(i,k-1)) * 0.5
            zq = (q(i,k)+q(i,k-1)) * 0.5
c
c           calculer Qs et dQs/dT:
c
            IF (thermcep) THEN
              zdelta = MAX(0.,SIGN(1.,RTT-zt))
              zcvm5 = R5LES*RLVTT/RCPD/(1.0+RVTMP2*zq)*(1.-zdelta) 
     .            + R5IES*RLSTT/RCPD/(1.0+RVTMP2*zq)*zdelta
              zqs = R2ES * FOEEW(zt,zdelta) / pplay(i,k)
              zqs = MIN(0.5,zqs)
              zcor = 1./(1.-RETV*zqs)
              zqs = zqs*zcor
              zdqs = FOEDE(zt,zdelta,zcvm5,zqs,zcor)
            ELSE
              IF (zt .LT. t_coup) THEN
                 zqs = qsats(zt) / pplay(i,k)
                 zdqs = dqsats(zt,zqs)
              ELSE
                 zqs = qsatl(zt) / pplay(i,k)
                 zdqs = dqsatl(zt,zqs)
              ENDIF
            ENDIF
c
c           calculer la fraction nuageuse (processus humide):
c
            zfr = (zq+ratqs*zq-zqs) / (2.0*ratqs*zq)
            zfr = MAX(0.0,MIN(1.0,zfr))
            IF (.NOT.richum) zfr = 0.0
c
c           calculer le nombre de Richardson:
c
            IF (tvirtu) THEN
            ztvd =( t(i,k)
     .             + zdphi/RCPD/(1.+RVTMP2*zq)
     .              *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) )
     .            )*(1.+RETV*q(i,k))
            ztvu =( t(i,k-1)
     .             - zdphi/RCPD/(1.+RVTMP2*zq)
     .              *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) )
     .            )*(1.+RETV*q(i,k-1))
            zri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))
            zri(i) = zri(i)
     .             + zmgeom(i)*zmgeom(i)/RG*gamt(k)
     .               *(paprs(i,k)/101325.0)**RKAPPA
     .               /(zdu2*0.5*(ztvd+ztvu))
c
            ELSE ! calcul de Ridchardson compatible LMD5
c
            zri(i) =(RCPD*(t(i,k)-t(i,k-1))
     .              -RD*0.5*(t(i,k)+t(i,k-1))/paprs(i,k)
     .               *(pplay(i,k)-pplay(i,k-1))
     .              )*zmgeom(i)/(zdu2*0.5*RCPD*(t(i,k-1)+t(i,k)))
            zri(i) = zri(i) +
     .             zmgeom(i)*zmgeom(i)*gamt(k)/RG
cSB     .             /(paprs(i,k)/101325.0)**RKAPPA
     .             *(paprs(i,k)/101325.0)**RKAPPA
     .             /(zdu2*0.5*(t(i,k-1)+t(i,k)))
            ENDIF
c
c           finalement, les coefficients d'echange sont obtenus:
c
            zcdn=SQRT(zdu2) / zmgeom(i) * 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).LT.0.0) THEN  ! situation instable
               zscf = ((zgeop(i,k)/zgeop(i,k-1))**(1./3.)-1.)**3
     .                / (zmgeom(i)/RG)**3 / (zgeop(i,k-1)/RG)
               zscf = SQRT(-zri(i)*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)*zscfm)
               pcfh(i,k)=zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh)
            ELSE ! situation stable
               zscf=SQRT(1.+cd*zri(i))
               pcfm(i,k)=zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf)
               pcfh(i,k)=zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf)
            ENDIF
          ELSE
            zl2(i)=(mixlen*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))/ric, kstable))
            pcfm(i,k)= zl2(i)* pcfm(i,k)
            pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
          ENDIF
      ENDDO
      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
      RETURN
      END

      SUBROUTINE coefkz2(nsrf, knon, paprs, pplay,t,
     .                  pcfm, pcfh)
      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 nsrf-----input-I- indicateur de la nature du sol
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 "dimensions.h"
#include "dimphy.h"
#include "YOMCST.h"
#include "indicesol.h"
#include "iniprint.h"
c
c Arguments:
c
      INTEGER knon, nsrf
      REAL paprs(klon,klev+1), pplay(klon,klev)
      REAL t(klon,klev)
c
      REAL pcfm(klon,klev), pcfh(klon,klev)
c
c Quelques constantes et options:
c
      REAL prandtl
      PARAMETER (prandtl=0.4)
      REAL kstable
      PARAMETER (kstable=0.002)
ccc      PARAMETER (kstable=0.001)
      REAL mixlen ! constante controlant longueur de melange
      PARAMETER (mixlen=35.0)
      REAL seuil ! au-dela l'inversion est consideree trop faible
      PARAMETER (seuil=-0.02)
ccc      PARAMETER (seuil=-0.04)
ccc      PARAMETER (seuil=-0.06)
ccc      PARAMETER (seuil=-0.09)
c
c Variables locales:
c
      INTEGER i, k, invb(knon)
      REAL zl2(knon)
      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
         zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1))
     .          - RD * 0.5*(t(i,k)+t(i,k+1))/RCPD/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
c Introduire une diffusion:
c
      DO k = 2, klev
      DO i = 1, knon
cIM cf FH/GK   IF ( (nsrf.NE.is_oce) .OR.  ! si ce n'est pas sur l'ocean
cIM cf FH/GK  .     (invb(i).EQ.klev) .OR. ! s'il n'y a pas d'inversion
      !IM cf JLD/ GKtest TERkz2
      ! IF ( (nsrf.EQ.is_ter) .OR.  ! si on est sur la terre
      ! fin GKtest
      IF ( (nsrf.EQ.is_oce) .AND.  ! si on est sur ocean et si 
     .     ( (invb(i).EQ.klev) .OR.      ! s'il n'y a pas d'inversion
     .     (zdthmin(i).GT.seuil) ) )THEN ! si l'inversion est trop faible
         zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,klev+1))
     .                       /(paprs(i,2)-paprs(i,klev+1)) ))**2
         pcfm(i,k)= zl2(i)* kstable
         pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
      ENDIF
      ENDDO
      ENDDO
c
      RETURN
      END
      SUBROUTINE calbeta(dtime,indice,knon,snow,qsol,
     .                    vbeta,vcal,vdif)
      IMPLICIT none
c======================================================================
c Auteur(s): Z.X. Li (LMD/CNRS) (adaptation du GCM du LMD)
c date: 19940414
c======================================================================
c
c Calculer quelques parametres pour appliquer la couche limite
c ------------------------------------------------------------
#include "dimensions.h"
#include "dimphy.h"
#include "YOMCST.h"
#include "indicesol.h"
#include "iniprint.h"
      REAL tau_gl ! temps de relaxation pour la glace de mer
ccc      PARAMETER (tau_gl=86400.0*30.0)
      PARAMETER (tau_gl=86400.0*5.0)
      REAL mx_eau_sol
      PARAMETER (mx_eau_sol=150.0)
c
      REAL calsol, calsno, calice ! epaisseur du sol: 0.15 m
      PARAMETER (calsol=1.0/(2.5578E+06*0.15))
      PARAMETER (calsno=1.0/(2.3867E+06*0.15))
      PARAMETER (calice=1.0/(5.1444E+06*0.15))
C
      INTEGER i
c
      REAL dtime
      REAL snow(klon), qsol(klon)
      INTEGER indice, knon
C
      REAL vbeta(klon)
      REAL vcal(klon)
      REAL vdif(klon)
C

      IF (indice.EQ.is_oce) THEN
      DO i = 1, knon
          vcal(i)   = 0.0
          vbeta(i)  = 1.0
          vdif(i) = 0.0
      ENDDO
      ENDIF
c
      IF (indice.EQ.is_sic) THEN
      DO i = 1, knon
          vcal(i) = calice
          IF (snow(i) .GT. 0.0) vcal(i) = calsno
          vbeta(i)  = 1.0
          vdif(i) = 1.0/tau_gl
ccc          vdif(i) = calice/tau_gl ! c'etait une erreur
      ENDDO
      ENDIF
c
      IF (indice.EQ.is_ter) THEN
      DO i = 1, knon
          vcal(i) = calsol
          IF (snow(i) .GT. 0.0) vcal(i) = calsno
          vbeta(i)  = MIN(2.0*qsol(i)/mx_eau_sol, 1.0)
          vdif(i) = 0.0
      ENDDO
      ENDIF
c
      IF (indice.EQ.is_lic) THEN
      DO i = 1, knon
          vcal(i) = calice
          IF (snow(i) .GT. 0.0) vcal(i) = calsno
          vbeta(i)  = 1.0
          vdif(i) = 0.0
      ENDDO
      ENDIF
c
      RETURN
      END
C======================================================================
      SUBROUTINE nonlocal(knon, paprs, pplay,
     .                    tsol,beta,u,v,t,q,
     .                    cd_h, cd_m, pcfh, pcfm, cgh, cgq)
      IMPLICIT none
c======================================================================
c Laurent Li (LMD/CNRS), le 30 septembre 1998
c Couche limite non-locale. Adaptation du code du CCM3.
c Code non teste, donc a ne pas utiliser.
c======================================================================
c Nonlocal scheme that determines eddy diffusivities based on a
c diagnosed boundary layer height and a turbulent velocity scale.
c Also countergradient effects for heat and moisture are included.
c
c For more information, see Holtslag, A.A.M., and B.A. Boville, 1993:
c Local versus nonlocal boundary-layer diffusion in a global climate
c model. J. of Climate, vol. 6, 1825-1842.
c======================================================================
#include "dimensions.h"
#include "dimphy.h"
#include "YOMCST.h"
#include "iniprint.h"
c
c Arguments:
c
      INTEGER knon ! nombre de points a calculer
      REAL tsol(klon) ! temperature du sol (K)
      REAL beta(klon) ! efficacite d'evaporation (entre 0 et 1)
      REAL paprs(klon,klev+1) ! pression a inter-couche (Pa)
      REAL pplay(klon,klev)   ! pression au milieu de couche (Pa)
      REAL u(klon,klev) ! vitesse U (m/s)
      REAL v(klon,klev) ! vitesse V (m/s)
      REAL t(klon,klev) ! temperature (K)
      REAL q(klon,klev) ! vapeur d'eau (kg/kg)
      REAL cd_h(klon) ! coefficient de friction au sol pour chaleur
      REAL cd_m(klon) ! coefficient de friction au sol pour vitesse
c
      INTEGER isommet
      PARAMETER (isommet=klev)
      REAL vk
      PARAMETER (vk=0.40)
      REAL ricr
      PARAMETER (ricr=0.4)
      REAL fak
      PARAMETER (fak=8.5)
      REAL fakn
      PARAMETER (fakn=7.2)
      REAL onet
      PARAMETER (onet=1.0/3.0)
      REAL t_coup
      PARAMETER(t_coup=273.15)
      REAL zkmin
      PARAMETER (zkmin=0.01)
      REAL betam
      PARAMETER (betam=15.0)
      REAL betah
      PARAMETER (betah=15.0)
      REAL betas
      PARAMETER (betas=5.0)
      REAL sffrac
      PARAMETER (sffrac=0.1)
      REAL binm
      PARAMETER (binm=betam*sffrac)
      REAL binh
      PARAMETER (binh=betah*sffrac)
      REAL ccon
      PARAMETER (ccon=fak*sffrac*vk)
c
      REAL z(klon,klev)
      REAL pcfm(klon,klev), pcfh(klon,klev)
c
      INTEGER i, k
      REAL zxt, zxq, zxu, zxv, zxmod, taux, tauy
      REAL zx_alf1, zx_alf2 ! parametres pour extrapolation
      REAL khfs(klon)       ! surface kinematic heat flux [mK/s]
      REAL kqfs(klon)       ! sfc kinematic constituent flux [m/s]
      REAL heatv(klon)      ! surface virtual heat flux
      REAL ustar(klon)
      REAL rino(klon,klev)  ! bulk Richardon no. from level to ref lev
      LOGICAL unstbl(klon)  ! pts w/unstbl pbl (positive virtual ht flx)
      LOGICAL stblev(klon)  ! stable pbl with levels within pbl
      LOGICAL unslev(klon)  ! unstbl pbl with levels within pbl
      LOGICAL unssrf(klon)  ! unstb pbl w/lvls within srf pbl lyr
      LOGICAL unsout(klon)  ! unstb pbl w/lvls in outer pbl lyr
      LOGICAL check(klon)   ! True=>chk if Richardson no.>critcal
      REAL pblh(klon)
      REAL cgh(klon,2:klev) ! counter-gradient term for heat [K/m]
      REAL cgq(klon,2:klev) ! counter-gradient term for constituents
      REAL cgs(klon,2:klev) ! counter-gradient star (cg/flux)
      REAL obklen(klon)
      REAL ztvd, ztvu, zdu2
      REAL therm(klon)      ! thermal virtual temperature excess
      REAL phiminv(klon)    ! inverse phi function for momentum
      REAL phihinv(klon)    ! inverse phi function for heat
      REAL wm(klon)         ! turbulent velocity scale for momentum
      REAL fak1(klon)       ! k*ustar*pblh
      REAL fak2(klon)       ! k*wm*pblh
      REAL fak3(klon)       ! fakn*wstr/wm
      REAL pblk(klon)       ! level eddy diffusivity for momentum
      REAL pr(klon)         ! Prandtl number for eddy diffusivities
      REAL zl(klon)         ! zmzp / Obukhov length
      REAL zh(klon)         ! zmzp / pblh
      REAL zzh(klon)        ! (1-(zmzp/pblh))**2
      REAL wstr(klon)       ! w*, convective velocity scale
      REAL zm(klon)         ! current level height
      REAL zp(klon)         ! current level height + one level up
      REAL zcor, zdelta, zcvm5, zxqs
      REAL fac, pblmin, zmzp, term
c
#include "YOETHF.h"
#include "FCTTRE.h"
c
c Initialisation
c
      DO i = 1, klon
         pcfh(i,1) = cd_h(i)
         pcfm(i,1) = cd_m(i)
      ENDDO
      DO k = 2, klev
      DO i = 1, klon
         pcfh(i,k) = zkmin
         pcfm(i,k) = zkmin
         cgs(i,k) = 0.0
         cgh(i,k) = 0.0
         cgq(i,k) = 0.0
      ENDDO
      ENDDO
c
c Calculer les hauteurs de chaque couche
c
      DO i = 1, knon
         z(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))
     .               * (paprs(i,1)-pplay(i,1)) / RG
      ENDDO
      DO k = 2, klev
      DO i = 1, knon
         z(i,k) = z(i,k-1)
     .              + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)
     .                   * (pplay(i,k-1)-pplay(i,k)) / RG
      ENDDO
      ENDDO
c
      DO i = 1, knon
         IF (thermcep) THEN
           zdelta=MAX(0.,SIGN(1.,RTT-tsol(i)))
           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
           zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q(i,1))
           zxqs= r2es * FOEEW(tsol(i),zdelta)/paprs(i,1)
           zxqs=MIN(0.5,zxqs)
           zcor=1./(1.-retv*zxqs)
           zxqs=zxqs*zcor
         ELSE
           IF (tsol(i).LT.t_coup) THEN
              zxqs = qsats(tsol(i)) / paprs(i,1)
           ELSE
              zxqs = qsatl(tsol(i)) / paprs(i,1)
           ENDIF
         ENDIF
        zx_alf1 = 1.0
        zx_alf2 = 1.0 - zx_alf1
        zxt = (t(i,1)+z(i,1)*RG/RCPD/(1.+RVTMP2*q(i,1)))
     .        *(1.+RETV*q(i,1))*zx_alf1
     .      + (t(i,2)+z(i,2)*RG/RCPD/(1.+RVTMP2*q(i,2)))
     .        *(1.+RETV*q(i,2))*zx_alf2
        zxu = u(i,1)*zx_alf1+u(i,2)*zx_alf2
        zxv = v(i,1)*zx_alf1+v(i,2)*zx_alf2
        zxq = q(i,1)*zx_alf1+q(i,2)*zx_alf2
        zxmod = 1.0+SQRT(zxu**2+zxv**2)
        khfs(i) = (tsol(i)*(1.+RETV*q(i,1))-zxt) *zxmod*cd_h(i)
        kqfs(i) = (zxqs-zxq) *zxmod*cd_h(i) * beta(i)
        heatv(i) = khfs(i) + 0.61*zxt*kqfs(i)
        taux = zxu *zxmod*cd_m(i)
        tauy = zxv *zxmod*cd_m(i)
        ustar(i) = SQRT(taux**2+tauy**2)
        ustar(i) = MAX(SQRT(ustar(i)),0.01)
      ENDDO
c
      DO i = 1, knon
         rino(i,1) = 0.0
         check(i) = .TRUE.
         pblh(i) = z(i,1)
         obklen(i) = -t(i,1)*ustar(i)**3/(RG*vk*heatv(i))
      ENDDO

C
C PBL height calculation:
C Search for level of pbl. Scan upward until the Richardson number between
C the first level and the current level exceeds the "critical" value.
C
      fac = 100.0
      DO k = 1, isommet
      DO i = 1, knon
      IF (check(i)) THEN
         zdu2 = (u(i,k)-u(i,1))**2+(v(i,k)-v(i,1))**2+fac*ustar(i)**2
         zdu2 = max(zdu2,1.0e-20)
         ztvd =(t(i,k)+z(i,k)*0.5*RG/RCPD/(1.+RVTMP2*q(i,k)))
     .         *(1.+RETV*q(i,k))
         ztvu =(t(i,1)-z(i,k)*0.5*RG/RCPD/(1.+RVTMP2*q(i,1)))
     .         *(1.+RETV*q(i,1))
         rino(i,k) = (z(i,k)-z(i,1))*RG*(ztvd-ztvu)
     .               /(zdu2*0.5*(ztvd+ztvu))
         IF (rino(i,k).GE.ricr) THEN
           pblh(i) = z(i,k-1) + (z(i,k-1)-z(i,k)) *
     .              (ricr-rino(i,k-1))/(rino(i,k-1)-rino(i,k))
           check(i) = .FALSE.
         ENDIF
      ENDIF
      ENDDO
      ENDDO

C
C Set pbl height to maximum value where computation exceeds number of
C layers allowed
C
      DO i = 1, knon
        if (check(i)) pblh(i) = z(i,isommet)
      ENDDO
C
C Improve estimate of pbl height for the unstable points.
C Find unstable points (sensible heat flux is upward):
C
      DO i = 1, knon
      IF (heatv(i) .GT. 0.) THEN
        unstbl(i) = .TRUE.
        check(i) = .TRUE.
      ELSE
        unstbl(i) = .FALSE.
        check(i) = .FALSE.
      ENDIF
      ENDDO
C
C For the unstable case, compute velocity scale and the
C convective temperature excess:
C
      DO i = 1, knon
      IF (check(i)) THEN
        phiminv(i) = (1.-binm*pblh(i)/obklen(i))**onet
        wm(i)= ustar(i)*phiminv(i)
        therm(i) = heatv(i)*fak/wm(i)
        rino(i,1) = 0.0
      ENDIF
      ENDDO
C
C Improve pblh estimate for unstable conditions using the
C convective temperature excess:
C
      DO k = 1, isommet
      DO i = 1, knon
      IF (check(i)) THEN
         zdu2 = (u(i,k)-u(i,1))**2+(v(i,k)-v(i,1))**2+fac*ustar(i)**2
         zdu2 = max(zdu2,1.0e-20)
         ztvd =(t(i,k)+z(i,k)*0.5*RG/RCPD/(1.+RVTMP2*q(i,k)))
     .         *(1.+RETV*q(i,k))
         ztvu =(t(i,1)+therm(i)-z(i,k)*0.5*RG/RCPD/(1.+RVTMP2*q(i,1)))
     .         *(1.+RETV*q(i,1))
         rino(i,k) = (z(i,k)-z(i,1))*RG*(ztvd-ztvu)
     .               /(zdu2*0.5*(ztvd+ztvu))
         IF (rino(i,k).GE.ricr) THEN
           pblh(i) = z(i,k-1) + (z(i,k-1)-z(i,k)) *
     .              (ricr-rino(i,k-1))/(rino(i,k-1)-rino(i,k))
           check(i) = .FALSE.
         ENDIF
      ENDIF
      ENDDO
      ENDDO
C
C Set pbl height to maximum value where computation exceeds number of
C layers allowed
C
      DO i = 1, knon
        if (check(i)) pblh(i) = z(i,isommet)
      ENDDO
C
C Points for which pblh exceeds number of pbl layers allowed;
C set to maximum
C
      DO i = 1, knon
        IF (check(i)) pblh(i) = z(i,isommet)
      ENDDO
C
C PBL height must be greater than some minimum mechanical mixing depth
C Several investigators have proposed minimum mechanical mixing depth
C relationships as a function of the local friction velocity, u*.  We
C make use of a linear relationship of the form h = c u* where c=700.
C The scaling arguments that give rise to this relationship most often
C represent the coefficient c as some constant over the local coriolis
C parameter.  Here we make use of the experimental results of Koracin
C and Berkowicz (1988) [BLM, Vol 43] for wich they recommend 0.07/f
C where f was evaluated at 39.5 N and 52 N.  Thus we use a typical mid
C latitude value for f so that c = 0.07/f = 700.
C
      DO i = 1, knon
        pblmin  = 700.0*ustar(i)
        pblh(i) = MAX(pblh(i),pblmin)
      ENDDO
C
C pblh is now available; do preparation for diffusivity calculation:
C
      DO i = 1, knon
        pblk(i) = 0.0
        fak1(i) = ustar(i)*pblh(i)*vk
C
C Do additional preparation for unstable cases only, set temperature
C and moisture perturbations depending on stability.
C
        IF (unstbl(i)) THEN
          zxt=(t(i,1)-z(i,1)*0.5*RG/RCPD/(1.+RVTMP2*q(i,1)))
     .         *(1.+RETV*q(i,1))
          phiminv(i) = (1. - binm*pblh(i)/obklen(i))**onet
          phihinv(i) = sqrt(1. - binh*pblh(i)/obklen(i))
          wm(i)      = ustar(i)*phiminv(i)
          fak2(i)    = wm(i)*pblh(i)*vk
          wstr(i)    = (heatv(i)*RG*pblh(i)/zxt)**onet
          fak3(i)    = fakn*wstr(i)/wm(i)
        ENDIF
      ENDDO

C Main level loop to compute the diffusivities and
C counter-gradient terms:
C
      DO 1000 k = 2, isommet
C
C Find levels within boundary layer:
C
        DO i = 1, knon
          unslev(i) = .FALSE.
          stblev(i) = .FALSE.
          zm(i) = z(i,k-1)
          zp(i) = z(i,k)
          IF (zkmin.EQ.0.0 .AND. zp(i).GT.pblh(i)) zp(i) = pblh(i)
          IF (zm(i) .LT. pblh(i)) THEN
            zmzp = 0.5*(zm(i) + zp(i))
            zh(i) = zmzp/pblh(i)
            zl(i) = zmzp/obklen(i)
            zzh(i) = 0.
            IF (zh(i).LE.1.0) zzh(i) = (1. - zh(i))**2
C
C stblev for points zm < plbh and stable and neutral
C unslev for points zm < plbh and unstable
C
            IF (unstbl(i)) THEN
              unslev(i) = .TRUE.
            ELSE
              stblev(i) = .TRUE.
            ENDIF
          ENDIF
        ENDDO
C
C Stable and neutral points; set diffusivities; counter-gradient
C terms zero for stable case:
C
        DO i = 1, knon
          IF (stblev(i)) THEN
            IF (zl(i).LE.1.) THEN
              pblk(i) = fak1(i)*zh(i)*zzh(i)/(1. + betas*zl(i))
            ELSE
              pblk(i) = fak1(i)*zh(i)*zzh(i)/(betas + zl(i))
            ENDIF
            pcfm(i,k) = pblk(i)
            pcfh(i,k) = pcfm(i,k)
          ENDIF
        ENDDO
C
C unssrf, unstable within surface layer of pbl
C unsout, unstable within outer   layer of pbl
C
        DO i = 1, knon
          unssrf(i) = .FALSE.
          unsout(i) = .FALSE.
          IF (unslev(i)) THEN
            IF (zh(i).lt.sffrac) THEN
              unssrf(i) = .TRUE.
            ELSE
              unsout(i) = .TRUE.
            ENDIF
          ENDIF
        ENDDO
C
C Unstable for surface layer; counter-gradient terms zero
C
        DO i = 1, knon
          IF (unssrf(i)) THEN
            term = (1. - betam*zl(i))**onet
            pblk(i) = fak1(i)*zh(i)*zzh(i)*term
            pr(i) = term/sqrt(1. - betah*zl(i))
          ENDIF
        ENDDO
C
C Unstable for outer layer; counter-gradient terms non-zero:
C
        DO i = 1, knon
          IF (unsout(i)) THEN
            pblk(i) = fak2(i)*zh(i)*zzh(i)
            cgs(i,k) = fak3(i)/(pblh(i)*wm(i))
            cgh(i,k) = khfs(i)*cgs(i,k)
            pr(i) = phiminv(i)/phihinv(i) + ccon*fak3(i)/fak
            cgq(i,k) = kqfs(i)*cgs(i,k)
          ENDIF
        ENDDO
C
C For all unstable layers, set diffusivities
C
        DO i = 1, knon
        IF (unslev(i)) THEN
            pcfm(i,k) = pblk(i)
            pcfh(i,k) = pblk(i)/pr(i)
        ENDIF
        ENDDO
 1000 continue           ! end of level loop

      RETURN
      END
