Context Navigation

mp2m_globals.F90 @ 3986

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

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

Download in other formats:

Original Format