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

Last change on this file since 5218 was 5153, checked in by abarral, 5 months ago

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