source: LMDZ5/branches/LF-private/libf/phylmd/hbtm.F @ 4990

Last change on this file since 4990 was 1818, checked in by idelkadi, 11 years ago

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