source: LMDZ4/branches/pre_V3/libf/phylmd/hbtm.F @ 5440

Last change on this file since 5440 was 674, checked in by lmdzadmin, 19 years ago

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