Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

hbtm_mod.F90 @ 4109

Last change on this file since 4109 was 3774, checked in by lguez, 4 years ago

Bug fix: encapsulate hbtm in a module

Fixes bug introduced in commit [3772]: therm was declared as an
assumed-shape dummy argument, which requires an explicit interface.

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:keywords set to Author Date Id Revision

File size: 29.3 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format