source: LMDZ6/branches/Amaury_dev/libf/phylmd/hbtm_mod.F90 @ 5137

Last change on this file since 5137 was 5119, checked in by abarral, 15 months ago

enforce PRIVATE by default in several modules, expose PUBLIC as needed
move eigen.f90 to obsolete/
(lint) aslong the way
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 29.3 KB
Line 
1module hbtm_mod
2
3  IMPLICIT NONE
4
5CONTAINS
6
7  SUBROUTINE hbtm(knon, paprs, pplay, t2m, t10m, q2m, q10m, ustar, wstar, &
8       flux_t, flux_q, u, v, t, q, pblh, cape, eauliq, ctei, pblt, therm, &
9       trmb1, trmb2, trmb3, plcl)
10    USE dimphy
11
12    ! ***************************************************************
13    ! *                                                             *
14    ! * HBTM2   D'apres Holstag&Boville et Troen&Mahrt              *
15    ! *                 JAS 47              BLM                     *
16    ! * Algorithme These Anne Mathieu                               *
17    ! * Critere d'Entrainement Peter Duynkerke (JAS 50)             *
18    ! * written by  : Anne MATHIEU & Alain LAHELLEC, 22/11/99       *
19    ! * features : implem. exces Mathieu                            *
20    ! ***************************************************************
21    ! * mods : decembre 99 passage th a niveau plus bas. voir fixer *
22    ! * la prise du th a z/Lambda = -.2 (max Ray)                   *
23    ! * Autre algo : entrainement ~ Theta+v =cste mais comment=>The?*
24    ! * on peut fixer q a .7qsat(cf non adiab)=>T2 et The2          *
25    ! * voir aussi //KE pblh = niveau The_e ou l = env.             *
26    ! ***************************************************************
27    ! * fin therm a la HBTM passage a forme Mathieu 12/09/2001      *
28    ! ***************************************************************
29    ! *
30
31
32    ! AM Fev 2003
33    ! Adaptation a LMDZ version couplee
34
35    ! Pour le moment on fait passer en argument les grdeurs de surface :
36    ! flux, t,q2m, t,q10m, on va utiliser systematiquement les grdeurs a 2m ms
37    ! on garde la possibilite de changer si besoin est (jusqu'a present la
38    ! forme de HB avec le 1er niveau modele etait conservee)
39
40
41
42
43
44    include "YOMCST.h"
45    REAL rlvcp, reps
46    ! Arguments:
47
48    INTEGER knon ! nombre de points a calculer
49    ! AM
50    REAL t2m(klon), t10m(klon) ! temperature a 2 et 10m
51    REAL q2m(klon), q10m(klon) ! q a 2 et 10m
52    REAL ustar(klon)
53    REAL wstar(klon) ! w*, convective velocity scale
54    REAL paprs(klon, klev+1) ! pression a inter-couche (Pa)
55    REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
56    REAL flux_t(klon, klev), flux_q(klon, klev) ! Flux
57    REAL u(klon, klev) ! vitesse U (m/s)
58    REAL v(klon, klev) ! vitesse V (m/s)
59    REAL t(klon, klev) ! temperature (K)
60    REAL q(klon, klev) ! vapeur d'eau (kg/kg)
61    ! AM      REAL cd_h(klon) ! coefficient de friction au sol pour chaleur
62    ! AM      REAL cd_m(klon) ! coefficient de friction au sol pour vitesse
63
64    INTEGER isommet
65    ! um      PARAMETER (isommet=klev) ! limite max sommet pbl
66    REAL, PARAMETER :: vk = 0.35 ! Von Karman => passer a .41 ! cf U.Olgstrom
67    REAL, PARAMETER :: ricr = 0.4
68    REAL, PARAMETER :: fak = 8.5 ! b calcul du Prandtl et de dTetas
69    REAL, PARAMETER :: fakn = 7.2 ! a
70    REAL, PARAMETER :: onet = 1.0/3.0
71    REAL, PARAMETER :: t_coup = 273.15
72    REAL, PARAMETER :: zkmin = 0.01
73    REAL, PARAMETER :: betam = 15.0 ! pour Phim / h dans la S.L stable
74    REAL, PARAMETER :: betah = 15.0
75
76    REAL, PARAMETER :: betas = 5.0
77    ! Phit dans la S.L. stable (mais 2 formes / z/OBL<>1
78
79    REAL, PARAMETER :: sffrac = 0.1 ! S.L. = z/h < .1
80    REAL, PARAMETER :: usmin = 1.E-12
81    REAL, PARAMETER :: binm = betam*sffrac
82    REAL, PARAMETER :: binh = betah*sffrac
83    REAL, PARAMETER :: ccon = fak*sffrac*vk
84    REAL, PARAMETER :: b1 = 70., b2 = 20.
85    REAL, PARAMETER :: zref = 2. ! Niveau de ref a 2m peut eventuellement
86    ! etre choisi a 10m
87    REAL q_star, t_star
88    REAL b212, b2sr ! Lambert correlations T' q' avec T* q*
89
90    REAL z(klon, klev)
91    ! AM      REAL pcfm(klon,klev), pcfh(klon,klev)
92    INTEGER i, k, j
93    REAL zxt
94    ! AM      REAL zxt, zxq, zxu, zxv, zxmod, taux, tauy
95    ! AM      REAL zx_alf1, zx_alf2 ! parametres pour extrapolation
96    REAL khfs(klon) ! surface kinematic heat flux [mK/s]
97    REAL kqfs(klon) ! sfc kinematic constituent flux [m/s]
98    REAL heatv(klon) ! surface virtual heat flux
99    REAL rhino(klon, klev) ! bulk Richardon no. mais en Theta_v
100    LOGICAL unstbl(klon) ! pts w/unstbl pbl (positive virtual ht flx)
101    LOGICAL stblev(klon) ! stable pbl with levels within pbl
102    LOGICAL unslev(klon) ! unstbl pbl with levels within pbl
103    LOGICAL unssrf(klon) ! unstb pbl w/lvls within srf pbl lyr
104    LOGICAL unsout(klon) ! unstb pbl w/lvls in outer pbl lyr
105    LOGICAL check(klon) ! True=>chk if Richardson no.>critcal
106    LOGICAL omegafl(klon) ! flag de prolongerment cape pour pt Omega
107    REAL pblh(klon)
108    REAL pblt(klon)
109    REAL plcl(klon)
110    ! AM      REAL cgh(klon,2:klev) ! counter-gradient term for heat [K/m]
111    ! AM      REAL cgq(klon,2:klev) ! counter-gradient term for constituents
112    ! AM      REAL cgs(klon,2:klev) ! counter-gradient star (cg/flux)
113    REAL unsobklen(klon) ! Monin-Obukhov lengh
114    ! AM      REAL ztvd, ztvu,
115    REAL zdu2
116    REAL, INTENT(OUT):: therm(:) ! (klon) thermal virtual temperature excess
117    REAL trmb1(klon), trmb2(klon), trmb3(klon)
118    ! Algorithme thermique
119    REAL s(klon, klev) ! [P/Po]^Kappa milieux couches
120    REAL th_th(klon) ! potential temperature of thermal
121    REAL the_th(klon) ! equivalent potential temperature of thermal
122    REAL qt_th(klon) ! total water  of thermal
123    REAL tbef(klon) ! T thermique niveau precedent
124    REAL qsatbef(klon)
125    LOGICAL zsat(klon) ! le thermique est sature
126    REAL cape(klon) ! Cape du thermique
127    REAL kape(klon) ! Cape locale
128    REAL eauliq(klon) ! Eau liqu integr du thermique
129    REAL ctei(klon) ! Critere d'instab d'entrainmt des nuages de CL
130    REAL the1, the2, aa, bb, zthvd, zthvu, xintpos, qqsat
131    ! IM 091204 BEG
132    REAL a1, a2, a3
133    ! IM 091204 END
134    REAL xhis, rnum, denom, th1, th2, thv1, thv2, ql2
135    REAL dqsat_dt, qsat2, qt1, q2, t1, t2, xnull, delt_the
136    REAL delt_qt, delt_2, quadsat, spblh, reduc
137
138    REAL phiminv(klon) ! inverse phi function for momentum
139    REAL phihinv(klon) ! inverse phi function for heat
140    REAL wm(klon) ! turbulent velocity scale for momentum
141    REAL fak1(klon) ! k*ustar*pblh
142    REAL fak2(klon) ! k*wm*pblh
143    REAL fak3(klon) ! fakn*wstar/wm
144    REAL pblk(klon) ! level eddy diffusivity for momentum
145    REAL pr(klon) ! Prandtl number for eddy diffusivities
146    REAL zl(klon) ! zmzp / Obukhov length
147    REAL zh(klon) ! zmzp / pblh
148    REAL zzh(klon) ! (1-(zmzp/pblh))**2
149    REAL zm(klon) ! current level height
150    REAL zp(klon) ! current level height + one level up
151    REAL zcor, zdelta, zcvm5
152    ! AM      REAL zxqs
153    REAL fac, pblmin, zmzp, term
154
155    include "YOETHF.h"
156    include "FCTTRE.h"
157
158
159
160    ! initialisations (Anne)
161    isommet = klev
162    th_th(:) = 0.
163    q_star = 0
164    t_star = 0
165    therm = 0.
166
167    b212 = sqrt(b1*b2)
168    b2sr = sqrt(b2)
169
170    ! ============================================================
171    ! Fonctions thermo implicites
172    ! ============================================================
173    ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
174    ! Tetens : pression partielle de vap d'eau e_sat(T)
175    ! =================================================
176    ! ++ e_sat(T) = r2*exp( r3*(T-Tf)/(T-r4) ) id a r2*FOEWE
177    ! ++ avec :
178    ! ++ Tf = 273.16 K  (Temp de fusion de la glace)
179    ! ++ r2 = 611.14 Pa
180    ! ++ r3 = 17.269 (liquide) 21.875 (solide) adim
181    ! ++ r4 = 35.86             7.66           Kelvin
182    ! ++  q_sat = eps*e_sat/(p-(1-eps)*e_sat)
183    ! ++ derivé :
184    ! ++ =========
185    ! ++                   r3*(Tf-r4)*q_sat(T,p)
186    ! ++ d_qsat_dT = --------------------------------
187    ! ++             (T-r4)^2*( 1-(1-eps)*e_sat(T)/p )
188    ! ++ pour zcvm5=Lv, c'est FOEDE
189    ! ++ Rq :(1.-REPS)*esarg/Parg id a RETV*Qsat
190    ! ------------------------------------------------------------------
191
192    ! Initialisation
193    rlvcp = rlvtt/rcpd
194    reps = rd/rv
195
196
197    ! DO i = 1, klon
198    ! pcfh(i,1) = cd_h(i)
199    ! pcfm(i,1) = cd_m(i)
200    ! ENDDO
201    ! DO k = 2, klev
202    ! DO i = 1, klon
203    ! pcfh(i,k) = zkmin
204    ! pcfm(i,k) = zkmin
205    ! cgs(i,k) = 0.0
206    ! cgh(i,k) = 0.0
207    ! cgq(i,k) = 0.0
208    ! ENDDO
209    ! ENDDO
210
211    ! Calculer les hauteurs de chaque couche
212    ! (geopotentielle Int_dp/ro = Int_[Rd.T.dp/p] z = geop/g)
213    ! pourquoi ne pas utiliser Phi/RG ?
214    DO i = 1, knon
215       z(i, 1) = rd*t(i, 1)/(0.5*(paprs(i,1)+pplay(i,1)))&
216            *(paprs(i,1)-pplay(i,1))/rg
217       s(i, 1) = (pplay(i,1)/paprs(i,1))**rkappa
218    END DO
219    ! s(k) = [pplay(k)/ps]^kappa
220    ! + + + + + + + + + pplay  <-> s(k)   t  dp=pplay(k-1)-pplay(k)
221
222    ! -----------------  paprs <-> sig(k)
223
224    ! + + + + + + + + + pplay  <-> s(k-1)
225
226
227    ! + + + + + + + + + pplay  <-> s(1)   t  dp=paprs-pplay   z(1)
228
229    ! -----------------  paprs <-> sig(1)
230
231    DO k = 2, klev
232       DO i = 1, knon
233          z(i, k) = z(i, k-1) + rd*0.5*(t(i,k-1)+t(i,k))/paprs(i, k)&
234               *(pplay(i,k-1)-pplay(i,k))/rg
235          s(i, k) = (pplay(i,k)/paprs(i,1))**rkappa
236       END DO
237    END DO
238    ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
239    ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
240    ! +++  Determination des grandeurs de surface  +++++++++++++++++++++
241    ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
242    ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
243    DO i = 1, knon
244       ! AM         IF (thermcep) THEN
245       ! AM           zdelta=MAX(0.,SIGN(1.,RTT-tsol(i)))
246       ! zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
247       ! zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q(i,1))
248       ! AM           zxqs= r2es * FOEEW(tsol(i),zdelta)/paprs(i,1)
249       ! AM           zxqs=MIN(0.5,zxqs)
250       ! AM           zcor=1./(1.-retv*zxqs)
251       ! AM           zxqs=zxqs*zcor
252       ! AM         ELSE
253       ! AM           IF (tsol(i).LT.t_coup) THEN
254       ! AM              zxqs = qsats(tsol(i)) / paprs(i,1)
255       ! AM           ELSE
256       ! AM              zxqs = qsatl(tsol(i)) / paprs(i,1)
257       ! AM           ENDIF
258       ! AM         ENDIF
259       ! niveau de reference bulk; mais ici, c,a pourrait etre le niveau de ref
260       ! du thermique
261       ! AM        zx_alf1 = 1.0
262       ! AM        zx_alf2 = 1.0 - zx_alf1
263       ! AM        zxt = (t(i,1)+z(i,1)*RG/RCPD/(1.+RVTMP2*q(i,1)))
264       ! AM     .        *(1.+RETV*q(i,1))*zx_alf1
265       ! AM     .      + (t(i,2)+z(i,2)*RG/RCPD/(1.+RVTMP2*q(i,2)))
266       ! AM     .        *(1.+RETV*q(i,2))*zx_alf2
267       ! AM        zxu = u(i,1)*zx_alf1+u(i,2)*zx_alf2
268       ! AM        zxv = v(i,1)*zx_alf1+v(i,2)*zx_alf2
269       ! AM        zxq = q(i,1)*zx_alf1+q(i,2)*zx_alf2
270       ! AM
271       ! AMAM           zxu = u10m(i)
272       ! AMAM           zxv = v10m(i)
273       ! AMAM           zxmod = 1.0+SQRT(zxu**2+zxv**2)
274       ! AM Niveau de ref choisi a 2m
275       zxt = t2m(i)
276
277       ! ***************************************************
278       ! attention, il doit s'agir de <w'theta'>
279       ! ;Calcul de tcls virtuel et de w'theta'virtuel
280       ! ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
281       ! tcls=tcls*(1+.608*qcls)
282
283       ! ;Pour avoir w'theta',
284       ! ; il faut diviser par ro.Cp
285       ! Cp=Cpd*(1+0.84*qcls)
286       ! fcs=fcs/(ro_surf*Cp)
287       ! ;On transforme w'theta' en w'thetav'
288       ! Lv=(2.501-0.00237*(tcls-273.15))*1.E6
289       ! xle=xle/(ro_surf*Lv)
290       ! fcsv=fcs+.608*xle*tcls
291       ! ***************************************************
292       ! AM        khfs(i) = (tsol(i)*(1.+RETV*q(i,1))-zxt) *zxmod*cd_h(i)
293       ! AM        kqfs(i) = (zxqs-zxq) *zxmod*cd_h(i) * beta(i)
294       ! AM
295       ! dif khfs est deja w't'_v / heatv(i) = khfs(i) + RETV*zxt*kqfs(i)
296       ! AM calcule de Ro = paprs(i,1)/Rd zxt
297       ! AM convention >0 vers le bas ds lmdz
298       khfs(i) = -flux_t(i, 1)*zxt*rd/(rcpd*paprs(i,1))
299       kqfs(i) = -flux_q(i, 1)*zxt*rd/(paprs(i,1))
300       ! AM   verifier que khfs et kqfs sont bien de la forme w'l'
301       heatv(i) = khfs(i) + 0.608*zxt*kqfs(i)
302       ! a comparer aussi aux sorties de clqh : flux_T/RoCp et flux_q/RoLv
303       ! AM        heatv(i) = khfs(i)
304       ! AM ustar est en entree
305       ! AM        taux = zxu *zxmod*cd_m(i)
306       ! AM        tauy = zxv *zxmod*cd_m(i)
307       ! AM        ustar(i) = SQRT(taux**2+tauy**2)
308       ! AM        ustar(i) = MAX(SQRT(ustar(i)),0.01)
309       ! Theta et qT du thermique sans exces (interpolin vers surf)
310       ! chgt de niveau du thermique (jeudi 30/12/1999)
311       ! (interpolation lineaire avant integration phi_h)
312       ! AM        qT_th(i) = zxqs*beta(i) + 4./z(i,1)*(q(i,1)-zxqs*beta(i))
313       ! AM        qT_th(i) = max(qT_th(i),q(i,1))
314       qt_th(i) = q2m(i)
315       ! n The_th restera la Theta du thermique sans exces jusqu'a 2eme calcul
316       ! n reste a regler convention P) pour Theta
317       ! The_th(i) = tsol(i) + 4./z(i,1)*(t(i,1)-tsol(i))
318       ! -                      + RLvCp*qT_th(i)
319       ! AM        Th_th(i) = tsol(i) + 4./z(i,1)*(t(i,1)-tsol(i))
320       th_th(i) = t2m(i)
321    END DO
322
323    DO i = 1, knon
324       rhino(i, 1) = 0.0 ! Global Richardson
325       check(i) = .TRUE.
326       pblh(i) = z(i, 1) ! on initialise pblh a l'altitude du 1er niveau
327       plcl(i) = 6000.
328       ! Lambda = -u*^3 / (alpha.g.kvon.<w'Theta'v>
329       unsobklen(i) = -rg*vk*heatv(i)/(t(i,1)*max(ustar(i),usmin)**3)
330       trmb1(i) = 0.
331       trmb2(i) = 0.
332       trmb3(i) = 0.
333    END DO
334
335
336    ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
337    ! PBL height calculation:
338    ! Search for level of pbl. Scan upward until the Richardson number between
339    ! the first level and the current level exceeds the "critical" value.
340    ! (bonne idee Nu de separer le Ric et l'exces de temp du thermique)
341    ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
342    fac = 100.0
343    DO k = 2, isommet
344       DO i = 1, knon
345          IF (check(i)) THEN
346             ! pourquoi / niveau 1 (au lieu du sol) et le terme en u*^2 ?
347             ! test     zdu2 =
348             ! (u(i,k)-u(i,1))**2+(v(i,k)-v(i,1))**2+fac*ustar(i)**2
349             zdu2 = u(i, k)**2 + v(i, k)**2
350             zdu2 = max(zdu2, 1.0E-20)
351             ! Theta_v environnement
352             zthvd = t(i, k)/s(i, k)*(1.+retv*q(i,k))
353
354             ! therm Theta_v sans exces (avec hypothese fausse de H&B, sinon,
355             ! passer par Theta_e et virpot)
356             ! zthvu=t(i,1)/s(i,1)*(1.+RETV*q(i,1))
357             ! AM         zthvu = Th_th(i)*(1.+RETV*q(i,1))
358             zthvu = th_th(i)*(1.+retv*qt_th(i))
359             ! Le Ri par Theta_v
360             ! AM         rhino(i,k) = (z(i,k)-z(i,1))*RG*(zthvd-zthvu)
361             ! AM     .               /(zdu2*0.5*(zthvd+zthvu))
362             ! AM On a nveau de ref a 2m ???
363             rhino(i, k) = (z(i,k)-zref)*rg*(zthvd-zthvu)/(zdu2*0.5&
364                  *(zthvd+zthvu))
365
366             IF (rhino(i,k)>=ricr) THEN
367                pblh(i) = z(i, k-1) + (z(i,k-1)-z(i,k))*(ricr-rhino(i,k-1))&
368                     /(rhino(i,k-1)-rhino(i,k))
369                ! test04
370                pblh(i) = pblh(i) + 100.
371                pblt(i) = t(i, k-1) + (t(i,k)-t(i,k-1))*(pblh(i)-z(i,k-1))&
372                     /(z(i,k)- z(i,k-1))
373                check(i) = .FALSE.
374             END IF
375          END IF
376       END DO
377    END DO
378
379
380    ! Set pbl height to maximum value where computation exceeds number of
381    ! layers allowed
382
383    DO i = 1, knon
384       IF (check(i)) pblh(i) = z(i, isommet)
385    END DO
386
387    ! Improve estimate of pbl height for the unstable points.
388    ! Find unstable points (sensible heat flux is upward):
389
390    DO i = 1, knon
391       IF (heatv(i)>0.) THEN
392          unstbl(i) = .TRUE.
393          check(i) = .TRUE.
394       ELSE
395          unstbl(i) = .FALSE.
396          check(i) = .FALSE.
397       END IF
398    END DO
399
400    ! For the unstable case, compute velocity scale and the
401    ! convective temperature excess:
402
403    DO i = 1, knon
404       IF (check(i)) THEN
405          phiminv(i) = (1.-binm*pblh(i)*unsobklen(i))**onet
406          ! ***************************************************
407          ! Wm ? et W* ? c'est la formule pour z/h < .1
408          ! ;Calcul de w* ;;
409          ! ;;;;;;;;;;;;;;;;
410          ! w_star=((g/tcls)*fcsv*z(ind))^(1/3.) [ou prendre la premiere approx
411          ! de h)
412          ! ;; CALCUL DE wm ;;
413          ! ;;;;;;;;;;;;;;;;;;
414          ! ; Ici on considerera que l'on est dans la couche de surf jusqu'a
415          ! 100 m
416          ! ; On prend svt couche de surface=0.1*h mais on ne connait pas h
417          ! ;;;;;;;;;;;Dans la couche de surface
418          ! if (z(ind) le 20) then begin
419          ! Phim=(1.-15.*(z(ind)/L))^(-1/3.)
420          ! wm=u_star/Phim
421          ! ;;;;;;;;;;;En dehors de la couche de surface
422          ! END IF ELSE IF (z(ind) gt 20) then begin
423          ! wm=(u_star^3+c1*w_star^3)^(1/3.)
424          ! END IF
425          ! ***************************************************
426          wm(i) = ustar(i)*phiminv(i)
427          ! ===================================================================
428          ! valeurs de Dominique Lambert de la campagne SEMAPHORE :
429          ! <T'^2> = 100.T*^2; <q'^2> = 20.q*^2 a 10m
430          ! <Tv'^2> = (1+1.2q).100.T* + 1.2Tv.sqrt(20*100).T*.q* +
431          ! (.608*Tv)^2*20.q*^2;
432          ! et dTetavS = sqrt(<Tv'^2>) ainsi calculee.
433          ! avec : T*=<w'T'>_s/w* et q*=<w'q'>/w*
434          ! !!! on peut donc utiliser w* pour les fluctuations <-> Lambert
435          ! (leur corellation pourrait dependre de beta par ex)
436          ! if fcsv(i,j) gt 0 then begin
437          ! dTetavs=b1*(1.+2.*.608*q_10(i,j))*(fcs(i,j)/wm(i,j))^2+$
438          ! (.608*Thetav_10(i,j))^2*b2*(xle(i,j)/wm(i,j))^2+$
439          ! 2.*.608*thetav_10(i,j)*sqrt(b1*b2)*(xle(i,j)/wm(i,j))*(fcs(i,j)
440          ! /wm(i,j))
441          ! dqs=b2*(xle(i,j)/wm(i,j))^2
442          ! theta_s(i,j)=thetav_10(i,j)+sqrt(dTetavs)
443          ! q_s(i,j)=q_10(i,j)+sqrt(dqs)
444          ! END IF else begin
445          ! Theta_s(i,j)=thetav_10(i,j)
446          ! q_s(i,j)=q_10(i,j)
447          ! endelse
448          ! ===================================================================
449
450          ! HBTM        therm(i) = heatv(i)*fak/wm(i)
451          ! forme Mathieu :
452          q_star = kqfs(i)/wm(i)
453          t_star = khfs(i)/wm(i)
454          ! IM 091204 BEG
455          IF (1==0) THEN
456             IF (t_star<0. .OR. q_star<0.) THEN
457                PRINT *, 'i t_star q_star khfs kqfs wm', i, t_star, q_star, &
458                     khfs(i), kqfs(i), wm(i)
459             END IF
460          END IF
461          ! IM 091204 END
462          ! AM Nveau cde ref 2m =>
463          ! AM        therm(i) = sqrt( b1*(1.+2.*RETV*q(i,1))*t_star**2
464          ! AM     +             + (RETV*T(i,1))**2*b2*q_star**2
465          ! AM     +             + 2.*RETV*T(i,1)*b212*q_star*t_star
466          ! AM     +                 )
467          ! IM 091204 BEG
468          a1 = b1*(1.+2.*retv*qt_th(i))*t_star**2
469          a2 = (retv*th_th(i))**2*b2*q_star*q_star
470          a3 = 2.*retv*th_th(i)*b212*q_star*t_star
471          aa = a1 + a2 + a3
472          IF (1==0) THEN
473             IF (aa<0.) THEN
474                PRINT *, 'i a1 a2 a3 aa', i, a1, a2, a3, aa
475                PRINT *, 'i qT_th Th_th t_star q_star RETV b1 b2 b212', i, &
476                     qt_th(i), th_th(i), t_star, q_star, retv, b1, b2, b212
477             END IF
478          END IF
479          ! IM 091204 END
480          therm(i) = sqrt(b1*(1.+2.*retv*qt_th(i))*t_star**2+(retv*th_th( &
481               i))**2*b2*q_star*q_star &  ! IM 101204  +             +
482                                ! 2.*RETV*Th_th(i)*b212*q_star*t_star
483               +max(0.,2.*retv*th_th(i)*b212*q_star*t_star))
484
485          ! Theta et qT du thermique (forme H&B) avec exces
486          ! (attention, on ajoute therm(i) qui est virtuelle ...)
487          ! pourquoi pas sqrt(b1)*t_star ?
488          ! dqs = b2sr*kqfs(i)/wm(i)
489          qt_th(i) = qt_th(i) + b2sr*q_star
490          ! new on differre le calcul de Theta_e
491          ! The_th(i) = The_th(i) + therm(i) + RLvCp*qT_th(i)
492          ! ou:    The_th(i) = The_th(i) + sqrt(b1)*khfs(i)/wm(i) +
493          ! RLvCp*qT_th(i)
494          rhino(i, 1) = 0.0
495       END IF
496    END DO
497
498    ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
499    ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
500    ! ++ Improve pblh estimate for unstable conditions using the +++++++
501    ! ++          convective temperature excess :                +++++++
502    ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
503    ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
504
505    DO k = 2, isommet
506       DO i = 1, knon
507          IF (check(i)) THEN
508             ! test     zdu2 =
509             ! (u(i,k)-u(i,1))**2+(v(i,k)-v(i,1))**2+fac*ustar(i)**2
510             zdu2 = u(i, k)**2 + v(i, k)**2
511             zdu2 = max(zdu2, 1.0E-20)
512             ! Theta_v environnement
513             zthvd = t(i, k)/s(i, k)*(1.+retv*q(i,k))
514
515             ! et therm Theta_v (avec hypothese de constance de H&B,
516             ! zthvu=(t(i,1)+therm(i))/s(i,1)*(1.+RETV*q(i,1))
517             zthvu = th_th(i)*(1.+retv*qt_th(i)) + therm(i)
518
519
520             ! Le Ri par Theta_v
521             ! AM Niveau de ref 2m
522             ! AM         rhino(i,k) = (z(i,k)-z(i,1))*RG*(zthvd-zthvu)
523             ! AM     .               /(zdu2*0.5*(zthvd+zthvu))
524             rhino(i, k) = (z(i,k)-zref)*rg*(zthvd-zthvu)/(zdu2*0.5&
525                  *(zthvd+zthvu))
526
527
528             IF (rhino(i,k)>=ricr) THEN
529                pblh(i) = z(i, k-1) + (z(i,k-1)-z(i,k))*(ricr-rhino(i,k-1))&
530                     /(rhino(i,k-1)-rhino(i,k))
531                ! test04
532                pblh(i) = pblh(i) + 100.
533                pblt(i) = t(i, k-1) + (t(i,k)-t(i,k-1))*(pblh(i)-z(i,k-1))&
534                     /(z(i,k)- z(i,k-1))
535                check(i) = .FALSE.
536                ! IM 170305 BEG
537                IF (1==0) THEN
538                   ! debug print -120;34       -34-        58 et    0;26 wamp
539                   IF (i==950 .OR. i==192 .OR. i==624 .OR. i==118) THEN
540                      PRINT *, ' i,Th_th,Therm,qT :', i, th_th(i), therm(i), &
541                           qt_th(i)
542                      q_star = kqfs(i)/wm(i)
543                      t_star = khfs(i)/wm(i)
544                      PRINT *, 'q* t*, b1,b2,b212 ', q_star, t_star, &
545                           b1*(1.+2.*retv*qt_th(i))*t_star**2, &
546                           (retv*th_th(i))**2*b2*q_star**2, 2.*retv*th_th(i)&
547                           *b212*q_star *t_star
548                      PRINT *, 'zdu2 ,100.*ustar(i)**2', zdu2, fac*ustar(i)**2
549                   END IF
550                END IF !(1.EQ.0) THEN
551                ! IM 170305 END
552                ! q_star = kqfs(i)/wm(i)
553                ! t_star = khfs(i)/wm(i)
554                ! trmb1(i) = b1*(1.+2.*RETV*q(i,1))*t_star**2
555                ! trmb2(i) = (RETV*T(i,1))**2*b2*q_star**2
556                ! Omega now   trmb3(i) = 2.*RETV*T(i,1)*b212*q_star*t_star
557             END IF
558          END IF
559       END DO
560    END DO
561
562    ! Set pbl height to maximum value where computation exceeds number of
563    ! layers allowed
564
565    DO i = 1, knon
566       IF (check(i)) pblh(i) = z(i, isommet)
567    END DO
568
569    ! PBL height must be greater than some minimum mechanical mixing depth
570    ! Several investigators have proposed minimum mechanical mixing depth
571    ! relationships as a function of the local friction velocity, u*.  We
572    ! make use of a linear relationship of the form h = c u* where c=700.
573    ! The scaling arguments that give rise to this relationship most often
574    ! represent the coefficient c as some constant over the local coriolis
575    ! parameter.  Here we make use of the experimental results of Koracin
576    ! and Berkowicz (1988) [BLM, Vol 43] for wich they recommend 0.07/f
577    ! where f was evaluated at 39.5 N and 52 N.  Thus we use a typical mid
578    ! latitude value for f so that c = 0.07/f = 700.
579
580    DO i = 1, knon
581       pblmin = 700.0*ustar(i)
582       pblh(i) = max(pblh(i), pblmin)
583       ! par exemple :
584       pblt(i) = t(i, 2) + (t(i,3)-t(i,2))*(pblh(i)-z(i,2))/(z(i,3)-z(i,2))
585    END DO
586
587    ! ********************************************************************
588    ! pblh is now available; do preparation for diffusivity calculation :
589    ! ********************************************************************
590    DO i = 1, knon
591       check(i) = .TRUE.
592       zsat(i) = .FALSE.
593       ! omegafl utilise pour prolongement CAPE
594       omegafl(i) = .FALSE.
595       cape(i) = 0.
596       kape(i) = 0.
597       eauliq(i) = 0.
598       ctei(i) = 0.
599       pblk(i) = 0.0
600       fak1(i) = ustar(i)*pblh(i)*vk
601
602       ! Do additional preparation for unstable cases only, set temperature
603       ! and moisture perturbations depending on stability.
604       ! *** Rq: les formule sont prises dans leur forme CS ***
605       IF (unstbl(i)) THEN
606          ! AM Niveau de ref du thermique
607          ! AM          zxt=(t(i,1)-z(i,1)*0.5*RG/RCPD/(1.+RVTMP2*q(i,1)))
608          ! AM     .         *(1.+RETV*q(i,1))
609          zxt = (th_th(i)-zref*0.5*rg/rcpd/(1.+rvtmp2*qt_th(i)))* &
610               (1.+retv*qt_th(i))
611          phiminv(i) = (1.-binm*pblh(i)*unsobklen(i))**onet
612          phihinv(i) = sqrt(1.-binh*pblh(i)*unsobklen(i))
613          wm(i) = ustar(i)*phiminv(i)
614          fak2(i) = wm(i)*pblh(i)*vk
615          wstar(i) = (heatv(i)*rg*pblh(i)/zxt)**onet
616          fak3(i) = fakn*wstar(i)/wm(i)
617       ELSE
618          wstar(i) = 0.
619       END IF
620       ! Computes Theta_e for thermal (all cases : to be modified)
621       ! attention ajout therm(i) = virtuelle
622       the_th(i) = th_th(i) + therm(i) + rlvcp*qt_th(i)
623       ! ou:    The_th(i) = Th_th(i) + sqrt(b1)*khfs(i)/wm(i) + RLvCp*qT_th(i)
624    END DO
625
626    ! Main level loop to compute the diffusivities and
627    ! counter-gradient terms:
628
629    DO k = 2, isommet
630
631       ! Find levels within boundary layer:
632
633       DO i = 1, knon
634          unslev(i) = .FALSE.
635          stblev(i) = .FALSE.
636          zm(i) = z(i, k-1)
637          zp(i) = z(i, k)
638          IF (zkmin==0.0 .AND. zp(i)>pblh(i)) zp(i) = pblh(i)
639          IF (zm(i)<pblh(i)) THEN
640             zmzp = 0.5*(zm(i)+zp(i))
641             ! debug
642             ! if (i.EQ.1864) THEN
643             ! PRINT*,'i,pblh(1864),obklen(1864)',i,pblh(i),obklen(i)
644             ! END IF
645
646             zh(i) = zmzp/pblh(i)
647             zl(i) = zmzp*unsobklen(i)
648             zzh(i) = 0.
649             IF (zh(i)<=1.0) zzh(i) = (1.-zh(i))**2
650
651             ! stblev for points zm < plbh and stable and neutral
652             ! unslev for points zm < plbh and unstable
653
654             IF (unstbl(i)) THEN
655                unslev(i) = .TRUE.
656             ELSE
657                stblev(i) = .TRUE.
658             END IF
659          END IF
660       END DO
661       ! PRINT*,'fin calcul niveaux'
662
663       ! Stable and neutral points; set diffusivities; counter-gradient
664       ! terms zero for stable case:
665
666       DO i = 1, knon
667          IF (stblev(i)) THEN
668             IF (zl(i)<=1.) THEN
669                pblk(i) = fak1(i)*zh(i)*zzh(i)/(1.+betas*zl(i))
670             ELSE
671                pblk(i) = fak1(i)*zh(i)*zzh(i)/(betas+zl(i))
672             END IF
673             ! pcfm(i,k) = pblk(i)
674             ! pcfh(i,k) = pcfm(i,k)
675          END IF
676       END DO
677
678       ! unssrf, unstable within surface layer of pbl
679       ! unsout, unstable within outer   layer of pbl
680
681       DO i = 1, knon
682          unssrf(i) = .FALSE.
683          unsout(i) = .FALSE.
684          IF (unslev(i)) THEN
685             IF (zh(i)<sffrac) THEN
686                unssrf(i) = .TRUE.
687             ELSE
688                unsout(i) = .TRUE.
689             END IF
690          END IF
691       END DO
692
693       ! Unstable for surface layer; counter-gradient terms zero
694
695       DO i = 1, knon
696          IF (unssrf(i)) THEN
697             term = (1.-betam*zl(i))**onet
698             pblk(i) = fak1(i)*zh(i)*zzh(i)*term
699             pr(i) = term/sqrt(1.-betah*zl(i))
700          END IF
701       END DO
702       ! PRINT*,'fin counter-gradient terms zero'
703
704       ! Unstable for outer layer; counter-gradient terms non-zero:
705
706       DO i = 1, knon
707          IF (unsout(i)) THEN
708             pblk(i) = fak2(i)*zh(i)*zzh(i)
709             ! cgs(i,k) = fak3(i)/(pblh(i)*wm(i))
710             ! cgh(i,k) = khfs(i)*cgs(i,k)
711             pr(i) = phiminv(i)/phihinv(i) + ccon*fak3(i)/fak
712             ! cgq(i,k) = kqfs(i)*cgs(i,k)
713          END IF
714       END DO
715       ! PRINT*,'fin counter-gradient terms non zero'
716
717       ! For all unstable layers, compute diffusivities and ctrgrad ter m
718
719       ! DO i = 1, knon
720       ! IF (unslev(i)) THEN
721       ! pcfm(i,k) = pblk(i)
722       ! pcfh(i,k) = pblk(i)/pr(i)
723       ! etc cf original
724       ! ENDIF
725       ! ENDDO
726
727       ! For all layers, compute integral info and CTEI
728
729       DO i = 1, knon
730          IF (check(i) .OR. omegafl(i)) THEN
731             IF (.NOT. zsat(i)) THEN
732                ! Th2 = The_th(i) - RLvCp*qT_th(i)
733                th2 = th_th(i)
734                t2 = th2*s(i, k)
735                ! thermodyn functions
736                zdelta = max(0., sign(1.,rtt-t2))
737                qqsat = r2es*foeew(t2, zdelta)/pplay(i, k)
738                qqsat = min(0.5, qqsat)
739                zcor = 1./(1.-retv*qqsat)
740                qqsat = qqsat*zcor
741
742                IF (qqsat<qt_th(i)) THEN
743                   ! on calcule lcl
744                   IF (k==2) THEN
745                      plcl(i) = z(i, k)
746                   ELSE
747                      plcl(i) = z(i, k-1) + (z(i,k-1)-z(i,k))*(qt_th(i)&
748                           -qsatbef(i))/(qsatbef(i)-qqsat)
749                   END IF
750                   zsat(i) = .TRUE.
751                   tbef(i) = t2
752                END IF
753
754                qsatbef(i) = qqsat ! bug dans la version orig ???
755             END IF
756             ! amn ???? cette ligne a deja ete faite normalement ?
757          END IF
758          ! PRINT*,'hbtm2 i,k=',i,k
759       END DO
760    END DO ! end of level loop
761    ! IM 170305 BEG
762    IF (1==0) THEN
763       PRINT *, 'hbtm2  ok'
764    END IF !(1.EQ.0) THEN
765    ! IM 170305 END
766
767  END SUBROUTINE hbtm
768
769END MODULE hbtm_mod
Note: See TracBrowser for help on using the repository browser.