source: LMDZ4/trunk/libf/phylmd/coef_diff_turb_mod.F90 @ 862

!
! $Header$
!
MODULE coef_diff_turb_mod
!
! This module contains some procedures for calculation of the coefficients of the
! turbulent diffusion in the atmosphere and coefficients for turbulent diffusion 
! at surface(cdrag)
!
  USE dimphy
  
  IMPLICIT NONE
  
CONTAINS
!
!****************************************************************************************
!
  SUBROUTINE coef_diff_turb(dtime, nsrf, knon, ni, &
       ypaprs, ypplay, yu, yv, yq, yt, yts, yrugos, yqsurf, &
       ycoefm, ycoefh)
! 
! Calculate coefficients(ycoefm, ycoefh) for turbulent diffusion in the 
! atmosphere and coefficients for turbulent diffusion at surface - cdrag
! (ycoefm(:,1), ycoefh(:,1)).
!
!
! Input arguments
!****************************************************************************************
    REAL, INTENT(IN)                           :: dtime
    INTEGER, INTENT(IN)                        :: nsrf, knon
    INTEGER, DIMENSION(klon), INTENT(IN)       :: ni
    REAL, DIMENSION(klon,klev+1), INTENT(IN)   :: ypaprs
    REAL, DIMENSION(klon,klev), INTENT(IN)     :: ypplay
    REAL, DIMENSION(klon,klev), INTENT(IN)     :: yu, yv
    REAL, DIMENSION(klon,klev), INTENT(IN)     :: yq, yt
    REAL, DIMENSION(klon), INTENT(IN)          :: yts, yrugos, yqsurf
  
! Output arguments
!****************************************************************************************
    REAL, DIMENSION(klon,klev), INTENT(OUT)    :: ycoefh
    REAL, DIMENSION(klon,klev), INTENT(OUT)    :: ycoefm

! Local variables with attribute SAVE
!****************************************************************************************
    REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE  :: q2save
    !$OMP THREADPRIVATE(q2save)
    LOGICAL, SAVE                              :: first_call=.TRUE.
    !$OMP THREADPRIVATE(first_call)    

! Other local variables
!****************************************************************************************
    INTEGER                                    :: k, i, j
    REAL, DIMENSION(klon,klev)                 :: ycoefm0, ycoefh0, yzlay, yteta
    REAL, DIMENSION(klon,klev+1)               :: yzlev, yq2, q2diag, ykmm, ykmn, ykmq
    REAL, DIMENSION(klon)                      :: y_cd_h, y_cd_m
    REAL, DIMENSION(klon)                      :: yustar

! Include
!****************************************************************************************
    INCLUDE "clesphys.h"
    INCLUDE "indicesol.h"
    INCLUDE "iniprint.h"
    INCLUDE "compbl.h"
    INCLUDE "YOETHF.h"
    INCLUDE "YOMCST.h"

!****************************************************************************************
! Start calculation
! - Initilalize output variables 
!****************************************************************************************

    ycoefm(:,:) = 0.0
    ycoefh(:,:) = 0.0


!****************************************************************************************    
! Methode 1 :
!
!****************************************************************************************    

    CALL coefkz(nsrf, knon, ypaprs, ypplay, &
         ksta, ksta_ter, &
         yts, yrugos, yu, yv, yt, yq, &
         yqsurf, &
         ycoefm, ycoefh)
  
!****************************************************************************************
! Methode 2 :
!
!****************************************************************************************

    IF (iflag_pbl.EQ.1) THEN
       CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, &
            ycoefm0, ycoefh0)

       DO k = 1, klev
          DO i = 1, knon
             ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
             ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
          ENDDO
       ENDDO
    ENDIF

!****************************************************************************************
! IM cf JLD : on seuille ycoefm et ycoefh
!
!****************************************************************************************
    IF (nsrf.EQ.is_oce) THEN
       DO j=1,knon
          ycoefm(j,1)=MIN(ycoefm(j,1),cdmmax)
          ycoefh(j,1)=MIN(ycoefh(j,1),cdhmax)
       ENDDO
    ENDIF
  
!****************************************************************************************  
! Calcul d'une diffusion minimale pour les conditions tres stables
!
!****************************************************************************************
    IF (ok_kzmin) THEN
       CALL coefkzmin(knon,ypaprs,ypplay,yu,yv,yt,yq,ycoefm, &
            ycoefm0,ycoefh0)
       
       DO k = 1, klev
          DO i = 1, knon
             ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
             ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
          ENDDO
       ENDDO
       
    ENDIF

  
!****************************************************************************************
! MELLOR ET YAMADA adapte a Mars Richard Fournier et Frederic Hourdin
! Methode 3 :
!****************************************************************************************

    IF (iflag_pbl.GE.3) THEN

       IF (first_call) THEN
          PRINT*, "Using method MELLOR&YAMADA"
          ! NB! q2save could/should be read and written from (re)startphy.nc
          ALLOCATE(q2save(klon,klev+1,nbsrf))
          q2save(:,:,:) = 1.e-8
          first_call=.FALSE.
       END IF

       yzlay(1:knon,1)= &
            RD*yt(1:knon,1)/(0.5*(ypaprs(1:knon,1)+ypplay(1:knon,1))) &
            *(ypaprs(1:knon,1)-ypplay(1:knon,1))/RG
       DO k=2,klev
          DO i = 1, knon
             yzlay(i,k)= &
                  yzlay(i,k-1)+RD*0.5*(yt(i,k-1)+yt(i,k)) &
                  /ypaprs(i,k)*(ypplay(i,k-1)-ypplay(i,k))/RG
          END DO
       END DO

       DO k=1,klev
          DO i = 1, knon
             yteta(i,k)= &
                  yt(i,k)*(ypaprs(i,1)/ypplay(i,k))**RKAPPA &
                  *(1.+0.61*yq(i,k))
          END DO
       END DO

       yzlev(1:knon,1)=0.
       yzlev(1:knon,klev+1)=2.*yzlay(1:knon,klev)-yzlay(1:knon,klev-1)
       DO k=2,klev
          DO i = 1, knon
             yzlev(i,k)=0.5*(yzlay(i,k)+yzlay(i,k-1))
          END DO
       END DO

! Compresse q2save
       DO k = 1, klev+1
          DO j = 1, knon
             i = ni(j)
             yq2(j,k)=q2save(i,k,nsrf)
          ENDDO
       ENDDO
         
!   Bug introduit volontairement pour converger avec les resultats
!   du papier sur les thermiques.
       IF (1.EQ.1) THEN
          y_cd_m(1:knon) = ycoefm(1:knon,1)
          y_cd_h(1:knon) = ycoefh(1:knon,1)
       ELSE
          y_cd_h(1:knon) = ycoefm(1:knon,1)
          y_cd_m(1:knon) = ycoefh(1:knon,1)
       ENDIF
       
       CALL ustarhb(knon,yu,yv,y_cd_m, yustar)
     
       IF (prt_level > 9) THEN
          WRITE(lunout,*) 'USTAR = ',yustar
       ENDIF
         
!   iflag_pbl peut etre utilise comme longuer de melange
       IF (iflag_pbl.GE.11) THEN
          CALL vdif_kcay(knon,dtime,RG,RD,ypaprs,yt, &
               yzlev,yzlay,yu,yv,yteta, &
               y_cd_m,yq2,q2diag,ykmm,ykmn,yustar, &
               iflag_pbl)
       ELSE
          CALL yamada4(knon,dtime,RG,RD,ypaprs,yt, &
               yzlev,yzlay,yu,yv,yteta, &
               y_cd_m,yq2,ykmm,ykmn,ykmq,yustar, &
               iflag_pbl)
       ENDIF
       
       ycoefm(1:knon,1)=y_cd_m(1:knon)
       ycoefh(1:knon,1)=y_cd_h(1:knon)
       ycoefm(1:knon,2:klev)=ykmm(1:knon,2:klev)
       ycoefh(1:knon,2:klev)=ykmn(1:knon,2:klev)
                
! update q2save with yq2
       DO j=1,knon
          DO k=1,klev+1
             i=ni(j)
             q2save(i,k,nsrf) = yq2(j,k)
          END DO
       END DO
       
    ENDIF !(iflag_pbl.ge.3)

  END SUBROUTINE coef_diff_turb
!
!****************************************************************************************
!
  SUBROUTINE coefkz(nsrf, knon, paprs, pplay, &
       ksta, ksta_ter, &
       ts, rugos, &
       u,v,t,q, &
       qsurf, &
       pcfm, pcfh)
    
!======================================================================
! Auteur(s) F. Hourdin, M. Forichon, Z.X. Li (LMD/CNRS) date: 19930922
!           (une version strictement identique a l'ancien modele)
! Objet: calculer le coefficient du frottement du sol (Cdrag) et les
!        coefficients d'echange turbulent dans l'atmosphere.
! Arguments:
! nsrf-----input-I- indicateur de la nature du sol
! knon-----input-I- nombre de points a traiter
! paprs----input-R- pregssion a chaque intercouche (en Pa)
! pplay----input-R- pression au milieu de chaque couche (en Pa)
! ts-------input-R- temperature du sol (en Kelvin)
! rugos----input-R- longeur de rugosite (en m)
! u--------input-R- vitesse u
! v--------input-R- vitesse v
! t--------input-R- temperature (K)
! q--------input-R- vapeur d'eau (kg/kg)
!
! pcfm-----output-R- coefficients a calculer (vitesse)
! pcfh-----output-R- coefficients a calculer (chaleur et humidite)
!======================================================================
    INCLUDE "YOETHF.h"
    INCLUDE "FCTTRE.h"
    INCLUDE "iniprint.h"
    INCLUDE "indicesol.h"
    INCLUDE "compbl.h"
    INCLUDE "YOMCST.h"
!
! Arguments:
!
    INTEGER, INTENT(IN)                      :: knon, nsrf
    REAL, DIMENSION(klon), INTENT(IN)        :: ts
    REAL, DIMENSION(klon,klev+1), INTENT(IN) ::  paprs
    REAL, DIMENSION(klon,klev), INTENT(IN)   ::  pplay
    REAL, DIMENSION(klon,klev), INTENT(IN)   :: u, v, t, q
    REAL, DIMENSION(klon), INTENT(IN)        :: rugos

    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: pcfm, pcfh


! Local variables

! numero de couche du sommet de la couche limite
    INTEGER, DIMENSION(klon)    :: itop
!
! Quelques constantes et options:
!
    REAL cepdu2, ckap, cb, cc, cd, clam
    PARAMETER (cepdu2 =(0.1)**2)
    PARAMETER (CKAP=0.4)
    PARAMETER (cb=5.0)
    PARAMETER (cc=5.0)
    PARAMETER (cd=5.0)
    PARAMETER (clam=160.0)
    REAL ratqs ! largeur de distribution de vapeur d'eau
    PARAMETER (ratqs=0.05)
    LOGICAL richum ! utilise le nombre de Richardson humide
    PARAMETER (richum=.TRUE.)
    REAL ric ! nombre de Richardson critique
    PARAMETER(ric=0.4)
    REAL prandtl
    PARAMETER (prandtl=0.4)
    REAL kstable ! diffusion minimale (situation stable)
    ! GKtest
    ! PARAMETER (kstable=1.0e-10)
    REAL ksta, ksta_ter
!IM: 261103     REAL kstable_ter, kstable_sinon
!IM: 211003 cf GK   PARAMETER (kstable_ter = 1.0e-6)
!IM: 261103     PARAMETER (kstable_ter = 1.0e-8)
!IM: 261103   PARAMETER (kstable_ter = 1.0e-10)
!IM: 261103   PARAMETER (kstable_sinon = 1.0e-10)
    ! fin GKtest
    REAL mixlen ! constante controlant longueur de melange
    PARAMETER (mixlen=35.0)
    INTEGER isommet ! le sommet de la couche limite
    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.)

!
! Variables locales:

    INTEGER i, k !IM 120704
    REAL zgeop(klon,klev)
    REAL zmgeom(klon)
    REAL zri(klon)
    REAL zl2(klon)
    
    REAL u1(klon), v1(klon), t1(klon), q1(klon), z1(klon)
    REAL pcfm1(klon), pcfh1(klon)

    REAL zdphi, zdu2, ztvd, ztvu, zcdn
    REAL zscf
    REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs
    REAL z2geomf, zalh2, zalm2, zscfh, zscfm
    REAL t_coup
    PARAMETER (t_coup=273.15)
!IM
    LOGICAL check
    PARAMETER (check=.FALSE.)
!
! contre-gradient pour la chaleur sensible: Kelvin/metre
    REAL gamt(2:klev)
    REAL qsurf(klon) 

    LOGICAL, SAVE :: appel1er
    !$OMP THREADPRIVATE(appel1er)
!
! Fonctions thermodynamiques et fonctions d'instabilite
    REAL fsta, fins, x
    LOGICAL zxli ! utiliser un jeu de fonctions simples
    PARAMETER (zxli=.FALSE.)

    fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))
    fins(x) = SQRT(1.0-18.0*x)

    DATA appel1er /.TRUE./


    isommet=klev
      
    IF (appel1er) THEN
       IF (prt_level > 9) THEN
          WRITE(lunout,*)'coefkz, opt_ec:', opt_ec
          WRITE(lunout,*)'coefkz, richum:', richum
          IF (richum) WRITE(lunout,*)'coefkz, ratqs:', ratqs
          WRITE(lunout,*)'coefkz, isommet:', isommet
          WRITE(lunout,*)'coefkz, tvirtu:', tvirtu
          appel1er = .FALSE.
       ENDIF
    ENDIF
!
! Initialiser les sorties
!
    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
    
!
! Prescrire la valeur de contre-gradient
!
    IF (iflag_pbl.EQ.1) THEN
       DO k = 3, klev
          gamt(k) = -1.0E-03
       ENDDO
       gamt(2) = -2.5E-03
    ELSE
       DO k = 2, klev
          gamt(k) = 0.0
       ENDDO
    ENDIF
!IM cf JLD/ GKtest
    IF ( nsrf .NE. is_oce ) THEN
!IM 261103     kstable = kstable_ter
       kstable = ksta_ter
    ELSE
!IM 261103     kstable = kstable_sinon
       kstable = ksta
    ENDIF
!IM cf JLD/ GKtest fin

!
! Calculer les geopotentiels de chaque couche
!
    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

!
! Calculer le frottement au sol (Cdrag)
!
    DO i = 1, knon
       u1(i) = u(i,1)
       v1(i) = v(i,1)
       t1(i) = t(i,1)
       q1(i) = q(i,1)
       z1(i) = zgeop(i,1)
    ENDDO

    CALL clcdrag(klon, knon, nsrf, zxli, &
         u1, v1, t1, q1, z1, &
         ts, qsurf, rugos, &
         pcfm1, pcfh1) 
!IM               ts, qsurf, rugos,

    DO i = 1, knon
       pcfm(i,1)=pcfm1(i)
       pcfh(i,1)=pcfh1(i)
    ENDDO
!
! Calculer les coefficients turbulents dans l'atmosphere
!
    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

!
! Calculer Qs et dQs/dT:
!
          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
!
!           calculer la fraction nuageuse (processus humide):
!
          zfr = (zq+ratqs*zq-zqs) / (2.0*ratqs*zq)
          zfr = MAX(0.0,MIN(1.0,zfr))
          IF (.NOT.richum) zfr = 0.0
!
!           calculer le nombre de Richardson:
!
          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))

          ELSE ! calcul de Ridchardson compatible LMD5

             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 &
                  *(paprs(i,k)/101325.0)**RKAPPA &
                  /(zdu2*0.5*(t(i,k-1)+t(i,k)))
          ENDIF
!
!           finalement, les coefficients d'echange sont obtenus:
!
          zcdn=SQRT(zdu2) / zmgeom(i) * RG

          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

!
! Au-dela du sommet, pas de diffusion turbulente:
!
    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
      
  END SUBROUTINE coefkz
!
!****************************************************************************************
!
  SUBROUTINE coefkz2(nsrf, knon, paprs, pplay,t, &
       pcfm, pcfh)

!======================================================================
! J'introduit un peu de diffusion sauf dans les endroits
! ou une forte inversion est presente
! On peut dire qu'il represente la convection peu profonde
!
! Arguments:
! nsrf-----input-I- indicateur de la nature du sol
! knon-----input-I- nombre de points a traiter
! paprs----input-R- pression a chaque intercouche (en Pa)
! pplay----input-R- pression au milieu de chaque couche (en Pa)
! t--------input-R- temperature (K)
!
! pcfm-----output-R- coefficients a calculer (vitesse)
! pcfh-----output-R- coefficients a calculer (chaleur et humidite)
!======================================================================
!
! Arguments:
!
    INTEGER, INTENT(IN)                       :: knon, nsrf
    REAL, DIMENSION(klon, klev+1), INTENT(IN) ::  paprs
    REAL, DIMENSION(klon, klev), INTENT(IN)   ::  pplay
    REAL, DIMENSION(klon, klev), INTENT(IN)   :: t(klon,klev)
    
    REAL, DIMENSION(klon, klev), INTENT(OUT)  :: pcfm, pcfh
!
! Quelques constantes et options:
!
    REAL prandtl
    PARAMETER (prandtl=0.4)
    REAL kstable
    PARAMETER (kstable=0.002)
!    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)
!    PARAMETER (seuil=-0.04)
!    PARAMETER (seuil=-0.06)
!    PARAMETER (seuil=-0.09)

!
! Variables locales:
!
    INTEGER i, k, invb(knon)
    REAL zl2(knon)
    REAL zdthmin(knon), zdthdp

    INCLUDE "indicesol.h"
    INCLUDE "YOMCST.h"
!
! Initialiser les sorties
!
    DO k = 1, klev
       DO i = 1, knon
          pcfm(i,k) = 0.0
          pcfh(i,k) = 0.0
       ENDDO
    ENDDO

!
! Chercher la zone d'inversion forte
!
    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

!
! Introduire une diffusion:
!
    IF ( nsrf.EQ.is_oce ) THEN
       DO k = 2, klev
          DO i = 1, knon
!IM cf FH/GK   IF ( (nsrf.NE.is_oce) .OR.  ! si ce n'est pas sur l'ocean
!IM cf FH/GK  .     (invb(i).EQ.klev) .OR. ! s'il n'y a pas d'inversion
      !IM cf JLD/ GKtest TERkz2
      ! IF ( (nsrf.EQ.is_ter) .OR.  ! si on est sur la terre
      ! fin GKtest


! s'il n'y a pas d'inversion ou si l'inversion est trop faible
!          IF ( (nsrf.EQ.is_oce) .AND. &
             IF ( (invb(i).EQ.klev) .OR. (zdthmin(i).GT.seuil) ) THEN 
                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
    ENDIF

  END SUBROUTINE coefkz2
!
!****************************************************************************************
!
END MODULE coef_diff_turb_mod