/[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 49 by guez, Wed Aug 24 11:43:14 2011 UTC trunk/Sources/phylmd/hbtm.f revision 186 by guez, Mon Mar 21 15:36:26 2016 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(knon, paprs, pplay, t2m, q2m, ustar, flux_t, flux_q, u, v, &
8         flux_q, u, v, t, q, pblh, cape, EauLiq, ctei, pblT, therm, trmb1, &         t, q, pblh, cape, EauLiq, ctei, pblT, therm, trmb1, trmb2, trmb3, 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
     ! 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      USE dimphy, ONLY: klev, klon, max      ! Algorithme th\'ese Anne Mathieu. Crit\'ere d'entra\^inement
14      USE suphec_m, ONLY: rcpd, rd, retv, rg, rkappa, rlvtt, rtt, rv      ! 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      USE yoethf_m, ONLY: r2es, rvtmp2
34      USE fcttre, ONLY: foeew      USE fcttre, ONLY: foeew
35    
     REAL RLvCp, REPS  
36      ! Arguments:      ! Arguments:
37    
38      ! nombre de points a calculer      ! nombre de points a calculer
39      INTEGER, intent(in):: knon      INTEGER, intent(in):: knon
40    
     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)  
41      ! pression a inter-couche (Pa)      ! pression a inter-couche (Pa)
42      REAL paprs(klon, klev+1)      REAL, intent(in):: paprs(klon, klev+1)
43      ! pression au milieu de couche (Pa)      ! pression au milieu de couche (Pa)
44      REAL pplay(klon, klev)      REAL, intent(in):: pplay(klon, klev)
45        REAL, intent(in):: t2m(klon) ! temperature a 2 m
46        ! q a 2 et 10m
47        REAL, intent(in):: q2m(klon)
48        REAL, intent(in):: ustar(klon)
49      ! Flux      ! Flux
50      REAL flux_t(klon, klev), flux_q(klon, klev)      REAL, intent(in):: flux_t(klon, klev), flux_q(klon, klev)
51      ! vitesse U (m/s)      ! vitesse U (m/s)
52      REAL u(klon, klev)      REAL, intent(in):: u(klon, klev)
53      ! vitesse V (m/s)      ! vitesse V (m/s)
54      REAL v(klon, klev)      REAL, intent(in):: v(klon, klev)
55      ! temperature (K)      ! temperature (K)
56      REAL t(klon, klev)      REAL, intent(in):: t(klon, klev)
57      ! vapeur d'eau (kg/kg)      ! vapeur d'eau (kg/kg)
58      REAL q(klon, klev)      REAL, intent(in):: q(klon, klev)
59    
60        REAL, intent(out):: pblh(:) ! (knon)
61        ! Cape du thermique
62        REAL Cape(klon)
63        ! Eau liqu integr du thermique
64        REAL EauLiq(klon)
65        ! Critere d'instab d'entrainmt des nuages de
66        REAL ctei(klon)
67        REAL pblT(klon)
68        ! thermal virtual temperature excess
69        REAL therm(klon)
70        REAL trmb1(klon), trmb2(klon), trmb3(klon)
71        REAL plcl(klon)
72    
73        ! Local:
74    
75      INTEGER isommet      INTEGER isommet
76      ! limite max sommet pbl      ! limite max sommet pbl
# Line 123  contains Line 135  contains
135      REAL rhino(klon, klev)      REAL rhino(klon, klev)
136      ! pts w/unstbl pbl (positive virtual ht flx)      ! pts w/unstbl pbl (positive virtual ht flx)
137      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)  
138      LOGICAL check(klon) ! Richardson number > critical      LOGICAL check(klon) ! Richardson number > critical
139      ! flag de prolongerment cape pour pt Omega      ! flag de prolongerment cape pour pt Omega
140      LOGICAL omegafl(klon)      LOGICAL omegafl(klon)
     REAL pblh(klon)  
     REAL pblT(klon)  
     REAL plcl(klon)  
141    
142      ! Monin-Obukhov lengh      ! Monin-Obukhov lengh
143      REAL obklen(klon)      REAL obklen(klon)
144    
145      REAL zdu2      REAL zdu2
     ! thermal virtual temperature excess  
     REAL therm(klon)  
     REAL trmb1(klon), trmb2(klon), trmb3(klon)  
146      ! Algorithme thermique      ! Algorithme thermique
147      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)  
148      ! total water of thermal      ! total water of thermal
149      REAL qT_th(klon)      REAL qT_th(klon)
150      ! T thermique niveau precedent      ! T thermique niveau precedent
     REAL Tbef(klon)  
151      REAL qsatbef(klon)      REAL qsatbef(klon)
152      ! le thermique est sature      ! le thermique est sature
153      LOGICAL Zsat(klon)      LOGICAL Zsat(klon)
154      ! Cape du thermique      REAL zthvd, zthvu, qqsat
155      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  
156    
157      ! inverse phi function for momentum      ! inverse phi function for momentum
158      REAL phiminv(klon)      REAL phiminv(klon)
     ! inverse phi function for heat  
     REAL phihinv(klon)  
159      ! turbulent velocity scale for momentum      ! turbulent velocity scale for momentum
160      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)  
161      ! current level height + one level up      ! current level height + one level up
162      REAL zp(klon)      REAL zp(klon)
163      REAL zcor, zdelta, zcvm5      REAL zcor
164    
165      REAL fac, pblmin, zmzp, term      REAL pblmin
166    
167      !-----------------------------------------------------------------      !-----------------------------------------------------------------
168    
# Line 211  contains Line 173  contains
173      b212=sqrt(b1*b2)      b212=sqrt(b1*b2)
174      b2sr=sqrt(b2)      b2sr=sqrt(b2)
175    
     ! Initialisation  
     RLvCp = RLVTT/RCPD  
     REPS = RD/RV  
   
176      ! Calculer les hauteurs de chaque couche      ! Calculer les hauteurs de chaque couche
177      ! (geopotentielle Int_dp/ro = Int_[Rd.T.dp/p] z = geop/g)      ! (geopotentielle Int_dp/ro = Int_[Rd.T.dp/p] z = geop/g)
178      ! pourquoi ne pas utiliser Phi/RG ?      ! pourquoi ne pas utiliser Phi/RG ?
# Line 246  contains Line 204  contains
204    
205         ! convention >0 vers le bas ds lmdz         ! convention >0 vers le bas ds lmdz
206         khfs(i) = - flux_t(i, 1)*zxt*Rd / (RCPD*paprs(i, 1))         khfs(i) = - flux_t(i, 1)*zxt*Rd / (RCPD*paprs(i, 1))
207         kqfs(i) = - flux_q(i, 1)*zxt*Rd / (paprs(i, 1))         kqfs(i) = - flux_q(i, 1)*zxt*Rd / paprs(i, 1)
208         ! verifier que khfs et kqfs sont bien de la forme w'l'         ! verifier que khfs et kqfs sont bien de la forme w'l'
209         heatv(i) = khfs(i) + 0.608*zxt*kqfs(i)         heatv(i) = khfs(i) + 0.608*zxt*kqfs(i)
210         ! 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 274  contains Line 232  contains
232      ! until the Richardson number between the first level and the      ! until the Richardson number between the first level and the
233      ! current level exceeds the "critical" value.  (bonne idee Nu de      ! current level exceeds the "critical" value.  (bonne idee Nu de
234      ! separer le Ric et l'exces de temp du thermique)      ! separer le Ric et l'exces de temp du thermique)
     fac = 100.  
235      DO k = 2, isommet      DO k = 2, isommet
236         DO i = 1, knon         DO i = 1, knon
237            IF (check(i)) THEN            IF (check(i)) THEN
# Line 339  contains Line 296  contains
296            q_star = kqfs(i)/wm(i)            q_star = kqfs(i)/wm(i)
297            t_star = khfs(i)/wm(i)            t_star = khfs(i)/wm(i)
298    
           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  
   
299            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 &
300                 + (RETV*T2m(i))**2*b2*q_star**2 &                 + (RETV*T2m(i))**2*b2*q_star**2 &
301                 + max(0., 2.*RETV*T2m(i)*b212*q_star*t_star))                 + max(0., 2.*RETV*T2m(i)*b212*q_star*t_star))
# Line 419  contains Line 371  contains
371         ! omegafl utilise pour prolongement CAPE         ! omegafl utilise pour prolongement CAPE
372         omegafl(i) = .FALSE.         omegafl(i) = .FALSE.
373         Cape(i) = 0.         Cape(i) = 0.
        Kape(i) = 0.  
374         EauLiq(i) = 0.         EauLiq(i) = 0.
375         CTEI(i) = 0.         CTEI(i) = 0.
        pblk(i) = 0.0  
        fak1(i) = ustar(i)*pblh(i)*vk  
376    
377         ! Do additional preparation for unstable cases only, set temperature         ! Do additional preparation for unstable cases only, set temperature
378         ! and moisture perturbations depending on stability.         ! and moisture perturbations depending on stability.
# Line 433  contains Line 382  contains
382            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))) &
383                 *(1.+RETV*qT_th(i))                 *(1.+RETV*qT_th(i))
384            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))  
385            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)  
386         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)  
387      ENDDO      ENDDO
388    
389      ! Main level loop to compute the diffusivities and      ! Main level loop to compute the diffusivities and
# Line 449  contains Line 391  contains
391      loop_level: DO k = 2, isommet      loop_level: DO k = 2, isommet
392         ! Find levels within boundary layer:         ! Find levels within boundary layer:
393         DO i = 1, knon         DO i = 1, knon
           unslev(i) = .FALSE.  
           stblev(i) = .FALSE.  
           zm(i) = z(i, k-1)  
394            zp(i) = z(i, k)            zp(i) = z(i, k)
395            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  
396         ENDDO         ENDDO
397    
398         ! For all layers, compute integral info and CTEI         ! For all layers, compute integral info and CTEI
# Line 520  contains Line 401  contains
401               if (.not. Zsat(i)) then               if (.not. Zsat(i)) then
402                  T2 = T2m(i) * s(i, k)                  T2 = T2m(i) * s(i, k)
403                  ! thermodyn functions                  ! thermodyn functions
404                  zdelta=MAX(0., SIGN(1., RTT - T2))                  qqsat= r2es * FOEEW(T2, RTT >= T2) / pplay(i, k)
                 qqsat= r2es * FOEEW(T2, zdelta) / pplay(i, k)  
405                  qqsat=MIN(0.5, qqsat)                  qqsat=MIN(0.5, qqsat)
406                  zcor=1./(1.-retv*qqsat)                  zcor=1./(1.-retv*qqsat)
407                  qqsat=qqsat*zcor                  qqsat=qqsat*zcor
# Line 535  contains Line 415  contains
415                             * (qT_th(i)-qsatbef(i)) / (qsatbef(i)-qqsat)                             * (qT_th(i)-qsatbef(i)) / (qsatbef(i)-qqsat)
416                     endif                     endif
417                     Zsat(i) = .true.                     Zsat(i) = .true.
                    Tbef(i) = T2  
418                  endif                  endif
419               endif               endif
420               qsatbef(i) = qqsat               qsatbef(i) = qqsat
Legend:
Removed from v.49  
changed lines
  Added in v.186
  ViewVC Help
Powered by ViewVC 1.1.21