Context Navigation

lmdz_lscp_tools.f90 @ 5279

Last change on this file since 5279 was 5274, checked in by abarral, 9 months ago
Replace yomcst.h by existing module
File size: 36.2 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	USE yomcst_mod_h, ONLY: RPI, RCLUM, RHPLA, RKBOL, RNAVO &
522	, RDAY, REA, REPSM, RSIYEA, RSIDAY, ROMEGA &
523	, R_ecc, R_peri, R_incl &
524	, RA, RG, R1SA &
525	, RSIGMA &
526	, R, RMD, RMV, RD, RV, RCPD &
527	, RMO3, RMCO2, RMC, RMCH4, RMN2O, RMCFC11, RMCFC12 &
528	, RCPV, RCVD, RCVV, RKAPPA, RETV, eps_w &
529	, RCW, RCS &
530	, RLVTT, RLSTT, RLMLT, RTT, RATM &
531	, RESTT, RALPW, RBETW, RGAMW, RALPS, RBETS, RGAMS &
532	, RALPD, RBETD, RGAMD
533	IMPLICIT NONE
534
535
536	include "YOETHF.h"
537	include "FCTTRE.h"
538
539	INTEGER, INTENT(IN) :: klon ! number of horizontal grid points
540	REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature in K
541	REAL, INTENT(IN), DIMENSION(klon) :: qtot ! total specific water in kg/kg
542	REAL, INTENT(IN), DIMENSION(klon) :: pressure ! pressure in Pa
543	REAL, INTENT(IN) :: tref ! reference temperature in K
544	LOGICAL, INTENT(IN) :: flagth ! flag for qsat calculation for thermals
545	INTEGER, INTENT(IN) :: phase
546	! phase: 0=depend on temperature sign (temp>tref -> liquid, temp<tref, solid)
547	! 1=liquid
548	! 2=solid
549
550	REAL, INTENT(OUT), DIMENSION(klon) :: qs ! saturation specific humidity [kg/kg]
551	REAL, INTENT(OUT), DIMENSION(klon) :: dqs ! derivation of saturation specific humidity wrt T
552
553	REAL delta, cor, cvm5
554	INTEGER i
555
556	DO i=1,klon
557
558	IF (phase .EQ. 1) THEN
559	delta=0.
560	ELSEIF (phase .EQ. 2) THEN
561	delta=1.
562	ELSE
563	delta=MAX(0.,SIGN(1.,tref-temp(i)))
564	ENDIF
565
566	IF (flagth) THEN
567	cvm5=R5LES(1.-delta) + R5IESdelta
568	ELSE
569	cvm5 = R5LESRLVTT(1.-delta) + R5IESRLSTTdelta
570	cvm5 = cvm5 /RCPD/(1.0+RVTMP2*(qtot(i)))
571	ENDIF
572
573	qs(i)= R2ES*FOEEW(temp(i),delta)/pressure(i)
574	qs(i)=MIN(0.5,qs(i))
575	cor=1./(1.-RETV*qs(i))
576	qs(i)=qs(i)*cor
577	dqs(i)= FOEDE(temp(i),delta,cvm5,qs(i),cor)
578
579	END DO
580
581	END SUBROUTINE CALC_QSAT_ECMWF
582	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
583
584
585	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
586	SUBROUTINE CALC_GAMMASAT(klon,temp,qtot,pressure,gammasat,dgammasatdt)
587
588	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
589	! programme that calculates the gammasat parameter that determines the
590	! homogeneous condensation thresholds for cold (<0oC) clouds
591	! condensation at q>gammasat*qsat
592	! Etienne Vignon, March 2021
593	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
594
595	use lmdz_lscp_ini, only: iflag_gammasat, temp_nowater, RTT
596	use lmdz_lscp_ini, only: a_homofreez, b_homofreez, delta_hetfreez
597
598	IMPLICIT NONE
599
600
601	INTEGER, INTENT(IN) :: klon ! number of horizontal grid points
602	REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature in K
603	REAL, INTENT(IN), DIMENSION(klon) :: qtot ! total specific water in kg/kg
604
605	REAL, INTENT(IN), DIMENSION(klon) :: pressure ! pressure in Pa
606
607	REAL, INTENT(OUT), DIMENSION(klon) :: gammasat ! coefficient to multiply qsat with to calculate saturation
608	REAL, INTENT(OUT), DIMENSION(klon) :: dgammasatdt ! derivative of gammasat wrt temperature
609
610	REAL, DIMENSION(klon) :: qsi,qsl,dqsl,dqsi
611	REAL f_homofreez, fac
612
613	INTEGER i
614
615	CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,1,.false.,qsl,dqsl)
616	CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,2,.false.,qsi,dqsi)
617
618	DO i = 1, klon
619
620	IF ( temp(i) .GE. RTT ) THEN
621	! warm clouds: condensation at saturation wrt liquid
622	gammasat(i) = 1.
623	dgammasatdt(i) = 0.
624
625	ELSE
626	! cold clouds: qsi > qsl
627
628	! homogeneous freezing of aerosols, according to
629	! Koop, 2000 and Ren and MacKenzie, 2005 (QJRMS)
630	! 'Cirrus regime'
631	! if f_homofreez > qsl / qsi, liquid nucleation
632	! if f_homofreez < qsl / qsi, homogeneous freezing of aerosols
633	! Note: f_homofreez = qsl / qsi for temp ~= -38degC
634	f_homofreez = a_homofreez - temp(i) / b_homofreez
635
636	IF ( iflag_gammasat .GE. 3 ) THEN
637	! condensation at homogeneous freezing threshold for temp < -38 degC
638	! condensation at liquid saturation for temp > -38 degC
639	IF ( f_homofreez .LE. qsl(i) / qsi(i) ) THEN
640	gammasat(i) = f_homofreez
641	dgammasatdt(i) = - 1. / b_homofreez
642	ELSE
643	gammasat(i) = qsl(i) / qsi(i)
644	dgammasatdt(i) = ( dqsl(i) * qsi(i) - dqsi(i) * qsl(i) ) / qsi(i) / qsi(i)
645	ENDIF
646
647	ELSEIF ( iflag_gammasat .EQ. 2 ) THEN
648	! condensation at homogeneous freezing threshold for temp < -38 degC
649	! condensation at a threshold linearly decreasing between homogeneous
650	! freezing and ice saturation for -38 degC < temp < temp_nowater
651	! condensation at ice saturation for temp > temp_nowater
652	! If temp_nowater = 235.15 K, this is equivalent to iflag_gammasat = 1
653	IF ( f_homofreez .LE. qsl(i) / qsi(i) ) THEN
654	gammasat(i) = f_homofreez
655	dgammasatdt(i) = - 1. / b_homofreez
656	ELSEIF ( temp(i) .LE. temp_nowater ) THEN
657	! Here, we assume that f_homofreez = qsl / qsi for temp = -38 degC = 235.15 K
658	dgammasatdt(i) = ( a_homofreez - 235.15 / b_homofreez - 1. ) &
659	/ ( 235.15 - temp_nowater )
660	gammasat(i) = dgammasatdt(i) * ( temp(i) - temp_nowater ) + 1.
661	ELSE
662	gammasat(i) = 1.
663	dgammasatdt(i) = 0.
664	ENDIF
665
666	ELSEIF ( iflag_gammasat .EQ. 1 ) THEN
667	! condensation at homogeneous freezing threshold for temp < -38 degC
668	! condensation at ice saturation for temp > -38 degC
669	IF ( f_homofreez .LE. qsl(i) / qsi(i) ) THEN
670	gammasat(i) = f_homofreez
671	dgammasatdt(i) = - 1. / b_homofreez
672	ELSE
673	gammasat(i) = 1.
674	dgammasatdt(i) = 0.
675	ENDIF
676
677	ELSE
678	! condensation at ice saturation for temp < -38 degC
679	! condensation at ice saturation for temp > -38 degC
680	gammasat(i) = 1.
681	dgammasatdt(i) = 0.
682
683	ENDIF
684
685	! Note that the delta_hetfreez parameter allows to linearly decrease the
686	! condensation threshold between the calculated threshold and the ice saturation
687	! for delta_hetfreez = 1, the threshold is the calculated condensation threshold
688	! for delta_hetfreez = 0, the threshold is the ice saturation
689	gammasat(i) = ( 1. - delta_hetfreez ) + delta_hetfreez * gammasat(i)
690	dgammasatdt(i) = delta_hetfreez * dgammasatdt(i)
691
692	ENDIF
693
694	END DO
695
696
697	END SUBROUTINE CALC_GAMMASAT
698	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
699
700
701	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
702	SUBROUTINE DISTANCE_TO_CLOUD_TOP(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D,temp_cltop)
703	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
704
705	USE lmdz_lscp_ini, ONLY : rd,rg,tresh_cl
706
707	IMPLICIT NONE
708
709	INTEGER, INTENT(IN) :: klon,klev !number of horizontal and vertical grid points
710	INTEGER, INTENT(IN) :: k ! vertical index
711	REAL, INTENT(IN), DIMENSION(klon,klev) :: temp ! temperature in K
712	REAL, INTENT(IN), DIMENSION(klon,klev) :: pplay ! pressure middle layer in Pa
713	REAL, INTENT(IN), DIMENSION(klon,klev+1) :: paprs ! pressure interfaces in Pa
714	REAL, INTENT(IN), DIMENSION(klon,klev) :: rneb ! cloud fraction
715
716	REAL, INTENT(OUT), DIMENSION(klon) :: distcltop1D ! distance from cloud top
717	REAL, INTENT(OUT), DIMENSION(klon) :: temp_cltop ! temperature of cloud top
718
719	REAL dzlay(klon,klev)
720	REAL zlay(klon,klev)
721	REAL dzinterf
722	INTEGER i,k_top, kvert
723	LOGICAL bool_cl
724
725
726	DO i=1,klon
727	! Initialization height middle of first layer
728	dzlay(i,1) = Rd * temp(i,1) / rg * log(paprs(i,1)/paprs(i,2))
729	zlay(i,1) = dzlay(i,1)/2
730
731	DO kvert=2,klev
732	IF (kvert.EQ.klev) THEN
733	dzlay(i,kvert) = 2(rd temp(i,kvert) / rg * log(paprs(i,kvert)/pplay(i,kvert)))
734	ELSE
735	dzlay(i,kvert) = rd * temp(i,kvert) / rg * log(paprs(i,kvert)/paprs(i,kvert+1))
736	ENDIF
737	dzinterf = rd * temp(i,kvert) / rg * log(pplay(i,kvert-1)/pplay(i,kvert))
738	zlay(i,kvert) = zlay(i,kvert-1) + dzinterf
739	ENDDO
740	ENDDO
741
742	DO i=1,klon
743	k_top = k
744	IF (rneb(i,k) .LE. tresh_cl) THEN
745	bool_cl = .FALSE.
746	ELSE
747	bool_cl = .TRUE.
748	ENDIF
749
750	DO WHILE ((bool_cl) .AND. (k_top .LE. klev))
751	! find cloud top
752	IF (rneb(i,k_top) .GT. tresh_cl) THEN
753	k_top = k_top + 1
754	ELSE
755	bool_cl = .FALSE.
756	k_top = k_top - 1
757	ENDIF
758	ENDDO
759	k_top=min(k_top,klev)
760
761	!dist to top is dist between current layer and layer of cloud top (from middle to middle) + dist middle to
762	!interf for layer of cloud top
763	distcltop1D(i) = zlay(i,k_top) - zlay(i,k) + dzlay(i,k_top)/2
764	temp_cltop(i) = temp(i,k_top)
765	ENDDO ! klon
766
767	END SUBROUTINE DISTANCE_TO_CLOUD_TOP
768	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
769
770	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
771	FUNCTION GAMMAINC ( p, x )
772
773	!*****************************************************************************80
774	!
775	!! GAMMAINC computes the regularized lower incomplete Gamma Integral
776	!
777	! Modified:
778	!
779	! 20 January 2008
780	!
781	! Author:
782	!
783	! Original FORTRAN77 version by B Shea.
784	! FORTRAN90 version by John Burkardt.
785	!
786	! Reference:
787	!
788	! B Shea,
789	! Algorithm AS 239:
790	! Chi-squared and Incomplete Gamma Integral,
791	! Applied Statistics,
792	! Volume 37, Number 3, 1988, pages 466-473.
793	!
794	! Parameters:
795	!
796	! Input, real X, P, the parameters of the incomplete
797	! gamma ratio. 0 <= X, and 0 < P.
798	!
799	! Output, real GAMMAINC, the value of the incomplete
800	! Gamma integral.
801	!
802	IMPLICIT NONE
803
804	REAL A
805	REAL AN
806	REAL ARG
807	REAL B
808	REAL C
809	REAL, PARAMETER :: ELIMIT = - 88.0E+00
810	REAL GAMMAINC
811	REAL, PARAMETER :: OFLO = 1.0E+37
812	REAL P
813	REAL, PARAMETER :: PLIMIT = 1000.0E+00
814	REAL PN1
815	REAL PN2
816	REAL PN3
817	REAL PN4
818	REAL PN5
819	REAL PN6
820	REAL RN
821	REAL, PARAMETER :: TOL = 1.0E-14
822	REAL X
823	REAL, PARAMETER :: XBIG = 1.0E+08
824
825	GAMMAINC = 0.0E+00
826
827	IF ( X == 0.0E+00 ) THEN
828	GAMMAINC = 0.0E+00
829	RETURN
830	END IF
831	!
832	! IF P IS LARGE, USE A NORMAL APPROXIMATION.
833	!
834	IF ( PLIMIT < P ) THEN
835
836	PN1 = 3.0E+00 * SQRT ( P ) * ( ( X / P )**( 1.0E+00 / 3.0E+00 ) &
837	+ 1.0E+00 / ( 9.0E+00 * P ) - 1.0E+00 )
838
839	GAMMAINC = 0.5E+00 * ( 1. + ERF ( PN1 ) )
840	RETURN
841
842	END IF
843	!
844	! IF X IS LARGE SET GAMMAD = 1.
845	!
846	IF ( XBIG < X ) THEN
847	GAMMAINC = 1.0E+00
848	RETURN
849	END IF
850	!
851	! USE PEARSON'S SERIES EXPANSION.
852	! (NOTE THAT P IS NOT LARGE ENOUGH TO FORCE OVERFLOW IN ALOGAM).
853	!
854	IF ( X <= 1.0E+00 .OR. X < P ) THEN
855
856	ARG = P * LOG ( X ) - X - LOG_GAMMA ( P + 1.0E+00 )
857	C = 1.0E+00
858	GAMMAINC = 1.0E+00
859	A = P
860
861	DO
862
863	A = A + 1.0E+00
864	C = C * X / A
865	GAMMAINC = GAMMAINC + C
866
867	IF ( C <= TOL ) THEN
868	EXIT
869	END IF
870
871	END DO
872
873	ARG = ARG + LOG ( GAMMAINC )
874
875	IF ( ELIMIT <= ARG ) THEN
876	GAMMAINC = EXP ( ARG )
877	ELSE
878	GAMMAINC = 0.0E+00
879	END IF
880	!
881	! USE A CONTINUED FRACTION EXPANSION.
882	!
883	ELSE
884
885	ARG = P * LOG ( X ) - X - LOG_GAMMA ( P )
886	A = 1.0E+00 - P
887	B = A + X + 1.0E+00
888	C = 0.0E+00
889	PN1 = 1.0E+00
890	PN2 = X
891	PN3 = X + 1.0E+00
892	PN4 = X * B
893	GAMMAINC = PN3 / PN4
894
895	DO
896
897	A = A + 1.0E+00
898	B = B + 2.0E+00
899	C = C + 1.0E+00
900	AN = A * C
901	PN5 = B * PN3 - AN * PN1
902	PN6 = B * PN4 - AN * PN2
903
904	IF ( PN6 /= 0.0E+00 ) THEN
905
906	RN = PN5 / PN6
907
908	IF ( ABS ( GAMMAINC - RN ) <= MIN ( TOL, TOL * RN ) ) THEN
909	EXIT
910	END IF
911
912	GAMMAINC = RN
913
914	END IF
915
916	PN1 = PN3
917	PN2 = PN4
918	PN3 = PN5
919	PN4 = PN6
920	!
921	! RE-SCALE TERMS IN CONTINUED FRACTION IF TERMS ARE LARGE.
922	!
923	IF ( OFLO <= ABS ( PN5 ) ) THEN
924	PN1 = PN1 / OFLO
925	PN2 = PN2 / OFLO
926	PN3 = PN3 / OFLO
927	PN4 = PN4 / OFLO
928	END IF
929
930	END DO
931
932	ARG = ARG + LOG ( GAMMAINC )
933
934	IF ( ELIMIT <= ARG ) THEN
935	GAMMAINC = 1.0E+00 - EXP ( ARG )
936	ELSE
937	GAMMAINC = 1.0E+00
938	END IF
939
940	END IF
941
942	RETURN
943	END FUNCTION GAMMAINC
944	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
945
946	END MODULE lmdz_lscp_tools
947
948

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

Download in other formats:

Original Format