source: LMDZ4/branches/LMDZ4_V2_patch/libf/phylmd/hbtm.F @ 5065

Last change on this file since 5065 was 644, checked in by Laurent Fairhead, 19 years ago

Synchronisation avec tous les diagnostiques de Ionela IM
Inclusion du slab ocean IM
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.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       b212=sqrt(b1*b2)
171       b2sr=sqrt(b2)
172c
173C ============================================================
174C     Fonctions thermo implicites
175C ============================================================
176c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
177c Tetens : pression partielle de vap d'eau e_sat(T)
178c =================================================
179C++ e_sat(T) = r2*exp( r3*(T-Tf)/(T-r4) ) id a r2*FOEWE
180C++ avec :
181C++ Tf = 273.16 K  (Temp de fusion de la glace)
182C++ r2 = 611.14 Pa
183C++ r3 = 17.269 (liquide) 21.875 (solide) adim
184C++ r4 = 35.86             7.66           Kelvin
185C++  q_sat = eps*e_sat/(p-(1-eps)*e_sat)
186C++ derivée :
187C++ =========
188C++                   r3*(Tf-r4)*q_sat(T,p)
189C++ d_qsat_dT = --------------------------------
190C++             (T-r4)^2*( 1-(1-eps)*e_sat(T)/p )
191c++ pour zcvm5=Lv, c'est FOEDE
192c++ Rq :(1.-REPS)*esarg/Parg id a RETV*Qsat
193C     ------------------------------------------------------------------
194c
195c Initialisation
196      RLvCp = RLVTT/RCPD
197      REPS  = RD/RV
198
199c
200c      DO i = 1, klon
201c         pcfh(i,1) = cd_h(i)
202c         pcfm(i,1) = cd_m(i)
203c      ENDDO
204c      DO k = 2, klev
205c      DO i = 1, klon
206c         pcfh(i,k) = zkmin
207c         pcfm(i,k) = zkmin
208c         cgs(i,k) = 0.0
209c         cgh(i,k) = 0.0
210c         cgq(i,k) = 0.0
211c      ENDDO
212c      ENDDO
213c
214c Calculer les hauteurs de chaque couche
215c (geopotentielle Int_dp/ro = Int_[Rd.T.dp/p] z = geop/g)
216c  pourquoi ne pas utiliser Phi/RG ?
217      DO i = 1, knon
218         z(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))
219     .               * (paprs(i,1)-pplay(i,1)) / RG
220         s(i,1) = (pplay(i,1)/paprs(i,1))**RKappa
221      ENDDO
222c                                 s(k) = [pplay(k)/ps]^kappa
223c    + + + + + + + + + pplay  <-> s(k)   t  dp=pplay(k-1)-pplay(k)
224c
225c    -----------------  paprs <-> sig(k)
226c
227c    + + + + + + + + + pplay  <-> s(k-1)
228c
229c
230c    + + + + + + + + + pplay  <-> s(1)   t  dp=paprs-pplay   z(1)
231c
232c    -----------------  paprs <-> sig(1)
233c
234      DO k = 2, klev
235      DO i = 1, knon
236         z(i,k) = z(i,k-1)
237     .              + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)
238     .                   * (pplay(i,k-1)-pplay(i,k)) / RG
239         s(i,k) = (pplay(i,k)/paprs(i,1))**RKappa
240      ENDDO
241      ENDDO
242c  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
243c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
244c +++  Determination des grandeurs de surface  +++++++++++++++++++++
245c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
246c  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
247      DO i = 1, knon
248cAM         IF (thermcep) THEN
249cAM           zdelta=MAX(0.,SIGN(1.,RTT-tsol(i)))
250c           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
251c           zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q(i,1))
252cAM           zxqs= r2es * FOEEW(tsol(i),zdelta)/paprs(i,1)
253cAM           zxqs=MIN(0.5,zxqs)
254cAM           zcor=1./(1.-retv*zxqs)
255cAM           zxqs=zxqs*zcor
256cAM         ELSE
257cAM           IF (tsol(i).LT.t_coup) THEN
258cAM              zxqs = qsats(tsol(i)) / paprs(i,1)
259cAM           ELSE
260cAM              zxqs = qsatl(tsol(i)) / paprs(i,1)
261cAM           ENDIF
262cAM         ENDIF
263c niveau de reference bulk; mais ici, c,a pourrait etre le niveau de ref du thermique
264cAM        zx_alf1 = 1.0
265cAM        zx_alf2 = 1.0 - zx_alf1
266cAM        zxt = (t(i,1)+z(i,1)*RG/RCPD/(1.+RVTMP2*q(i,1)))
267cAM     .        *(1.+RETV*q(i,1))*zx_alf1
268cAM     .      + (t(i,2)+z(i,2)*RG/RCPD/(1.+RVTMP2*q(i,2)))
269cAM     .        *(1.+RETV*q(i,2))*zx_alf2
270cAM        zxu = u(i,1)*zx_alf1+u(i,2)*zx_alf2
271cAM        zxv = v(i,1)*zx_alf1+v(i,2)*zx_alf2
272cAM        zxq = q(i,1)*zx_alf1+q(i,2)*zx_alf2
273cAM     
274cAMAM           zxu = u10m(i)
275cAMAM           zxv = v10m(i)
276cAMAM           zxmod = 1.0+SQRT(zxu**2+zxv**2)
277cAM Niveau de ref choisi a 2m
278        zxt = t2m(i)
279
280c ***************************************************
281c attention, il doit s'agir de <w'theta'>
282c   ;Calcul de tcls virtuel et de w'theta'virtuel
283c   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
284c   tcls=tcls*(1+.608*qcls)
285c
286c   ;Pour avoir w'theta',
287c   ; il faut diviser par ro.Cp
288c   Cp=Cpd*(1+0.84*qcls)
289c   fcs=fcs/(ro_surf*Cp)
290c   ;On transforme w'theta' en w'thetav'
291c   Lv=(2.501-0.00237*(tcls-273.15))*1.E6
292c   xle=xle/(ro_surf*Lv)
293c   fcsv=fcs+.608*xle*tcls
294c ***************************************************
295cAM        khfs(i) = (tsol(i)*(1.+RETV*q(i,1))-zxt) *zxmod*cd_h(i)
296cAM        kqfs(i) = (zxqs-zxq) *zxmod*cd_h(i) * beta(i)
297cAM
298cdif khfs est deja w't'_v / heatv(i) = khfs(i) + RETV*zxt*kqfs(i)
299cAM calcule de Ro = paprs(i,1)/Rd zxt
300cAM convention >0 vers le bas ds lmdz
301        khfs(i) = - flux_t(i,1)*zxt*Rd / (RCPD*paprs(i,1))
302        kqfs(i) = - flux_q(i,1)*zxt*Rd / (paprs(i,1))
303cAM   verifier que khfs et kqfs sont bien de la forme w'l'
304        heatv(i) = khfs(i) + 0.608*zxt*kqfs(i)
305c a comparer aussi aux sorties de clqh : flux_T/RoCp et flux_q/RoLv
306cAM        heatv(i) = khfs(i)
307cAM ustar est en entree
308cAM        taux = zxu *zxmod*cd_m(i)
309cAM        tauy = zxv *zxmod*cd_m(i)
310cAM        ustar(i) = SQRT(taux**2+tauy**2)
311cAM        ustar(i) = MAX(SQRT(ustar(i)),0.01)
312c Theta et qT du thermique sans exces (interpolin vers surf)
313c chgt de niveau du thermique (jeudi 30/12/1999)
314c (interpolation lineaire avant integration phi_h)
315cAM        qT_th(i) = zxqs*beta(i) + 4./z(i,1)*(q(i,1)-zxqs*beta(i))
316cAM        qT_th(i) = max(qT_th(i),q(i,1))
317        qT_th(i) = q2m(i)
318cn The_th restera la Theta du thermique sans exces jusqu'a 2eme calcul
319cn reste a regler convention P) pour Theta
320c        The_th(i) = tsol(i) + 4./z(i,1)*(t(i,1)-tsol(i))
321c     -                      + RLvCp*qT_th(i)
322cAM        Th_th(i) = tsol(i) + 4./z(i,1)*(t(i,1)-tsol(i))
323        Th_th(i) = t2m(i)
324      ENDDO
325c
326      DO i = 1, knon
327         rhino(i,1) = 0.0   ! Global Richardson
328         check(i) = .TRUE.
329         pblh(i) = z(i,1)   ! on initialise pblh a l'altitude du 1er niveau
330         plcl(i) = 6000.
331c Lambda = -u*^3 / (alpha.g.kvon.<w'Theta'v>
332         obklen(i) = -t(i,1)*ustar(i)**3/(RG*vk*heatv(i))
333         trmb1(i)   = 0.
334         trmb2(i)   = 0.
335         trmb3(i) = 0.
336      ENDDO
337
338C
339c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
340C PBL height calculation:
341C Search for level of pbl. Scan upward until the Richardson number between
342C the first level and the current level exceeds the "critical" value.
343C (bonne idee Nu de separer le Ric et l'exces de temp du thermique)
344c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
345      fac = 100.0
346      DO k = 2, isommet
347      DO i = 1, knon
348      IF (check(i)) THEN
349! pourquoi / niveau 1 (au lieu du sol) et le terme en u*^2 ?
350ctest     zdu2 = (u(i,k)-u(i,1))**2+(v(i,k)-v(i,1))**2+fac*ustar(i)**2
351         zdu2 = u(i,k)**2+v(i,k)**2
352         zdu2 = max(zdu2,1.0e-20)
353c Theta_v environnement
354         zthvd=t(i,k)/s(i,k)*(1.+RETV*q(i,k))
355c
356c therm Theta_v sans exces (avec hypothese fausse de H&B, sinon,
357c passer par Theta_e et virpot)
358c         zthvu=t(i,1)/s(i,1)*(1.+RETV*q(i,1))
359cAM         zthvu = Th_th(i)*(1.+RETV*q(i,1))
360         zthvu = Th_th(i)*(1.+RETV*qT_th(i))
361c  Le Ri par Theta_v
362cAM         rhino(i,k) = (z(i,k)-z(i,1))*RG*(zthvd-zthvu)
363cAM     .               /(zdu2*0.5*(zthvd+zthvu))
364cAM On a nveau de ref a 2m ???
365         rhino(i,k) = (z(i,k)-zref)*RG*(zthvd-zthvu)
366     .               /(zdu2*0.5*(zthvd+zthvu))
367c
368         IF (rhino(i,k).GE.ricr) THEN
369           pblh(i) = z(i,k-1) + (z(i,k-1)-z(i,k)) *
370     .              (ricr-rhino(i,k-1))/(rhino(i,k-1)-rhino(i,k))
371c test04
372           pblh(i) = pblh(i) + 100.
373           pblT(i) = t(i,k-1) + (t(i,k)-t(i,k-1)) *
374     .              (pblh(i)-z(i,k-1))/(z(i,k)-z(i,k-1))
375           check(i) = .FALSE.
376         ENDIF
377      ENDIF
378      ENDDO
379      ENDDO
380
381C
382C Set pbl height to maximum value where computation exceeds number of
383C layers allowed
384C
385      DO i = 1, knon
386        if (check(i)) pblh(i) = z(i,isommet)
387      ENDDO
388C
389C Improve estimate of pbl height for the unstable points.
390C Find unstable points (sensible heat flux is upward):
391C
392      DO i = 1, knon
393      IF (heatv(i) .GT. 0.) THEN
394        unstbl(i) = .TRUE.
395        check(i) = .TRUE.
396      ELSE
397        unstbl(i) = .FALSE.
398        check(i) = .FALSE.
399      ENDIF
400      ENDDO
401C
402C For the unstable case, compute velocity scale and the
403C convective temperature excess:
404C
405      DO i = 1, knon
406      IF (check(i)) THEN
407        phiminv(i) = (1.-binm*pblh(i)/obklen(i))**onet
408c ***************************************************
409c Wm ? et W* ? c'est la formule pour z/h < .1
410c   ;Calcul de w* ;;
411c   ;;;;;;;;;;;;;;;;
412c   w_star=((g/tcls)*fcsv*z(ind))^(1/3.) [ou prendre la premiere approx de h)
413c   ;; CALCUL DE wm ;;
414c   ;;;;;;;;;;;;;;;;;;
415c   ; Ici on considerera que l'on est dans la couche de surf jusqu'a 100m
416c   ; On prend svt couche de surface=0.1*h mais on ne connait pas h
417c   ;;;;;;;;;;;Dans la couche de surface
418c   if (z(ind) le 20) then begin
419c   Phim=(1.-15.*(z(ind)/L))^(-1/3.)
420c   wm=u_star/Phim
421c   ;;;;;;;;;;;En dehors de la couche de surface
422c   endif else if (z(ind) gt 20) then begin
423c   wm=(u_star^3+c1*w_star^3)^(1/3.)
424c   endif
425c ***************************************************
426        wm(i)= ustar(i)*phiminv(i)
427c======================================================================
428cvaleurs de Dominique Lambert de la campagne SEMAPHORE :
429c <T'^2> = 100.T*^2; <q'^2> = 20.q*^2 a 10m
430c <Tv'^2> = (1+1.2q).100.T* + 1.2Tv.sqrt(20*100).T*.q* + (.608*Tv)^2*20.q*^2;
431c et dTetavS = sqrt(<Tv'^2>) ainsi calculee.
432c avec : T*=<w'T'>_s/w* et q*=<w'q'>/w*
433c !!! on peut donc utiliser w* pour les fluctuations <-> Lambert
434c(leur corellation pourrait dependre de beta par ex)
435c  if fcsv(i,j) gt 0 then begin
436c    dTetavs=b1*(1.+2.*.608*q_10(i,j))*(fcs(i,j)/wm(i,j))^2+$
437c    (.608*Thetav_10(i,j))^2*b2*(xle(i,j)/wm(i,j))^2+$
438c    2.*.608*thetav_10(i,j)*sqrt(b1*b2)*(xle(i,j)/wm(i,j))*(fcs(i,j)/wm(i,j))
439c    dqs=b2*(xle(i,j)/wm(i,j))^2
440c    theta_s(i,j)=thetav_10(i,j)+sqrt(dTetavs)
441c    q_s(i,j)=q_10(i,j)+sqrt(dqs)
442c  endif else begin
443c    Theta_s(i,j)=thetav_10(i,j)
444c    q_s(i,j)=q_10(i,j)
445c  endelse
446c======================================================================
447c
448cHBTM        therm(i) = heatv(i)*fak/wm(i)
449c forme Mathieu :
450        q_star = kqfs(i)/wm(i)
451        t_star = khfs(i)/wm(i)
452cIM 091204 BEG
453        IF(1.EQ.0) THEN
454        IF(t_star.LT.0..OR.q_star.LT.0.) THEN
455          print*,'i t_star q_star khfs kqfs wm',i,t_star,q_star,
456     $    khfs(i),kqfs(i),wm(i)
457        ENDIF
458        ENDIF
459cIM 091204 END
460cAM Nveau cde ref 2m =>
461cAM        therm(i) = sqrt( b1*(1.+2.*RETV*q(i,1))*t_star**2
462cAM     +             + (RETV*T(i,1))**2*b2*q_star**2
463cAM     +             + 2.*RETV*T(i,1)*b212*q_star*t_star
464cAM     +                 )
465cIM 091204 BEG
466        a1=b1*(1.+2.*RETV*qT_th(i))*t_star**2
467        a2=(RETV*Th_th(i))**2*b2*q_star**2
468        a3=2.*RETV*Th_th(i)*b212*q_star*t_star
469        aa=a1+a2+a3
470        IF(1.EQ.0) THEN
471        IF (aa.LT.0.) THEN
472         print*,'i a1 a2 a3 aa',i,a1,a2,a3,aa
473         print*,'i qT_th Th_th t_star q_star RETV b1 b2 b212',
474     $   i,qT_th(i),Th_th(i),t_star,q_star,RETV,b1,b2,b212
475        ENDIF
476        ENDIF
477cIM 091204 END
478        therm(i) = sqrt( b1*(1.+2.*RETV*qT_th(i))*t_star**2
479     +             + (RETV*Th_th(i))**2*b2*q_star**2
480cIM 101204  +             + 2.*RETV*Th_th(i)*b212*q_star*t_star
481     +             + max(0.,2.*RETV*Th_th(i)*b212*q_star*t_star)
482     +                 )
483c
484c Theta et qT du thermique (forme H&B) avec exces
485c (attention, on ajoute therm(i) qui est virtuelle ...)
486c pourquoi pas sqrt(b1)*t_star ?
487c        dqs = b2sr*kqfs(i)/wm(i)
488        qT_th(i) = qT_th(i)  + b2sr*q_star
489cnew on differre le calcul de Theta_e
490c        The_th(i) = The_th(i) + therm(i) + RLvCp*qT_th(i)
491c ou:    The_th(i) = The_th(i) + sqrt(b1)*khfs(i)/wm(i) + RLvCp*qT_th(i)
492        rhino(i,1) = 0.0
493      ENDIF
494      ENDDO
495C
496c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
497c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
498C ++ Improve pblh estimate for unstable conditions using the +++++++
499C ++          convective temperature excess :                +++++++
500c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
501c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
502C
503      DO k = 2, isommet
504      DO i = 1, knon
505      IF (check(i)) THEN
506ctest     zdu2 = (u(i,k)-u(i,1))**2+(v(i,k)-v(i,1))**2+fac*ustar(i)**2
507         zdu2 = u(i,k)**2+v(i,k)**2
508         zdu2 = max(zdu2,1.0e-20)
509c Theta_v environnement
510         zthvd=t(i,k)/s(i,k)*(1.+RETV*q(i,k))
511c
512c et therm Theta_v (avec hypothese de constance de H&B,
513c         zthvu=(t(i,1)+therm(i))/s(i,1)*(1.+RETV*q(i,1))
514         zthvu = Th_th(i)*(1.+RETV*qT_th(i)) + therm(i)
515
516c
517c  Le Ri par Theta_v
518cAM Niveau de ref 2m
519cAM         rhino(i,k) = (z(i,k)-z(i,1))*RG*(zthvd-zthvu)
520cAM     .               /(zdu2*0.5*(zthvd+zthvu))
521         rhino(i,k) = (z(i,k)-zref)*RG*(zthvd-zthvu)
522     .               /(zdu2*0.5*(zthvd+zthvu))
523c
524c
525         IF (rhino(i,k).GE.ricr) THEN
526           pblh(i) = z(i,k-1) + (z(i,k-1)-z(i,k)) *
527     .              (ricr-rhino(i,k-1))/(rhino(i,k-1)-rhino(i,k))
528c test04
529           pblh(i) = pblh(i) + 100.
530           pblT(i) = t(i,k-1) + (t(i,k)-t(i,k-1)) *
531     .              (pblh(i)-z(i,k-1))/(z(i,k)-z(i,k-1))
532           check(i) = .FALSE.
533cIM 170305 BEG
534      IF(1.EQ.0) THEN
535c debug print -120;34       -34-        58 et    0;26 wamp
536      if (i.eq.950.or.i.eq.192.or.i.eq.624.or.i.eq.118) then
537            print*,' i,Th_th,Therm,qT :',i,Th_th(i),therm(i),qT_th(i)
538            q_star = kqfs(i)/wm(i)
539            t_star = khfs(i)/wm(i)
540            print*,'q* t*, b1,b2,b212 ',q_star,t_star
541     -            , b1*(1.+2.*RETV*qT_th(i))*t_star**2
542     -            , (RETV*Th_th(i))**2*b2*q_star**2
543     -            , 2.*RETV*Th_th(i)*b212*q_star*t_star
544            print*,'zdu2 ,100.*ustar(i)**2',zdu2 ,fac*ustar(i)**2
545      endif
546      ENDIF !(1.EQ.0) THEN
547cIM 170305 END
548c             q_star = kqfs(i)/wm(i)
549c             t_star = khfs(i)/wm(i)
550c             trmb1(i) = b1*(1.+2.*RETV*q(i,1))*t_star**2
551c             trmb2(i) = (RETV*T(i,1))**2*b2*q_star**2
552c Omega now   trmb3(i) = 2.*RETV*T(i,1)*b212*q_star*t_star
553         ENDIF
554      ENDIF
555      ENDDO
556      ENDDO
557C
558C Set pbl height to maximum value where computation exceeds number of
559C layers allowed
560C
561      DO i = 1, knon
562        if (check(i)) pblh(i) = z(i,isommet)
563      ENDDO
564C
565C PBL height must be greater than some minimum mechanical mixing depth
566C Several investigators have proposed minimum mechanical mixing depth
567C relationships as a function of the local friction velocity, u*.  We
568C make use of a linear relationship of the form h = c u* where c=700.
569C The scaling arguments that give rise to this relationship most often
570C represent the coefficient c as some constant over the local coriolis
571C parameter.  Here we make use of the experimental results of Koracin
572C and Berkowicz (1988) [BLM, Vol 43] for wich they recommend 0.07/f
573C where f was evaluated at 39.5 N and 52 N.  Thus we use a typical mid
574C latitude value for f so that c = 0.07/f = 700.
575C
576      DO i = 1, knon
577        pblmin  = 700.0*ustar(i)
578        pblh(i) = MAX(pblh(i),pblmin)
579c par exemple :
580        pblT(i) = t(i,2) + (t(i,3)-t(i,2)) *
581     .              (pblh(i)-z(i,2))/(z(i,3)-z(i,2))
582      ENDDO
583
584C ********************************************************************
585C  pblh is now available; do preparation for diffusivity calculation :
586C ********************************************************************
587      DO i = 1, knon
588        check(i) = .TRUE.
589        Zsat(i)   = .FALSE.
590c omegafl utilise pour prolongement CAPE
591        omegafl(i) = .FALSE.
592        Cape(i)   = 0.
593        Kape(i)   = 0.
594        EauLiq(i) = 0.
595        CTEI(i)   = 0.
596        pblk(i) = 0.0
597        fak1(i) = ustar(i)*pblh(i)*vk
598C
599C Do additional preparation for unstable cases only, set temperature
600C and moisture perturbations depending on stability.
601C *** Rq: les formule sont prises dans leur forme CS ***
602        IF (unstbl(i)) THEN
603cAM Niveau de ref du thermique
604cAM          zxt=(t(i,1)-z(i,1)*0.5*RG/RCPD/(1.+RVTMP2*q(i,1)))
605cAM     .         *(1.+RETV*q(i,1))
606          zxt=(Th_th(i)-zref*0.5*RG/RCPD/(1.+RVTMP2*qT_th(i)))
607     .         *(1.+RETV*qT_th(i))
608          phiminv(i) = (1. - binm*pblh(i)/obklen(i))**onet
609          phihinv(i) = sqrt(1. - binh*pblh(i)/obklen(i))
610          wm(i)      = ustar(i)*phiminv(i)
611          fak2(i)    = wm(i)*pblh(i)*vk
612          wstr(i)    = (heatv(i)*RG*pblh(i)/zxt)**onet
613          fak3(i)    = fakn*wstr(i)/wm(i)
614        ENDIF
615c Computes Theta_e for thermal (all cases : to be modified)
616c   attention ajout therm(i) = virtuelle
617        The_th(i) = Th_th(i) + therm(i) + RLvCp*qT_th(i)
618c ou:    The_th(i) = Th_th(i) + sqrt(b1)*khfs(i)/wm(i) + RLvCp*qT_th(i)
619      ENDDO
620
621C Main level loop to compute the diffusivities and
622C counter-gradient terms:
623C
624      DO 1000 k = 2, isommet
625C
626C Find levels within boundary layer:
627C
628        DO i = 1, knon
629          unslev(i) = .FALSE.
630          stblev(i) = .FALSE.
631          zm(i) = z(i,k-1)
632          zp(i) = z(i,k)
633          IF (zkmin.EQ.0.0 .AND. zp(i).GT.pblh(i)) zp(i) = pblh(i)
634          IF (zm(i) .LT. pblh(i)) THEN
635            zmzp = 0.5*(zm(i) + zp(i))
636C debug
637c          if (i.EQ.1864) then
638c             print*,'i,pblh(1864),obklen(1864)',i,pblh(i),obklen(i)
639c          endif
640
641            zh(i) = zmzp/pblh(i)
642            zl(i) = zmzp/obklen(i)
643            zzh(i) = 0.
644            IF (zh(i).LE.1.0) zzh(i) = (1. - zh(i))**2
645C
646C stblev for points zm < plbh and stable and neutral
647C unslev for points zm < plbh and unstable
648C
649            IF (unstbl(i)) THEN
650              unslev(i) = .TRUE.
651            ELSE
652              stblev(i) = .TRUE.
653            ENDIF
654          ENDIF
655        ENDDO
656c        print*,'fin calcul niveaux'
657C
658C Stable and neutral points; set diffusivities; counter-gradient
659C terms zero for stable case:
660C
661        DO i = 1, knon
662          IF (stblev(i)) THEN
663            IF (zl(i).LE.1.) THEN
664              pblk(i) = fak1(i)*zh(i)*zzh(i)/(1. + betas*zl(i))
665            ELSE
666              pblk(i) = fak1(i)*zh(i)*zzh(i)/(betas + zl(i))
667            ENDIF
668c            pcfm(i,k) = pblk(i)
669c            pcfh(i,k) = pcfm(i,k)
670          ENDIF
671        ENDDO
672C
673C unssrf, unstable within surface layer of pbl
674C unsout, unstable within outer   layer of pbl
675C
676        DO i = 1, knon
677          unssrf(i) = .FALSE.
678          unsout(i) = .FALSE.
679          IF (unslev(i)) THEN
680            IF (zh(i).lt.sffrac) THEN
681              unssrf(i) = .TRUE.
682            ELSE
683              unsout(i) = .TRUE.
684            ENDIF
685          ENDIF
686        ENDDO
687C
688C Unstable for surface layer; counter-gradient terms zero
689C
690        DO i = 1, knon
691          IF (unssrf(i)) THEN
692            term = (1. - betam*zl(i))**onet
693            pblk(i) = fak1(i)*zh(i)*zzh(i)*term
694            pr(i) = term/sqrt(1. - betah*zl(i))
695          ENDIF
696        ENDDO
697c        print*,'fin counter-gradient terms zero'
698C
699C Unstable for outer layer; counter-gradient terms non-zero:
700C
701        DO i = 1, knon
702          IF (unsout(i)) THEN
703            pblk(i) = fak2(i)*zh(i)*zzh(i)
704c            cgs(i,k) = fak3(i)/(pblh(i)*wm(i))
705c            cgh(i,k) = khfs(i)*cgs(i,k)
706            pr(i) = phiminv(i)/phihinv(i) + ccon*fak3(i)/fak
707c            cgq(i,k) = kqfs(i)*cgs(i,k)
708          ENDIF
709        ENDDO
710c        print*,'fin counter-gradient terms non zero'
711C
712C For all unstable layers, compute diffusivities and ctrgrad ter m
713C
714c        DO i = 1, knon
715c        IF (unslev(i)) THEN
716c            pcfm(i,k) = pblk(i)
717c            pcfh(i,k) = pblk(i)/pr(i)
718c etc cf original
719c        ENDIF
720c        ENDDO
721C
722C For all layers, compute integral info and CTEI
723C
724        DO i = 1, knon
725        if (check(i).or.omegafl(i)) then
726          if (.not.Zsat(i)) then
727c            Th2 = The_th(i) - RLvCp*qT_th(i)
728            Th2 = Th_th(i)
729            T2 = Th2*s(i,k)
730c thermodyn functions
731            zdelta=MAX(0.,SIGN(1.,RTT-T2))
732            qqsat= r2es * FOEEW(T2,zdelta)/pplay(i,k)
733            qqsat=MIN(0.5,qqsat)
734            zcor=1./(1.-retv*qqsat)
735            qqsat=qqsat*zcor
736c
737            if (qqsat.lt.qT_th(i)) then
738c on calcule lcl
739              if (k.eq.2) then
740                plcl(i) = z(i,k)
741              else
742                plcl(i) =  z(i,k-1) + (z(i,k-1)-z(i,k)) *
743     .                 (qT_th(i)-qsatbef(i))/(qsatbef(i)-qqsat)
744              endif
745              Zsat(i) = .true.
746              Tbef(i) = T2
747            endif
748c
749          endif
750          qsatbef(i) = qqsat
751camn ???? cette ligne a deja ete faite normalement ?
752        endif
753c            print*,'hbtm2 i,k=',i,k
754        ENDDO
755 1000 continue           ! end of level loop
756cIM 170305 BEG
757        IF(1.EQ.0) THEN
758            print*,'hbtm2  ok'
759        ENDIF !(1.EQ.0) THEN
760cIM 170305 END
761      RETURN
762      END
Note: See TracBrowser for help on using the repository browser.