Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

hbtm_mod.f90

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

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

Download in other formats:

Original Format