source: LMDZ4/branches/LMDZ4_AR5/libf/phylmd/hbtm.F @ 1500

Last change on this file since 1500 was 1279, checked in by Laurent Fairhead, 14 years ago

Merged LMDZ4-dev branch changes r1241:1278 into the trunk
Running trunk and LMDZ4-dev in LMDZOR configuration on local
machine (sequential) and SX8 (4-proc) yields identical results
(restart and restartphy are identical binarily)
Log history from r1241 to r1278 is available by switching to
source:LMDZ4/branches/LMDZ4-dev-20091210
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 27.3 KB
Line 
1!
2! $Header$
3!
4
5      SUBROUTINE HBTM(knon, paprs, pplay,
6     .                t2m,t10m,q2m,q10m,ustar,
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 paprs(klon,klev+1) ! pression a inter-couche (Pa)
57      REAL pplay(klon,klev)   ! pression au milieu de couche (Pa)
58      REAL flux_t(klon,klev), flux_q(klon,klev)     ! Flux
59      REAL u(klon,klev) ! vitesse U (m/s)
60      REAL v(klon,klev) ! vitesse V (m/s)
61      REAL t(klon,klev) ! temperature (K)
62      REAL q(klon,klev) ! vapeur d'eau (kg/kg)
63cAM      REAL cd_h(klon) ! coefficient de friction au sol pour chaleur
64cAM      REAL cd_m(klon) ! coefficient de friction au sol pour vitesse
65c
66      INTEGER isommet
67cum      PARAMETER (isommet=klev) ! limite max sommet pbl
68      REAL vk
69      PARAMETER (vk=0.35)     ! Von Karman => passer a .41 ! cf U.Olgstrom
70      REAL ricr
71      PARAMETER (ricr=0.4)
72      REAL fak
73      PARAMETER (fak=8.5)     ! b calcul du Prandtl et de dTetas
74      REAL fakn
75      PARAMETER (fakn=7.2)    ! a
76      REAL onet
77      PARAMETER (onet=1.0/3.0)
78      REAL t_coup
79      PARAMETER(t_coup=273.15)
80      REAL zkmin
81      PARAMETER (zkmin=0.01)
82      REAL betam
83      PARAMETER (betam=15.0)  ! pour Phim / h dans la S.L stable
84      REAL betah
85      PARAMETER (betah=15.0)
86      REAL betas
87      PARAMETER (betas=5.0)   ! Phit dans la S.L. stable (mais 2 formes / z/OBL<>1
88      REAL sffrac
89      PARAMETER (sffrac=0.1)  ! S.L. = z/h < .1
90      REAL binm
91      PARAMETER (binm=betam*sffrac)
92      REAL binh
93      PARAMETER (binh=betah*sffrac)
94      REAL ccon
95      PARAMETER (ccon=fak*sffrac*vk)
96c
97      REAL q_star,t_star
98      REAL b1,b2,b212,b2sr     ! Lambert correlations T' q' avec T* q*
99      PARAMETER (b1=70.,b2=20.)
100c
101      REAL z(klon,klev)
102cAM      REAL pcfm(klon,klev), pcfh(klon,klev)
103cAM
104      REAL zref
105      PARAMETER (zref=2.)    ! Niveau de ref a 2m peut eventuellement
106c                              etre choisi a 10m
107cMA
108c
109      INTEGER i, k, j
110      REAL zxt
111cAM      REAL zxt, zxq, zxu, zxv, zxmod, taux, tauy
112cAM      REAL zx_alf1, zx_alf2 ! parametres pour extrapolation
113      REAL khfs(klon)       ! surface kinematic heat flux [mK/s]
114      REAL kqfs(klon)       ! sfc kinematic constituent flux [m/s]
115      REAL heatv(klon)      ! surface virtual heat flux
116      REAL rhino(klon,klev) ! bulk Richardon no. mais en Theta_v
117      LOGICAL unstbl(klon)  ! pts w/unstbl pbl (positive virtual ht flx)
118      LOGICAL stblev(klon)  ! stable pbl with levels within pbl
119      LOGICAL unslev(klon)  ! unstbl pbl with levels within pbl
120      LOGICAL unssrf(klon)  ! unstb pbl w/lvls within srf pbl lyr
121      LOGICAL unsout(klon)  ! unstb pbl w/lvls in outer pbl lyr
122      LOGICAL check(klon)   ! True=>chk if Richardson no.>critcal
123      LOGICAL omegafl(klon) ! flag de prolongerment cape pour pt Omega
124      REAL pblh(klon)
125      REAL pblT(klon)
126      REAL plcl(klon)
127cAM      REAL cgh(klon,2:klev) ! counter-gradient term for heat [K/m]
128cAM      REAL cgq(klon,2:klev) ! counter-gradient term for constituents
129cAM      REAL cgs(klon,2:klev) ! counter-gradient star (cg/flux)
130      REAL obklen(klon)     ! Monin-Obukhov lengh
131cAM      REAL ztvd, ztvu,
132      REAL zdu2
133      REAL therm(klon)      ! thermal virtual temperature excess
134      REAL trmb1(klon),trmb2(klon),trmb3(klon)
135C  Algorithme thermique
136      REAL s(klon,klev)     ! [P/Po]^Kappa milieux couches
137      REAL Th_th(klon)      ! potential temperature of thermal
138      REAL The_th(klon)     ! equivalent potential temperature of thermal
139      REAL qT_th(klon)      ! total water  of thermal
140      REAL Tbef(klon)       ! T thermique niveau precedent
141      REAL qsatbef(klon)
142      LOGICAL Zsat(klon)    ! le thermique est sature
143      REAL Cape(klon)       ! Cape du thermique
144      REAL Kape(klon)       ! Cape locale
145      REAL EauLiq(klon)     ! Eau liqu integr du thermique
146      REAL ctei(klon)       ! Critere d'instab d'entrainmt des nuages de CL
147      REAL the1,the2,aa,bb,zthvd,zthvu,xintpos,qqsat
148cIM 091204 BEG
149      REAL a1,a2,a3
150cIM 091204 END
151      REAL xhis,rnum,denom,th1,th2,thv1,thv2,ql2
152      REAL dqsat_dt,qsat2,qT1,q2,t1,t2,xnull,delt_the
153      REAL delt_qt,delt_2,quadsat,spblh,reduc
154c
155      REAL phiminv(klon)    ! inverse phi function for momentum
156      REAL phihinv(klon)    ! inverse phi function for heat
157      REAL wm(klon)         ! turbulent velocity scale for momentum
158      REAL fak1(klon)       ! k*ustar*pblh
159      REAL fak2(klon)       ! k*wm*pblh
160      REAL fak3(klon)       ! fakn*wstr/wm
161      REAL pblk(klon)       ! level eddy diffusivity for momentum
162      REAL pr(klon)         ! Prandtl number for eddy diffusivities
163      REAL zl(klon)         ! zmzp / Obukhov length
164      REAL zh(klon)         ! zmzp / pblh
165      REAL zzh(klon)        ! (1-(zmzp/pblh))**2
166      REAL wstr(klon)       ! w*, convective velocity scale
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          wstr(i)    = (heatv(i)*RG*pblh(i)/zxt)**onet
628          fak3(i)    = fakn*wstr(i)/wm(i)
629        ENDIF
630c Computes Theta_e for thermal (all cases : to be modified)
631c   attention ajout therm(i) = virtuelle
632        The_th(i) = Th_th(i) + therm(i) + RLvCp*qT_th(i)
633c ou:    The_th(i) = Th_th(i) + sqrt(b1)*khfs(i)/wm(i) + RLvCp*qT_th(i)
634      ENDDO
635
636C Main level loop to compute the diffusivities and
637C counter-gradient terms:
638C
639      DO 1000 k = 2, isommet
640C
641C Find levels within boundary layer:
642C
643        DO i = 1, knon
644          unslev(i) = .FALSE.
645          stblev(i) = .FALSE.
646          zm(i) = z(i,k-1)
647          zp(i) = z(i,k)
648          IF (zkmin.EQ.0.0 .AND. zp(i).GT.pblh(i)) zp(i) = pblh(i)
649          IF (zm(i) .LT. pblh(i)) THEN
650            zmzp = 0.5*(zm(i) + zp(i))
651C debug
652c          if (i.EQ.1864) then
653c             print*,'i,pblh(1864),obklen(1864)',i,pblh(i),obklen(i)
654c          endif
655
656            zh(i) = zmzp/pblh(i)
657            zl(i) = zmzp/obklen(i)
658            zzh(i) = 0.
659            IF (zh(i).LE.1.0) zzh(i) = (1. - zh(i))**2
660C
661C stblev for points zm < plbh and stable and neutral
662C unslev for points zm < plbh and unstable
663C
664            IF (unstbl(i)) THEN
665              unslev(i) = .TRUE.
666            ELSE
667              stblev(i) = .TRUE.
668            ENDIF
669          ENDIF
670        ENDDO
671c        print*,'fin calcul niveaux'
672C
673C Stable and neutral points; set diffusivities; counter-gradient
674C terms zero for stable case:
675C
676        DO i = 1, knon
677          IF (stblev(i)) THEN
678            IF (zl(i).LE.1.) THEN
679              pblk(i) = fak1(i)*zh(i)*zzh(i)/(1. + betas*zl(i))
680            ELSE
681              pblk(i) = fak1(i)*zh(i)*zzh(i)/(betas + zl(i))
682            ENDIF
683c            pcfm(i,k) = pblk(i)
684c            pcfh(i,k) = pcfm(i,k)
685          ENDIF
686        ENDDO
687C
688C unssrf, unstable within surface layer of pbl
689C unsout, unstable within outer   layer of pbl
690C
691        DO i = 1, knon
692          unssrf(i) = .FALSE.
693          unsout(i) = .FALSE.
694          IF (unslev(i)) THEN
695            IF (zh(i).lt.sffrac) THEN
696              unssrf(i) = .TRUE.
697            ELSE
698              unsout(i) = .TRUE.
699            ENDIF
700          ENDIF
701        ENDDO
702C
703C Unstable for surface layer; counter-gradient terms zero
704C
705        DO i = 1, knon
706          IF (unssrf(i)) THEN
707            term = (1. - betam*zl(i))**onet
708            pblk(i) = fak1(i)*zh(i)*zzh(i)*term
709            pr(i) = term/sqrt(1. - betah*zl(i))
710          ENDIF
711        ENDDO
712c        print*,'fin counter-gradient terms zero'
713C
714C Unstable for outer layer; counter-gradient terms non-zero:
715C
716        DO i = 1, knon
717          IF (unsout(i)) THEN
718            pblk(i) = fak2(i)*zh(i)*zzh(i)
719c            cgs(i,k) = fak3(i)/(pblh(i)*wm(i))
720c            cgh(i,k) = khfs(i)*cgs(i,k)
721            pr(i) = phiminv(i)/phihinv(i) + ccon*fak3(i)/fak
722c            cgq(i,k) = kqfs(i)*cgs(i,k)
723          ENDIF
724        ENDDO
725c        print*,'fin counter-gradient terms non zero'
726C
727C For all unstable layers, compute diffusivities and ctrgrad ter m
728C
729c        DO i = 1, knon
730c        IF (unslev(i)) THEN
731c            pcfm(i,k) = pblk(i)
732c            pcfh(i,k) = pblk(i)/pr(i)
733c etc cf original
734c        ENDIF
735c        ENDDO
736C
737C For all layers, compute integral info and CTEI
738C
739        DO i = 1, knon
740        if (check(i).or.omegafl(i)) then
741          if (.not.Zsat(i)) then
742c            Th2 = The_th(i) - RLvCp*qT_th(i)
743            Th2 = Th_th(i)
744            T2 = Th2*s(i,k)
745c thermodyn functions
746            zdelta=MAX(0.,SIGN(1.,RTT-T2))
747            qqsat= r2es * FOEEW(T2,zdelta)/pplay(i,k)
748            qqsat=MIN(0.5,qqsat)
749            zcor=1./(1.-retv*qqsat)
750            qqsat=qqsat*zcor
751c
752            if (qqsat.lt.qT_th(i)) then
753c on calcule lcl
754              if (k.eq.2) then
755                plcl(i) = z(i,k)
756              else
757                plcl(i) =  z(i,k-1) + (z(i,k-1)-z(i,k)) *
758     .                 (qT_th(i)-qsatbef(i))/(qsatbef(i)-qqsat)
759              endif
760              Zsat(i) = .true.
761              Tbef(i) = T2
762            endif
763c
764            qsatbef(i) = qqsat    ! bug dans la version orig ???
765          endif
766camn ???? cette ligne a deja ete faite normalement ?
767        endif
768c            print*,'hbtm2 i,k=',i,k
769        ENDDO
770 1000 continue           ! end of level loop
771cIM 170305 BEG
772        IF(1.EQ.0) THEN
773            print*,'hbtm2  ok'
774        ENDIF !(1.EQ.0) THEN
775cIM 170305 END
776      RETURN
777      END
Note: See TracBrowser for help on using the repository browser.