Context Navigation

mp2m_clouds.F90

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

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

Download in other formats:

Original Format