/[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/Sources/phylmd/hbtm.f revision 188 by guez, Tue Mar 22 16:31:39 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)  
   
     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        ! 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      ! 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 69  contains Line 80  contains
80      PARAMETER (vk=0.35)      PARAMETER (vk=0.35)
81      REAL ricr      REAL ricr
82      PARAMETER (ricr=0.4)      PARAMETER (ricr=0.4)
     REAL fak  
     ! b calcul du Prandtl et de dTetas  
     PARAMETER (fak=8.5)  
     REAL fakn  
83      ! a      ! a
     PARAMETER (fakn=7.2)  
84      REAL onet      REAL onet
85      PARAMETER (onet=1.0/3.0)      PARAMETER (onet=1.0/3.0)
     REAL t_coup  
     PARAMETER(t_coup=273.15)  
86      REAL zkmin      REAL zkmin
87      PARAMETER (zkmin=0.01)      PARAMETER (zkmin=0.01)
88      REAL betam      REAL betam
89      ! pour Phim / h dans la S.L stable      ! pour Phim / h dans la S.L stable
90      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)  
91      ! z/OBL<>1      ! z/OBL<>1
92      REAL sffrac      REAL sffrac
93      ! S.L. = z/h < .1      ! S.L. = z/h < .1
94      PARAMETER (sffrac=0.1)      PARAMETER (sffrac=0.1)
95      REAL binm      REAL binm
96      PARAMETER (binm=betam*sffrac)      PARAMETER (binm=betam*sffrac)
     REAL binh  
     PARAMETER (binh=betah*sffrac)  
     REAL ccon  
     PARAMETER (ccon=fak*sffrac*vk)  
97    
98      REAL q_star, t_star      REAL q_star, t_star
99      ! Lambert correlations T' q' avec T* q*      ! Lambert correlations T' q' avec T* q*
# Line 124  contains Line 119  contains
119      REAL rhino(klon, klev)      REAL rhino(klon, klev)
120      ! pts w/unstbl pbl (positive virtual ht flx)      ! pts w/unstbl pbl (positive virtual ht flx)
121      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)  
122      LOGICAL check(klon) ! Richardson number > critical      LOGICAL check(klon) ! Richardson number > critical
123      ! flag de prolongerment cape pour pt Omega      ! flag de prolongerment cape pour pt Omega
124      LOGICAL omegafl(klon)      LOGICAL omegafl(klon)
     REAL pblh(klon)  
     REAL pblT(klon)  
     REAL plcl(klon)  
125    
126      ! Monin-Obukhov lengh      ! Monin-Obukhov lengh
127      REAL obklen(klon)      REAL obklen(klon)
128    
129      REAL zdu2      REAL zdu2
     ! thermal virtual temperature excess  
     REAL therm(klon)  
     REAL trmb1(klon), trmb2(klon), trmb3(klon)  
130      ! Algorithme thermique      ! Algorithme thermique
131      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)  
132      ! total water of thermal      ! total water of thermal
133      REAL qT_th(klon)      REAL qT_th(klon)
134      ! T thermique niveau precedent      ! T thermique niveau precedent
     REAL Tbef(klon)  
135      REAL qsatbef(klon)      REAL qsatbef(klon)
136      ! le thermique est sature      ! le thermique est sature
137      LOGICAL Zsat(klon)      LOGICAL Zsat(klon)
138      ! Cape du thermique      REAL zthvd, zthvu, qqsat
139      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  
140    
141      ! inverse phi function for momentum      ! inverse phi function for momentum
142      REAL phiminv(klon)      REAL phiminv(klon)
     ! inverse phi function for heat  
     REAL phihinv(klon)  
143      ! turbulent velocity scale for momentum      ! turbulent velocity scale for momentum
144      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)  
145      ! current level height + one level up      ! current level height + one level up
146      REAL zp(klon)      REAL zp(klon)
147      REAL zcor, zdelta, zcvm5      REAL zcor
148    
149      REAL fac, pblmin, zmzp, term      REAL pblmin
150    
151      !-----------------------------------------------------------------      !-----------------------------------------------------------------
152    
# Line 212  contains Line 157  contains
157      b212=sqrt(b1*b2)      b212=sqrt(b1*b2)
158      b2sr=sqrt(b2)      b2sr=sqrt(b2)
159    
     ! Initialisation  
     RLvCp = RLVTT/RCPD  
     REPS = RD/RV  
   
160      ! Calculer les hauteurs de chaque couche      ! Calculer les hauteurs de chaque couche
161      ! (geopotentielle Int_dp/ro = Int_[Rd.T.dp/p] z = geop/g)      ! (geopotentielle Int_dp/ro = Int_[Rd.T.dp/p] z = geop/g)
162      ! pourquoi ne pas utiliser Phi/RG ?      ! pourquoi ne pas utiliser Phi/RG ?
# Line 247  contains Line 188  contains
188    
189         ! convention >0 vers le bas ds lmdz         ! convention >0 vers le bas ds lmdz
190         khfs(i) = - flux_t(i, 1)*zxt*Rd / (RCPD*paprs(i, 1))         khfs(i) = - flux_t(i, 1)*zxt*Rd / (RCPD*paprs(i, 1))
191         kqfs(i) = - flux_q(i, 1)*zxt*Rd / (paprs(i, 1))         kqfs(i) = - flux_q(i, 1)*zxt*Rd / paprs(i, 1)
192         ! verifier que khfs et kqfs sont bien de la forme w'l'         ! verifier que khfs et kqfs sont bien de la forme w'l'
193         heatv(i) = khfs(i) + 0.608*zxt*kqfs(i)         heatv(i) = khfs(i) + 0.608*zxt*kqfs(i)
194         ! 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 275  contains Line 216  contains
216      ! until the Richardson number between the first level and the      ! until the Richardson number between the first level and the
217      ! current level exceeds the "critical" value.  (bonne idee Nu de      ! current level exceeds the "critical" value.  (bonne idee Nu de
218      ! separer le Ric et l'exces de temp du thermique)      ! separer le Ric et l'exces de temp du thermique)
     fac = 100.  
219      DO k = 2, isommet      DO k = 2, isommet
220         DO i = 1, knon         DO i = 1, knon
221            IF (check(i)) THEN            IF (check(i)) THEN
# Line 340  contains Line 280  contains
280            q_star = kqfs(i)/wm(i)            q_star = kqfs(i)/wm(i)
281            t_star = khfs(i)/wm(i)            t_star = khfs(i)/wm(i)
282    
           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  
   
283            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 &
284                 + (RETV*T2m(i))**2*b2*q_star**2 &                 + (RETV*T2m(i))**2*b2*q_star**2 &
285                 + 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 355  contains
355         ! omegafl utilise pour prolongement CAPE         ! omegafl utilise pour prolongement CAPE
356         omegafl(i) = .FALSE.         omegafl(i) = .FALSE.
357         Cape(i) = 0.         Cape(i) = 0.
        Kape(i) = 0.  
358         EauLiq(i) = 0.         EauLiq(i) = 0.
359         CTEI(i) = 0.         CTEI(i) = 0.
        pblk(i) = 0.0  
        fak1(i) = ustar(i)*pblh(i)*vk  
360    
361         ! Do additional preparation for unstable cases only, set temperature         ! Do additional preparation for unstable cases only, set temperature
362         ! and moisture perturbations depending on stability.         ! and moisture perturbations depending on stability.
# Line 434  contains Line 366  contains
366            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))) &
367                 *(1.+RETV*qT_th(i))                 *(1.+RETV*qT_th(i))
368            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))  
369            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)  
370         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)  
371      ENDDO      ENDDO
372    
373      ! Main level loop to compute the diffusivities and      ! Main level loop to compute the diffusivities and
374      ! counter-gradient terms:      ! counter-gradient terms:
375      DO k = 2, isommet      loop_level: DO k = 2, isommet
376         ! Find levels within boundary layer:         ! Find levels within boundary layer:
377         DO i = 1, knon         DO i = 1, knon
           unslev(i) = .FALSE.  
           stblev(i) = .FALSE.  
           zm(i) = z(i, k-1)  
378            zp(i) = z(i, k)            zp(i) = z(i, k)
379            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  
380         ENDDO         ENDDO
381    
382         ! For all layers, compute integral info and CTEI         ! For all layers, compute integral info and CTEI
383         DO i = 1, knon         DO i = 1, knon
384            if (check(i).or.omegafl(i)) then            if (check(i) .or. omegafl(i)) then
385               if (.not.Zsat(i)) then               if (.not. Zsat(i)) then
386                  T2 = T2m(i) * s(i, k)                  T2 = T2m(i) * s(i, k)
387                  ! thermodyn functions                  ! thermodyn functions
388                  zdelta=MAX(0., SIGN(1., RTT - T2))                  qqsat= r2es * FOEEW(T2, RTT >= T2) / pplay(i, k)
                 qqsat= r2es * FOEEW(T2, zdelta) / pplay(i, k)  
389                  qqsat=MIN(0.5, qqsat)                  qqsat=MIN(0.5, qqsat)
390                  zcor=1./(1.-retv*qqsat)                  zcor=1./(1.-retv*qqsat)
391                  qqsat=qqsat*zcor                  qqsat=qqsat*zcor
# Line 536  contains Line 399  contains
399                             * (qT_th(i)-qsatbef(i)) / (qsatbef(i)-qqsat)                             * (qT_th(i)-qsatbef(i)) / (qsatbef(i)-qqsat)
400                     endif                     endif
401                     Zsat(i) = .true.                     Zsat(i) = .true.
                    Tbef(i) = T2  
402                  endif                  endif
403               endif               endif
404               qsatbef(i) = qqsat               qsatbef(i) = qqsat
405               ! cette ligne a deja ete faite normalement ?               ! cette ligne a deja ete faite normalement ?
406            endif            endif
407         ENDDO         ENDDO
408      end DO      end DO loop_level
409    
410    END SUBROUTINE HBTM    END SUBROUTINE HBTM
411    
Legend:
Removed from v.37  
changed lines
  Added in v.188
  ViewVC Help
Powered by ViewVC 1.1.21