Context Navigation

lmdz_lscp_tools.F90 @ 5171

Last change on this file since 5171 was 5170, checked in by evignon, 10 months ago
petites correction dans la nouvelle routine de partitionnement de phase des nuages de Lea Raillard
File size: 30.0 KB

Line
1	MODULE lmdz_lscp_tools
2
3	IMPLICIT NONE
4
5	CONTAINS
6
7	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
8	SUBROUTINE FALLICE_VELOCITY(klon,iwc,temp,rho,pres,ptconv,velo)
9
10	! Ref:
11	! Stubenrauch, C. J., Bonazzola, M.,
12	! Protopapadaki, S. E., & Musat, I. (2019).
13	! New cloud system metrics to assess bulk
14	! ice cloud schemes in a GCM. Journal of
15	! Advances in Modeling Earth Systems, 11,
16	! 3212–3234. https://doi.org/10.1029/2019MS001642
17
18	use lmdz_lscp_ini, only: iflag_vice, ffallv_con, ffallv_lsc
19	use lmdz_lscp_ini, only: cice_velo, dice_velo
20
21	IMPLICIT NONE
22
23	INTEGER, INTENT(IN) :: klon
24	REAL, INTENT(IN), DIMENSION(klon) :: iwc ! specific ice water content [kg/m3]
25	REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature [K]
26	REAL, INTENT(IN), DIMENSION(klon) :: rho ! dry air density [kg/m3]
27	REAL, INTENT(IN), DIMENSION(klon) :: pres ! air pressure [Pa]
28	LOGICAL, INTENT(IN), DIMENSION(klon) :: ptconv ! convective point [-]
29
30	REAL, INTENT(OUT), DIMENSION(klon) :: velo ! fallspeed velocity of crystals [m/s]
31
32
33	INTEGER i
34	REAL logvm,iwcg,tempc,phpa,fallv_tun
35	REAL m2ice, m2snow, vmice, vmsnow
36	REAL aice, bice, asnow, bsnow
37
38
39	DO i=1,klon
40
41	IF (ptconv(i)) THEN
42	fallv_tun=ffallv_con
43	ELSE
44	fallv_tun=ffallv_lsc
45	ENDIF
46
47	tempc=temp(i)-273.15 ! celcius temp
48	iwcg=MAX(iwc(i)*1000.,1E-3) ! iwc in g/m3. We set a minimum value to prevent from division by 0
49	phpa=pres(i)/100. ! pressure in hPa
50
51	IF (iflag_vice .EQ. 1) THEN
52	! so-called 'empirical parameterization' in Stubenrauch et al. 2019
53	if (tempc .GE. -60.0) then
54
55	logvm= -0.0000414122tempctempclog(iwcg)-0.00538922tempc*log(iwcg) &
56	-0.0516344log(iwcg)+0.00216078tempc + 1.9714
57	velo(i)=exp(logvm)
58	else
59	velo(i)=65.0(iwcg0.2)(150./phpa)**0.15
60	endif
61
62	velo(i)=fallv_tun*velo(i)/100.0 ! from cm/s to m/s
63
64	ELSE IF (iflag_vice .EQ. 2) THEN
65	! so called PSDM empirical coherent bulk ice scheme in Stubenrauch et al. 2019
66	aice=0.587
67	bice=2.45
68	asnow=0.0444
69	bsnow=2.1
70
71	m2ice=((iwcg0.001/aice)/(exp(13.6-bice7.76+0.479bice2) &
72	exp((-0.0361+bice0.0151+0.00149bice*2)tempc))) &
73	*(1./(0.807+bice0.00581+0.0457bice*2))
74
75	vmice=100.1042.4exp(13.6-(bice+1)7.76+0.479(bice+1.)*2)exp((-0.0361+ &
76	(bice+1.)0.0151+0.00149(bice+1.)*2)tempc) &
77	(m2ice(0.807+(bice+1.)0.00581+0.0457(bice+1.)2))/(iwcg0.001/aice)
78
79
80	vmice=vmice((1000./phpa)*0.2)
81
82	m2snow=((iwcg0.001/asnow)/(exp(13.6-bsnow7.76+0.479bsnow2) &
83	exp((-0.0361+bsnow0.0151+0.00149bsnow*2)tempc))) &
84	*(1./(0.807+bsnow0.00581+0.0457bsnow*2))
85
86
87	vmsnow=100.14.3exp(13.6-(bsnow+.416)7.76+0.479(bsnow+.416)**2)&
88	exp((-0.0361+(bsnow+.416)0.0151+0.00149(bsnow+.416)2)tempc)&
89	(m2snow(0.807+(bsnow+.416)0.00581+0.0457(bsnow+.416)2))/(iwcg0.001/asnow)
90
91	vmsnow=vmsnow((1000./phpa)*0.35)
92	velo(i)=fallv_tun*min(vmsnow,vmice)/100. ! to m/s
93
94	ELSE
95	! By default, fallspeed velocity of ice crystals according to Heymsfield & Donner 1990
96	velo(i) = fallv_tuncice_velo((iwcg/1000.)**dice_velo)
97	ENDIF
98	ENDDO
99
100	END SUBROUTINE FALLICE_VELOCITY
101	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
102
103	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
104	SUBROUTINE ICEFRAC_LSCP(klon, temp, iflag_ice_thermo, distcltop, temp_cltop, icefrac, dicefracdT)
105	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
106
107	! Compute the ice fraction 1-xliq (see e.g.
108	! Doutriaux-Boucher & Quaas 2004, section 2.2.)
109	! as a function of temperature
110	! see also Fig 3 of Madeleine et al. 2020, JAMES
111	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
112
113
114	USE print_control_mod, ONLY: lunout, prt_level
115	USE lmdz_lscp_ini, ONLY: t_glace_min, t_glace_max, exposant_glace, iflag_t_glace
116	USE lmdz_lscp_ini, ONLY : RTT, dist_liq, temp_nowater
117
118	IMPLICIT NONE
119
120
121	INTEGER, INTENT(IN) :: klon ! number of horizontal grid points
122	REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature
123	REAL, INTENT(IN), DIMENSION(klon) :: distcltop ! distance to cloud top
124	REAL, INTENT(IN), DIMENSION(klon) :: temp_cltop ! temperature of cloud top
125	INTEGER, INTENT(IN) :: iflag_ice_thermo
126	REAL, INTENT(OUT), DIMENSION(klon) :: icefrac
127	REAL, INTENT(OUT), DIMENSION(klon) :: dicefracdT
128
129
130	INTEGER i
131	REAL liqfrac_tmp, dicefrac_tmp
132	REAL Dv, denomdep,beta,qsi,dqsidt
133	LOGICAL ice_thermo
134
135	CHARACTER (len = 20) :: modname = 'lscp_tools'
136	CHARACTER (len = 80) :: abort_message
137
138	IF ((iflag_t_glace.LT.2)) THEN !.OR. (iflag_t_glace.GT.6)) THEN
139	abort_message = 'lscp cannot be used if iflag_t_glace<2 or >6'
140	CALL abort_physic(modname,abort_message,1)
141	ENDIF
142
143	IF (.NOT.((iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3))) THEN
144	abort_message = 'lscp cannot be used without ice thermodynamics'
145	CALL abort_physic(modname,abort_message,1)
146	ENDIF
147
148
149	DO i=1,klon
150
151	! old function with sole dependence upon temperature
152	IF (iflag_t_glace .EQ. 2) THEN
153	liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
154	liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
155	icefrac(i) = (1.0-liqfrac_tmp)**exposant_glace
156	IF (icefrac(i) .GT.0.) THEN
157	dicefracdT(i)= exposant_glace * (icefrac(i)**(exposant_glace-1.)) &
158	/ (t_glace_min - t_glace_max)
159	ENDIF
160
161	IF ((icefrac(i).EQ.0).OR.(icefrac(i).EQ.1)) THEN
162	dicefracdT(i)=0.
163	ENDIF
164
165	ENDIF
166
167	! function of temperature used in CMIP6 physics
168	IF (iflag_t_glace .EQ. 3) THEN
169	liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
170	liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
171	icefrac(i) = 1.0-liqfrac_tmp**exposant_glace
172	IF ((icefrac(i) .GT.0.) .AND. (liqfrac_tmp .GT. 0.)) THEN
173	dicefracdT(i)= exposant_glace * ((liqfrac_tmp)**(exposant_glace-1.)) &
174	/ (t_glace_min - t_glace_max)
175	ELSE
176	dicefracdT(i)=0.
177	ENDIF
178	ENDIF
179
180	! for iflag_t_glace .GE. 4, the liquid fraction depends upon temperature at cloud top
181	! and then decreases with decreasing height
182
183	!with linear function of temperature at cloud top
184	IF (iflag_t_glace .EQ. 4) THEN
185	liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
186	liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
187	icefrac(i) = MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.)
188	dicefrac_tmp = - temp(i)/(t_glace_max-t_glace_min)
189	dicefracdT(i) = dicefrac_tmp*exp(-distcltop(i)/dist_liq)
190	IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN
191	dicefracdT(i) = 0.
192	ENDIF
193	ENDIF
194
195	! with CMIP6 function of temperature at cloud top
196	IF ((iflag_t_glace .EQ. 5) .OR. (iflag_t_glace .EQ. 7)) THEN
197	liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
198	liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
199	liqfrac_tmp = liqfrac_tmp**exposant_glace
200	icefrac(i) = MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.)
201	IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN
202	dicefracdT(i) = 0.
203	ELSE
204	dicefracdT(i) = exposant_glace((liqfrac_tmp)*(exposant_glace-1.))/(t_glace_min- t_glace_max) &
205	*exp(-distcltop(i)/dist_liq)
206	ENDIF
207	ENDIF
208
209	! with modified function of temperature at cloud top
210	! to get largere values around 260 K, works well with t_glace_min = 241K
211	IF (iflag_t_glace .EQ. 6) THEN
212	IF (temp(i) .GT. t_glace_max) THEN
213	liqfrac_tmp = 1.
214	ELSE
215	liqfrac_tmp = -((temp(i)-t_glace_max) / (t_glace_max-t_glace_min))**2+1.
216	ENDIF
217	liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
218	icefrac(i) = MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.)
219	IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN
220	dicefracdT(i) = 0.
221	ELSE
222	dicefracdT(i) = 2*((temp(i)-t_glace_max) / (t_glace_max-t_glace_min))/(t_glace_max-t_glace_min) &
223	*exp(-distcltop(i)/dist_liq)
224	ENDIF
225	ENDIF
226
227	! if temperature of cloud top <-40°C,
228	IF (iflag_t_glace .GE. 4) THEN
229	IF ((temp_cltop(i) .LE. temp_nowater) .AND. (temp(i) .LE. t_glace_max)) THEN
230	icefrac(i) = 1.
231	dicefracdT(i) = 0.
232	ENDIF
233	ENDIF
234
235
236	ENDDO ! klon
237	RETURN
238
239	END SUBROUTINE ICEFRAC_LSCP
240	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
241
242	SUBROUTINE ICEFRAC_LSCP_TURB(klon, dtime, temp, pplay, paprsdn, paprsup, qice_ini, snowcld, qtot_incl, cldfra, tke, tke_dissip, qliq, qvap_cld, qice, icefrac, dicefracdT, cldfraliq, sigma2_icefracturb, mean_icefracturb)
243	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
244	! Compute the liquid, ice and vapour content (+ice fraction) based
245	! on turbulence (see Fields 2014, Furtado 2016, Raillard 2025)
246	! L.Raillard (30/08/24)
247	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
248
249
250	USE lmdz_lscp_ini, ONLY : prt_level, lunout
251	USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RV, RPI
252	USE lmdz_lscp_ini, ONLY : seuil_neb, temp_nowater
253	USE lmdz_lscp_ini, ONLY : tau_mixenv, lmix_mpc, naero5, gamma_snwretro, gamma_taud, capa_crystal
254	USE lmdz_lscp_ini, ONLY : eps
255
256	IMPLICIT NONE
257
258	INTEGER, INTENT(IN) :: klon !--number of horizontal grid points
259	REAL, INTENT(IN) :: dtime !--time step [s]
260
261	REAL, INTENT(IN), DIMENSION(klon) :: temp !--temperature
262	REAL, INTENT(IN), DIMENSION(klon) :: pplay !--pressure in the middle of the layer [Pa]
263	REAL, INTENT(IN), DIMENSION(klon) :: paprsdn !--pressure at the bottom interface of the layer [Pa]
264	REAL, INTENT(IN), DIMENSION(klon) :: paprsup !--pressure at the top interface of the layer [Pa]
265	REAL, INTENT(IN), DIMENSION(klon) :: qtot_incl !--specific total cloud water content, in-cloud content [kg/kg]
266	REAL, INTENT(IN), DIMENSION(klon) :: cldfra !--cloud fraction in gridbox [-]
267	REAL, INTENT(IN), DIMENSION(klon) :: tke !--turbulent kinetic energy [m2/s2]
268	REAL, INTENT(IN), DIMENSION(klon) :: tke_dissip !--TKE dissipation [m2/s3]
269
270	REAL, INTENT(IN), DIMENSION(klon) :: qice_ini !--initial specific ice content gridbox-mean [kg/kg]
271	REAL, INTENT(IN), DIMENSION(klon) :: snowcld
272	REAL, INTENT(OUT), DIMENSION(klon) :: qliq !--specific liquid content gridbox-mean [kg/kg]
273	REAL, INTENT(OUT), DIMENSION(klon) :: qvap_cld !--specific cloud vapor content, gridbox-mean [kg/kg]
274	REAL, INTENT(OUT), DIMENSION(klon) :: qice !--specific ice content gridbox-mean [kg/kg]
275	REAL, INTENT(OUT), DIMENSION(klon) :: icefrac !--fraction of ice in condensed water [-]
276	REAL, INTENT(OUT), DIMENSION(klon) :: dicefracdT
277
278	REAL, INTENT(OUT), DIMENSION(klon) :: cldfraliq !--fraction of cldfra which is liquid only
279	REAL, INTENT(OUT), DIMENSION(klon) :: sigma2_icefracturb !--Temporary
280	REAL, INTENT(OUT), DIMENSION(klon) :: mean_icefracturb !--Temporary
281
282	REAL, DIMENSION(klon) :: qzero, qsatl, dqsatl, qsati, dqsati !--specific humidity saturation values
283	INTEGER :: i
284
285	REAL :: qvap_incl, qice_incl, qliq_incl, qiceini_incl !--In-cloud specific quantities [kg/kg]
286	REAL :: qsnowcld_incl
287	!REAL :: capa_crystal !--Capacitance of ice crystals [-]
288	REAL :: water_vapor_diff !--Water-vapour diffusion coefficient in air [m2/s] (function of T&P)
289	REAL :: air_thermal_conduct !--Thermal conductivity of air [J/m/K/s] (function of T)
290	REAL :: C0 !--Lagrangian structure function [-]
291	REAL :: tau_mixingenv
292	REAL :: tau_dissipturb
293	REAL :: invtau_phaserelax
294	REAL :: sigma2_pdf, mean_pdf
295	REAL :: ai, bi, B0
296	REAL :: sursat_iceliq
297	REAL :: sursat_env
298	REAL :: liqfra_max
299	REAL :: sursat_iceext
300	REAL :: nb_crystals !--number concentration of ice crystals [#/m3]
301	REAL :: moment1_PSD !--1st moment of ice PSD
302	REAL :: N0_PSD, lambda_PSD !--parameters of the exponential PSD
303
304	REAL :: rho_ice !--ice density [kg/m3]
305	REAL :: cldfra1D
306	REAL :: deltaz, rho_air
307	REAL :: psati !--saturation vapor pressure wrt i [Pa]
308
309	C0 = 10. !--value assumed in Field2014
310	rho_ice = 950.
311	sursat_iceext = -0.1
312	!capa_crystal = 1. !r_ice
313	qzero(:) = 0.
314	cldfraliq(:) = 0.
315	icefrac(:) = 0.
316	dicefracdT(:) = 0.
317
318	sigma2_icefracturb(:) = 0.
319	mean_icefracturb(:) = 0.
320
321	!--wrt liquid water
322	CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,1,.false.,qsatl(:),dqsatl(:))
323	!--wrt ice
324	CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,2,.false.,qsati(:),dqsati(:))
325
326
327	DO i=1,klon
328
329
330	rho_air = pplay(i) / temp(i) / RD
331	!deltaz = ( paprsdn(i) - paprsup(i) ) / RG / rho_air(i)
332	! because cldfra is intent in, but can be locally modified due to test
333	cldfra1D = cldfra(i)
334	IF (cldfra(i) .LE. 0.) THEN
335	qvap_cld(i) = 0.
336	qliq(i) = 0.
337	qice(i) = 0.
338	cldfraliq(i) = 0.
339	icefrac(i) = 0.
340	dicefracdT(i) = 0.
341
342	! If there is a cloud
343	ELSE
344	IF (cldfra(i) .GE. 1.0) THEN
345	cldfra1D = 1.0
346	END IF
347
348	! T>0°C, no ice allowed
349	IF ( temp(i) .GE. RTT ) THEN
350	qvap_cld(i) = qsatl(i) * cldfra1D
351	qliq(i) = MAX(0.0,qtot_incl(i)-qsatl(i)) * cldfra1D
352	qice(i) = 0.
353	cldfraliq(i) = 1.
354	icefrac(i) = 0.
355	dicefracdT(i) = 0.
356
357	! T<-38°C, no liquid allowed
358	ELSE IF ( temp(i) .LE. temp_nowater) THEN
359	qvap_cld(i) = qsati(i) * cldfra1D
360	qliq(i) = 0.
361	qice(i) = MAX(0.0,qtot_incl(i)-qsati(i)) * cldfra1D
362	cldfraliq(i) = 0.
363	icefrac(i) = 1.
364	dicefracdT(i) = 0.
365
366	! MPC temperature
367	ELSE
368	! Not enough TKE
369	IF ( tke_dissip(i) .LE. eps ) THEN
370	qvap_cld(i) = qsati(i) * cldfra1D
371	qliq(i) = 0.
372	qice(i) = MAX(0.,qtot_incl(i)-qsati(i)) * cldfra1D
373	cldfraliq(i) = 0.
374	icefrac(i) = 1.
375	dicefracdT(i) = 0.
376
377	! Enough TKE
378	ELSE
379	print*,"MOUCHOIRACTIVE"
380	!---------------------------------------------------------
381	!-- ICE SUPERSATURATION PDF
382	!---------------------------------------------------------
383	!--If -38°C< T <0°C and there is enough turbulence,
384	!--we compute the cloud liquid properties with a Gaussian PDF
385	!--of ice supersaturation F(Si) (Field2014, Furtado2016).
386	!--Parameters of the PDF are function of turbulence and
387	!--microphysics/existing ice.
388
389	sursat_iceliq = qsatl(i)/qsati(i) - 1.
390	psati = qsati(i) * pplay(i) / (RD/RV)
391
392	!-------------- MICROPHYSICAL TERMS --------------
393	!--We assume an exponential ice PSD whose parameters
394	!--are computed following Morrison&Gettelman 2008
395	!--Ice number density is assumed equals to INP density
396	!--which is a function of temperature (DeMott 2010)
397	!--bi and B0 are microphysical function characterizing
398	!--vapor/ice interactions
399	!--tau_phase_relax is the typical time of vapor deposition
400	!--onto ice crystals
401
402	qiceini_incl = qice_ini(i) / cldfra1D
403	qsnowcld_incl = snowcld(i) * RG * dtime / ( paprsdn(i) - paprsup(i) ) / cldfra1D
404	sursat_env = max(0., (qtot_incl(i) - qiceini_incl)/qsati(i) - 1.)
405	IF ( qiceini_incl .GT. eps ) THEN
406	nb_crystals = 1.e3 * 5.94e-5 * ( RTT - temp(i) )*3.33 naero5*(0.0264(RTT-temp(i))+0.0033)
407	lambda_PSD = ( (RPIrho_icenb_crystals) / (rho_air(qiceini_incl + gamma_snwretro qsnowcld_incl)) ) ** (1./3.)
408	N0_PSD = nb_crystals * lambda_PSD
409	moment1_PSD = N0_PSD/lambda_PSD**2
410	ELSE
411	moment1_PSD = 0.
412	ENDIF
413
414	!--Formulae for air thermal conductivity and water vapor diffusivity
415	!--comes respectively from Beard and Pruppacher (1971)
416	!--and Hall and Pruppacher (1976)
417
418	air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184
419	water_vapor_diff = 2.111e-5 ( temp(i) / RTT )*1.94 ( 101325 / pplay(i) )
420
421	bi = 1./((qsati(i)+qsatl(i))/2.) + RLSTT2 / RCPD / RV / temp(i)2
422	B0 = 4. * RPI * capa_crystal * 1. / ( RLSTT2 / air_thermal_conduct / RV / temp(i)2 &
423	+ RV * temp(i) / psati / water_vapor_diff )
424
425	invtau_phaserelax = (bi * B0 * moment1_PSD )
426
427	! Old way of estimating moment1 : spherical crystals + monodisperse PSD
428	! nb_crystals = rho_air * qiceini_incl / ( 4. / 3. * RPI * r_ice*3. rho_ice )
429	! moment1_PSD = nb_crystals * r_ice
430
431	!----------------- TURBULENT SOURCE/SINK TERMS -----------------
432	!--Tau_mixingenv is the time needed to homogeneize the parcel
433	!--with its environment by turbulent diffusion over the parcel
434	!--length scale
435	!--if lmix_mpc <0, tau_mixigenv value is prescribed
436	!--else tau_mixigenv value is derived from tke_dissip and lmix_mpc
437	!--Tau_dissipturb is the time needed turbulence to decay due to
438	!--viscosity
439
440	ai = RG / RD / temp(i) * ( RD * RLSTT / RCPD / RV / temp(i) - 1. )
441	IF ( lmix_mpc .GT. 0 ) THEN
442	tau_mixingenv = ( lmix_mpc2. / tke_dissip(i) )(1./3.)
443	ELSE
444	tau_mixingenv = tau_mixenv
445	ENDIF
446
447	tau_dissipturb = gamma_taud * 2. * 2./3. * tke(i) / tke_dissip(i) / C0
448
449	!--------------------- PDF COMPUTATIONS ---------------------
450	!--Formulae for sigma2_pdf (variance), mean of PDF in Furtado2016
451	!--cloud liquid fraction and in-cloud liquid content are given
452	!--by integrating resp. F(Si) and Si*F(Si)
453	!--Liquid is limited by the available water vapor trough a
454	!--maximal liquid fraction
455
456	liqfra_max = MAX(0., (MIN (1.,( qtot_incl(i) - qiceini_incl - qsati(i) * (1 + sursat_iceext ) ) / ( qsatl(i) - qsati(i) ) ) ) )
457	sigma2_pdf = 1./2. * ( ai*2 ) 2./3. * tke(i) * tau_dissipturb / ( invtau_phaserelax + 1./tau_mixingenv )
458	mean_pdf = sursat_env * 1./tau_mixingenv / ( invtau_phaserelax + 1./tau_mixingenv )
459	cldfraliq(i) = 0.5 * (1. - erf( ( sursat_iceliq - mean_pdf) / (SQRT(2.* sigma2_pdf) ) ) )
460	IF (cldfraliq(i) .GT. liqfra_max) THEN
461	cldfraliq(i) = liqfra_max
462	ENDIF
463
464	qliq_incl = qsati(i) * SQRT(sigma2_pdf) / SQRT(2.RPI) EXP( -1.(sursat_iceliq - mean_pdf)2. / (2.sigma2_pdf) ) &
465	- qsati(i) * cldfraliq(i) * (sursat_iceliq - mean_pdf )
466
467	sigma2_icefracturb(i)= sigma2_pdf
468	mean_icefracturb(i) = mean_pdf
469
470	!------------ SPECIFIC VAPOR CONTENT AND WATER CONSERVATION ------------
471
472	IF ( (qliq_incl .LE. eps) .OR. (cldfraliq(i) .LE. eps) ) THEN
473	qliq_incl = 0.
474	cldfraliq(i) = 0.
475	END IF
476
477	!--Specific humidity is the max between qsati and the weighted mean between
478	!--qv in MPC patches and qv in ice-only parts. We assume that MPC parts are
479	!--always at qsatl and ice-only parts slightly subsaturated (qsati*sursat_iceext+1)
480	!--The whole cloud can therefore be supersaturated but never subsaturated.
481
482	qvap_incl = MAX(qsati(i), ( 1. - cldfraliq(i) ) * (sursat_iceext + 1.) * qsati(i) + cldfraliq(i) * qsatl(i) )
483
484
485	IF ( qvap_incl .GE. qtot_incl(i) ) THEN
486	qvap_incl = qsati(i)
487	qliq_incl = qtot_incl(i) - qvap_incl
488	qice_incl = 0.
489
490	ELSEIF ( (qvap_incl + qliq_incl) .GE. qtot_incl(i) ) THEN
491	qliq_incl = MAX(0.0,qtot_incl(i) - qvap_incl)
492	qice_incl = 0.
493	ELSE
494	qice_incl = qtot_incl(i) - qvap_incl - qliq_incl
495	END IF
496
497	qvap_cld(i) = qvap_incl * cldfra1D
498	qliq(i) = qliq_incl * cldfra1D
499	qice(i) = qice_incl * cldfra1D
500	icefrac(i) = qice(i) / ( qice(i) + qliq(i) )
501	dicefracdT(i) = 0.
502	!print*,'MPC turb'
503
504	END IF ! Enough TKE
505
506	END IF ! ! MPC temperature
507
508	END IF ! cldfra
509
510	ENDDO ! klon
511	END SUBROUTINE ICEFRAC_LSCP_TURB
512	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
513
514
515	SUBROUTINE CALC_QSAT_ECMWF(klon,temp,qtot,pressure,tref,phase,flagth,qs,dqs)
516	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
517	! Calculate qsat following ECMWF method
518	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
519
520
521	IMPLICIT NONE
522
523	include "YOMCST.h"
524	include "YOETHF.h"
525	include "FCTTRE.h"
526
527	INTEGER, INTENT(IN) :: klon ! number of horizontal grid points
528	REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature in K
529	REAL, INTENT(IN), DIMENSION(klon) :: qtot ! total specific water in kg/kg
530	REAL, INTENT(IN), DIMENSION(klon) :: pressure ! pressure in Pa
531	REAL, INTENT(IN) :: tref ! reference temperature in K
532	LOGICAL, INTENT(IN) :: flagth ! flag for qsat calculation for thermals
533	INTEGER, INTENT(IN) :: phase
534	! phase: 0=depend on temperature sign (temp>tref -> liquid, temp<tref, solid)
535	! 1=liquid
536	! 2=solid
537
538	REAL, INTENT(OUT), DIMENSION(klon) :: qs ! saturation specific humidity [kg/kg]
539	REAL, INTENT(OUT), DIMENSION(klon) :: dqs ! derivation of saturation specific humidity wrt T
540
541	REAL delta, cor, cvm5
542	INTEGER i
543
544	DO i=1,klon
545
546	IF (phase .EQ. 1) THEN
547	delta=0.
548	ELSEIF (phase .EQ. 2) THEN
549	delta=1.
550	ELSE
551	delta=MAX(0.,SIGN(1.,tref-temp(i)))
552	ENDIF
553
554	IF (flagth) THEN
555	cvm5=R5LES(1.-delta) + R5IESdelta
556	ELSE
557	cvm5 = R5LESRLVTT(1.-delta) + R5IESRLSTTdelta
558	cvm5 = cvm5 /RCPD/(1.0+RVTMP2*(qtot(i)))
559	ENDIF
560
561	qs(i)= R2ES*FOEEW(temp(i),delta)/pressure(i)
562	qs(i)=MIN(0.5,qs(i))
563	cor=1./(1.-RETV*qs(i))
564	qs(i)=qs(i)*cor
565	dqs(i)= FOEDE(temp(i),delta,cvm5,qs(i),cor)
566
567	END DO
568
569	END SUBROUTINE CALC_QSAT_ECMWF
570	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
571
572
573	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
574	SUBROUTINE CALC_GAMMASAT(klon,temp,qtot,pressure,gammasat,dgammasatdt)
575
576	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
577	! programme that calculates the gammasat parameter that determines the
578	! homogeneous condensation thresholds for cold (<0oC) clouds
579	! condensation at q>gammasat*qsat
580	! Etienne Vignon, March 2021
581	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
582
583	use lmdz_lscp_ini, only: iflag_gammasat, t_glace_min, RTT
584
585	IMPLICIT NONE
586
587
588	INTEGER, INTENT(IN) :: klon ! number of horizontal grid points
589	REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature in K
590	REAL, INTENT(IN), DIMENSION(klon) :: qtot ! total specific water in kg/kg
591
592	REAL, INTENT(IN), DIMENSION(klon) :: pressure ! pressure in Pa
593
594	REAL, INTENT(OUT), DIMENSION(klon) :: gammasat ! coefficient to multiply qsat with to calculate saturation
595	REAL, INTENT(OUT), DIMENSION(klon) :: dgammasatdt ! derivative of gammasat wrt temperature
596
597	REAL, DIMENSION(klon) :: qsi,qsl,dqsl,dqsi
598	REAL fcirrus, fac
599	REAL, PARAMETER :: acirrus=2.349
600	REAL, PARAMETER :: bcirrus=259.0
601
602	INTEGER i
603
604	CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,1,.false.,qsl,dqsl)
605	CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,2,.false.,qsi,dqsi)
606
607	DO i=1,klon
608
609	IF (temp(i) .GE. RTT) THEN
610	! warm clouds: condensation at saturation wrt liquid
611	gammasat(i)=1.
612	dgammasatdt(i)=0.
613
614	ELSEIF ((temp(i) .LT. RTT) .AND. (temp(i) .GT. t_glace_min)) THEN
615
616	IF (iflag_gammasat .GE. 2) THEN
617	gammasat(i)=qsl(i)/qsi(i)
618	dgammasatdt(i)=(dqsl(i)qsi(i)-dqsi(i)qsl(i))/qsi(i)/qsi(i)
619	ELSE
620	gammasat(i)=1.
621	dgammasatdt(i)=0.
622	ENDIF
623
624	ELSE
625
626	IF (iflag_gammasat .GE.1) THEN
627	! homogeneous freezing of aerosols, according to
628	! Koop, 2000 and Karcher 2008, QJRMS
629	! 'Cirrus regime'
630	fcirrus=acirrus-temp(i)/bcirrus
631	IF (fcirrus .GT. qsl(i)/qsi(i)) THEN
632	gammasat(i)=qsl(i)/qsi(i)
633	dgammasatdt(i)=(dqsl(i)qsi(i)-dqsi(i)qsl(i))/qsi(i)/qsi(i)
634	ELSE
635	gammasat(i)=fcirrus
636	dgammasatdt(i)=-1.0/bcirrus
637	ENDIF
638
639	ELSE
640
641	gammasat(i)=1.
642	dgammasatdt(i)=0.
643
644	ENDIF
645
646	ENDIF
647
648	END DO
649
650
651	END SUBROUTINE CALC_GAMMASAT
652	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
653
654
655	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
656	SUBROUTINE DISTANCE_TO_CLOUD_TOP(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D,temp_cltop)
657	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
658
659	USE lmdz_lscp_ini, ONLY : rd,rg,tresh_cl
660
661	IMPLICIT NONE
662
663	INTEGER, INTENT(IN) :: klon,klev !number of horizontal and vertical grid points
664	INTEGER, INTENT(IN) :: k ! vertical index
665	REAL, INTENT(IN), DIMENSION(klon,klev) :: temp ! temperature in K
666	REAL, INTENT(IN), DIMENSION(klon,klev) :: pplay ! pressure middle layer in Pa
667	REAL, INTENT(IN), DIMENSION(klon,klev+1) :: paprs ! pressure interfaces in Pa
668	REAL, INTENT(IN), DIMENSION(klon,klev) :: rneb ! cloud fraction
669
670	REAL, INTENT(OUT), DIMENSION(klon) :: distcltop1D ! distance from cloud top
671	REAL, INTENT(OUT), DIMENSION(klon) :: temp_cltop ! temperature of cloud top
672
673	REAL dzlay(klon,klev)
674	REAL zlay(klon,klev)
675	REAL dzinterf
676	INTEGER i,k_top, kvert
677	LOGICAL bool_cl
678
679
680	DO i=1,klon
681	! Initialization height middle of first layer
682	dzlay(i,1) = Rd * temp(i,1) / rg * log(paprs(i,1)/paprs(i,2))
683	zlay(i,1) = dzlay(i,1)/2
684
685	DO kvert=2,klev
686	IF (kvert.EQ.klev) THEN
687	dzlay(i,kvert) = 2(rd temp(i,kvert) / rg * log(paprs(i,kvert)/pplay(i,kvert)))
688	ELSE
689	dzlay(i,kvert) = rd * temp(i,kvert) / rg * log(paprs(i,kvert)/paprs(i,kvert+1))
690	ENDIF
691	dzinterf = rd * temp(i,kvert) / rg * log(pplay(i,kvert-1)/pplay(i,kvert))
692	zlay(i,kvert) = zlay(i,kvert-1) + dzinterf
693	ENDDO
694	ENDDO
695
696	DO i=1,klon
697	k_top = k
698	IF (rneb(i,k) .LE. tresh_cl) THEN
699	bool_cl = .FALSE.
700	ELSE
701	bool_cl = .TRUE.
702	ENDIF
703
704	DO WHILE ((bool_cl) .AND. (k_top .LE. klev))
705	! find cloud top
706	IF (rneb(i,k_top) .GT. tresh_cl) THEN
707	k_top = k_top + 1
708	ELSE
709	bool_cl = .FALSE.
710	k_top = k_top - 1
711	ENDIF
712	ENDDO
713	k_top=min(k_top,klev)
714
715	!dist to top is dist between current layer and layer of cloud top (from middle to middle) + dist middle to
716	!interf for layer of cloud top
717	distcltop1D(i) = zlay(i,k_top) - zlay(i,k) + dzlay(i,k_top)/2
718	temp_cltop(i) = temp(i,k_top)
719	ENDDO ! klon
720
721	END SUBROUTINE DISTANCE_TO_CLOUD_TOP
722	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
723
724	END MODULE lmdz_lscp_tools
725
726

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

Download in other formats:

Original Format