Context Navigation

lmdz_lscp_tools.F90 @ 5201

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

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

Download in other formats:

Original Format