Context Navigation

lmdz_lscp_tools.F90 @ 5169

Last change on this file since 5169 was 5007, checked in by evignon, 12 months ago
ajout de la nouvelle paramétrisation du partitonnement de phase dans les nuages de phase mixte de Lea Raillard La parametrisation s'active avec iflag_icefrac=1 et est fondé sur la theorie de creation et maintien de sursaturation en atmosphere turbulente avec ou sans presence de cristaux de glace
File size: 29.7 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)
246	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
247
248
249	USE lmdz_lscp_ini, ONLY : prt_level, lunout
250	USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RV, RPI
251	USE lmdz_lscp_ini, ONLY : seuil_neb, temp_nowater
252	USE lmdz_lscp_ini, ONLY : tau_mixenv, lmix_mpc, naero5, gamma_snwretro, gamma_taud, capa_crystal
253	USE lmdz_lscp_ini, ONLY : eps
254
255	IMPLICIT NONE
256
257	INTEGER, INTENT(IN) :: klon !--number of horizontal grid points
258	REAL, INTENT(IN) :: dtime !--time step [s]
259
260	REAL, INTENT(IN), DIMENSION(klon) :: temp !--temperature
261	REAL, INTENT(IN), DIMENSION(klon) :: pplay !--pressure in the middle of the layer [Pa]
262	REAL, INTENT(IN), DIMENSION(klon) :: paprsdn !--pressure at the bottom interface of the layer [Pa]
263	REAL, INTENT(IN), DIMENSION(klon) :: paprsup !--pressure at the top interface of the layer [Pa]
264	REAL, INTENT(IN), DIMENSION(klon) :: qtot_incl !--specific total cloud water content, in-cloud content [kg/kg]
265	REAL, INTENT(IN), DIMENSION(klon) :: cldfra !--cloud fraction in gridbox [-]
266	REAL, INTENT(IN), DIMENSION(klon) :: tke !--turbulent kinetic energy [m2/s2]
267	REAL, INTENT(IN), DIMENSION(klon) :: tke_dissip !--TKE dissipation [m2/s3]
268
269	REAL, INTENT(IN), DIMENSION(klon) :: qice_ini !--initial specific ice content gridbox-mean [kg/kg]
270	REAL, INTENT(IN), DIMENSION(klon) :: snowcld
271	REAL, INTENT(OUT), DIMENSION(klon) :: qliq !--specific liquid content gridbox-mean [kg/kg]
272	REAL, INTENT(OUT), DIMENSION(klon) :: qvap_cld !--specific cloud vapor content, gridbox-mean [kg/kg]
273	REAL, INTENT(OUT), DIMENSION(klon) :: qice !--specific ice content gridbox-mean [kg/kg]
274	REAL, INTENT(OUT), DIMENSION(klon) :: icefrac !--fraction of ice in condensed water [-]
275	REAL, INTENT(OUT), DIMENSION(klon) :: dicefracdT
276
277	REAL, INTENT(OUT), DIMENSION(klon) :: cldfraliq !--fraction of cldfra which is liquid only
278	REAL, INTENT(OUT), DIMENSION(klon) :: sigma2_icefracturb !--Temporary
279	REAL, INTENT(OUT), DIMENSION(klon) :: mean_icefracturb !--Temporary
280
281	REAL, DIMENSION(klon) :: qzero, qsatl, dqsatl, qsati, dqsati !--specific humidity saturation values
282	INTEGER :: i
283
284	REAL :: qvap_incl, qice_incl, qliq_incl, qiceini_incl !--In-cloud specific quantities [kg/kg]
285	REAL :: qsnowcld_incl
286	!REAL :: capa_crystal !--Capacitance of ice crystals [-]
287	REAL :: water_vapor_diff !--Water-vapour diffusion coefficient in air [m2/s] (function of T&P)
288	REAL :: air_thermal_conduct !--Thermal conductivity of air [J/m/K/s] (function of T)
289	REAL :: C0 !--Lagrangian structure function [-]
290	REAL :: tau_mixingenv
291	REAL :: tau_dissipturb
292	REAL :: invtau_phaserelax
293	REAL :: sigma2_pdf, mean_pdf
294	REAL :: ai, bi, B0
295	REAL :: sursat_iceliq
296	REAL :: sursat_env
297	REAL :: liqfra_max
298	REAL :: sursat_iceext
299	REAL :: nb_crystals !--number concentration of ice crystals [#/m3]
300	REAL :: moment1_PSD !--1st moment of ice PSD
301	REAL :: N0_PSD, lambda_PSD !--parameters of the exponential PSD
302
303	REAL :: rho_ice !--ice density [kg/m3]
304	REAL :: cldfra1D
305	REAL :: deltaz, rho_air
306	REAL :: psati !--saturation vapor pressure wrt i [Pa]
307
308	C0 = 10. !--value assumed in Field2014
309	rho_ice = 950.
310	sursat_iceext = -0.1
311	!capa_crystal = 1. !r_ice
312	qzero(:) = 0.
313	cldfraliq(:) = 0.
314	icefrac(:) = 0.
315	dicefracdT(:) = 0.
316
317	sigma2_icefracturb(:) = 0.
318	mean_icefracturb(:) = 0.
319
320	!--wrt liquid water
321	CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,1,.false.,qsatl(:),dqsatl(:))
322	!--wrt ice
323	CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,2,.false.,qsati(:),dqsati(:))
324
325
326	DO i=1,klon
327
328
329	rho_air = pplay(i) / temp(i) / RD
330	!deltaz = ( paprsdn(i) - paprsup(i) ) / RG / rho_air(i)
331	! because cldfra is intent in, but can be locally modified due to test
332	cldfra1D = cldfra(i)
333	IF (cldfra(i) .LE. 0.) THEN
334	qvap_cld(i) = 0.
335	qliq(i) = 0.
336	qice(i) = 0.
337	cldfraliq(i) = 0.
338	icefrac(i) = 0.
339	dicefracdT(i) = 0.
340
341	! If there is a cloud
342	ELSE
343	IF (cldfra(i) .GE. 1.0) THEN
344	cldfra1D = 1.0
345	END IF
346
347	! T>0°C, no ice allowed
348	IF ( temp(i) .GE. RTT ) THEN
349	qvap_cld(i) = qsatl(i) * cldfra1D
350	qliq(i) = MAX(0.0,qtot_incl(i)-qsatl(i)) * cldfra1D
351	qice(i) = 0.
352	cldfraliq(i) = 1.
353	icefrac(i) = 0.
354	dicefracdT(i) = 0.
355
356	! T<-38°C, no liquid allowed
357	ELSE IF ( temp(i) .LE. temp_nowater) THEN
358	qvap_cld(i) = qsati(i) * cldfra1D
359	qliq(i) = 0.
360	qice(i) = MAX(0.0,qtot_incl(i)-qsati(i)) * cldfra1D
361	cldfraliq(i) = 0.
362	icefrac(i) = 1.
363	dicefracdT(i) = 0.
364
365	! MPC temperature
366	ELSE
367	! Not enough TKE
368	IF ( tke_dissip(i) .LE. eps ) THEN
369	qvap_cld(i) = qsati(i) * cldfra1D
370	qliq(i) = 0.
371	qice(i) = MAX(0.,qtot_incl(i)-qsati(i)) * cldfra1D
372	cldfraliq(i) = 0.
373	icefrac(i) = 1.
374	dicefracdT(i) = 0.
375
376	! Enough TKE
377	ELSE
378	!---------------------------------------------------------
379	!-- ICE SUPERSATURATION PDF
380	!---------------------------------------------------------
381	!--If -38°C< T <0°C and there is enough turbulence,
382	!--we compute the cloud liquid properties with a Gaussian PDF
383	!--of ice supersaturation F(Si) (Field2014, Furtado2016).
384	!--Parameters of the PDF are function of turbulence and
385	!--microphysics/existing ice.
386
387	sursat_iceliq = qsatl(i)/qsati(i) - 1.
388	psati = qsati(i) * pplay(i) / (RD/RV)
389
390	!-------------- MICROPHYSICAL TERMS --------------
391	!--We assume an exponential ice PSD whose parameters
392	!--are computed following Morrison&Gettelman 2008
393	!--Ice number density is assumed equals to INP density
394	!--which is a function of temperature (DeMott 2010)
395	!--bi and B0 are microphysical function characterizing
396	!--vapor/ice interactions
397	!--tau_phase_relax is the typical time of vapor deposition
398	!--onto ice crystals
399
400	qiceini_incl = qice_ini(i) / cldfra1D
401	qsnowcld_incl = snowcld(i) * RG * dtime / ( paprsdn(i) - paprsup(i) ) / cldfra1D
402	sursat_env = max(0., (qtot_incl(i) - qiceini_incl)/qsati(i) - 1.)
403	IF ( qiceini_incl .GT. eps ) THEN
404	nb_crystals = 1.e3 * 5.94e-5 * ( RTT - temp(i) )*3.33 naero5*(0.0264(RTT-temp(i))+0.0033)
405	lambda_PSD = ( (RPIrho_icenb_crystals24.) / (6.(qiceini_incl + gamma_snwretro * qsnowcld_incl)) ) ** (1./3.)
406	N0_PSD = nb_crystals * lambda_PSD
407	moment1_PSD = N0_PSD/2./lambda_PSD**2
408	ELSE
409	moment1_PSD = 0.
410	ENDIF
411
412	!--Formulae for air thermal conductivity and water vapor diffusivity
413	!--comes respectively from Beard and Pruppacher (1971)
414	!--and Hall and Pruppacher (1976)
415
416	air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184
417	water_vapor_diff = 2.111e-5 ( temp(i) / RTT )*1.94 ( 101325 / pplay(i) )
418
419	bi = 1./((qsati(i)+qsatl(i))/2.) + RLSTT2 / RCPD / RV / temp(i)2
420	B0 = 4. * RPI * capa_crystal * 1. / ( RLSTT2 / air_thermal_conduct / RV / temp(i)2 &
421	+ RV * temp(i) / psati / water_vapor_diff )
422
423	invtau_phaserelax = (bi * B0 * moment1_PSD )
424
425	! Old way of estimating moment1 : spherical crystals + monodisperse PSD
426	! nb_crystals = rho_air * qiceini_incl / ( 4. / 3. * RPI * r_ice*3. rho_ice )
427	! moment1_PSD = nb_crystals * r_ice
428
429	!----------------- TURBULENT SOURCE/SINK TERMS -----------------
430	!--Tau_mixingenv is the time needed to homogeneize the parcel
431	!--with its environment by turbulent diffusion over the parcel
432	!--length scale
433	!--if lmix_mpc <0, tau_mixigenv value is prescribed
434	!--else tau_mixigenv value is derived from tke_dissip and lmix_mpc
435	!--Tau_dissipturb is the time needed turbulence to decay due to
436	!--viscosity
437
438	ai = RG / RD / temp(i) * ( RD * RLSTT / RCPD / RV / temp(i) - 1. )
439	IF ( lmix_mpc .GT. 0 ) THEN
440	tau_mixingenv = ( lmix_mpc2. / tke_dissip(i) )(1./3.)
441	ELSE
442	tau_mixingenv = tau_mixenv
443	ENDIF
444
445	tau_dissipturb = gamma_taud * 2. * 2./3. * tke(i) / tke_dissip(i) / C0
446
447	!--------------------- PDF COMPUTATIONS ---------------------
448	!--Formulae for sigma2_pdf (variance), mean of PDF in Furtado2016
449	!--cloud liquid fraction and in-cloud liquid content are given
450	!--by integrating resp. F(Si) and Si*F(Si)
451	!--Liquid is limited by the available water vapor trough a
452	!--maximal liquid fraction
453
454	liqfra_max = MAX(0., (MIN (1.,( qtot_incl(i) - qiceini_incl - qsati(i) * (1 + sursat_iceext ) ) / ( qsatl(i) - qsati(i) ) ) ) )
455	sigma2_pdf = 1./2. * ( ai*2 ) 2./3. * tke(i) * tau_dissipturb / ( invtau_phaserelax + 1./tau_mixingenv )
456	mean_pdf = sursat_env * 1./tau_mixingenv / ( invtau_phaserelax + 1./tau_mixingenv )
457	cldfraliq(i) = 0.5 * (1. - erf( ( sursat_iceliq - mean_pdf) / (SQRT(2.* sigma2_pdf) ) ) )
458	!IF (cldfraliq(i) .GT. liqfra_max) THEN
459	! cldfraliq(i) = liqfra_max
460	!ENDIF
461
462	qliq_incl = qsati(i) * SQRT(sigma2_pdf) / SQRT(2.RPI) EXP( -1.(sursat_iceliq - mean_pdf)2. / (2.sigma2_pdf) ) &
463	- qsati(i) * cldfraliq(i) * (sursat_iceliq - mean_pdf )
464
465	sigma2_icefracturb(i)= sigma2_pdf
466	mean_icefracturb(i) = mean_pdf
467	!------------ ICE AMOUNT AND WATER CONSERVATION ------------
468
469	IF ( (qliq_incl .LE. eps) .OR. (cldfraliq(i) .LE. eps) ) THEN
470	qliq_incl = 0.
471	cldfraliq(i) = 0.
472	END IF
473
474	!--Choice for in-cloud vapor :
475	!--1.Weighted mean between qvap in MPC parts and in ice-only parts
476	!--2.Always at ice saturation
477	qvap_incl = MAX(qsati(i), ( 1. - cldfraliq(i) ) * (sursat_iceext + 1.) * qsati(i) + cldfraliq(i) * qsatl(i) )
478
479	IF ( qvap_incl .GE. qtot_incl(i) ) THEN
480	qvap_incl = qsati(i)
481	qliq_incl = qtot_incl(i) - qvap_incl
482	qice_incl = 0.
483
484	ELSEIF ( (qvap_incl + qliq_incl) .GE. qtot_incl(i) ) THEN
485	qliq_incl = MAX(0.0,qtot_incl(i) - qvap_incl)
486	qice_incl = 0.
487	ELSE
488	qice_incl = qtot_incl(i) - qvap_incl - qliq_incl
489	END IF
490
491	qvap_cld(i) = qvap_incl * cldfra1D
492	qliq(i) = qliq_incl * cldfra1D
493	qice(i) = qice_incl * cldfra1D
494	icefrac(i) = qice(i) / ( qice(i) + qliq(i) )
495	dicefracdT(i) = 0.
496	!print*,'MPC turb'
497
498	END IF ! Enough TKE
499
500	END IF ! ! MPC temperature
501
502	END IF ! cldfra
503
504	ENDDO ! klon
505	END SUBROUTINE ICEFRAC_LSCP_TURB
506	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
507
508
509	SUBROUTINE CALC_QSAT_ECMWF(klon,temp,qtot,pressure,tref,phase,flagth,qs,dqs)
510	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
511	! Calculate qsat following ECMWF method
512	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
513
514
515	IMPLICIT NONE
516
517	include "YOMCST.h"
518	include "YOETHF.h"
519	include "FCTTRE.h"
520
521	INTEGER, INTENT(IN) :: klon ! number of horizontal grid points
522	REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature in K
523	REAL, INTENT(IN), DIMENSION(klon) :: qtot ! total specific water in kg/kg
524	REAL, INTENT(IN), DIMENSION(klon) :: pressure ! pressure in Pa
525	REAL, INTENT(IN) :: tref ! reference temperature in K
526	LOGICAL, INTENT(IN) :: flagth ! flag for qsat calculation for thermals
527	INTEGER, INTENT(IN) :: phase
528	! phase: 0=depend on temperature sign (temp>tref -> liquid, temp<tref, solid)
529	! 1=liquid
530	! 2=solid
531
532	REAL, INTENT(OUT), DIMENSION(klon) :: qs ! saturation specific humidity [kg/kg]
533	REAL, INTENT(OUT), DIMENSION(klon) :: dqs ! derivation of saturation specific humidity wrt T
534
535	REAL delta, cor, cvm5
536	INTEGER i
537
538	DO i=1,klon
539
540	IF (phase .EQ. 1) THEN
541	delta=0.
542	ELSEIF (phase .EQ. 2) THEN
543	delta=1.
544	ELSE
545	delta=MAX(0.,SIGN(1.,tref-temp(i)))
546	ENDIF
547
548	IF (flagth) THEN
549	cvm5=R5LES(1.-delta) + R5IESdelta
550	ELSE
551	cvm5 = R5LESRLVTT(1.-delta) + R5IESRLSTTdelta
552	cvm5 = cvm5 /RCPD/(1.0+RVTMP2*(qtot(i)))
553	ENDIF
554
555	qs(i)= R2ES*FOEEW(temp(i),delta)/pressure(i)
556	qs(i)=MIN(0.5,qs(i))
557	cor=1./(1.-RETV*qs(i))
558	qs(i)=qs(i)*cor
559	dqs(i)= FOEDE(temp(i),delta,cvm5,qs(i),cor)
560
561	END DO
562
563	END SUBROUTINE CALC_QSAT_ECMWF
564	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
565
566
567	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
568	SUBROUTINE CALC_GAMMASAT(klon,temp,qtot,pressure,gammasat,dgammasatdt)
569
570	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
571	! programme that calculates the gammasat parameter that determines the
572	! homogeneous condensation thresholds for cold (<0oC) clouds
573	! condensation at q>gammasat*qsat
574	! Etienne Vignon, March 2021
575	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
576
577	use lmdz_lscp_ini, only: iflag_gammasat, t_glace_min, RTT
578
579	IMPLICIT NONE
580
581
582	INTEGER, INTENT(IN) :: klon ! number of horizontal grid points
583	REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature in K
584	REAL, INTENT(IN), DIMENSION(klon) :: qtot ! total specific water in kg/kg
585
586	REAL, INTENT(IN), DIMENSION(klon) :: pressure ! pressure in Pa
587
588	REAL, INTENT(OUT), DIMENSION(klon) :: gammasat ! coefficient to multiply qsat with to calculate saturation
589	REAL, INTENT(OUT), DIMENSION(klon) :: dgammasatdt ! derivative of gammasat wrt temperature
590
591	REAL, DIMENSION(klon) :: qsi,qsl,dqsl,dqsi
592	REAL fcirrus, fac
593	REAL, PARAMETER :: acirrus=2.349
594	REAL, PARAMETER :: bcirrus=259.0
595
596	INTEGER i
597
598	CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,1,.false.,qsl,dqsl)
599	CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,2,.false.,qsi,dqsi)
600
601	DO i=1,klon
602
603	IF (temp(i) .GE. RTT) THEN
604	! warm clouds: condensation at saturation wrt liquid
605	gammasat(i)=1.
606	dgammasatdt(i)=0.
607
608	ELSEIF ((temp(i) .LT. RTT) .AND. (temp(i) .GT. t_glace_min)) THEN
609
610	IF (iflag_gammasat .GE. 2) THEN
611	gammasat(i)=qsl(i)/qsi(i)
612	dgammasatdt(i)=(dqsl(i)qsi(i)-dqsi(i)qsl(i))/qsi(i)/qsi(i)
613	ELSE
614	gammasat(i)=1.
615	dgammasatdt(i)=0.
616	ENDIF
617
618	ELSE
619
620	IF (iflag_gammasat .GE.1) THEN
621	! homogeneous freezing of aerosols, according to
622	! Koop, 2000 and Karcher 2008, QJRMS
623	! 'Cirrus regime'
624	fcirrus=acirrus-temp(i)/bcirrus
625	IF (fcirrus .GT. qsl(i)/qsi(i)) THEN
626	gammasat(i)=qsl(i)/qsi(i)
627	dgammasatdt(i)=(dqsl(i)qsi(i)-dqsi(i)qsl(i))/qsi(i)/qsi(i)
628	ELSE
629	gammasat(i)=fcirrus
630	dgammasatdt(i)=-1.0/bcirrus
631	ENDIF
632
633	ELSE
634
635	gammasat(i)=1.
636	dgammasatdt(i)=0.
637
638	ENDIF
639
640	ENDIF
641
642	END DO
643
644
645	END SUBROUTINE CALC_GAMMASAT
646	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
647
648
649	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
650	SUBROUTINE DISTANCE_TO_CLOUD_TOP(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D,temp_cltop)
651	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
652
653	USE lmdz_lscp_ini, ONLY : rd,rg,tresh_cl
654
655	IMPLICIT NONE
656
657	INTEGER, INTENT(IN) :: klon,klev !number of horizontal and vertical grid points
658	INTEGER, INTENT(IN) :: k ! vertical index
659	REAL, INTENT(IN), DIMENSION(klon,klev) :: temp ! temperature in K
660	REAL, INTENT(IN), DIMENSION(klon,klev) :: pplay ! pressure middle layer in Pa
661	REAL, INTENT(IN), DIMENSION(klon,klev+1) :: paprs ! pressure interfaces in Pa
662	REAL, INTENT(IN), DIMENSION(klon,klev) :: rneb ! cloud fraction
663
664	REAL, INTENT(OUT), DIMENSION(klon) :: distcltop1D ! distance from cloud top
665	REAL, INTENT(OUT), DIMENSION(klon) :: temp_cltop ! temperature of cloud top
666
667	REAL dzlay(klon,klev)
668	REAL zlay(klon,klev)
669	REAL dzinterf
670	INTEGER i,k_top, kvert
671	LOGICAL bool_cl
672
673
674	DO i=1,klon
675	! Initialization height middle of first layer
676	dzlay(i,1) = Rd * temp(i,1) / rg * log(paprs(i,1)/paprs(i,2))
677	zlay(i,1) = dzlay(i,1)/2
678
679	DO kvert=2,klev
680	IF (kvert.EQ.klev) THEN
681	dzlay(i,kvert) = 2(rd temp(i,kvert) / rg * log(paprs(i,kvert)/pplay(i,kvert)))
682	ELSE
683	dzlay(i,kvert) = rd * temp(i,kvert) / rg * log(paprs(i,kvert)/paprs(i,kvert+1))
684	ENDIF
685	dzinterf = rd * temp(i,kvert) / rg * log(pplay(i,kvert-1)/pplay(i,kvert))
686	zlay(i,kvert) = zlay(i,kvert-1) + dzinterf
687	ENDDO
688	ENDDO
689
690	DO i=1,klon
691	k_top = k
692	IF (rneb(i,k) .LE. tresh_cl) THEN
693	bool_cl = .FALSE.
694	ELSE
695	bool_cl = .TRUE.
696	ENDIF
697
698	DO WHILE ((bool_cl) .AND. (k_top .LE. klev))
699	! find cloud top
700	IF (rneb(i,k_top) .GT. tresh_cl) THEN
701	k_top = k_top + 1
702	ELSE
703	bool_cl = .FALSE.
704	k_top = k_top - 1
705	ENDIF
706	ENDDO
707	k_top=min(k_top,klev)
708
709	!dist to top is dist between current layer and layer of cloud top (from middle to middle) + dist middle to
710	!interf for layer of cloud top
711	distcltop1D(i) = zlay(i,k_top) - zlay(i,k) + dzlay(i,k_top)/2
712	temp_cltop(i) = temp(i,k_top)
713	ENDDO ! klon
714
715	END SUBROUTINE DISTANCE_TO_CLOUD_TOP
716	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
717
718	END MODULE lmdz_lscp_tools
719
720

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

Download in other formats:

Original Format