Context Navigation

lmdz_lscp_tools.f90 @ 5756

Last change on this file since 5756 was 5727, checked in by evignon, 11 days ago
legere modification de lscp pour enlever un warning obsolete sur la parameterisation de phase
File size: 38.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
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 or 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
243	SUBROUTINE ICEFRAC_LSCP_TURB(klon, dtime, pticefracturb, temp, pplay, paprsdn, paprsup, wvel, qice_ini, snowcld, qtot_incl, cldfra, tke, &
244	tke_dissip, sursat_e, invtau_e, qliq, qvap_cld, qice, icefrac, dicefracdT, cldfraliq, sigma2_icefracturb, mean_icefracturb)
245	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
246	! Compute the liquid, ice and vapour content (+ice fraction) based
247	! on turbulence (see Fields 2014, Furtado 2016, Raillard 2025)
248	! L.Raillard (23/09/24)
249	! E.Vignon (03/2025) : additional elements for treatment of convective
250	! boundary layer clouds
251	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
252
253
254	USE lmdz_lscp_ini, ONLY : prt_level, lunout
255	USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RV, RPI
256	USE lmdz_lscp_ini, ONLY : seuil_neb, temp_nowater
257	USE lmdz_lscp_ini, ONLY : naero5, gamma_snwretro, gamma_taud, capa_crystal, rho_ice
258	USE lmdz_lscp_ini, ONLY : eps, snow_fallspeed
259
260	IMPLICIT NONE
261
262	INTEGER, INTENT(IN) :: klon !--number of horizontal grid points
263	REAL, INTENT(IN) :: dtime !--time step [s]
264	LOGICAL, INTENT(IN), DIMENSION(klon) :: pticefracturb !--grid points concerned by this routine
265	REAL, INTENT(IN), DIMENSION(klon) :: temp !--temperature
266	REAL, INTENT(IN), DIMENSION(klon) :: pplay !--pressure in the middle of the layer [Pa]
267	REAL, INTENT(IN), DIMENSION(klon) :: paprsdn !--pressure at the bottom interface of the layer [Pa]
268	REAL, INTENT(IN), DIMENSION(klon) :: paprsup !--pressure at the top interface of the layer [Pa]
269	REAL, INTENT(IN), DIMENSION(klon) :: wvel !--vertical velocity [m/s]
270	REAL, INTENT(IN), DIMENSION(klon) :: qtot_incl !--specific total cloud water in-cloud content [kg/kg]
271	REAL, INTENT(IN), DIMENSION(klon) :: cldfra !--cloud fraction in gridbox [-]
272	REAL, INTENT(IN), DIMENSION(klon) :: tke !--turbulent kinetic energy [m2/s2]
273	REAL, INTENT(IN), DIMENSION(klon) :: tke_dissip !--TKE dissipation [m2/s3]
274
275	REAL, INTENT(IN), DIMENSION(klon) :: qice_ini !--initial specific ice content gridbox-mean [kg/kg]
276	REAL, INTENT(IN), DIMENSION(klon) :: snowcld !--in-cloud snowfall flux [kg/m2/s]
277	REAL, INTENT(IN), DIMENSION(klon) :: sursat_e !--environment supersaturation [-]
278	REAL, INTENT(IN), DIMENSION(klon) :: invtau_e !--inverse time-scale of mixing with environment [s-1]
279	REAL, INTENT(OUT), DIMENSION(klon) :: qliq !--specific liquid content gridbox-mean [kg/kg]
280	REAL, INTENT(OUT), DIMENSION(klon) :: qvap_cld !--specific cloud vapor content, gridbox-mean [kg/kg]
281	REAL, INTENT(OUT), DIMENSION(klon) :: qice !--specific ice content gridbox-mean [kg/kg]
282
283	REAL, INTENT(INOUT), DIMENSION(klon) :: icefrac !--fraction of ice in condensed water [-]
284	REAL, INTENT(INOUT), DIMENSION(klon) :: dicefracdT
285
286	REAL, INTENT(OUT), DIMENSION(klon) :: cldfraliq !--fraction of cldfra where liquid [-]
287	REAL, INTENT(OUT), DIMENSION(klon) :: sigma2_icefracturb !--Sigma2 of the ice supersaturation PDF [-]
288	REAL, INTENT(OUT), DIMENSION(klon) :: mean_icefracturb !--Mean of the ice supersaturation PDF [-]
289
290	REAL, DIMENSION(klon) :: qzero, qsatl, dqsatl, qsati, dqsati !--specific humidity saturation values
291	INTEGER :: i
292
293	REAL :: qvap_incl, qice_incl, qliq_incl, qiceini_incl !--In-cloud specific quantities [kg/kg]
294	REAL :: water_vapor_diff !--Water-vapour diffusion coeff in air (f(T,P)) [m2/s]
295	REAL :: air_thermal_conduct !--Thermal conductivity of air (f(T)) [J/m/K/s]
296	REAL :: C0 !--Lagrangian structure function [-]
297	REAL :: tau_dissipturb
298	REAL :: invtau_phaserelax
299	REAL :: sigma2_pdf
300	REAL :: ai, bi, B0
301	REAL :: sursat_iceliq
302	REAL :: sursat_equ
303	REAL :: liqfra_max
304	REAL :: sursat_iceext
305	REAL :: nb_crystals !--number concentration of ice crystals [#/m3]
306	REAL :: moment1_PSD !--1st moment of ice PSD
307	REAL :: N0_PSD, lambda_PSD !--parameters of the exponential PSD
308
309	REAL :: cldfra1D
310	REAL :: rho_air
311	REAL :: psati !--saturation vapor pressure wrt ice [Pa]
312
313
314	REAL :: tempvig1, tempvig2
315
316	tempvig1 = -21.06 + RTT
317	tempvig2 = -30.35 + RTT
318	C0 = 10. !--value assumed in Field2014
319	sursat_iceext = -0.1
320	qzero(:) = 0.
321	cldfraliq(:) = 0.
322	qliq(:) = 0.
323	qice(:) = 0.
324	qvap_cld(:) = 0.
325	sigma2_icefracturb(:) = 0.
326	mean_icefracturb(:) = 0.
327
328	!--wrt liquid
329	CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,1,.false.,qsatl(:),dqsatl(:))
330	!--wrt ice
331	CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,2,.false.,qsati(:),dqsati(:))
332
333
334	DO i=1,klon
335	rho_air = pplay(i) / temp(i) / RD
336	! because cldfra is intent in, but can be locally modified due to test
337	cldfra1D = cldfra(i)
338	! activate param for concerned grid points and for cloudy conditions
339	IF ((pticefracturb(i)) .AND. (cldfra(i) .GT. 0.)) THEN
340	IF (cldfra(i) .GE. 1.0) THEN
341	cldfra1D = 1.0
342	END IF
343
344	! T>0°C, no ice allowed
345	IF ( temp(i) .GE. RTT ) THEN
346	qvap_cld(i) = qsatl(i) * cldfra1D
347	qliq(i) = MAX(0.0,qtot_incl(i)-qsatl(i)) * cldfra1D
348	qice(i) = 0.
349	cldfraliq(i) = 1.
350	icefrac(i) = 0.
351	dicefracdT(i) = 0.
352
353	! T<-38°C, no liquid allowed
354	ELSE IF ( temp(i) .LE. temp_nowater) THEN
355	qvap_cld(i) = qsati(i) * cldfra1D
356	qliq(i) = 0.
357	qice(i) = MAX(0.0,qtot_incl(i)-qsati(i)) * cldfra1D
358	cldfraliq(i) = 0.
359	icefrac(i) = 1.
360	dicefracdT(i) = 0.
361
362
363	!---------------------------------------------------------
364	!-- MIXED PHASE TEMPERATURE REGIME
365	!---------------------------------------------------------
366	!--In the mixed phase regime (-38°C< T <0°C) we distinguish
367	!--3 possible subcases.
368	!--1. No pre-existing ice
369	!--2A. Pre-existing ice and no turbulence
370	!--2B. Pre-existing ice and turbulence
371	ELSE
372
373	! gamma_snwretro controls the contribution of snowflakes to the negative feedback
374	! note that for reasons related to inetarctions with the condensation iteration in lscp_main
375	! we consider here the mean snowflake concentration in the mesh (not the in-cloud concentration)
376	! when poprecip is active, it will be worth testing considering the incloud fraction, dividing
377	! by znebprecipcld
378	! qiceini_incl = qice_ini(i) / cldfra1D + &
379	! gamma_snwretro * snowcld(i) * RG * dtime / ( paprsdn(i) - paprsup(i) )
380	! assuming constant snowfall velocity
381	qiceini_incl = qice_ini(i) / cldfra1D + gamma_snwretro * snowcld(i) / pplay(i) * RD * temp(i) / snow_fallspeed
382
383	!--1. No preexisting ice and no mixing with environment: if vertical motion, fully liquid
384	!--cloud else fully iced cloud
385	IF ( (qiceini_incl .LT. eps) .AND. (invtau_e(i) .LT. eps) ) THEN
386	IF ( (wvel(i) .GT. eps) .OR. (tke(i) .GT. eps) ) THEN
387	qvap_cld(i) = qsatl(i) * cldfra1D
388	qliq(i) = MAX(0.,qtot_incl(i)-qsatl(i)) * cldfra1D
389	qice(i) = 0.
390	cldfraliq(i) = 1.
391	icefrac(i) = 0.
392	dicefracdT(i) = 0.
393	ELSE
394	qvap_cld(i) = qsati(i) * cldfra1D
395	qliq(i) = 0.
396	qice(i) = MAX(0.,qtot_incl(i)-qsati(i)) * cldfra1D
397	cldfraliq(i) = 0.
398	icefrac(i) = 1.
399	dicefracdT(i) = 0.
400	ENDIF
401
402
403	!--2. Pre-existing ice and/or mixing with environment:computation of ice properties for
404	!--feedback
405	ELSE
406
407	sursat_iceliq = qsatl(i)/qsati(i) - 1.
408	psati = qsati(i) * pplay(i) / (RD/RV)
409
410	!--We assume an exponential ice PSD whose parameters
411	!--are computed following Morrison&Gettelman 2008
412	!--Ice number density is assumed equals to INP density
413	!--which is for naero5>0 a function of temperature (DeMott 2010)
414	!--bi and B0 are microphysical function characterizing
415	!--vapor/ice interactions
416	!--tau_phase_relax is the typical time of vapor deposition
417	!--onto ice crystals
418
419	!--For naero5<=0 INP density is derived from the empirical fit
420	!--from MARCUS campaign from Vignon 2021
421	!--/!\ Note that option is very specific and should be use for
422	!--the Southern Ocean and the Antarctic
423
424	IF (naero5 .LE. 0) THEN
425	IF ( temp(i) .GT. tempvig1 ) THEN
426	nb_crystals = 1.e3 * 10*(-0.14(temp(i)-tempvig1) - 2.88)
427	ELSE IF ( temp(i) .GT. tempvig2 ) THEN
428	nb_crystals = 1.e3 * 10*(-0.31(temp(i)-tempvig1) - 2.88)
429	ELSE
430	nb_crystals = 1.e3 * 10**(0.)
431	ENDIF
432	ELSE
433	nb_crystals = 1.e3 * 5.94e-5 * ( RTT - temp(i) )*3.33 naero5*(0.0264(RTT-temp(i))+0.0033)
434	ENDIF
435	lambda_PSD = ( (RPIrho_icenb_crystals) / (rho_air * MAX(qiceini_incl , eps) ) ) ** (1./3.)
436	N0_PSD = nb_crystals * lambda_PSD
437	moment1_PSD = N0_PSD/lambda_PSD**2
438
439	!--Formulae for air thermal conductivity and water vapor diffusivity
440	!--comes respectively from Beard and Pruppacher (1971)
441	!--and Hall and Pruppacher (1976)
442
443	air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184
444	water_vapor_diff = 2.111e-5 ( temp(i) / RTT )*1.94 ( 101325 / pplay(i) )
445
446	bi = 1./((qsati(i)+qsatl(i))/2.) + RLSTT2 / RCPD / RV / temp(i)2
447	B0 = 4. * RPI * capa_crystal * 1. / ( RLSTT2 / air_thermal_conduct / RV / temp(i)2 &
448	+ RV * temp(i) / psati / water_vapor_diff )
449	invtau_phaserelax = bi * B0 * moment1_PSD
450
451	ai = RG / RD / temp(i) * ( RD * RLSTT / RCPD / RV / temp(i) - 1. )
452	sursat_equ = (ai * wvel(i) + sursat_e(i)*invtau_e(i)) / (invtau_phaserelax + invtau_e(i))
453	! as sursaturation is by definition lower than -1 and
454	! because local supersaturation > 1 are never found in the atmosphere
455
456	!--2A. No TKE : stationnary binary solution depending on vertical velocity and mixing with env.
457	! If Sequ > Siw liquid cloud, else ice cloud
458	IF ( tke_dissip(i) .LE. eps ) THEN
459	sigma2_icefracturb(i)= 0.
460	mean_icefracturb(i) = sursat_equ
461	IF (sursat_equ .GT. sursat_iceliq) THEN
462	qvap_cld(i) = qsatl(i) * cldfra1D
463	qliq(i) = MAX(0.,qtot_incl(i)-qsatl(i)) * cldfra1D
464	qice(i) = 0.
465	cldfraliq(i) = 1.
466	icefrac(i) = 0.
467	dicefracdT(i) = 0.
468	ELSE
469	qvap_cld(i) = qsati(i) * cldfra1D
470	qliq(i) = 0.
471	qice(i) = MAX(0.,qtot_incl(i)-qsati(i)) * cldfra1D
472	cldfraliq(i) = 0.
473	icefrac(i) = 1.
474	dicefracdT(i) = 0.
475	ENDIF
476
477	!--2B. TKE and ice : ice supersaturation PDF
478	!--we compute the cloud liquid properties with a Gaussian PDF
479	!--of ice supersaturation F(Si) (Field2014, Furtado2016).
480	!--Parameters of the PDF are function of turbulence and
481	!--microphysics/existing ice.
482	ELSE
483
484	!--Tau_dissipturb is the time needed for turbulence to decay
485	!--due to viscosity
486	tau_dissipturb = gamma_taud * 2. * 2./3. * tke(i) / tke_dissip(i) / C0
487
488	!--------------------- PDF COMPUTATIONS ---------------------
489	!--Formulae for sigma2_pdf (variance), mean of PDF in Raillard2025
490	!--cloud liquid fraction and in-cloud liquid content are given
491	!--by integrating resp. F(Si) and Si*F(Si)
492	!--Liquid is limited by the available water vapor trough a
493	!--maximal liquid fraction
494	!--qice_ini(i) / cldfra1D = qiceincld without precip
495
496	liqfra_max = MAX(0., (MIN (1.,( qtot_incl(i) - (qice_ini(i) / cldfra1D) - qsati(i) * (1 + sursat_iceext ) ) / ( qsatl(i) - qsati(i) ) ) ) )
497	sigma2_pdf = 1./2. * ( ai*2 ) 2./3. * tke(i) * tau_dissipturb / (invtau_phaserelax + invtau_e(i))
498	! sursat ranges between -1 and 1, so we prevent sigma2 so exceed 1
499	cldfraliq(i) = 0.5 * (1. - erf( ( sursat_iceliq - sursat_equ) / (SQRT(2.* sigma2_pdf) ) ) )
500	IF (cldfraliq(i) .GT. liqfra_max) THEN
501	cldfraliq(i) = liqfra_max
502	ENDIF
503
504	qliq_incl = qsati(i) * SQRT(sigma2_pdf) / SQRT(2.RPI) EXP( -1.(sursat_iceliq - sursat_equ)2. / (2.sigma2_pdf) ) &
505	- qsati(i) * cldfraliq(i) * (sursat_iceliq - sursat_equ )
506
507	sigma2_icefracturb(i)= sigma2_pdf
508	mean_icefracturb(i) = sursat_equ
509
510	!------------ SPECIFIC VAPOR CONTENT AND WATER CONSERVATION ------------
511
512	IF ( (qliq_incl .LE. eps) .OR. (cldfraliq(i) .LE. eps) ) THEN
513	qliq_incl = 0.
514	cldfraliq(i) = 0.
515	END IF
516
517	!--Specific humidity is the max between qsati and the weighted mean between
518	!--qv in MPC patches and qv in ice-only parts. We assume that MPC parts are
519	!--always at qsatl and ice-only parts slightly subsaturated (qsati*sursat_iceext+1)
520	!--The whole cloud can therefore be supersaturated but never subsaturated.
521
522	qvap_incl = MAX(qsati(i), ( 1. - cldfraliq(i) ) * (sursat_iceext + 1.) * qsati(i) + cldfraliq(i) * qsatl(i) )
523
524	IF ( qvap_incl .GE. qtot_incl(i) ) THEN
525	qvap_incl = qsati(i)
526	qliq_incl = MAX(0.0,qtot_incl(i) - qvap_incl)
527	qice_incl = 0.
528
529	ELSEIF ( (qvap_incl + qliq_incl) .GE. qtot_incl(i) ) THEN
530	qliq_incl = MAX(0.0,qtot_incl(i) - qvap_incl)
531	qice_incl = 0.
532	ELSE
533	qice_incl = qtot_incl(i) - qvap_incl - qliq_incl
534	END IF
535
536	qvap_cld(i) = qvap_incl * cldfra1D
537	qliq(i) = qliq_incl * cldfra1D
538	qice(i) = qice_incl * cldfra1D
539	IF ((qice(i)+qliq(i)) .GT. 0.) THEN
540	icefrac(i) = qice(i) / ( qice(i) + qliq(i) )
541	ELSE
542	icefrac(i) = 1. ! to keep computation of qsat wrt ice in condensation loop in lmdz_lscp_main
543	ENDIF
544	dicefracdT(i) = 0.
545
546	END IF ! Enough TKE
547
548	END IF ! End qini
549
550	END IF ! ! MPC temperature
551
552	END IF ! pticefracturb and cldfra
553
554	ENDDO ! klon
555	END SUBROUTINE ICEFRAC_LSCP_TURB
556	!
557	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
558
559
560	SUBROUTINE CALC_QSAT_ECMWF(klon,temp,qtot,pressure,tref,phase,flagth,qs,dqs)
561	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
562	! Calculate qsat following ECMWF method
563	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
564
565
566	USE yoethf_mod_h
567	USE yomcst_mod_h
568	IMPLICIT NONE
569
570
571	include "FCTTRE.h"
572
573	INTEGER, INTENT(IN) :: klon ! number of horizontal grid points
574	REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature in K
575	REAL, INTENT(IN), DIMENSION(klon) :: qtot ! total specific water in kg/kg
576	REAL, INTENT(IN), DIMENSION(klon) :: pressure ! pressure in Pa
577	REAL, INTENT(IN) :: tref ! reference temperature in K
578	LOGICAL, INTENT(IN) :: flagth ! flag for qsat calculation for thermals
579	INTEGER, INTENT(IN) :: phase
580	! phase: 0=depend on temperature sign (temp>tref -> liquid, temp<tref, solid)
581	! 1=liquid
582	! 2=solid
583
584	REAL, INTENT(OUT), DIMENSION(klon) :: qs ! saturation specific humidity [kg/kg]
585	REAL, INTENT(OUT), DIMENSION(klon) :: dqs ! derivation of saturation specific humidity wrt T
586
587	REAL delta, cor, cvm5
588	INTEGER i
589
590	DO i=1,klon
591
592	IF (phase .EQ. 1) THEN
593	delta=0.
594	ELSEIF (phase .EQ. 2) THEN
595	delta=1.
596	ELSE
597	delta=MAX(0.,SIGN(1.,tref-temp(i)))
598	ENDIF
599
600	IF (flagth) THEN
601	cvm5=R5LES(1.-delta) + R5IESdelta
602	ELSE
603	cvm5 = R5LESRLVTT(1.-delta) + R5IESRLSTTdelta
604	cvm5 = cvm5 /RCPD/(1.0+RVTMP2*(qtot(i)))
605	ENDIF
606
607	qs(i)= R2ES*FOEEW(temp(i),delta)/pressure(i)
608	qs(i)=MIN(0.5,qs(i))
609	cor=1./(1.-RETV*qs(i))
610	qs(i)=qs(i)*cor
611	dqs(i)= FOEDE(temp(i),delta,cvm5,qs(i),cor)
612
613	END DO
614
615	END SUBROUTINE CALC_QSAT_ECMWF
616	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
617
618
619	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
620	SUBROUTINE CALC_GAMMASAT(klon,temp,qtot,pressure,gammasat,dgammasatdt)
621
622	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
623	! programme that calculates the gammasat parameter that determines the
624	! homogeneous condensation thresholds for cold (<0oC) clouds
625	! condensation at q>gammasat*qsat
626	! Etienne Vignon, March 2021
627	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
628
629	use lmdz_lscp_ini, only: iflag_gammasat, temp_nowater, RTT
630	use lmdz_lscp_ini, only: a_homofreez, b_homofreez, delta_hetfreez
631
632	IMPLICIT NONE
633
634
635	INTEGER, INTENT(IN) :: klon ! number of horizontal grid points
636	REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature in K
637	REAL, INTENT(IN), DIMENSION(klon) :: qtot ! total specific water in kg/kg
638
639	REAL, INTENT(IN), DIMENSION(klon) :: pressure ! pressure in Pa
640
641	REAL, INTENT(OUT), DIMENSION(klon) :: gammasat ! coefficient to multiply qsat with to calculate saturation
642	REAL, INTENT(OUT), DIMENSION(klon) :: dgammasatdt ! derivative of gammasat wrt temperature
643
644	REAL, DIMENSION(klon) :: qsi,qsl,dqsl,dqsi
645	REAL f_homofreez, fac
646
647	INTEGER i
648
649	CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,1,.false.,qsl,dqsl)
650	CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,2,.false.,qsi,dqsi)
651
652	DO i = 1, klon
653
654	IF ( temp(i) .GE. RTT ) THEN
655	! warm clouds: condensation at saturation wrt liquid
656	gammasat(i) = 1.
657	dgammasatdt(i) = 0.
658
659	ELSE
660	! cold clouds: qsi > qsl
661
662	! homogeneous freezing of aerosols, according to
663	! Koop, 2000 and Ren and MacKenzie, 2005 (QJRMS)
664	! 'Cirrus regime'
665	! if f_homofreez > qsl / qsi, liquid nucleation
666	! if f_homofreez < qsl / qsi, homogeneous freezing of aerosols
667	! Note: f_homofreez = qsl / qsi for temp ~= -38degC
668	f_homofreez = a_homofreez - temp(i) / b_homofreez
669
670	IF ( iflag_gammasat .GE. 3 ) THEN
671	! condensation at homogeneous freezing threshold for temp < -38 degC
672	! condensation at liquid saturation for temp > -38 degC
673	IF ( f_homofreez .LE. qsl(i) / qsi(i) ) THEN
674	gammasat(i) = f_homofreez
675	dgammasatdt(i) = - 1. / b_homofreez
676	ELSE
677	gammasat(i) = qsl(i) / qsi(i)
678	dgammasatdt(i) = ( dqsl(i) * qsi(i) - dqsi(i) * qsl(i) ) / qsi(i) / qsi(i)
679	ENDIF
680
681	ELSEIF ( iflag_gammasat .EQ. 2 ) THEN
682	! condensation at homogeneous freezing threshold for temp < -38 degC
683	! condensation at a threshold linearly decreasing between homogeneous
684	! freezing and ice saturation for -38 degC < temp < temp_nowater
685	! condensation at ice saturation for temp > temp_nowater
686	! If temp_nowater = 235.15 K, this is equivalent to iflag_gammasat = 1
687	IF ( f_homofreez .LE. qsl(i) / qsi(i) ) THEN
688	gammasat(i) = f_homofreez
689	dgammasatdt(i) = - 1. / b_homofreez
690	ELSEIF ( temp(i) .LE. temp_nowater ) THEN
691	! Here, we assume that f_homofreez = qsl / qsi for temp = -38 degC = 235.15 K
692	dgammasatdt(i) = ( a_homofreez - 235.15 / b_homofreez - 1. ) &
693	/ ( 235.15 - temp_nowater )
694	gammasat(i) = dgammasatdt(i) * ( temp(i) - temp_nowater ) + 1.
695	ELSE
696	gammasat(i) = 1.
697	dgammasatdt(i) = 0.
698	ENDIF
699
700	ELSEIF ( iflag_gammasat .EQ. 1 ) THEN
701	! condensation at homogeneous freezing threshold for temp < -38 degC
702	! condensation at ice saturation for temp > -38 degC
703	IF ( f_homofreez .LE. qsl(i) / qsi(i) ) THEN
704	gammasat(i) = f_homofreez
705	dgammasatdt(i) = - 1. / b_homofreez
706	ELSE
707	gammasat(i) = 1.
708	dgammasatdt(i) = 0.
709	ENDIF
710
711	ELSE
712	! condensation at ice saturation for temp < -38 degC
713	! condensation at ice saturation for temp > -38 degC
714	gammasat(i) = 1.
715	dgammasatdt(i) = 0.
716
717	ENDIF
718
719	! Note that the delta_hetfreez parameter allows to linearly decrease the
720	! condensation threshold between the calculated threshold and the ice saturation
721	! for delta_hetfreez = 1, the threshold is the calculated condensation threshold
722	! for delta_hetfreez = 0, the threshold is the ice saturation
723	gammasat(i) = ( 1. - delta_hetfreez ) + delta_hetfreez * gammasat(i)
724	dgammasatdt(i) = delta_hetfreez * dgammasatdt(i)
725
726	ENDIF
727
728	END DO
729
730
731	END SUBROUTINE CALC_GAMMASAT
732	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
733
734
735	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
736	SUBROUTINE DISTANCE_TO_CLOUD_TOP(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D,temp_cltop)
737	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
738
739	USE lmdz_lscp_ini, ONLY : rd,rg,tresh_cl
740
741	IMPLICIT NONE
742
743	INTEGER, INTENT(IN) :: klon,klev !number of horizontal and vertical grid points
744	INTEGER, INTENT(IN) :: k ! vertical index
745	REAL, INTENT(IN), DIMENSION(klon,klev) :: temp ! temperature in K
746	REAL, INTENT(IN), DIMENSION(klon,klev) :: pplay ! pressure middle layer in Pa
747	REAL, INTENT(IN), DIMENSION(klon,klev+1) :: paprs ! pressure interfaces in Pa
748	REAL, INTENT(IN), DIMENSION(klon,klev) :: rneb ! cloud fraction
749
750	REAL, INTENT(OUT), DIMENSION(klon) :: distcltop1D ! distance from cloud top
751	REAL, INTENT(OUT), DIMENSION(klon) :: temp_cltop ! temperature of cloud top
752
753	REAL dzlay(klon,klev)
754	REAL zlay(klon,klev)
755	REAL dzinterf
756	INTEGER i,k_top, kvert
757	LOGICAL bool_cl
758
759
760	DO i=1,klon
761	! Initialization height middle of first layer
762	dzlay(i,1) = Rd * temp(i,1) / rg * log(paprs(i,1)/paprs(i,2))
763	zlay(i,1) = dzlay(i,1)/2
764
765	DO kvert=2,klev
766	IF (kvert.EQ.klev) THEN
767	dzlay(i,kvert) = 2(rd temp(i,kvert) / rg * log(paprs(i,kvert)/pplay(i,kvert)))
768	ELSE
769	dzlay(i,kvert) = rd * temp(i,kvert) / rg * log(paprs(i,kvert)/paprs(i,kvert+1))
770	ENDIF
771	dzinterf = rd * temp(i,kvert) / rg * log(pplay(i,kvert-1)/pplay(i,kvert))
772	zlay(i,kvert) = zlay(i,kvert-1) + dzinterf
773	ENDDO
774	ENDDO
775
776	DO i=1,klon
777	k_top = k
778	IF (rneb(i,k) .LE. tresh_cl) THEN
779	bool_cl = .FALSE.
780	ELSE
781	bool_cl = .TRUE.
782	ENDIF
783
784	DO WHILE ((bool_cl) .AND. (k_top .LE. klev))
785	! find cloud top
786	IF (rneb(i,k_top) .GT. tresh_cl) THEN
787	k_top = k_top + 1
788	ELSE
789	bool_cl = .FALSE.
790	k_top = k_top - 1
791	ENDIF
792	ENDDO
793	k_top=min(k_top,klev)
794
795	!dist to top is dist between current layer and layer of cloud top (from middle to middle) + dist middle to
796	!interf for layer of cloud top
797	distcltop1D(i) = zlay(i,k_top) - zlay(i,k) + dzlay(i,k_top)/2
798	temp_cltop(i) = temp(i,k_top)
799	ENDDO ! klon
800
801	END SUBROUTINE DISTANCE_TO_CLOUD_TOP
802	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
803
804	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
805	FUNCTION GAMMAINC ( p, x )
806
807	!*****************************************************************************80
808	!
809	!! GAMMAINC computes the regularized lower incomplete Gamma Integral
810	!
811	! Modified:
812	!
813	! 20 January 2008
814	!
815	! Author:
816	!
817	! Original FORTRAN77 version by B Shea.
818	! FORTRAN90 version by John Burkardt.
819	!
820	! Reference:
821	!
822	! B Shea,
823	! Algorithm AS 239:
824	! Chi-squared and Incomplete Gamma Integral,
825	! Applied Statistics,
826	! Volume 37, Number 3, 1988, pages 466-473.
827	!
828	! Parameters:
829	!
830	! Input, real X, P, the parameters of the incomplete
831	! gamma ratio. 0 <= X, and 0 < P.
832	!
833	! Output, real GAMMAINC, the value of the incomplete
834	! Gamma integral.
835	!
836	IMPLICIT NONE
837
838	REAL A
839	REAL AN
840	REAL ARG
841	REAL B
842	REAL C
843	REAL, PARAMETER :: ELIMIT = - 88.0E+00
844	REAL GAMMAINC
845	REAL, PARAMETER :: OFLO = 1.0E+37
846	REAL P
847	REAL, PARAMETER :: PLIMIT = 1000.0E+00
848	REAL PN1
849	REAL PN2
850	REAL PN3
851	REAL PN4
852	REAL PN5
853	REAL PN6
854	REAL RN
855	REAL, PARAMETER :: TOL = 1.0E-14
856	REAL X
857	REAL, PARAMETER :: XBIG = 1.0E+08
858
859	GAMMAINC = 0.0E+00
860
861	IF ( X == 0.0E+00 ) THEN
862	GAMMAINC = 0.0E+00
863	RETURN
864	END IF
865	!
866	! IF P IS LARGE, USE A NORMAL APPROXIMATION.
867	!
868	IF ( PLIMIT < P ) THEN
869
870	PN1 = 3.0E+00 * SQRT ( P ) * ( ( X / P )**( 1.0E+00 / 3.0E+00 ) &
871	+ 1.0E+00 / ( 9.0E+00 * P ) - 1.0E+00 )
872
873	GAMMAINC = 0.5E+00 * ( 1. + ERF ( PN1 ) )
874	RETURN
875
876	END IF
877	!
878	! IF X IS LARGE SET GAMMAD = 1.
879	!
880	IF ( XBIG < X ) THEN
881	GAMMAINC = 1.0E+00
882	RETURN
883	END IF
884	!
885	! USE PEARSON'S SERIES EXPANSION.
886	! (NOTE THAT P IS NOT LARGE ENOUGH TO FORCE OVERFLOW IN ALOGAM).
887	!
888	IF ( X <= 1.0E+00 .OR. X < P ) THEN
889
890	ARG = P * LOG ( X ) - X - LOG_GAMMA ( P + 1.0E+00 )
891	C = 1.0E+00
892	GAMMAINC = 1.0E+00
893	A = P
894
895	DO
896
897	A = A + 1.0E+00
898	C = C * X / A
899	GAMMAINC = GAMMAINC + C
900
901	IF ( C <= TOL ) THEN
902	EXIT
903	END IF
904
905	END DO
906
907	ARG = ARG + LOG ( GAMMAINC )
908
909	IF ( ELIMIT <= ARG ) THEN
910	GAMMAINC = EXP ( ARG )
911	ELSE
912	GAMMAINC = 0.0E+00
913	END IF
914	!
915	! USE A CONTINUED FRACTION EXPANSION.
916	!
917	ELSE
918
919	ARG = P * LOG ( X ) - X - LOG_GAMMA ( P )
920	A = 1.0E+00 - P
921	B = A + X + 1.0E+00
922	C = 0.0E+00
923	PN1 = 1.0E+00
924	PN2 = X
925	PN3 = X + 1.0E+00
926	PN4 = X * B
927	GAMMAINC = PN3 / PN4
928
929	DO
930
931	A = A + 1.0E+00
932	B = B + 2.0E+00
933	C = C + 1.0E+00
934	AN = A * C
935	PN5 = B * PN3 - AN * PN1
936	PN6 = B * PN4 - AN * PN2
937
938	IF ( PN6 /= 0.0E+00 ) THEN
939
940	RN = PN5 / PN6
941
942	IF ( ABS ( GAMMAINC - RN ) <= MIN ( TOL, TOL * RN ) ) THEN
943	EXIT
944	END IF
945
946	GAMMAINC = RN
947
948	END IF
949
950	PN1 = PN3
951	PN2 = PN4
952	PN3 = PN5
953	PN4 = PN6
954	!
955	! RE-SCALE TERMS IN CONTINUED FRACTION IF TERMS ARE LARGE.
956	!
957	IF ( OFLO <= ABS ( PN5 ) ) THEN
958	PN1 = PN1 / OFLO
959	PN2 = PN2 / OFLO
960	PN3 = PN3 / OFLO
961	PN4 = PN4 / OFLO
962	END IF
963
964	END DO
965
966	ARG = ARG + LOG ( GAMMAINC )
967
968	IF ( ELIMIT <= ARG ) THEN
969	GAMMAINC = 1.0E+00 - EXP ( ARG )
970	ELSE
971	GAMMAINC = 1.0E+00
972	END IF
973
974	END IF
975
976	RETURN
977	END FUNCTION GAMMAINC
978	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
979
980	END MODULE lmdz_lscp_tools
981
982

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

Download in other formats:

Original Format