Context Navigation

mp2m_clouds.F90 @ 4096

Last change on this file since 4096 was 4096, checked in by debatzbr, 2 weeks ago
Pluto PCM: ccn moment first layer correction (from revision 4084 - Titan PCM, CP). BBT
File size: 56.2 KB

Line
1	MODULE MP2M_CLOUDS
2	!============================================================================
3	!
4	! Purpose
5	! -------
6	! Clouds microphysics module.
7	!
8	! This module contains all definitions of the microphysics processes related to
9	! clouds (nucleation, condensation, sedimentation).
10	!
11	! The interface methods always use the global variables defined in [[mm_globals(module)]] when values
12	! (temperature, pressure, moments...) over the vertical grid are required.
13	! Consequently, all these functions only deal with output arguments which are most of the time the
14	! tendencies of relevant variables on the atmospheric column.
15	!
16	! @Warning
17	! The tendencies returned by the method are always defined over the vertical grid
18	! from __TOP__ to __GROUND__.
19	!
20	! The module also contains ten methods:
21	! - mm_cloud_microphysics \| Evolution of moments tracers through clouds microphysics processes.
22	! - mm_cloud_nucond \| Get moments tendencies through nucleation/condensation/evaporation.
23	! * nc_esp \| Get moments tendencies through nucleation/condensation/evaporation of a given condensible species.
24	! * hetnucl_rate_aer \| Get heterogeneous nucleation rate on aerosols.
25	! * growth_rate \| Get growth rate through condensation/evaporation process.
26	! - mm_cloud_sedimentation \| Compute the tendency of clouds related moments through sedimentation process.
27	! * exchange \| Compute the exchange matrix.
28	! * getnzs \| Compute displacement of a cell under sedimentation process.
29	! * wsettle \| Compute the settling velocity of a spherical particle.
30	! * get_mass_flux \| Get the mass flux of clouds related moment through sedimention.
31	! - mm_cloud_coagulation \| Compute the tendency of clouds related moments through coagulation process.
32	!
33	! Authors
34	! -------
35	! B. de Batz de Trenquelléon (10/2025)
36	!
37	!============================================================================
38
39	USE MP2M_MPREC
40	USE MP2M_GLOBALS
41	USE MP2M_METHODS
42	USE MP2M_CLOUDS_METHODS
43	IMPLICIT NONE
44
45	PRIVATE
46
47	PUBLIC :: mm_cloud_microphysics, mm_cloud_nucond, mm_cloud_sedimentation
48
49	CONTAINS
50
51	!============================================================================
52	! CLOUDS MICROPHYSICS INTERFACE SUBROUTINE
53	!============================================================================
54
55	SUBROUTINE mm_cloud_microphysics(Hdm0as,Hdm3as,Hdm0af,Hdm3af,&
56	Cdm0as,Cdm3as,Cdm0af,Cdm3af,&
57	Cdm0ccn,Cdm3ccn,Cdm3ice,Cdmugas)
58	!! Get the evolution of moments tracers through cloud microphysics processes.
59	!! The subroutine is a wrapper to the cloud microphysics methods. It computes the tendencies of moments
60	!! tracers for nucleation, condensation and sedimentation processes for the atmospheric column.
61	!!
62	!! @note
63	!! Both __dm3ice__ and __dmugas__ are 2D-array with the vertical layers in first dimension and the number
64	!! of ice components in the second.
65	!!
66	! Tendency of the 0th and 3rd order moment of the aerosols (spherical mode) through haze microphysics processes.
67	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: Hdm0as ! (m-3)
68	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: Hdm3as ! (m3.m-3)
69	! Tendency of the 0th and 3rd order moment of the aerosols (fractal mode) through haze microphysics processes.
70	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: Hdm0af ! (m-3)
71	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: Hdm3af ! (m3.m-3)
72	! Tendency of the 0th and 3rd order moment of the aerosols (spherical mode) through cloud microphysics processes.
73	REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: Cdm0as ! (m-3)
74	REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: Cdm3as ! (m3.m-3)
75	! Tendency of the 0th and 3rd order moment of the aerosols (fractal mode) through cloud microphysics processes.
76	REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: Cdm0af ! (m-3)
77	REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: Cdm3af ! (m3.m-3)
78	! Tendency of the 0th and 3rd order moment of the ccn distribution through cloud microphysics processes.
79	REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: Cdm0ccn ! (m-3)
80	REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: Cdm3ccn ! (m3.m-3)
81	! Tendencies of the 3rd order moments of each ice components through cloud microphysics processes.
82	REAL(kind=mm_wp), DIMENSION(:,:), INTENT(inout) :: Cdm3ice ! (m3.m-3)
83	! Tendencies of each condensible gaz species through cloud microphysics processes.
84	REAL(kind=mm_wp), DIMENSION(:,:), INTENT(inout) :: Cdmugas ! (mol.mol-1)
85
86	! Local variables
87	!~~~~~~~~~~~~~~~~
88	INTEGER :: i
89	! Updated tracers with haze and cloud related processes tendencies
90	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: m0as,m3as
91	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: m0af,m3af
92	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: m0ccn,m3ccn
93	REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE :: m3ice,mugas
94	! Temporary tendencies through nucleation/condensation/sublimation
95	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: nCdm0as,nCdm3as
96	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: nCdm0af,nCdm3af
97	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: nCdm0ccn,nCdm3ccn
98	REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE :: nCdm3ice,nCdmugas
99	! Temporary tendencies through coagulation
100	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: cCdm0ccn,cCdm3ccn
101	! Temporary tendencies through sedimentation (usefull for diagnostics)
102	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: zCdm0ccn,zCdm3ccn
103	REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE :: zCdm3ice
104
105	! Initialization:
106	!~~~~~~~~~~~~~~~~
107	ALLOCATE(m0as(mm_nla),m3as(mm_nla),m0af(mm_nla),m3af(mm_nla))
108	ALLOCATE(m0ccn(mm_nla),m3ccn(mm_nla),m3ice(mm_nla,mm_nesp),mugas(mm_nla,mm_nesp))
109	ALLOCATE(nCdm0as(mm_nla),nCdm3as(mm_nla),nCdm0af(mm_nla),nCdm3af(mm_nla))
110	ALLOCATE(nCdm0ccn(mm_nla),nCdm3ccn(mm_nla),nCdm3ice(mm_nla,mm_nesp),nCdmugas(mm_nla,mm_nesp))
111	ALLOCATE(cCdm0ccn(mm_nla),cCdm3ccn(mm_nla))
112	ALLOCATE(zCdm0ccn(mm_nla),zCdm3ccn(mm_nla),zCdm3ice(mm_nla,mm_nesp))
113
114	m0as(:) = mm_m0aer_s + Hdm0as ; m3as(:) = mm_m3aer_s + Hdm3as
115	m0af(:) = mm_m0aer_f + Hdm0af ; m3af(:) = mm_m3aer_f + Hdm3af
116	m0ccn(:) = mm_m0ccn ; m3ccn(:) = mm_m3ccn
117	m3ice(:,:) = mm_m3ice ; mugas(:,:) = mm_gas
118
119	Cdm0as(:) = 0._mm_wp ; Cdm3as(:) = 0._mm_wp
120	Cdm0af(:) = 0._mm_wp ; Cdm3af(:) = 0._mm_wp
121	Cdm0ccn(:) = 0._mm_wp ; Cdm3ccn(:) = 0._mm_wp
122	Cdm3ice(:,:) = 0._mm_wp ; Cdmugas(:,:) = 0._mm_wp
123
124	nCdm0as(:) = 0._mm_wp ; nCdm3as(:) = 0._mm_wp
125	nCdm0af(:) = 0._mm_wp ; nCdm3af(:) = 0._mm_wp
126	nCdm0ccn(:) = 0._mm_wp ; nCdm3ccn(:) = 0._mm_wp
127	nCdm3ice(:,:) = 0._mm_wp ; nCdmugas(:,:) = 0._mm_wp
128	cCdm0ccn(:) = 0._mm_wp ; cCdm3ccn(:) = 0._mm_wp
129	zCdm0ccn(:) = 0._mm_wp ; zCdm3ccn(:) = 0._mm_wp
130	zCdm3ice(:,:) = 0._mm_wp
131
132	! Calls nucleation/condensation:
133	! And update saturation ratio, nucleation rate and growth rate diagnostic.
134	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
135	call mm_cloud_nucond(m0as,m3as,m0af,m3af,m0ccn,m3ccn,m3ice,mugas,&
136	nCdm0as,nCdm3as,nCdm0af,nCdm3af,&
137	nCdm0ccn,nCdm3ccn,nCdm3ice,nCdmugas,mm_gas_sat,mm_nrate,mm_grate)
138
139	! Update tracers
140	m0as = m0as + nCdm0as ; m3as = m3as + nCdm3as
141	m0af = m0af + nCdm0af ; m3af = m3af + nCdm3af
142	m0ccn = m0ccn + nCdm0ccn ; m3ccn = m3ccn + nCdm3ccn
143	m3ice = m3ice + nCdm3ice ; mugas = mugas + nCdmugas
144
145	where(m3as > 0._mm_wp .AND. m0as > 0._mm_wp)
146	mm_rcs = mm_get_rcs(m0as,m3as)
147	elsewhere
148	mm_rcs = 0._mm_wp
149	endwhere
150	where(m3af > 0._mm_wp .AND. m0af > 0._mm_wp)
151	mm_rcf = mm_get_rcf(m0af,m3af)
152	elsewhere
153	mm_rcf = 0._mm_wp
154	endwhere
155	call mm_cloud_properties(m0ccn,m3ccn,m3ice,mm_drad,mm_drho)
156
157	! Calls coagulation:
158	!~~~~~~~~~~~~~~~~~~~
159	call mm_cloud_coagulation(m0ccn,m3ccn,m3ice,cCdm0ccn,cCdm3ccn)
160
161	! Update tracers
162	m0ccn = m0ccn + cCdm0ccn ; m3ccn = m3ccn + cCdm3ccn
163	call mm_cloud_properties(m0ccn,m3ccn,m3ice,mm_drad,mm_drho)
164
165	! Calls sedimentation:
166	!~~~~~~~~~~~~~~~~~~~~~
167	call mm_cloud_sedimentation(m0ccn,m3ccn,m3ice,zCdm0ccn,zCdm3ccn,zCdm3ice)
168
169	! Update tracers
170	m0ccn = m0ccn + zCdm0ccn ; m3ccn = m3ccn + zCdm3ccn
171	m3ice = m3ice + zCdm3ice
172	call mm_cloud_properties(m0ccn,m3ccn,m3ice,mm_drad,mm_drho)
173
174	! Computes diagnostics:
175	!~~~~~~~~~~~~~~~~~~~~~~
176	! Settling velocity [m.s-1] of clouds (ccn and ices)
177	mm_cld_vsed(:) = wsettle(mm_play,mm_temp,mm_zlay,mm_rhoair,mm_drho,mm_drad)
178
179	! Flux [kg.m-2.s-1] and precipitation [kg.m-2] of ccn
180	mm_ccn_flux(:) = get_mass_flux(mm_rhoaer,m3ccn(:))
181	mm_ccn_prec = SUM(zCdm3ccn(:)mm_dzlevmm_rhoaer)
182
183	! Flux [kg.m-2.s-1] and precipitation [kg.m-2] of ices
184	DO i = 1, mm_nesp
185	mm_ice_fluxes(:,i) = get_mass_flux(mm_xESPS(i)%rho_s,m3ice(:,i))
186	mm_ice_prec(i) = SUM(zCdm3ice(:,i)mm_dzlevmm_xESPS(i)%rho_s)
187	ENDDO
188
189	! Updates tendencies:
190	!~~~~~~~~~~~~~~~~~~~~
191	Cdm0as = nCdm0as ; Cdm3as = nCdm3as
192	Cdm0af = nCdm0af ; Cdm3af = nCdm3af
193	Cdm0ccn = nCdm0ccn + cCdm0ccn + zCdm0ccn ; Cdm3ccn = nCdm3ccn + cCdm3ccn + zCdm3ccn
194	Cdm3ice = nCdm3ice + zCdm3ice ; Cdmugas = nCdmugas
195
196	END SUBROUTINE mm_cloud_microphysics
197
198	!============================================================================
199	! NUCLEATION/CONDENSATION PROCESS RELATED METHODS
200	!============================================================================
201
202	SUBROUTINE mm_cloud_nucond(m0as,m3as,m0af,m3af,m0ccn,m3ccn,m3ice,mugas,&
203	dm0s,dm3s,dm0f,dm3f,dm0n,dm3n,dm3i,dmugas,&
204	gassat,nrate,grate)
205	!! Get moments tendencies through nucleation and condensation/evaporation.
206	!!
207	!! The method is a wrapper of [[mm_clouds(module):nc_esp(subroutine)]] which computes the
208	!! tendencies of tracers for all the condensible species given in the vector __xESPS__.
209	!!
210	!! Aerosols and CCN distribution evolution depends on the ice components:
211	!! - For nucleation only creation of CCN can occur.
212	!! - For condensation/evaporation only loss of CCN can occur.
213	!!
214	!! @note
215	!! We use the simple following rule to compute the variation of CCN and aerosols:
216	!! The global variation of CCN (and thus aerosols) is determined from the most intense activity
217	!! of the ice components.
218	!!
219	!! @warning
220	!! __xESPS__, __m3i__ and __gazes__ must share the same indexing.
221	!! For example if xESPS(i) corresponds to CH4 properties then m3i(i) must be
222	!! the total volume of CH4 ice and gazs(i) its vapor mole fraction.
223	!!
224	! 0th and 3rd order moment of the aerosols (spherical mode) component.
225	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m0as ! (m-3)
226	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m3as ! (m3.m-3)
227	! 0th and 3rd order moment of the aerosols (fractal mode) component.
228	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m0af ! (m-3)
229	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m3af ! (m3.m-3)
230	! 0th and 3rd order moment of the cloud condensation nuclei component.
231	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m0ccn ! (m-3)
232	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m3ccn ! (m3.m-3)
233	! 3rd order moment of ices and gases component.
234	REAL(kind=mm_wp), DIMENSION(:,:), INTENT(in) :: m3ice ! (m-3)
235	REAL(kind=mm_wp), DIMENSION(:,:), INTENT(in) :: mugas ! (m3.m-3)
236	! Tendency of the 0th and 3rd order moment of the aerosols (spherical mode) through cloud microphysics processes.
237	REAL(kind=mm_wp), DIMENSION(:), INTENT(out) :: dm0s ! (m-3)
238	REAL(kind=mm_wp), DIMENSION(:), INTENT(out) :: dm3s ! (m3.m-3)
239	! Tendency of the 0th and 3rd order moment of the aerosols (fractal mode) through cloud microphysics processes.
240	REAL(kind=mm_wp), DIMENSION(:), INTENT(out) :: dm0f ! (m-3)
241	REAL(kind=mm_wp), DIMENSION(:), INTENT(out) :: dm3f ! (m3.m-3)
242	! Tendency of the 0th and 3rd order moment of the ccn distribution through cloud microphysics processes.
243	REAL(kind=mm_wp), DIMENSION(:), INTENT(out) :: dm0n ! (m-3)
244	REAL(kind=mm_wp), DIMENSION(:), INTENT(out) :: dm3n ! (m3.m-3)
245	! Tendencies of the 3rd order moments of each ice components through cloud microphysics processes.
246	REAL(kind=mm_wp), DIMENSION(:,:), INTENT(out) :: dm3i ! (m3.m-3)
247	! Tendencies of each condensible gaz species through cloud microphysics processes.
248	REAL(kind=mm_wp), DIMENSION(:,:), INTENT(out) :: dmugas ! (mol.mol-1)
249	! Saturation ratio of each condensible species.
250	REAL(kind=mm_wp), DIMENSION(:,:), INTENT(out) :: gassat
251	! Nucleation rate values of each condensible species (m-2.s-1).
252	REAL(kind=mm_wp), DIMENSION(:,:),INTENT(out) :: nrate
253	! Growth rate values of each condensible species (m2.s-1).
254	REAL(kind=mm_wp), DIMENSION(:,:),INTENT(out) :: grate
255
256	! Local variables:
257	!~~~~~~~~~~~~~~~~~
258	INTEGER :: i,idx
259	REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE :: zdm0s,zdm3s
260	REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE :: zdm0f,zdm3f
261	REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE :: zdm0n,zdm3n
262
263	! Initialization:
264	!~~~~~~~~~~~~~~~~
265	ALLOCATE(zdm0s(mm_nla,mm_nesp),zdm3s(mm_nla,mm_nesp))
266	ALLOCATE(zdm0f(mm_nla,mm_nesp),zdm3f(mm_nla,mm_nesp))
267	ALLOCATE(zdm0n(mm_nla,mm_nesp),zdm3n(mm_nla,mm_nesp))
268
269	zdm0s(:,:) = 0._mm_wp ; zdm3s(:,:) = 0._mm_wp
270	zdm0f(:,:) = 0._mm_wp ; zdm3f(:,:) = 0._mm_wp
271	zdm0n(:,:) = 0._mm_wp ; zdm3n(:,:) = 0._mm_wp
272
273	dm0s(:) = 0._mm_wp ; dm3s(:) = 0._mm_wp
274	dm0f(:) = 0._mm_wp ; dm3f(:) = 0._mm_wp
275	dm0n(:) = 0._mm_wp ; dm3n(:) = 0._mm_wp
276	dm3i(:,:) = 0._mm_wp ; dmugas(:,:) = 0._mm_wp
277
278	! Calls nucleation/condensation for each species:
279	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
280	DO i = 1, mm_nesp
281	CALL nc_esp(mm_xESPS(i),mugas(:,i),m3ice(:,i), &
282	m0as(:),m3as(:),m0af(:),m3af(:),m0ccn(:),m3ccn(:),&
283	zdm0s(:,i),zdm3s(:,i),zdm0f(:,i),zdm3f(:,i), &
284	zdm0n(:,i),zdm3n(:,i),dm3i(:,i),dmugas(:,i), &
285	gassat(:,i),nrate(:,i),grate(:,i))
286	ENDDO
287
288	! Retrieve the index of the maximum tendency of CCN where ice variation is not null:
289	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
290	DO i = 1, mm_nla
291	idx = MAXLOC(zdm0n(i,:),DIM=1,MASK=(dm3i(i,:) /= 0._mm_wp .OR. m3ice(i,:)+dm3i(i,:) >= 0._mm_wp))
292
293	IF (idx == 0) THEN
294	dm0s(i) = 0._mm_wp
295	dm3s(i) = 0._mm_wp
296	dm0f(i) = 0._mm_wp
297	dm3f(i) = 0._mm_wp
298	dm0n(i) = 0._mm_wp
299	dm3n(i) = 0._mm_wp
300
301	ELSE
302	IF (mm_debug) THEN
303	WRITE(,'((a),I2.2,(a),ES10.3,(a))') "[MM_DEBUG - mm_cloud_nucond] Z(",i,") = ", mm_play(i)1e2, &
304	" mbar: Max aer/ccn exchange variation due to species: "//TRIM(mm_xESPS(idx)%name)
305	ENDIF
306	dm0s(i) = zdm0s(i,idx)
307	dm3s(i) = zdm3s(i,idx)
308	dm0f(i) = zdm0f(i,idx)
309	dm3f(i) = zdm3f(i,idx)
310	dm0n(i) = zdm0n(i,idx)
311	dm3n(i) = zdm3n(i,idx)
312	ENDIF
313	ENDDO
314	END SUBROUTINE mm_cloud_nucond
315
316
317	SUBROUTINE nc_esp(xESP,Xvap,Xm3ice,m0as,m3as,m0af,m3af,m0ccn,m3ccn,&
318	dm0as,dm3as,dm0af,dm3af,dm0ccn,dm3ccn,Xdm3ice,Xdvap,Xsat,Xnrate,Xgrate)
319	!! Get moments tendencies through nucleation/condensation/evaporation of a given condensible species.
320	!! The method computes the global tendencies of the aerosols, ccn and ice moments through cloud
321	!! microphysics processes (nucleation & condensation/evaporation).
322	!!
323	!! @warning
324	!! Input quantities __m0aer__,__m3aer__,__m0ccn__,__m3ccn__,__m3ice__ are assumed to be in (X.m-3),
325	!! where X is the unit of the moment that is, a number for M0 and a volume for M3.
326	!! __Xvap__ must be expressed in term of molar fraction.
327	!!
328	!! Implicit scheme:
329	!! For nucleation we have the following equations:
330	!! (1) dM_aer(k)/dt = - dM_ccn(k)/dt
331	!! (2) dM_aer(k)/dt = - (4 * pi * Jhet) / rm * M_aer(k+3)
332	!! = - (4 * pi * Jhet) / rm * alpha(k+3)/alpha(k) * rc*3 M_aer(k)
333	!! With:
334	!! CST_M(k) = (4 * pi * Jhet) / rm * alpha(k+3)/alpha(k) * rc*3 dt
335	!! We solve:
336	!! (3) M_aer(k)[t+dt] = (1 / (1+CST_M(k))) * M_aer(k)[t]
337	!! Then, from (2):
338	!! (4) M_ccn(k)[t+dt] = M_ccn(k)[t] + (CST_M(k)/(1+CST_M(k))) * M_aer(k)[t]
339	!!
340	! Condensate species properties.
341	TYPE(mm_esp), INTENT(in) :: xESP
342	! Gas species molar fraction on the vertical grid from __TOP__ to __GROUND__.
343	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: Xvap ! (mol.mol-1)
344	! 3rd order moment of the ice component.
345	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: Xm3ice ! (m3.m-3)
346	! 0th and 3rd order moment of the aerosols (spherical mode) component.
347	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m0as ! (m-3)
348	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m3as ! (m3.m-3)
349	! 0th and 3rd order moment of the aerosols (fractal mode) component.
350	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m0af ! (m-3)
351	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m3af ! (m3.m-3)
352	! 0th and 3rd order moment of the cloud condensation nuclei component.
353	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m0ccn ! (m-3)
354	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m3ccn ! (m3.m-3)
355	! Tendency of the 0th and 3rd order moment of the aerosols (spherical mode) through cloud microphysics processes.
356	REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: dm0as ! (m-3)
357	REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: dm3as ! (m3.m-3)
358	! Tendency of the 0th and 3rd order moment of the aerosols (fractal mode) through cloud microphysics processes.
359	REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: dm0af ! (m-3)
360	REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: dm3af ! (m3.m-3)
361	! Tendency of the 0th and 3rd order moment of the ccn through cloud microphysics processes.
362	REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: dm0ccn ! (m-3)
363	REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: dm3ccn ! (m3.m-3)
364	! Tendency of the 3rd order moment of the ice component.
365	REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: Xdm3ice ! (m3.m-3)
366	! Tendency of gas species.
367	REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: Xdvap ! (mol.mol-1)
368	! Saturation ratio values on the vertical grid.
369	REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: Xsat
370	! Nucleation rate values on the vertical grid for the species X.
371	REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: Xnrate ! (m-2.s-1)
372	! Growth rate values on the vertical grid for the species X.
373	REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: Xgrate ! (m2.s-1)
374
375	! Local variables:
376	!~~~~~~~~~~~~~~~~~
377	INTEGER :: i, nucl_k
378	REAL(kind=mm_wp) :: bef,aft
379	REAL(kind=mm_wp), DIMENSION(SIZE(Xvap)) :: s_m0as,s_m3as,s_m0af,s_m3af,s_m0ccn,s_m3ccn,s_Xm3ice,s_Xvap
380	REAL(kind=mm_wp), DIMENSION(SIZE(Xvap)) :: z_m0as,z_m3as,z_m0af,z_m3af,z_m0ccn,z_m3ccn,z_Xm3ice,z_Xvap
381	REAL(kind=mm_wp), DIMENSION(SIZE(Xvap)) :: drad,sig,qsat,pX
382	REAL(kind=mm_wp), DIMENSION(SIZE(Xvap)) :: Jhet_aers,cst_m0aers,cst_m3aers
383	REAL(kind=mm_wp), DIMENSION(SIZE(Xvap)) :: Jhet_aerf,cst_m0aerf,cst_m3aerf
384	REAL(kind=mm_wp), DIMENSION(SIZE(Xvap)) :: S_eq,g_rate
385	REAL(kind=mm_wp), DIMENSION(SIZE(Xvap)) :: ctot,up,down,newvap
386
387	REAL(kind=mm_wp), PARAMETER :: athird = 1._mm_wp / 3._mm_wp
388	REAL(kind=mm_wp), PARAMETER :: pifac = (4._mm_wp * mm_pi) / 3._mm_wp
389
390	! Initialization:
391	!~~~~~~~~~~~~~~~~
392	! Types of Nucleation
393	! Default value (0 = Jaerf, 1 = Jaers, else = Jaers + Jaerf)
394	nucl_k = 2
395
396	! Copy input argument and convert units X.m-3 to X.kg-1 or X.mol-1 to X.kg-1.
397	! @warning: we must update the aerosol tracers through haze microphysics processes.
398
399	! s_XXX is the initial converted value saved:
400	s_m0as = m0as / mm_rhoair ; s_m3as = m3as / mm_rhoair
401	s_m0af = m0af / mm_rhoair ; s_m3af = m3af / mm_rhoair
402	s_m0ccn = m0ccn / mm_rhoair
403	s_m3ccn = m3ccn / mm_rhoair
404	s_Xm3ice = Xm3ice / mm_rhoair
405	s_Xvap = Xvap * xESP%fmol2fmas
406
407	! z_XXX is our working copy:
408	z_m0as = s_m0as ; z_m3as = s_m3as
409	z_m0af = s_m0af ; z_m3af = s_m3af
410	z_m0ccn = s_m0ccn ; z_m3ccn = s_m3ccn
411	z_Xm3ice = s_Xm3ice
412	z_Xvap = s_Xvap
413
414	! Initialize local variables:
415	bef = 0._mm_wp ; aft = 0._mm_wp
416	drad(:) = 0._mm_wp
417	sig(:) = 0._mm_wp ; qsat(:) = 0._mm_wp ; pX(:) = 0._mm_wp
418	Jhet_aers(:) = 0._mm_wp ; cst_m0aers(:) = 0._mm_wp ; cst_m3aers(:) = 0._mm_wp
419	Jhet_aerf(:) = 0._mm_wp ; cst_m0aerf(:) = 0._mm_wp ; cst_m3aerf(:) = 0._mm_wp
420	S_eq(:) = 0._mm_wp ; g_rate(:) = 0._mm_wp
421	ctot(:) = 0._mm_wp ; up(:) = 0._mm_wp ; down(:) = 0._mm_wp
422	newvap(:) = 0._mm_wp
423
424	! Get a copy of drop radius:
425	drad(:) = mm_drad(:)
426
427	! Surface tension (N.m-1):
428	!~~~~~~~~~~~~~~~~~~~~~~~~~
429	sig = mm_sigX(mm_temp,xESP)
430
431	! Species mass mixing ratio at saturation (kg.kg-1):
432	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
433	qsat = mm_qsatX(mm_temp,mm_play,xESP) * xESP%fmol2fmas
434
435	! Partial pressure of species:
436	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
437	pX = Xvap * mm_play
438
439	! Saturation ratio:
440	!~~~~~~~~~~~~~~~~~~
441	where (qsat > 0._mm_wp)
442	Xsat = z_Xvap / qsat
443	elsewhere
444	Xsat = 0._mm_wp
445	endwhere
446
447	!~~~~~~~~~~~~~~~~~~~~~
448	! GETS NUCLEATION RATE
449	!~~~~~~~~~~~~~~~~~~~~~
450
451	! Gets heterogeneous nucleation rate:
452	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
453	if (nucl_k == 1) then
454	! Spherical aerosols:
455	!~~~~~~~~~~~~~~~~~~~~
456	call hetnucl_rate_aer(mm_rcs,mm_temp,xESP,pX,Xsat,Jhet_aers)
457
458	! We set the drop radius to the spherical radius.
459	! Doing so will prevent a nasty bug to occur later when ice volume is updated!
460	where (drad <= mm_drad_min .AND. Jhet_aers > 0._mm_wp)
461	drad = mm_rcs
462	endwhere
463
464	else if (nucl_k == 0) then
465	! Fractal aerosols (ccn radius is the monomer):
466	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
467	call hetnucl_rate_aer((/(mm_rm, i=1,mm_nla)/),mm_temp,xESP,pX,Xsat,Jhet_aerf)
468
469	! ==> We set the drop radius to the fractal radius.
470	! Doing so will prevent a nasty bug to occur later when ice volume is updated!
471	where (drad <= mm_drad_min .AND. Jhet_aerf > 0._mm_wp)
472	drad = mm_rm
473	endwhere
474
475	else
476	! Spherical aerosols + Fractal aerosols:
477	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
478	call hetnucl_rate_aer(mm_rcs,mm_temp,xESP,pX,Xsat,Jhet_aers)
479	call hetnucl_rate_aer((/(mm_rm, i=1,mm_nla)/),mm_temp,xESP,pX,Xsat,Jhet_aerf)
480
481	! ==> We set the drop radius to the minimum radius.
482	! Doing so will prevent a nasty bug to occur later when ice volume is updated!
483	where (drad <= mm_drad_min .AND. (Jhet_aers > 0._mm_wp .OR. Jhet_aerf > 0._mm_wp))
484	drad = MIN(mm_rcs,mm_rm)
485	endwhere
486
487	endif
488
489	! Diagnostics
490	Xnrate = Jhet_aers + Jhet_aerf
491
492	! /!\ IMPORTANT:
493	! Update CCN and aerosols moment from nucleation now !
494	! Doing so should prevent a nasty bug that occurs if we want to generate clouds from scratch
495	! (i.e. a "dry" atmosphere without any clouds tracers already present).
496	!
497	! In such case, we do not produce ice variation on the first call of the method,
498	! at most only CCN are produced (i.e. dm3i == 0, dm3n != 0)
499	! But the rules for computing the global tendencies in mm_cloud_nucond state that the global
500	! variation for CCN is due to the most active species exchange.
501
502	if (nucl_k == 1) then
503	! Spherical aerosols:
504	!~~~~~~~~~~~~~~~~~~~~
505	cst_m0aers = (4._mm_wpmm_piJhet_aers) / mm_rm * (mm_alpha_s(3._mm_wp)/mm_alpha_s(0._mm_wp)mm_rcs3) mm_dt
506	cst_m3aers = (4._mm_wpmm_piJhet_aers) / mm_rm * (mm_alpha_s(6._mm_wp)/mm_alpha_s(3._mm_wp)mm_rcs3) mm_dt
507
508	z_m0as = (1._mm_wp/(1._mm_wp+cst_m0aers)) * z_m0as
509	z_m3as = (1._mm_wp/(1._mm_wp+cst_m3aers)) * z_m3as
510
511	where (z_m0as <= 0._mm_wp .OR. z_m3as <= 0._mm_wp)
512	z_m0as = 0._mm_wp
513	z_m3as = 0._mm_wp
514	z_m0ccn = z_m0ccn + s_m0as
515	z_m3ccn = z_m3ccn + s_m3as
516	elsewhere
517	z_m0ccn = z_m0ccn + cst_m0aers*z_m0as
518	z_m3ccn = z_m3ccn + cst_m3aers*z_m3as
519	endwhere
520
521	else if (nucl_k == 0) then
522	! Fractal aerosols
523	!~~~~~~~~~~~~~~~~~
524	cst_m0aerf = (4._mm_wpmm_piJhet_aerf) / mm_rm * (mm_alpha_f(3._mm_wp)/mm_alpha_f(0._mm_wp)mm_rcf3) mm_dt
525	cst_m3aerf = (4._mm_wpmm_piJhet_aerf) / mm_rm * (mm_alpha_f(6._mm_wp)/mm_alpha_f(3._mm_wp)mm_rcf3) mm_dt
526
527	z_m0af = (1._mm_wp/(1._mm_wp+cst_m0aerf)) * z_m0af
528	z_m3af = (1._mm_wp/(1._mm_wp+cst_m3aerf)) * z_m3af
529
530	where (z_m0af <= 0._mm_wp .OR. z_m3af <= 0._mm_wp)
531	z_m0af = 0._mm_wp
532	z_m3af = 0._mm_wp
533	z_m0ccn = z_m0ccn + s_m0af
534	z_m3ccn = z_m3ccn + s_m3af
535	elsewhere
536	z_m0ccn = z_m0ccn + cst_m0aerf*z_m0af
537	z_m3ccn = z_m3ccn + cst_m3aerf*z_m3af
538	endwhere
539
540	else
541	! Spherical aerosols + Fractal aerosols:
542	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
543	cst_m0aers = (4._mm_wpmm_piJhet_aers) / mm_rm * (mm_alpha_s(3._mm_wp)/mm_alpha_s(0._mm_wp)mm_rcs3) mm_dt
544	cst_m3aers = (4._mm_wpmm_piJhet_aers) / mm_rm * (mm_alpha_s(6._mm_wp)/mm_alpha_s(3._mm_wp)mm_rcs3) mm_dt
545
546	z_m0as = (1._mm_wp/(1._mm_wp+cst_m0aers)) * z_m0as
547	z_m3as = (1._mm_wp/(1._mm_wp+cst_m3aers)) * z_m3as
548
549	where (z_m0as <= 0._mm_wp .OR. z_m3as <= 0._mm_wp)
550	z_m0as = 0._mm_wp
551	z_m3as = 0._mm_wp
552	z_m0ccn = z_m0ccn + s_m0as
553	z_m3ccn = z_m3ccn + s_m3as
554	elsewhere
555	z_m0ccn = z_m0ccn + cst_m0aers*z_m0as
556	z_m3ccn = z_m3ccn + cst_m3aers*z_m3as
557	endwhere
558
559	cst_m0aerf = (4._mm_wpmm_piJhet_aerf) / mm_rm * (mm_alpha_f(3._mm_wp)/mm_alpha_f(0._mm_wp)mm_rcf3) mm_dt
560	cst_m3aerf = (4._mm_wpmm_piJhet_aerf) / mm_rm * (mm_alpha_f(6._mm_wp)/mm_alpha_f(3._mm_wp)mm_rcf3) mm_dt
561
562	z_m0af = (1._mm_wp/(1._mm_wp+cst_m0aerf)) * z_m0af
563	z_m3af = (1._mm_wp/(1._mm_wp+cst_m3aerf)) * z_m3af
564
565	where (z_m0af <= 0._mm_wp .OR. z_m3af <= 0._mm_wp)
566	z_m0af = 0._mm_wp
567	z_m3af = 0._mm_wp
568	z_m0ccn = z_m0ccn + s_m0af
569	z_m3ccn = z_m3ccn + s_m3af
570	elsewhere
571	z_m0ccn = z_m0ccn + cst_m0aerf*z_m0af
572	z_m3ccn = z_m3ccn + cst_m3aerf*z_m3af
573	endwhere
574
575	endif
576
577	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
578	! GETS CONDENSATION/EVAPORATION RATE
579	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
580
581	! Equilibrium saturation near the drop:
582	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
583	S_eq = exp(min((2._mm_wpsigxESP%masmol) / (mm_rgasmm_tempxESP%rho_s*drad),30._mm_wp))
584
585	! Gets growth rate:
586	!~~~~~~~~~~~~~~~~~~
587	call growth_rate(mm_temp,mm_play,pX/Xsat,xESP,S_eq,drad,g_rate)
588
589	ctot = z_Xvap + (z_Xm3ice * xESP%rho_s)
590	up = z_Xvap + mm_dt * g_rate * 4._mm_wp * mm_pi * xESP%rho_s * drad * S_eq * z_m0ccn
591	down = 1._mm_wp + mm_dt * g_rate * 4._mm_wp * mm_pi * xESP%rho_s * drad / qsat * z_m0ccn
592
593	! Gets new vapor X species mass mixing ratio:
594	! Cannot be greater than the total gas + ice and lower than nothing.
595	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
596	newvap = max(min(up/down,ctot),0._mm_wp)
597
598	! Gets "true" growth rate:
599	!~~~~~~~~~~~~~~~~~~~~~~~~~
600	g_rate = g_rate * (newvap/qsat - S_eq)
601
602	! Diagnostics
603	where(z_Xm3ice <= 0._mm_wp .AND. g_rate <= 0._mm_wp)
604	Xgrate = 0._mm_wp
605	elsewhere
606	Xgrate = g_rate
607	endwhere
608
609	! Update ice volume through condensation/evaporation:
610	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
611	DO i = 1, mm_nla
612	! Check for the specific case: no ice and sublimation
613	IF (z_Xm3ice(i) <= 0._mm_wp .AND. g_rate(i) <= 0._mm_wp) THEN
614	z_Xm3ice(i) = 0._mm_wp
615
616	ELSE
617	z_Xm3ice(i) = z_Xm3ice(i) + (4._mm_wp * mm_pi * drad(i) * g_rate(i) * z_m0ccn(i) * mm_dt)
618
619	! Check if there is ice left in the ccn.
620	! @note: 0.8/0.2 for now... need to be improve!
621	IF (z_Xm3ice(i) <= 0._mm_wp) THEN
622	IF (nucl_k == 1) THEN
623	z_m0as(i) = z_m0as(i) + z_m0ccn(i)
624	z_m3as(i) = z_m3as(i) + z_m3ccn(i)
625	ELSE IF (nucl_k == 0) THEN
626	z_m0af(i) = z_m0af(i) + z_m0ccn(i)
627	z_m3af(i) = z_m3af(i) + z_m3ccn(i)
628	ELSE
629	z_m0as(i) = z_m0as(i) + 0.8_mm_wp * z_m0ccn(i)
630	z_m3as(i) = z_m3as(i) + 0.8_mm_wp * z_m3ccn(i)
631	z_m0af(i) = z_m0af(i) + 0.2_mm_wp * z_m0ccn(i)
632	z_m3af(i) = z_m3af(i) + 0.2_mm_wp * z_m3ccn(i)
633	ENDIF
634	z_m0ccn(i) = 0._mm_wp
635	z_m3ccn(i) = 0._mm_wp
636	z_Xm3ice(i) = 0._mm_wp
637	ENDIF
638	ENDIF
639	ENDDO
640
641	! Sanity check:
642	!~~~~~~~~~~~~~~
643	IF (mm_debug) THEN
644	DO i = 1, mm_nla
645	bef = s_m0as(i) + s_m0af(i) + s_m0ccn(i)
646	aft = z_m0as(i) + z_m0af(i) + z_m0ccn(i)
647
648	IF (ABS(bef-aft)/bef > 1e-10_mm_wp) THEN
649	WRITE(*,'((a),I2.2,(a),ES20.12,(a),ES20.12)') &
650	"[MM_DEBUG - nc_esp] nc_esp("//TRIM(xESP%name)//"): M0 not conserved (z=",i,")",bef," <-> ",aft
651	ENDIF
652
653	bef = s_m3as(i) + s_m3af(i) + s_m3ccn(i)
654	aft = z_m3as(i) + z_m3af(i) + z_m3ccn(i)
655
656	IF (ABS(bef-aft)/bef > 1e-10_mm_wp) THEN
657	WRITE(*,'((a),I2.2,(a),ES20.12,(a),ES20.12)') &
658	"[MM_DEBUG - nc_esp] nc_esp("//TRIM(xESP%name)//"): M3 not conserved (z=",i,")",bef," <-> ",aft
659	ENDIF
660	ENDDO
661	ENDIF
662
663	! Compute tendencies (in X.m-3 or X.mol-1):
664	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
665	dm0as = (z_m0as - s_m0as) * mm_rhoair ; dm3as = (z_m3as - s_m3as) * mm_rhoair
666	dm0af = (z_m0af - s_m0af) * mm_rhoair ; dm3af = (z_m3af - s_m3af) * mm_rhoair
667	dm0ccn = (z_m0ccn - s_m0ccn) * mm_rhoair ; dm3ccn = (z_m3ccn - s_m3ccn) * mm_rhoair
668
669	Xdm3ice = (z_Xm3ice - s_Xm3ice) * mm_rhoair
670	Xdvap = - (z_Xm3ice - s_Xm3ice) * xESP%rho_s / xESP%fmol2fmas
671
672	END SUBROUTINE nc_esp
673
674	SUBROUTINE hetnucl_rate_aer(rccn,temp,xESP,pvp,sat,rate)
675	!! Get heterogeneous nucleation rate on aerosols.
676	!! The method computes the heterogeneous nucleation rate for the given species on a aerosols of size __rccn__.
677	!! Except __xESP__, all arguments are vectors of the same size (vertical grid).
678	!!
679	REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: rccn ! Radius of the cloud condensation nuclei (m).
680	REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: temp ! Temperature (K).
681	TYPE(mm_esp), INTENT(in) :: xESP ! Species properties.
682	REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: pvp ! Partial vapor pressure of X species (Pa).
683	REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: sat ! Saturation ratio of given species.
684	REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: rate ! The nucleation rate (m-2.s-1).
685
686	! Local variables:
687	!~~~~~~~~~~~~~~~~~
688	INTEGER :: i
689	REAL(kind=mm_wp) :: S, T, r
690	REAL(kind=mm_wp) :: sig,nX,rstar
691	REAL(kind=mm_wp) :: x,fsh,deltaFstar
692	REAL(kind=mm_wp) :: deltaF, gstar,zeldov
693
694	! Initialization:
695	!~~~~~~~~~~~~~~~~
696	rate(:) = 0._mm_wp
697
698	! Activation condition:
699	!~~~~~~~~~~~~~~~~~~~~~~
700	DO i = 1, SIZE(rccn)
701	S = sat(i)
702
703	IF (S > 1._mm_wp) THEN
704	T = temp(i)
705	r = rccn(i)
706
707	sig = mm_sigX(T,xESP)
708	nX = pvp(i) / (mm_kboltz*T)
709	rstar = (2._mm_wpsigxESP%vol) / (mm_kboltzTdlog(S))
710
711	! Curvature radius and shape factor
712	x = r / rstar
713	fsh = mm_fshape(xESP%mteta,x)
714	deltaFstar = (4._mm_wpmm_pi/3._mm_wp) sig * (rstar*2.) fsh
715
716	deltaF = MIN(MAX((2.xESP%fdes - xESP%fdif - deltaFstar) / (mm_kboltzT),-100._mm_wp),100._mm_wp)
717
718	IF (deltaF > -100._mm_wp) THEN
719	gstar = ((4._mm_wpmm_pi/3._mm_wp) (rstar**3)) / (xESP%vol)
720	zeldov = dsqrt(deltaFstar / (3._mm_wpmm_pimm_kboltzT(gstar**2)))
721
722	! Heterogeneous nucleation rate
723	rate(i) = ((zeldovmm_kboltzT) / (fshxESP%nusxESP%mas)) * (nXrstar)2._mm_wp dexp(deltaF)
724	ENDIF ! (deltaF > -100._mm_wp)
725	ENDIF ! (S > 1._mm_wp)
726	ENDDO ! SIZE(rccn)
727
728	RETURN
729
730	END SUBROUTINE hetnucl_rate_aer
731
732	SUBROUTINE growth_rate(temp,pres,pXsat,xESP,S_eq,drad,rate)
733	!! Get growth rate through condensation/evaporation process.
734	!! The method computes the growth rate a drop through condensation/evaporation processes:
735	!! Except __xESP__ which is a scalar, all arguments are vectors of the same size (vertical grid).
736	!!
737	REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: temp ! Temperature (K).
738	REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: pres ! Pressure level (Pa).
739	REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: pXsat ! Saturation vapor pressure of species (Pa).
740	TYPE(mm_esp), INTENT(in) :: xESP ! Specie properties.
741	REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: S_eq ! Equilibrium saturation near the drop.
742	REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: drad ! Drop radius (m).
743	REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: rate ! Growth rate (m2.s-1).
744
745	! Local variables:
746	!~~~~~~~~~~~~~~~~~
747	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: k, L, knu, slf
748	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: Dv, Rc, Rd
749
750	ALLOCATE(k(SIZE(temp)),L(SIZE(temp)),knu(SIZE(temp)),slf(SIZE(temp)))
751	ALLOCATE(Dv(SIZE(temp)),Rc(SIZE(temp)),Rd(SIZE(temp)))
752
753	! Initialization:
754	!~~~~~~~~~~~~~~~~
755	k(:) = 0._mm_wp ; L(:) = 0._mm_wp ; knu(:) = 0._mm_wp ; slf(:) = 0._mm_wp
756	Dv(:) = 0._mm_wp ; Rc(:) = 0._mm_wp ; Rd(:) = 0._mm_wp
757
758	! N2 (air) Thermal conductivity:
759	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
760	k(:) = (2.857e-2_mm_wptemp - 0.5428_mm_wp) 4.184e-3_mm_wp
761
762	! Gas mean free path:
763	!~~~~~~~~~~~~~~~~~~~~
764	L(:) = mm_lambda_air(temp,pres)
765
766	! The knudsen number of the drop:
767	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
768	knu(:) = L(:) / drad(:)
769
770	! Slip flow correction:
771	!~~~~~~~~~~~~~~~~~~~~~~
772	slf(:) = knu(:) * (1.333_mm_wp + 0.71_mm_wp/knu(:)) / (1._mm_wp + 1._mm_wp/knu(:))
773
774	! Molecular diffusivity coefficient of each species:
775	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
776	Dv(:) = (1._mm_wp/3._mm_wp) * dsqrt(8._mm_wpmm_rgastemp(:) / (mm_pixESP%masmol)) &
777	mm_kboltztemp(:) / (mm_pi pres(:) * (mm_air_rad+xESP%ray)*2 dsqrt(1._mm_wp+xESP%fmol2fmas))
778
779	! Transitional regime:
780	!~~~~~~~~~~~~~~~~~~~~~
781	Dv(:) = Dv(:) / (1._mm_wp + slf(:))
782
783	! Latent heat resistance coefficient Rc:
784	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
785	Rc(:) = (mm_LheatX(temp(:),xESP)*2 xESP%rho_s * xESP%masmol) / (k(:) * mm_rgas * temp(:)**2)
786
787	! Molecular diffusion resistance coefficient Rd:
788	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
789	Rd(:) = (mm_rgas * temp(:) * xESP%rho_s) / (Dv(:) * pXsat(:) * xESP%masmol)
790
791	! Growth rate: rdr/dt = rate * (S - S_eq):
792	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
793	rate(:) = 1._mm_wp / (S_eq(:)*Rc(:) + Rd(:))
794
795	RETURN
796
797	END SUBROUTINE growth_rate
798
799	!============================================================================
800	! COAGULATION PROCESS RELATED METHODS
801	!============================================================================
802
803	SUBROUTINE mm_cloud_coagulation(m0ccn,m3ccn,m3ice,dm0n,dm3n)
804	!! Get the evolution of the cloud droplets moments vertical column due to coagulation process.
805	!!
806	! 0th and 3rd order moment of the cloud condensation nuclei component.
807	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m0ccn ! (m-3)
808	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m3ccn ! (m3.m-3)
809	! 3rd order moment of ices and gases component.
810	REAL(kind=mm_wp), DIMENSION(:,:), INTENT(in) :: m3ice ! (m-3)
811	! Tendency of the 0th order moment of the droplets size-distribution over a time step (m-3).
812	REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm0n
813	! Tendency of the 3rd order moment of the droplets size-distribution (m3.m-3).
814	REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3n
815
816	! Local variables:
817	!~~~~~~~~~~~~~~~~~
818	INTEGER :: i
819	REAL(kind=mm_wp), DIMENSION(SIZE(dm0n)) :: Gij
820
821	! 0. Initialization:
822	!~~~~~~~~~~~~~~~~~~~
823	dm0n(:) = 0._mm_wp ; dm3n(:) = 0._mm_wp
824
825	Gij(:) = 0._mm_wp
826
827	! 1. Molecular's case:
828	!~~~~~~~~~~~~~~~~~~~~~
829	! Types of Coagulation
830	where (mm_drad > 0._mm_wp .and. mm_drho > 0._mm_wp)
831	! Pre-factor:
832	Gij = sqrt((12._mm_wp * mm_kboltz * mm_btemp) / (mm_pi*2 mm_drad*3 mm_drho))
833
834	! Tendencies of the moments:
835	dm0n = - mm_pi * Gij * m0ccn*2 mm_drad*2 (mm_alpha_ccn(1._mm_wp)*2 + mm_alpha_ccn(2._mm_wp)) mm_dt
836	dm3n = 0._mm_wp
837	endwhere
838
839	! 2. Transitory's case:
840	!~~~~~~~~~~~~~~~~~~~~~~
841	!where (mm_drad > 0._mm_wp .and. mm_drho > 0._mm_wp)
842	! ! Pre-factor:
843	! Gij(:) = sqrt((6._mm_wp * mm_kboltz * mm_btemp(:)) / mm_drho(:))
844	!
845	! ! Tendencies of the moments:
846	! dm0n = - mm_get_btk(3,0) * Gij * m0ccn*2 mm_drad*(1._mm_wp/2._mm_wp) &
847	! (mm_alpha_ccn(1._mm_wp/2._mm_wp) + &
848	! mm_alpha_ccn(-3._mm_wp/2._mm_wp)*mm_alpha_ccn(2._mm_wp) + &
849	! 2._mm_wpmm_alpha_ccn(-1._mm_wp/2._mm_wp)mm_alpha_ccn(1._mm_wp))
850	! dm3n = 0._mm_wp
851	!endwhere
852
853	RETURN
854	END SUBROUTINE mm_cloud_coagulation
855
856	!============================================================================
857	! SEDIMENTATION PROCESS RELATED METHODS
858	!============================================================================
859
860	SUBROUTINE mm_cloud_sedimentation(m0ccn,m3ccn,m3ice,dm0n,dm3n,dm3i)
861	!! Compute the tendency of clouds related moments through sedimentation process.
862	!!
863	!! The method computes the tendencies of moments related to cloud microphysics through
864	!! sedimentation process. The algorithm used here differs from
865	!! [[mm_haze_sedimentation(subroutine)]] as all moments settle with the same
866	!! terminal velocity which is computed with the average drop radius of the size distribution.
867	!! We simply compute an exchange matrix that stores the new positions of each cells through
868	!! sedimentation process and then computes the matrix
869	!! product with input moments values to get final tendencies.
870	!!
871	! 0th and 3rd order moment of the cloud condensation nuclei component.
872	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m0ccn ! (m-3)
873	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m3ccn ! (m3.m-3)
874	! 3rd order moment of ices and gases component.
875	REAL(kind=mm_wp), DIMENSION(:,:), INTENT(in) :: m3ice ! (m-3)
876	! Tendency of the 0th order moment of the ccn distribution (m-3).
877	REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm0n
878	! Tendency of the 3rd order moment of the ccn distribution (m3.m-3).
879	REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3n
880	! Tendencies of the 3rd order moment of each ice component of the cloud (m3.m-3).
881	REAL(kind=mm_wp), INTENT(out), DIMENSION(:,:) :: dm3i
882
883	! Local variables:
884	!~~~~~~~~~~~~~~~~~
885	INTEGER :: i, nm
886	REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE :: moms_i, moms_f, chg_matrix
887
888	! Initialization:
889	!~~~~~~~~~~~~~~~~
890	dm0n(:) = 0._mm_wp ; dm3n(:) = 0._mm_wp ; dm3i(:,:) = 0._mm_wp
891
892	nm = 2 + mm_nesp
893	ALLOCATE(moms_i(mm_nla,nm),moms_f(mm_nla,nm),chg_matrix(mm_nla,mm_nla))
894
895	moms_i(:,1) = m0ccn * mm_dzlev
896	moms_i(:,2) = m3ccn * mm_dzlev
897	DO i = 1, mm_nesp
898	moms_i(:,2+i) = m3ice(:,i) * mm_dzlev
899	ENDDO
900
901	! Computes exchange matrix:
902	!~~~~~~~~~~~~~~~~~~~~~~~~~~
903	CALL exchange(mm_drad,mm_drho,mm_dt,chg_matrix)
904
905	! Computes final moments values:
906	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
907	DO i = 1, nm
908	moms_f(:,i) = MATMUL(chg_matrix,moms_i(:,i))
909	ENDDO
910
911	! Computes tendencies (converted in X.m-3):
912	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
913	dm0n = (moms_f(:,1) - moms_i(:,1)) / mm_dzlev
914	dm3n = (moms_f(:,2) - moms_i(:,2)) / mm_dzlev
915	DO i = 1, mm_nesp
916	dm3i(:,i) = (moms_f(:,2+i) - moms_i(:,2+i)) / mm_dzlev
917	ENDDO
918
919	RETURN
920	END SUBROUTINE mm_cloud_sedimentation
921
922	SUBROUTINE exchange(rad,rhog,dt,matrix)
923	!! Compute the exchange matrix.
924	!!
925	!! The subroutine computes the matrix exchange used by [[mm_cloud_sedimentation(subroutine)]]
926	!! to compute moments tendencies through sedimentation process. Both rad and rhog must be vector with relevant
927	!! values over the atmospheric vertical structure.
928	!! matrix is square 2D-array with same dimension size than rad.
929	!!
930	! Cloud drop radius over the atmospheric vertical structure (m).
931	REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: rad
932	! Cloud drop density over the atmospheric vertical structure (kg.m-3).
933	REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: rhog
934	! Timestep (s).
935	REAL(kind=mm_wp), INTENT(in) :: dt
936	! The output _exchange matrix.
937	REAL(kind=mm_wp), INTENT(out) :: matrix(:,:)
938
939	! Local variables:
940	!~~~~~~~~~~~~~~~~~
941	INTEGER :: nz,i,j,jj,jinf,jsup
942	REAL(kind=mm_wp) :: zni,znip1,xf,xft,xcnt
943	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: puit
944	REAL(kind=mm_wp) :: cpte,cpte2
945	REAL(kind=mm_wp) :: zsurf
946
947	INTEGER, PARAMETER:: ichx = 1
948
949	! Initialization:
950	!~~~~~~~~~~~~~~~~
951	matrix = 0._mm_wp
952	nz = SIZE(rad)
953	zsurf = mm_zlev(nz+1)
954	ALLOCATE(puit(nz))
955
956	! Compute exchange matrix:
957	!~~~~~~~~~~~~~~~~~~~~~~~~~
958	DO i=1,nz
959	puit(i) = 0._mm_wp
960	xcnt = 0._mm_wp
961	! Computes drop move (i.e. its new positions)
962	CALL getnzs(ichx,i,rad,rhog,dt,zni,znip1)
963
964	! Peculiar case : Ground level precipitation [znip1 < zsurf && (zni < zsurf \|\| zni > zsurf)]
965	! - complete precipitation [ znip1 <= 0 && zni <= 0 ] :
966	IF(zni <= zsurf .and. znip1 <= zsurf) THEN
967	xft=0._mm_wp
968	xf=1._mm_wp
969	xcnt=xcnt+xf
970	puit(i)=puit(i)+xf
971	ENDIF
972	! - partial precipitation [ znip1 <= zsurf && zni > zsurf ] :
973	IF (zni > zsurf .and. znip1 <= zsurf) THEN
974	xft=(zni-zsurf)/(zni-znip1)
975	xf=(1.-xft)
976	xcnt=xcnt+xf
977	puit(i)=puit(i)+xf
978	ENDIF
979
980	! General case : no ground precipitation [ znip1 > zsurf && zni > zsurf ]
981	IF (zni > zsurf .and. znip1 > zsurf) THEN
982	xft = 1._mm_wp ! on a la totalite de la case
983	xf = 0._mm_wp
984	xcnt=xcnt+xf
985	puit(i)=puit(i)+xf
986	ENDIF
987	IF (zni < znip1) THEN
988	WRITE(*,'("[EXCHANGES] WARNING, missing this case :",2(2X,ES10.3))') zni, znip1
989	ENDIF
990
991	! Fix minimum level to the ground
992	znip1 = MAX(znip1,zsurf)
993	zni = MAX(zni,zsurf)
994
995	! Locate new "drop" position in the verical grid
996	jsup=nz+1
997	jinf=nz+1
998	DO j=1,nz
999	IF (zni<=mm_zlev(j).and.zni>=mm_zlev(j+1)) jsup=j
1000	IF (znip1<=mm_zlev(j).and.znip1>=mm_zlev(j+1)) jinf=j
1001	ENDDO
1002
1003	! Volume is out of range: (all drops have touched the ground!)
1004	! Note: cannot happen here, it has been treated previously :)
1005	IF (jsup>=nz+1.and.jinf==jsup) THEN
1006	WRITE(*,'(a)') "[EXCHANGE] speaking: The impossible happened !"
1007	call EXIT(666)
1008	ENDIF
1009
1010	! Volume inside a single level
1011	IF (jsup==jinf.and.jsup<=nz) THEN
1012	xf=1._mm_wp
1013	xcnt=xcnt+xft*xf
1014	matrix(jinf,i)=matrix(jinf,i)+xft*xf
1015	ENDIF
1016
1017	! Volume over 2 levels
1018	IF (jinf==jsup+1) THEN
1019	xf=(zni-mm_zlev(jinf))/(zni-znip1)
1020	xcnt=xcnt+xf*xft
1021	IF(jsup<=nz) THEN
1022	matrix(jsup,i)=matrix(jsup,i)+xft*xf
1023	ENDIF
1024	xf=(mm_zlev(jinf)-znip1)/(zni-znip1)
1025	xcnt=xcnt+xf*xft
1026	IF (jinf<=nz) THEN
1027	matrix(jinf,i)=matrix(jinf,i)+xft*xf
1028	ENDIF
1029	ENDIF
1030
1031	! Volume over 3 or more levels
1032	IF (jinf > jsup+1) THEN
1033	! first cell
1034	xf=(zni-mm_zlev(jsup+1))/(zni-znip1)
1035	xcnt=xcnt+xf*xft
1036	matrix(jsup,i)=matrix(jsup,i)+xft*xf
1037	! last cell
1038	xf=(mm_zlev(jinf)-znip1)/(zni-znip1)
1039	xcnt=xcnt+xf*xft
1040	matrix(jinf,i)=matrix(jinf,i)+xft*xf
1041	! other :)
1042	DO jj=jsup+1,jinf-1
1043	xf=(mm_zlev(jj)-mm_zlev(jj+1))/(zni-znip1)
1044	xcnt=xcnt+xf*xft
1045	matrix(jj,i)=matrix(jj,i)+xft*xf
1046	ENDDO
1047	ENDIF
1048	ENDDO
1049
1050	! Checking if everything is alright if debug enabled...
1051	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1052	IF (mm_debug) THEN
1053	cpte=0._mm_wp ; cpte2=0._mm_wp
1054	DO j=1,nz
1055	DO jj=1,nz
1056	cpte=cpte+matrix(jj,j)
1057	ENDDO
1058	ENDDO
1059	cpte2=cpte+sum(puit)
1060	IF (abs(cpte2-nz)>1.e-4_mm_wp) THEN
1061	WRITE(*,'("[MM_DEBUG - exchange] Tx expl (/nz):",2(2X,ES10.3))') cpte,cpte2
1062	ENDIF
1063	ENDIF
1064
1065	RETURN
1066	END SUBROUTINE exchange
1067
1068	SUBROUTINE getnzs(ichx,idx,rad,rho,dt,zni,zns)
1069	!! Compute displacement of a cell under sedimentation process.
1070	!!
1071	!! The method computes the new position of a drop cell through sedimentation process.
1072	!! New positions are returned in zni and zns ouptut arguments.
1073	!!
1074	! Velocity extrapolation control flag (0 for linear, 1 for exponential -preferred -).
1075	INTEGER, INTENT(in) :: ichx
1076	! Initial position of the drop (subscript of vertical layers vectors).
1077	INTEGER, INTENT(in) :: idx
1078	! Cloud drop radius over the atmospheric vertical structure (m).
1079	REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: rad
1080	! Cloud drop density over the atmospheric vertical structure (kg.m-3).
1081	REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: rho
1082	! Timestep (s).
1083	REAL(kind=mm_wp), INTENT(in) :: dt
1084	! Final layer upper boundary position (m).
1085	REAL(kind=mm_wp), INTENT(out) :: zni
1086	! Final layer lower boundary position (m).
1087	REAL(kind=mm_wp), INTENT(out) :: zns
1088
1089	! Local variables:
1090	!~~~~~~~~~~~~~~~~~
1091	INTEGER :: i,nz
1092	REAL(kind=mm_wp) :: ws,wi,w,zi,zs,rhoair_mid
1093	REAL(kind=mm_wp) :: alpha,argexp,v0,arg1,arg2
1094	REAL(kind=mm_wp), PARAMETER :: es = 30._mm_wp
1095
1096	nz = SIZE(rad)
1097
1098	! Linear extrapolation of velocity:
1099	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1100	IF (ichx==0) THEN
1101	! Velocity upper interface
1102	ws = wsettle(mm_plev(idx),mm_btemp(idx),mm_zlev(idx),mm_rhoair(idx),rho(idx),rad(idx))
1103	IF (idx==nz) THEN
1104	! Veloctity center layer
1105	rhoair_mid = (mm_rhoair(idx-1) + mm_rhoair(idx)) / 2._mm_wp
1106	wi = wsettle(mm_play(idx),mm_temp(idx),mm_zlay(idx),rhoair_mid,rho(idx),rad(idx))
1107	ELSEIF(idx<nz) THEN
1108	! Velocity lower interface
1109	wi = wsettle(mm_plev(idx+1),mm_btemp(idx+1),mm_zlev(idx+1),mm_rhoair(idx+1), &
1110	rho(idx+1),rad(idx+1))
1111	ELSE
1112	WRITE(*,'(a)') "[ERROR - getnzs] This is the fatal error..."
1113	WRITE(*,'(a)') "Index is higher than number of levels"
1114	call EXIT(111)
1115	ENDIF
1116	w = (ws+wi)/2._mm_wp
1117	zni = mm_zlev(idx)-w*dt
1118	zns = mm_zlev(idx)-mm_dzlev(idx)-w*dt
1119	RETURN
1120
1121	! Exponential extrapolation of velocity:
1122	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1123	ELSEIF(ichx==1) THEN
1124	zs = mm_zlev(idx)
1125	ws = wsettle(mm_plev(idx),mm_btemp(idx),zs,mm_rhoair(idx),rho(idx),rad(idx))
1126	zi=mm_zlay(idx)
1127	IF (idx == 1) THEN
1128	rhoair_mid = mm_rhoair(idx)
1129	ELSE
1130	rhoair_mid = (mm_rhoair(idx-1) + mm_rhoair(idx)) / 2._mm_wp
1131	ENDIF
1132	wi = wsettle(mm_play(idx),mm_temp(idx),zi,rhoair_mid,rho(idx),rad(idx))
1133	! ws & wi must be different !
1134	IF(dabs(wi-ws)/wi <= 1.e-3_mm_wp) wi=ws/1.001_mm_wp
1135	IF (wi/=0._mm_wp.AND.(ws/wi)>0._mm_wp.AND.(zs-zi)/=0._mm_wp) alpha = dlog(ws/wi)/(zs-zi) ! Alpha < 0 if wi > ws
1136	! -es < argexp < es
1137	argexp=MAX(MIN(alpha*zs,es),-es)
1138	v0 = ws/dexp(argexp)
1139	arg1=1._mm_wp+v0alphadexp(argexp)*dt
1140	argexp=MAX(MIN(alpha*(mm_zlev(idx)-mm_dzlev(idx)),es),-es)
1141	arg2=1._mm_wp+v0alphadexp(argexp)*dt
1142	IF (arg1<=0._mm_wp.OR.arg2<=0._mm_wp) THEN
1143	! correct velocity
1144	! divides the velocity argument in arg1 and arg2 :
1145	! argX=1+alphav0exp(alphaz)dt <==> argX-1=alphav0exp(alphaz)dt
1146	! ==> argX' = 1 + (argX-1)/2 <==> argX' = (1+argX)/2.
1147	DO i=1,25
1148	IF (arg1<=0._mm_wp.OR.arg2<=0._mm_wp) THEN
1149	IF (mm_debug) &
1150	WRITE(*,'((a),I2.2,(a))') "[MM_DEBUG - getnzs] Must adjust velocity (iter:",i,"/25)"
1151	arg1=(arg1+1._mm_wp)/2._mm_wp ; arg2=(arg2+1._mm_wp)/2._mm_wp
1152	ELSE
1153	EXIT
1154	ENDIF
1155	ENDDO
1156	! We have to stop
1157	IF (i>25) THEN
1158	WRITE(*,'(a)')"[ERROR - getnzs] Cannot adjust velocities !"
1159	call EXIT(111)
1160	ENDIF
1161	ENDIF
1162	zni = mm_zlev(idx)-dlog(arg1)/alpha
1163	zns = mm_zlev(idx)-mm_dzlev(idx)-dlog(arg2)/alpha
1164	RETURN
1165	ENDIF ! end of ichx
1166	END SUBROUTINE getnzs
1167
1168	ELEMENTAL FUNCTION wsettle(p,t,z,rhoair,rho,rad) RESULT(w)
1169	!! Compute the settling velocity of a spherical particle.
1170	!!
1171	!! The method computes the effective settling velocity of spherical particle of radius rad.
1172	!!
1173	REAL(kind=mm_wp), INTENT(in) :: p ! The pressure level (Pa).
1174	REAL(kind=mm_wp), INTENT(in) :: t ! The temperature (K).
1175	REAL(kind=mm_wp), INTENT(in) :: z ! The altitude level (m).
1176	REAL(kind=mm_wp), INTENT(in) :: rhoair ! Density of air (kg.m-3).
1177	REAL(kind=mm_wp), INTENT(in) :: rho ! Density of the particle (kg.m-3).
1178	REAL(kind=mm_wp), INTENT(in) :: rad ! Radius of the particle (m).
1179	REAL(kind=mm_wp) :: w ! Settling velocity (m.s-1).
1180
1181	! Local variables:
1182	REAL(kind=mm_wp), PARAMETER :: wmax = 30.0_mm_wp ! Maximal settling velocity (m.s-1)
1183	REAL(kind=mm_wp), PARAMETER :: mrcoef = 0.74_mm_wp ! Molecular reflexion coefficient
1184
1185	REAL(kind=mm_wp) :: thermal_w
1186	REAL(kind=mm_wp) :: kn, Fc, Us
1187
1188	! Molecular's case
1189	!~~~~~~~~~~~~~~~~~
1190	! Thermal velocity
1191	thermal_w = sqrt((8._mm_wp * mm_kboltz * t) / (mm_pi * mm_air_mmas))
1192
1193	! Computes settling velocity
1194	! Types of Sedimentation
1195	w = mrcoef * mm_effg(z) * (rho/rhoair) * rad / thermal_w
1196	!w = mrcoef * mm_effg(z) * (mm_rhoaer/rhoair) * mm_rm / thermal_w
1197
1198	! Hydrodynamical's case
1199	! In fact: transitory case which is the general case
1200	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1201	! Knudsen number
1202	!kn = (mm_kboltz * t) / (p * 4._mm_wp * sqrt(2._mm_wp) * mm_pi * mm_air_rad**2) / rad
1203
1204	! Computes slip-flow correction
1205	!Fc = (1._mm_wp + 1.257_mm_wpkn + 0.4_mm_wpkn*dexp(-1.1_mm_wp/kn))
1206
1207	! Computes Stokes settling velocity
1208	!Us = (2._mm_wp * rad*2 rho * mm_effg(z)) / (9._mm_wp * mm_eta_air(t))
1209
1210	! Cunningham-Millikan correction
1211	!w = Us * Fc
1212
1213	! Velocity limit: drop deformation (Lorenz 1993)
1214	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1215	!w = 1._mm_wp / ((1._mm_wp / w) + (1._mm_wp / wmax))
1216
1217	RETURN
1218	END FUNCTION wsettle
1219
1220	FUNCTION get_mass_flux(rho,M3) RESULT(flx)
1221	!! Get the mass flux of clouds related moment through sedimention.
1222	!!
1223	!! @warning
1224	!! The method is only valid for cloud moments (i.e. ice or ccn).
1225	!! It calls [[wsettle(function)]] that compute the mean settling velocity of a cloud drop.
1226	!!
1227	!! @note
1228	!! The computed flux is always positive.
1229	!!
1230	! Tracer density (kg.m-3).
1231	REAL(kind=mm_wp), INTENT(in) :: rho
1232	! Vertical profile of the total volume of tracer from __TOP__ to __GROUND__ (m3.m-3).
1233	REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: M3
1234	! Mass sedimentation fluxes at each layer from __TOP__ to __GROUND__ (kg.m-2.s-1).
1235	REAL(kind=mm_wp), DIMENSION(SIZE(M3)) :: flx
1236
1237	! Local variables
1238	INTEGER :: i
1239	REAL(kind=mm_wp), DIMENSION(SIZE(M3)) :: w_cld
1240	REAL(kind=mm_wp), SAVE :: pifac = (4._mm_wp * mm_pi) / 3._mm_wp
1241
1242	! Cloud drop settling velocity
1243	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1244	w_cld = wsettle(mm_play,mm_temp,mm_zlay,mm_rhoair,mm_drho,mm_drad)
1245
1246	! Computes flux through sedimention
1247	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1248	flx = rho * pifac * M3 * w_cld
1249
1250	RETURN
1251	END FUNCTION get_mass_flux
1252
1253	END MODULE MP2M_CLOUDS

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

Download in other formats:

Original Format