/[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/libf/phylmd/hbtm.f90 revision 37 by guez, Tue Dec 21 15:45:48 2010 UTC trunk/phylmd/Interface_surf/hbtm.f revision 286 by guez, Tue Jul 24 15:22:48 2018 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)  
   
     use dimens_m  
     use dimphy  
     use YOMCST  
     use yoethf  
     use fcttre  
9    
10      ! D'apres Holstag & Boville et Troen & Mahrt      ! D'apr\'es Holstag et Boville et Troen et Mahrt
11      ! JAS 47 BLM      ! JAS 47 BLM
     ! Algorithme thèse Anne Mathieu  
     ! Critère d'entraînement Peter Duynkerke (JAS 50)  
     ! written by: Anne MATHIEU and Alain LAHELLEC, 22nd November 1999  
     ! features : implem. exces Mathieu  
   
     ! modifications : decembre 99 passage th a niveau plus bas. voir fixer  
     ! la prise du th a z/Lambda = -.2 (max Ray)  
     ! Autre algo : entrainement ~ Theta+v =cste mais comment=>The?  
     ! on peut fixer q a .7 qsat (cf. non adiabatique) => T2 et The2  
     ! voir aussi //KE pblh = niveau The_e ou l = env.  
   
     ! fin therm a la HBTM passage a forme Mathieu 12/09/2001  
   
     ! Adaptation a LMDZ version couplee  
     ! Pour le moment on fait passer en argument les grandeurs de surface :  
     ! flux, t, q2m, t, q10m, on va utiliser systematiquement les grandeurs a 2m  
     ! mais on garde la possibilité de changer si besoin est (jusqu'à présent  
     ! la forme de HB avec le 1er niveau modele etait conservee)  
12    
13      REAL RLvCp, REPS      ! Algorithme th\'ese Anne Mathieu. Crit\'ere d'entra\^inement
14      ! Arguments:      ! Peter Duynkerke (JAS 50).  Written by: Anne MATHIEU and Alain
15        ! LAHELLEC, 22nd November 1999.
16    
17        ! Modifications : d\'ecembre 99 passage th \`a niveau plus bas. Voir fixer
18        ! la prise du th \`a z/Lambda = -.2 (max Ray)
19        ! Autre algorithme : entra\^inement ~ Theta + v =constante
20        !  mais comment ? The ?
21        ! 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    
24        ! Adaptation \`a LMDZ version coupl\'ee. Pour le moment on fait
25        ! passer en argument les grandeurs de surface : flux, t, q2m. On
26        ! va utiliser syst\'ematiquement les grandeurs \`a 2 m mais on
27        ! garde la possibilit\'e de changer si besoin (jusqu'\`a pr\'esent
28        ! la forme de HB avec le premier niveau mod\`ele \'etait
29        ! conserv\'ee).
30    
31        USE dimphy, ONLY: klev, klon
32        USE suphec_m, ONLY: rcpd, rd, retv, rg, rkappa, rtt
33        USE yoethf_m, ONLY: r2es, rvtmp2
34        USE fcttre, ONLY: foeew
35    
36      ! nombre de points a calculer      ! Arguments:
     INTEGER, intent(in):: knon  
37    
     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)  
48        REAL, intent(in):: u(klon, klev) ! vitesse U (m/s)
49        REAL, intent(in):: v(klon, klev) ! vitesse V (m/s)
50    
51      ! temperature (K)      ! temperature (K)
52      REAL t(klon, klev)      REAL, intent(in):: t(klon, klev)
53      ! vapeur d'eau (kg/kg)      ! vapeur d'eau (kg/kg)
54      REAL q(klon, klev)      REAL, intent(in):: q(klon, klev)
55    
56        REAL, intent(out):: pblh(:) ! (knon)
57        ! Cape du thermique
58        REAL Cape(klon)
59        ! Eau liqu integr du thermique
60        REAL EauLiq(klon)
61        ! Critere d'instab d'entrainmt des nuages de
62        REAL ctei(klon)
63        REAL pblT(klon)
64        ! thermal virtual temperature excess
65        REAL therm(klon)
66        REAL plcl(klon)
67    
68        ! Local:
69        
70        INTEGER knon ! nombre de points a calculer
71      INTEGER isommet      INTEGER isommet
72      ! limite max sommet pbl      ! limite max sommet pbl
73      PARAMETER (isommet=klev)      PARAMETER (isommet=klev)
# Line 69  contains Line 76  contains
76      PARAMETER (vk=0.35)      PARAMETER (vk=0.35)
77      REAL ricr      REAL ricr
78      PARAMETER (ricr=0.4)      PARAMETER (ricr=0.4)
     REAL fak  
     ! b calcul du Prandtl et de dTetas  
     PARAMETER (fak=8.5)  
     REAL fakn  
79      ! a      ! a
     PARAMETER (fakn=7.2)  
80      REAL onet      REAL onet
81      PARAMETER (onet=1.0/3.0)      PARAMETER (onet=1.0/3.0)
     REAL t_coup  
     PARAMETER(t_coup=273.15)  
82      REAL zkmin      REAL zkmin
83      PARAMETER (zkmin=0.01)      PARAMETER (zkmin=0.01)
84      REAL betam      REAL betam
85      ! pour Phim / h dans la S.L stable      ! pour Phim / h dans la S.L stable
86      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)  
87      ! z/OBL<>1      ! z/OBL<>1
88      REAL sffrac      REAL sffrac
89      ! S.L. = z/h < .1      ! S.L. = z/h < .1
90      PARAMETER (sffrac=0.1)      PARAMETER (sffrac=0.1)
91      REAL binm      REAL binm
92      PARAMETER (binm=betam*sffrac)      PARAMETER (binm=betam*sffrac)
     REAL binh  
     PARAMETER (binh=betah*sffrac)  
     REAL ccon  
     PARAMETER (ccon=fak*sffrac*vk)  
93    
94      REAL q_star, t_star      REAL q_star, t_star
95      ! Lambert correlations T' q' avec T* q*      ! Lambert correlations T' q' avec T* q*
# Line 124  contains Line 115  contains
115      REAL rhino(klon, klev)      REAL rhino(klon, klev)
116      ! pts w/unstbl pbl (positive virtual ht flx)      ! pts w/unstbl pbl (positive virtual ht flx)
117      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)  
118      LOGICAL check(klon) ! Richardson number > critical      LOGICAL check(klon) ! Richardson number > critical
119      ! flag de prolongerment cape pour pt Omega      ! flag de prolongerment cape pour pt Omega
120      LOGICAL omegafl(klon)      LOGICAL omegafl(klon)
     REAL pblh(klon)  
     REAL pblT(klon)  
     REAL plcl(klon)  
121    
122      ! Monin-Obukhov lengh      ! Monin-Obukhov lengh
123      REAL obklen(klon)      REAL obklen(klon)
124    
125      REAL zdu2      REAL zdu2
     ! thermal virtual temperature excess  
     REAL therm(klon)  
     REAL trmb1(klon), trmb2(klon), trmb3(klon)  
126      ! Algorithme thermique      ! Algorithme thermique
127      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)  
128      ! total water of thermal      ! total water of thermal
129      REAL qT_th(klon)      REAL qT_th(klon)
130      ! T thermique niveau precedent      ! T thermique niveau precedent
     REAL Tbef(klon)  
131      REAL qsatbef(klon)      REAL qsatbef(klon)
132      ! le thermique est sature      ! le thermique est sature
133      LOGICAL Zsat(klon)      LOGICAL Zsat(klon)
134      ! Cape du thermique      REAL zthvd, zthvu, qqsat
135      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  
136    
137      ! inverse phi function for momentum      ! inverse phi function for momentum
138      REAL phiminv(klon)      REAL phiminv(klon)
     ! inverse phi function for heat  
     REAL phihinv(klon)  
139      ! turbulent velocity scale for momentum      ! turbulent velocity scale for momentum
140      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)  
141      ! current level height + one level up      ! current level height + one level up
142      REAL zp(klon)      REAL zp(klon)
143      REAL zcor, zdelta, zcvm5      REAL zcor
144    
145      REAL fac, pblmin, zmzp, term      REAL pblmin
146    
147      !-----------------------------------------------------------------      !-----------------------------------------------------------------
148    
149        knon = size(pblh)
150    
151      ! initialisations      ! initialisations
152      q_star = 0      q_star = 0
153      t_star = 0      t_star = 0
# Line 212  contains Line 155  contains
155      b212=sqrt(b1*b2)      b212=sqrt(b1*b2)
156      b2sr=sqrt(b2)      b2sr=sqrt(b2)
157    
     ! Initialisation  
     RLvCp = RLVTT/RCPD  
     REPS = RD/RV  
   
158      ! Calculer les hauteurs de chaque couche      ! Calculer les hauteurs de chaque couche
159      ! (geopotentielle Int_dp/ro = Int_[Rd.T.dp/p] z = geop/g)      ! (geopotentielle Int_dp/ro = Int_[Rd.T.dp/p] z = geop/g)
160      ! pourquoi ne pas utiliser Phi/RG ?      ! pourquoi ne pas utiliser Phi/RG ?
# Line 246  contains Line 185  contains
185         zxt = t2m(i)         zxt = t2m(i)
186    
187         ! convention >0 vers le bas ds lmdz         ! convention >0 vers le bas ds lmdz
188         khfs(i) = - flux_t(i, 1)*zxt*Rd / (RCPD*paprs(i, 1))         khfs(i) = - flux_t(i)*zxt*Rd / (RCPD*paprs(i, 1))
189         kqfs(i) = - flux_q(i, 1)*zxt*Rd / (paprs(i, 1))         kqfs(i) = - flux_q(i)*zxt*Rd / paprs(i, 1)
190         ! verifier que khfs et kqfs sont bien de la forme w'l'         ! verifier que khfs et kqfs sont bien de la forme w'l'
191         heatv(i) = khfs(i) + 0.608*zxt*kqfs(i)         heatv(i) = khfs(i) + 0.608*zxt*kqfs(i)
192         ! 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 266  contains Line 205  contains
205         plcl(i) = 6000.         plcl(i) = 6000.
206         ! Lambda = -u*^3 / (alpha.g.kvon.<w'Theta'v>         ! Lambda = -u*^3 / (alpha.g.kvon.<w'Theta'v>
207         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.  
208      ENDDO      ENDDO
209    
210      ! PBL height calculation: Search for level of pbl. Scan upward      ! PBL height calculation: Search for level of pbl. Scan upward
211      ! until the Richardson number between the first level and the      ! until the Richardson number between the first level and the
212      ! current level exceeds the "critical" value.  (bonne idee Nu de      ! current level exceeds the "critical" value.  (bonne idee Nu de
213      ! separer le Ric et l'exces de temp du thermique)      ! separer le Ric et l'exces de temp du thermique)
     fac = 100.  
214      DO k = 2, isommet      DO k = 2, isommet
215         DO i = 1, knon         DO i = 1, knon
216            IF (check(i)) THEN            IF (check(i)) THEN
# Line 340  contains Line 275  contains
275            q_star = kqfs(i)/wm(i)            q_star = kqfs(i)/wm(i)
276            t_star = khfs(i)/wm(i)            t_star = khfs(i)/wm(i)
277    
           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  
   
278            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 &
279                 + (RETV*T2m(i))**2*b2*q_star**2 &                 + (RETV*T2m(i))**2*b2*q_star**2 &
280                 + max(0., 2.*RETV*T2m(i)*b212*q_star*t_star))                 + max(0., 2.*RETV*T2m(i)*b212*q_star*t_star))
# Line 420  contains Line 350  contains
350         ! omegafl utilise pour prolongement CAPE         ! omegafl utilise pour prolongement CAPE
351         omegafl(i) = .FALSE.         omegafl(i) = .FALSE.
352         Cape(i) = 0.         Cape(i) = 0.
        Kape(i) = 0.  
353         EauLiq(i) = 0.         EauLiq(i) = 0.
354         CTEI(i) = 0.         CTEI(i) = 0.
        pblk(i) = 0.0  
        fak1(i) = ustar(i)*pblh(i)*vk  
355    
356         ! Do additional preparation for unstable cases only, set temperature         ! Do additional preparation for unstable cases only, set temperature
357         ! and moisture perturbations depending on stability.         ! and moisture perturbations depending on stability.
# Line 434  contains Line 361  contains
361            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))) &
362                 *(1.+RETV*qT_th(i))                 *(1.+RETV*qT_th(i))
363            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))  
364            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)  
365         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)  
366      ENDDO      ENDDO
367    
368      ! Main level loop to compute the diffusivities and      ! Main level loop to compute the diffusivities and
369      ! counter-gradient terms:      ! counter-gradient terms:
370      DO k = 2, isommet      loop_level: DO k = 2, isommet
371         ! Find levels within boundary layer:         ! Find levels within boundary layer:
372         DO i = 1, knon         DO i = 1, knon
           unslev(i) = .FALSE.  
           stblev(i) = .FALSE.  
           zm(i) = z(i, k-1)  
373            zp(i) = z(i, k)            zp(i) = z(i, k)
374            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  
375         ENDDO         ENDDO
376    
377         ! For all layers, compute integral info and CTEI         ! For all layers, compute integral info and CTEI
378         DO i = 1, knon         DO i = 1, knon
379            if (check(i).or.omegafl(i)) then            if (check(i) .or. omegafl(i)) then
380               if (.not.Zsat(i)) then               if (.not. Zsat(i)) then
381                  T2 = T2m(i) * s(i, k)                  T2 = T2m(i) * s(i, k)
382                  ! thermodyn functions                  ! thermodyn functions
383                  zdelta=MAX(0., SIGN(1., RTT - T2))                  qqsat= r2es * FOEEW(T2, RTT >= T2) / pplay(i, k)
                 qqsat= r2es * FOEEW(T2, zdelta) / pplay(i, k)  
384                  qqsat=MIN(0.5, qqsat)                  qqsat=MIN(0.5, qqsat)
385                  zcor=1./(1.-retv*qqsat)                  zcor=1./(1.-retv*qqsat)
386                  qqsat=qqsat*zcor                  qqsat=qqsat*zcor
# Line 536  contains Line 394  contains
394                             * (qT_th(i)-qsatbef(i)) / (qsatbef(i)-qqsat)                             * (qT_th(i)-qsatbef(i)) / (qsatbef(i)-qqsat)
395                     endif                     endif
396                     Zsat(i) = .true.                     Zsat(i) = .true.
                    Tbef(i) = T2  
397                  endif                  endif
398               endif               endif
399               qsatbef(i) = qqsat               qsatbef(i) = qqsat
400               ! cette ligne a deja ete faite normalement ?               ! cette ligne a deja ete faite normalement ?
401            endif            endif
402         ENDDO         ENDDO
403      end DO      end DO loop_level
404    
405    END SUBROUTINE HBTM    END SUBROUTINE HBTM
406    
Legend:
Removed from v.37  
changed lines
  Added in v.286
  ViewVC Help
Powered by ViewVC 1.1.21