/[lmdze]/trunk/phylmd/Interface_surf/hbtm.f90
ViewVC logotype

Diff of /trunk/phylmd/Interface_surf/hbtm.f90

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/Sources/phylmd/hbtm.f revision 134 by guez, Wed Apr 29 15:47:56 2015 UTC trunk/phylmd/Interface_surf/hbtm.f90 revision 346 by guez, Mon Dec 9 20:15:29 2019 UTC
# Line 4  module HBTM_m Line 4  module HBTM_m
4    
5  contains  contains
6    
7    SUBROUTINE HBTM(knon, paprs, pplay, t2m, t10m, q2m, q10m, ustar, flux_t, &    SUBROUTINE HBTM(paprs, pplay, t2m, q2m, ustar, flux_t, flux_q, u, v, t, q, &
8         flux_q, u, v, t, q, pblh, cape, EauLiq, ctei, pblT, therm, trmb1, &         pblh, cape, EauLiq, ctei, pblT, therm, plcl)
        trmb2, trmb3, plcl)  
9    
10      ! D'après Holstag et Boville et Troen et Mahrt      ! D'apr\'es Holstag et Boville et Troen et Mahrt
11      ! JAS 47 BLM      ! JAS 47 BLM
12      ! Algorithme thèse Anne Mathieu  
13      ! Critère d'entraînement Peter Duynkerke (JAS 50)      ! Algorithme th\'ese Anne Mathieu. Crit\'ere d'entra\^inement
14      ! written by: Anne MATHIEU and Alain LAHELLEC, 22nd November 1999      ! Peter Duynkerke (JAS 50).  Written by: Anne MATHIEU and Alain
15      ! features : implem. exces Mathieu      ! LAHELLEC, 22nd November 1999.
16    
17      ! modifications : decembre 99 passage th a niveau plus bas. voir fixer      ! Modifications : d\'ecembre 99 passage th \`a niveau plus bas. Voir fixer
18      ! la prise du th a z/Lambda = -.2 (max Ray)      ! la prise du th \`a z/Lambda = -.2 (max Ray)
19      ! Autre algo : entrainement ~ Theta+v =cste mais comment=>The?      ! Autre algorithme : entra\^inement ~ Theta + v =constante
20      ! on peut fixer q a .7 qsat (cf. non adiabatique) => T2 et The2      !  mais comment ? The ?
21      ! voir aussi //KE pblh = niveau The_e ou l = env.      ! On peut fixer q \`a 0.7 qsat (cf. non adiabatique) d'où T2 et The2.
22        ! Voir aussi KE pblh = niveau The_e ou l = env.
23      ! fin therm a la HBTM passage a forme Mathieu 12/09/2001  
24        ! Adaptation \`a LMDZ version coupl\'ee. Pour le moment on fait
25      ! Adaptation a LMDZ version couplee      ! passer en argument les grandeurs de surface : flux, t, q2m. On
26      ! Pour le moment on fait passer en argument les grandeurs de surface :      ! va utiliser syst\'ematiquement les grandeurs \`a 2 m mais on
27      ! flux, t, q2m, t, q10m, on va utiliser systematiquement les grandeurs a 2m      ! garde la possibilit\'e de changer si besoin (jusqu'\`a pr\'esent
28      ! mais on garde la possibilité de changer si besoin est (jusqu'à présent      ! la forme de HB avec le premier niveau mod\`ele \'etait
29      ! la forme de HB avec le 1er niveau modele etait conservee)      ! conserv\'ee).
30    
31      USE dimphy, ONLY: klev, klon      USE dimphy, ONLY: klev, klon
32      USE suphec_m, ONLY: rcpd, rd, retv, rg, rkappa, rlvtt, rtt, rv      USE suphec_m, ONLY: rcpd, rd, retv, rg, rkappa, rtt
33      USE yoethf_m, ONLY: r2es, rvtmp2      USE yoethf_m, ONLY: r2es, rvtmp2
34      USE fcttre, ONLY: foeew      USE fcttre, ONLY: foeew
35    
     REAL RLvCp, REPS  
36      ! Arguments:      ! Arguments:
37    
     ! nombre de points a calculer  
     INTEGER, intent(in):: knon  
   
     REAL, intent(in):: t2m(klon) ! temperature a 2 m  
     real t10m(klon) ! temperature a 10 m  
     ! q a 2 et 10m  
     REAL q2m(klon), q10m(klon)  
     REAL ustar(klon)  
38      ! pression a inter-couche (Pa)      ! pression a inter-couche (Pa)
39      REAL paprs(klon, klev+1)      REAL, intent(in):: paprs(klon, klev+1)
40      ! pression au milieu de couche (Pa)      ! pression au milieu de couche (Pa)
41      REAL pplay(klon, klev)      REAL, intent(in):: pplay(klon, klev)
42      ! Flux      REAL, intent(in):: t2m(klon) ! temperature a 2 m
43      REAL flux_t(klon, klev), flux_q(klon, klev)      ! q a 2 et 10m
44      ! vitesse U (m/s)      REAL, intent(in):: q2m(klon)
45      REAL u(klon, klev)      REAL, intent(in):: ustar(:) ! (knon)
46      ! vitesse V (m/s)      REAL, intent(in):: flux_t(:), flux_q(:) ! (knon) flux à la surface
47      REAL v(klon, klev)      REAL, intent(in):: u(:, :) ! (knon, klev) vitesse U (m/s)
48      ! temperature (K)      REAL, intent(in):: v(:, :) ! (knon, klev) vitesse V (m/s)
49      REAL t(klon, klev)      REAL, intent(in):: t(:, :) ! (knon, klev) temperature (K)
50      ! vapeur d'eau (kg/kg)      REAL, intent(in):: q(:, :) ! (knon, klev) vapeur d'eau (kg/kg)
     REAL q(klon, klev)  
51    
52        REAL, intent(out):: pblh(:) ! (knon)
53        ! Cape du thermique
54        REAL Cape(klon)
55        ! Eau liqu integr du thermique
56        REAL EauLiq(klon)
57        ! Critere d'instab d'entrainmt des nuages de
58        REAL ctei(klon)
59        REAL pblT(klon)
60        ! thermal virtual temperature excess
61        REAL therm(klon)
62        REAL plcl(klon)
63    
64        ! Local:
65        
66        INTEGER knon ! nombre de points a calculer
67      INTEGER isommet      INTEGER isommet
68      ! limite max sommet pbl      ! limite max sommet pbl
69      PARAMETER (isommet=klev)      PARAMETER (isommet=klev)
# Line 68  contains Line 72  contains
72      PARAMETER (vk=0.35)      PARAMETER (vk=0.35)
73      REAL ricr      REAL ricr
74      PARAMETER (ricr=0.4)      PARAMETER (ricr=0.4)
     REAL fak  
     ! b calcul du Prandtl et de dTetas  
     PARAMETER (fak=8.5)  
     REAL fakn  
75      ! a      ! a
     PARAMETER (fakn=7.2)  
76      REAL onet      REAL onet
77      PARAMETER (onet=1.0/3.0)      PARAMETER (onet=1.0/3.0)
     REAL t_coup  
     PARAMETER(t_coup=273.15)  
78      REAL zkmin      REAL zkmin
79      PARAMETER (zkmin=0.01)      PARAMETER (zkmin=0.01)
80      REAL betam      REAL betam
81      ! pour Phim / h dans la S.L stable      ! pour Phim / h dans la S.L stable
82      PARAMETER (betam=15.0)      PARAMETER (betam=15.0)
     REAL betah  
     PARAMETER (betah=15.0)  
     REAL betas  
     ! Phit dans la S.L. stable (mais 2 formes /  
     PARAMETER (betas=5.0)  
83      ! z/OBL<>1      ! z/OBL<>1
84      REAL sffrac      REAL sffrac
85      ! S.L. = z/h < .1      ! S.L. = z/h < .1
86      PARAMETER (sffrac=0.1)      PARAMETER (sffrac=0.1)
87      REAL binm      REAL binm
88      PARAMETER (binm=betam*sffrac)      PARAMETER (binm=betam*sffrac)
     REAL binh  
     PARAMETER (binh=betah*sffrac)  
     REAL ccon  
     PARAMETER (ccon=fak*sffrac*vk)  
89    
90      REAL q_star, t_star      REAL q_star, t_star
91      ! Lambert correlations T' q' avec T* q*      ! Lambert correlations T' q' avec T* q*
# Line 123  contains Line 111  contains
111      REAL rhino(klon, klev)      REAL rhino(klon, klev)
112      ! pts w/unstbl pbl (positive virtual ht flx)      ! pts w/unstbl pbl (positive virtual ht flx)
113      LOGICAL unstbl(klon)      LOGICAL unstbl(klon)
     ! stable pbl with levels within pbl  
     LOGICAL stblev(klon)  
     ! unstbl pbl with levels within pbl  
     LOGICAL unslev(klon)  
     ! unstb pbl w/lvls within srf pbl lyr  
     LOGICAL unssrf(klon)  
     ! unstb pbl w/lvls in outer pbl lyr  
     LOGICAL unsout(klon)  
114      LOGICAL check(klon) ! Richardson number > critical      LOGICAL check(klon) ! Richardson number > critical
115      ! flag de prolongerment cape pour pt Omega      ! flag de prolongerment cape pour pt Omega
116      LOGICAL omegafl(klon)      LOGICAL omegafl(klon)
     REAL pblh(klon)  
     REAL pblT(klon)  
     REAL plcl(klon)  
117    
118      ! Monin-Obukhov lengh      ! Monin-Obukhov lengh
119      REAL obklen(klon)      REAL obklen(klon)
120    
121      REAL zdu2      REAL zdu2
     ! thermal virtual temperature excess  
     REAL therm(klon)  
     REAL trmb1(klon), trmb2(klon), trmb3(klon)  
122      ! Algorithme thermique      ! Algorithme thermique
123      REAL s(klon, klev) ! [P/Po]^Kappa milieux couches      REAL s(klon, klev) ! [P/Po]^Kappa milieux couches
     ! equivalent potential temperature of therma  
     REAL The_th(klon)  
124      ! total water of thermal      ! total water of thermal
125      REAL qT_th(klon)      REAL qT_th(klon)
126      ! T thermique niveau precedent      ! T thermique niveau precedent
     REAL Tbef(klon)  
127      REAL qsatbef(klon)      REAL qsatbef(klon)
128      ! le thermique est sature      ! le thermique est sature
129      LOGICAL Zsat(klon)      LOGICAL Zsat(klon)
130      ! Cape du thermique      REAL zthvd, zthvu, qqsat
131      REAL Cape(klon)      REAL t2
     ! Cape locale  
     REAL Kape(klon)  
     ! Eau liqu integr du thermique  
     REAL EauLiq(klon)  
     ! Critere d'instab d'entrainmt des nuages de  
     REAL ctei(klon)  
     REAL the1, the2, aa, zthvd, zthvu, xintpos, qqsat  
     REAL a1, a2, a3  
     REAL xhis, rnum, th1, thv1, thv2, ql2  
     REAL qsat2, qT1, q2, t1, t2, xnull  
     REAL quadsat, spblh, reduc  
132    
133      ! inverse phi function for momentum      ! inverse phi function for momentum
134      REAL phiminv(klon)      REAL phiminv(klon)
     ! inverse phi function for heat  
     REAL phihinv(klon)  
135      ! turbulent velocity scale for momentum      ! turbulent velocity scale for momentum
136      REAL wm(klon)      REAL wm(klon)
     ! k*ustar*pblh  
     REAL fak1(klon)  
     ! k*wm*pblh  
     REAL fak2(klon)  
     ! fakn*wstr/wm  
     REAL fak3(klon)  
     ! level eddy diffusivity for momentum  
     REAL pblk(klon)  
     ! Prandtl number for eddy diffusivities  
     REAL pr(klon)  
     ! zmzp / Obukhov length  
     REAL zl(klon)  
     ! zmzp / pblh  
     REAL zh(klon)  
     ! (1-(zmzp/pblh))**2  
     REAL zzh(klon)  
     ! w*, convective velocity scale  
     REAL wstr(klon)  
     ! current level height  
     REAL zm(klon)  
137      ! current level height + one level up      ! current level height + one level up
138      REAL zp(klon)      REAL zp(klon)
139      REAL zcor, zcvm5      REAL zcor
140    
141      REAL fac, pblmin, zmzp, term      REAL pblmin
142    
143      !-----------------------------------------------------------------      !-----------------------------------------------------------------
144    
145        knon = size(pblh)
146    
147      ! initialisations      ! initialisations
148      q_star = 0      q_star = 0
149      t_star = 0      t_star = 0
# Line 211  contains Line 151  contains
151      b212=sqrt(b1*b2)      b212=sqrt(b1*b2)
152      b2sr=sqrt(b2)      b2sr=sqrt(b2)
153    
     ! Initialisation  
     RLvCp = RLVTT/RCPD  
     REPS = RD/RV  
   
154      ! Calculer les hauteurs de chaque couche      ! Calculer les hauteurs de chaque couche
155      ! (geopotentielle Int_dp/ro = Int_[Rd.T.dp/p] z = geop/g)      ! (geopotentielle Int_dp/ro = Int_[Rd.T.dp/p] z = geop/g)
156      ! pourquoi ne pas utiliser Phi/RG ?      ! pourquoi ne pas utiliser Phi/RG ?
# Line 245  contains Line 181  contains
181         zxt = t2m(i)         zxt = t2m(i)
182    
183         ! convention >0 vers le bas ds lmdz         ! convention >0 vers le bas ds lmdz
184         khfs(i) = - flux_t(i, 1)*zxt*Rd / (RCPD*paprs(i, 1))         khfs(i) = - flux_t(i)*zxt*Rd / (RCPD*paprs(i, 1))
185         kqfs(i) = - flux_q(i, 1)*zxt*Rd / (paprs(i, 1))         kqfs(i) = - flux_q(i)*zxt*Rd / paprs(i, 1)
186         ! verifier que khfs et kqfs sont bien de la forme w'l'         ! verifier que khfs et kqfs sont bien de la forme w'l'
187         heatv(i) = khfs(i) + 0.608*zxt*kqfs(i)         heatv(i) = khfs(i) + 0.608*zxt*kqfs(i)
188         ! a comparer aussi aux sorties de clqh : flux_T/RoCp et flux_q/RoLv         ! a comparer aussi aux sorties de clqh : flux_T/RoCp et flux_q/RoLv
# Line 265  contains Line 201  contains
201         plcl(i) = 6000.         plcl(i) = 6000.
202         ! Lambda = -u*^3 / (alpha.g.kvon.<w'Theta'v>         ! Lambda = -u*^3 / (alpha.g.kvon.<w'Theta'v>
203         obklen(i) = -t(i, 1)*ustar(i)**3/(RG*vk*heatv(i))         obklen(i) = -t(i, 1)*ustar(i)**3/(RG*vk*heatv(i))
        trmb1(i) = 0.  
        trmb2(i) = 0.  
        trmb3(i) = 0.  
204      ENDDO      ENDDO
205    
206      ! PBL height calculation: Search for level of pbl. Scan upward      ! PBL height calculation: Search for level of pbl. Scan upward
207      ! until the Richardson number between the first level and the      ! until the Richardson number between the first level and the
208      ! current level exceeds the "critical" value.  (bonne idee Nu de      ! current level exceeds the "critical" value.  (bonne idee Nu de
209      ! separer le Ric et l'exces de temp du thermique)      ! separer le Ric et l'exces de temp du thermique)
     fac = 100.  
210      DO k = 2, isommet      DO k = 2, isommet
211         DO i = 1, knon         DO i = 1, knon
212            IF (check(i)) THEN            IF (check(i)) THEN
# Line 339  contains Line 271  contains
271            q_star = kqfs(i)/wm(i)            q_star = kqfs(i)/wm(i)
272            t_star = khfs(i)/wm(i)            t_star = khfs(i)/wm(i)
273    
           a1=b1*(1.+2.*RETV*qT_th(i))*t_star**2  
           a2=(RETV*T2m(i))**2*b2*q_star**2  
           a3=2.*RETV*T2m(i)*b212*q_star*t_star  
           aa=a1+a2+a3  
   
274            therm(i) = sqrt( b1*(1.+2.*RETV*qT_th(i))*t_star**2 &            therm(i) = sqrt( b1*(1.+2.*RETV*qT_th(i))*t_star**2 &
275                 + (RETV*T2m(i))**2*b2*q_star**2 &                 + (RETV*T2m(i))**2*b2*q_star**2 &
276                 + max(0., 2.*RETV*T2m(i)*b212*q_star*t_star))                 + max(0., 2.*RETV*T2m(i)*b212*q_star*t_star))
# Line 352  contains Line 279  contains
279            ! (attention, on ajoute therm(i) qui est virtuelle ...)            ! (attention, on ajoute therm(i) qui est virtuelle ...)
280            ! pourquoi pas sqrt(b1)*t_star ?            ! pourquoi pas sqrt(b1)*t_star ?
281            qT_th(i) = qT_th(i) + b2sr*q_star            qT_th(i) = qT_th(i) + b2sr*q_star
282            ! new on differre le calcul de Theta_e            ! new on diff\`ere le calcul de Theta_e
283            rhino(i, 1) = 0.            rhino(i, 1) = 0.
284         ENDIF         ENDIF
285      ENDDO      ENDDO
# Line 419  contains Line 346  contains
346         ! omegafl utilise pour prolongement CAPE         ! omegafl utilise pour prolongement CAPE
347         omegafl(i) = .FALSE.         omegafl(i) = .FALSE.
348         Cape(i) = 0.         Cape(i) = 0.
        Kape(i) = 0.  
349         EauLiq(i) = 0.         EauLiq(i) = 0.
350         CTEI(i) = 0.         CTEI(i) = 0.
        pblk(i) = 0.0  
        fak1(i) = ustar(i)*pblh(i)*vk  
351    
352         ! Do additional preparation for unstable cases only, set temperature         ! Do additional preparation for unstable cases only, set temperature
353         ! and moisture perturbations depending on stability.         ! and moisture perturbations depending on stability.
# Line 433  contains Line 357  contains
357            zxt=(T2m(i)-zref*0.5*RG/RCPD/(1.+RVTMP2*qT_th(i))) &            zxt=(T2m(i)-zref*0.5*RG/RCPD/(1.+RVTMP2*qT_th(i))) &
358                 *(1.+RETV*qT_th(i))                 *(1.+RETV*qT_th(i))
359            phiminv(i) = (1. - binm*pblh(i)/obklen(i))**onet            phiminv(i) = (1. - binm*pblh(i)/obklen(i))**onet
           phihinv(i) = sqrt(1. - binh*pblh(i)/obklen(i))  
360            wm(i) = ustar(i)*phiminv(i)            wm(i) = ustar(i)*phiminv(i)
           fak2(i) = wm(i)*pblh(i)*vk  
           wstr(i) = (heatv(i)*RG*pblh(i)/zxt)**onet  
           fak3(i) = fakn*wstr(i)/wm(i)  
361         ENDIF         ENDIF
        ! Computes Theta_e for thermal (all cases : to be modified)  
        ! attention ajout therm(i) = virtuelle  
        The_th(i) = T2m(i) + therm(i) + RLvCp*qT_th(i)  
362      ENDDO      ENDDO
363    
364      ! Main level loop to compute the diffusivities and      ! Main level loop to compute the diffusivities and
# Line 449  contains Line 366  contains
366      loop_level: DO k = 2, isommet      loop_level: DO k = 2, isommet
367         ! Find levels within boundary layer:         ! Find levels within boundary layer:
368         DO i = 1, knon         DO i = 1, knon
           unslev(i) = .FALSE.  
           stblev(i) = .FALSE.  
           zm(i) = z(i, k-1)  
369            zp(i) = z(i, k)            zp(i) = z(i, k)
370            IF (zkmin == 0. .AND. zp(i) > pblh(i)) zp(i) = pblh(i)            IF (zkmin == 0. .AND. zp(i) > pblh(i)) zp(i) = pblh(i)
           IF (zm(i) < pblh(i)) THEN  
              zmzp = 0.5*(zm(i) + zp(i))  
              zh(i) = zmzp/pblh(i)  
              zl(i) = zmzp/obklen(i)  
              zzh(i) = 0.  
              IF (zh(i) <= 1.) zzh(i) = (1. - zh(i))**2  
   
              ! stblev for points zm < plbh and stable and neutral  
              ! unslev for points zm < plbh and unstable  
              IF (unstbl(i)) THEN  
                 unslev(i) = .TRUE.  
              ELSE  
                 stblev(i) = .TRUE.  
              ENDIF  
           ENDIF  
        ENDDO  
   
        ! Stable and neutral points; set diffusivities; counter-gradient  
        ! terms zero for stable case:  
        DO i = 1, knon  
           IF (stblev(i)) THEN  
              IF (zl(i) <= 1.) THEN  
                 pblk(i) = fak1(i)*zh(i)*zzh(i)/(1. + betas*zl(i))  
              ELSE  
                 pblk(i) = fak1(i)*zh(i)*zzh(i)/(betas + zl(i))  
              ENDIF  
           ENDIF  
        ENDDO  
   
        ! unssrf, unstable within surface layer of pbl  
        ! unsout, unstable within outer layer of pbl  
        DO i = 1, knon  
           unssrf(i) = .FALSE.  
           unsout(i) = .FALSE.  
           IF (unslev(i)) THEN  
              IF (zh(i) < sffrac) THEN  
                 unssrf(i) = .TRUE.  
              ELSE  
                 unsout(i) = .TRUE.  
              ENDIF  
           ENDIF  
        ENDDO  
   
        ! Unstable for surface layer; counter-gradient terms zero  
        DO i = 1, knon  
           IF (unssrf(i)) THEN  
              term = (1. - betam*zl(i))**onet  
              pblk(i) = fak1(i)*zh(i)*zzh(i)*term  
              pr(i) = term/sqrt(1. - betah*zl(i))  
           ENDIF  
        ENDDO  
   
        ! Unstable for outer layer; counter-gradient terms non-zero:  
        DO i = 1, knon  
           IF (unsout(i)) THEN  
              pblk(i) = fak2(i)*zh(i)*zzh(i)  
              pr(i) = phiminv(i)/phihinv(i) + ccon*fak3(i)/fak  
           ENDIF  
371         ENDDO         ENDDO
372    
373         ! For all layers, compute integral info and CTEI         ! For all layers, compute integral info and CTEI
# Line 534  contains Line 390  contains
390                             * (qT_th(i)-qsatbef(i)) / (qsatbef(i)-qqsat)                             * (qT_th(i)-qsatbef(i)) / (qsatbef(i)-qqsat)
391                     endif                     endif
392                     Zsat(i) = .true.                     Zsat(i) = .true.
                    Tbef(i) = T2  
393                  endif                  endif
394               endif               endif
395               qsatbef(i) = qqsat               qsatbef(i) = qqsat
Legend:
Removed from v.134  
changed lines
  Added in v.346
  ViewVC Help
Powered by ViewVC 1.1.21