Context Navigation

lmdz_lscp_tools.f90 @ 5779

Last change on this file since 5779 was 5717, checked in by aborella, 5 months ago
Merge with trunk r5653
File size: 39.4 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,ptpronclds,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, ffallv_issr
19	use lmdz_lscp_ini, only: cice_velo, dice_velo
20	use lmdz_lscp_ini, only: ok_bug_ice_fallspeed, eps
21
22	IMPLICIT NONE
23
24	INTEGER, INTENT(IN) :: klon
25	REAL, INTENT(IN), DIMENSION(klon) :: iwc ! specific ice water content [kg/m3]
26	REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature [K]
27	REAL, INTENT(IN), DIMENSION(klon) :: rho ! dry air density [kg/m3]
28	REAL, INTENT(IN), DIMENSION(klon) :: pres ! air pressure [Pa]
29	LOGICAL, INTENT(IN), DIMENSION(klon) :: ptconv ! convective point [-]
30	LOGICAL, INTENT(IN), DIMENSION(klon) :: ptpronclds! prognostic clouds point [-]
31
32	REAL, INTENT(OUT), DIMENSION(klon) :: velo ! fallspeed velocity of crystals [m/s]
33
34
35	INTEGER i
36	REAL logvm,iwcg,tempc,phpa,fallv_tun
37	REAL m2ice, m2snow, vmice, vmsnow
38	REAL aice, bice, asnow, bsnow
39
40
41	DO i=1,klon
42
43	IF (ptconv(i)) THEN
44	fallv_tun=ffallv_con
45	ELSEIF (ptpronclds(i)) THEN
46	fallv_tun=ffallv_issr
47	ELSE
48	fallv_tun=ffallv_lsc
49	ENDIF
50
51	tempc=temp(i)-273.15 ! celcius temp
52	IF ( ok_bug_ice_fallspeed ) THEN
53	iwcg=MAX(iwc(i)*1000.,1E-3) ! iwc in g/m3. We set a minimum value to prevent from division by 0
54	ELSE
55	! AB the threshold is way too high, we reduce it
56	iwcg=MAX(iwc(i)*1000.,eps) ! iwc in g/m3. We set a minimum value to prevent from division by 0
57	ENDIF
58	phpa=pres(i)/100. ! pressure in hPa
59
60	IF (iflag_vice .EQ. 1) THEN
61	! so-called 'empirical parameterization' in Stubenrauch et al. 2019
62	if (tempc .GE. -60.0) then
63
64	logvm= -0.0000414122tempctempclog(iwcg)-0.00538922tempc*log(iwcg) &
65	-0.0516344log(iwcg)+0.00216078tempc + 1.9714
66	velo(i)=exp(logvm)
67	else
68	velo(i)=65.0(iwcg0.2)(150./phpa)**0.15
69	endif
70
71	velo(i)=fallv_tun*velo(i)/100.0 ! from cm/s to m/s
72
73	ELSE IF (iflag_vice .EQ. 2) THEN
74	! so called PSDM empirical coherent bulk ice scheme in Stubenrauch et al. 2019
75	aice=0.587
76	bice=2.45
77	asnow=0.0444
78	bsnow=2.1
79
80	m2ice=((iwcg0.001/aice)/(exp(13.6-bice7.76+0.479bice2) &
81	exp((-0.0361+bice0.0151+0.00149bice*2)tempc))) &
82	*(1./(0.807+bice0.00581+0.0457bice*2))
83
84	vmice=100.1042.4exp(13.6-(bice+1)7.76+0.479(bice+1.)*2)exp((-0.0361+ &
85	(bice+1.)0.0151+0.00149(bice+1.)*2)tempc) &
86	(m2ice(0.807+(bice+1.)0.00581+0.0457(bice+1.)2))/(iwcg0.001/aice)
87
88
89	vmice=vmice((1000./phpa)*0.2)
90
91	m2snow=((iwcg0.001/asnow)/(exp(13.6-bsnow7.76+0.479bsnow2) &
92	exp((-0.0361+bsnow0.0151+0.00149bsnow*2)tempc))) &
93	*(1./(0.807+bsnow0.00581+0.0457bsnow*2))
94
95
96	vmsnow=100.14.3exp(13.6-(bsnow+.416)7.76+0.479(bsnow+.416)**2)&
97	exp((-0.0361+(bsnow+.416)0.0151+0.00149(bsnow+.416)2)tempc)&
98	(m2snow(0.807+(bsnow+.416)0.00581+0.0457(bsnow+.416)2))/(iwcg0.001/asnow)
99
100	vmsnow=vmsnow((1000./phpa)*0.35)
101	velo(i)=fallv_tun*min(vmsnow,vmice)/100. ! to m/s
102
103	ELSE
104	! By default, fallspeed velocity of ice crystals according to Heymsfield & Donner 1990
105	velo(i) = fallv_tuncice_velo((iwcg/1000.)**dice_velo)
106	ENDIF
107	ENDDO
108
109	END SUBROUTINE FALLICE_VELOCITY
110	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
111
112	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
113	SUBROUTINE ICEFRAC_LSCP(klon, temp, iflag_ice_thermo, distcltop, temp_cltop, icefrac, dicefracdT)
114	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
115
116	! Compute the ice fraction 1-xliq (see e.g.
117	! Doutriaux-Boucher & Quaas 2004, section 2.2.)
118	! as a function of temperature
119	! see also Fig 3 of Madeleine et al. 2020, JAMES
120	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
121
122
123	USE print_control_mod, ONLY: lunout, prt_level
124	USE lmdz_lscp_ini, ONLY: t_glace_min, t_glace_max, exposant_glace, iflag_t_glace
125	USE lmdz_lscp_ini, ONLY : RTT, dist_liq, temp_nowater
126
127	IMPLICIT NONE
128
129
130	INTEGER, INTENT(IN) :: klon ! number of horizontal grid points
131	REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature
132	REAL, INTENT(IN), DIMENSION(klon) :: distcltop ! distance to cloud top
133	REAL, INTENT(IN), DIMENSION(klon) :: temp_cltop ! temperature of cloud top
134	INTEGER, INTENT(IN) :: iflag_ice_thermo
135	REAL, INTENT(OUT), DIMENSION(klon) :: icefrac
136	REAL, INTENT(OUT), DIMENSION(klon) :: dicefracdT
137
138
139	INTEGER i
140	REAL liqfrac_tmp, dicefrac_tmp
141	REAL Dv, denomdep,beta,qsi,dqsidt
142	LOGICAL ice_thermo
143
144	CHARACTER (len = 20) :: modname = 'lscp_tools'
145	CHARACTER (len = 80) :: abort_message
146
147	IF ((iflag_t_glace.LT.2)) THEN !.OR. (iflag_t_glace.GT.6)) THEN
148	abort_message = 'lscp cannot be used if iflag_t_glace<2 or >6'
149	CALL abort_physic(modname,abort_message,1)
150	ENDIF
151
152	IF (.NOT.((iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3))) THEN
153	abort_message = 'lscp cannot be used without ice thermodynamics'
154	CALL abort_physic(modname,abort_message,1)
155	ENDIF
156
157
158	DO i=1,klon
159
160	! old function with sole dependence upon temperature
161	IF (iflag_t_glace .EQ. 2) THEN
162	liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
163	liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
164	icefrac(i) = (1.0-liqfrac_tmp)**exposant_glace
165	IF (icefrac(i) .GT.0.) THEN
166	dicefracdT(i)= exposant_glace * (icefrac(i)**(exposant_glace-1.)) &
167	/ (t_glace_min - t_glace_max)
168	ENDIF
169
170	IF ((icefrac(i).EQ.0).OR.(icefrac(i).EQ.1)) THEN
171	dicefracdT(i)=0.
172	ENDIF
173
174	ENDIF
175
176	! function of temperature used in CMIP6 physics
177	IF (iflag_t_glace .EQ. 3) THEN
178	liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
179	liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
180	icefrac(i) = 1.0-liqfrac_tmp**exposant_glace
181	IF ((icefrac(i) .GT.0.) .AND. (liqfrac_tmp .GT. 0.)) THEN
182	dicefracdT(i)= exposant_glace * ((liqfrac_tmp)**(exposant_glace-1.)) &
183	/ (t_glace_min - t_glace_max)
184	ELSE
185	dicefracdT(i)=0.
186	ENDIF
187	ENDIF
188
189	! for iflag_t_glace .GE. 4, the liquid fraction depends upon temperature at cloud top
190	! and then decreases with decreasing height
191
192	!with linear function of temperature at cloud top
193	IF (iflag_t_glace .EQ. 4) THEN
194	liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
195	liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
196	icefrac(i) = MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.)
197	dicefrac_tmp = - temp(i)/(t_glace_max-t_glace_min)
198	dicefracdT(i) = dicefrac_tmp*exp(-distcltop(i)/dist_liq)
199	IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN
200	dicefracdT(i) = 0.
201	ENDIF
202	ENDIF
203
204	! with CMIP6 function of temperature at cloud top
205	IF ((iflag_t_glace .EQ. 5) .OR. (iflag_t_glace .EQ. 7)) THEN
206	liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)
207	liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
208	liqfrac_tmp = liqfrac_tmp**exposant_glace
209	icefrac(i) = MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.)
210	IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN
211	dicefracdT(i) = 0.
212	ELSE
213	dicefracdT(i) = exposant_glace((liqfrac_tmp)*(exposant_glace-1.))/(t_glace_min- t_glace_max) &
214	*exp(-distcltop(i)/dist_liq)
215	ENDIF
216	ENDIF
217
218	! with modified function of temperature at cloud top
219	! to get largere values around 260 K, works well with t_glace_min = 241K
220	IF (iflag_t_glace .EQ. 6) THEN
221	IF (temp(i) .GT. t_glace_max) THEN
222	liqfrac_tmp = 1.
223	ELSE
224	liqfrac_tmp = -((temp(i)-t_glace_max) / (t_glace_max-t_glace_min))**2+1.
225	ENDIF
226	liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
227	icefrac(i) = MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.)
228	IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN
229	dicefracdT(i) = 0.
230	ELSE
231	dicefracdT(i) = 2*((temp(i)-t_glace_max) / (t_glace_max-t_glace_min))/(t_glace_max-t_glace_min) &
232	*exp(-distcltop(i)/dist_liq)
233	ENDIF
234	ENDIF
235
236	! if temperature or temperature of cloud top <-40°C,
237	IF (iflag_t_glace .GE. 4) THEN
238	IF ((temp_cltop(i) .LE. temp_nowater) .AND. (temp(i) .LE. t_glace_max)) THEN
239	icefrac(i) = 1.
240	dicefracdT(i) = 0.
241	ENDIF
242	ENDIF
243
244
245	ENDDO ! klon
246	RETURN
247
248	END SUBROUTINE ICEFRAC_LSCP
249	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
250
251
252	SUBROUTINE ICEFRAC_LSCP_TURB(klon, dtime, pticefracturb, temp, pplay, paprsdn, paprsup, wvel, qice_ini, snowcld, snowfracld, qtot_incl, cldfra, tke, &
253	tke_dissip, sursat_e, invtau_e, qliq, qvap_cld, qice, icefrac, dicefracdT, cldfraliq, sigma2_icefracturb, mean_icefracturb)
254	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
255	! Compute the liquid, ice and vapour content (+ice fraction) based
256	! on turbulence (see Fields 2014, Furtado 2016, Raillard 2025)
257	! L.Raillard (23/09/24)
258	! E.Vignon (03/2025) : additional elements for treatment of convective
259	! boundary layer clouds
260	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
261
262
263	USE lmdz_lscp_ini, ONLY : prt_level, lunout
264	USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RV, RPI
265	USE lmdz_lscp_ini, ONLY : seuil_neb, temp_nowater
266	USE lmdz_lscp_ini, ONLY : naero5, gamma_snwretro, gamma_taud, capa_crystal, rho_ice
267	USE lmdz_lscp_ini, ONLY : eps, snow_fallspeed
268
269	IMPLICIT NONE
270
271	INTEGER, INTENT(IN) :: klon !--number of horizontal grid points
272	REAL, INTENT(IN) :: dtime !--time step [s]
273	LOGICAL, INTENT(IN), DIMENSION(klon) :: pticefracturb !--grid points concerned by this routine
274	REAL, INTENT(IN), DIMENSION(klon) :: temp !--temperature
275	REAL, INTENT(IN), DIMENSION(klon) :: pplay !--pressure in the middle of the layer [Pa]
276	REAL, INTENT(IN), DIMENSION(klon) :: paprsdn !--pressure at the bottom interface of the layer [Pa]
277	REAL, INTENT(IN), DIMENSION(klon) :: paprsup !--pressure at the top interface of the layer [Pa]
278	REAL, INTENT(IN), DIMENSION(klon) :: wvel !--vertical velocity [m/s]
279	REAL, INTENT(IN), DIMENSION(klon) :: qtot_incl !--specific total cloud water in-cloud content [kg/kg]
280	REAL, INTENT(IN), DIMENSION(klon) :: cldfra !--cloud fraction in gridbox [-]
281	REAL, INTENT(IN), DIMENSION(klon) :: tke !--turbulent kinetic energy [m2/s2]
282	REAL, INTENT(IN), DIMENSION(klon) :: tke_dissip !--TKE dissipation [m2/s3]
283
284	REAL, INTENT(IN), DIMENSION(klon) :: qice_ini !--initial specific ice content gridbox-mean [kg/kg]
285	REAL, INTENT(IN), DIMENSION(klon) :: snowcld !--in-cloud snowfall flux [kg/m2/s]
286	REAL, INTENT(IN), DIMENSION(klon) :: snowfracld !--cloudy precip fraction [-]
287	REAL, INTENT(IN), DIMENSION(klon) :: sursat_e !--environment supersaturation [-]
288	REAL, INTENT(IN), DIMENSION(klon) :: invtau_e !--inverse time-scale of mixing with environment [s-1]
289	REAL, INTENT(OUT), DIMENSION(klon) :: qliq !--specific liquid content gridbox-mean [kg/kg]
290	REAL, INTENT(OUT), DIMENSION(klon) :: qvap_cld !--specific cloud vapor content, gridbox-mean [kg/kg]
291	REAL, INTENT(OUT), DIMENSION(klon) :: qice !--specific ice content gridbox-mean [kg/kg]
292
293	REAL, INTENT(INOUT), DIMENSION(klon) :: icefrac !--fraction of ice in condensed water [-]
294	REAL, INTENT(INOUT), DIMENSION(klon) :: dicefracdT
295
296	REAL, INTENT(OUT), DIMENSION(klon) :: cldfraliq !--fraction of cldfra where liquid [-]
297	REAL, INTENT(OUT), DIMENSION(klon) :: sigma2_icefracturb !--Sigma2 of the ice supersaturation PDF [-]
298	REAL, INTENT(OUT), DIMENSION(klon) :: mean_icefracturb !--Mean of the ice supersaturation PDF [-]
299
300	REAL, DIMENSION(klon) :: qzero, qsatl, dqsatl, qsati, dqsati !--specific humidity saturation values
301	INTEGER :: i
302
303	REAL :: qvap_incl, qice_incl, qliq_incl, qiceini_incl !--In-cloud specific quantities [kg/kg]
304	REAL :: water_vapor_diff !--Water-vapour diffusion coeff in air (f(T,P)) [m2/s]
305	REAL :: air_thermal_conduct !--Thermal conductivity of air (f(T)) [J/m/K/s]
306	REAL :: C0 !--Lagrangian structure function [-]
307	REAL :: tau_dissipturb
308	REAL :: invtau_phaserelax
309	REAL :: sigma2_pdf
310	REAL :: ai, bi, B0
311	REAL :: sursat_iceliq
312	REAL :: sursat_equ
313	REAL :: liqfra_max
314	REAL :: sursat_iceext
315	REAL :: nb_crystals !--number concentration of ice crystals [#/m3]
316	REAL :: moment1_PSD !--1st moment of ice PSD
317	REAL :: N0_PSD, lambda_PSD !--parameters of the exponential PSD
318
319	REAL :: cldfra1D
320	REAL :: rho_air
321	REAL :: psati !--saturation vapor pressure wrt ice [Pa]
322
323
324	REAL :: tempvig1, tempvig2
325
326	tempvig1 = -21.06 + RTT
327	tempvig2 = -30.35 + RTT
328	C0 = 10. !--value assumed in Field2014
329	sursat_iceext = -0.1
330	qzero(:) = 0.
331	cldfraliq(:) = 0.
332	qliq(:) = 0.
333	qice(:) = 0.
334	qvap_cld(:) = 0.
335	sigma2_icefracturb(:) = 0.
336	mean_icefracturb(:) = 0.
337
338	!--wrt liquid
339	CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,1,.false.,qsatl(:),dqsatl(:))
340	!--wrt ice
341	CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,2,.false.,qsati(:),dqsati(:))
342
343
344	DO i=1,klon
345	rho_air = pplay(i) / temp(i) / RD
346	! because cldfra is intent in, but can be locally modified due to test
347	cldfra1D = cldfra(i)
348	! activate param for concerned grid points and for cloudy conditions
349	IF ((pticefracturb(i)) .AND. (cldfra(i) .GT. 0.)) THEN
350	IF (cldfra(i) .GE. 1.0) THEN
351	cldfra1D = 1.0
352	END IF
353
354	! T>0°C, no ice allowed
355	IF ( temp(i) .GE. RTT ) THEN
356	qvap_cld(i) = qsatl(i) * cldfra1D
357	qliq(i) = MAX(0.0,qtot_incl(i)-qsatl(i)) * cldfra1D
358	qice(i) = 0.
359	cldfraliq(i) = 1.
360	icefrac(i) = 0.
361	dicefracdT(i) = 0.
362
363	! T<-38°C, no liquid allowed
364	ELSE IF ( temp(i) .LE. temp_nowater) THEN
365	qvap_cld(i) = qsati(i) * cldfra1D
366	qliq(i) = 0.
367	qice(i) = MAX(0.0,qtot_incl(i)-qsati(i)) * cldfra1D
368	cldfraliq(i) = 0.
369	icefrac(i) = 1.
370	dicefracdT(i) = 0.
371
372
373	!---------------------------------------------------------
374	!-- MIXED PHASE TEMPERATURE REGIME
375	!---------------------------------------------------------
376	!--In the mixed phase regime (-38°C< T <0°C) we distinguish
377	!--3 possible subcases.
378	!--1. No pre-existing ice
379	!--2A. Pre-existing ice and no turbulence
380	!--2B. Pre-existing ice and turbulence
381	ELSE
382
383	! gamma_snwretro controls the contribution of snowflakes to the negative feedback
384	! note that for reasons related to inetarctions with the condensation iteration in lscp_main
385	! we consider here the mean snowflake concentration in the mesh (not the in-cloud concentration)
386	! when poprecip is active, it will be worth testing considering the incloud fraction, dividing
387	! by snowfracld
388	! qiceini_incl = qice_ini(i) / cldfra1D + &
389	! gamma_snwretro * snowcld(i) * RG * dtime / ( paprsdn(i) - paprsup(i) )
390	! assuming constant snowfall velocity
391	qiceini_incl = qice_ini(i) / cldfra1D + gamma_snwretro * snowcld(i) / pplay(i) * RD * temp(i) / snow_fallspeed
392
393	!--1. No preexisting ice and no mixing with environment: if vertical motion, fully liquid
394	!--cloud else fully iced cloud
395	IF ( (qiceini_incl .LT. eps) .AND. (invtau_e(i) .LT. eps) ) THEN
396	IF ( (wvel(i) .GT. eps) .OR. (tke(i) .GT. eps) ) THEN
397	qvap_cld(i) = qsatl(i) * cldfra1D
398	qliq(i) = MAX(0.,qtot_incl(i)-qsatl(i)) * cldfra1D
399	qice(i) = 0.
400	cldfraliq(i) = 1.
401	icefrac(i) = 0.
402	dicefracdT(i) = 0.
403	ELSE
404	qvap_cld(i) = qsati(i) * cldfra1D
405	qliq(i) = 0.
406	qice(i) = MAX(0.,qtot_incl(i)-qsati(i)) * cldfra1D
407	cldfraliq(i) = 0.
408	icefrac(i) = 1.
409	dicefracdT(i) = 0.
410	ENDIF
411
412
413	!--2. Pre-existing ice and/or mixing with environment:computation of ice properties for
414	!--feedback
415	ELSE
416
417	sursat_iceliq = qsatl(i)/qsati(i) - 1.
418	psati = qsati(i) * pplay(i) / (RD/RV)
419
420	!--We assume an exponential ice PSD whose parameters
421	!--are computed following Morrison&Gettelman 2008
422	!--Ice number density is assumed equals to INP density
423	!--which is for naero5>0 a function of temperature (DeMott 2010)
424	!--bi and B0 are microphysical function characterizing
425	!--vapor/ice interactions
426	!--tau_phase_relax is the typical time of vapor deposition
427	!--onto ice crystals
428
429	!--For naero5<=0 INP density is derived from the empirical fit
430	!--from MARCUS campaign from Vignon 2021
431	!--/!\ Note that option is very specific and should be use for
432	!--the Southern Ocean and the Antarctic
433
434	IF (naero5 .LE. 0) THEN
435	IF ( temp(i) .GT. tempvig1 ) THEN
436	nb_crystals = 1.e3 * 10*(-0.14(temp(i)-tempvig1) - 2.88)
437	ELSE IF ( temp(i) .GT. tempvig2 ) THEN
438	nb_crystals = 1.e3 * 10*(-0.31(temp(i)-tempvig1) - 2.88)
439	ELSE
440	nb_crystals = 1.e3 * 10**(0.)
441	ENDIF
442	ELSE
443	nb_crystals = 1.e3 * 5.94e-5 * ( RTT - temp(i) )*3.33 naero5*(0.0264(RTT-temp(i))+0.0033)
444	ENDIF
445	lambda_PSD = ( (RPIrho_icenb_crystals) / (rho_air * MAX(qiceini_incl , eps) ) ) ** (1./3.)
446	N0_PSD = nb_crystals * lambda_PSD
447	moment1_PSD = N0_PSD/lambda_PSD**2
448
449	!--Formulae for air thermal conductivity and water vapor diffusivity
450	!--comes respectively from Beard and Pruppacher (1971)
451	!--and Hall and Pruppacher (1976)
452
453	air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184
454	water_vapor_diff = 2.111e-5 ( temp(i) / RTT )*1.94 ( 101325 / pplay(i) )
455
456	bi = 1./((qsati(i)+qsatl(i))/2.) + RLSTT2 / RCPD / RV / temp(i)2
457	B0 = 4. * RPI * capa_crystal * 1. / ( RLSTT2 / air_thermal_conduct / RV / temp(i)2 &
458	+ RV * temp(i) / psati / water_vapor_diff )
459	invtau_phaserelax = bi * B0 * moment1_PSD
460
461	ai = RG / RD / temp(i) * ( RD * RLSTT / RCPD / RV / temp(i) - 1. )
462	sursat_equ = (ai * wvel(i) + sursat_e(i)*invtau_e(i)) / (invtau_phaserelax + invtau_e(i))
463	! as sursaturation is by definition lower than -1 and
464	! because local supersaturation > 1 are never found in the atmosphere
465
466	!--2A. No TKE : stationnary binary solution depending on vertical velocity and mixing with env.
467	! If Sequ > Siw liquid cloud, else ice cloud
468	IF ( tke_dissip(i) .LE. eps ) THEN
469	sigma2_icefracturb(i)= 0.
470	mean_icefracturb(i) = sursat_equ
471	IF (sursat_equ .GT. sursat_iceliq) THEN
472	qvap_cld(i) = qsatl(i) * cldfra1D
473	qliq(i) = MAX(0.,qtot_incl(i)-qsatl(i)) * cldfra1D
474	qice(i) = 0.
475	cldfraliq(i) = 1.
476	icefrac(i) = 0.
477	dicefracdT(i) = 0.
478	ELSE
479	qvap_cld(i) = qsati(i) * cldfra1D
480	qliq(i) = 0.
481	qice(i) = MAX(0.,qtot_incl(i)-qsati(i)) * cldfra1D
482	cldfraliq(i) = 0.
483	icefrac(i) = 1.
484	dicefracdT(i) = 0.
485	ENDIF
486
487	!--2B. TKE and ice : ice supersaturation PDF
488	!--we compute the cloud liquid properties with a Gaussian PDF
489	!--of ice supersaturation F(Si) (Field2014, Furtado2016).
490	!--Parameters of the PDF are function of turbulence and
491	!--microphysics/existing ice.
492	ELSE
493
494	!--Tau_dissipturb is the time needed for turbulence to decay
495	!--due to viscosity
496	tau_dissipturb = gamma_taud * 2. * 2./3. * tke(i) / tke_dissip(i) / C0
497
498	!--------------------- PDF COMPUTATIONS ---------------------
499	!--Formulae for sigma2_pdf (variance), mean of PDF in Raillard2025
500	!--cloud liquid fraction and in-cloud liquid content are given
501	!--by integrating resp. F(Si) and Si*F(Si)
502	!--Liquid is limited by the available water vapor trough a
503	!--maximal liquid fraction
504	!--qice_ini(i) / cldfra1D = qiceincld without precip
505
506	liqfra_max = MAX(0., (MIN (1.,( qtot_incl(i) - (qice_ini(i) / cldfra1D) - qsati(i) * (1 + sursat_iceext ) ) / ( qsatl(i) - qsati(i) ) ) ) )
507	sigma2_pdf = 1./2. * ( ai*2 ) 2./3. * tke(i) * tau_dissipturb / (invtau_phaserelax + invtau_e(i))
508	! sursat ranges between -1 and 1, so we prevent sigma2 so exceed 1
509	cldfraliq(i) = 0.5 * (1. - erf( ( sursat_iceliq - sursat_equ) / (SQRT(2.* sigma2_pdf) ) ) )
510	IF (cldfraliq(i) .GT. liqfra_max) THEN
511	cldfraliq(i) = liqfra_max
512	ENDIF
513
514	qliq_incl = qsati(i) * SQRT(sigma2_pdf) / SQRT(2.RPI) EXP( -1.(sursat_iceliq - sursat_equ)2. / (2.sigma2_pdf) ) &
515	- qsati(i) * cldfraliq(i) * (sursat_iceliq - sursat_equ )
516
517	sigma2_icefracturb(i)= sigma2_pdf
518	mean_icefracturb(i) = sursat_equ
519
520	!------------ SPECIFIC VAPOR CONTENT AND WATER CONSERVATION ------------
521
522	IF ( (qliq_incl .LE. eps) .OR. (cldfraliq(i) .LE. eps) ) THEN
523	qliq_incl = 0.
524	cldfraliq(i) = 0.
525	END IF
526
527	!--Specific humidity is the max between qsati and the weighted mean between
528	!--qv in MPC patches and qv in ice-only parts. We assume that MPC parts are
529	!--always at qsatl and ice-only parts slightly subsaturated (qsati*sursat_iceext+1)
530	!--The whole cloud can therefore be supersaturated but never subsaturated.
531
532	qvap_incl = MAX(qsati(i), ( 1. - cldfraliq(i) ) * (sursat_iceext + 1.) * qsati(i) + cldfraliq(i) * qsatl(i) )
533
534	IF ( qvap_incl .GE. qtot_incl(i) ) THEN
535	qvap_incl = qsati(i)
536	qliq_incl = MAX(0.0,qtot_incl(i) - qvap_incl)
537	qice_incl = 0.
538
539	ELSEIF ( (qvap_incl + qliq_incl) .GE. qtot_incl(i) ) THEN
540	qliq_incl = MAX(0.0,qtot_incl(i) - qvap_incl)
541	qice_incl = 0.
542	ELSE
543	qice_incl = qtot_incl(i) - qvap_incl - qliq_incl
544	END IF
545
546	qvap_cld(i) = qvap_incl * cldfra1D
547	qliq(i) = qliq_incl * cldfra1D
548	qice(i) = qice_incl * cldfra1D
549	IF ((qice(i)+qliq(i)) .GT. 0.) THEN
550	icefrac(i) = qice(i) / ( qice(i) + qliq(i) )
551	ELSE
552	icefrac(i) = 1. ! to keep computation of qsat wrt ice in condensation loop in lmdz_lscp_main
553	ENDIF
554	dicefracdT(i) = 0.
555
556	END IF ! Enough TKE
557
558	END IF ! End qini
559
560	END IF ! ! MPC temperature
561
562	END IF ! pticefracturb and cldfra
563
564	ENDDO ! klon
565	END SUBROUTINE ICEFRAC_LSCP_TURB
566	!
567	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
568
569
570	SUBROUTINE CALC_QSAT_ECMWF(klon,temp,qtot,pressure,tref,phase,flagth,qs,dqs)
571	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
572	! Calculate qsat following ECMWF method
573	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
574
575
576	USE yoethf_mod_h
577	USE yomcst_mod_h
578	IMPLICIT NONE
579
580
581	include "FCTTRE.h"
582
583	INTEGER, INTENT(IN) :: klon ! number of horizontal grid points
584	REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature in K
585	REAL, INTENT(IN), DIMENSION(klon) :: qtot ! total specific water in kg/kg
586	REAL, INTENT(IN), DIMENSION(klon) :: pressure ! pressure in Pa
587	REAL, INTENT(IN) :: tref ! reference temperature in K
588	LOGICAL, INTENT(IN) :: flagth ! flag for qsat calculation for thermals
589	INTEGER, INTENT(IN) :: phase
590	! phase: 0=depend on temperature sign (temp>tref -> liquid, temp<tref, solid)
591	! 1=liquid
592	! 2=solid
593
594	REAL, INTENT(OUT), DIMENSION(klon) :: qs ! saturation specific humidity [kg/kg]
595	REAL, INTENT(OUT), DIMENSION(klon) :: dqs ! derivation of saturation specific humidity wrt T
596
597	REAL delta, cor, cvm5
598	INTEGER i
599
600	DO i=1,klon
601
602	IF (phase .EQ. 1) THEN
603	delta=0.
604	ELSEIF (phase .EQ. 2) THEN
605	delta=1.
606	ELSE
607	delta=MAX(0.,SIGN(1.,tref-temp(i)))
608	ENDIF
609
610	IF (flagth) THEN
611	cvm5=R5LES(1.-delta) + R5IESdelta
612	ELSE
613	cvm5 = R5LESRLVTT(1.-delta) + R5IESRLSTTdelta
614	cvm5 = cvm5 /RCPD/(1.0+RVTMP2*(qtot(i)))
615	ENDIF
616
617	qs(i)= R2ES*FOEEW(temp(i),delta)/pressure(i)
618	qs(i)=MIN(0.5,qs(i))
619	cor=1./(1.-RETV*qs(i))
620	qs(i)=qs(i)*cor
621	dqs(i)= FOEDE(temp(i),delta,cvm5,qs(i),cor)
622
623	END DO
624
625	END SUBROUTINE CALC_QSAT_ECMWF
626	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
627
628
629	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
630	SUBROUTINE CALC_GAMMASAT(klon,temp,qtot,pressure,gammasat,dgammasatdt)
631
632	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
633	! programme that calculates the gammasat parameter that determines the
634	! homogeneous condensation thresholds for cold (<0oC) clouds
635	! condensation at q>gammasat*qsat
636	! Etienne Vignon, March 2021
637	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
638
639	use lmdz_lscp_ini, only: iflag_gammasat, temp_nowater, RTT
640	use lmdz_lscp_ini, only: a_homofreez, b_homofreez, delta_hetfreez
641
642	IMPLICIT NONE
643
644
645	INTEGER, INTENT(IN) :: klon ! number of horizontal grid points
646	REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature in K
647	REAL, INTENT(IN), DIMENSION(klon) :: qtot ! total specific water in kg/kg
648
649	REAL, INTENT(IN), DIMENSION(klon) :: pressure ! pressure in Pa
650
651	REAL, INTENT(OUT), DIMENSION(klon) :: gammasat ! coefficient to multiply qsat with to calculate saturation
652	REAL, INTENT(OUT), DIMENSION(klon) :: dgammasatdt ! derivative of gammasat wrt temperature
653
654	REAL, DIMENSION(klon) :: qsi,qsl,dqsl,dqsi
655	REAL f_homofreez, fac
656
657	INTEGER i
658
659	CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,1,.false.,qsl,dqsl)
660	CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,2,.false.,qsi,dqsi)
661
662	DO i = 1, klon
663
664	IF ( temp(i) .GE. RTT ) THEN
665	! warm clouds: condensation at saturation wrt liquid
666	gammasat(i) = 1.
667	dgammasatdt(i) = 0.
668
669	ELSE
670	! cold clouds: qsi > qsl
671
672	! homogeneous freezing of aerosols, according to
673	! Koop, 2000 and Ren and MacKenzie, 2005 (QJRMS)
674	! 'Cirrus regime'
675	! if f_homofreez > qsl / qsi, liquid nucleation
676	! if f_homofreez < qsl / qsi, homogeneous freezing of aerosols
677	! Note: f_homofreez = qsl / qsi for temp ~= -38degC
678	f_homofreez = a_homofreez - temp(i) / b_homofreez
679
680	IF ( iflag_gammasat .GE. 3 ) THEN
681	! condensation at homogeneous freezing threshold for temp < -38 degC
682	! condensation at liquid saturation for temp > -38 degC
683	IF ( f_homofreez .LE. qsl(i) / qsi(i) ) THEN
684	gammasat(i) = f_homofreez
685	dgammasatdt(i) = - 1. / b_homofreez
686	ELSE
687	gammasat(i) = qsl(i) / qsi(i)
688	dgammasatdt(i) = ( dqsl(i) * qsi(i) - dqsi(i) * qsl(i) ) / qsi(i) / qsi(i)
689	ENDIF
690
691	ELSEIF ( iflag_gammasat .EQ. 2 ) THEN
692	! condensation at homogeneous freezing threshold for temp < -38 degC
693	! condensation at a threshold linearly decreasing between homogeneous
694	! freezing and ice saturation for -38 degC < temp < temp_nowater
695	! condensation at ice saturation for temp > temp_nowater
696	! If temp_nowater = 235.15 K, this is equivalent to iflag_gammasat = 1
697	IF ( f_homofreez .LE. qsl(i) / qsi(i) ) THEN
698	gammasat(i) = f_homofreez
699	dgammasatdt(i) = - 1. / b_homofreez
700	ELSEIF ( temp(i) .LE. temp_nowater ) THEN
701	! Here, we assume that f_homofreez = qsl / qsi for temp = -38 degC = 235.15 K
702	dgammasatdt(i) = ( a_homofreez - 235.15 / b_homofreez - 1. ) &
703	/ ( 235.15 - temp_nowater )
704	gammasat(i) = dgammasatdt(i) * ( temp(i) - temp_nowater ) + 1.
705	ELSE
706	gammasat(i) = 1.
707	dgammasatdt(i) = 0.
708	ENDIF
709
710	ELSEIF ( iflag_gammasat .EQ. 1 ) THEN
711	! condensation at homogeneous freezing threshold for temp < -38 degC
712	! condensation at ice saturation for temp > -38 degC
713	IF ( f_homofreez .LE. qsl(i) / qsi(i) ) THEN
714	gammasat(i) = f_homofreez
715	dgammasatdt(i) = - 1. / b_homofreez
716	ELSE
717	gammasat(i) = 1.
718	dgammasatdt(i) = 0.
719	ENDIF
720
721	ELSE
722	! condensation at ice saturation for temp < -38 degC
723	! condensation at ice saturation for temp > -38 degC
724	gammasat(i) = 1.
725	dgammasatdt(i) = 0.
726
727	ENDIF
728
729	! Note that the delta_hetfreez parameter allows to linearly decrease the
730	! condensation threshold between the calculated threshold and the ice saturation
731	! for delta_hetfreez = 1, the threshold is the calculated condensation threshold
732	! for delta_hetfreez = 0, the threshold is the ice saturation
733	gammasat(i) = ( 1. - delta_hetfreez ) + delta_hetfreez * gammasat(i)
734	dgammasatdt(i) = delta_hetfreez * dgammasatdt(i)
735
736	ENDIF
737
738	END DO
739
740
741	END SUBROUTINE CALC_GAMMASAT
742	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
743
744
745	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
746	SUBROUTINE DISTANCE_TO_CLOUD_TOP(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D,temp_cltop)
747	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
748
749	USE lmdz_lscp_ini, ONLY : rd,rg,tresh_cl
750
751	IMPLICIT NONE
752
753	INTEGER, INTENT(IN) :: klon,klev !number of horizontal and vertical grid points
754	INTEGER, INTENT(IN) :: k ! vertical index
755	REAL, INTENT(IN), DIMENSION(klon,klev) :: temp ! temperature in K
756	REAL, INTENT(IN), DIMENSION(klon,klev) :: pplay ! pressure middle layer in Pa
757	REAL, INTENT(IN), DIMENSION(klon,klev+1) :: paprs ! pressure interfaces in Pa
758	REAL, INTENT(IN), DIMENSION(klon,klev) :: rneb ! cloud fraction
759
760	REAL, INTENT(OUT), DIMENSION(klon) :: distcltop1D ! distance from cloud top
761	REAL, INTENT(OUT), DIMENSION(klon) :: temp_cltop ! temperature of cloud top
762
763	REAL dzlay(klon,klev)
764	REAL zlay(klon,klev)
765	REAL dzinterf
766	INTEGER i,k_top, kvert
767	LOGICAL bool_cl
768
769
770	DO i=1,klon
771	! Initialization height middle of first layer
772	dzlay(i,1) = Rd * temp(i,1) / rg * log(paprs(i,1)/paprs(i,2))
773	zlay(i,1) = dzlay(i,1)/2
774
775	DO kvert=2,klev
776	IF (kvert.EQ.klev) THEN
777	dzlay(i,kvert) = 2(rd temp(i,kvert) / rg * log(paprs(i,kvert)/pplay(i,kvert)))
778	ELSE
779	dzlay(i,kvert) = rd * temp(i,kvert) / rg * log(paprs(i,kvert)/paprs(i,kvert+1))
780	ENDIF
781	dzinterf = rd * temp(i,kvert) / rg * log(pplay(i,kvert-1)/pplay(i,kvert))
782	zlay(i,kvert) = zlay(i,kvert-1) + dzinterf
783	ENDDO
784	ENDDO
785
786	DO i=1,klon
787	k_top = k
788	IF (rneb(i,k) .LE. tresh_cl) THEN
789	bool_cl = .FALSE.
790	ELSE
791	bool_cl = .TRUE.
792	ENDIF
793
794	DO WHILE ((bool_cl) .AND. (k_top .LE. klev))
795	! find cloud top
796	IF (rneb(i,k_top) .GT. tresh_cl) THEN
797	k_top = k_top + 1
798	ELSE
799	bool_cl = .FALSE.
800	k_top = k_top - 1
801	ENDIF
802	ENDDO
803	k_top=min(k_top,klev)
804
805	!dist to top is dist between current layer and layer of cloud top (from middle to middle) + dist middle to
806	!interf for layer of cloud top
807	distcltop1D(i) = zlay(i,k_top) - zlay(i,k) + dzlay(i,k_top)/2
808	temp_cltop(i) = temp(i,k_top)
809	ENDDO ! klon
810
811	END SUBROUTINE DISTANCE_TO_CLOUD_TOP
812	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
813
814	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
815	FUNCTION GAMMAINC ( p, x )
816
817	!*****************************************************************************80
818	!
819	!! GAMMAINC computes the regularized lower incomplete Gamma Integral
820	!
821	! Modified:
822	!
823	! 20 January 2008
824	!
825	! Author:
826	!
827	! Original FORTRAN77 version by B Shea.
828	! FORTRAN90 version by John Burkardt.
829	!
830	! Reference:
831	!
832	! B Shea,
833	! Algorithm AS 239:
834	! Chi-squared and Incomplete Gamma Integral,
835	! Applied Statistics,
836	! Volume 37, Number 3, 1988, pages 466-473.
837	!
838	! Parameters:
839	!
840	! Input, real X, P, the parameters of the incomplete
841	! gamma ratio. 0 <= X, and 0 < P.
842	!
843	! Output, real GAMMAINC, the value of the incomplete
844	! Gamma integral.
845	!
846	IMPLICIT NONE
847
848	REAL A
849	REAL AN
850	REAL ARG
851	REAL B
852	REAL C
853	REAL, PARAMETER :: ELIMIT = - 88.0E+00
854	REAL GAMMAINC
855	REAL, PARAMETER :: OFLO = 1.0E+37
856	REAL P
857	REAL, PARAMETER :: PLIMIT = 1000.0E+00
858	REAL PN1
859	REAL PN2
860	REAL PN3
861	REAL PN4
862	REAL PN5
863	REAL PN6
864	REAL RN
865	REAL, PARAMETER :: TOL = 1.0E-14
866	REAL X
867	REAL, PARAMETER :: XBIG = 1.0E+08
868
869	GAMMAINC = 0.0E+00
870
871	IF ( X == 0.0E+00 ) THEN
872	GAMMAINC = 0.0E+00
873	RETURN
874	END IF
875	!
876	! IF P IS LARGE, USE A NORMAL APPROXIMATION.
877	!
878	IF ( PLIMIT < P ) THEN
879
880	PN1 = 3.0E+00 * SQRT ( P ) * ( ( X / P )**( 1.0E+00 / 3.0E+00 ) &
881	+ 1.0E+00 / ( 9.0E+00 * P ) - 1.0E+00 )
882
883	GAMMAINC = 0.5E+00 * ( 1. + ERF ( PN1 ) )
884	RETURN
885
886	END IF
887	!
888	! IF X IS LARGE SET GAMMAD = 1.
889	!
890	IF ( XBIG < X ) THEN
891	GAMMAINC = 1.0E+00
892	RETURN
893	END IF
894	!
895	! USE PEARSON'S SERIES EXPANSION.
896	! (NOTE THAT P IS NOT LARGE ENOUGH TO FORCE OVERFLOW IN ALOGAM).
897	!
898	IF ( X <= 1.0E+00 .OR. X < P ) THEN
899
900	ARG = P * LOG ( X ) - X - LOG_GAMMA ( P + 1.0E+00 )
901	C = 1.0E+00
902	GAMMAINC = 1.0E+00
903	A = P
904
905	DO
906
907	A = A + 1.0E+00
908	C = C * X / A
909	GAMMAINC = GAMMAINC + C
910
911	IF ( C <= TOL ) THEN
912	EXIT
913	END IF
914
915	END DO
916
917	ARG = ARG + LOG ( GAMMAINC )
918
919	IF ( ELIMIT <= ARG ) THEN
920	GAMMAINC = EXP ( ARG )
921	ELSE
922	GAMMAINC = 0.0E+00
923	END IF
924	!
925	! USE A CONTINUED FRACTION EXPANSION.
926	!
927	ELSE
928
929	ARG = P * LOG ( X ) - X - LOG_GAMMA ( P )
930	A = 1.0E+00 - P
931	B = A + X + 1.0E+00
932	C = 0.0E+00
933	PN1 = 1.0E+00
934	PN2 = X
935	PN3 = X + 1.0E+00
936	PN4 = X * B
937	GAMMAINC = PN3 / PN4
938
939	DO
940
941	A = A + 1.0E+00
942	B = B + 2.0E+00
943	C = C + 1.0E+00
944	AN = A * C
945	PN5 = B * PN3 - AN * PN1
946	PN6 = B * PN4 - AN * PN2
947
948	IF ( PN6 /= 0.0E+00 ) THEN
949
950	RN = PN5 / PN6
951
952	IF ( ABS ( GAMMAINC - RN ) <= MIN ( TOL, TOL * RN ) ) THEN
953	EXIT
954	END IF
955
956	GAMMAINC = RN
957
958	END IF
959
960	PN1 = PN3
961	PN2 = PN4
962	PN3 = PN5
963	PN4 = PN6
964	!
965	! RE-SCALE TERMS IN CONTINUED FRACTION IF TERMS ARE LARGE.
966	!
967	IF ( OFLO <= ABS ( PN5 ) ) THEN
968	PN1 = PN1 / OFLO
969	PN2 = PN2 / OFLO
970	PN3 = PN3 / OFLO
971	PN4 = PN4 / OFLO
972	END IF
973
974	END DO
975
976	ARG = ARG + LOG ( GAMMAINC )
977
978	IF ( ELIMIT <= ARG ) THEN
979	GAMMAINC = 1.0E+00 - EXP ( ARG )
980	ELSE
981	GAMMAINC = 1.0E+00
982	END IF
983
984	END IF
985
986	RETURN
987	END FUNCTION GAMMAINC
988	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
989
990	END MODULE lmdz_lscp_tools
991
992

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

Download in other formats:

Original Format