source: lmdz_wrf/trunk/WRFV3/lmdz/hbtm.F90 @ 1531

Last change on this file since 1531 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 28.6 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
13!c***************************************************************
14!c*                                                             *
15!c* HBTM2   D'apres Holstag&Boville et Troen&Mahrt              *
16!c*                 JAS 47              BLM                     *
17!c* Algorithme These Anne Mathieu                               *
18!c* Critere d'Entrainement Peter Duynkerke (JAS 50)             *
19!c* written by  : Anne MATHIEU & Alain LAHELLEC, 22/11/99       *
20!c* features : implem. exces Mathieu                            *
21!c***************************************************************
22!c* mods : decembre 99 passage th a niveau plus bas. voir fixer *
23!c* la prise du th a z/Lambda = -.2 (max Ray)                   *
24!c* Autre algo : entrainement ~ Theta+v =cste mais comment=>The?*
25!c* on peut fixer q a .7qsat(cf non adiab)=>T2 et The2          *
26!c* voir aussi //KE pblh = niveau The_e ou l = env.             *
27!c***************************************************************
28!c* fin therm a la HBTM passage a forme Mathieu 12/09/2001      *
29!c***************************************************************
30!c*
31!c
32!c
33!cAM Fev 2003
34!c Adaptation a LMDZ version couplee
35!c
36!c Pour le moment on fait passer en argument les grdeurs de surface :
37!c flux, t,q2m, t,q10m, on va utiliser systematiquement les grdeurs a 2m ms
38!c on garde la possibilite de changer si besoin est (jusqu'a present la
39!c forme de HB avec le 1er niveau modele etait conservee)
40!c
41!c
42!c
43!c
44!c
45!cym#include "dimensions.h"
46!cym#include "dimphy.h"
47#include "YOMCST.h"
48      REAL RLvCp, REPS
49!c Arguments:
50!c
51      INTEGER knon ! nombre de points a calculer
52!cAM
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)
63!cAM      REAL cd_h(klon) ! coefficient de friction au sol pour chaleur
64!cAM      REAL cd_m(klon) ! coefficient de friction au sol pour vitesse
65!c
66      INTEGER isommet
67!cum      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)
96!c
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.)
100!c
101      REAL z(klon,klev)
102!cAM      REAL pcfm(klon,klev), pcfh(klon,klev)
103!cAM
104      REAL zref
105      PARAMETER (zref=2.)    ! Niveau de ref a 2m peut eventuellement
106!c                              etre choisi a 10m
107!cMA
108!c
109      INTEGER i, k, j
110      REAL zxt
111!cAM      REAL zxt, zxq, zxu, zxv, zxmod, taux, tauy
112!cAM      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)
127!cAM      REAL cgh(klon,2:klev) ! counter-gradient term for heat [K/m]
128!cAM      REAL cgq(klon,2:klev) ! counter-gradient term for constituents
129!cAM      REAL cgs(klon,2:klev) ! counter-gradient star (cg/flux)
130      REAL obklen(klon)     ! Monin-Obukhov lengh
131!cAM      REAL ztvd, ztvu,
132      REAL zdu2
133      REAL therm(klon)      ! thermal virtual temperature excess
134      REAL trmb1(klon),trmb2(klon),trmb3(klon)
135!C  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
148!IM 091204 BEG
149      REAL a1,a2,a3
150!IM 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
154!c
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
170!cAM      REAL zxqs
171      REAL fac, pblmin, zmzp, term
172!c
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)
187!c
188!C ============================================================
189!C     Fonctions thermo implicites
190!C ============================================================
191!c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
192!c Tetens : pression partielle de vap d'eau e_sat(T)
193!c =================================================
194!C++ e_sat(T) = r2*exp( r3*(T-Tf)/(T-r4) ) id a r2*FOEWE
195!C++ avec :
196!C++ Tf = 273.16 K  (Temp de fusion de la glace)
197!C++ r2 = 611.14 Pa
198!C++ r3 = 17.269 (liquide) 21.875 (solide) adim
199!C++ r4 = 35.86             7.66           Kelvin
200!C++  q_sat = eps*e_sat/(p-(1-eps)*e_sat)
201!C++ derivᅵ :
202!C++ =========
203!C++                   r3*(Tf-r4)*q_sat(T,p)
204!C++ d_qsat_dT = --------------------------------
205!C++             (T-r4)^2*( 1-(1-eps)*e_sat(T)/p )
206!c++ pour zcvm5=Lv, c'est FOEDE
207!c++ Rq :(1.-REPS)*esarg/Parg id a RETV*Qsat
208!C     ------------------------------------------------------------------
209!c
210!c Initialisation
211      RLvCp = RLVTT/RCPD
212      REPS  = RD/RV
213
214!c
215!c      DO i = 1, klon
216!c         pcfh(i,1) = cd_h(i)
217!c         pcfm(i,1) = cd_m(i)
218!c      ENDDO
219!c      DO k = 2, klev
220!c      DO i = 1, klon
221!c         pcfh(i,k) = zkmin
222!c         pcfm(i,k) = zkmin
223!c         cgs(i,k) = 0.0
224!c         cgh(i,k) = 0.0
225!c         cgq(i,k) = 0.0
226!c      ENDDO
227!c      ENDDO
228!c
229!c Calculer les hauteurs de chaque couche
230!c (geopotentielle Int_dp/ro = Int_[Rd.T.dp/p] z = geop/g)
231!c  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
237!c                                 s(k) = [pplay(k)/ps]^kappa
238!c    + + + + + + + + + pplay  <-> s(k)   t  dp=pplay(k-1)-pplay(k)
239!c
240!c    -----------------  paprs <-> sig(k)
241!c
242!c    + + + + + + + + + pplay  <-> s(k-1)
243!c
244!c
245!c    + + + + + + + + + pplay  <-> s(1)   t  dp=paprs-pplay   z(1)
246!c
247!c    -----------------  paprs <-> sig(1)
248!c
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
257!c  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
258!c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
259!c +++  Determination des grandeurs de surface  +++++++++++++++++++++
260!c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
261!c  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
262      DO i = 1, knon
263!cAM         IF (thermcep) THEN
264!cAM           zdelta=MAX(0.,SIGN(1.,RTT-tsol(i)))
265!c           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
266!c           zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q(i,1))
267!cAM           zxqs= r2es * FOEEW(tsol(i),zdelta)/paprs(i,1)
268!cAM           zxqs=MIN(0.5,zxqs)
269!cAM           zcor=1./(1.-retv*zxqs)
270!cAM           zxqs=zxqs*zcor
271!cAM         ELSE
272!cAM           IF (tsol(i).LT.t_coup) THEN
273!cAM              zxqs = qsats(tsol(i)) / paprs(i,1)
274!cAM           ELSE
275!cAM              zxqs = qsatl(tsol(i)) / paprs(i,1)
276!cAM           ENDIF
277!cAM         ENDIF
278!c niveau de reference bulk; mais ici, c,a pourrait etre le niveau de ref du thermique
279!cAM        zx_alf1 = 1.0
280!cAM        zx_alf2 = 1.0 - zx_alf1
281!cAM        zxt = (t(i,1)+z(i,1)*RG/RCPD/(1.+RVTMP2*q(i,1)))
282!cAM     .        *(1.+RETV*q(i,1))*zx_alf1
283!cAM     .      + (t(i,2)+z(i,2)*RG/RCPD/(1.+RVTMP2*q(i,2)))
284!cAM     .        *(1.+RETV*q(i,2))*zx_alf2
285!cAM        zxu = u(i,1)*zx_alf1+u(i,2)*zx_alf2
286!cAM        zxv = v(i,1)*zx_alf1+v(i,2)*zx_alf2
287!cAM        zxq = q(i,1)*zx_alf1+q(i,2)*zx_alf2
288!cAM     
289!cAMAM           zxu = u10m(i)
290!cAMAM           zxv = v10m(i)
291!cAMAM           zxmod = 1.0+SQRT(zxu**2+zxv**2)
292!cAM Niveau de ref choisi a 2m
293        zxt = t2m(i)
294
295!c ***************************************************
296!c attention, il doit s'agir de <w'theta'>
297!c   ;Calcul de tcls virtuel et de w'theta'virtuel
298!c   ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
299!c   tcls=tcls*(1+.608*qcls)
300!c
301!c   ;Pour avoir w'theta',
302!c   ; il faut diviser par ro.Cp
303!c   Cp=Cpd*(1+0.84*qcls)
304!c   fcs=fcs/(ro_surf*Cp)
305!c   ;On transforme w'theta' en w'thetav'
306!c   Lv=(2.501-0.00237*(tcls-273.15))*1.E6
307!c   xle=xle/(ro_surf*Lv)
308!c   fcsv=fcs+.608*xle*tcls
309!c ***************************************************
310!cAM        khfs(i) = (tsol(i)*(1.+RETV*q(i,1))-zxt) *zxmod*cd_h(i)
311!cAM        kqfs(i) = (zxqs-zxq) *zxmod*cd_h(i) * beta(i)
312!cAM
313!cdif khfs est deja w't'_v / heatv(i) = khfs(i) + RETV*zxt*kqfs(i)
314!cAM calcule de Ro = paprs(i,1)/Rd zxt
315!cAM 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))
318!cAM   verifier que khfs et kqfs sont bien de la forme w'l'
319        heatv(i) = khfs(i) + 0.608*zxt*kqfs(i)
320!c a comparer aussi aux sorties de clqh : flux_T/RoCp et flux_q/RoLv
321!cAM        heatv(i) = khfs(i)
322!cAM ustar est en entree
323!cAM        taux = zxu *zxmod*cd_m(i)
324!cAM        tauy = zxv *zxmod*cd_m(i)
325!cAM        ustar(i) = SQRT(taux**2+tauy**2)
326!cAM        ustar(i) = MAX(SQRT(ustar(i)),0.01)
327!c Theta et qT du thermique sans exces (interpolin vers surf)
328!c chgt de niveau du thermique (jeudi 30/12/1999)
329!c (interpolation lineaire avant integration phi_h)
330!cAM        qT_th(i) = zxqs*beta(i) + 4./z(i,1)*(q(i,1)-zxqs*beta(i))
331!cAM        qT_th(i) = max(qT_th(i),q(i,1))
332        qT_th(i) = q2m(i)
333!cn The_th restera la Theta du thermique sans exces jusqu'a 2eme calcul
334!cn reste a regler convention P) pour Theta
335!c        The_th(i) = tsol(i) + 4./z(i,1)*(t(i,1)-tsol(i))
336!c     -                      + RLvCp*qT_th(i)
337!cAM        Th_th(i) = tsol(i) + 4./z(i,1)*(t(i,1)-tsol(i))
338        Th_th(i) = t2m(i)
339      ENDDO
340!c
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.
346!c 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
353!C
354!c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
355!C PBL height calculation:
356!C Search for level of pbl. Scan upward until the Richardson number between
357!C the first level and the current level exceeds the "critical" value.
358!C (bonne idee Nu de separer le Ric et l'exces de temp du thermique)
359!c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
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 ?
365!ctest     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)
368!c Theta_v environnement
369         zthvd=t(i,k)/s(i,k)*(1.+RETV*q(i,k))
370!c
371!c therm Theta_v sans exces (avec hypothese fausse de H&B, sinon,
372!c passer par Theta_e et virpot)
373!c         zthvu=t(i,1)/s(i,1)*(1.+RETV*q(i,1))
374!cAM         zthvu = Th_th(i)*(1.+RETV*q(i,1))
375         zthvu = Th_th(i)*(1.+RETV*qT_th(i))
376!c  Le Ri par Theta_v
377!cAM         rhino(i,k) = (z(i,k)-z(i,1))*RG*(zthvd-zthvu)
378!cAM     .               /(zdu2*0.5*(zthvd+zthvu))
379!cAM On a nveau de ref a 2m ???
380         rhino(i,k) = (z(i,k)-zref)*RG*(zthvd-zthvu)                                 &
381       &               /(zdu2*0.5*(zthvd+zthvu))
382!c
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))
386!c 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
396!C
397!C Set pbl height to maximum value where computation exceeds number of
398!C layers allowed
399!C
400      DO i = 1, knon
401        if (check(i)) pblh(i) = z(i,isommet)
402      ENDDO
403!C
404!C Improve estimate of pbl height for the unstable points.
405!C Find unstable points (sensible heat flux is upward):
406!C
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
416!C
417!C For the unstable case, compute velocity scale and the
418!C convective temperature excess:
419!C
420      DO i = 1, knon
421      IF (check(i)) THEN
422        phiminv(i) = (1.-binm*pblh(i)/obklen(i))**onet
423!c ***************************************************
424!c Wm ? et W* ? c'est la formule pour z/h < .1
425!c   ;Calcul de w* ;;
426!c   ;;;;;;;;;;;;;;;;
427!c   w_star=((g/tcls)*fcsv*z(ind))^(1/3.) [ou prendre la premiere approx de h)
428!c   ;; CALCUL DE wm ;;
429!c   ;;;;;;;;;;;;;;;;;;
430!c   ; Ici on considerera que l'on est dans la couche de surf jusqu'a 100m
431!c   ; On prend svt couche de surface=0.1*h mais on ne connait pas h
432!c   ;;;;;;;;;;;Dans la couche de surface
433!c   if (z(ind) le 20) then begin
434!c   Phim=(1.-15.*(z(ind)/L))^(-1/3.)
435!c   wm=u_star/Phim
436!c   ;;;;;;;;;;;En dehors de la couche de surface
437!c   endif else if (z(ind) gt 20) then begin
438!c   wm=(u_star^3+c1*w_star^3)^(1/3.)
439!c   endif
440!c ***************************************************
441        wm(i)= ustar(i)*phiminv(i)
442!c======================================================================
443!cvaleurs de Dominique Lambert de la campagne SEMAPHORE :
444!c <T'^2> = 100.T*^2; <q'^2> = 20.q*^2 a 10m
445!c <Tv'^2> = (1+1.2q).100.T* + 1.2Tv.sqrt(20*100).T*.q* + (.608*Tv)^2*20.q*^2;
446!c et dTetavS = sqrt(<Tv'^2>) ainsi calculee.
447!c avec : T*=<w'T'>_s/w* et q*=<w'q'>/w*
448!c !!! on peut donc utiliser w* pour les fluctuations <-> Lambert
449!c(leur corellation pourrait dependre de beta par ex)
450!c  if fcsv(i,j) gt 0 then begin
451!c    dTetavs=b1*(1.+2.*.608*q_10(i,j))*(fcs(i,j)/wm(i,j))^2+$
452!c    (.608*Thetav_10(i,j))^2*b2*(xle(i,j)/wm(i,j))^2+$
453!c    2.*.608*thetav_10(i,j)*sqrt(b1*b2)*(xle(i,j)/wm(i,j))*(fcs(i,j)/wm(i,j))
454!c    dqs=b2*(xle(i,j)/wm(i,j))^2
455!c    theta_s(i,j)=thetav_10(i,j)+sqrt(dTetavs)
456!c    q_s(i,j)=q_10(i,j)+sqrt(dqs)
457!c  endif else begin
458!c    Theta_s(i,j)=thetav_10(i,j)
459!c    q_s(i,j)=q_10(i,j)
460!c  endelse
461!c======================================================================
462!c
463!cHBTM        therm(i) = heatv(i)*fak/wm(i)
464!c forme Mathieu :
465        q_star = kqfs(i)/wm(i)
466        t_star = khfs(i)/wm(i)
467!IM 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
474!IM 091204 END
475!cAM Nveau cde ref 2m =>
476!cAM        therm(i) = sqrt( b1*(1.+2.*RETV*q(i,1))*t_star**2
477!cAM     +             + (RETV*T(i,1))**2*b2*q_star**2
478!cAM     +             + 2.*RETV*T(i,1)*b212*q_star*t_star
479!cAM     +                 )
480!IM 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
492!IM 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                           &
495!IM 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       &                 )
498!c
499!c Theta et qT du thermique (forme H&B) avec exces
500!c (attention, on ajoute therm(i) qui est virtuelle ...)
501!c pourquoi pas sqrt(b1)*t_star ?
502!c        dqs = b2sr*kqfs(i)/wm(i)
503        qT_th(i) = qT_th(i)  + b2sr*q_star
504!cnew on differre le calcul de Theta_e
505!c        The_th(i) = The_th(i) + therm(i) + RLvCp*qT_th(i)
506!c 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
510!C
511!c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
512!c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
513!C ++ Improve pblh estimate for unstable conditions using the +++++++
514!C ++          convective temperature excess :                +++++++
515!c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
516!c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
517!C
518      DO k = 2, isommet
519      DO i = 1, knon
520      IF (check(i)) THEN
521!ctest     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)
524!c Theta_v environnement
525         zthvd=t(i,k)/s(i,k)*(1.+RETV*q(i,k))
526!c
527!c et therm Theta_v (avec hypothese de constance de H&B,
528!c         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
531!c
532!c  Le Ri par Theta_v
533!cAM Niveau de ref 2m
534!cAM         rhino(i,k) = (z(i,k)-z(i,1))*RG*(zthvd-zthvu)
535!cAM     .               /(zdu2*0.5*(zthvd+zthvu))
536         rhino(i,k) = (z(i,k)-zref)*RG*(zthvd-zthvu)                                 &
537       &               /(zdu2*0.5*(zthvd+zthvu))
538!c
539!c
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))
543!c 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.
548!IM 170305 BEG
549      IF(1.EQ.0) THEN
550!c 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
562!IM 170305 END
563!c             q_star = kqfs(i)/wm(i)
564!c             t_star = khfs(i)/wm(i)
565!c             trmb1(i) = b1*(1.+2.*RETV*q(i,1))*t_star**2
566!c             trmb2(i) = (RETV*T(i,1))**2*b2*q_star**2
567!c Omega now   trmb3(i) = 2.*RETV*T(i,1)*b212*q_star*t_star
568         ENDIF
569      ENDIF
570      ENDDO
571      ENDDO
572!C
573!C Set pbl height to maximum value where computation exceeds number of
574!C layers allowed
575!C
576      DO i = 1, knon
577        if (check(i)) pblh(i) = z(i,isommet)
578      ENDDO
579!C
580!C PBL height must be greater than some minimum mechanical mixing depth
581!C Several investigators have proposed minimum mechanical mixing depth
582!C relationships as a function of the local friction velocity, u*.  We
583!C make use of a linear relationship of the form h = c u* where c=700.
584!C The scaling arguments that give rise to this relationship most often
585!C represent the coefficient c as some constant over the local coriolis
586!C parameter.  Here we make use of the experimental results of Koracin
587!C and Berkowicz (1988) [BLM, Vol 43] for wich they recommend 0.07/f
588!C where f was evaluated at 39.5 N and 52 N.  Thus we use a typical mid
589!C latitude value for f so that c = 0.07/f = 700.
590!C
591      DO i = 1, knon
592        pblmin  = 700.0*ustar(i)
593        pblh(i) = MAX(pblh(i),pblmin)
594!c 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
599!C ********************************************************************
600!C  pblh is now available; do preparation for diffusivity calculation :
601!C ********************************************************************
602      DO i = 1, knon
603        check(i) = .TRUE.
604        Zsat(i)   = .FALSE.
605!c 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
613!C
614!C Do additional preparation for unstable cases only, set temperature
615!C and moisture perturbations depending on stability.
616!C *** Rq: les formule sont prises dans leur forme CS ***
617        IF (unstbl(i)) THEN
618!cAM Niveau de ref du thermique
619!cAM          zxt=(t(i,1)-z(i,1)*0.5*RG/RCPD/(1.+RVTMP2*q(i,1)))
620!cAM     .         *(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
630!c Computes Theta_e for thermal (all cases : to be modified)
631!c   attention ajout therm(i) = virtuelle
632        The_th(i) = Th_th(i) + therm(i) + RLvCp*qT_th(i)
633!c ou:    The_th(i) = Th_th(i) + sqrt(b1)*khfs(i)/wm(i) + RLvCp*qT_th(i)
634      ENDDO
635
636!C Main level loop to compute the diffusivities and
637!C counter-gradient terms:
638!C
639      DO 1000 k = 2, isommet
640!C
641!C Find levels within boundary layer:
642!C
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))
651!C debug
652!c          if (i.EQ.1864) then
653!c             print*,'i,pblh(1864),obklen(1864)',i,pblh(i),obklen(i)
654!c          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
660!C
661!C stblev for points zm < plbh and stable and neutral
662!C unslev for points zm < plbh and unstable
663!C
664            IF (unstbl(i)) THEN
665              unslev(i) = .TRUE.
666            ELSE
667              stblev(i) = .TRUE.
668            ENDIF
669          ENDIF
670        ENDDO
671!c        print*,'fin calcul niveaux'
672!C
673!C Stable and neutral points; set diffusivities; counter-gradient
674!C terms zero for stable case:
675!C
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
683!c            pcfm(i,k) = pblk(i)
684!c            pcfh(i,k) = pcfm(i,k)
685          ENDIF
686        ENDDO
687!C
688!C unssrf, unstable within surface layer of pbl
689!C unsout, unstable within outer   layer of pbl
690!C
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
702!C
703!C Unstable for surface layer; counter-gradient terms zero
704!C
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
712!c        print*,'fin counter-gradient terms zero'
713!C
714!C Unstable for outer layer; counter-gradient terms non-zero:
715!C
716        DO i = 1, knon
717          IF (unsout(i)) THEN
718            pblk(i) = fak2(i)*zh(i)*zzh(i)
719!c            cgs(i,k) = fak3(i)/(pblh(i)*wm(i))
720!c            cgh(i,k) = khfs(i)*cgs(i,k)
721            pr(i) = phiminv(i)/phihinv(i) + ccon*fak3(i)/fak
722!c            cgq(i,k) = kqfs(i)*cgs(i,k)
723          ENDIF
724        ENDDO
725!c        print*,'fin counter-gradient terms non zero'
726!C
727!C For all unstable layers, compute diffusivities and ctrgrad ter m
728!C
729!c        DO i = 1, knon
730!c        IF (unslev(i)) THEN
731!c            pcfm(i,k) = pblk(i)
732!c            pcfh(i,k) = pblk(i)/pr(i)
733!c etc cf original
734!c        ENDIF
735!c        ENDDO
736!C
737!C For all layers, compute integral info and CTEI
738!C
739        DO i = 1, knon
740        if (check(i).or.omegafl(i)) then
741          if (.not.Zsat(i)) then
742!c            Th2 = The_th(i) - RLvCp*qT_th(i)
743            Th2 = Th_th(i)
744            T2 = Th2*s(i,k)
745!c 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
751!c
752            if (qqsat.lt.qT_th(i)) then
753!c 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
763!c
764            qsatbef(i) = qqsat    ! bug dans la version orig ???
765          endif
766!camn ???? cette ligne a deja ete faite normalement ?
767        endif
768!c            print*,'hbtm2 i,k=',i,k
769        ENDDO
770 1000 continue           ! end of level loop
771!IM 170305 BEG
772        IF(1.EQ.0) THEN
773            print*,'hbtm2  ok'
774        ENDIF !(1.EQ.0) THEN
775!IM 170305 END
776      RETURN
777      END
Note: See TracBrowser for help on using the repository browser.