!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  USE dimphy, only : klev,klon

  USE phys_local_var_mod, only: tke_dissip,wprime

  USE yamada_ini_mod, only : new_yamada4,yamada4_num,hboville

  USE yamada_ini_mod, only : prt_level, lunout,pbl_lmixmin_alpha,b1,kap,viscom,viscoh

  USE yamada_ini_mod, only : ric, yun,ydeux,lmixmin

  IMPLICIT NONE

  ! ************************************************************************************************

  !

  ! yamada4: subroutine qui calcule le transfert turbulent avec une fermeture d'ordre 2 ou 2.5

  !

  ! Reference: Simulation of nocturnal drainage flows by a q2l Turbulence Closure Model

  !            Yamada T.

  !            J. Atmos. Sci, 40, 91-106, 1983

  !

  !************************************************************************************************

  ! Input :

  !'======

  ! ni: indice horizontal sur la grille de base, non restreinte

  ! nsrf: type de surface

  ! ngrid: nombre de mailles concern??es sur l'horizontal

  ! dt : pas de temps

  ! g  : g

  ! rconst: constante de l'air sec

  ! zlev : altitude a chaque niveau (interface inferieure de la couche

  ! de meme indice)

  ! zlay : altitude au centre de chaque couche

  ! u,v : vitesse au centre de chaque couche

  ! (en entree : la valeur au debut du pas de temps)

  ! teta : temperature potentielle virtuelle au centre de chaque couche

  ! (en entree : la valeur au debut du pas de temps)

  ! cd : cdrag pour la quantit?? de mouvement

  ! (en entree : la valeur au debut du pas de temps)

  ! ustar: vitesse de friction calcul??e par une formule diagnostique

  ! iflag_pbl: flag pour choisir des options du sch??ma de turbulence

  !

  !             iflag_pbl doit valoir entre 6 et 9

  !             l=6, on prend  systematiquement une longueur d'equilibre

  !             iflag_pbl=6 : MY 2.0

  !             iflag_pbl=7 : MY 2.0.Fournier

  !             iflag_pbl=8/9 : MY 2.5

  !             iflag_pbl=8 with special obsolete treatments for convergence

  !             with Cmpi5 NPv3.1 simulations

  !             iflag_pbl=10/11 :  New scheme M2 and N2 explicit and dissiptation exact

  !             iflag_pbl=12 = 11 with vertical diffusion off q2

  !

  !             2013/04/01 (FH hourdin@lmd.jussieu.fr)

  !             Correction for very stable PBLs (iflag_pbl=10 and 11)

  !             iflag_pbl=8 converges numerically with NPv3.1

  !             iflag_pbl=11 -> the model starts with NP from start files created by ce0l

  !                          -> the model can run with longer time-steps.

  !             2016/11/30 (EV etienne.vignon@lmd.ipsl.fr)

  !               On met tke (=q2/2) en entr??e plut??t que q2

  !               On corrige l'update de la tke

  !             2020/10/01 (EV)

  !               On ajoute la dissipation de la TKE en diagnostique de sortie

  !

  ! Inpout/Output :

  !==============

  ! tke : tke au bas de chaque couche

  ! (en entree : la valeur au debut du pas de temps)

  ! (en sortie : la valeur a la fin du pas de temps)

  ! Outputs:

  !==========

  ! km : diffusivite turbulente de quantite de mouvement (au bas de chaque

  ! couche)

  ! (en sortie : la valeur a la fin du pas de temps)

  ! kn : diffusivite turbulente des scalaires (au bas de chaque couche)

  ! (en sortie : la valeur a la fin du pas de temps)

  !

  !.......................................................................

  !=======================================================================

  ! Declarations:

  !=======================================================================

  ! Inputs/Outputs

  !----------------

  REAL dt, g, rconst

  REAL plev(klon, klev+1), temp(klon, klev)

  REAL ustar(klon)

  REAL zlev(klon, klev+1)

  REAL zlay(klon, klev)

  REAL u(klon, klev)

  REAL v(klon, klev)

  REAL teta(klon, klev)

  REAL cd(klon)

  REAL tke(klon, klev+1)

  REAL kn(klon, klev+1)

  REAL km(klon, klev+1)

  INTEGER iflag_pbl, ngrid, nsrf

  INTEGER ni(klon)

!FC

  REAL drgpro(klon,klev)

  ! Local

  !-------

  REAL kq(klon, klev+1)

  INTEGER nlay, nlev

  INTEGER ig, jg, k

  REAL ri, zrif, zalpha, zsm, zsn

  REAL m2cstat, mcstat, kmcstat

  INTEGER iter

  REAL :: disseff

  REAL,SAVE :: rifc

  !$OMP THREADPRIVATE(rifc)

  REAL,SAVE :: seuilsm, seuilalpha

  !$OMP THREADPRIVATE(seuilsm, seuilalpha)

  REAL frif, falpha, fsm

  REAL rino(klon, klev+1), smyam(klon, klev), styam(klon, klev), &

    lyam(klon, klev), knyam(klon, klev), w2yam(klon, klev), t2yam(klon, klev)

  CHARACTER (len = 20) :: modname = 'yamada4'

  CHARACTER (len = 80) :: abort_message

  ! Fonctions utilis??es

  !--------------------

  frif(ri) = 0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156))

  falpha(ri) = 1.318*(0.2231-ri)/(0.2341-ri)

  fsm(ri) = 1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri))

! Corrections et reglages issus du travail de these d'Etienne Vignon.

    ELSE

    ENDIF

!===============================================================================

! Flags, tests et ??valuations de constantes

!===============================================================================

! On utilise ou non la routine de Holstalg Boville pour les cas tres stables

  END IF

!========================================================================

! Calcul des increments verticaux

!=========================================================================

! Attention: zlev n'est pas declare a nlev

  END DO

    END DO

  END DO

  END DO

    END DO

  END DO

  END DO

!=========================================================================

! Richardson number and stability functions

!=========================================================================

! initialize arrays:

!------------------------------------------------------------

      m2(ig, k) = ((u(ig,k)-u(ig,k-1))**2+(v(ig,k)-v(ig, &

      ELSE

      END IF

else

      ELSE

      END IF

end if

    END DO

  END DO

  !=======================================================================

  !     DIFFERENT TYPES  DE SCHEMA de  YAMADA

  !=======================================================================

  ! On commence par calculer q2 a partir de la tke

      ENDDO

  ELSE

      ENDDO

  ENDIF

! ====================================================================

! Computing the mixing length

! ====================================================================

  !--------------

  ! Yamada 2.0

  !--------------

    END DO

  !------------------

  ! Yamada 2.Fournier

  !------------------

    ! Calcul de l,  km, au pas precedent

    !....................................

      END DO

    END DO

     ! Puis on ecrit la valeur de q qui annule l'equation de m supposee en q3

     !.........................................................................

          kmcstat = 1.E+0/mcstat*(unsdz(ig,k)*kmpre(ig,k+1)*mpre(ig,k+1)+ &

            unsdz(ig,k-1)*cd(ig)*(sqrt(u(ig,3)**2+v(ig,3)**2)-mcstat/unsdzdec &

            (ig,k)-mpre(ig,k+1)/unsdzdec(ig,k+1))**2)/(unsdz(ig,k)+unsdz(ig,k &

        ELSE

          kmcstat = 1.E+0/mcstat*(unsdz(ig,k)*kmpre(ig,k+1)*mpre(ig,k+1)+ &

            unsdz(ig,k-1)*kmpre(ig,k-1)*mpre(ig,k-1))/ &

        END IF

      END DO

    END DO

    ! ------------------------

    ! Yamada 2.5 a la Didi

    !-------------------------

    ! Calcul de l, km, au pas precedent

    !....................................

        END IF

        qpre = sqrt(q2(ig,k))

        ELSE

        END IF

        ! else ! iflag_pbl=9

        ! if (aa(ig,k)*qpre.gt.0.9) then

        ! q2(ig,k)=(qpre*10.)**2

        ! else

        ! q2(ig,k)=(qpre/(1.-aa(ig,k)*qpre))**2

        ! endif

        ! endif

      END DO

    END DO

!      dissip(ig,k)=min(max(((sqrt(q2(ig,k)))**3)/(b1*l(ig,k)),1.E-12),1.E4)

     ENDDO

    ENDDO

         tkeprov= tkeprov*                           &

           &  (tkeprov+dt*(shear(ig,k)+max(0.,buoy(ig,k))))/ &

        ENDDO

       ENDDO

         ! En cas stable, on traite la flotabilite comme la

         ! dissipation, en supposant que buoy/q2^3 est constant.

         ! Puis on prend la solution exacte

        ENDDO

       ENDDO

         ! En cas stable, on traite la flotabilite comme la

         ! dissipation, en supposant que buoy/q2^3 est constant.

         ! Puis on prend la solution exacte

        ENDDO

       ENDDO

         tkeprov= tkeprov*                           &

           &  tkeprov/ &

         ! En cas stable, on traite la flotabilite comme la

         ! dissipation, en supposant que buoy/q2^3 est constant.

         ! Puis on prend la solution exacte

        ENDDO

       ENDDO

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      !! Attention, yamada4_num=5 est inexacte car néglige les termes de flottabilité

      !! en conditions instables

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!FC on ajoute la dissipation due aux arbres

! on prend en compte la tke cree par les arbres

         tkeprov= (shear(ig,k)+ &

         ! En cas stable, on traite la flotabilite comme la

         ! dissipation, en supposant que buoy/q2^3 est constant.

         ! Puis on prend la solution exacte

        ENDDO

       ENDDO

         ! En cas stable, on traite la flotabilite comme la

         ! dissipation, en supposant que dissipeff/TKE est constant.

         ! Puis on prend la solution exacte

         ! With drag and dissipation from high vegetation (EV & FC, 05/10/2020)

        ENDDO

       ENDDO

    ENDIF

      ENDDO

    ENDDO

   ELSE

!     q2(1:ngrid, k) = q2(1:ngrid, k) + dt*km(1:ngrid, k)*m2(1:ngrid, k)*(1.-rif(1:ngrid,k))

!     q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(2*l(1:ngrid,k)*b1))

    END DO

  ENDIF

  ELSE

  END IF ! Fin du cas 8

  ! ====================================================================

  ! Calcul des coefficients de melange

  ! ====================================================================

    END DO

  END DO

  !====================================================================

  ! Transport diffusif vertical de la TKE par la TKE

  !====================================================================

    ! initialize near-surface and top-layer mixing coefficients

    !...........................................................

    ! Transport diffusif vertical de la TKE.

    !.......................................

  END IF

  !====================================================================

  ! Traitement particulier pour les cas tres stables, introduction d'une

  ! longueur de m??lange minimale

  !====================================================================

  !

  ! Reference: Local versus Nonlocal boundary-layer diffusion in a global climate model

  !            Holtslag A.A.M. and Boville B.A.

  !            J. Clim., 6, 1825-1842, 1993

  END IF

  END DO

  IF (1==1) THEN

          ELSE

            kmin = -1. ! kmin n'est utilise que pour les SL stables.

          END IF

 ! la longueur de melange est suposee etre l= kap z

 ! K=l q Sm d'ou q2=(K/l Sm)**2

          END IF

        END DO

      END DO

    ELSE

          ELSE

            kmin = -1. ! kmin n'est utilise que pour les SL stables.

          END IF

            kn(ig, k) = kmin

            km(ig, k) = kmin

            kq(ig, k) = kmin

 ! la longueur de melange est suposee etre l= kap z

 ! K=l q Sm d'ou q2=(K/l Sm)**2

          END IF

        END DO

      END DO

    END IF

  END IF

 END IF ! hboville

! Ajout d'une viscosite moleculaire

  END IF !(prt_level>1) THEN

 !======================================================

 ! Estimations de w'2 et T'2 d'apres Abdela et McFarlane

 !======================================================

 !

 ! Reference: A New Second-Order Turbulence Closure Scheme for the Planetary Boundary Layer

 !            Abdella K and McFarlane N

 !            J. Atmos. Sci., 54, 1850-1867, 1997

  ! Diagnostique pour stokage

  !..........................

  IF (1==0) THEN

    rino = rif

    smyam(1:ngrid, 1) = 0.

    styam(1:ngrid, 1) = 0.

    lyam(1:ngrid, 1) = 0.

    knyam(1:ngrid, 1) = 0.

    w2yam(1:ngrid, 1) = 0.

    t2yam(1:ngrid, 1) = 0.

    smyam(1:ngrid, 2:klev) = sm(1:ngrid, 2:klev)

    styam(1:ngrid, 2:klev) = sm(1:ngrid, 2:klev)*alpha(1:ngrid, 2:klev)

    lyam(1:ngrid, 2:klev) = l(1:ngrid, 2:klev)

    knyam(1:ngrid, 2:klev) = kn(1:ngrid, 2:klev)

  ! Calcul de w'2 et T'2

  !.......................

    w2yam(1:ngrid, 2:klev) = q2(1:ngrid, 2:klev)*0.24 + &

      lyam(1:ngrid, 2:klev)*5.17*kn(1:ngrid, 2:klev)*n2(1:ngrid, 2:klev)/ &

      sqrt(q2(1:ngrid,2:klev))

    t2yam(1:ngrid, 2:klev) = 9.1*kn(1:ngrid, 2:klev)* &

      dtetadz(1:ngrid, 2:klev)**2/sqrt(q2(1:ngrid,2:klev))* &

      lyam(1:ngrid, 2:klev)

  END IF

!============================================================================

! Mise a jour de la tke

!============================================================================

     ENDDO

  ELSE

     ENDDO

  ENDIF

!============================================================================

! Diagnostique de la dissipation et vitesse verticale

!============================================================================

! Diagnostics

    ENDDO

 ENDDO

!=============================================================================

END SUBROUTINE yamada4

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  USE dimphy, only : klev,klon

  IMPLICIT NONE

!    vdif_q2: subroutine qui calcule la diffusion de la TKE par la TKE

!             avec un schema implicite en temps avec

!             inversion d'un syst??me tridiagonal

!

!     Reference: Description of the interface with the surface and

!                the computation of the turbulet diffusion in LMDZ

!                Technical note on LMDZ

!                Dufresne, J-L, Ghattas, J. and Grandpeix, J-Y

!

!============================================================================

! Declarations

!============================================================================

  REAL plev(klon, klev+1)

  REAL temp(klon, klev)

  REAL timestep

  REAL gravity, rconst

  REAL kmy(klon, klev+1)

  REAL q2(klon, klev+1)

  INTEGER ngrid

  INTEGER i, k

!=========================================================================

! Calcul

!=========================================================================

      kstar(i, k) = 0.125*(kmy(i,k+1)+kmy(i,k))*zz*zz/ &

    END DO

  END DO

    END DO

  END DO

  END DO

      denom(i, k) = deltap(i, k) + (1.-beta(i,k+1))*kstar(i, k) + &

    END DO

  END DO

  ! Si on recalcule q2(1)

  !.......................

  IF (1==0) THEN

    DO i = 1, ngrid

      denom(i, 1) = deltap(i, 1) + (1-beta(i,2))*kstar(i, 1)

      q2(i, 1) = (q2(i,1)*deltap(i,1)+kstar(i,1)*alpha(i,2))/denom(i, 1)

    END DO

  END IF

    END DO

  END DO

END SUBROUTINE vdif_q2

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

   USE dimphy, only : klev,klon

  IMPLICIT NONE

! vdif_q2e: subroutine qui calcule la diffusion de TKE par la TKE

!           avec un schema explicite en temps

!====================================================

! Declarations

!====================================================

  REAL plev(klon, klev+1)

  REAL temp(klon, klev)

  REAL timestep

  REAL gravity, rconst

  REAL kmy(klon, klev+1)

  REAL q2(klon, klev+1)

  REAL denom(klon, klev+1), alpha(klon, klev+1), beta(klon, klev+1)

  INTEGER ngrid

  INTEGER i, k

!==================================================

! Calcul

!==================================================

      kstar(i, k) = 0.125*(kmy(i,k+1)+kmy(i,k))*zz*zz/ &

    END DO

  END DO

    END DO

  END DO

  END DO

      q2(i, k) = q2(i, k) + (kstar(i,k)*(q2(i,k+1)-q2(i, &

    END DO

  END DO

    q2(i, klev+1) = q2(i, klev+1) + (-kstar(i,klev)*(q2(i,klev+1)-q2(i, &

  END DO

END SUBROUTINE vdif_q2e

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  USE dimphy, only : klev,klon

  USE yamada_ini_mod, only : l0

  USE phys_state_var_mod, only: zstd, zsig, zmea

  USE phys_local_var_mod, only: l_mixmin, l_mix

  USE yamada_ini_mod, only : kap, kapb

 ! zstd: ecart type de la'altitud e sous-maille

 ! zmea: altitude moyenne sous maille

 ! zsig: pente moyenne de le maille

  USE geometry_mod, only: cell_area

  ! aire_cell: aire de la maille

  IMPLICIT NONE

!*************************************************************************

! Subrourine qui calcule la longueur de m??lange dans le sch??ma de turbulence

! avec la formule de Blackadar

! Calcul d'un  minimum en fonction de l'orographie sous-maille:

! L'id??e est la suivante: plus il y a de relief, plus il y a du m??lange

! induit par les circulations meso et submeso ??chelles.

!

! References: * The vertical distribution of wind and turbulent exchange in a neutral atmosphere

!               Blackadar A.K.

!               J. Geophys. Res., 64, No 8, 1962

!

!             * An evaluation of neutral and convective planetary boundary-layer parametrisations relative

!               to large eddy simulations

!               Ayotte K et al

!               Boundary Layer Meteorology, 79, 131-175, 1996

!

!

!             * Local Similarity in the Stable Boundary Layer and Mixing length Approaches: consistency of concepts

!               Van de Wiel B.J.H et al

!               Boundary-Lay Meteorol, 128, 103-166, 2008

!

!

! Histoire:

!----------

! * premi??re r??daction, Etienne et Frederic, 09/06/2016

!

! ***********************************************************************

!==================================================================

! Declarations

!==================================================================

! Inputs

!-------

 INTEGER            ni(klon)           ! indice sur la grille original (non restreinte)

 INTEGER            nsrf               ! Type de surface

 INTEGER            ngrid              ! Nombre de points concern??s sur l'horizontal

 INTEGER            iflag_pbl          ! Choix du sch??ma de turbulence

 REAL               pbl_lmixmin_alpha  ! on active ou non le calcul de la longueur de melange minimum en fonction du relief

 REAL               lmixmin            ! Minimum absolu de la longueur de m??lange

 REAL               zlay(klon, klev)   ! altitude du centre de la couche

 REAL               zlev(klon, klev+1) ! atitude de l'interface inf??rieure de la couche

 REAL               u(klon, klev)      ! vitesse du vent zonal

 REAL               v(klon, klev)      ! vitesse du vent meridional

 REAL               q2(klon, klev+1)   ! energie cin??tique turbulente

 REAL               n2(klon, klev+1)   ! frequence de Brunt-Vaisala

!In/out

!-------

! Outputs

!---------

 REAL               lmix(klon, klev+1)    ! Longueur de melange

! Local

!-------

 INTEGER  ig,jg, k

 REAL zq

 REAL fl, zzz, zl0, zq2, zn2

 REAL famorti, zzzz, zh_oro, zhlim

 REAL xcell

 REAL zstdslope(klon)

 REAL lmax

 REAL l2strat, l2neutre, extent

!===============================================================

! Fonctions utiles

!===============================================================

! Calcul de l suivant la formule de Blackadar 1962 adapt??e par Ayotte 1996

!..........................................................................

 fl(zzz, zl0, zq2, zn2) = max(min(l0(ig)*kap*zlev(ig, &

    k)/(kap*zlev(ig,k)+l0(ig)),0.5*sqrt(q2(ig,k))/sqrt( &

    max(n2(ig,k),1.E-10))), 1.E-5)

! Fonction d'amortissement de la turbulence au dessus de la montagne

! On augmente l'amortissement en diminuant la valeur de hlim (extent) dans le code

!.....................................................................

 famorti(zzzz, zh_oro, zhlim)=(-1.)*ATAN((zzzz-zh_oro)/(zhlim-zh_oro))*2./3.1416+1.