SUBROUTINE clmain(dtime,pctsrf,t,q,u,v,
     .                  soil_model,ts,soilcap,soilflux,
     .                  paprs,pplay,radsol,snow,qsol,
     .                  xlat, rugos,
     .                  d_t,d_q,d_u,d_v,d_ts,
     .                  flux_t,flux_q,flux_u,flux_v,cdragh,cdragm,
     .                  rugmer, dflux_t,dflux_q,
     .                  zcoefh,zu1,zv1)
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).
 
      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 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 capsol---input-R- inversion de l'effective capacite du sol (J/m2/K)
c beta-----input-R- coefficient de l'evaporation reelle (0 a 1)
c dif_grnd-input-R- coeff. de diffusion (chaleur) vers le sol profond
c xlat-----input-R- latitude en degree
c rugos----input-R- longeur de rugosite (en 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 rugmer---output-R- longeur de rugosite sur mer (m)
c dflux_t derive du flux sensible
c dflux_q derive du flux latent
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
c======================================================================
#include "dimensions.h"
#include "dimphy.h"
#include "indicesol.h"
c
      LOGICAL soil_model
c
      REAL dtime
      REAL t(klon,klev), q(klon,klev)
      REAL u(klon,klev), v(klon,klev)
      REAL paprs(klon,klev+1), pplay(klon,klev), radsol(klon)
      REAL xlat(klon)
      REAL d_t(klon, klev), d_q(klon, klev)
      REAL d_u(klon, klev), d_v(klon, klev)
      REAL flux_t(klon,klev), flux_q(klon,klev)
      REAL dflux_t(klon), dflux_q(klon)
      REAL flux_u(klon,klev), flux_v(klon,klev)
      REAL rugmer(klon)
      REAL cdragh(klon), cdragm(klon)
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 qsol(klon,nbsrf)
      REAL rugos(klon,nbsrf)
cAA
      REAL zcoefh(klon,klev)
      REAL zu1(klon)
      REAL zv1(klon)
cAA
c======================================================================
      EXTERNAL clqh, clvent, coefkz, calbeta, cltrac
c======================================================================
      REAL yts(klon), yrugos(klon), ypct(klon)
      REAL ycal(klon), ybeta(klon), ydif(klon)
      REAL yu1(klon), yv1(klon)
      REAL yrugm(klon), yrads(klon)
      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 ygamt(klon,2:klev) ! contre-gradient pour temperature
      REAL ygamq(klon,2:klev) ! contre-gradient pour humidite
      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
      LOGICAL contreg
      PARAMETER (contreg=.TRUE.)
c
      LOGICAL ok_nonloc
      PARAMETER (ok_nonloc=.FALSE.)
      REAL y_cd_h(klon), y_cd_m(klon)
      REAL ycoefm0(klon,klev), ycoefh0(klon,klev)
c
#include "YOMCST.h"
      REAL u1lay(klon), v1lay(klon)
      REAL delp(klon,klev)
      REAL capsol(klon), beta(klon), dif_grnd(klon)
      REAL cal(klon)
      REAL soilcap(klon,nbsrf), soilflux(klon,nbsrf)
      REAL totalflu(klon)
      INTEGER i, k, nsrf 
cAA   INTEGER it
      INTEGER ni(klon), knon, j
c======================================================================
      REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.
c======================================================================
      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
      DO nsrf = 1, nbsrf
      DO i = 1, klon
         d_ts(i,nsrf) = 0.0
      ENDDO
      ENDDO
      DO k = 1, klev
      DO i = 1, klon
         d_t(i,k) = 0.0
         d_q(i,k) = 0.0
         flux_t(i,k) = 0.0
         flux_q(i,k) = 0.0
         d_u(i,k) = 0.0
         d_v(i,k) = 0.0
         flux_u(i,k) = 0.0
         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
      DO 99999 nsrf = 1, nbsrf
c
c prescrire les proprietes du sol:
      CALL calbeta(dtime,nsrf,snow,qsol, beta,capsol,dif_grnd)
      IF (.NOT.soil_model) THEN
         DO i = 1, klon
            cal(i) = RCPD * capsol(i)
            totalflu(i) = radsol(i)
         ENDDO
      ELSE
         DO i = 1, klon
            totalflu(i) = soilflux(i,nsrf) + radsol(i)
            IF (nsrf.EQ.is_oce) THEN
               cal(i) = 0.0
            ELSE
               cal(i) = RCPD / soilcap(i,nsrf)
            ENDIF
         ENDDO
      ENDIF
c
c chercher les indices:
      DO j = 1, klon
         ni(j) = 0
      ENDDO
      knon = 0
      DO i = 1, klon
      IF (pctsrf(i,nsrf).GT.epsfra) THEN
         knon = knon + 1
         ni(knon) = i
      ENDIF
      ENDDO
c
      IF (knon.EQ.0) GOTO 99999
      DO j = 1, knon
      i = ni(j)
        ypct(j) = pctsrf(i,nsrf)
        yts(j) = ts(i,nsrf)
        yrugos(j) = rugos(i,nsrf)
        yu1(j) = u1lay(i)
        yv1(j) = v1lay(i)
        yrads(j) = totalflu(i)
        ycal(j) = cal(i)
        ybeta(j) = beta(i)
        ydif(j) = dif_grnd(i)
        ypaprs(j,klev+1) = paprs(i,klev+1)
      ENDDO
      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
cAA      IF (itr.GE.1) THEN
cAA      DO it = 1, itr
cAA        DO k = 1, klev
cAA        DO j = 1, knon
cAA        i = ni(j)
cAA          ytr(j,k,it) = tr(i,k,it)
cAA        ENDDO
cAA        ENDDO
cAA        DO j = 1, knon
cAA        i = ni(j)
cAA          yflxsrf(j,it) = flux_surf(i,it)
cAA        ENDDO
cAA      ENDDO
cAA      ENDIF
c
c calculer Cdrag et les coefficients d'echange
      CALL coefkz(nsrf, knon, ypaprs, ypplay,
     .            yts, yrugos, yu, yv, yt, yq,
     .            ycoefm, ycoefh)
      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
c
c parametrisation non-locale:
      IF (ok_nonloc) THEN
         DO i = 1, knon
            y_cd_h(i) = ycoefh(i,1)
            y_cd_m(i) = ycoefm(i,1)
         ENDDO
         CALL nonlocal(knon, ypaprs, ypplay,
     .        yts,ybeta,yu,yv,yt,yq,
     .        y_cd_h, y_cd_m, ycoefm0, ycoefh0, ygamt, ygamq)
         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
      ELSE
         IF (.NOT.contreg) THEN
            DO k = 2, klev
            DO i = 1, knon
               ygamq(i,k) = 0.0
               ygamt(i,k) = 0.0
            ENDDO
            ENDDO
         ELSE
         DO k = 3, klev
         DO i = 1, knon
            ygamq(i,k) = 0.0
            ygamt(i,k) = -1.0E-03
         ENDDO
         ENDDO
         DO i = 1, knon
            ygamq(i,2) = 0.0
            ygamt(i,2) = -2.5E-03
         ENDDO
         ENDIF
      ENDIF
c
c calculer la diffusion de "q" et de "h"
      CALL clqh(knon, dtime, yu1, yv1,
     e          ycoefh,yt,yq,yts,ypaprs,ypplay,ydelp,yrads,
     e          ycal,ybeta,ydif,ygamt,ygamq,
     s          y_d_t, y_d_q, y_d_ts, 
     s          y_flux_t, y_flux_q, y_dflux_t, y_dflux_q)
c
c calculer la diffusion des vitesses "u" et "v"
      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
c calculer la longueur de rugosite sur ocean
      IF (nsrf.EQ.is_oce) THEN
      DO j = 1, knon
         yrugm(j) = 0.018*ycoefm(j,1) * (yu1(j)**2+yv1(j)**2)/RG
         yrugm(j) = MAX(1.5e-05,yrugm(j))
      ENDDO
      ENDIF
c
cAA MAINTENANT DANS PHYTRAC
cAAc calculer la diffusion des traceurs
cAA      IF (itr.GE.1) THEN
cAA      DO it = 1, itr
cAA      CALL cltrac(knon,dtime,ycoefh, yt, ytr(1,1,it), yflxsrf(1,it),
cAA     e            ypaprs, ypplay, ydelp,
cAA     s            y_d_tr(1,1,it))
cAA      ENDDO
cAA      ENDIF
c
      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
         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)
         y_flux_t(j,k) = y_flux_t(j,k) * ypct(j)
         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)
         y_flux_u(j,k) = y_flux_u(j,k) * ypct(j)
         y_flux_v(j,k) = y_flux_v(j,k) * ypct(j)
      ENDDO
      ENDDO
c
      DO j = 1, knon
      i = ni(j)
         d_ts(i,nsrf) = y_d_ts(j)
         rugmer(i) = yrugm(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)
      ENDDO
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)
         flux_t(i,k) = flux_t(i,k) + y_flux_t(j,k)
         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)
         flux_u(i,k) = flux_u(i,k) + y_flux_u(j,k)
         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
cAA      IF (itr.GE.1) THEN
cAA      DO it = 1, itr
cAA      DO k = 1, klev
cAA      DO j = 1, knon
cAA         y_d_tr(j,k,it) = y_d_tr(j,k,it) * ypct(j)
cAA      ENDDO
cAA      ENDDO
cAA      ENDDO
cAA      DO j = 1, knon
cAA      i = ni(j)
cAA      DO it = 1, itr
cAA      DO k = 1, klev
cAA         d_tr(i,k,it) = d_tr(i,k,it) + y_d_tr(j,k,it)
cAA      ENDDO
cAA      ENDDO
cAA      ENDDO
cAA      ENDIF
c
99999 CONTINUE
c
      RETURN
      END
      SUBROUTINE clqh(knon,dtime,u1lay,v1lay,coef,t,q,ts,paprs,pplay,
     e                delp,radsol,cal,beta,dif_grnd, gamt,gamq,
     s                d_t, d_q, d_ts, flux_t, flux_q,dflux_s,dflux_l)
      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"
c Arguments:
      INTEGER knon
      REAL dtime              ! intervalle du temps (s)
      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 cal(klon)          ! Cp/cal indique la capacite calorifique
c                               surfacique du sol
      REAL beta(klon)         ! evap. reelle / evapotranspiration
      REAL dif_grnd(klon)     ! coeff. diffusion vers le sol profond
      REAL t(klon,klev)       ! temperature (K)
      REAL q(klon,klev)       ! humidite specifique (kg/kg)
      REAL ts(klon)           ! temperature du sol (K)
      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
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
      REAL dflux_g(klon) ! derivee du flux du sol profond dF/dTs
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_a
      REAL zx_b
      REAL zx_qs
      REAL zx_dq_s_dh
      REAL zx_h_grnd
      REAL zx_cq0(klon)
      REAL zx_dq0(klon)
      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 zx_q_0(klon)
      REAL zx_h_ts(klon)
      REAL zx_sl(klon)
      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======================================================================
C Variables intermediaires pour le calcul des fluxs a la surface
      real zx_mh(klon), zx_nh(klon), zx_oh(klon)
      real zx_mq(klon), zx_nq(klon), zx_oq(klon)
      real zx_k1(klon), zx_dq_s_dt(klon)
      real zx_qsat(klon)
c======================================================================
      REAL zcor, zdelta, zcvm5
#include "YOMCST.h"
#include "YOETHF.h"
#include "FCTTRE.h"
c======================================================================
      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
      DO i = 1, knon
         zx_coef(i,1) = coef(i,1)
     .                 * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2))
     .                 * pplay(i,1)/(RD*t(i,1))
      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
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
c======================================================================
c
c zx_qs = qsat en kg/kg
c
      DO i = 1, knon
         IF (thermcep) THEN
           zdelta=MAX(0.,SIGN(1.,rtt-local_ts(i)))
           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
           zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q(i,1))
           zx_qs= r2es * FOEEW(local_ts(i),zdelta)/paprs(i,1)
           zx_qs=MIN(0.5,zx_qs)
           zcor=1./(1.-retv*zx_qs)
           zx_qs=zx_qs*zcor
           zx_dq_s_dh = FOEDE(local_ts(i),zdelta,zcvm5,zx_qs,zcor)
     .                 /RLVTT / zx_pkh(i,1)
         ELSE
           IF (local_ts(i).LT.t_coup) THEN
              zx_qs = qsats(local_ts(i)) / paprs(i,1)
              zx_dq_s_dh = dqsats(local_ts(i),zx_qs)/RLVTT
     .                    / zx_pkh(i,1)
           ELSE
              zx_qs = qsatl(local_ts(i)) / paprs(i,1)
              zx_dq_s_dh = dqsatl(local_ts(i),zx_qs)/RLVTT
     .                    / zx_pkh(i,1)
           ENDIF
         ENDIF

c        zx_dq_s_dh = 0.
         zx_dq_s_dt(i) = RCPD * zx_pkh(i,1) * zx_dq_s_dh
         zx_qsat(i) = zx_qs
c
      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
         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
         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
         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 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
         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 === Calcul de la temperature de surface ===
C 
C zx_sl = chaleur latente d'evaporation ou de sublimation
C
      do i = 1, knon
        zx_sl(i) = RLVTT
        if (local_ts(i) .LT. RTT) zx_sl(i) = RLSTT
        zx_k1(i) = zx_coef(i,1)
      enddo
      do i = 1, knon
C Q
        zx_oq(i) = 1. - (beta(i) * zx_k1(i) * zx_dq(i,1) * dtime)
        zx_mq(i) = beta(i) * zx_k1(i) * 
     .              (zx_cq(i,1) - zx_qsat(i) 
     .                             + zx_dq_s_dt(i) * local_ts(i))
     .              / zx_oq(i)
        zx_nq(i) = beta(i) * zx_k1(i) * (-1. * zx_dq_s_dt(i))
     .                                                / zx_oq(i)
c H
        zx_oh(i) = 1. - (zx_k1(i) * zx_dh(i,1) * dtime)
        zx_mh(i) = zx_k1(i) * zx_ch(i,1) / zx_oh(i)
        zx_nh(i) = - (zx_k1(i) * RCPD * zx_pkh(i,1))/ zx_oh(i)
c Tsurface
        local_ts(i) = (ts(i) + cal(i)/(RCPD * zx_pkh(i,1)) * dtime *
     .                   (radsol(i) + zx_mh(i) + zx_sl(i) * zx_mq(i)) 
     .                       + dif_grnd(i) * t_grnd * dtime)/
     .                ( 1. - dtime * cal(i)/(RCPD * zx_pkh(i,1)) * (
     .                       zx_nh(i) + zx_sl(i) * zx_nq(i))  
     .                     + dtime * dif_grnd(i))
        zx_h_ts(i) = local_ts(i) * RCPD * zx_pkh(i,1)
        d_ts(i) = local_ts(i) - ts(i)
        zx_q_0(i) = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
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)
        flux_q(i,1) = zx_mq(i) + zx_nq(i) * local_ts(i) 
        flux_t(i,1) = zx_mh(i) + zx_nh(i) * local_ts(i)
c Derives des flux dF/dTs (W m-2 K-1):
        dflux_s(i) = zx_nh(i)
        dflux_l(i) = (zx_sl(i) * zx_nq(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"
      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,
     .                  ts, rugos,
     .                  u,v,t,q,
     .                  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 xlat-----input-R- latitude en degree
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"
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)
      REAL xlat(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.35)
      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)
      PARAMETER (kstable=1.0e-10)
      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.)
      LOGICAL contreg ! utiliser le contre-gradient dans Ri
      PARAMETER (contreg=.TRUE.)
c
c Variables locales:
c
      INTEGER i, k
      REAL zgeop(klon,klev)
      REAL zmgeom(klon)
      REAL zri(klon)
      REAL zl2(klon)
      REAL zcfm1(klon), zcfm2(klon)
      REAL zcfh1(klon), zcfh2(klon)
      REAL zdphi, zdu2, ztvd, ztvu, ztsolv, zcdn
      REAL zscf, zucf, zcr
      REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs
      REAL z2geomf, zalh2, zalm2, zscfh, zscfm
      REAL t_coup
      PARAMETER (t_coup=273.15)
c
c contre-gradient pour la chaleur sensible: Kelvin/metre
      REAL gamt(2:klev)
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
         PRINT*, 'coefkz, opt_ec:', opt_ec
         PRINT*, 'coefkz, richum:', richum
         IF (richum) PRINT*, 'coefkz, ratqs:', ratqs
         PRINT*, 'coefkz, isommet:', isommet
         PRINT*, 'coefkz, tvirtu:', tvirtu
         appel1er = .FALSE.
      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 (.NOT.contreg) THEN
         DO k = 2, klev
            gamt(k) = 0.0
         ENDDO
      ELSE
         DO k = 3, klev
            gamt(k) = -1.0E-03
         ENDDO
         gamt(2) = -2.5E-03
      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 le frottement au sol (Cdrag)
c
      DO i = 1, knon
         zdu2=max(cepdu2,u(i,1)**2+v(i,1)**2)
         zdphi=zgeop(i,1)
         ztsolv = ts(i) * (1.0+RETV*q(i,1)) ! qsol approx = q(i,1)
         ztvd=(t(i,1)+zdphi/RCPD/(1.+RVTMP2*q(i,1)))
     .       *(1.+RETV*q(i,1))
         zri(i)=zgeop(i,1)*(ztvd-ztsolv)/(zdu2*ztvd)
         zcdn = (ckap/log(1.+zgeop(i,1)/(RG*rugos(i))))**2
         IF (zri(i) .ge. 0.) THEN ! situation stable
           IF (.NOT.zxli) THEN
           zscf=SQRT(1.+cd*ABS(zri(i)))
           zcfm1(i) = zcdn/(1.+2.0*cb*zri(i)/zscf)
           zcfh1(i) = zcdn/(1.+3.0*cb*zri(i)*zscf)
           pcfm(i,1) = zcfm1(i)
           pcfh(i,1) = zcfh1(i)
           ELSE
           pcfm(i,1) = zcdn* fsta(zri(i))
           pcfh(i,1) = zcdn* fsta(zri(i))
           ENDIF
         ELSE ! situation instable
           IF (.NOT.zxli) THEN
           zucf=1./(1.+3.0*cb*cc*zcdn*SQRT(ABS(zri(i))
     .              *(1.0+zgeop(i,1)/(RG*rugos(i)))))
           zcfm2(i) = zcdn*(1.-2.0*cb*zri(i)*zucf)
           zcfh2(i) = zcdn*(1.-3.0*cb*zri(i)*zucf)
           pcfm(i,1) = zcfm2(i)
           pcfh(i,1) = zcfh2(i)
           ELSE
           pcfm(i,1) = zcdn* fins(zri(i))
           pcfh(i,1) = zcdn* fins(zri(i))
           ENDIF
           zcr = (0.0016/(zcdn*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.)
           IF(nsrf.EQ.is_oce)pcfh(i,1)=zcdn*(1.0+zcr**1.25)**(1./1.25)
         ENDIF
      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"
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
      IF ( (nsrf.NE.is_oce) .OR.  ! si ce n'est pas sur l'ocean
     .     (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,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"
      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,nbsrf), qsol(klon,nbsrf)
      INTEGER indice
C
      REAL vbeta(klon)
      REAL vcal(klon)
      REAL vdif(klon)
C
      IF (indice.EQ.is_oce) THEN
      DO i = 1, klon
          vcal(i)   = 0.0
          vbeta(i)  = 1.0
          vdif(i) = 0.0
      ENDDO
      ENDIF
c
      IF (indice.EQ.is_sic) THEN
      DO i = 1, klon
          vcal(i) = calice
          IF (snow(i,is_sic) .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, klon
          vcal(i) = calsol
          IF (snow(i,is_ter) .GT. 0.0) vcal(i) = calsno
          vbeta(i)  = MIN(2.0*qsol(i,is_ter)/mx_eau_sol, 1.0)
          vdif(i) = 0.0
      ENDDO
      ENDIF
c
      IF (indice.EQ.is_lic) THEN
      DO i = 1, klon
          vcal(i) = calice
          IF (snow(i,is_lic) .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"
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.35)
      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