Context Navigation

mp2m_globals.F90 @ 4076

Last change on this file since 4076 was 4073, checked in by debatzbr, 3 weeks ago
Pluto PCM: Improvement of Ice Microphysics + some cleans. Allows for nucleation on spherical, fractal, or both aerosol types. Add coagulation of ice particles. BBT
File size: 66.4 KB

Line
1	MODULE MP2M_GLOBALS
2	!============================================================================
3	!
4	! Purpose
5	! -------
6	! Parameters and global variables module.
7	! Defines and initialize all the parameters and global variables that are
8	! common to all other modules of the library.
9	!
10	! It is separated in two parts :
11	! - Main parameters and global saved variables. Most of these variables should be initialized
12	! once and should hold the same value during run-time. These variables are completly public and
13	! initialized by [[mm_globals(module):mm_global_init(interface)]] method.
14	! - The second part defines a set of vectors that defines the vertical structure of the atmosphere.
15	! Each time a new atmospheric column has to be computed (either on a new timestep or on a new couple
16	! of longitude/latitude), these vectors should be intialized with new values by calling
17	! [[mm_globals(module):mm_column_init(function)]] method.
18	!
19	! @note
20	! All the vectors that represent the vertical structure of the atmosphere (altitude, pressure and
21	! temperature...) are oriented from the __TOP__ of the atmosphere to the __GROUND__.
22	!
23	! Global variables overview:
24	! - Protected variables
25	! - Control flags
26	! - Related free parameters
27	! - Physical constants
28	! - Free parameters
29	! - Miscellaneous parameters
30	! - Vertical structure part
31	!
32	! The module also contains twenty-one methods:
33	! - mm_global_init_0
34	! - mm_global_init_1
35	! - mm_column_init
36	! - mm_aerosols_init
37	! - mm_clouds_init
38	! - mm_alpha_s, mm_alpha_f, mm_alpha_ccn
39	! - mm_set_moments_thresholds
40	! - mm_get_rcs, mm_get_rcf
41	! - mm_set_moments_cld_thresholds
42	! - cldprop_sc, cldprop_ve
43	! - read_esp
44	! - mm_effg
45	! - mm_dump_parameters
46	! - check_r1, check_l1, check_i1, check_s1
47	!
48	! Authors
49	! -------
50	! B. de Batz de Trenquelléon, J. Burgalat (11/2024)
51	!
52	!============================================================================
53	USE MP2M_MPREC
54	! From lint
55	USE LINT_DATASETS
56	! From swift
57	USE SWIFT_CFGPARSE
58	USE SWIFT_STRING_OP
59	USE SWIFT_ERRORS
60	IMPLICIT NONE
61
62	PUBLIC
63
64	PRIVATE :: cldprop_sc,cldprop_ve,read_esp,check_r1,check_i1,check_l1,check_s1
65
66	! ~~~~~~~~~~~~~~~~~~~
67	! Protected variables
68	! ~~~~~~~~~~~~~~~~~~~
69	! The following variables are read-only outside this module.
70	! One must call the afferent subroutine to update them.
71
72	! Haze model:
73	!~~~~~~~~~~~~
74	! Initialization control flags (cannot be updated)
75	PROTECTED :: mm_ini,mm_ini_col,mm_ini_aer
76	! Model parameters (mm_global_init)
77	PROTECTED :: mm_dt,mm_rhoaer,mm_df,mm_rm,mm_call_CH4hazeprod,mm_p_prod,mm_rc_prod,mm_tx_prod,mm_rpla,mm_g0,mm_rb2ra
78	! Atmospheric vertical structure (mm_column_init)
79	PROTECTED :: mm_nla,mm_nle,mm_zlay,mm_zlev,mm_play,mm_plev,mm_temp,mm_rhoair,mm_btemp,mm_dzlev,mm_dzlay
80	! Moments parameters (mm_aerosols_init)
81	PROTECTED :: mm_m0aer_s, mm_m3aer_s, mm_m0aer_f, mm_m3aer_f
82	! Thresholds parameters
83	PROTECTED :: mm_m0as_min, mm_m3as_min, mm_rcs_min, mm_m0af_min, mm_m3af_min, mm_rcf_min
84
85	! Cloud model:
86	!~~~~~~~~~~~~~
87	! Initialization control flags (cannot be updated)
88	PROTECTED :: mm_ini_cld
89	! Condensible species parameters (mm_global_init)
90	PROTECTED :: mm_nesp, mm_spcname, mm_xESPS
91	! Condensible species parameters (mm_clouds_init)
92	PROTECTED :: mm_gas
93	! Moments parameters (mm_clouds_init)
94	PROTECTED :: mm_m0ccn, mm_m3ccn, mm_m3ice
95	! Thresholds parameters
96	PROTECTED :: mm_m0ccn_min, mm_m3ccn_min, mm_m3cld_min, mm_drad_min, mm_drad_max
97
98	! ~~~~~~~~~~~~~
99	! Control flags
100	! ~~~~~~~~~~~~~
101	! Haze model:
102	!~~~~~~~~~~~~
103	! Enable/Disable haze production.
104	LOGICAL, SAVE :: mm_call_hazeprod = .true.
105	! Enable/Disable haze production from CH4 photolysis.
106	LOGICAL, SAVE :: mm_call_CH4hazeprod = .true.
107	! Enable/Disable haze sedimentation.
108	LOGICAL, SAVE :: mm_call_hazesed = .true.
109	! Enable/Disable haze coagulation.
110	LOGICAL, SAVE :: mm_call_hazecoag = .true.
111	! Force all aerosols moments to fall at M0 settling velocity.
112	LOGICAL, SAVE :: mm_wsed_m0 = .false.
113	! Force all aerosols moments to fall at M3 settling velocity.
114	LOGICAL, SAVE :: mm_wsed_m3 = .false.
115	! Enable/Disable spherical probability transfert.
116	LOGICAL, SAVE :: mm_call_ps2s = .true.
117	! Enable/Disable aerosol electric charge correction.
118	LOGICAL, SAVE :: mm_call_qe = .true.
119
120	! Cloud model:
121	!~~~~~~~~~~~~~
122	! Enable/Disable clouds microphysics.
123	LOGICAL, SAVE :: mm_call_clouds = .false.
124
125	! Enable/Disable QnD debug mode (can be used for devel).
126	LOGICAL, SAVE :: mm_debug = .false.
127	! Enable/Disable log mode (for configuration only).
128	LOGICAL, SAVE :: mm_log = .false.
129	! Enable/Disable effective G for computations.
130	LOGICAL, SAVE :: mm_use_effg = .true.
131
132	! Initialization control flag [[mm_globals(module):mm_global_init(interface)]].
133	LOGICAL, PUBLIC, SAVE :: mm_ini = .false.
134	! Initialization control flag [[mm_globals(module):mm_column_init(function)]].
135	LOGICAL, PUBLIC, SAVE :: mm_ini_col = .false.
136	! Initialization control flag [[mm_globals(module):mm_aerosols_init(function)]].
137	LOGICAL, PUBLIC, SAVE :: mm_ini_aer = .false.
138	! Initialization control flag [[mm_globals(module):mm_clouds_init(function)]].
139	LOGICAL, PUBLIC, SAVE :: mm_ini_cld = .false.
140
141	! ~~~~~~~~~~~~~~~~~~~~~~~
142	! Related free parameters
143	! ~~~~~~~~~~~~~~~~~~~~~~~
144	! No mode interaction for coagulation (i.e. no coagulation at all).
145	INTEGER, PARAMETER :: mm_coag_no = 0
146	! SS mode interaction for coagulation.
147	INTEGER, PARAMETER :: mm_coag_ss = 1
148	! SF mode interaction for coagulation.
149	INTEGER, PARAMETER :: mm_coag_sf = 2
150	! FF mode interaction for coagulation.
151	INTEGER, PARAMETER :: mm_coag_ff = 4
152	! Default interactions to activate (all by default).
153	INTEGER, SAVE :: mm_coag_choice = mm_coag_ss+mm_coag_sf+mm_coag_ff
154
155	! ~~~~~~~~~~~~~~~~~~
156	! Physical constants
157	! ~~~~~~~~~~~~~~~~~~
158	! Pi number.
159	REAL(kind=mm_wp), PARAMETER :: mm_pi = 4._mm_wp*atan(1._mm_wp)
160	! Avogadro number (mol-1).
161	REAL(kind=mm_wp), PARAMETER :: mm_navo = 6.0221367e23_mm_wp
162	! Boltzmann constant (J.K-1).
163	REAL(kind=mm_wp), PARAMETER :: mm_kboltz = 1.3806488e-23_mm_wp
164	! Perfect gas constant (J.mol-1.K-1).
165	REAL(kind=mm_wp), PARAMETER :: mm_rgas = mm_kboltz * mm_navo
166	! Approximated slip-flow correction coefficient.
167	REAL(kind=mm_wp), PARAMETER :: mm_akn = 1.591_mm_wp
168
169	! ~~~~~~~~~~~~~~~
170	! Free parameters
171	! ~~~~~~~~~~~~~~~
172	! Haze model:
173	!~~~~~~~~~~~~
174	! Spherical aerosol production pressure level (Pa).
175	REAL(kind=mm_wp), SAVE :: mm_p_prod = 1.e-2_mm_wp
176	! Spherical aerosol production rate (kg.m-2.s-1).
177	REAL(kind=mm_wp), SAVE :: mm_tx_prod = 9.8e-14_mm_wp
178	! Spherical aerosol equivalent radius production (m)
179	REAL(kind=mm_wp), SAVE :: mm_rc_prod = 1.e-9_mm_wp
180	! Monomer radius (m).
181	REAL(kind=mm_wp), SAVE :: mm_rm = 1.e-8_mm_wp
182	! Fractal dimension of fractal aerosols.
183	REAL(kind=mm_wp), SAVE :: mm_df = 2._mm_wp
184	! Aerosol density (kg.m-3).
185	REAL(kind=mm_wp), SAVE :: mm_rhoaer = 8.e2_mm_wp
186
187	! Total number of aerosols minimum threshold for the spherical mode.
188	REAL(kind=mm_wp), SAVE :: mm_m0as_min = 1.e-8_mm_wp
189	! Total volume of aerosols minimum threshold for the spherical mode.
190	REAL(kind=mm_wp), SAVE :: mm_m3as_min = 1.e-35_mm_wp
191	! Characteristic radius minimum threshold for the spherical mode.
192	REAL(kind=mm_wp), SAVE :: mm_rcs_min = 1.e-9_mm_wp
193	! Total number of aerosols minimum threshold for the fractal mode.
194	REAL(kind=mm_wp), SAVE :: mm_m0af_min = 1.e-8_mm_wp
195	! Total volume of aerosols minimum threshold for the fractal mode.
196	REAL(kind=mm_wp), SAVE :: mm_m3af_min = 1.e-35_mm_wp
197	! Characteristic radius minimum threshold for the fractal mode.
198	REAL(kind=mm_wp), SAVE :: mm_rcf_min = 1.e-9_mm_wp
199
200	! Cloud model:
201	!~~~~~~~~~~~~~
202	! Total number of cloud condensation nuclei minimum threshold.
203	REAL(kind=mm_wp), SAVE :: mm_m0ccn_min = 1.e-8_mm_wp
204	! Total volume of cloud condensation nuclei minimum threshold.
205	REAL(kind=mm_wp), SAVE :: mm_m3ccn_min = 1.e-35_mm_wp
206	! Total volume of cloud drop minimum threshold.
207	REAL(kind=mm_wp), SAVE :: mm_m3cld_min = 1.e-35_mm_wp
208	! Characteristic cloud drop radius minimum threshold.
209	REAL(kind=mm_wp), SAVE :: mm_drad_min = 1.e-9_mm_wp
210	! Characteristic cloud drop radius Maximum threshold.
211	REAL(kind=mm_wp), SAVE :: mm_drad_max = 1.e-3_mm_wp
212
213	! Planet radius (m) and gravity acceleration (m.s-2).
214	! WARNING: initialization for Pluto.
215	REAL(kind=mm_wp), SAVE :: mm_rpla = 1187000._mm_wp
216	REAL(kind=mm_wp), SAVE :: mm_g0 = 0.617_mm_wp
217	! Air molecules mean radius (m), molar mass (kg.mol-1), and molecular mass (kg).
218	! WARNING: initialization for N2.
219	REAL(kind=mm_wp), SAVE :: mm_air_rad = 1.75e-10_mm_wp
220	REAL(kind=mm_wp), SAVE :: mm_air_mmol = 28.e-3_mm_wp
221	REAL(kind=mm_wp), SAVE :: mm_air_mmas = 28.e-3_mm_wp / mm_navo
222	! Microphysical time-step (s).
223	REAL(kind=mm_wp), SAVE :: mm_dt = 180._mm_wp
224
225	! ~~~~~~~~~~~~~~~~~~~~~~~~
226	! Miscellaneous parameters
227	! ~~~~~~~~~~~~~~~~~~~~~~~~
228
229	! Bulk to apparent radius (Haze model)
230	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
231	! Bulk to apparent radius conversion pre-factor (m^X).
232	!
233	! With r_a = r_b^(3/Df) . r_m^((Df-3)/(Df))
234	! Then rb2ra = r_m^((Df-3)/(Df))
235	REAL(kind=mm_wp), SAVE :: mm_rb2ra = 1._mm_wp
236
237	! Inter-moment relation (Haze model)
238	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
239	! Alpha function parameters.
240	! It stores the parameters of the inter-moments relation functions.
241	!
242	! The inter-moments relation function is represented by the sum of exponential quadratic expressions:
243	! alpha(k) = Sum_{i=1}^{n} exp(a_i.k^2 + bi.k^2 + c_i)
244	TYPE, PUBLIC :: aprm
245	! Quadratic coefficients of the quadratic expressions.
246	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: a
247	! Linear coefficients of the quadratic expressions.
248	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: b
249	! Free term of the quadratic expressions.
250	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: c
251	END TYPE
252	! Inter-moment relation set of parameters for the spherical mode.
253	TYPE(aprm), PUBLIC, SAVE :: mm_asp
254	! Inter-moment relation set of parameters for the fractal mode.
255	TYPE(aprm), PUBLIC, SAVE :: mm_afp
256
257	! Transfert probabilities (S --> F) (Haze model)
258	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
259	! Data set for linear interpolation of transfert probability (M0/CO).
260	TYPE(dset1d), PUBLIC, SAVE, TARGET :: mm_pco0p
261	! Data set for linear interpolation of transfert probability (M3/CO).
262	TYPE(dset1d), PUBLIC, SAVE, TARGET :: mm_pco3p
263	! Data set for linear interpolation of transfert probability (M0/FM).
264	TYPE(dset1d), PUBLIC, SAVE, TARGET :: mm_pfm0p
265	! Data set for linear interpolation of transfert probability (M3/FM).
266	TYPE(dset1d), PUBLIC, SAVE, TARGET :: mm_pfm3p
267
268	! Mean electric correction (Haze model)
269	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
270	! Data set for Q_SF(M0).
271	TYPE(dset2d), PUBLIC, SAVE, TARGET :: mm_qbsf0
272	! Extended values of [[mm_gcm(module):mm_qbsf0(variable)]] dataset.
273	REAL(kind=mm_wp), PUBLIC, SAVE, DIMENSION(2,2) :: mm_qbsf0_e
274	! Data set for Q_SF(M3).
275	TYPE(dset2d), PUBLIC, SAVE, TARGET :: mm_qbsf3
276	! Extended values of [[mm_gcm(module):mm_qbsf3(variable)]] dataset.
277	REAL(kind=mm_wp), PUBLIC, SAVE, DIMENSION(2,2) :: mm_qbsf3_e
278	! Data set for Q_FF(M0).
279	TYPE(dset2d), PUBLIC, SAVE, TARGET :: mm_qbff0
280	! Extended values of [[mm_gcm(module):mm_qbff0(variable)]] dataset.
281	REAL(kind=mm_wp), PUBLIC, SAVE, DIMENSION(2,2) :: mm_qbff0_e
282
283	! btk coefficients (Haze model)
284	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
285	! Coefficients for Free-molecular regime kernel approximation (b_0(t)).
286	REAL(kind=mm_wp), PUBLIC, SAVE, DIMENSION(5) :: mm_bt0 = (/1._mm_wp,1._mm_wp,1._mm_wp,1._mm_wp,1._mm_wp/)
287	! Coefficients for Free-molecular regime kernel approximation (b_3(t)).
288	REAL(kind=mm_wp), PUBLIC, SAVE, DIMENSION(5) :: mm_bt3 = (/1._mm_wp,1._mm_wp,1._mm_wp,1._mm_wp,1._mm_wp/)
289
290	! Chemical species properties (Cloud model)
291	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
292	TYPE, PUBLIC :: mm_esp
293	!! Cloud related chemical species properties.
294	!! This derived type is used in thermodynamic methods related to cloud microphysics.
295	!!
296	CHARACTER(LEN=10) :: name !! Specie name.
297	REAL(kind=mm_wp) :: mas !! Molecular weight (kg).
298	REAL(kind=mm_wp) :: vol !! Molecular volume (m3).
299	REAL(kind=mm_wp) :: ray !! Molecular radius (m).
300	REAL(kind=mm_wp) :: masmol !! Molar mass (kg.mol-1).
301	REAL(kind=mm_wp) :: rho_l !! Liquid density (kg.m-3).
302	REAL(kind=mm_wp) :: rho_s !! Ice density (kg.m-3).
303	REAL(kind=mm_wp) :: Tc !! Critical temperature (K).
304	REAL(kind=mm_wp) :: pc !! Critical pressure (Bar).
305	REAL(kind=mm_wp) :: Tb !! Boiling point temperature (K).
306	REAL(kind=mm_wp) :: w !! Acentric factor.
307	REAL(kind=mm_wp) :: a0_sat !! Saturation equation A0 coefficient.
308	REAL(kind=mm_wp) :: a1_sat !! Saturation equation A1 coefficient.
309	REAL(kind=mm_wp) :: a2_sat !! saturation equation A2 coefficient.
310	REAL(kind=mm_wp) :: a3_sat !! Saturation equation A3 coefficient.
311	REAL(kind=mm_wp) :: a4_sat !! Saturation equation A4 coefficient.
312	REAL(kind=mm_wp) :: a5_sat !! Saturation equation A5 coefficient.
313	REAL(kind=mm_wp) :: a6_sat !! Saturation equation A6 coefficient.
314	REAL(kind=mm_wp) :: mteta !! Wettability.
315	REAL(kind=mm_wp) :: fdes !! Desorption energy (J).
316	REAL(kind=mm_wp) :: fdif !! Surface diffusion energy (J).
317	REAL(kind=mm_wp) :: nus !! Jump frequency (s-1).
318	REAL(kind=mm_wp) :: fmol2fmas !! Molar fraction to mass fraction coefficient = masmol(X)/masmol(AIR)
319	END TYPE mm_esp
320
321	! Name of condensible species.
322	CHARACTER(len=30), DIMENSION(:), ALLOCATABLE, SAVE :: mm_spcname
323	! Total number of clouds condensible species.
324	INTEGER, SAVE :: mm_nesp = -1
325	! Clouds chemical species properties.
326	TYPE(mm_esp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_xESPS
327
328	! ~~~~~~~~~~~~~~~~~~~~~~~
329	! Vertical structure part
330	! ~~~~~~~~~~~~~~~~~~~~~~~
331	! Number of vertical layers.
332	INTEGER, SAVE :: mm_nla = -1
333	! Number of vertical levels.
334	INTEGER, SAVE :: mm_nle = -1
335
336	! Altitude layers (m).
337	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_zlay
338	! Altitude levels (m).
339	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_zlev
340	!> Pressure layers (Pa).
341	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_play
342	! Pressure levels (Pa).
343	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_plev
344	! Temperature vertical profile (K).
345	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_temp
346	! Air density vertical profile (kg.m-3).
347	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_rhoair
348	! Temperature vertical profil at interfaces (K).
349	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_btemp
350
351	! Atmospheric levels thickness (m).
352	! @note: __mm_dzlev__ is defined on the total number of layers and actually
353	! corresponds to the thickness of a given layer.
354	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_dzlev
355	! Atmospheric layers thickness (m).
356	! @note: __mm_dzlay__ is defined on the total number of layers. The last
357	! value of __mm_dzlay__ is set to twice the altitude of the ground layer.
358	! This value corresponds to the thickness between the center of the
359	! __GROUND__ layer and below the surface.
360	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_dzlay
361
362	! Haze model:
363	! ~~~~~~~~~~~
364	! Spherical mode - 0th order moment (m-3).
365	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m0aer_s
366	! Spherical mode - 3rd order moment (m3.m-3).
367	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m3aer_s
368	! Spherical mode - characteristic radius (m).
369	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_rcs
370	! Fractal mode - 0th order moment (m-3).
371	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m0aer_f
372	! Fractal mode - 3rd order moment (m3.m-3).
373	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m3aer_f
374	! Fractal mode - characteristic radius (m).
375	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_rcf
376
377	! Spherical aerosol precipitation (kg.m-2).
378	REAL(kind=mm_wp), SAVE :: mm_aers_prec = 0._mm_wp
379	! Fractal aerosol precipitation (kg.m-2).
380	REAL(kind=mm_wp), SAVE :: mm_aerf_prec = 0._mm_wp
381
382	! Spherical mode (M0) settling velocity (m.s-1).
383	! @note: This variable is always negative.
384	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m0as_vsed
385	! Spherical mode (M3) settling velocity (m.s-1).
386	! @note: This variable is always negative.
387	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m3as_vsed
388	! Fractal mode (M0) settling velocity (m.s-1).
389	! @note: This variable is always negative.
390	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m0af_vsed
391	! Fractal mode (M3) settling velocity (m.s-1).
392	! @note: This variable is always negative.
393	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m3af_vsed
394
395	! Spherical aerosol mass fluxes (kg.m-2.s-1).
396	! @note: This variable is always negative.
397	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_aer_s_flux
398	! Fractal aerosol mass fluxes (kg.m-2.s-1).
399	! @note: This variable is always negative.
400	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_aer_f_flux
401
402	! Cloud model:
403	! ~~~~~~~~~~~~
404	! Cloud condensation nuclei - 0th order moment (m-3).
405	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m0ccn
406	! Cloud condensation nuclei - 3rd order moment (m3.m-3).
407	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m3ccn
408	! Mean drop radius (m).
409	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_drad
410	! Mean Drop density (kg.m-3).
411	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_drho
412
413	! Ice components - 3rd order moments (m3.m-3).
414	! It is a 2D array with the vertical layers in first dimension, and
415	! the number of ice components in the second.
416	REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: mm_m3ice
417	! Condensible species molar fraction (mol.mol-1).
418	! It is a 2D array with the vertical layers in first dimension, and
419	! the number of condensible species in the second.
420	REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: mm_gas
421
422	! CCN precipitation (kg.m-2).
423	REAL(kind=mm_wp), SAVE :: mm_ccn_prec = 0._mm_wp
424	! Ice components precipitation (kg.m-2).
425	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_ice_prec
426
427	! CCN (and ices) settling velocity (m.s-1).
428	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_cld_vsed
429
430	! CCN mass fluxes (kg.m-2.s-1).
431	REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_ccn_flux
432	! Ice components sedimentation fluxes (kg.m-2.s-1).
433	REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: mm_ice_fluxes
434
435	! Condensible species saturation ratio.
436	REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: mm_gas_sat
437	! Condensible components nucleation rates (m-2.s-1).
438	REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: mm_nrate
439	! Condensible components growth rates (m2.s-1).
440	REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: mm_grate
441
442	! --- OPENMP ---------------
443	! All variables related to column computations should be private to each thread
444	!$OMP THREADPRIVATE(mm_ini_col,mm_ini_aer,mm_ini_cld)
445	!$OMP THREADPRIVATE(mm_zlay,mm_zlev,mm_play,mm_plev,mm_temp,mm_rhoair,mm_btemp,mm_dzlev,mm_dzlay)
446	!$OMP THREADPRIVATE(mm_m0aer_s,mm_m3aer_s,mm_m0aer_f,mm_m3aer_f)
447	!$OMP THREADPRIVATE(mm_m0ccn,mm_m3ccn,mm_m3ice,mm_gas)
448	!$OMP THREADPRIVATE(mm_rcs,mm_rcf,mm_drad,mm_drho)
449	!$OMP THREADPRIVATE(mm_m0as_vsed,mm_m3as_vsed,mm_m0af_vsed,mm_m3af_vsed,mm_cld_vsed)
450	!$OMP THREADPRIVATE(mm_aer_s_flux,mm_aer_f_flux,mm_ccn_prec,mm_ice_prec,mm_ccn_flux,mm_ice_fluxes,mm_gas_sat,mm_nrate,mm_grate)
451	!$OMP THREADPRIVATE(mm_m0as_min,mm_m3as_min,mm_rcs_min,mm_m0af_min,mm_m3af_min,mm_rcf_min)
452	!$OMP THREADPRIVATE(mm_m0ccn_min,mm_m3ccn_min,mm_m3cld_min,mm_drad_min,mm_drad_max)
453	!$OMP THREADPRIVATE(mm_nla,mm_nle)
454	! --------------------------
455
456	! Interface to global initialization.
457	! The method performs the global initialization of the model.
458	INTERFACE mm_global_init
459	MODULE PROCEDURE mm_global_init_0,mm_global_init_1
460	END INTERFACE mm_global_init
461
462	! Check an option from the configuration system.
463	! The method checks for an option in the configuration system and optionally
464	! set a default value if the option is not found. This is an overloaded method
465	! that can take in input either a floating point, integer, logical or string
466	! option value.
467	INTERFACE mm_check_opt
468	MODULE PROCEDURE check_r1,check_i1,check_l1,check_s1
469	END INTERFACE mm_check_opt
470
471	! Interface to cloud properties methods.
472	! The method computes clouds properties (mean drop radius and denstity) from their afferent
473	! moments. It is overloaded to compute properties at a single level or over all the vertical
474	! atmospheric structure.
475	INTERFACE mm_cloud_properties
476	MODULE PROCEDURE cldprop_sc,cldprop_ve
477	END INTERFACE mm_cloud_properties
478
479
480	CONTAINS
481
482
483	!============================================================================
484	! INITIALIZATION METHODS
485	!============================================================================
486
487	FUNCTION mm_global_init_0(dt,df,rm,rho_aer,call_CH4hazeprod,p_prod,tx_prod,rc_prod, &
488	rplanet,g0,air_rad,air_mmol, &
489	coag_interactions,call_hazeprod,call_hazesed,call_hazecoag,&
490	force_wsed_to_m0,force_wsed_to_m3, &
491	m0as_min,rcs_min,m0af_min,rcf_min, &
492	clouds,spcfile,m0ccn_min,drad_min,debug) RESULT(err)
493	!! Initialize global parameters of the model.
494	!!
495	!! The function initializes all the global parameters of the model from direct input.
496	!! Store input values in global variables
497	!!
498	!! @note
499	!! If the method fails to initialize parameters the model should be aborted as the
500	!! global variables of the model will not be correctly setup.
501	!!
502	! Microphysical timestep (s).
503	REAL(kind=mm_wp), INTENT(in) :: dt
504	! Fractal dimension of fractal aerosol (-).
505	REAL(kind=mm_wp), INTENT(in) :: df
506	! Monomer radius (m).
507	REAL(kind=mm_wp), INTENT(in) :: rm
508	! Aerosol density(kg.m-3).
509	REAL(kind=mm_wp), INTENT(in) :: rho_aer
510	! Enable/Disable production from CH4 photolysis.
511	LOGICAL, INTENT(in) :: call_CH4hazeprod
512	! Aerosol production pressure level (Pa).
513	REAL(kind=mm_wp), INTENT(in) :: p_prod
514	! Spherical aerosol production rate kg.m-2.s-1).
515	REAL(kind=mm_wp), INTENT(in) :: tx_prod
516	! Spherical aerosol characteristic radius of production (m).
517	REAL(kind=mm_wp), INTENT(in) :: rc_prod
518	! Planet radius (m).
519	REAL(kind=mm_wp), INTENT(in) :: rplanet
520	! Planet gravity acceleration at ground level (m.s-2).
521	REAL(kind=mm_wp), INTENT(in) :: g0
522	! Radius of air molecules (m).
523	REAL(kind=mm_wp), INTENT(in) :: air_rad
524	! Molar mass of air molecules (kg.mol-1).
525	REAL(kind=mm_wp), INTENT(in) :: air_mmol
526
527	! Coagulation interactions process control flag.
528	INTEGER, INTENT(in) :: coag_interactions
529	! Haze production process control flag.
530	LOGICAL, INTENT(in) :: call_hazeprod
531	! Haze sedimentation process control flag.
532	LOGICAL, INTENT(in) :: call_hazesed
533	! Haze coagulation process control flag.
534	LOGICAL, INTENT(in) :: call_hazecoag
535	! Force __all__ aerosols moments to fall at M0 settling velocity.
536	LOGICAL, INTENT(in) :: force_wsed_to_m0
537	! Force __all__ aerosols moments to fall at M3 settling velocity
538	LOGICAL, INTENT(in) :: force_wsed_to_m3
539
540	! Minimum threshold for M0 of the spherical mode (m-3).
541	REAL(kind=mm_wp), INTENT(in) :: m0as_min
542	! Minimum threshold for the characteristic radius of the spherical mode (m).
543	REAL(kind=mm_wp), INTENT(in) :: rcs_min
544	! Minimum threshold for M0 of the factal mode (m-3).
545	REAL(kind=mm_wp), INTENT(in) :: m0af_min
546	! Minimum threshold for the characteristic radius of the fractal mode (m).
547	REAL(kind=mm_wp), INTENT(in) :: rcf_min
548
549	! Clouds microphysics control flag.
550	LOGICAL, INTENT(in) :: clouds
551	! Clouds microphysics condensible species properties file.
552	CHARACTER(len=*), INTENT(in) :: spcfile
553	! Minimum threshold for M0 of cloud condensation nuceli (m-3).
554	REAL(kind=mm_wp), INTENT(in) :: m0ccn_min
555	! Minimum threshold for the cloud drop radius (m).
556	REAL(kind=mm_wp), INTENT(in) :: drad_min
557
558	! Debug mode control flag.
559	LOGICAL, INTENT(in) :: debug
560
561	! Error status of the function.
562	TYPE(error) :: err
563
564	! Local variables:
565	INTEGER :: i
566	TYPE(cfgparser) :: cp
567	CHARACTER(len=st_slen), DIMENSION(:), ALLOCATABLE :: species
568
569	err = noerror
570
571	! Sanity check:
572	IF (mm_ini) THEN
573	err = error("mm_global_init: YAMMS global initialization already performed !",-1)
574	RETURN
575	ENDIF
576
577	! Free parameters:
578	mm_dt = dt
579	mm_df = df
580	mm_rm = rm
581	mm_rhoaer = rho_aer
582	mm_call_CH4hazeprod = call_CH4hazeprod
583	mm_p_prod = p_prod
584	mm_tx_prod = tx_prod
585	mm_rc_prod = rc_prod
586	mm_rpla = rplanet
587	mm_g0 = g0
588	mm_air_rad = air_rad
589	mm_air_mmol = air_mmol
590	mm_air_mmas = air_mmol / mm_navo
591
592	! Microphysical processes - Haze:
593	mm_coag_choice = coag_interactions
594	IF (mm_coag_choice < 0 .OR. mm_coag_choice > 7) THEN
595	err = error("mm_global_init: Invalid choice for coagulation interactions activation",-1)
596	RETURN
597	ENDIF
598	mm_call_hazeprod = call_hazeprod
599	mm_call_hazesed = call_hazesed
600	mm_call_hazecoag = call_hazecoag
601	mm_wsed_m0 = force_wsed_to_m0
602	mm_wsed_m3 = force_wsed_to_m3
603
604	! Moment threshold flags:
605	mm_m0as_min = MAX(0._mm_wp,m0as_min)
606	mm_rcs_min = MAX(1.e-9_mm_wp,rcs_min)
607	mm_m0af_min = MAX(0._mm_wp,m0af_min)
608	mm_rcf_min = MAX(mm_rm,rcf_min)
609
610	! Computes M3 thresholds from user-defined thresholds:
611	mm_m3as_min = mm_m0as_minmm_alpha_s(3._mm_wp) mm_rcs_min**3._mm_wp
612	mm_m3af_min = mm_m0af_minmm_alpha_f(3._mm_wp) mm_rcf_min**3._mm_wp
613
614	! Microphysical processes - Clouds:
615	mm_call_clouds = clouds
616
617	! Check clouds microphysics species file
618	IF (mm_call_clouds) THEN
619	IF (LEN_TRIM(spcfile) == 0) THEN
620	err = error("mm_global_init: No species properties file given",-1)
621	RETURN
622	ENDIF
623
624	! Reads species properties configuration file
625	err = cfg_read_config(cp,TRIM(spcfile))
626	IF (err /= 0) THEN
627	write(,) err
628	RETURN
629	ENDIF
630
631	! Reads used species
632	err = cfg_get_value(cp,"used_species",species)
633	IF (err /= 0) THEN
634	err = error("mm_global_init: cannot retrieve 'used_species' values",-1)
635	RETURN
636	ENDIF
637
638	mm_nesp = SIZE(species)
639	ALLOCATE(mm_spcname(mm_nesp),mm_xESPS(mm_nesp))
640
641	! Reads used species properties
642	DO i=1,mm_nesp
643	mm_spcname(i) = TRIM(species(i))
644	IF(.NOT.cfg_has_section(cp,TRIM(mm_spcname(i)))) THEN
645	err = error("mm_global_init: Cannot find "//TRIM(mm_spcname(i))//" properties",-1)
646	RETURN
647	ELSE
648	err = read_esp(cp,TRIM(mm_spcname(i)),mm_xESPS(i))
649	! Compute conversion factor: mol.mol-1 => kg.kg-1
650	mm_xESPS(i)%fmol2fmas = mm_xESPS(i)%masmol / mm_air_mmol
651	IF (err/=0) THEN
652	err = error("mm_global_init: "//TRIM(mm_spcname(i))//": "//TRIM(err%msg),-1)
653	RETURN
654	ENDIF
655	ENDIF
656	ENDDO
657
658	! Moment threshold flags:
659	mm_m0ccn_min = MAX(0._mm_wp,m0ccn_min)
660	mm_drad_min = MAX(1.e-9_mm_wp,drad_min)
661
662	! Computes M3 thresholds:
663	mm_m3ccn_min = mm_m0ccn_minmm_alpha_s(3._mm_wp) mm_rcs_min**3._mm_wp
664	mm_m3cld_min = mm_m0ccn_min * mm_drad_min**3._mm_wp
665	ENDIF ! end of mm_call_clouds
666
667	! Debug mode:
668	mm_debug = debug
669
670	! Computes conversion factor for bulk to apparent radius:
671	mm_rb2ra = mm_rm**((mm_df-3._mm_wp)/mm_df)
672
673	! Sanity check for settling velocity:
674	IF (mm_wsed_m0 .AND. mm_wsed_m3) THEN
675	err = error("'wsed_m0' and 'wsed_m3' options are mutually exclusive",-1)
676	ENDIF
677
678	! End of initialization:
679	mm_ini = err == noerror
680
681	END FUNCTION mm_global_init_0
682
683
684	FUNCTION mm_global_init_1(cfg) RESULT(err)
685	!! Set global configuration from a configuration file.
686	!!
687	!! @note:
688	!! See [[mm_globals(module):mm_global_init_0(function)]].
689	!!
690
691	! Configuration file path.
692	TYPE(cfgparser), INTENT(in) :: cfg
693	! Error status of the function.
694	TYPE(error) :: err
695
696	! Local variables:
697	INTEGER :: i
698	TYPE(cfgparser) :: spccfg
699	CHARACTER(len=st_slen) :: spcpath
700	CHARACTER(len=st_slen), DIMENSION(:), ALLOCATABLE :: species
701
702	err = noerror
703
704	! Sanity check:
705	IF (mm_ini) THEN
706	err = error("mm_global_init: YAMMS global initialization already performed !",-1)
707	RETURN
708	ENDIF
709
710	! Free parameters:
711	err = mm_check_opt(cfg_get_value(cfg,"timestep",mm_dt),mm_dt,wlog=mm_log)
712	IF (err/=0) RETURN
713	err = mm_check_opt(cfg_get_value(cfg,"df",mm_df),mm_df,wlog=mm_log)
714	IF (err/=0) RETURN
715	err = mm_check_opt(cfg_get_value(cfg,"rm",mm_rm),mm_rm,wlog=mm_log)
716	IF (err/=0) RETURN
717	err = mm_check_opt(cfg_get_value(cfg,"rho_aer",mm_rhoaer),mm_rhoaer,wlog=mm_log)
718	IF (err/=0) RETURN
719	err = mm_check_opt(cfg_get_value(cfg,"call_haze_prod",mm_call_CH4hazeprod),mm_call_CH4hazeprod,wlog=mm_log)
720	IF (err/=0) RETURN
721	err = mm_check_opt(cfg_get_value(cfg,"p_prod",mm_p_prod),mm_p_prod,wlog=mm_log)
722	IF (err/=0) RETURN
723	err = mm_check_opt(cfg_get_value(cfg,"tx_prod",mm_tx_prod),mm_tx_prod,wlog=mm_log)
724	IF (err/=0) RETURN
725	err = mm_check_opt(cfg_get_value(cfg,"rc_prod",mm_rc_prod),mm_rc_prod,wlog=mm_log)
726	IF (err/=0) RETURN
727	err = mm_check_opt(cfg_get_value(cfg,"planet_radius",mm_rpla),mm_rpla,wlog=mm_log)
728	IF (err/=0) RETURN
729	err = mm_check_opt(cfg_get_value(cfg,"g0",mm_g0),mm_g0,wlog=mm_log)
730	IF (err/=0) RETURN
731	err = mm_check_opt(cfg_get_value(cfg,"air_radius",mm_air_rad),mm_air_rad,wlog=mm_log)
732	IF (err/=0) RETURN
733	err = mm_check_opt(cfg_get_value(cfg,"air_molarmass",mm_air_mmol),mm_air_mmol,wlog=mm_log)
734	IF (err/=0) RETURN
735	err = mm_check_opt(cfg_get_value(cfg,"air_molecularmass",mm_air_mmas),mm_air_mmas,wlog=mm_log)
736	IF (err/=0) RETURN
737
738	! Microphysical processes - Haze:
739	err = mm_check_opt(cfg_get_value(cfg,"haze_coag_interactions",mm_coag_choice),mm_coag_choice,wlog=mm_log)
740	IF (err/=0) RETURN
741	IF (mm_coag_choice < 0 .OR. mm_coag_choice > 7) THEN
742	err = error("mm_global_init: Invalid choice for coagulation interactions activation",-1)
743	RETURN
744	ENDIF
745	err = mm_check_opt(cfg_get_value(cfg,"haze_production",mm_call_hazeprod),mm_call_hazeprod,wlog=mm_log)
746	IF (err/=0) RETURN
747	err = mm_check_opt(cfg_get_value(cfg,"haze_sedimentation",mm_call_hazesed),mm_call_hazesed,wlog=mm_log)
748	IF (err/=0) RETURN
749	err = mm_check_opt(cfg_get_value(cfg,"haze_coagulation",mm_call_hazecoag),mm_call_hazecoag,wlog=mm_log)
750	IF (err/=0) RETURN
751	err = mm_check_opt(cfg_get_value(cfg,"wsed_m0",mm_wsed_m0),mm_wsed_m0,wlog=mm_log)
752	IF (err/=0) RETURN
753	err = mm_check_opt(cfg_get_value(cfg,"wsed_m3",mm_wsed_m3),mm_wsed_m3,wlog=mm_log)
754	IF (err/=0) RETURN
755
756	! Moment threshold flags:
757	err = mm_check_opt(cfg_get_value(cfg,"m0as_min",mm_m0as_min),mm_m0as_min,wlog=mm_log)
758	IF (err/=0) RETURN
759	err = mm_check_opt(cfg_get_value(cfg,"rcs_min",mm_rcs_min),mm_rcs_min,wlog=mm_log)
760	IF (err/=0) RETURN
761	err = mm_check_opt(cfg_get_value(cfg,"m0af_min",mm_m0af_min),mm_m0af_min,wlog=mm_log)
762	IF (err/=0) RETURN
763	err = mm_check_opt(cfg_get_value(cfg,"rcf_min",mm_rcf_min),mm_rcf_min,wlog=mm_log)
764	IF (err/=0) RETURN
765
766	! Microphysical processes - Clouds:
767	err = mm_check_opt(cfg_get_value(cfg,"clouds_microphysics",mm_call_clouds),mm_call_clouds,wlog=mm_log)
768	IF (err/=0) RETURN
769
770	! Check clouds microphysics input
771	IF (mm_call_clouds) THEN
772	! Gets species property file path
773	err = cfg_get_value(cfg,'species_cfg',spcpath) ; IF (err /= 0) RETURN
774	! Reads species properties configuration file
775	err = cfg_read_config(spccfg,trim(spcpath)) ; IF (err /= 0) RETURN
776	err = cfg_get_value(spccfg,"used_species",species)
777	IF (err /= 0) THEN
778	err = error("mm_global_init: cannot retrieve 'used_species' values",-1)
779	RETURN
780	ENDIF
781	mm_nesp = SIZE(species)
782	ALLOCATE(mm_spcname(mm_nesp),mm_xESPS(mm_nesp))
783	DO i=1,mm_nesp
784	mm_spcname(i) = TRIM(species(i))
785	IF (.NOT.cfg_has_section(spccfg,TRIM(mm_spcname(i)))) THEN
786	err = error("mm_global_init: Cannot find "//TRIM(mm_spcname(i))//" properties",-1)
787	RETURN
788	ELSE
789	err = read_esp(spccfg,TRIM(mm_spcname(i)),mm_xESPS(i))
790	! Compute conversion factor: mol.mol-1 => kg.kg-1
791	mm_xESPS(i)%fmol2fmas = mm_xESPS(i)%masmol / mm_air_mmol
792	IF (err/=0) THEN
793	err = error(TRIM(mm_spcname(i))//": "//TRIM(err%msg),-2)
794	RETURN
795	ENDIF
796	ENDIF
797	ENDDO
798	ENDIF
799
800	! Debug mode:
801	err = mm_check_opt(cfg_get_value(cfg,"debug",mm_debug),mm_debug,wlog=mm_log)
802	IF (err/=0) RETURN
803
804	! Computes M3 thresholds from user-defined thresholds:
805	mm_m0as_min = MAX(0._mm_wp,mm_m0as_min)
806	mm_rcs_min = MAX(1.e-10_mm_wp,mm_rcs_min)
807	mm_m0af_min = MAX(0._mm_wp,mm_m0af_min)
808	mm_rcf_min = MAX(mm_rm,mm_rcf_min)
809	mm_m3as_min = mm_m0as_minmm_alpha_s(3._mm_wp) mm_rcs_min**3._mm_wp
810	mm_m3af_min = mm_m0af_minmm_alpha_f(3._mm_wp) mm_rcf_min**3._mm_wp
811
812	! Computes conversion factor for bulk to apparent radius:
813	mm_rb2ra = mm_rm**((mm_df-3._mm_wp)/mm_df)
814
815	! Sanity check for settling velocity:
816	IF (mm_wsed_m0 .AND. mm_wsed_m3) THEN
817	err = error("'wsed_m0' and 'wsed_m3' options are mutually exclusive",-1)
818	ENDIF
819
820	! End of initialization:
821	mm_ini = err == noerror
822
823	END FUNCTION mm_global_init_1
824
825
826	FUNCTION mm_column_init(plev,zlev,play,zlay,temp) RESULT(err)
827	!! Initialize vertical atmospheric fields.
828	!!
829	!! This subroutine initializes vertical fields needed by the microphysics:
830	!! 1. Save reversed input field into "local" array
831	!! 2. Compute thicknesses layers and levels
832	!! 3. Interpolate temperature at levels
833	!!
834	!! @warning
835	!! The method should be called whenever the vertical structure of the atmosphere changes.
836	!! All the input vectors should be defined from __GROUND__ to __TOP__ of the atmosphere,
837	!! otherwise nasty things will occur in computations.
838	!!
839	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: plev ! Pressure levels (Pa).
840	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: zlev ! Altitude levels (m).
841	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: play ! Pressure layers (Pa).
842	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: zlay ! Altitude at the center of each layer (m).
843	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: temp ! Temperature at the center of each layer (K).
844	TYPE(error) :: err ! Error status of the function.
845	INTEGER :: i
846
847	err = noerror
848	mm_ini_col = .false.
849
850	! Global initialization must be done before:
851	IF (.NOT.mm_ini) THEN
852	err = error("mm_column_init: Global initialization not done yet",-1)
853	RETURN
854	ENDIF
855
856	! Check number of vertical layers:
857	IF (mm_nla < 0) THEN
858	mm_nla = SIZE(play)
859	ELSE
860	IF (mm_nla /= SIZE(play)) THEN
861	err = error("mm_column_init: mm_nla cannot be modified dynamically within the run",-1)
862	RETURN
863	ENDIF
864	ENDIF
865
866	! Check number of vertical levels:
867	IF (mm_nle < 0) THEN
868	mm_nle = SIZE(plev)
869	ELSE
870	IF (mm_nle /= SIZE(plev)) THEN
871	err = error("mm_column_init: mm_nle cannot be modified dynamically within the run",-1)
872	RETURN
873	ENDIF
874	ENDIF
875
876	! Sanity check:
877	IF (mm_nla+1 /= mm_nle) THEN
878	err = error("mm_column_init: Inconsistent number of layers/levels",-1)
879	RETURN
880	ENDIF
881
882	! Allocates if required:
883	IF (.NOT.ALLOCATED(mm_plev)) ALLOCATE(mm_plev(mm_nle))
884	IF (.NOT.ALLOCATED(mm_zlev)) ALLOCATE(mm_zlev(mm_nle))
885	IF (.NOT.ALLOCATED(mm_play)) ALLOCATE(mm_play(mm_nla))
886	IF (.NOT.ALLOCATED(mm_zlay)) ALLOCATE(mm_zlay(mm_nla))
887	IF (.NOT.ALLOCATED(mm_temp)) ALLOCATE(mm_temp(mm_nla))
888	IF (.NOT.ALLOCATED(mm_btemp)) ALLOCATE(mm_btemp(mm_nle))
889	IF (.NOT.ALLOCATED(mm_dzlev)) ALLOCATE(mm_dzlev(mm_nla))
890	IF (.NOT.ALLOCATED(mm_dzlay)) ALLOCATE(mm_dzlay(mm_nla))
891	IF (.NOT.ALLOCATED(mm_rhoair)) ALLOCATE(mm_rhoair(mm_nla))
892
893	! Saves reversed input vectors:
894	mm_plev = plev(mm_nle:1:-1)
895	mm_zlev = zlev(mm_nle:1:-1)
896	mm_play = play(mm_nla:1:-1)
897	mm_zlay = zlay(mm_nla:1:-1)
898	mm_temp = temp(mm_nla:1:-1)
899
900	! Computes temperature vertical profil at interfaces:
901	mm_btemp(2:mm_nla) = (mm_temp(1:mm_nla-1) + mm_temp(2:mm_nla)) / 2._mm_wp
902	mm_btemp(1) = mm_temp(1)
903	mm_btemp(mm_nle) = mm_temp(mm_nla) + (mm_temp(mm_nla) - mm_temp(mm_nla-1)) / 2._mm_wp
904
905	! Computes atmospheric levels thickness:
906	mm_dzlev(1:mm_nla) = mm_zlev(1:mm_nle-1)-mm_zlev(2:mm_nle)
907
908	! Computes atmospheric layers thickness :
909	mm_dzlay(1:mm_nla-1) = mm_zlay(1:mm_nla-1)-mm_zlay(2:mm_nla)
910	mm_dzlay(mm_nla) = mm_dzlay(mm_nla-1)
911
912	! Hydrostatic equilibrium:
913	mm_rhoair(1:mm_nla) = (mm_plev(2:mm_nle)-mm_plev(1:mm_nla)) / (mm_effg(mm_zlay)*mm_dzlev)
914
915	! Write out profiles for debug and log:
916	IF (mm_log.AND.mm_debug) THEN
917	WRITE(*,'(a)') '# TEMP PLAY ZLAY DZLAY RHOAIR'
918	DO i=1,mm_nla
919	WRITE(*,'(5(ES15.7,2X))') mm_temp(i),mm_play(i),mm_zlay(i),mm_dzlay(i), mm_rhoair(i)
920	ENDDO
921	WRITE(*,'(a)') '# TEMP PLEV ZLEV DZLEV'
922	DO i=1,mm_nle
923	IF (i /= mm_nle) THEN
924	WRITE(*,'(4(ES15.7,2X))') mm_btemp(i),mm_plev(i),mm_zlev(i),mm_dzlev(i)
925	ELSE
926	WRITE(*,'(3(ES15.7,2X))') mm_btemp(i),mm_plev(i),mm_zlev(i)
927	ENDIF
928	ENDDO
929	ENDIF
930
931	! End of initialization:
932	mm_ini_col = .true.
933
934	RETURN
935	END FUNCTION mm_column_init
936
937
938	FUNCTION mm_aerosols_init(m0aer_s,m3aer_s,m0aer_f,m3aer_f) RESULT(err)
939	!! Initialize aerosol tracers vertical grid.
940	!!
941	!! The subroutine initializes aerosols microphysics tracers columns. It allocates variables if
942	!! required and stores input vectors in reversed order. It also computes the characteristic radii
943	!! of each mode.
944	!!
945	!! @warning
946	!! The method should be called after mm_global_init and mm_column_init. Moreover, it should be called
947	!! whenever the vertical structure of the atmosphere changes.
948	!! All the input vectors should be defined from __GROUND__ to __TOP__ of the atmosphere,
949	!! otherwise nasty things will occur in computations.
950	!!
951	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m0aer_s ! 0th order moment of the spherical mode (m-2).
952	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m3aer_s ! 3rd order moment of the spherical mode (m3.m-2).
953	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m0aer_f ! 0th order moment of the fractal mode (m-2).
954	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m3aer_f ! 3rd order moment of the fractal mode (m3.m-2).
955	TYPE(error) :: err ! Error status of the function.
956
957	err = noerror
958
959	! Global initialization must be done before:
960	IF (.NOT.mm_ini) THEN
961	err = error("mm_aerosols_init: Global initialization not done yet",-1) ; RETURN
962	ENDIF
963
964	! Column initialization must be done before:
965	IF (.NOT.mm_ini_col) THEN
966	err = error("mm_aerosols_init: Column initialization not done yet",-1) ; RETURN
967	ENDIF
968
969	! Sanity check:
970	IF (SIZE(m0aer_s) /= mm_nla) THEN
971	err = error("mm_aerosols_init: Invalid size for input arrays",-1) ; RETURN
972	ENDIF
973
974	! Allocate variable if required:
975	IF (.NOT.ALLOCATED(mm_m0aer_s)) ALLOCATE(mm_m0aer_s(mm_nla))
976	IF (.NOT.ALLOCATED(mm_m3aer_s)) ALLOCATE(mm_m3aer_s(mm_nla))
977	IF (.NOT.ALLOCATED(mm_m0aer_f)) ALLOCATE(mm_m0aer_f(mm_nla))
978	IF (.NOT.ALLOCATED(mm_m3aer_f)) ALLOCATE(mm_m3aer_f(mm_nla))
979	IF (.NOT.ALLOCATED(mm_rcs)) ALLOCATE(mm_rcs(mm_nla))
980	IF (.NOT.ALLOCATED(mm_rcf)) ALLOCATE(mm_rcf(mm_nla))
981
982	! Allocate memory for diagnostics if required:
983	IF (.NOT.ALLOCATED(mm_m0as_vsed)) THEN
984	ALLOCATE(mm_m0as_vsed(mm_nla))
985	mm_m0as_vsed(:) = 0._mm_wp
986	ENDIF
987	IF (.NOT.ALLOCATED(mm_m3as_vsed)) THEN
988	ALLOCATE(mm_m3as_vsed(mm_nla))
989	mm_m3as_vsed(:) = 0._mm_wp
990	ENDIF
991	IF (.NOT.ALLOCATED(mm_m0af_vsed)) THEN
992	ALLOCATE(mm_m0af_vsed(mm_nla))
993	mm_m0af_vsed(:) = 0._mm_wp
994	ENDIF
995	IF (.NOT.ALLOCATED(mm_m3af_vsed)) THEN
996	ALLOCATE(mm_m3af_vsed(mm_nla))
997	mm_m3af_vsed(:) = 0._mm_wp
998	ENDIF
999	IF (.NOT.ALLOCATED(mm_aer_s_flux)) THEN
1000	ALLOCATE(mm_aer_s_flux(mm_nla))
1001	mm_aer_s_flux(:) = 0._mm_wp
1002	ENDIF
1003	IF (.NOT.ALLOCATED(mm_aer_f_flux)) THEN
1004	ALLOCATE(mm_aer_f_flux(mm_nla))
1005	mm_aer_f_flux(:) = 0._mm_wp
1006	ENDIF
1007
1008	! Initialization of aerosol tracers:
1009	! @note: mm_dzlev is already from top to ground
1010	mm_m0aer_s = m0aer_s(mm_nla:1:-1) / mm_dzlev(:)
1011	mm_m3aer_s = m3aer_s(mm_nla:1:-1) / mm_dzlev(:)
1012	mm_m0aer_f = m0aer_f(mm_nla:1:-1) / mm_dzlev(:)
1013	mm_m3aer_f = m3aer_f(mm_nla:1:-1) / mm_dzlev(:)
1014
1015	! Setup threshold (useful for preventing bugs):
1016	call mm_set_moments_thresholds()
1017
1018	! Initialization of spherical aerosol characteristic radii:
1019	WHERE(mm_m3aer_s > 0._mm_wp .AND. mm_m0aer_s > 0._mm_wp)
1020	mm_rcs = mm_get_rcs(mm_m0aer_s,mm_m3aer_s)
1021	ELSEWHERE
1022	mm_rcs = 0._mm_wp
1023	ENDWHERE
1024
1025	! Initialization of fractal aerosol characteristic radii:
1026	WHERE(mm_m3aer_f > 0._mm_wp .AND. mm_m0aer_f > 0._mm_wp)
1027	mm_rcf = mm_get_rcf(mm_m0aer_f,mm_m3aer_f)
1028	ELSEWHERE
1029	mm_rcf = 0._mm_wp
1030	ENDWHERE
1031
1032	! End of initialization:
1033	mm_ini_aer = .true.
1034
1035	END FUNCTION mm_aerosols_init
1036
1037	FUNCTION mm_clouds_init(m0ccn,m3ccn,m3ice,gas) RESULT(err)
1038	!! Initialize clouds tracers vertical grid.
1039	!!
1040	!! The subroutine initializes cloud microphysics tracers columns. It allocates variables if
1041	!! required and stores input vectors in reversed order. It also computes the mean drop radius
1042	!! and density and allocates diagnostic vectors.
1043	!!
1044	!! @note
1045	!! All the input arguments should be defined from __GROUND__ to __TOP__.
1046	!!
1047	!! @warning
1048	!! [[mm_global_init(interface)]] and [[mm_column_init(function)]] must have been called at least once before this method is called.
1049	!! Moreover, this method should be after each call of [[mm_column_init(function)]] to reflect changes in the
1050	!! vertical atmospheric structure.
1051	!!
1052	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m0ccn ! 0th order moment of the CCN distribution (m-2).
1053	REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m3ccn ! 3rd order moment of the CCN distribution (m3.m-2).
1054	REAL(kind=mm_wp), DIMENSION(:,:), INTENT(in) :: m3ice ! 3rd order moments of the ice components (m3.m-2).
1055	REAL(kind=mm_wp), DIMENSION(:,:), INTENT(in) :: gas ! Condensible species gas molar fraction (mol.mol-1).
1056
1057	INTEGER :: i
1058	TYPE(error) :: err ! Error status of the function.
1059
1060	err = noerror
1061	IF (.NOT.mm_ini) THEN
1062	err = error("Global initialization not done yet",-8)
1063	RETURN
1064	ENDIF
1065
1066	IF (.NOT.mm_call_clouds) THEN
1067	IF (mm_debug) WRITE(*,'(a)') "[MM_DEBUG - mm_clouds_init] Cloud microphysic is not enabled..."
1068	RETURN
1069	ENDIF
1070
1071	! Allocate variable if required:
1072	IF (.NOT.ALLOCATED(mm_m0ccn)) ALLOCATE(mm_m0ccn(mm_nla))
1073	IF (.NOT.ALLOCATED(mm_m3ccn)) ALLOCATE(mm_m3ccn(mm_nla))
1074	IF (.NOT.ALLOCATED(mm_m3ice)) ALLOCATE(mm_m3ice(mm_nla,mm_nesp))
1075	IF (.NOT.ALLOCATED(mm_gas)) ALLOCATE(mm_gas(mm_nla,mm_nesp))
1076	IF (.NOT.ALLOCATED(mm_drad)) ALLOCATE(mm_drad(mm_nla))
1077	IF (.NOT.ALLOCATED(mm_drho)) ALLOCATE(mm_drho(mm_nla))
1078
1079	! Allocate memory for diagnostics:
1080	mm_ccn_prec = 0._mm_wp
1081	IF (.NOT.ALLOCATED(mm_ice_prec)) THEN
1082	ALLOCATE(mm_ice_prec(mm_nesp)) ; mm_ice_prec(:) = 0._mm_wp
1083	ENDIF
1084	IF (.NOT.ALLOCATED(mm_cld_vsed)) THEN
1085	ALLOCATE(mm_cld_vsed(mm_nla)) ; mm_cld_vsed(:) = 0._mm_wp
1086	ENDIF
1087	IF (.NOT.ALLOCATED(mm_ccn_flux)) THEN
1088	ALLOCATE(mm_ccn_flux(mm_nla)) ; mm_ccn_flux(:) = 0._mm_wp
1089	ENDIF
1090	IF (.NOT.ALLOCATED(mm_ice_fluxes)) THEN
1091	ALLOCATE(mm_ice_fluxes(mm_nla,mm_nesp)) ; mm_ice_fluxes(:,:) = 0._mm_wp
1092	ENDIF
1093	IF (.NOT.ALLOCATED(mm_gas_sat)) THEN
1094	ALLOCATE(mm_gas_sat(mm_nla,mm_nesp)) ; mm_gas_sat(:,:) = 0._mm_wp
1095	ENDIF
1096	IF (.NOT.ALLOCATED(mm_nrate)) THEN
1097	ALLOCATE(mm_nrate(mm_nla,mm_nesp)) ; mm_nrate(:,:) = 0._mm_wp
1098	ENDIF
1099	IF (.NOT.ALLOCATED(mm_grate)) THEN
1100	ALLOCATE(mm_grate(mm_nla,mm_nesp)) ; mm_grate(:,:) = 0._mm_wp
1101	ENDIF
1102
1103	! /!\ mm_dzlev already from top to ground
1104	mm_m0ccn = m0ccn(mm_nla:1:-1) / mm_dzlev(:)
1105	mm_m3ccn = m3ccn(mm_nla:1:-1) / mm_dzlev(:)
1106	DO i = 1, mm_nesp
1107	mm_m3ice(:,i) = m3ice(mm_nla:1:-1,i) / mm_dzlev(:)
1108	mm_gas(:,i) = gas(mm_nla:1:-1,i)
1109	ENDDO
1110
1111	! Setup threshold:
1112	call mm_set_moments_cld_thresholds()
1113
1114	! Drop mean radius and density:
1115	call mm_cloud_properties(mm_m0ccn,mm_m3ccn,mm_m3ice,mm_drad,mm_drho)
1116
1117	! End of initialization:
1118	mm_ini_cld = .true.
1119
1120	END FUNCTION mm_clouds_init
1121
1122
1123	!============================================================================
1124	! INTER-MOMENT RELATION METHODS
1125	!============================================================================
1126
1127	PURE FUNCTION mm_alpha_s(k) RESULT (res)
1128	!! Inter-moment relation for spherical aerosols size distribution law.
1129	!! Mk / M0 = rc^k . alpha(k)
1130	!!
1131	REAL(kind=mm_wp), INTENT(in) :: k ! k order of the moment.
1132	REAL(kind=mm_wp) :: sigma ! Standard deviation.
1133	REAL(kind=mm_wp) :: res ! Alpha value.
1134
1135	! Titan's case
1136	!~~~~~~~~~~~~~
1137	! res = SUM(dexp(mm_asp%ak2 + mm_asp%bk + mm_asp%c))
1138
1139	! Pluto's case
1140	!~~~~~~~~~~~~~
1141	sigma = 0.2_mm_wp
1142	res = exp(k*2 sigma**2 / 2._mm_wp)
1143
1144	RETURN
1145	END FUNCTION mm_alpha_s
1146
1147
1148	PURE FUNCTION mm_alpha_f(k) RESULT (res)
1149	!! Inter-moment relation for fractal aerosols size distribution law.
1150	!! Mk / M0 = rc^k . alpha(k)
1151	!!
1152	REAL(kind=mm_wp), INTENT(in) :: k ! k order of the moment.
1153	REAL(kind=mm_wp) :: sigma ! Standard deviation.
1154	REAL(kind=mm_wp) :: res ! Alpha value.
1155
1156	! Titan's case
1157	!~~~~~~~~~~~~~
1158	! res = SUM(dexp(mm_afp%ak2 + mm_afp%bk + mm_afp%c))
1159
1160	! Pluto's case
1161	!~~~~~~~~~~~~~
1162	sigma = 0.35_mm_wp
1163	res = exp(k*2 sigma**2 / 2._mm_wp)
1164
1165	RETURN
1166	END FUNCTION mm_alpha_f
1167
1168	PURE FUNCTION mm_alpha_ccn(k) RESULT (res)
1169	!! Inter-moment relation for fractal aerosols size distribution law.
1170	!! Mk / M0 = rc^k . alpha(k)
1171	!!
1172	REAL(kind=mm_wp), INTENT(in) :: k ! k order of the moment.
1173	REAL(kind=mm_wp) :: sigma ! Standard deviation.
1174	REAL(kind=mm_wp) :: res ! Alpha value.
1175
1176	! Pluto's case
1177	!~~~~~~~~~~~~~
1178	sigma = 0.2_mm_wp
1179	res = exp(k*2 sigma**2 / 2._mm_wp)
1180
1181	RETURN
1182	END FUNCTION mm_alpha_ccn
1183
1184
1185	!============================================================================
1186	! HAZE RELATED METHODS
1187	!============================================================================
1188
1189	SUBROUTINE mm_set_moments_thresholds()
1190	!! Apply minimum threshold for the aerosols moments.
1191	!!
1192	!! The method resets moments (for both modes and orders, 0 and 3) values to zero if
1193	!! their current value is below the minimum threholds.
1194	!!
1195	INTEGER :: i
1196	DO i=1,mm_nla
1197	IF ((mm_m0aer_s(i) < mm_m0as_min) .OR. (mm_m3aer_s(i) < mm_m3as_min)) THEN
1198	mm_m0aer_s(i) = 0._mm_wp
1199	mm_m3aer_s(i) = 0._mm_wp
1200	ENDIF
1201	IF ((mm_m0aer_f(i) < mm_m0af_min) .OR. (mm_m3aer_f(i) < mm_m3af_min)) THEN
1202	mm_m0aer_f(i) = 0._mm_wp
1203	mm_m3aer_f(i) = 0._mm_wp
1204	ENDIF
1205	ENDDO
1206	END SUBROUTINE mm_set_moments_thresholds
1207
1208
1209	ELEMENTAL FUNCTION mm_get_rcs(m0,m3) RESULT(res)
1210	!! Get the characteristic radius for the spherical aerosols size distribution.
1211	!!
1212	!! The method computes the characteristic radius of the spherical aerosol size distribution
1213	!! law according to its moments and its inter-moments relation.
1214	!!
1215	REAL(kind=mm_wp), INTENT(in) :: m0 ! 0th order moment
1216	REAL(kind=mm_wp), INTENT(in) :: m3 ! 3rd order moment
1217	REAL(kind=mm_wp) :: res ! Radius
1218	res = (m3 / (m0mm_alpha_s(3._mm_wp)))*(1._mm_wp/3._mm_wp)
1219	END FUNCTION mm_get_rcs
1220
1221
1222	ELEMENTAL FUNCTION mm_get_rcf(m0,m3) RESULT(res)
1223	!! Get the characteristic radius for the fractal aerosols size distribution.
1224	!!
1225	!! The method computes the characteristic radius of the fractal aerosol size distribution
1226	!! law according to its moments and its inter-moments relation.
1227	!!
1228	REAL(kind=mm_wp), INTENT(in) :: m0 ! 0th order moment
1229	REAL(kind=mm_wp), INTENT(in) :: m3 ! 3rd order moment
1230	REAL(kind=mm_wp) :: res ! Radius
1231	res = (m3 / (m0mm_alpha_f(3._mm_wp)))*(1._mm_wp/3._mm_wp)
1232	END FUNCTION mm_get_rcf
1233
1234
1235	!============================================================================
1236	! CLOUD RELATED METHODS
1237	!============================================================================
1238
1239	SUBROUTINE mm_set_moments_cld_thresholds()
1240	!! Apply minimum threshold for the cloud drop moments.
1241	!! The method resets moments values to zero if their current value is below the minimum threholds.
1242	!!
1243	INTEGER :: i, j
1244	REAL(kind=mm_wp) :: m3cld = 0._mm_wp
1245
1246	DO i = 1, mm_nla
1247	m3cld = mm_m3ccn(i)
1248	DO j = 1, mm_nesp
1249	m3cld = m3cld + mm_m3ice(i,j)
1250	ENDDO
1251
1252	IF ((mm_m0ccn(i) < mm_m0ccn_min) .OR. (mm_m3ccn(i) < mm_m3ccn_min) .OR. (m3cld < mm_m3cld_min)) THEN
1253	mm_m0ccn(i) = 0._mm_wp
1254	mm_m3ccn(i) = 0._mm_wp
1255	DO j = 1, mm_nesp
1256	mm_m3ice(i,j) = 0._mm_wp
1257	ENDDO
1258	ENDIF
1259	ENDDO
1260	END SUBROUTINE mm_set_moments_cld_thresholds
1261
1262
1263	SUBROUTINE cldprop_sc(m0ccn,m3ccn,m3ice,drad,drho)
1264	!! Get cloud drop properties (scalar).
1265	!! The method computes the mean radius and mean density of cloud drops.
1266	!!
1267	!! @warning
1268	!! If __drad__ is greater than __drmax__ it is automatically set to __drmax__, but computation of
1269	!! __drho__ remains unmodified. So __drho__ is not correct in that case!
1270	!!
1271	REAL(kind=mm_wp), INTENT(in) :: m0ccn ! 0th order moment of the ccn
1272	REAL(kind=mm_wp), INTENT(in) :: m3ccn ! 3rd order moment of the ccn
1273	REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: m3ice ! 3rd order moments of each ice component
1274	REAL(kind=mm_wp), INTENT(out) :: drad ! Output mean drop radius
1275	REAL(kind=mm_wp), INTENT(out), OPTIONAL :: drho ! Optional output mean drop density
1276
1277	REAL(kind=mm_wp) :: Ntot, Vtot, Wtot
1278	REAL(kind=mm_wp), PARAMETER :: athird = 1._mm_wp / 3._mm_wp
1279	REAL(kind=mm_wp), PARAMETER :: pifac = (4._mm_wp * mm_pi) / 3._mm_wp
1280
1281	! Set to zero :
1282	drad = 0._mm_wp
1283	IF (PRESENT(drho)) drho = 0._mm_wp
1284
1285	! Initialization :
1286	Ntot = m0ccn
1287	Vtot = m3ccn + SUM(m3ice)
1288	Wtot = m3ccnmm_rhoaer + SUM(m3icemm_xESPS(:)%rho_s)
1289
1290	IF (Ntot <= mm_m0ccn_min .OR. Vtot <= mm_m3cld_min) THEN
1291	drad = 0._mm_wp
1292	IF (PRESENT(drho)) drho = 0._mm_wp
1293	ELSE
1294	drad = (Vtot / Ntot)**athird
1295	drad = MAX(MIN(drad,mm_drad_max),mm_drad_min)
1296	IF (PRESENT(drho)) drho = Wtot / Vtot
1297	ENDIF
1298
1299	RETURN
1300	END SUBROUTINE cldprop_sc
1301
1302	SUBROUTINE cldprop_ve(m0ccn,m3ccn,m3ice,drad,drho)
1303	!! Get cloud drop properties (vector).
1304	!!
1305	!! The method performs the same computations than [[cldprop_sc(subroutine)]]
1306	!! but for the entire vertical atmospheric structure.
1307	!!
1308	REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: m0ccn ! 0th order moment of the ccn.
1309	REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: m3ccn ! 3rd order moment of the ccn.
1310	REAL(kind=mm_wp), INTENT(in), DIMENSION(:,:) :: m3ice ! 3rd order moments of each ice component.
1311	REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: drad ! Output mean drop radius.
1312	REAL(kind=mm_wp), INTENT(out), DIMENSION(:), OPTIONAL :: drho ! Optional output mean drop density.
1313
1314	INTEGER :: i
1315
1316	IF (PRESENT(drho)) THEN
1317	DO i = 1, SIZE(m0ccn)
1318	call cldprop_sc(m0ccn(i),m3ccn(i),m3ice(i,:),drad(i),drho(i))
1319	ENDDO
1320	ELSE
1321	DO i = 1, SIZE(m0ccn)
1322	call cldprop_sc(m0ccn(i),m3ccn(i),m3ice(i,:),drad(i))
1323	ENDDO
1324	ENDIF
1325
1326	RETURN
1327	END SUBROUTINE cldprop_ve
1328
1329
1330	FUNCTION read_esp(parser,sec,pp) RESULT (err)
1331	!! Read and store [[mm_esp(type)]] parameters.
1332	!!
1333	TYPE(cfgparser), INTENT(in) :: parser ! Configuration parser.
1334	CHARACTER(len=*), INTENT(in) :: sec ! Name of the species.
1335	TYPE(mm_esp), INTENT(out) :: pp ! [[mm_esp(type)]] object that stores the parameters.
1336
1337	TYPE(error) :: err ! Error status of the function.
1338
1339	err = cfg_get_value(parser,TRIM(sec)//'/name',pp%name) ; IF (err /= 0) RETURN
1340	err = cfg_get_value(parser,TRIM(sec)//'/mas',pp%mas) ; IF (err /= 0) RETURN
1341	err = cfg_get_value(parser,TRIM(sec)//'/vol',pp%vol) ; IF (err /= 0) RETURN
1342	err = cfg_get_value(parser,TRIM(sec)//'/ray',pp%ray) ; IF (err /= 0) RETURN
1343	err = cfg_get_value(parser,TRIM(sec)//'/masmol',pp%masmol) ; IF (err /= 0) RETURN
1344	err = cfg_get_value(parser,TRIM(sec)//'/rho_l',pp%rho_l) ; IF (err /= 0) RETURN
1345	err = cfg_get_value(parser,TRIM(sec)//'/rho_s',pp%rho_s) ; IF (err /= 0) RETURN
1346	err = cfg_get_value(parser,TRIM(sec)//'/Tc',pp%Tc) ; IF (err /= 0) RETURN
1347	err = cfg_get_value(parser,TRIM(sec)//'/pc',pp%pc) ; IF (err /= 0) RETURN
1348	err = cfg_get_value(parser,TRIM(sec)//'/Tb',pp%Tb) ; IF (err /= 0) RETURN
1349	err = cfg_get_value(parser,TRIM(sec)//'/w',pp%w) ; IF (err /= 0) RETURN
1350	err = cfg_get_value(parser,TRIM(sec)//'/a0_sat',pp%a0_sat) ; IF (err /= 0) RETURN
1351	err = cfg_get_value(parser,TRIM(sec)//'/a1_sat',pp%a1_sat) ; IF (err /= 0) RETURN
1352	err = cfg_get_value(parser,TRIM(sec)//'/a2_sat',pp%a2_sat) ; IF (err /= 0) RETURN
1353	err = cfg_get_value(parser,TRIM(sec)//'/a3_sat',pp%a3_sat) ; IF (err /= 0) RETURN
1354	err = cfg_get_value(parser,TRIM(sec)//'/a4_sat',pp%a4_sat) ; IF (err /= 0) RETURN
1355	err = cfg_get_value(parser,TRIM(sec)//'/a5_sat',pp%a5_sat) ; IF (err /= 0) RETURN
1356	err = cfg_get_value(parser,TRIM(sec)//'/a6_sat',pp%a6_sat) ; IF (err /= 0) RETURN
1357	err = cfg_get_value(parser,TRIM(sec)//'/mteta',pp%mteta) ; IF (err /= 0) RETURN
1358	err = cfg_get_value(parser,TRIM(sec)//'/fdes',pp%fdes) ; IF (err /= 0) RETURN
1359	err = cfg_get_value(parser,TRIM(sec)//'/fdif',pp%fdif) ; IF (err /= 0) RETURN
1360	err = cfg_get_value(parser,TRIM(sec)//'/nus',pp%nus) ; IF (err /= 0) RETURN
1361
1362	RETURN
1363	END FUNCTION read_esp
1364
1365
1366	!============================================================================
1367	! MISCELLANEOUS METHODS
1368	!============================================================================
1369
1370	ELEMENTAL FUNCTION mm_effg(z) RESULT(effg)
1371	!! Compute effective gravitational acceleration.
1372	!!
1373	REAL(kind=mm_wp), INTENT(in) :: z ! Altitude (m)
1374	REAL(kind=mm_wp) :: effg ! Effective gravitational acceleration (m.s-2)
1375	effg = mm_g0
1376	IF (mm_use_effg) THEN
1377	effg = effg * (mm_rpla/(mm_rpla+z))**2
1378	ENDIF
1379	RETURN
1380	END FUNCTION mm_effg
1381
1382
1383	SUBROUTINE mm_dump_parameters()
1384	!! Dump global parameters on stdout.
1385	!!
1386	WRITE(*,'(a)') "========= YAMMS PARAMETERS ============"
1387	WRITE(*,'(a,a)') "mm_fp_precision : ", mm_wp_s
1388	WRITE(*,'(a,L2)') "mm_debug : ", mm_debug
1389	WRITE(*,'(a)') "---------------------------------------"
1390	WRITE(*,'(a)') "Microphysical control flags"
1391
1392	! Haze production:
1393	WRITE(*,'(a,L2)') "mm_call_hazeprod : ", mm_call_hazeprod
1394	WRITE(*,'(a,ES14.7)') "mm_rc_prod (m) : ", mm_rc_prod
1395	WRITE(*,'(a,L2)') "mm_call_CH4hazeprod : ", mm_call_CH4hazeprod
1396	IF (.NOT. mm_call_CH4hazeprod) THEN
1397	WRITE(*,'(a,ES14.7)') " --> mm_p_prod (Pa) : ", mm_p_prod
1398	WRITE(*,'(a,ES14.7)') " --> mm_tx_prod (kg.m-2.s-1): ", mm_tx_prod
1399
1400	ENDIF
1401
1402	! Haze coagulation:
1403	WRITE(*,'(a,L2)') "mm_call_hazecoag : ", mm_call_hazecoag
1404	IF (mm_call_hazecoag) THEN
1405	WRITE(*,'(a,I2.2)') " --> mm_coag_interactions : ", mm_coag_choice
1406	ENDIF
1407
1408	! Haze sedimentation:
1409	WRITE(*,'(a,L2)') "mm_call_hazesed : ", mm_call_hazesed
1410	IF (mm_call_hazesed) THEN
1411	WRITE(*,'(a,L2)') " --> mm_wsed_m0 : ", mm_wsed_m0
1412	WRITE(*,'(a,L2)') " --> mm_wsed_m3 : ", mm_wsed_m3
1413	ENDIF
1414	WRITE(*,'(a)') "---------------------------------------"
1415
1416	! Haze threshold:
1417	WRITE(*,'(a)') "Spherical aerosol thresholds"
1418	WRITE(*,'(a,ES14.7)') " mm_m0as_min (m-3) : ", mm_m0as_min
1419	WRITE(*,'(a,ES14.7)') " mm_rcs_min (m) : ", mm_rcs_min
1420	WRITE(*,'(a)') "Fractal aerosol thresholds"
1421	WRITE(*,'(a,ES14.7)') " mm_m0af_min (m-3) : ", mm_m0af_min
1422	WRITE(*,'(a,ES14.7)') " mm_rcf_min (m) : ", mm_rcf_min
1423	WRITE(*,'(a)') "---------------------------------------"
1424
1425	! Clouds related:
1426	WRITE(*,'(a,L2)') "mm_call_clouds : ", mm_call_clouds
1427	IF (mm_call_clouds) THEN
1428	WRITE(*,'(a)') " Thresholds clouds drop"
1429	WRITE(*,'(a,ES14.7)') " --> mm_m0ccn_min : ", mm_m0ccn_min
1430	WRITE(*,'(a,ES14.7)') " --> mm_drad_min : ", mm_drad_min
1431	WRITE(*,'(a,ES14.7)') " --> mm_drad_max : ", mm_drad_max
1432	ENDIF
1433	WRITE(*,'(a)') "---------------------------------------"
1434
1435	! Free parameters:
1436	WRITE(*,'(a)') "Free parameters"
1437	WRITE(*,'(a,ES14.7)') "mm_rm (m) : ", mm_rm
1438	WRITE(*,'(a,ES14.7)') "mm_df (-) : ", mm_df
1439	WRITE(*,'(a,ES14.7)') "mm_rhoaer (kg.m-3) : ", mm_rhoaer
1440	WRITE(*,'(a,ES14.7)') "mm_rplanete (m) : ", mm_rpla
1441	WRITE(*,'(a,ES14.7)') "mm_g0 (m.s-2) : ", mm_g0
1442	WRITE(*,'(a,ES14.7)') "mm_air_rad (m) : ", mm_air_rad
1443	WRITE(*,'(a,ES14.7)') "mm_air_mmol (kg.mol-1) : ", mm_air_mmol
1444	WRITE(*,'(a,ES14.7)') "mm_air_mmas (kg) : ", mm_air_mmas
1445	WRITE(*,'(a,ES14.7)') "mm_dt (s) : ", mm_dt
1446	IF (mm_nla > -1) THEN
1447	WRITE(*,'(a,I3.3)') "mm_nla : ", mm_nla
1448	ELSE
1449	WRITE(*,'(a)') "mm_nla : not initialized yet"
1450	ENDIF
1451	WRITE(*,'(a)') "======================================="
1452	END SUBROUTINE mm_dump_parameters
1453
1454
1455	! =========================================================================
1456	! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1457	! CONFIGURATION PARSER checking methods
1458	! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1459	! =========================================================================
1460
1461	FUNCTION check_r1(err,var,def,wlog) RESULT(ret)
1462	!! Check an option value (float).
1463	!!
1464	!! The method checks an option value and optionally set a default value, __def__ to initialize
1465	!! __var__ on error if given.
1466	!!
1467	TYPE(error), INTENT(in) :: err ! Error object from value getter.
1468	REAL(kind=mm_wp), INTENT(inout) :: var ! Input/output option value.
1469	REAL(kind=mm_wp), INTENT(in), OPTIONAL :: def ! Default value to set.
1470	LOGICAL, INTENT(in), OPTIONAL :: wlog ! .true. to print warning/error message.
1471	TYPE(error) :: ret ! Input error.
1472	CHARACTER(len=*), PARAMETER :: defmsg = '... Using default value: '
1473	LOGICAL :: zlog
1474	ret = err
1475	zlog = .false. ; IF (PRESENT(wlog)) zlog = wlog
1476	IF (err == 0) RETURN
1477	IF (PRESENT(def)) THEN
1478	var = def
1479	IF (zlog) WRITE(*,'(a,a,a)') error_to_string(err,'',.true.),defmsg,to_string(var)
1480	ret = noerror
1481	ELSE
1482	IF (zlog) WRITE(*,'(a)') error_to_string(err,'',.true.)
1483	ENDIF
1484	END FUNCTION check_r1
1485
1486	FUNCTION check_l1(err,var,def,wlog) RESULT(ret)
1487	!! Check an option value (logical).
1488	!!
1489	!! The method checks an option value and optionally set a default value, __def__ to initialize
1490	!! __var__ on error if given.
1491	!!
1492	TYPE(error), INTENT(in) :: err ! Error object from value getter.
1493	LOGICAL, INTENT(inout) :: var ! Input/output option value.
1494	LOGICAL, INTENT(in), OPTIONAL :: def ! Default value to set.
1495	LOGICAL, INTENT(in), OPTIONAL :: wlog ! .true. to print warning/error message.
1496	TYPE(error) :: ret ! Input error.
1497	CHARACTER(len=*), PARAMETER :: defmsg = '... Using default value: '
1498	LOGICAL :: zlog
1499	ret = err
1500	zlog = .false. ; IF (PRESENT(wlog)) zlog = wlog
1501	IF (err == 0) RETURN
1502	IF (PRESENT(def)) THEN
1503	var = def
1504	IF (zlog) WRITE(*,'(a,a,a)') error_to_string(err,'',.true.),defmsg,to_string(var)
1505	ret = noerror
1506	ELSE
1507	IF (zlog) WRITE(*,'(a)') error_to_string(err,'',.true.)
1508	ENDIF
1509	END FUNCTION check_l1
1510
1511	FUNCTION check_i1(err,var,def,wlog) RESULT(ret)
1512	!! Check an option value (integer).
1513	!!
1514	!! The method checks an option value and optionally set a default value, __def__ to initialize
1515	!! __var__ on error if given.
1516	!!
1517	TYPE(error), INTENT(in) :: err ! Error object from value getter.
1518	INTEGER, INTENT(inout) :: var ! Input/output option value.
1519	INTEGER, INTENT(in), OPTIONAL :: def ! Default value to set.
1520	LOGICAL, INTENT(in), OPTIONAL :: wlog ! .true. to print warning/error message.
1521	TYPE(error) :: ret ! Input error.
1522	CHARACTER(len=*), PARAMETER :: defmsg = '... Using default value: '
1523	LOGICAL :: zlog
1524	ret = err
1525	zlog = .false. ; IF (PRESENT(wlog)) zlog = wlog
1526	IF (err == 0) RETURN
1527	IF (PRESENT(def)) THEN
1528	var = def
1529	IF (zlog) WRITE(*,'(a,a,a)') error_to_string(err,'',.true.),defmsg,to_string(var)
1530	ret = noerror
1531	ELSE
1532	IF (zlog) WRITE(*,'(a)') error_to_string(err,'',.true.)
1533	ENDIF
1534	END FUNCTION check_i1
1535
1536	FUNCTION check_s1(err,var,def,wlog) RESULT(ret)
1537	!! Check an option value (string).
1538	!!
1539	!! The method checks an option value and optionally set a default value, __def__ to initialize
1540	!! __var__ on error if given.
1541	!!
1542	TYPE(error), INTENT(in) :: err ! Error object from value getter.
1543	CHARACTER(len=*), INTENT(inout) :: var ! Input/output option value.
1544	CHARACTER(len=*), INTENT(in), OPTIONAL :: def ! Default value to set.
1545	LOGICAL, INTENT(in), OPTIONAL :: wlog ! .true. to print warning/error message.
1546	TYPE(error) :: ret ! Input error.
1547	CHARACTER(len=*), PARAMETER :: defmsg = '... Using default value: '
1548	LOGICAL :: zlog
1549	ret = err
1550	zlog = .false. ; IF (PRESENT(wlog)) zlog = wlog
1551	IF (err == 0) RETURN
1552	IF (PRESENT(def)) THEN
1553	var = TRIM(def)
1554	IF (zlog) WRITE(*,'(a,a,a)') error_to_string(err,'',.true.),defmsg,var
1555	ret = noerror
1556	ELSE
1557	IF (zlog) WRITE(*,'(a)') error_to_string(err,'')
1558	ENDIF
1559	RETURN
1560	END FUNCTION check_s1
1561
1562	END MODULE MP2M_GLOBALS

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

Download in other formats:

Original Format