c
c $Header$
c

      SUBROUTINE clmain(dtime,itap,date0,pctsrf,
     .                  t,q,u,v,
     .                  jour, rmu0,
     .                  ok_veget, ocean, npas, nexca, ts,
     .                  soil_model,ftsoil,
     .                  paprs,pplay,radsol,snow,qsol,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,
     .                  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).
      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
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$$$ PB ajout pour soil
#include "dimsoil.h"
c
      REAL dtime
      real date0
      integer itap
      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 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)
      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
      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 qsol(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)
      REAL sollw(klon), solsw(klon), 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
      REAL ftsoil(klon,nsoilmx,nbsrf)
      REAL ytsoil(klon,nsoilmx)
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), yqsol(klon), yagesno(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
      LOGICAL contreg
      PARAMETER (contreg=.TRUE.)
c
      LOGICAL ok_nonloc
      PARAMETER (ok_nonloc=.FALSE.)
      REAL ycoefm0(klon,klev), ycoefh0(klon,klev)
c
#include "YOMCST.h"
      REAL u1lay(klon), v1lay(klon)
      REAL delp(klon,klev)
      REAL totalflu(klon)
      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/.false./
      LOGICAL debugindex
      SAVE debugindex
      DATA debugindex/.false./
#include "temps.h"
      
      IF (first_appel) THEN
          first_appel=.false.
!
! initialisation sorties netcdf
!
          CALL ymds2ju(anne_ini, 1, 1, 0.0, zjulian)
          zjulian = zjulian + day_ini
          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,jjm+1,zx_lat,1,iim,1,jjm
     $        +1, 0,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
      yqsol = 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
      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.

      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
       totalflu = radsol

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
c      write(*,*)'CLMAIN, nsrf, knon =',nsrf, knon
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)
        ysnow(j) = snow(i,nsrf)
        yqsol(j) = qsol(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)
c$$$        ysolsw(j) = solsw(i)
       ysolsw(j) = (1 - albe(i,nsrf)) 
     $      /(1 - pctsrf(i,is_ter) * albe(i,is_ter) 
     $      - pctsrf(i, is_lic) *albe(i,is_lic) 
     $      - pctsrf(i, is_oce) *albe(i,is_oce)
     $      - pctsrf(i, is_sic) *albe(i,is_sic) 
     $      ) * solsw(i)
        ysollw(j) = sollw(i)
        ysollwdown(j) = sollwdown(i)
        yrugos(j) = rugos(i,nsrf)
        yrugoro(j) = rugoro(i)
        yu1(j) = u1lay(i)
        yv1(j) = v1lay(i)
c$$$        yrads(j) = totalflu(i)
        yrads(j) = (1 - albe(i,nsrf)) 
     $      /(1 - pctsrf(i,is_ter) * albe(i,is_ter) 
     $      - pctsrf(i, is_lic) *albe(i,is_lic) 
     $      - pctsrf(i, is_oce) *albe(i,is_oce)
     $      - pctsrf(i, is_sic) *albe(i,is_sic) 
     $      ) * solsw(i) + sollw(i)
        ypaprs(j,klev+1) = paprs(i,klev+1)
      END DO
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,
     .            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
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 pour le couplage
      ytaux = y_flux_u(:,1)
      ytauy = y_flux_v(:,1)

c calculer la diffusion de "q" et de "h"
      CALL clqh(dtime, itap, date0,jour, debut,lafin,
     e          rlon, rlat, cufi, cvfi,
     e          knon, nsrf, ni, pctsrf,
     e          soil_model, ytsoil,
     e          ok_veget, ocean, npas, nexca,
     e          rmu0, yrugos, yrugoro,
     e          yu1, yv1, ycoefh,
     e          yt,yq,yts,ypaprs,ypplay,
     e          ydelp,yrads,yalb, yalblw, ysnow, yqsol, 
     e          yrain_f, ysnow_f, yfder, ytaux, ytauy, 
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)
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
         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.
      qsol(:, 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)
         qsol(i,nsrf) = yqsol(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(i)
          endif	
         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
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
99999 CONTINUE
c
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,
     e                ok_veget, ocean, npas, nexca,
     e                rmu0, rugos, rugoro,
     e                u1lay,v1lay,coef,
     e                t,q,ts,paprs,pplay,
     e                delp,radsol,albedo,alblw,snow,qsol, 
     e                precip_rain, precip_snow, fder, taux, tauy,
     $                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)

      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"
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 qsol(klon)         ! humidite de la surface
      real precip_rain(klon), precip_snow(klon)
      REAL agesno(klon)
      REAL rugoro(klon)
      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 
      character*6 ocean
      integer npas, nexca

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
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======================================================================
      logical contreg
      parameter (contreg=.true.)
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)

! 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

      if (.not. contreg) then
        do k = 2, klev
          do i = 1, knon
            gamq(i,k) = 0.0
            gamt(i,k) = 0.0
          enddo
        enddo
      else
        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
      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
         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 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)
        swdown(1:knon) = swnet(1:knon)
c      enddo
c En attendant mieux
      ccanopy = 365.

      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,
     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, rugos, rugoro,
     e albedo, snow, qsol,
     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)


      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"
      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 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)
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 essai qsurf
      real qsurf(klon)
      real friv, frih
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

      do i = 1, knon
        qsurf(i) = qsatl(ts(i))/paprs(i,1)
      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
      CALL clcdrag(knon, nsrf, is_oce, zxli, u, v, t, q, zgeop,
     .             ts, qsurf, rugos, pcfm, pcfh, zcdn, zri)
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 clcdrag(knon, nsrf, is_oce, zxli, 
     .                   u, v, t, q, zgeop, 
     .                   ts, qsurf, rugos,
     .                   pcfm, pcfh, zcdn, zri)
c ================================================================= c
c Objet : calcul cdrags pour le moment et les flux chaleur sensible, latente (pcfm,pcfh)
c         et du nombre de Richardson zri
c ================================================================= c
      IMPLICIT NONE
#include "dimensions.h"
#include "dimphy.h"
#include "YOMCST.h"
#include "indicesol.h"
c
      INTEGER knon, nsrf, is_oce
      REAL ts(klon), qsurf(klon)
      REAL u(klon,klev), v(klon,klev), t(klon,klev), q(klon,klev)
      REAL zgeop(klon,klev)
      REAL rugos(klon), zri(klon)
c
      REAL zcdn
      REAL pcfm(klon,klev), pcfh(klon,klev)
c
c Quelques constantes et options:
c
      REAL ckap, cb, cc, cd, cepdu2
      PARAMETER (ckap=0.35)
      PARAMETER (cb=5.0)
      PARAMETER (cc=5.0)
      PARAMETER (cd=5.0)
      PARAMETER (cepdu2 =(0.1)**2)
c
c Variables locales
      INTEGER i
      REAL zdu2, zdphi, ztsolv, ztvd, zscf, zucf, zcr
      REAL friv, frih
      REAL zcfm1(klon), zcfm2(klon)
      REAL zcfh1(klon), zcfh2(klon)
c
c Fonctions thermodynamiques et fonctions d'instabilite
      REAL fsta, fins, x
      LOGICAL zxli
      fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))
      fins(x) = SQRT(1.0-18.0*x)
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)
c        ztsolv = ts(i) * (1.0+RETV*q(i,1)) ! qsol approx = q(i,1)
         ztsolv = ts(i) * (1.0+RETV*qsurf(i))
         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)))
           FRIV = AMAX1(1. / (1.+2.*CB*zri(i)/ZSCF), 0.1)
!           zcfm1(i) = zcdn/(1.+2.0*cb*zri(i)/ zscf)
           zcfm1(i) = zcdn * FRIV
           FRIH = AMAX1(1./ (1.+3.*CB*zri(i)*ZSCF), 0.1 )
!           zcfh1(i) = zcdn/(1.+3.0*cb*zri(i)*zscf)
           zcfh1(i) = zcdn * FRIH
           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*amax1((1.-2.0*cb*zri(i)*zucf),0.1)
           zcfh2(i) = zcdn*amax1((1.-3.0*cb*zri(i)*zucf),0.1)
           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
      END DO
      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,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"
      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"
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