Context Navigation

lmdz_lscp_tools.F90 @ 5117

Last change on this file since 5117 was 5117, checked in by abarral, 4 months ago
rename modules properly lmdz_* move some unused files to obsolete/ (lint) uppercase fortran keywords
File size: 29.6 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 == 1) THEN
52	! so-called 'empirical parameterization' in Stubenrauch et al. 2019
53	IF (tempc >= -60.0) THEN
54	logvm= -0.0000414122tempctempclog(iwcg)-0.00538922tempc*log(iwcg) &
55	-0.0516344log(iwcg)+0.00216078tempc + 1.9714
56	velo(i)=exp(logvm)
57	else
58	velo(i)=65.0(iwcg0.2)(150./phpa)**0.15
59	endif
60
61	velo(i)=fallv_tun*velo(i)/100.0 ! from cm/s to m/s
62
63	ELSE IF (iflag_vice == 2) THEN
64	! so called PSDM empirical coherent bulk ice scheme in Stubenrauch et al. 2019
65	aice=0.587
66	bice=2.45
67	asnow=0.0444
68	bsnow=2.1
69
70	m2ice=((iwcg0.001/aice)/(exp(13.6-bice7.76+0.479bice2) &
71	exp((-0.0361+bice0.0151+0.00149bice*2)tempc))) &
72	*(1./(0.807+bice0.00581+0.0457bice*2))
73
74	vmice=100.1042.4exp(13.6-(bice+1)7.76+0.479(bice+1.)*2)exp((-0.0361+ &
75	(bice+1.)0.0151+0.00149(bice+1.)*2)tempc) &
76	(m2ice(0.807+(bice+1.)0.00581+0.0457(bice+1.)2))/(iwcg0.001/aice)
77
78
79	vmice=vmice((1000./phpa)*0.2)
80
81	m2snow=((iwcg0.001/asnow)/(exp(13.6-bsnow7.76+0.479bsnow2) &
82	exp((-0.0361+bsnow0.0151+0.00149bsnow*2)tempc))) &
83	*(1./(0.807+bsnow0.00581+0.0457bsnow*2))
84
85
86	vmsnow=100.14.3exp(13.6-(bsnow+.416)7.76+0.479(bsnow+.416)**2)&
87	exp((-0.0361+(bsnow+.416)0.0151+0.00149(bsnow+.416)2)tempc)&
88	(m2snow(0.807+(bsnow+.416)0.00581+0.0457(bsnow+.416)2))/(iwcg0.001/asnow)
89
90	vmsnow=vmsnow((1000./phpa)*0.35)
91	velo(i)=fallv_tun*min(vmsnow,vmice)/100. ! to m/s
92
93	ELSE
94	! By default, fallspeed velocity of ice crystals according to Heymsfield & Donner 1990
95	velo(i) = fallv_tuncice_velo((iwcg/1000.)**dice_velo)
96	ENDIF
97	ENDDO
98
99	END SUBROUTINE FALLICE_VELOCITY
100	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
101
102	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
103	SUBROUTINE ICEFRAC_LSCP(klon, temp, iflag_ice_thermo, distcltop, temp_cltop, icefrac, dicefracdT)
104	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
105
106	! Compute the ice fraction 1-xliq (see e.g.
107	! Doutriaux-Boucher & Quaas 2004, section 2.2.)
108	! as a function of temperature
109	! see also Fig 3 of Madeleine et al. 2020, JAMES
110	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
111
112
113	USE lmdz_print_control, ONLY: lunout, prt_level
114	USE lmdz_lscp_ini, ONLY: t_glace_min, t_glace_max, exposant_glace, iflag_t_glace
115	USE lmdz_lscp_ini, ONLY: RTT, dist_liq, temp_nowater
116	USE lmdz_abort_physic, ONLY: abort_physic
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<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 == 1).OR.(iflag_ice_thermo >= 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 == 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) >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)==0).OR.(icefrac(i)==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 == 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) >0.) .AND. (liqfrac_tmp > 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 == 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 <=0) .OR. (liqfrac_tmp >= 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 == 5) .OR. (iflag_t_glace == 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 <=0) .OR. (liqfrac_tmp >= 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 == 6) THEN
212	IF (temp(i) > 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 <=0) .OR. (liqfrac_tmp >= 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 >= 4) THEN
229	IF ((temp_cltop(i) <= temp_nowater) .AND. (temp(i) <= t_glace_max)) THEN
230	icefrac(i) = 1.
231	dicefracdT(i) = 0.
232	ENDIF
233	ENDIF
234
235
236	ENDDO ! klon
237
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) <= 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) >= 1.0) THEN
344	cldfra1D = 1.0
345	END IF
346
347	! T>0°C, no ice allowed
348	IF ( temp(i) >= 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) <= 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) <= 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 > 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 > 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 <= eps) .OR. (cldfraliq(i) <= 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 >= 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) >= 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 == 1) THEN
541	delta=0.
542	ELSEIF (phase == 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) >= RTT) THEN
604	! warm clouds: condensation at saturation wrt liquid
605	gammasat(i)=1.
606	dgammasatdt(i)=0.
607
608	ELSEIF ((temp(i) < RTT) .AND. (temp(i) > t_glace_min)) THEN
609
610	IF (iflag_gammasat >= 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 >=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 > 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==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) <= tresh_cl) THEN
693	bool_cl = .FALSE.
694	ELSE
695	bool_cl = .TRUE.
696	ENDIF
697
698	DO WHILE ((bool_cl) .AND. (k_top <= klev))
699	! find cloud top
700	IF (rneb(i,k_top) > 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