source: LMDZ5/trunk/libf/phylmd/hbtm.F @ 1983

Last change on this file since 1983 was 1943, checked in by idelkadi, 11 years ago

Changes concerinng the Hourdin et al. 2002 version of the thermals and
Ayotte cases

calltherm.F90

call to thermcell_2002 modified :
1) iflag_thermals changed from 1 to >= 1000 ; iflag_thermals-1000
controls sub options.
2) thermals w and fraction in output files.

hbtm.F

Singularity in the 1/heat dependency of the Monin Obukov length L when
heat=0. Since 1/L is used rather than L, it is preferable to compute
directly L. There is a dependency in 1/u* then which is treated with a
threshold value.
(+ some cleaning in the syntax).

lmdz1d.F

If nday<0, -nday is used as the total number of time steps of the
simulation.
The option with imposed wtsurf and wqsurf read in lmdz1d.def was not
active anymore.
< IF(.NOT.ok_flux_surf) THEN
changed to

IF(.NOT.ok_flux_surf.or.max(abs(wtsurf),abs(wqsurf))>0.) THEN

before

fsens=-wtsurf*rcpd*rho(1)
flat=-wqsurf*rlvtt*rho(1)

phys_output_write.F90 et phys_local_var_mod.F90

Removing the d_u_ajsb contribution to duthe (same for dv).
Those tendencies are not computed by the model ...
< zx_tmp_fi3d(1:klon,1:klev)=d_u_ajs(1:klon,1:klev)/pdtphys - &
< d_u_ajsb(1:klon,1:klev)/pdtphys
---

zx_tmp_fi3d(1:klon,1:klev)=d_u_ajs(1:klon,1:klev)/pdtphys

! d_u_ajsb(1:klon,1:klev)/pdtphys

thermcell_dq.F90, thermcell_main.F90

Some cleaning

thermcell_old.F

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