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