Changeset 3774 for LMDZ6/trunk


Ignore:
Timestamp:
Oct 1, 2020, 5:25:34 PM (4 years ago)
Author:
lguez
Message:

Bug fix: encapsulate hbtm in a module

Fixes bug introduced in commit [3772]: therm was declared as an
assumed-shape dummy argument, which requires an explicit interface.

Location:
LMDZ6/trunk/libf/phylmd
Files:
1 edited
1 moved

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/hbtm_mod.F90

    r3773 r3774  
    1 
    2 ! $Header$
    3 
    4 
    5 SUBROUTINE hbtm(knon, paprs, pplay, t2m, t10m, q2m, q10m, ustar, wstar, &
    6     flux_t, flux_q, u, v, t, q, pblh, cape, eauliq, ctei, pblt, therm, trmb1, &
    7     trmb2, trmb3, plcl)
    8   USE dimphy
     1module hbtm_mod
     2
    93  IMPLICIT NONE
    104
    11   ! ***************************************************************
    12   ! *                                                             *
    13   ! * HBTM2   D'apres Holstag&Boville et Troen&Mahrt              *
    14   ! *                 JAS 47              BLM                     *
    15   ! * Algorithme These Anne Mathieu                               *
    16   ! * Critere d'Entrainement Peter Duynkerke (JAS 50)             *
    17   ! * written by  : Anne MATHIEU & Alain LAHELLEC, 22/11/99       *
    18   ! * features : implem. exces Mathieu                            *
    19   ! ***************************************************************
    20   ! * mods : decembre 99 passage th a niveau plus bas. voir fixer *
    21   ! * la prise du th a z/Lambda = -.2 (max Ray)                   *
    22   ! * Autre algo : entrainement ~ Theta+v =cste mais comment=>The?*
    23   ! * on peut fixer q a .7qsat(cf non adiab)=>T2 et The2          *
    24   ! * voir aussi //KE pblh = niveau The_e ou l = env.             *
    25   ! ***************************************************************
    26   ! * fin therm a la HBTM passage a forme Mathieu 12/09/2001      *
    27   ! ***************************************************************
    28   ! *
    29 
    30 
    31   ! AM Fev 2003
    32   ! Adaptation a LMDZ version couplee
    33 
    34   ! Pour le moment on fait passer en argument les grdeurs de surface :
    35   ! flux, t,q2m, t,q10m, on va utiliser systematiquement les grdeurs a 2m ms
    36   ! on garde la possibilite de changer si besoin est (jusqu'a present la
    37   ! forme de HB avec le 1er niveau modele etait conservee)
    38 
    39 
    40 
    41 
    42 
    43   include "YOMCST.h"
    44   REAL rlvcp, reps
    45   ! Arguments:
    46 
    47   INTEGER knon ! nombre de points a calculer
    48   ! AM
    49   REAL t2m(klon), t10m(klon) ! temperature a 2 et 10m
    50   REAL q2m(klon), q10m(klon) ! q a 2 et 10m
    51   REAL ustar(klon)
    52   REAL wstar(klon) ! w*, convective velocity scale
    53   REAL paprs(klon, klev+1) ! pression a inter-couche (Pa)
    54   REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
    55   REAL flux_t(klon, klev), flux_q(klon, klev) ! Flux
    56   REAL u(klon, klev) ! vitesse U (m/s)
    57   REAL v(klon, klev) ! vitesse V (m/s)
    58   REAL t(klon, klev) ! temperature (K)
    59   REAL q(klon, klev) ! vapeur d'eau (kg/kg)
    60   ! AM      REAL cd_h(klon) ! coefficient de friction au sol pour chaleur
    61   ! AM      REAL cd_m(klon) ! coefficient de friction au sol pour vitesse
    62 
    63   INTEGER isommet
    64   ! um      PARAMETER (isommet=klev) ! limite max sommet pbl
    65   REAL, PARAMETER :: vk = 0.35 ! Von Karman => passer a .41 ! cf U.Olgstrom
    66   REAL, PARAMETER :: ricr = 0.4
    67   REAL, PARAMETER :: fak = 8.5 ! b calcul du Prandtl et de dTetas
    68   REAL, PARAMETER :: fakn = 7.2 ! a
    69   REAL, PARAMETER :: onet = 1.0/3.0
    70   REAL, PARAMETER :: t_coup = 273.15
    71   REAL, PARAMETER :: zkmin = 0.01
    72   REAL, PARAMETER :: betam = 15.0 ! pour Phim / h dans la S.L stable
    73   REAL, PARAMETER :: betah = 15.0
    74   REAL, PARAMETER :: betas = 5.0 ! Phit dans la S.L. stable (mais 2 formes / z/OBL<>1
    75   REAL, PARAMETER :: sffrac = 0.1 ! S.L. = z/h < .1
    76   REAL, PARAMETER :: usmin = 1.E-12
    77   REAL, PARAMETER :: binm = betam*sffrac
    78   REAL, PARAMETER :: binh = betah*sffrac
    79   REAL, PARAMETER :: ccon = fak*sffrac*vk
    80   REAL, PARAMETER :: b1 = 70., b2 = 20.
    81   REAL, PARAMETER :: zref = 2. ! Niveau de ref a 2m peut eventuellement
    82   ! etre choisi a 10m
    83   REAL q_star, t_star
    84   REAL b212, b2sr ! Lambert correlations T' q' avec T* q*
    85 
    86   REAL z(klon, klev)
    87   ! AM      REAL pcfm(klon,klev), pcfh(klon,klev)
    88   INTEGER i, k, j
    89   REAL zxt
    90   ! AM      REAL zxt, zxq, zxu, zxv, zxmod, taux, tauy
    91   ! AM      REAL zx_alf1, zx_alf2 ! parametres pour extrapolation
    92   REAL khfs(klon) ! surface kinematic heat flux [mK/s]
    93   REAL kqfs(klon) ! sfc kinematic constituent flux [m/s]
    94   REAL heatv(klon) ! surface virtual heat flux
    95   REAL rhino(klon, klev) ! bulk Richardon no. mais en Theta_v
    96   LOGICAL unstbl(klon) ! pts w/unstbl pbl (positive virtual ht flx)
    97   LOGICAL stblev(klon) ! stable pbl with levels within pbl
    98   LOGICAL unslev(klon) ! unstbl pbl with levels within pbl
    99   LOGICAL unssrf(klon) ! unstb pbl w/lvls within srf pbl lyr
    100   LOGICAL unsout(klon) ! unstb pbl w/lvls in outer pbl lyr
    101   LOGICAL check(klon) ! True=>chk if Richardson no.>critcal
    102   LOGICAL omegafl(klon) ! flag de prolongerment cape pour pt Omega
    103   REAL pblh(klon)
    104   REAL pblt(klon)
    105   REAL plcl(klon)
    106   ! AM      REAL cgh(klon,2:klev) ! counter-gradient term for heat [K/m]
    107   ! AM      REAL cgq(klon,2:klev) ! counter-gradient term for constituents
    108   ! AM      REAL cgs(klon,2:klev) ! counter-gradient star (cg/flux)
    109   REAL unsobklen(klon) ! Monin-Obukhov lengh
    110   ! AM      REAL ztvd, ztvu,
    111   REAL zdu2
    112   REAL, intent(out):: therm(:) ! (klon) thermal virtual temperature excess
    113   REAL trmb1(klon), trmb2(klon), trmb3(klon)
    114   ! Algorithme thermique
    115   REAL s(klon, klev) ! [P/Po]^Kappa milieux couches
    116   REAL th_th(klon) ! potential temperature of thermal
    117   REAL the_th(klon) ! equivalent potential temperature of thermal
    118   REAL qt_th(klon) ! total water  of thermal
    119   REAL tbef(klon) ! T thermique niveau precedent
    120   REAL qsatbef(klon)
    121   LOGICAL zsat(klon) ! le thermique est sature
    122   REAL cape(klon) ! Cape du thermique
    123   REAL kape(klon) ! Cape locale
    124   REAL eauliq(klon) ! Eau liqu integr du thermique
    125   REAL ctei(klon) ! Critere d'instab d'entrainmt des nuages de CL
    126   REAL the1, the2, aa, bb, zthvd, zthvu, xintpos, qqsat
    127   ! IM 091204 BEG
    128   REAL a1, a2, a3
    129   ! IM 091204 END
    130   REAL xhis, rnum, denom, th1, th2, thv1, thv2, ql2
    131   REAL dqsat_dt, qsat2, qt1, q2, t1, t2, xnull, delt_the
    132   REAL delt_qt, delt_2, quadsat, spblh, reduc
    133 
    134   REAL phiminv(klon) ! inverse phi function for momentum
    135   REAL phihinv(klon) ! inverse phi function for heat
    136   REAL wm(klon) ! turbulent velocity scale for momentum
    137   REAL fak1(klon) ! k*ustar*pblh
    138   REAL fak2(klon) ! k*wm*pblh
    139   REAL fak3(klon) ! fakn*wstar/wm
    140   REAL pblk(klon) ! level eddy diffusivity for momentum
    141   REAL pr(klon) ! Prandtl number for eddy diffusivities
    142   REAL zl(klon) ! zmzp / Obukhov length
    143   REAL zh(klon) ! zmzp / pblh
    144   REAL zzh(klon) ! (1-(zmzp/pblh))**2
    145   REAL zm(klon) ! current level height
    146   REAL zp(klon) ! current level height + one level up
    147   REAL zcor, zdelta, zcvm5
    148   ! AM      REAL zxqs
    149   REAL fac, pblmin, zmzp, term
    150 
    151   include "YOETHF.h"
    152   include "FCTTRE.h"
    153 
    154 
    155 
    156   ! initialisations (Anne)
    157   isommet = klev
    158   th_th(:) = 0.
    159   q_star = 0
    160   t_star = 0
    161   therm = 0.
    162 
    163   b212 = sqrt(b1*b2)
    164   b2sr = sqrt(b2)
    165 
    166   ! ============================================================
    167   ! Fonctions thermo implicites
    168   ! ============================================================
    169   ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    170   ! Tetens : pression partielle de vap d'eau e_sat(T)
    171   ! =================================================
    172   ! ++ e_sat(T) = r2*exp( r3*(T-Tf)/(T-r4) ) id a r2*FOEWE
    173   ! ++ avec :
    174   ! ++ Tf = 273.16 K  (Temp de fusion de la glace)
    175   ! ++ r2 = 611.14 Pa
    176   ! ++ r3 = 17.269 (liquide) 21.875 (solide) adim
    177   ! ++ r4 = 35.86             7.66           Kelvin
    178   ! ++  q_sat = eps*e_sat/(p-(1-eps)*e_sat)
    179   ! ++ deriv� :
    180   ! ++ =========
    181   ! ++                   r3*(Tf-r4)*q_sat(T,p)
    182   ! ++ d_qsat_dT = --------------------------------
    183   ! ++             (T-r4)^2*( 1-(1-eps)*e_sat(T)/p )
    184   ! ++ pour zcvm5=Lv, c'est FOEDE
    185   ! ++ Rq :(1.-REPS)*esarg/Parg id a RETV*Qsat
    186   ! ------------------------------------------------------------------
    187 
    188   ! Initialisation
    189   rlvcp = rlvtt/rcpd
    190   reps = rd/rv
    191 
    192 
    193   ! DO i = 1, klon
    194   ! pcfh(i,1) = cd_h(i)
    195   ! pcfm(i,1) = cd_m(i)
    196   ! ENDDO
    197   ! DO k = 2, klev
    198   ! DO i = 1, klon
    199   ! pcfh(i,k) = zkmin
    200   ! pcfm(i,k) = zkmin
    201   ! cgs(i,k) = 0.0
    202   ! cgh(i,k) = 0.0
    203   ! cgq(i,k) = 0.0
    204   ! ENDDO
    205   ! ENDDO
    206 
    207   ! Calculer les hauteurs de chaque couche
    208   ! (geopotentielle Int_dp/ro = Int_[Rd.T.dp/p] z = geop/g)
    209   ! pourquoi ne pas utiliser Phi/RG ?
    210   DO i = 1, knon
    211     z(i, 1) = rd*t(i, 1)/(0.5*(paprs(i,1)+pplay(i,1)))*(paprs(i,1)-pplay(i,1) &
    212       )/rg
    213     s(i, 1) = (pplay(i,1)/paprs(i,1))**rkappa
    214   END DO
    215   ! s(k) = [pplay(k)/ps]^kappa
    216   ! + + + + + + + + + pplay  <-> s(k)   t  dp=pplay(k-1)-pplay(k)
    217 
    218   ! -----------------  paprs <-> sig(k)
    219 
    220   ! + + + + + + + + + pplay  <-> s(k-1)
    221 
    222 
    223   ! + + + + + + + + + pplay  <-> s(1)   t  dp=paprs-pplay   z(1)
    224 
    225   ! -----------------  paprs <-> sig(1)
    226 
    227   DO k = 2, klev
    228     DO i = 1, knon
    229       z(i, k) = z(i, k-1) + rd*0.5*(t(i,k-1)+t(i,k))/paprs(i, k)*(pplay(i,k-1 &
    230         )-pplay(i,k))/rg
    231       s(i, k) = (pplay(i,k)/paprs(i,1))**rkappa
    232     END DO
    233   END DO
    234   ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    235   ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    236   ! +++  Determination des grandeurs de surface  +++++++++++++++++++++
    237   ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    238   ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    239   DO i = 1, knon
    240     ! AM         IF (thermcep) THEN
    241     ! AM           zdelta=MAX(0.,SIGN(1.,RTT-tsol(i)))
    242     ! zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
    243     ! zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q(i,1))
    244     ! AM           zxqs= r2es * FOEEW(tsol(i),zdelta)/paprs(i,1)
    245     ! AM           zxqs=MIN(0.5,zxqs)
    246     ! AM           zcor=1./(1.-retv*zxqs)
    247     ! AM           zxqs=zxqs*zcor
    248     ! AM         ELSE
    249     ! AM           IF (tsol(i).LT.t_coup) THEN
    250     ! AM              zxqs = qsats(tsol(i)) / paprs(i,1)
    251     ! AM           ELSE
    252     ! AM              zxqs = qsatl(tsol(i)) / paprs(i,1)
    253     ! AM           ENDIF
    254     ! AM         ENDIF
    255     ! niveau de reference bulk; mais ici, c,a pourrait etre le niveau de ref
    256     ! du thermique
    257     ! AM        zx_alf1 = 1.0
    258     ! AM        zx_alf2 = 1.0 - zx_alf1
    259     ! AM        zxt = (t(i,1)+z(i,1)*RG/RCPD/(1.+RVTMP2*q(i,1)))
    260     ! AM     .        *(1.+RETV*q(i,1))*zx_alf1
    261     ! AM     .      + (t(i,2)+z(i,2)*RG/RCPD/(1.+RVTMP2*q(i,2)))
    262     ! AM     .        *(1.+RETV*q(i,2))*zx_alf2
    263     ! AM        zxu = u(i,1)*zx_alf1+u(i,2)*zx_alf2
    264     ! AM        zxv = v(i,1)*zx_alf1+v(i,2)*zx_alf2
    265     ! AM        zxq = q(i,1)*zx_alf1+q(i,2)*zx_alf2
     5contains
     6
     7  SUBROUTINE hbtm(knon, paprs, pplay, t2m, t10m, q2m, q10m, ustar, wstar, &
     8       flux_t, flux_q, u, v, t, q, pblh, cape, eauliq, ctei, pblt, therm, &
     9       trmb1, trmb2, trmb3, plcl)
     10    USE dimphy
     11
     12    ! ***************************************************************
     13    ! *                                                             *
     14    ! * HBTM2   D'apres Holstag&Boville et Troen&Mahrt              *
     15    ! *                 JAS 47              BLM                     *
     16    ! * Algorithme These Anne Mathieu                               *
     17    ! * Critere d'Entrainement Peter Duynkerke (JAS 50)             *
     18    ! * written by  : Anne MATHIEU & Alain LAHELLEC, 22/11/99       *
     19    ! * features : implem. exces Mathieu                            *
     20    ! ***************************************************************
     21    ! * mods : decembre 99 passage th a niveau plus bas. voir fixer *
     22    ! * la prise du th a z/Lambda = -.2 (max Ray)                   *
     23    ! * Autre algo : entrainement ~ Theta+v =cste mais comment=>The?*
     24    ! * on peut fixer q a .7qsat(cf non adiab)=>T2 et The2          *
     25    ! * voir aussi //KE pblh = niveau The_e ou l = env.             *
     26    ! ***************************************************************
     27    ! * fin therm a la HBTM passage a forme Mathieu 12/09/2001      *
     28    ! ***************************************************************
     29    ! *
     30
     31
     32    ! AM Fev 2003
     33    ! Adaptation a LMDZ version couplee
     34
     35    ! Pour le moment on fait passer en argument les grdeurs de surface :
     36    ! flux, t,q2m, t,q10m, on va utiliser systematiquement les grdeurs a 2m ms
     37    ! on garde la possibilite de changer si besoin est (jusqu'a present la
     38    ! forme de HB avec le 1er niveau modele etait conservee)
     39
     40
     41
     42
     43
     44    include "YOMCST.h"
     45    REAL rlvcp, reps
     46    ! Arguments:
     47
     48    INTEGER knon ! nombre de points a calculer
    26649    ! AM
    267     ! AMAM           zxu = u10m(i)
    268     ! AMAM           zxv = v10m(i)
    269     ! AMAM           zxmod = 1.0+SQRT(zxu**2+zxv**2)
    270     ! AM Niveau de ref choisi a 2m
    271     zxt = t2m(i)
    272 
    273     ! ***************************************************
    274     ! attention, il doit s'agir de <w'theta'>
    275     ! ;Calcul de tcls virtuel et de w'theta'virtuel
    276     ! ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    277     ! tcls=tcls*(1+.608*qcls)
    278 
    279     ! ;Pour avoir w'theta',
    280     ! ; il faut diviser par ro.Cp
    281     ! Cp=Cpd*(1+0.84*qcls)
    282     ! fcs=fcs/(ro_surf*Cp)
    283     ! ;On transforme w'theta' en w'thetav'
    284     ! Lv=(2.501-0.00237*(tcls-273.15))*1.E6
    285     ! xle=xle/(ro_surf*Lv)
    286     ! fcsv=fcs+.608*xle*tcls
    287     ! ***************************************************
    288     ! AM        khfs(i) = (tsol(i)*(1.+RETV*q(i,1))-zxt) *zxmod*cd_h(i)
    289     ! AM        kqfs(i) = (zxqs-zxq) *zxmod*cd_h(i) * beta(i)
    290     ! AM
    291     ! dif khfs est deja w't'_v / heatv(i) = khfs(i) + RETV*zxt*kqfs(i)
    292     ! AM calcule de Ro = paprs(i,1)/Rd zxt
    293     ! AM convention >0 vers le bas ds lmdz
    294     khfs(i) = -flux_t(i, 1)*zxt*rd/(rcpd*paprs(i,1))
    295     kqfs(i) = -flux_q(i, 1)*zxt*rd/(paprs(i,1))
    296     ! AM   verifier que khfs et kqfs sont bien de la forme w'l'
    297     heatv(i) = khfs(i) + 0.608*zxt*kqfs(i)
    298     ! a comparer aussi aux sorties de clqh : flux_T/RoCp et flux_q/RoLv
    299     ! AM        heatv(i) = khfs(i)
    300     ! AM ustar est en entree
    301     ! AM        taux = zxu *zxmod*cd_m(i)
    302     ! AM        tauy = zxv *zxmod*cd_m(i)
    303     ! AM        ustar(i) = SQRT(taux**2+tauy**2)
    304     ! AM        ustar(i) = MAX(SQRT(ustar(i)),0.01)
    305     ! Theta et qT du thermique sans exces (interpolin vers surf)
    306     ! chgt de niveau du thermique (jeudi 30/12/1999)
    307     ! (interpolation lineaire avant integration phi_h)
    308     ! AM        qT_th(i) = zxqs*beta(i) + 4./z(i,1)*(q(i,1)-zxqs*beta(i))
    309     ! AM        qT_th(i) = max(qT_th(i),q(i,1))
    310     qt_th(i) = q2m(i)
    311     ! n The_th restera la Theta du thermique sans exces jusqu'a 2eme calcul
    312     ! n reste a regler convention P) pour Theta
    313     ! The_th(i) = tsol(i) + 4./z(i,1)*(t(i,1)-tsol(i))
    314     ! -                      + RLvCp*qT_th(i)
    315     ! AM        Th_th(i) = tsol(i) + 4./z(i,1)*(t(i,1)-tsol(i))
    316     th_th(i) = t2m(i)
    317   END DO
    318 
    319   DO i = 1, knon
    320     rhino(i, 1) = 0.0 ! Global Richardson
    321     check(i) = .TRUE.
    322     pblh(i) = z(i, 1) ! on initialise pblh a l'altitude du 1er niveau
    323     plcl(i) = 6000.
    324     ! Lambda = -u*^3 / (alpha.g.kvon.<w'Theta'v>
    325     unsobklen(i) = -rg*vk*heatv(i)/(t(i,1)*max(ustar(i),usmin)**3)
    326     trmb1(i) = 0.
    327     trmb2(i) = 0.
    328     trmb3(i) = 0.
    329   END DO
    330 
    331 
    332   ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    333   ! PBL height calculation:
    334   ! Search for level of pbl. Scan upward until the Richardson number between
    335   ! the first level and the current level exceeds the "critical" value.
    336   ! (bonne idee Nu de separer le Ric et l'exces de temp du thermique)
    337   ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    338   fac = 100.0
    339   DO k = 2, isommet
    340     DO i = 1, knon
    341       IF (check(i)) THEN
    342         ! pourquoi / niveau 1 (au lieu du sol) et le terme en u*^2 ?
    343         ! test     zdu2 =
    344         ! (u(i,k)-u(i,1))**2+(v(i,k)-v(i,1))**2+fac*ustar(i)**2
    345         zdu2 = u(i, k)**2 + v(i, k)**2
    346         zdu2 = max(zdu2, 1.0E-20)
    347         ! Theta_v environnement
    348         zthvd = t(i, k)/s(i, k)*(1.+retv*q(i,k))
    349 
    350         ! therm Theta_v sans exces (avec hypothese fausse de H&B, sinon,
    351         ! passer par Theta_e et virpot)
    352         ! zthvu=t(i,1)/s(i,1)*(1.+RETV*q(i,1))
    353         ! AM         zthvu = Th_th(i)*(1.+RETV*q(i,1))
    354         zthvu = th_th(i)*(1.+retv*qt_th(i))
    355         ! Le Ri par Theta_v
    356         ! AM         rhino(i,k) = (z(i,k)-z(i,1))*RG*(zthvd-zthvu)
    357         ! AM     .               /(zdu2*0.5*(zthvd+zthvu))
    358         ! AM On a nveau de ref a 2m ???
    359         rhino(i, k) = (z(i,k)-zref)*rg*(zthvd-zthvu)/(zdu2*0.5*(zthvd+zthvu))
    360 
    361         IF (rhino(i,k)>=ricr) THEN
    362           pblh(i) = z(i, k-1) + (z(i,k-1)-z(i,k))*(ricr-rhino(i,k-1))/(rhino( &
    363             i,k-1)-rhino(i,k))
    364           ! test04
    365           pblh(i) = pblh(i) + 100.
    366           pblt(i) = t(i, k-1) + (t(i,k)-t(i,k-1))*(pblh(i)-z(i,k-1))/(z(i,k)- &
    367             z(i,k-1))
     50    REAL t2m(klon), t10m(klon) ! temperature a 2 et 10m
     51    REAL q2m(klon), q10m(klon) ! q a 2 et 10m
     52    REAL ustar(klon)
     53    REAL wstar(klon) ! w*, convective velocity scale
     54    REAL paprs(klon, klev+1) ! pression a inter-couche (Pa)
     55    REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
     56    REAL flux_t(klon, klev), flux_q(klon, klev) ! Flux
     57    REAL u(klon, klev) ! vitesse U (m/s)
     58    REAL v(klon, klev) ! vitesse V (m/s)
     59    REAL t(klon, klev) ! temperature (K)
     60    REAL q(klon, klev) ! vapeur d'eau (kg/kg)
     61    ! AM      REAL cd_h(klon) ! coefficient de friction au sol pour chaleur
     62    ! AM      REAL cd_m(klon) ! coefficient de friction au sol pour vitesse
     63
     64    INTEGER isommet
     65    ! um      PARAMETER (isommet=klev) ! limite max sommet pbl
     66    REAL, PARAMETER :: vk = 0.35 ! Von Karman => passer a .41 ! cf U.Olgstrom
     67    REAL, PARAMETER :: ricr = 0.4
     68    REAL, PARAMETER :: fak = 8.5 ! b calcul du Prandtl et de dTetas
     69    REAL, PARAMETER :: fakn = 7.2 ! a
     70    REAL, PARAMETER :: onet = 1.0/3.0
     71    REAL, PARAMETER :: t_coup = 273.15
     72    REAL, PARAMETER :: zkmin = 0.01
     73    REAL, PARAMETER :: betam = 15.0 ! pour Phim / h dans la S.L stable
     74    REAL, PARAMETER :: betah = 15.0
     75
     76    REAL, PARAMETER :: betas = 5.0
     77    ! Phit dans la S.L. stable (mais 2 formes / z/OBL<>1
     78
     79    REAL, PARAMETER :: sffrac = 0.1 ! S.L. = z/h < .1
     80    REAL, PARAMETER :: usmin = 1.E-12
     81    REAL, PARAMETER :: binm = betam*sffrac
     82    REAL, PARAMETER :: binh = betah*sffrac
     83    REAL, PARAMETER :: ccon = fak*sffrac*vk
     84    REAL, PARAMETER :: b1 = 70., b2 = 20.
     85    REAL, PARAMETER :: zref = 2. ! Niveau de ref a 2m peut eventuellement
     86    ! etre choisi a 10m
     87    REAL q_star, t_star
     88    REAL b212, b2sr ! Lambert correlations T' q' avec T* q*
     89
     90    REAL z(klon, klev)
     91    ! AM      REAL pcfm(klon,klev), pcfh(klon,klev)
     92    INTEGER i, k, j
     93    REAL zxt
     94    ! AM      REAL zxt, zxq, zxu, zxv, zxmod, taux, tauy
     95    ! AM      REAL zx_alf1, zx_alf2 ! parametres pour extrapolation
     96    REAL khfs(klon) ! surface kinematic heat flux [mK/s]
     97    REAL kqfs(klon) ! sfc kinematic constituent flux [m/s]
     98    REAL heatv(klon) ! surface virtual heat flux
     99    REAL rhino(klon, klev) ! bulk Richardon no. mais en Theta_v
     100    LOGICAL unstbl(klon) ! pts w/unstbl pbl (positive virtual ht flx)
     101    LOGICAL stblev(klon) ! stable pbl with levels within pbl
     102    LOGICAL unslev(klon) ! unstbl pbl with levels within pbl
     103    LOGICAL unssrf(klon) ! unstb pbl w/lvls within srf pbl lyr
     104    LOGICAL unsout(klon) ! unstb pbl w/lvls in outer pbl lyr
     105    LOGICAL check(klon) ! True=>chk if Richardson no.>critcal
     106    LOGICAL omegafl(klon) ! flag de prolongerment cape pour pt Omega
     107    REAL pblh(klon)
     108    REAL pblt(klon)
     109    REAL plcl(klon)
     110    ! AM      REAL cgh(klon,2:klev) ! counter-gradient term for heat [K/m]
     111    ! AM      REAL cgq(klon,2:klev) ! counter-gradient term for constituents
     112    ! AM      REAL cgs(klon,2:klev) ! counter-gradient star (cg/flux)
     113    REAL unsobklen(klon) ! Monin-Obukhov lengh
     114    ! AM      REAL ztvd, ztvu,
     115    REAL zdu2
     116    REAL, intent(out):: therm(:) ! (klon) thermal virtual temperature excess
     117    REAL trmb1(klon), trmb2(klon), trmb3(klon)
     118    ! Algorithme thermique
     119    REAL s(klon, klev) ! [P/Po]^Kappa milieux couches
     120    REAL th_th(klon) ! potential temperature of thermal
     121    REAL the_th(klon) ! equivalent potential temperature of thermal
     122    REAL qt_th(klon) ! total water  of thermal
     123    REAL tbef(klon) ! T thermique niveau precedent
     124    REAL qsatbef(klon)
     125    LOGICAL zsat(klon) ! le thermique est sature
     126    REAL cape(klon) ! Cape du thermique
     127    REAL kape(klon) ! Cape locale
     128    REAL eauliq(klon) ! Eau liqu integr du thermique
     129    REAL ctei(klon) ! Critere d'instab d'entrainmt des nuages de CL
     130    REAL the1, the2, aa, bb, zthvd, zthvu, xintpos, qqsat
     131    ! IM 091204 BEG
     132    REAL a1, a2, a3
     133    ! IM 091204 END
     134    REAL xhis, rnum, denom, th1, th2, thv1, thv2, ql2
     135    REAL dqsat_dt, qsat2, qt1, q2, t1, t2, xnull, delt_the
     136    REAL delt_qt, delt_2, quadsat, spblh, reduc
     137
     138    REAL phiminv(klon) ! inverse phi function for momentum
     139    REAL phihinv(klon) ! inverse phi function for heat
     140    REAL wm(klon) ! turbulent velocity scale for momentum
     141    REAL fak1(klon) ! k*ustar*pblh
     142    REAL fak2(klon) ! k*wm*pblh
     143    REAL fak3(klon) ! fakn*wstar/wm
     144    REAL pblk(klon) ! level eddy diffusivity for momentum
     145    REAL pr(klon) ! Prandtl number for eddy diffusivities
     146    REAL zl(klon) ! zmzp / Obukhov length
     147    REAL zh(klon) ! zmzp / pblh
     148    REAL zzh(klon) ! (1-(zmzp/pblh))**2
     149    REAL zm(klon) ! current level height
     150    REAL zp(klon) ! current level height + one level up
     151    REAL zcor, zdelta, zcvm5
     152    ! AM      REAL zxqs
     153    REAL fac, pblmin, zmzp, term
     154
     155    include "YOETHF.h"
     156    include "FCTTRE.h"
     157
     158
     159
     160    ! initialisations (Anne)
     161    isommet = klev
     162    th_th(:) = 0.
     163    q_star = 0
     164    t_star = 0
     165    therm = 0.
     166
     167    b212 = sqrt(b1*b2)
     168    b2sr = sqrt(b2)
     169
     170    ! ============================================================
     171    ! Fonctions thermo implicites
     172    ! ============================================================
     173    ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     174    ! Tetens : pression partielle de vap d'eau e_sat(T)
     175    ! =================================================
     176    ! ++ e_sat(T) = r2*exp( r3*(T-Tf)/(T-r4) ) id a r2*FOEWE
     177    ! ++ avec :
     178    ! ++ Tf = 273.16 K  (Temp de fusion de la glace)
     179    ! ++ r2 = 611.14 Pa
     180    ! ++ r3 = 17.269 (liquide) 21.875 (solide) adim
     181    ! ++ r4 = 35.86             7.66           Kelvin
     182    ! ++  q_sat = eps*e_sat/(p-(1-eps)*e_sat)
     183    ! ++ deriv� :
     184    ! ++ =========
     185    ! ++                   r3*(Tf-r4)*q_sat(T,p)
     186    ! ++ d_qsat_dT = --------------------------------
     187    ! ++             (T-r4)^2*( 1-(1-eps)*e_sat(T)/p )
     188    ! ++ pour zcvm5=Lv, c'est FOEDE
     189    ! ++ Rq :(1.-REPS)*esarg/Parg id a RETV*Qsat
     190    ! ------------------------------------------------------------------
     191
     192    ! Initialisation
     193    rlvcp = rlvtt/rcpd
     194    reps = rd/rv
     195
     196
     197    ! DO i = 1, klon
     198    ! pcfh(i,1) = cd_h(i)
     199    ! pcfm(i,1) = cd_m(i)
     200    ! ENDDO
     201    ! DO k = 2, klev
     202    ! DO i = 1, klon
     203    ! pcfh(i,k) = zkmin
     204    ! pcfm(i,k) = zkmin
     205    ! cgs(i,k) = 0.0
     206    ! cgh(i,k) = 0.0
     207    ! cgq(i,k) = 0.0
     208    ! ENDDO
     209    ! ENDDO
     210
     211    ! Calculer les hauteurs de chaque couche
     212    ! (geopotentielle Int_dp/ro = Int_[Rd.T.dp/p] z = geop/g)
     213    ! pourquoi ne pas utiliser Phi/RG ?
     214    DO i = 1, knon
     215       z(i, 1) = rd*t(i, 1)/(0.5*(paprs(i,1)+pplay(i,1)))&
     216            *(paprs(i,1)-pplay(i,1))/rg
     217       s(i, 1) = (pplay(i,1)/paprs(i,1))**rkappa
     218    END DO
     219    ! s(k) = [pplay(k)/ps]^kappa
     220    ! + + + + + + + + + pplay  <-> s(k)   t  dp=pplay(k-1)-pplay(k)
     221
     222    ! -----------------  paprs <-> sig(k)
     223
     224    ! + + + + + + + + + pplay  <-> s(k-1)
     225
     226
     227    ! + + + + + + + + + pplay  <-> s(1)   t  dp=paprs-pplay   z(1)
     228
     229    ! -----------------  paprs <-> sig(1)
     230
     231    DO k = 2, klev
     232       DO i = 1, knon
     233          z(i, k) = z(i, k-1) + rd*0.5*(t(i,k-1)+t(i,k))/paprs(i, k)&
     234               *(pplay(i,k-1)-pplay(i,k))/rg
     235          s(i, k) = (pplay(i,k)/paprs(i,1))**rkappa
     236       END DO
     237    END DO
     238    ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     239    ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     240    ! +++  Determination des grandeurs de surface  +++++++++++++++++++++
     241    ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     242    ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     243    DO i = 1, knon
     244       ! AM         IF (thermcep) THEN
     245       ! AM           zdelta=MAX(0.,SIGN(1.,RTT-tsol(i)))
     246       ! zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
     247       ! zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q(i,1))
     248       ! AM           zxqs= r2es * FOEEW(tsol(i),zdelta)/paprs(i,1)
     249       ! AM           zxqs=MIN(0.5,zxqs)
     250       ! AM           zcor=1./(1.-retv*zxqs)
     251       ! AM           zxqs=zxqs*zcor
     252       ! AM         ELSE
     253       ! AM           IF (tsol(i).LT.t_coup) THEN
     254       ! AM              zxqs = qsats(tsol(i)) / paprs(i,1)
     255       ! AM           ELSE
     256       ! AM              zxqs = qsatl(tsol(i)) / paprs(i,1)
     257       ! AM           ENDIF
     258       ! AM         ENDIF
     259       ! niveau de reference bulk; mais ici, c,a pourrait etre le niveau de ref
     260       ! du thermique
     261       ! AM        zx_alf1 = 1.0
     262       ! AM        zx_alf2 = 1.0 - zx_alf1
     263       ! AM        zxt = (t(i,1)+z(i,1)*RG/RCPD/(1.+RVTMP2*q(i,1)))
     264       ! AM     .        *(1.+RETV*q(i,1))*zx_alf1
     265       ! AM     .      + (t(i,2)+z(i,2)*RG/RCPD/(1.+RVTMP2*q(i,2)))
     266       ! AM     .        *(1.+RETV*q(i,2))*zx_alf2
     267       ! AM        zxu = u(i,1)*zx_alf1+u(i,2)*zx_alf2
     268       ! AM        zxv = v(i,1)*zx_alf1+v(i,2)*zx_alf2
     269       ! AM        zxq = q(i,1)*zx_alf1+q(i,2)*zx_alf2
     270       ! AM
     271       ! AMAM           zxu = u10m(i)
     272       ! AMAM           zxv = v10m(i)
     273       ! AMAM           zxmod = 1.0+SQRT(zxu**2+zxv**2)
     274       ! AM Niveau de ref choisi a 2m
     275       zxt = t2m(i)
     276
     277       ! ***************************************************
     278       ! attention, il doit s'agir de <w'theta'>
     279       ! ;Calcul de tcls virtuel et de w'theta'virtuel
     280       ! ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     281       ! tcls=tcls*(1+.608*qcls)
     282
     283       ! ;Pour avoir w'theta',
     284       ! ; il faut diviser par ro.Cp
     285       ! Cp=Cpd*(1+0.84*qcls)
     286       ! fcs=fcs/(ro_surf*Cp)
     287       ! ;On transforme w'theta' en w'thetav'
     288       ! Lv=(2.501-0.00237*(tcls-273.15))*1.E6
     289       ! xle=xle/(ro_surf*Lv)
     290       ! fcsv=fcs+.608*xle*tcls
     291       ! ***************************************************
     292       ! AM        khfs(i) = (tsol(i)*(1.+RETV*q(i,1))-zxt) *zxmod*cd_h(i)
     293       ! AM        kqfs(i) = (zxqs-zxq) *zxmod*cd_h(i) * beta(i)
     294       ! AM
     295       ! dif khfs est deja w't'_v / heatv(i) = khfs(i) + RETV*zxt*kqfs(i)
     296       ! AM calcule de Ro = paprs(i,1)/Rd zxt
     297       ! AM convention >0 vers le bas ds lmdz
     298       khfs(i) = -flux_t(i, 1)*zxt*rd/(rcpd*paprs(i,1))
     299       kqfs(i) = -flux_q(i, 1)*zxt*rd/(paprs(i,1))
     300       ! AM   verifier que khfs et kqfs sont bien de la forme w'l'
     301       heatv(i) = khfs(i) + 0.608*zxt*kqfs(i)
     302       ! a comparer aussi aux sorties de clqh : flux_T/RoCp et flux_q/RoLv
     303       ! AM        heatv(i) = khfs(i)
     304       ! AM ustar est en entree
     305       ! AM        taux = zxu *zxmod*cd_m(i)
     306       ! AM        tauy = zxv *zxmod*cd_m(i)
     307       ! AM        ustar(i) = SQRT(taux**2+tauy**2)
     308       ! AM        ustar(i) = MAX(SQRT(ustar(i)),0.01)
     309       ! Theta et qT du thermique sans exces (interpolin vers surf)
     310       ! chgt de niveau du thermique (jeudi 30/12/1999)
     311       ! (interpolation lineaire avant integration phi_h)
     312       ! AM        qT_th(i) = zxqs*beta(i) + 4./z(i,1)*(q(i,1)-zxqs*beta(i))
     313       ! AM        qT_th(i) = max(qT_th(i),q(i,1))
     314       qt_th(i) = q2m(i)
     315       ! n The_th restera la Theta du thermique sans exces jusqu'a 2eme calcul
     316       ! n reste a regler convention P) pour Theta
     317       ! The_th(i) = tsol(i) + 4./z(i,1)*(t(i,1)-tsol(i))
     318       ! -                      + RLvCp*qT_th(i)
     319       ! AM        Th_th(i) = tsol(i) + 4./z(i,1)*(t(i,1)-tsol(i))
     320       th_th(i) = t2m(i)
     321    END DO
     322
     323    DO i = 1, knon
     324       rhino(i, 1) = 0.0 ! Global Richardson
     325       check(i) = .TRUE.
     326       pblh(i) = z(i, 1) ! on initialise pblh a l'altitude du 1er niveau
     327       plcl(i) = 6000.
     328       ! Lambda = -u*^3 / (alpha.g.kvon.<w'Theta'v>
     329       unsobklen(i) = -rg*vk*heatv(i)/(t(i,1)*max(ustar(i),usmin)**3)
     330       trmb1(i) = 0.
     331       trmb2(i) = 0.
     332       trmb3(i) = 0.
     333    END DO
     334
     335
     336    ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     337    ! PBL height calculation:
     338    ! Search for level of pbl. Scan upward until the Richardson number between
     339    ! the first level and the current level exceeds the "critical" value.
     340    ! (bonne idee Nu de separer le Ric et l'exces de temp du thermique)
     341    ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     342    fac = 100.0
     343    DO k = 2, isommet
     344       DO i = 1, knon
     345          IF (check(i)) THEN
     346             ! pourquoi / niveau 1 (au lieu du sol) et le terme en u*^2 ?
     347             ! test     zdu2 =
     348             ! (u(i,k)-u(i,1))**2+(v(i,k)-v(i,1))**2+fac*ustar(i)**2
     349             zdu2 = u(i, k)**2 + v(i, k)**2
     350             zdu2 = max(zdu2, 1.0E-20)
     351             ! Theta_v environnement
     352             zthvd = t(i, k)/s(i, k)*(1.+retv*q(i,k))
     353
     354             ! therm Theta_v sans exces (avec hypothese fausse de H&B, sinon,
     355             ! passer par Theta_e et virpot)
     356             ! zthvu=t(i,1)/s(i,1)*(1.+RETV*q(i,1))
     357             ! AM         zthvu = Th_th(i)*(1.+RETV*q(i,1))
     358             zthvu = th_th(i)*(1.+retv*qt_th(i))
     359             ! Le Ri par Theta_v
     360             ! AM         rhino(i,k) = (z(i,k)-z(i,1))*RG*(zthvd-zthvu)
     361             ! AM     .               /(zdu2*0.5*(zthvd+zthvu))
     362             ! AM On a nveau de ref a 2m ???
     363             rhino(i, k) = (z(i,k)-zref)*rg*(zthvd-zthvu)/(zdu2*0.5&
     364                  *(zthvd+zthvu))
     365
     366             IF (rhino(i,k)>=ricr) THEN
     367                pblh(i) = z(i, k-1) + (z(i,k-1)-z(i,k))*(ricr-rhino(i,k-1))&
     368                     /(rhino(i,k-1)-rhino(i,k))
     369                ! test04
     370                pblh(i) = pblh(i) + 100.
     371                pblt(i) = t(i, k-1) + (t(i,k)-t(i,k-1))*(pblh(i)-z(i,k-1))&
     372                     /(z(i,k)- z(i,k-1))
     373                check(i) = .FALSE.
     374             END IF
     375          END IF
     376       END DO
     377    END DO
     378
     379
     380    ! Set pbl height to maximum value where computation exceeds number of
     381    ! layers allowed
     382
     383    DO i = 1, knon
     384       IF (check(i)) pblh(i) = z(i, isommet)
     385    END DO
     386
     387    ! Improve estimate of pbl height for the unstable points.
     388    ! Find unstable points (sensible heat flux is upward):
     389
     390    DO i = 1, knon
     391       IF (heatv(i)>0.) THEN
     392          unstbl(i) = .TRUE.
     393          check(i) = .TRUE.
     394       ELSE
     395          unstbl(i) = .FALSE.
    368396          check(i) = .FALSE.
    369         END IF
    370       END IF
    371     END DO
    372   END DO
    373 
    374 
    375   ! Set pbl height to maximum value where computation exceeds number of
    376   ! layers allowed
    377 
    378   DO i = 1, knon
    379     IF (check(i)) pblh(i) = z(i, isommet)
    380   END DO
    381 
    382   ! Improve estimate of pbl height for the unstable points.
    383   ! Find unstable points (sensible heat flux is upward):
    384 
    385   DO i = 1, knon
    386     IF (heatv(i)>0.) THEN
    387       unstbl(i) = .TRUE.
    388       check(i) = .TRUE.
    389     ELSE
    390       unstbl(i) = .FALSE.
    391       check(i) = .FALSE.
    392     END IF
    393   END DO
    394 
    395   ! For the unstable case, compute velocity scale and the
    396   ! convective temperature excess:
    397 
    398   DO i = 1, knon
    399     IF (check(i)) THEN
    400       phiminv(i) = (1.-binm*pblh(i)*unsobklen(i))**onet
    401       ! ***************************************************
    402       ! Wm ? et W* ? c'est la formule pour z/h < .1
    403       ! ;Calcul de w* ;;
    404       ! ;;;;;;;;;;;;;;;;
    405       ! w_star=((g/tcls)*fcsv*z(ind))^(1/3.) [ou prendre la premiere approx
    406       ! de h)
    407       ! ;; CALCUL DE wm ;;
    408       ! ;;;;;;;;;;;;;;;;;;
    409       ! ; Ici on considerera que l'on est dans la couche de surf jusqu'a 100m
    410       ! ; On prend svt couche de surface=0.1*h mais on ne connait pas h
    411       ! ;;;;;;;;;;;Dans la couche de surface
    412       ! if (z(ind) le 20) then begin
    413       ! Phim=(1.-15.*(z(ind)/L))^(-1/3.)
    414       ! wm=u_star/Phim
    415       ! ;;;;;;;;;;;En dehors de la couche de surface
    416       ! endif else if (z(ind) gt 20) then begin
    417       ! wm=(u_star^3+c1*w_star^3)^(1/3.)
    418       ! endif
    419       ! ***************************************************
    420       wm(i) = ustar(i)*phiminv(i)
    421       ! ======================================================================
    422       ! valeurs de Dominique Lambert de la campagne SEMAPHORE :
    423       ! <T'^2> = 100.T*^2; <q'^2> = 20.q*^2 a 10m
    424       ! <Tv'^2> = (1+1.2q).100.T* + 1.2Tv.sqrt(20*100).T*.q* +
    425       ! (.608*Tv)^2*20.q*^2;
    426       ! et dTetavS = sqrt(<Tv'^2>) ainsi calculee.
    427       ! avec : T*=<w'T'>_s/w* et q*=<w'q'>/w*
    428       ! !!! on peut donc utiliser w* pour les fluctuations <-> Lambert
    429       ! (leur corellation pourrait dependre de beta par ex)
    430       ! if fcsv(i,j) gt 0 then begin
    431       ! dTetavs=b1*(1.+2.*.608*q_10(i,j))*(fcs(i,j)/wm(i,j))^2+$
    432       ! (.608*Thetav_10(i,j))^2*b2*(xle(i,j)/wm(i,j))^2+$
    433       ! 2.*.608*thetav_10(i,j)*sqrt(b1*b2)*(xle(i,j)/wm(i,j))*(fcs(i,j)/wm(i,j))
    434       ! dqs=b2*(xle(i,j)/wm(i,j))^2
    435       ! theta_s(i,j)=thetav_10(i,j)+sqrt(dTetavs)
    436       ! q_s(i,j)=q_10(i,j)+sqrt(dqs)
    437       ! endif else begin
    438       ! Theta_s(i,j)=thetav_10(i,j)
    439       ! q_s(i,j)=q_10(i,j)
    440       ! endelse
    441       ! ======================================================================
    442 
    443       ! HBTM        therm(i) = heatv(i)*fak/wm(i)
    444       ! forme Mathieu :
    445       q_star = kqfs(i)/wm(i)
    446       t_star = khfs(i)/wm(i)
    447       ! IM 091204 BEG
    448       IF (1==0) THEN
    449         IF (t_star<0. .OR. q_star<0.) THEN
    450           PRINT *, 'i t_star q_star khfs kqfs wm', i, t_star, q_star, &
    451             khfs(i), kqfs(i), wm(i)
    452         END IF
    453       END IF
    454       ! IM 091204 END
    455       ! AM Nveau cde ref 2m =>
    456       ! AM        therm(i) = sqrt( b1*(1.+2.*RETV*q(i,1))*t_star**2
    457       ! AM     +             + (RETV*T(i,1))**2*b2*q_star**2
    458       ! AM     +             + 2.*RETV*T(i,1)*b212*q_star*t_star
    459       ! AM     +                 )
    460       ! IM 091204 BEG
    461       a1 = b1*(1.+2.*retv*qt_th(i))*t_star**2
    462       a2 = (retv*th_th(i))**2*b2*q_star*q_star
    463       a3 = 2.*retv*th_th(i)*b212*q_star*t_star
    464       aa = a1 + a2 + a3
    465       IF (1==0) THEN
    466         IF (aa<0.) THEN
    467           PRINT *, 'i a1 a2 a3 aa', i, a1, a2, a3, aa
    468           PRINT *, 'i qT_th Th_th t_star q_star RETV b1 b2 b212', i, &
    469             qt_th(i), th_th(i), t_star, q_star, retv, b1, b2, b212
    470         END IF
    471       END IF
    472       ! IM 091204 END
    473       therm(i) = sqrt(b1*(1.+2.*retv*qt_th(i))*t_star**2+(retv*th_th( &
    474         i))**2*b2*q_star*q_star &  ! IM 101204  +             +
    475                                    ! 2.*RETV*Th_th(i)*b212*q_star*t_star
    476         +max(0.,2.*retv*th_th(i)*b212*q_star*t_star))
    477 
    478       ! Theta et qT du thermique (forme H&B) avec exces
    479       ! (attention, on ajoute therm(i) qui est virtuelle ...)
    480       ! pourquoi pas sqrt(b1)*t_star ?
    481       ! dqs = b2sr*kqfs(i)/wm(i)
    482       qt_th(i) = qt_th(i) + b2sr*q_star
    483       ! new on differre le calcul de Theta_e
    484       ! The_th(i) = The_th(i) + therm(i) + RLvCp*qT_th(i)
    485       ! ou:    The_th(i) = The_th(i) + sqrt(b1)*khfs(i)/wm(i) +
    486       ! RLvCp*qT_th(i)
    487       rhino(i, 1) = 0.0
    488     END IF
    489   END DO
    490 
    491   ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    492   ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    493   ! ++ Improve pblh estimate for unstable conditions using the +++++++
    494   ! ++          convective temperature excess :                +++++++
    495   ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    496   ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    497 
    498   DO k = 2, isommet
    499     DO i = 1, knon
    500       IF (check(i)) THEN
    501         ! test     zdu2 =
    502         ! (u(i,k)-u(i,1))**2+(v(i,k)-v(i,1))**2+fac*ustar(i)**2
    503         zdu2 = u(i, k)**2 + v(i, k)**2
    504         zdu2 = max(zdu2, 1.0E-20)
    505         ! Theta_v environnement
    506         zthvd = t(i, k)/s(i, k)*(1.+retv*q(i,k))
    507 
    508         ! et therm Theta_v (avec hypothese de constance de H&B,
    509         ! zthvu=(t(i,1)+therm(i))/s(i,1)*(1.+RETV*q(i,1))
    510         zthvu = th_th(i)*(1.+retv*qt_th(i)) + therm(i)
    511 
    512 
    513         ! Le Ri par Theta_v
    514         ! AM Niveau de ref 2m
    515         ! AM         rhino(i,k) = (z(i,k)-z(i,1))*RG*(zthvd-zthvu)
    516         ! AM     .               /(zdu2*0.5*(zthvd+zthvu))
    517         rhino(i, k) = (z(i,k)-zref)*rg*(zthvd-zthvu)/(zdu2*0.5*(zthvd+zthvu))
    518 
    519 
    520         IF (rhino(i,k)>=ricr) THEN
    521           pblh(i) = z(i, k-1) + (z(i,k-1)-z(i,k))*(ricr-rhino(i,k-1))/(rhino( &
    522             i,k-1)-rhino(i,k))
    523           ! test04
    524           pblh(i) = pblh(i) + 100.
    525           pblt(i) = t(i, k-1) + (t(i,k)-t(i,k-1))*(pblh(i)-z(i,k-1))/(z(i,k)- &
    526             z(i,k-1))
    527           check(i) = .FALSE.
    528           ! IM 170305 BEG
     397       END IF
     398    END DO
     399
     400    ! For the unstable case, compute velocity scale and the
     401    ! convective temperature excess:
     402
     403    DO i = 1, knon
     404       IF (check(i)) THEN
     405          phiminv(i) = (1.-binm*pblh(i)*unsobklen(i))**onet
     406          ! ***************************************************
     407          ! Wm ? et W* ? c'est la formule pour z/h < .1
     408          ! ;Calcul de w* ;;
     409          ! ;;;;;;;;;;;;;;;;
     410          ! w_star=((g/tcls)*fcsv*z(ind))^(1/3.) [ou prendre la premiere approx
     411          ! de h)
     412          ! ;; CALCUL DE wm ;;
     413          ! ;;;;;;;;;;;;;;;;;;
     414          ! ; Ici on considerera que l'on est dans la couche de surf jusqu'a
     415          ! 100 m
     416          ! ; On prend svt couche de surface=0.1*h mais on ne connait pas h
     417          ! ;;;;;;;;;;;Dans la couche de surface
     418          ! if (z(ind) le 20) then begin
     419          ! Phim=(1.-15.*(z(ind)/L))^(-1/3.)
     420          ! wm=u_star/Phim
     421          ! ;;;;;;;;;;;En dehors de la couche de surface
     422          ! endif else if (z(ind) gt 20) then begin
     423          ! wm=(u_star^3+c1*w_star^3)^(1/3.)
     424          ! endif
     425          ! ***************************************************
     426          wm(i) = ustar(i)*phiminv(i)
     427          ! ===================================================================
     428          ! valeurs de Dominique Lambert de la campagne SEMAPHORE :
     429          ! <T'^2> = 100.T*^2; <q'^2> = 20.q*^2 a 10m
     430          ! <Tv'^2> = (1+1.2q).100.T* + 1.2Tv.sqrt(20*100).T*.q* +
     431          ! (.608*Tv)^2*20.q*^2;
     432          ! et dTetavS = sqrt(<Tv'^2>) ainsi calculee.
     433          ! avec : T*=<w'T'>_s/w* et q*=<w'q'>/w*
     434          ! !!! on peut donc utiliser w* pour les fluctuations <-> Lambert
     435          ! (leur corellation pourrait dependre de beta par ex)
     436          ! if fcsv(i,j) gt 0 then begin
     437          ! dTetavs=b1*(1.+2.*.608*q_10(i,j))*(fcs(i,j)/wm(i,j))^2+$
     438          ! (.608*Thetav_10(i,j))^2*b2*(xle(i,j)/wm(i,j))^2+$
     439          ! 2.*.608*thetav_10(i,j)*sqrt(b1*b2)*(xle(i,j)/wm(i,j))*(fcs(i,j)
     440          ! /wm(i,j))
     441          ! dqs=b2*(xle(i,j)/wm(i,j))^2
     442          ! theta_s(i,j)=thetav_10(i,j)+sqrt(dTetavs)
     443          ! q_s(i,j)=q_10(i,j)+sqrt(dqs)
     444          ! endif else begin
     445          ! Theta_s(i,j)=thetav_10(i,j)
     446          ! q_s(i,j)=q_10(i,j)
     447          ! endelse
     448          ! ===================================================================
     449
     450          ! HBTM        therm(i) = heatv(i)*fak/wm(i)
     451          ! forme Mathieu :
     452          q_star = kqfs(i)/wm(i)
     453          t_star = khfs(i)/wm(i)
     454          ! IM 091204 BEG
    529455          IF (1==0) THEN
    530             ! debug print -120;34       -34-        58 et    0;26 wamp
    531             IF (i==950 .OR. i==192 .OR. i==624 .OR. i==118) THEN
    532               PRINT *, ' i,Th_th,Therm,qT :', i, th_th(i), therm(i), qt_th(i)
    533               q_star = kqfs(i)/wm(i)
    534               t_star = khfs(i)/wm(i)
    535               PRINT *, 'q* t*, b1,b2,b212 ', q_star, t_star, &
    536                 b1*(1.+2.*retv*qt_th(i))*t_star**2, &
    537                 (retv*th_th(i))**2*b2*q_star**2, 2.*retv*th_th(i)*b212*q_star &
    538                 *t_star
    539               PRINT *, 'zdu2 ,100.*ustar(i)**2', zdu2, fac*ustar(i)**2
    540             END IF
    541           END IF !(1.EQ.0) THEN
    542           ! IM 170305 END
    543           ! q_star = kqfs(i)/wm(i)
    544           ! t_star = khfs(i)/wm(i)
    545           ! trmb1(i) = b1*(1.+2.*RETV*q(i,1))*t_star**2
    546           ! trmb2(i) = (RETV*T(i,1))**2*b2*q_star**2
    547           ! Omega now   trmb3(i) = 2.*RETV*T(i,1)*b212*q_star*t_star
    548         END IF
    549       END IF
    550     END DO
    551   END DO
    552 
    553   ! Set pbl height to maximum value where computation exceeds number of
    554   ! layers allowed
    555 
    556   DO i = 1, knon
    557     IF (check(i)) pblh(i) = z(i, isommet)
    558   END DO
    559 
    560   ! PBL height must be greater than some minimum mechanical mixing depth
    561   ! Several investigators have proposed minimum mechanical mixing depth
    562   ! relationships as a function of the local friction velocity, u*.  We
    563   ! make use of a linear relationship of the form h = c u* where c=700.
    564   ! The scaling arguments that give rise to this relationship most often
    565   ! represent the coefficient c as some constant over the local coriolis
    566   ! parameter.  Here we make use of the experimental results of Koracin
    567   ! and Berkowicz (1988) [BLM, Vol 43] for wich they recommend 0.07/f
    568   ! where f was evaluated at 39.5 N and 52 N.  Thus we use a typical mid
    569   ! latitude value for f so that c = 0.07/f = 700.
    570 
    571   DO i = 1, knon
    572     pblmin = 700.0*ustar(i)
    573     pblh(i) = max(pblh(i), pblmin)
    574     ! par exemple :
    575     pblt(i) = t(i, 2) + (t(i,3)-t(i,2))*(pblh(i)-z(i,2))/(z(i,3)-z(i,2))
    576   END DO
    577 
    578   ! ********************************************************************
    579   ! pblh is now available; do preparation for diffusivity calculation :
    580   ! ********************************************************************
    581   DO i = 1, knon
    582     check(i) = .TRUE.
    583     zsat(i) = .FALSE.
    584     ! omegafl utilise pour prolongement CAPE
    585     omegafl(i) = .FALSE.
    586     cape(i) = 0.
    587     kape(i) = 0.
    588     eauliq(i) = 0.
    589     ctei(i) = 0.
    590     pblk(i) = 0.0
    591     fak1(i) = ustar(i)*pblh(i)*vk
    592 
    593     ! Do additional preparation for unstable cases only, set temperature
    594     ! and moisture perturbations depending on stability.
    595     ! *** Rq: les formule sont prises dans leur forme CS ***
    596     IF (unstbl(i)) THEN
    597       ! AM Niveau de ref du thermique
    598       ! AM          zxt=(t(i,1)-z(i,1)*0.5*RG/RCPD/(1.+RVTMP2*q(i,1)))
    599       ! AM     .         *(1.+RETV*q(i,1))
    600       zxt = (th_th(i)-zref*0.5*rg/rcpd/(1.+rvtmp2*qt_th(i)))* &
    601         (1.+retv*qt_th(i))
    602       phiminv(i) = (1.-binm*pblh(i)*unsobklen(i))**onet
    603       phihinv(i) = sqrt(1.-binh*pblh(i)*unsobklen(i))
    604       wm(i) = ustar(i)*phiminv(i)
    605       fak2(i) = wm(i)*pblh(i)*vk
    606       wstar(i) = (heatv(i)*rg*pblh(i)/zxt)**onet
    607       fak3(i) = fakn*wstar(i)/wm(i)
    608     ELSE
    609       wstar(i) = 0.
    610     END IF
    611     ! Computes Theta_e for thermal (all cases : to be modified)
    612     ! attention ajout therm(i) = virtuelle
    613     the_th(i) = th_th(i) + therm(i) + rlvcp*qt_th(i)
    614     ! ou:    The_th(i) = Th_th(i) + sqrt(b1)*khfs(i)/wm(i) + RLvCp*qT_th(i)
    615   END DO
    616 
    617   ! Main level loop to compute the diffusivities and
    618   ! counter-gradient terms:
    619 
    620   DO k = 2, isommet
    621 
    622     ! Find levels within boundary layer:
    623 
    624     DO i = 1, knon
    625       unslev(i) = .FALSE.
    626       stblev(i) = .FALSE.
    627       zm(i) = z(i, k-1)
    628       zp(i) = z(i, k)
    629       IF (zkmin==0.0 .AND. zp(i)>pblh(i)) zp(i) = pblh(i)
    630       IF (zm(i)<pblh(i)) THEN
    631         zmzp = 0.5*(zm(i)+zp(i))
    632         ! debug
    633         ! if (i.EQ.1864) then
    634         ! print*,'i,pblh(1864),obklen(1864)',i,pblh(i),obklen(i)
    635         ! endif
    636 
    637         zh(i) = zmzp/pblh(i)
    638         zl(i) = zmzp*unsobklen(i)
    639         zzh(i) = 0.
    640         IF (zh(i)<=1.0) zzh(i) = (1.-zh(i))**2
    641 
    642         ! stblev for points zm < plbh and stable and neutral
    643         ! unslev for points zm < plbh and unstable
    644 
    645         IF (unstbl(i)) THEN
    646           unslev(i) = .TRUE.
    647         ELSE
    648           stblev(i) = .TRUE.
    649         END IF
    650       END IF
    651     END DO
    652     ! print*,'fin calcul niveaux'
    653 
    654     ! Stable and neutral points; set diffusivities; counter-gradient
    655     ! terms zero for stable case:
    656 
    657     DO i = 1, knon
    658       IF (stblev(i)) THEN
    659         IF (zl(i)<=1.) THEN
    660           pblk(i) = fak1(i)*zh(i)*zzh(i)/(1.+betas*zl(i))
    661         ELSE
    662           pblk(i) = fak1(i)*zh(i)*zzh(i)/(betas+zl(i))
    663         END IF
    664         ! pcfm(i,k) = pblk(i)
    665         ! pcfh(i,k) = pcfm(i,k)
    666       END IF
    667     END DO
    668 
    669     ! unssrf, unstable within surface layer of pbl
    670     ! unsout, unstable within outer   layer of pbl
    671 
    672     DO i = 1, knon
    673       unssrf(i) = .FALSE.
    674       unsout(i) = .FALSE.
    675       IF (unslev(i)) THEN
    676         IF (zh(i)<sffrac) THEN
    677           unssrf(i) = .TRUE.
    678         ELSE
    679           unsout(i) = .TRUE.
    680         END IF
    681       END IF
    682     END DO
    683 
    684     ! Unstable for surface layer; counter-gradient terms zero
    685 
    686     DO i = 1, knon
    687       IF (unssrf(i)) THEN
    688         term = (1.-betam*zl(i))**onet
    689         pblk(i) = fak1(i)*zh(i)*zzh(i)*term
    690         pr(i) = term/sqrt(1.-betah*zl(i))
    691       END IF
    692     END DO
    693     ! print*,'fin counter-gradient terms zero'
    694 
    695     ! Unstable for outer layer; counter-gradient terms non-zero:
    696 
    697     DO i = 1, knon
    698       IF (unsout(i)) THEN
    699         pblk(i) = fak2(i)*zh(i)*zzh(i)
    700         ! cgs(i,k) = fak3(i)/(pblh(i)*wm(i))
    701         ! cgh(i,k) = khfs(i)*cgs(i,k)
    702         pr(i) = phiminv(i)/phihinv(i) + ccon*fak3(i)/fak
    703         ! cgq(i,k) = kqfs(i)*cgs(i,k)
    704       END IF
    705     END DO
    706     ! print*,'fin counter-gradient terms non zero'
    707 
    708     ! For all unstable layers, compute diffusivities and ctrgrad ter m
    709 
    710     ! DO i = 1, knon
    711     ! IF (unslev(i)) THEN
    712     ! pcfm(i,k) = pblk(i)
    713     ! pcfh(i,k) = pblk(i)/pr(i)
    714     ! etc cf original
    715     ! ENDIF
    716     ! ENDDO
    717 
    718     ! For all layers, compute integral info and CTEI
    719 
    720     DO i = 1, knon
    721       IF (check(i) .OR. omegafl(i)) THEN
    722         IF (.NOT. zsat(i)) THEN
    723           ! Th2 = The_th(i) - RLvCp*qT_th(i)
    724           th2 = th_th(i)
    725           t2 = th2*s(i, k)
    726           ! thermodyn functions
    727           zdelta = max(0., sign(1.,rtt-t2))
    728           qqsat = r2es*foeew(t2, zdelta)/pplay(i, k)
    729           qqsat = min(0.5, qqsat)
    730           zcor = 1./(1.-retv*qqsat)
    731           qqsat = qqsat*zcor
    732 
    733           IF (qqsat<qt_th(i)) THEN
    734             ! on calcule lcl
    735             IF (k==2) THEN
    736               plcl(i) = z(i, k)
    737             ELSE
    738               plcl(i) = z(i, k-1) + (z(i,k-1)-z(i,k))*(qt_th(i)-qsatbef(i))/( &
    739                 qsatbef(i)-qqsat)
    740             END IF
    741             zsat(i) = .TRUE.
    742             tbef(i) = t2
    743           END IF
    744 
    745           qsatbef(i) = qqsat ! bug dans la version orig ???
    746         END IF
    747         ! amn ???? cette ligne a deja ete faite normalement ?
    748       END IF
    749       ! print*,'hbtm2 i,k=',i,k
    750     END DO
    751   END DO ! end of level loop
    752   ! IM 170305 BEG
    753   IF (1==0) THEN
    754     PRINT *, 'hbtm2  ok'
    755   END IF !(1.EQ.0) THEN
    756   ! IM 170305 END
    757   RETURN
    758 END SUBROUTINE hbtm
     456             IF (t_star<0. .OR. q_star<0.) THEN
     457                PRINT *, 'i t_star q_star khfs kqfs wm', i, t_star, q_star, &
     458                     khfs(i), kqfs(i), wm(i)
     459             END IF
     460          END IF
     461          ! IM 091204 END
     462          ! AM Nveau cde ref 2m =>
     463          ! AM        therm(i) = sqrt( b1*(1.+2.*RETV*q(i,1))*t_star**2
     464          ! AM     +             + (RETV*T(i,1))**2*b2*q_star**2
     465          ! AM     +             + 2.*RETV*T(i,1)*b212*q_star*t_star
     466          ! AM     +                 )
     467          ! IM 091204 BEG
     468          a1 = b1*(1.+2.*retv*qt_th(i))*t_star**2
     469          a2 = (retv*th_th(i))**2*b2*q_star*q_star
     470          a3 = 2.*retv*th_th(i)*b212*q_star*t_star
     471          aa = a1 + a2 + a3
     472          IF (1==0) THEN
     473             IF (aa<0.) THEN
     474                PRINT *, 'i a1 a2 a3 aa', i, a1, a2, a3, aa
     475                PRINT *, 'i qT_th Th_th t_star q_star RETV b1 b2 b212', i, &
     476                     qt_th(i), th_th(i), t_star, q_star, retv, b1, b2, b212
     477             END IF
     478          END IF
     479          ! IM 091204 END
     480          therm(i) = sqrt(b1*(1.+2.*retv*qt_th(i))*t_star**2+(retv*th_th( &
     481               i))**2*b2*q_star*q_star &  ! IM 101204  +             +
     482                                ! 2.*RETV*Th_th(i)*b212*q_star*t_star
     483               +max(0.,2.*retv*th_th(i)*b212*q_star*t_star))
     484
     485          ! Theta et qT du thermique (forme H&B) avec exces
     486          ! (attention, on ajoute therm(i) qui est virtuelle ...)
     487          ! pourquoi pas sqrt(b1)*t_star ?
     488          ! dqs = b2sr*kqfs(i)/wm(i)
     489          qt_th(i) = qt_th(i) + b2sr*q_star
     490          ! new on differre le calcul de Theta_e
     491          ! The_th(i) = The_th(i) + therm(i) + RLvCp*qT_th(i)
     492          ! ou:    The_th(i) = The_th(i) + sqrt(b1)*khfs(i)/wm(i) +
     493          ! RLvCp*qT_th(i)
     494          rhino(i, 1) = 0.0
     495       END IF
     496    END DO
     497
     498    ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     499    ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     500    ! ++ Improve pblh estimate for unstable conditions using the +++++++
     501    ! ++          convective temperature excess :                +++++++
     502    ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     503    ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     504
     505    DO k = 2, isommet
     506       DO i = 1, knon
     507          IF (check(i)) THEN
     508             ! test     zdu2 =
     509             ! (u(i,k)-u(i,1))**2+(v(i,k)-v(i,1))**2+fac*ustar(i)**2
     510             zdu2 = u(i, k)**2 + v(i, k)**2
     511             zdu2 = max(zdu2, 1.0E-20)
     512             ! Theta_v environnement
     513             zthvd = t(i, k)/s(i, k)*(1.+retv*q(i,k))
     514
     515             ! et therm Theta_v (avec hypothese de constance de H&B,
     516             ! zthvu=(t(i,1)+therm(i))/s(i,1)*(1.+RETV*q(i,1))
     517             zthvu = th_th(i)*(1.+retv*qt_th(i)) + therm(i)
     518
     519
     520             ! Le Ri par Theta_v
     521             ! AM Niveau de ref 2m
     522             ! AM         rhino(i,k) = (z(i,k)-z(i,1))*RG*(zthvd-zthvu)
     523             ! AM     .               /(zdu2*0.5*(zthvd+zthvu))
     524             rhino(i, k) = (z(i,k)-zref)*rg*(zthvd-zthvu)/(zdu2*0.5&
     525                  *(zthvd+zthvu))
     526
     527
     528             IF (rhino(i,k)>=ricr) THEN
     529                pblh(i) = z(i, k-1) + (z(i,k-1)-z(i,k))*(ricr-rhino(i,k-1))&
     530                     /(rhino(i,k-1)-rhino(i,k))
     531                ! test04
     532                pblh(i) = pblh(i) + 100.
     533                pblt(i) = t(i, k-1) + (t(i,k)-t(i,k-1))*(pblh(i)-z(i,k-1))&
     534                     /(z(i,k)- z(i,k-1))
     535                check(i) = .FALSE.
     536                ! IM 170305 BEG
     537                IF (1==0) THEN
     538                   ! debug print -120;34       -34-        58 et    0;26 wamp
     539                   IF (i==950 .OR. i==192 .OR. i==624 .OR. i==118) THEN
     540                      PRINT *, ' i,Th_th,Therm,qT :', i, th_th(i), therm(i), &
     541                           qt_th(i)
     542                      q_star = kqfs(i)/wm(i)
     543                      t_star = khfs(i)/wm(i)
     544                      PRINT *, 'q* t*, b1,b2,b212 ', q_star, t_star, &
     545                           b1*(1.+2.*retv*qt_th(i))*t_star**2, &
     546                           (retv*th_th(i))**2*b2*q_star**2, 2.*retv*th_th(i)&
     547                           *b212*q_star *t_star
     548                      PRINT *, 'zdu2 ,100.*ustar(i)**2', zdu2, fac*ustar(i)**2
     549                   END IF
     550                END IF !(1.EQ.0) THEN
     551                ! IM 170305 END
     552                ! q_star = kqfs(i)/wm(i)
     553                ! t_star = khfs(i)/wm(i)
     554                ! trmb1(i) = b1*(1.+2.*RETV*q(i,1))*t_star**2
     555                ! trmb2(i) = (RETV*T(i,1))**2*b2*q_star**2
     556                ! Omega now   trmb3(i) = 2.*RETV*T(i,1)*b212*q_star*t_star
     557             END IF
     558          END IF
     559       END DO
     560    END DO
     561
     562    ! Set pbl height to maximum value where computation exceeds number of
     563    ! layers allowed
     564
     565    DO i = 1, knon
     566       IF (check(i)) pblh(i) = z(i, isommet)
     567    END DO
     568
     569    ! PBL height must be greater than some minimum mechanical mixing depth
     570    ! Several investigators have proposed minimum mechanical mixing depth
     571    ! relationships as a function of the local friction velocity, u*.  We
     572    ! make use of a linear relationship of the form h = c u* where c=700.
     573    ! The scaling arguments that give rise to this relationship most often
     574    ! represent the coefficient c as some constant over the local coriolis
     575    ! parameter.  Here we make use of the experimental results of Koracin
     576    ! and Berkowicz (1988) [BLM, Vol 43] for wich they recommend 0.07/f
     577    ! where f was evaluated at 39.5 N and 52 N.  Thus we use a typical mid
     578    ! latitude value for f so that c = 0.07/f = 700.
     579
     580    DO i = 1, knon
     581       pblmin = 700.0*ustar(i)
     582       pblh(i) = max(pblh(i), pblmin)
     583       ! par exemple :
     584       pblt(i) = t(i, 2) + (t(i,3)-t(i,2))*(pblh(i)-z(i,2))/(z(i,3)-z(i,2))
     585    END DO
     586
     587    ! ********************************************************************
     588    ! pblh is now available; do preparation for diffusivity calculation :
     589    ! ********************************************************************
     590    DO i = 1, knon
     591       check(i) = .TRUE.
     592       zsat(i) = .FALSE.
     593       ! omegafl utilise pour prolongement CAPE
     594       omegafl(i) = .FALSE.
     595       cape(i) = 0.
     596       kape(i) = 0.
     597       eauliq(i) = 0.
     598       ctei(i) = 0.
     599       pblk(i) = 0.0
     600       fak1(i) = ustar(i)*pblh(i)*vk
     601
     602       ! Do additional preparation for unstable cases only, set temperature
     603       ! and moisture perturbations depending on stability.
     604       ! *** Rq: les formule sont prises dans leur forme CS ***
     605       IF (unstbl(i)) THEN
     606          ! AM Niveau de ref du thermique
     607          ! AM          zxt=(t(i,1)-z(i,1)*0.5*RG/RCPD/(1.+RVTMP2*q(i,1)))
     608          ! AM     .         *(1.+RETV*q(i,1))
     609          zxt = (th_th(i)-zref*0.5*rg/rcpd/(1.+rvtmp2*qt_th(i)))* &
     610               (1.+retv*qt_th(i))
     611          phiminv(i) = (1.-binm*pblh(i)*unsobklen(i))**onet
     612          phihinv(i) = sqrt(1.-binh*pblh(i)*unsobklen(i))
     613          wm(i) = ustar(i)*phiminv(i)
     614          fak2(i) = wm(i)*pblh(i)*vk
     615          wstar(i) = (heatv(i)*rg*pblh(i)/zxt)**onet
     616          fak3(i) = fakn*wstar(i)/wm(i)
     617       ELSE
     618          wstar(i) = 0.
     619       END IF
     620       ! Computes Theta_e for thermal (all cases : to be modified)
     621       ! attention ajout therm(i) = virtuelle
     622       the_th(i) = th_th(i) + therm(i) + rlvcp*qt_th(i)
     623       ! ou:    The_th(i) = Th_th(i) + sqrt(b1)*khfs(i)/wm(i) + RLvCp*qT_th(i)
     624    END DO
     625
     626    ! Main level loop to compute the diffusivities and
     627    ! counter-gradient terms:
     628
     629    DO k = 2, isommet
     630
     631       ! Find levels within boundary layer:
     632
     633       DO i = 1, knon
     634          unslev(i) = .FALSE.
     635          stblev(i) = .FALSE.
     636          zm(i) = z(i, k-1)
     637          zp(i) = z(i, k)
     638          IF (zkmin==0.0 .AND. zp(i)>pblh(i)) zp(i) = pblh(i)
     639          IF (zm(i)<pblh(i)) THEN
     640             zmzp = 0.5*(zm(i)+zp(i))
     641             ! debug
     642             ! if (i.EQ.1864) then
     643             ! print*,'i,pblh(1864),obklen(1864)',i,pblh(i),obklen(i)
     644             ! endif
     645
     646             zh(i) = zmzp/pblh(i)
     647             zl(i) = zmzp*unsobklen(i)
     648             zzh(i) = 0.
     649             IF (zh(i)<=1.0) zzh(i) = (1.-zh(i))**2
     650
     651             ! stblev for points zm < plbh and stable and neutral
     652             ! unslev for points zm < plbh and unstable
     653
     654             IF (unstbl(i)) THEN
     655                unslev(i) = .TRUE.
     656             ELSE
     657                stblev(i) = .TRUE.
     658             END IF
     659          END IF
     660       END DO
     661       ! print*,'fin calcul niveaux'
     662
     663       ! Stable and neutral points; set diffusivities; counter-gradient
     664       ! terms zero for stable case:
     665
     666       DO i = 1, knon
     667          IF (stblev(i)) THEN
     668             IF (zl(i)<=1.) THEN
     669                pblk(i) = fak1(i)*zh(i)*zzh(i)/(1.+betas*zl(i))
     670             ELSE
     671                pblk(i) = fak1(i)*zh(i)*zzh(i)/(betas+zl(i))
     672             END IF
     673             ! pcfm(i,k) = pblk(i)
     674             ! pcfh(i,k) = pcfm(i,k)
     675          END IF
     676       END DO
     677
     678       ! unssrf, unstable within surface layer of pbl
     679       ! unsout, unstable within outer   layer of pbl
     680
     681       DO i = 1, knon
     682          unssrf(i) = .FALSE.
     683          unsout(i) = .FALSE.
     684          IF (unslev(i)) THEN
     685             IF (zh(i)<sffrac) THEN
     686                unssrf(i) = .TRUE.
     687             ELSE
     688                unsout(i) = .TRUE.
     689             END IF
     690          END IF
     691       END DO
     692
     693       ! Unstable for surface layer; counter-gradient terms zero
     694
     695       DO i = 1, knon
     696          IF (unssrf(i)) THEN
     697             term = (1.-betam*zl(i))**onet
     698             pblk(i) = fak1(i)*zh(i)*zzh(i)*term
     699             pr(i) = term/sqrt(1.-betah*zl(i))
     700          END IF
     701       END DO
     702       ! print*,'fin counter-gradient terms zero'
     703
     704       ! Unstable for outer layer; counter-gradient terms non-zero:
     705
     706       DO i = 1, knon
     707          IF (unsout(i)) THEN
     708             pblk(i) = fak2(i)*zh(i)*zzh(i)
     709             ! cgs(i,k) = fak3(i)/(pblh(i)*wm(i))
     710             ! cgh(i,k) = khfs(i)*cgs(i,k)
     711             pr(i) = phiminv(i)/phihinv(i) + ccon*fak3(i)/fak
     712             ! cgq(i,k) = kqfs(i)*cgs(i,k)
     713          END IF
     714       END DO
     715       ! print*,'fin counter-gradient terms non zero'
     716
     717       ! For all unstable layers, compute diffusivities and ctrgrad ter m
     718
     719       ! DO i = 1, knon
     720       ! IF (unslev(i)) THEN
     721       ! pcfm(i,k) = pblk(i)
     722       ! pcfh(i,k) = pblk(i)/pr(i)
     723       ! etc cf original
     724       ! ENDIF
     725       ! ENDDO
     726
     727       ! For all layers, compute integral info and CTEI
     728
     729       DO i = 1, knon
     730          IF (check(i) .OR. omegafl(i)) THEN
     731             IF (.NOT. zsat(i)) THEN
     732                ! Th2 = The_th(i) - RLvCp*qT_th(i)
     733                th2 = th_th(i)
     734                t2 = th2*s(i, k)
     735                ! thermodyn functions
     736                zdelta = max(0., sign(1.,rtt-t2))
     737                qqsat = r2es*foeew(t2, zdelta)/pplay(i, k)
     738                qqsat = min(0.5, qqsat)
     739                zcor = 1./(1.-retv*qqsat)
     740                qqsat = qqsat*zcor
     741
     742                IF (qqsat<qt_th(i)) THEN
     743                   ! on calcule lcl
     744                   IF (k==2) THEN
     745                      plcl(i) = z(i, k)
     746                   ELSE
     747                      plcl(i) = z(i, k-1) + (z(i,k-1)-z(i,k))*(qt_th(i)&
     748                           -qsatbef(i))/(qsatbef(i)-qqsat)
     749                   END IF
     750                   zsat(i) = .TRUE.
     751                   tbef(i) = t2
     752                END IF
     753
     754                qsatbef(i) = qqsat ! bug dans la version orig ???
     755             END IF
     756             ! amn ???? cette ligne a deja ete faite normalement ?
     757          END IF
     758          ! print*,'hbtm2 i,k=',i,k
     759       END DO
     760    END DO ! end of level loop
     761    ! IM 170305 BEG
     762    IF (1==0) THEN
     763       PRINT *, 'hbtm2  ok'
     764    END IF !(1.EQ.0) THEN
     765    ! IM 170305 END
     766
     767  END SUBROUTINE hbtm
     768
     769end module hbtm_mod

  • LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90

    r3772 r3774  
    286286    USE carbon_cycle_mod,   ONLY : carbon_cycle_cpl, carbon_cycle_tr, level_coupling_esm
    287287    USE carbon_cycle_mod,   ONLY : co2_send, nbcf_out, fields_out, yfields_out, cfname_out
     288    use hbtm_mod, only: hbtm
    288289    USE indice_sol_mod
    289290    USE time_phylmdz_mod,   ONLY : day_ini,annee_ref,itau_phy
Note: See TracChangeset for help on using the changeset viewer.