source: LMDZ4/tags/V3_conv/libf/phylmd/hbtm.F @ 1386

Last change on this file since 1386 was 704, checked in by Laurent Fairhead, 18 years ago

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