source: trunk/LMDZ.TITAN/libf/muphytitan/mm_globals.f90 @ 1814

Last change on this file since 1814 was 1814, checked in by jvatant, 7 years ago

Correct string management within muphy for ifort
JVO
File size: 66.6 KB
Line 
1! Copyright 2013-2015 Université de Reims Champagne-Ardenne
2! Contributor: J. Burgalat (GSMA, URCA)
3! email of the author : jeremie.burgalat@univ-reims.fr
4!
5! This software is a computer program whose purpose is to compute
6! microphysics processes using a two-moments scheme.
7!
8! This library is governed by the CeCILL-B license under French law and
9! abiding by the rules of distribution of free software.  You can  use,
10! modify and/ or redistribute the software under the terms of the CeCILL-B
11! license as circulated by CEA, CNRS and INRIA at the following URL
12! "http://www.cecill.info".
13!
14! As a counterpart to the access to the source code and  rights to copy,
15! modify and redistribute granted by the license, users are provided only
16! with a limited warranty  and the software's author,  the holder of the
17! economic rights,  and the successive licensors  have only  limited
18! liability.
19!
20! In this respect, the user's attention is drawn to the risks associated
21! with loading,  using,  modifying and/or developing or reproducing the
22! software by the user in light of its specific status of free software,
23! that may mean  that it is complicated to manipulate,  and  that  also
24! therefore means  that it is reserved for developers  and  experienced
25! professionals having in-depth computer knowledge. Users are therefore
26! encouraged to load and test the software's suitability as regards their
27! requirements in conditions enabling the security of their systems and/or
28! data to be ensured and,  more generally, to use and operate it in the
29! same conditions as regards security.
30!
31! The fact that you are presently reading this means that you have had
32! knowledge of the CeCILL-B license and that you accept its terms.
33
34!! file: mm_globals.f90
35!! summary: Parameters and global variables module.
36!! author: J. Burgalat
37!! date: 2013-2015
38
39MODULE MM_GLOBALS
40  !! Parameters and global variables module.
41  !! 
42  !! # Module overview
43  !!
44  !! The module defines all the parameters and global variables that are common
45  !! to all other modules of the library.
46  !!
47  !! It is separated in two parts :
48  !!
49  !! - Main parameters and global saved variables. Most of these variables should
50  !!   be initialized once and should hold the same value during run-time. These
51  !!   variables are completly public and initialized by [[mm_globals(module):mm_global_init(interface)]]
52  !!   method.
53  !! - The second part defines a set of vectors that defines the vertical structure of the atmosphere.
54  !!   Each time a new atmospheric column has to be computed (either on a new timestep or on a new couple
55  !!   of longitude/latitude), these vectors should be intialized with new values by calling
56  !!   [[mm_globals(module):mm_column_init(function)]] method.
57  !!   This part is separated in two sets :
58  !!
59  !!   - The atmospheric structure with temperature, pressure levels and altitude definitions.
60  !!   - The vertical profiles of tracers with the moments of the two aerosols modes (both \(M_{0}\)
61  !!     and \(M_{3}\) for a total of 4 vectors), the _clouds_ microphysics moments tracers (i.e.
62  !!     \(M_{0}\) and \(M_{3}\) for the ccn and \(M_{3}\) for the ice components).
63  !!     Additionally, the module also stores intermediates variables of interest such as the
64  !!     characteristic radii of the aerosols modes, the mean drop radius and the drop density,
65  !!     the molar fraction of each condensible species (related to ice components) and some
66  !!     scalar variables that holds arrays sizes.
67  !!
68  !! @note
69  !! All the vectors that represent the vertical structure of the atmosphere (altitude, pressure and
70  !! temperature...) are oriented from the __TOP__ of the atmosphere to the __GROUND__.
71  !!
72  !! @note
73  !! The module also imports errors module from __FCCP__ library to get definitions of the error object
74  !! everywhere in the library ([[mm_globals(module)]] is always imported, except in [[mm_mprec(module)]]).
75  !!
76  !! # Global variables
77  !!
78  !! [[mm_globals(module)]] module contains the declaration of all global/common variable that are shared
79  !! by all other modules of the model. Except for few physical constant which are declared as parameters,
80  !! these variables are onlu SAVEd. They are initialized by [[mm_globals(module):mm_global_init(interface)]]
81  !! methods.
82  !! the following sections list all the global variables by category.
83  !!
84  !! ## Control flags
85  !!
86  !! | Name               | Description
87  !! | :----------------- | :-----------------
88  !! | mm_log             | Enable log mode (verbose)
89  !! | mm_w_haze_prod     | Enable/Disable haze production
90  !! | mm_w_haze_sed      | Enable/Disable haze sedimentation
91  !! | mm_w_haze_coag     | Enable/Disable haze coagulation
92  !! | mm_w_clouds        | Enable/Disable clouds microphysics
93  !! | mm_w_clouds_sed    | Enable/Disable clouds microphysics sedimentation
94  !! | mm_w_clouds_nucond | Enable/Disable clouds microphysics nucleation/condensation
95  !! | mm_wsed_m0         | Force all aerosols moments to fall at M0 settling velocity
96  !! | mm_wsed_m3         | Force all aerosols moments to fall at M3 settling velocity
97  !! | mm_no_fiadero_w    | Enable/Disable __Fiadero__ correction
98  !!
99  !! ### Related free parameters:
100  !!
101  !! | Name            | Description
102  !! | :-------------- | :-----------------
103  !! | mm_fiadero_min  | Minimum ratio for __Fiadero__'s correction
104  !! | mm_fiadero_max  | Maximum ratio for __Fiadero__'s correction
105  !! | mm_coag_choice  | Coagulation interaction activation flag. It should be a combination of [[mm_globals(module):mm_coag_no(variable)]], [[mm_globals(module):mm_coag_ss(variable)]], [[mm_globals(module):mm_coag_sf(variable)]] and [[mm_globals(module):mm_coag_ff(variable)]].
106  !!
107  !! ## Physical constants
108  !!
109  !! | Name      | Description
110  !! | :-------- | :-----------------
111  !! | mm_pi     | Pi number
112  !! | mm_navo   | Avogadro number
113  !! | mm_kboltz | Boltzmann constant (\(J.K^{-1}\))
114  !! | mm_rgas   | Perfect gas constant (\(J.mol^{-1}.K^{-1}\))
115  !! | mm_fdes   | Desorption energy (\(J\)) (nucleation)
116  !! | mm_fdif   | Surface diffusion energy (\(J\)) (nucleation)
117  !! | mm_fnus   | Jump frequency (\(s^{-1}\)) (nucleation)
118  !! | mm_akn    | Approximated slip-flow correction coefficient (
119  !!
120  !! ## Free parameters
121  !!
122  !! | Name        | Description
123  !! | :---------- | :-----------------
124  !! | mm_rhoaer   | Aerosol density (in \(kg.m^{-3}\))
125  !! | mm_df       | Fractal dimension
126  !! | mm_rm       | Monomer radius (in m)
127  !! | mm_p_prod   | Spherical aerosols production pressure level (Pa)
128  !! | mm_p_rcprod | Spherical aerosols equivalent radius production (m)
129  !! | mm_tx_prod  | Production rate of spherical aerosols (\(kg.m^{-2}.s^{-1}\))
130  !! | mm_d_prod   | Time-dependent sine wve pre-factor.
131  !! | mm_w_prod   | Angular frequency of the time-dependent production rate.
132  !! | mm_ne       | Electric charging of aerosols (\(e^{-}.m^{-1}\)) (unused)
133  !! | mm_rb2ra    | Bulk to apparent radius conversion pre-factor (\(m^X\))
134  !! | mm_rpla     | Planet radius (m)
135  !! | mm_g0       | Planet acceleration due to gravity constant (ground) (\(m.s^{-2}\))
136  !! | mm_air_rad  | Air molecules mean radius (m)
137  !! | mm_air_mmol | Air molecules molar mass (\(kg.mol^{-1}\))
138  !! | mm_dt       | Microphysic time step (s)
139  USE MM_MPREC
140  USE MM_INTERFACES
141  ! from swift
142  USE CFGPARSE
143  USE STRINGS
144  USE ERRORS
145  IMPLICIT NONE
146
147  PUBLIC
148
149  PRIVATE :: cldprop_sc,cldprop_ve,read_esp,check_r1,check_i1,check_l1,check_s1
150
151  ! Protected variables
152  ! the following variables are read-only outside this module.
153  ! One must call the afferent subroutine to update them.
154   
155  ! initialization control flags (cannot be updated)
156  PROTECTED :: mm_ini,mm_ini_col,mm_ini_aer,mm_ini_cld
157  ! model parameters (mm_global_init)
158  PROTECTED :: mm_dt,mm_rhoaer,mm_df,mm_rm,mm_p_prod,mm_rc_prod,mm_tx_prod,mm_rpla,mm_g0,mm_rb2ra
159  ! atmospheric vertical structure (mm_column_init)
160  PROTECTED :: mm_nla,mm_nle,mm_zlay,mm_zlev,mm_play,mm_plev,mm_temp,mm_rhoair,mm_btemp,mm_dzlev,mm_dzlay
161  ! Condensible species parameters (mm_global_init)
162  PROTECTED :: mm_nesp,mm_spcname,mm_xESPS
163  ! Moments parameters (mm_aerosols_init / mm_clouds_init)
164  PROTECTED :: mm_m0aer_s, mm_m3aer_s, mm_m0aer_f, mm_m3aer_f, mm_m0ccn, mm_m3ccn, mm_m3ice
165  ! Moments parameters (derived, are updated with moments parameters)
166  PROTECTED :: mm_rcs, mm_rcf, mm_drad, mm_drho
167
168  LOGICAL, SAVE :: mm_debug = .true.  !! Enable QnD debug mode.
169  LOGICAL, SAVE :: mm_log = .false.   !! Enable log mode.
170
171  LOGICAL, SAVE :: mm_w_haze_prod = .true. !! Enable/Disable haze production.
172  LOGICAL, SAVE :: mm_w_haze_sed = .true.  !! Enable/Disable haze sedimentation.
173  LOGICAL, SAVE :: mm_w_haze_coag = .true. !! Activate haze coagulation.
174
175  LOGICAL, SAVE :: mm_wsed_m0 = .false. !! Force all aerosols moments to fall at M0 settling velocity.
176  LOGICAL, SAVE :: mm_wsed_m3 = .false. !! Force all aerosols moments to fall at M3 settling velocity.
177
178  LOGICAL, SAVE :: mm_var_prod = .false. !! Time variation of production rate control flag.
179
180  !> Enable/Disable __Fiadero__'s correction.
181  !!
182  !! This flag enables/disables the __Fiadero__ correction alogrithm for fractal mode settling velocity
183  !! computation.
184  !!
185  !! @bug
186  !! Currently, the Fiadero correction creates instatibilities on the vertical structure. It seems to be
187  !! related to the coupling between the two moments. In order to reduce the instabilities, settling
188  !! velocity of moments are forced to be the same, see [[mm_globals(module):mm_wsed_m0(variable)]] and
189  !! [[mm_globals(module):mm_wsed_m3(variable)]]).
190  LOGICAL, SAVE          :: mm_no_fiadero_w = .false.
191
192  !> Minimum ratio for __Fiadero__ correction.
193  !!
194  !! When [[mm_globals(module):mm_no_fiadero_w(variable)]] is disabled, this variable defines the minimum
195  !! value of the moment's ratio between two adjacents vertical cells to be used within the correction.
196  REAL(kind=mm_wp), SAVE :: mm_fiadero_min  = 0.1_mm_wp
197
198  !> Maximum ratio for __Fiadero__ correction.
199  !!
200  !! When [[mm_globals(module):mm_no_fiadero_w(variable)]] is disabled, this variable defines the maximum
201  !! value of the moment's ratio between two adjacents vertical cells to be used within the correction.
202  REAL(kind=mm_wp), SAVE :: mm_fiadero_max  = 10._mm_wp
203
204  LOGICAL, SAVE :: mm_w_clouds = .true.       !! Enable/Disable clouds microphysics.
205  LOGICAL, SAVE :: mm_w_cloud_sed = .true.    !! Enable/Disable cloud sedimentation.
206  LOGICAL, SAVE :: mm_w_cloud_nucond = .true. !! Activate cloud nucleation/condensation.
207
208  INTEGER, PARAMETER :: mm_coag_no = 0 !! no mode interaction for coagulation (i.e. no coagulation at all).
209  INTEGER, PARAMETER :: mm_coag_ss = 1 !! SS mode interaction for coagulation.
210  INTEGER, PARAMETER :: mm_coag_sf = 2 !! SF mode interaction for coagulation.
211  INTEGER, PARAMETER :: mm_coag_ff = 4 !! FF mode interaction for coagulation.
212  !> Default interactions to activate (all by default).
213  INTEGER, SAVE      :: mm_coag_choice = mm_coag_ss+mm_coag_sf+mm_coag_ff
214
215  !> Pi number.
216  REAL(kind=mm_wp), PARAMETER :: mm_pi = 4._mm_wp*atan(1._mm_wp)
217  !> Avogadro number.
218  REAL(kind=mm_wp), PARAMETER :: mm_navo = 6.0221367e23_mm_wp
219  !> Boltzmann constant (\(J.K^{-1}\)).
220  REAL(kind=mm_wp), PARAMETER :: mm_kboltz = 1.3806488e-23_mm_wp
221  !> Perfect gas constant (\(J.mol^{-1}.K^{-1}\)).
222  REAL(kind=mm_wp), PARAMETER :: mm_rgas = mm_kboltz * mm_navo
223  !> Desorption energy (\(J\)) (nucleation).
224  REAL(kind=mm_wp), PARAMETER :: mm_fdes = 0.288e-19_mm_wp
225  !> Surface diffusion energy (\(J\)) (nucleation).
226  REAL(kind=mm_wp), PARAMETER :: mm_fdif = 0.288e-20_mm_wp
227  !> Jump frequency (\(s^{-1}\)) (nucleation).
228  REAL(kind=mm_wp), PARAMETER :: mm_nus = 1.e+13_mm_wp
229  !> Approximated slip-flow correction coefficient.
230  REAL(kind=mm_wp), PARAMETER :: mm_akn = 1.591_mm_wp
231
232  !> Aerosols density (\(kg.m^{-3}\)).
233  REAL(kind=mm_wp), SAVE :: mm_rhoaer = 1.e3_mm_wp
234
235  !> Fractal dimension of fractal aerosols.
236  REAL(kind=mm_wp), SAVE :: mm_df = 3._mm_wp
237
238  !> Monomer radius (m).
239  REAL(kind=mm_wp), SAVE :: mm_rm = 6.66e-8_mm_wp
240
241  !> Spherical aerosols production pressure level (Pa).
242  REAL(kind=mm_wp), SAVE :: mm_p_prod = 1._mm_wp
243
244  !> Spherical aerosols equivalent radius production (m)
245  REAL(kind=mm_wp), SAVE :: mm_rc_prod = 1.3101721857598102e-9_mm_wp
246
247  !> Production rate of spherical aerosols (\(kg.m^{-2}.s^{-1}\)).
248  REAL(kind=mm_wp), SAVE :: mm_tx_prod = 3.5e-13_mm_wp
249
250  !> Aerosol production delta if time variations is enabled (fraction).
251  REAL(kind=mm_wp), SAVE :: mm_d_prod  = 0.25_mm_wp
252
253  !> Aerosol production variations angular frequency if time variations is enabled (\(rad.s^{-1}\)).
254  REAL(kind=mm_wp), SAVE :: mm_w_prod  = 2.*mm_pi / (86400.*16.)
255
256
257  !> Electric charging of aerosols (\(e^{-}.m^{-1}\)).
258  REAL(kind=mm_wp), SAVE :: mm_ne = -15.e6_mm_wp
259
260  !> Bulk to apparent radius conversion pre-factor (\(m^{X}\)).
261  !!
262  !! It is initialized using [[mm_globals(module):mm_rm(variable)]] in
263  !! [[mm_globals(module):mm_global_init(interface)]] from the following equation:
264  !!
265  !! $$ r_{a} = r_{b}^{3/D_{f}}\times r_{m}^{\frac{D_{f}-3}{D_{f}}} $$
266  !!
267  !! Where \(r_{a}\) is the apparent radius, \(r_{b}\) the bulk radius and
268  !! \(rb2ra = r_{m}^{\frac{D_{f}-3}{D_{f}}}\) is the returned pre-factor
269  REAL(kind=mm_wp), SAVE :: mm_rb2ra = 1._mm_wp
270
271  !> Characteristic radius threshold.
272  REAL(kind=mm_wp), SAVE :: mm_rc_min = 1.e-200_mm_wp
273
274  !> Name of condensible species.
275  CHARACTER(len=30), DIMENSION(:), ALLOCATABLE, SAVE :: mm_spcname
276
277  TYPE, PUBLIC :: mm_esp
278    !! Cloud related chemical specie properties.
279    !!
280    !! This derived type is used in thermodynamic methods related to cloud microphysics.
281    !! Most of its fields represent parameters of equations from \cite{reid1986}.
282    CHARACTER(LEN=10) :: name      !! Specie name.
283    REAL(kind=mm_wp)  :: mas       !! Molecular weight (kg).
284    REAL(kind=mm_wp)  :: vol       !! Molecular volume (\(m^{3}\)).
285    REAL(kind=mm_wp)  :: ray       !! Molecular radius (m).
286    REAL(kind=mm_wp)  :: masmol    !! Molar mass (\(kg.mol^{-1}\)).
287    REAL(kind=mm_wp)  :: rho       !! density (liquid) (\(kg.m^{-3}\)).
288    REAL(kind=mm_wp)  :: tc        !! Critical temperature (K).
289    REAL(kind=mm_wp)  :: pc        !! Critical pressure (Bar).
290    REAL(kind=mm_wp)  :: tb        !! Boiling point temperature (K).
291    REAL(kind=mm_wp)  :: w         !! Acentric factor (--).
292    REAL(kind=mm_wp)  :: a_sat     !! Saturation equation A coefficient.
293    REAL(kind=mm_wp)  :: b_sat     !! Saturation equation B coefficient.
294    REAL(kind=mm_wp)  :: c_sat     !! saturation equation C coefficient.
295    REAL(kind=mm_wp)  :: d_sat     !! Saturation equation D coefficient.
296    REAL(kind=mm_wp)  :: mteta     !! Wettability.
297    REAL(kind=mm_wp)  :: tx_prod   !! Production rate.
298    REAL(kind=mm_wp)  :: fmol2fmas !! molar fraction to mass fraction coefficient.
299    ! = masmol(X)/masmol(AIR)
300  END TYPE mm_esp
301
302  !> Planet radius (m).
303  REAL(kind=mm_wp), SAVE                        :: mm_rpla     = 2575000._mm_wp
304  !> Planet acceleration due to gravity constant (ground) (\(m.s^{-2}\)).
305  REAL(kind=mm_wp), SAVE                        :: mm_g0       = 1.35_mm_wp
306  !> Air molecules mean radius (m).
307  REAL(kind=mm_wp), SAVE                        :: mm_air_rad  = 1.75e-10_mm_wp
308  !> Air molecules molar mass (\(kg.mol^{-1}\)).
309  REAL(kind=mm_wp), SAVE                        :: mm_air_mmol = 28e-3_mm_wp
310  !> Microphysic time step (s).
311  REAL(kind=mm_wp), SAVE                        :: mm_dt       = 5529.6_mm_wp
312  !> Model current time tracer (s).
313  REAL(kind=mm_wp), SAVE                        :: mm_ct       = 0.0
314  !> Total number of clouds condensible species.
315  INTEGER, SAVE                                 :: mm_nesp     = -1
316  !> Clouds chemical species properties.
317  TYPE(mm_esp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_xESPS
318
319  !------------------------
320  ! Vertical structure part
321  !------------------------
322
323  !> Number of vertical layers.
324  INTEGER, SAVE :: mm_nla = -1
325  !> Number of vertical levels.
326  INTEGER, SAVE :: mm_nle = -1
327
328  !> Altitude layers (m).
329  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_zlay
330  !> Altitude levels (m).
331  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_zlev
332  !> Pressure layers (Pa).
333  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_play
334  !> Pressure levels (Pa).
335  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_plev
336  !>  Temperature vertical profile (K).
337  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_temp
338  !>  Air density vertical profile (\(kg.m^{-3}\)).
339  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_rhoair
340  !> Temperature vertical profil at interfaces (K).
341  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_btemp
342
343  !> Atmospheric levels thickness (m).
344  !!
345  !! Atmospheric thickness between two adjacent levels (\(m\)) from the
346  !! __TOP__ to the __GROUND__.
347  !! @note __mm_dzlev__ is defined on the total number of layers and actually
348  !! corresponds to the thickness of a given layer.
349  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_dzlev
350
351  !> Atmospheric layers "thickness" (m).
352  !!
353  !! Atmospheric thickness between the center of two adjacent layers (\(m\))
354  !! from the __TOP__ to the __GROUND__.
355  !! @note
356  !! __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  !! @note This value corresponds to the thickness between the center of the
359  !! __GROUND__ layer and below the surface. It is arbitrary and not used.
360  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_dzlay
361
362  !> Spherical mode \(0^{th}\) order moment (\(m^{-3}\)).
363  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m0aer_s
364  !> Spherical mode \(3^{rd}\) order moment (\(m^{3}.m^{-3}\)).
365  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m3aer_s
366  !> Fractal mode \(0^{th}\) order moment (\(m^{-3}\)).
367  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m0aer_f
368  !> Fractal mode \(3^{rd}\) order moment (\(m^{3}.m^{-3}\)).
369  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m3aer_f
370  !> CCN \(0^{th}\) order moment (\(m^{-3}\)).
371  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m0ccn
372  !> CCN \(3^{rd}\) order moment (\(m^{3}.m^{-3}\)).
373  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m3ccn
374
375  !> Ice components 3rd order moments (\(m^{3}.m^{-3}\)).
376  !!
377  !! It is a 2D array with the vertical layers in first dimension, and the number of ice
378  !! components in the second.
379  !! @note
380  !! Both [[mm_globals(module):mm_m3ice(variable)]] and [[mm_globals(module):mm_gazs(variable)]]
381  !! share the same indexing (related to species order).
382  REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: mm_m3ice
383
384  !> Condensible species molar fraction (\(mol.mol^{-1}\)).
385  !!
386  !! It is a 2D array with the vertical layers in first dimension, and
387  !! the number of condensible species in the second.
388  !! @note
389  !! Both [[mm_globals(module):mm_m3ice(variable)]] and [[mm_globals(module):mm_gazs(variable)]]
390  !! share the same indexing (related to species order).
391  REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: mm_gazs
392
393  !> Spherical mode characteristic radius (m).
394  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_rcs
395  !> Fractal mode characteristic radius (m).
396  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_rcf
397  !> Mean Drop radius (m).
398  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_drad
399  !> Mean Drop density (\(kg.m^{-3}\)).
400  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_drho
401
402  !> Aerosols precipitations (m).
403  !!
404  !! Aerosols precipitations take into account both spherical and fractal modes.
405  !! It is updated in [[mm_haze(module):mm_haze_microphysics(subroutine)]].
406  REAL(kind=mm_wp), SAVE :: mm_aer_prec = 0._mm_wp
407
408  !> Spherical mode \(M_{0}\) settling velocity (\(m.s^{-1}\)).
409  !!
410  !! It is a vector with the vertical layers that contains the settling velocity for
411  !! the \(0^{th}\) order moment of the spherical mode.
412  !! It is updated in [[mm_haze(module):mm_haze_sedimentation(subroutine)]].
413  !! @note
414  !! This variable is always negative.
415  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m0as_vsed
416
417  !> Spherical mode \(M_{3}\) settling velocity (\(m.s^{-1}\)).
418  !!
419  !! It is a vector with the vertical layers that contains the settling velocity for the
420  !! \(3^{rd}\) order moment of the spherical mode.
421  !! It is updated in [[mm_haze(module):mm_haze_sedimentation(subroutine)]].
422  !! @note
423  !! This variable is always negative.
424  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m3as_vsed
425
426  !> Fractal mode \(M_{0}\) settling velocity (\(m.s^{-1}\)).
427  !!
428  !! It is a vector with the vertical layers that contains the settling velocity for the
429  !! \(0^{th}\) order moment of the fractal mode.
430  !! It is updated in [[mm_haze(module):mm_haze_sedimentation(subroutine)]].
431  !! @note
432  !! This variable is always negative.
433  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m0af_vsed
434
435  !> Fractal mode \(M_{3}\) settling velocity (\(m.s^{-1}\)).
436  !!
437  !! It is a vector with the vertical layers that contains the settling velocity for the
438  !! \(3^{rd}\) order moment of the fractal mode.
439  !! It is updated in [[mm_haze(module):mm_haze_sedimentation(subroutine)]].
440  !! @note
441  !! This variable is always negative.
442  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m3af_vsed
443
444  !> Spherical aerosol mass fluxes (\(kg.m^{-2}.s^{-1}\)).
445  !!
446  !! It is a vector with the vertical layers that contains the mass fluxes for spherical aerosols.
447  !! It is updated in [[mm_haze(module):mm_haze_sedimentation(subroutine)]].
448  !! @note
449  !! This variable is always negative.
450  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_aer_s_flux
451
452  !> Fractal aerosol mass fluxes (\(kg.m^{-2}.s^{-1}\)).
453  !!
454  !! It is a vector with the vertical layers that contains the mass fluxes for fractal aerosols
455  !! It is updated in [[mm_haze(module):mm_haze_sedimentation(subroutine)]].
456  !! @note
457  !! This variable is always negative.
458  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_aer_f_flux
459
460  !> CCN precipitations (m).
461  !! It is updated in [[mm_clouds(module):mm_cloud_microphysics(subroutine)]].
462  REAL(kind=mm_wp), SAVE :: mm_ccn_prec = 0._mm_wp
463
464  !> CCN mass fluxes (\(kg.m^{-2}.s^{-1}\)).
465  !!
466  !! It is a vector with the vertical layers that contains the
467  !! mass fluxes for CCN.
468  !! It is updated in [[mm_clouds(module):mm_cloud_microphysics(subroutine)]].
469  !! @note
470  !! This variable is always negative.
471  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_ccn_flux
472
473  !> Ice components precipitations (m).
474  !!
475  !! It is a vector of [[mm_globals(module):mm_nesp(variable)]] values which share the same indexing
476  !! than [[mm_globals(module):mm_m3ice(variable)]] and [[mm_globals(module):mm_gazs(variable)]].
477  !! It is updated in [[mm_clouds(module):mm_cloud_microphysics(subroutine)]].
478  !! @note
479  !! This variable is always negative.
480  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_ice_prec
481
482  !> Ice components sedimentation fluxes (\(kg.m^{-2}.s-1\)).
483  !!
484  !! It is a 2D-array with the vertical layers in first dimension and the number of ice components
485  !! in the second. It is updated in [[mm_clouds(module):mm_cloud_microphysics(subroutine)]].
486  !! @note
487  !! This variable is always negative.
488  REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: mm_ice_fluxes
489
490  !> Condensible species saturation ratio (--).
491  !!
492  !! It is a 2D-array with the vertical layers in first dimension and the number of condensible
493  !! species in the second.
494  !! It is updated in [[mm_clouds(module):mm_cloud_microphysics(subroutine)]].
495  REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: mm_gazs_sat
496
497  !> [[mm_globals(module):mm_global_init(interface)]] initialization control flag.
498  LOGICAL, PUBLIC, SAVE :: mm_ini     = .false.
499
500  !> [[mm_globals(module):mm_column_init(function)]] initialization control flag.
501  LOGICAL, PUBLIC, SAVE :: mm_ini_col = .false.
502
503  !> [[mm_globals(module):mm_aerosols_init(function)]] initialization control flag.
504  LOGICAL, PUBLIC, SAVE :: mm_ini_aer = .false.
505
506  !> [[mm_globals(module):mm_clouds_init(function)]] initialization control flag.
507  LOGICAL, PUBLIC, SAVE :: mm_ini_cld = .false.
508
509  !> Interface to cloud properties methods.
510  !!
511  !! The method computes clouds properties (mean drop radius and denstity) from their afferent
512  !! moments. It is overloaded to compute properties at a single level or over all the vertical
513  !! atmospheric structure.
514  INTERFACE mm_cloud_properties
515    MODULE PROCEDURE cldprop_sc,cldprop_ve
516  END INTERFACE
517
518  !> Interface to global initialization.
519  !!
520  !! The method performs the global initialization of the model.
521  !! @warning
522  !! If OpenMP is activated, this subroutine must be called in an $OMP SINGLE statement as it
523  !! initializes global variable that are not thread private.
524  !!
525  !! '''
526  !! !$OMP SINGLE
527  !! call mm_global_init(...)
528  !! !$OMP END SINGLE
529  INTERFACE mm_global_init
530    MODULE PROCEDURE mm_global_init_0,mm_global_init_1
531  END INTERFACE
532
533  !> Check an option from the configuration system.
534  !!
535  !! The method checks for an option in the configuration system and optionally
536  !! set a default value if the option is not found. This is an overloaded method
537  !! that can take in input either a floating point, integer, logical or string
538  !! option value.
539  INTERFACE mm_check_opt
540    MODULE PROCEDURE check_r1,check_i1,check_l1,check_s1
541  END INTERFACE
542
543  ! --- OPENMP ---------------
544  ! All variable related to column computations should be private to each thread
545  !
546  !$OMP THREADPRIVATE(mm_ini_col,mm_ini_aer,mm_ini_cld)
547  !$OMP THREADPRIVATE(mm_zlay,mm_zlev,mm_play,mm_plev,mm_temp,mm_rhoair,mm_btemp,mm_dzlev,mm_dzlay)
548  !$OMP THREADPRIVATE(mm_m0aer_s,mm_m3aer_s,mm_m0aer_f,mm_m3aer_f)
549  !$OMP THREADPRIVATE(mm_m0ccn,mm_m3ccn,mm_m3ice,mm_gazs)
550  !$OMP THREADPRIVATE(mm_rcs,mm_rcf,mm_drad,mm_drho)
551  !$OMP THREADPRIVATE(mm_aer_s_flux,mm_aer_f_flux,mm_ccn_flux,mm_ice_prec,mm_ice_fluxes,mm_gazs_sat)
552  !$OMP THREADPRIVATE(mm_m0as_vsed,mm_m3as_vsed,mm_m0af_vsed,mm_m3af_vsed)
553
554  !$OMP THREADPRIVATE(mm_nla,mm_nle)
555
556  ! --------------------------
557
558
559  CONTAINS
560
561  FUNCTION mm_global_init_0(dt,df,rm,rho_aer,p_prod,tx_prod,rc_prod,rplanet,g0, &
562                            air_rad,air_mmol,coag_interactions,clouds,spcfile,  &
563                            w_haze_prod,w_haze_sed,w_haze_coag,w_cloud_nucond,  &
564                            w_cloud_sed,force_wsed_to_m0,force_wsed_to_m3,      &
565                            no_fiadero,fiadero_min,fiadero_max) RESULT(err)
566    !! Initialize global parameters of the model.
567    !!
568    !! The function initializes all the global parameters of the model from direct input.
569    !! Boolean (and Fiadero) parameters are optional as they are rather testing parameters. Their
570    !! default values are suitable for production runs. 
571    !! @note
572    !! If the method fails to initialize parameters (i.e. returned error is not 0). Then the model
573    !! should probably be aborted as the global variables of the model will not be correctly setup.
574    !! @warning
575    !! If OpenMP is activated, this subroutine must be called in an $OMP SINGLE statement as it
576    !! initializes global variable that are not thread private.
577    !!
578    !! '''
579    !! !$OMP SINGLE
580    !! call mm_global_init_0(...)
581    !! !$OMP END SINGLE
582    REAL(kind=mm_wp), INTENT(in)           :: dt
583      !! Microphysics timestep in seconds.
584    REAL(kind=mm_wp), INTENT(in)           :: df
585      !! Fractal dimension of fractal aerosol.
586    REAL(kind=mm_wp), INTENT(in)           :: rm
587      !! Monomer radius in meter.
588    REAL(kind=mm_wp), INTENT(in)           :: rho_aer
589      !! Aerosol density in \(kg.m^{-3}\).
590    REAL(kind=mm_wp), INTENT(in)           :: p_prod
591      !!  Aerosol production pressure level in Pa.
592    REAL(kind=mm_wp), INTENT(in)           :: tx_prod
593      !! Spherical aerosol mode production rate in \(kg.m^{-2}.s^{-1}\).
594    REAL(kind=mm_wp), INTENT(in)           :: rc_prod
595      !! Spherical mode characteristic radius for production in meter.
596    REAL(kind=mm_wp), INTENT(in)           :: rplanet
597      !! Planet radius in meter
598    REAL(kind=mm_wp), INTENT(in)           :: g0
599      !! Planet gravity acceleration at ground level in \(m.s^{-2}\).
600    REAL(kind=mm_wp), INTENT(in)           :: air_rad
601      !! Air molecules mean radius in meter.
602    REAL(kind=mm_wp), INTENT(in)           :: air_mmol
603      !! Air molecules mean molar mass in \(kg.mol^{-1}\).
604    INTEGER, INTENT(in)                    :: coag_interactions
605      !! Coagulation interactions process control flag.
606    LOGICAL, INTENT(in)                    :: clouds
607      !! Clouds microphysics control flag.
608    CHARACTER(len=*), INTENT(in)           :: spcfile
609      !! Clouds microphysics condensible species properties file.
610    REAL(kind=mm_wp), INTENT(in), OPTIONAL :: fiadero_max
611      !! Maximum moment ratio threshold for Fiadero correction (default: 10.) .
612    REAL(kind=mm_wp), INTENT(in), OPTIONAL :: fiadero_min
613      !! Minimum moment ratio threshold for Fiadero correction (default: 0.1).
614    LOGICAL, INTENT(in), OPTIONAL          :: w_haze_prod
615      !! Haze microphysics production process control flag (default: T).
616    LOGICAL, INTENT(in), OPTIONAL          :: w_haze_sed
617      !! Haze microphysics sedimentation process control flag (default: T).
618    LOGICAL, INTENT(in), OPTIONAL          :: w_haze_coag
619      !! Haze microphysics coagulation process control flag (default: T).
620    LOGICAL, INTENT(in), OPTIONAL          :: w_cloud_sed
621      !! Cloud microphysics nucleation/conensation process control flag (default: __clouds__ value).
622    LOGICAL, INTENT(in), OPTIONAL          :: w_cloud_nucond
623      !! Cloud microphysics production process control flag (default: __clouds__ value).
624    LOGICAL, INTENT(in), OPTIONAL          :: no_fiadero
625      !! Disable Fiadero correction for haze sedimentation process (default: F).
626    LOGICAL, INTENT(in), OPTIONAL          :: force_wsed_to_m0
627      !! force __all__ aerosols moments to fall at M0 settling velocity (default: T).
628    LOGICAL, INTENT(in), OPTIONAL          :: force_wsed_to_m3
629      !! Force __all__ aerosols moments to fall at M3 settling velocity (default: F).
630    TYPE(error) :: err
631      !! Error status of the function.
632    INTEGER                                           :: i
633    TYPE(cfgparser)                                   :: cp
634    CHARACTER(len=st_slen)                            :: spcpath
635    CHARACTER(len=:), ALLOCATABLE                     :: defmsg
636    CHARACTER(len=st_slen), DIMENSION(:), ALLOCATABLE :: species
637    REAL(kind=mm_wp)                                  :: zfiamin,zfiamax
638    LOGICAL                                           :: zwhp,zwhs,zwhc,zwcs,zwcn,znofia, &
639                                                         zwstom0,zwstom3
640
641    zwhp = .true. ; zwhs = .true. ; zwhc = .true.
642    zwcs = clouds ; zwcn = clouds
643    znofia = .false. ; zfiamin = 0.1_mm_wp ; zfiamax = 10._mm_wp
644    zwstom0 = .true. ; zwstom3 = .false.
645    err = noerror
646    IF (mm_ini) THEN
647      err = error("mm_global_init: YAMMS global initialization already performed !",-1)
648      RETURN
649    ENDIF
650
651    ! Store options values in global variables...
652    mm_df          = df
653    mm_rm          = rm
654    mm_rb2ra       = mm_rm**((mm_df-3._mm_wp)/mm_df) ! conversion factor for bulk -> fractal radius
655    mm_rhoaer      = rho_aer
656    mm_p_prod      = p_prod
657    mm_tx_prod     = tx_prod
658    mm_rc_prod     = rc_prod
659    mm_rpla        = rplanet
660    mm_g0          = g0
661    mm_dt          = dt
662    mm_air_rad     = mm_air_rad
663    mm_air_mmol    = air_mmol
664    mm_coag_choice = coag_interactions
665    ! check coagulation interactions choice
666    IF (mm_coag_choice < 0 .OR. mm_coag_choice > 7) THEN
667      err = error("mm_global_init: Invalid choice for coagulation interactions activation",-1)
668      RETURN
669    ENDIF
670
671    mm_w_clouds = clouds
672
673    ! Check clouds microphysics species file
674    ! (only if clouds is activated)
675    IF (mm_w_clouds) THEN
676      IF (LEN_TRIM(spcfile) == 0) THEN
677        err = error("mm_global_init: No species properties file given",-1)
678        RETURN
679      ENDIF
680      ! Reads species properties configuration file 
681      err = cfg_read_config(cp,TRIM(spcfile)) ; IF (err /= 0) RETURN
682      err = cfg_get_value(cp,"used_species",species)
683      IF (err /= 0) THEN
684        err = error("mm_global_init: cannot retrieve 'used_species' values",-1)
685        RETURN
686      ENDIF
687      ! Now attempts to find species properties !!!
688      mm_nesp = SIZE(species)
689      ALLOCATE(mm_spcname(mm_nesp),mm_xESPS(mm_nesp))
690      DO i=1,mm_nesp
691        mm_spcname(i) = str_to_lower(species(i))
692        IF(.NOT.cfg_has_section(cp,TRIM(mm_spcname(i)))) THEN
693          err = error("mm_global_init: Cannot find "//TRIM(mm_spcname(i))//" properties",-1)
694          RETURN
695        ELSE
696          err = read_esp(cp,TRIM(mm_spcname(i)),mm_xESPS(i))
697          ! compute conversion factor: mol.mol-1 => kg.kg-1
698          mm_xESPS(i)%fmol2fmas = mm_xESPS(i)%masmol / mm_air_mmol
699          IF (err/=0) THEN
700            err = error("mm_global_init: "//TRIM(mm_spcname(i))//": "//TRIM(err%msg),-1)
701            RETURN
702          ENDIF
703        ENDIF
704      ENDDO
705    ENDIF
706
707    ! optional flags
708    ! haze control flags
709    IF (PRESENT(w_haze_prod)) THEN
710      mm_w_haze_prod = w_haze_prod
711    ELSE
712      mm_w_haze_prod = zwhp
713      call printw("mm_haze_production",to_string(mm_w_haze_prod))
714    ENDIF
715    IF (PRESENT(w_haze_sed)) THEN
716      mm_w_haze_sed = w_haze_sed
717    ELSE
718      mm_w_haze_sed = zwhs
719      call printw("mm_haze_sedimentation",to_string(mm_w_haze_sed))
720    ENDIF
721    IF (PRESENT(w_haze_coag)) THEN
722      mm_w_haze_coag = w_haze_coag
723    ELSE
724      mm_w_haze_coag = zwhc
725      call printw("mm_haze_coagulation",to_string(mm_w_haze_coag))
726    ENDIF
727    IF (PRESENT(force_wsed_to_m0)) THEN
728      mm_wsed_m0 = force_wsed_to_m0
729    ELSE
730      mm_wsed_m0 = zwstom0
731      call printw("mm_wsed_m0",to_string(mm_wsed_m0))
732    ENDIF
733    IF (PRESENT(force_wsed_to_m3)) THEN
734      mm_wsed_m3 = force_wsed_to_m3
735    ELSE
736      mm_wsed_m3 = zwstom3
737      call printw("mm_wsed_m3",to_string(mm_wsed_m3))
738    ENDIF
739    IF (PRESENT(no_fiadero)) THEN
740      mm_no_fiadero_w = no_fiadero
741    ELSE
742      mm_no_fiadero_w = znofia
743      call printw("mm_no_fiadero",to_string(mm_no_fiadero_w))
744    ENDIF
745    IF (PRESENT(fiadero_min)) THEN
746      mm_fiadero_min = fiadero_min
747    ELSE
748      mm_fiadero_min = zfiamin
749      call printw("mm_fiadero_min",to_string(mm_fiadero_min))
750    ENDIF
751    IF (PRESENT(fiadero_max)) THEN
752      mm_fiadero_max = fiadero_max
753    ELSE
754      mm_fiadero_max = zfiamax
755      call printw("mm_fiadero_max",to_string(mm_fiadero_max))
756    ENDIF
757    ! clouds control flags
758    IF (mm_w_clouds) THEN
759      IF (PRESENT(w_cloud_sed)) THEN
760        mm_w_cloud_sed = w_cloud_sed
761      ELSE
762        mm_w_cloud_sed = zwcs
763        call printw("mm_cloud_sed",to_string(mm_w_cloud_sed))
764      ENDIF
765      IF (PRESENT(w_cloud_nucond)) THEN
766        mm_w_cloud_nucond = w_cloud_nucond
767      ELSE
768        mm_w_cloud_nucond = zwcs
769        call printw("mm_cloud_nucond",to_string(mm_w_cloud_nucond))
770      ENDIF
771    ENDIF
772
773    ! check w sed flags
774    err = noerror
775    ! special check for settling velocity
776    IF (mm_wsed_m0 .AND. mm_wsed_m3) THEN
777      err = error("'wsed_m0' and 'wsed_m3' options are mutually exclusive",-1)
778    ENDIF
779    mm_ini = err == noerror
780
781    CONTAINS
782
783    SUBROUTINE printw(string,value)
784      !! Print a warning message.
785      CHARACTER(len=*), INTENT(in) :: string !! Name of the option.
786      CHARACTER(len=*), INTENT(in) :: value  !! (string) Value of the option.
787      IF (mm_log) &
788      WRITE(*,'(a,a,a)') "warning: Parameter "//string//"not given... Using default value: "//value
789    END SUBROUTINE printw
790  END FUNCTION mm_global_init_0
791
792  FUNCTION mm_global_init_1(cfg) RESULT(err)
793    !! Set global configuration from a configuration file.
794    !!
795    !! See [[mm_globals(module):mm_global_init_0(function)]].
796    TYPE(cfgparser), INTENT(in) :: cfg  !! Configuration file path.
797    TYPE(error) :: err                  !! Error status of the function.
798    INTEGER                                           :: i
799    TYPE(cfgparser)                                   :: spccfg
800    CHARACTER(len=st_slen)                            :: spcpath
801    CHARACTER(len=:), ALLOCATABLE                     :: defmsg
802    CHARACTER(len=st_slen), DIMENSION(:), ALLOCATABLE :: species
803    REAL(kind=mm_wp)                                  :: zfiamin,zfiamax
804    LOGICAL                                           :: zwhp,zwhs,zwhc,zwcs,zwcn,znofia, &
805                                                         zwstom0,zwstom3
806
807    err = noerror
808
809    IF (mm_ini) THEN
810      err = error("mm_global_init: YAMMS global initialization already performed !",-1)
811      RETURN
812    ENDIF
813
814    ! MP2M mandatory parameters
815    err = mm_check_opt(cfg_get_value(cfg,"df",mm_df),mm_df,wlog=mm_log)
816    IF (err/=0) RETURN
817    err = mm_check_opt(cfg_get_value(cfg,"rm",mm_rm),mm_rm,wlog=mm_log)
818    IF (err/=0) RETURN
819    err = mm_check_opt(cfg_get_value(cfg,"rho_aer",mm_rhoaer),mm_rhoaer,wlog=mm_log)
820    IF (err/=0) RETURN
821    err = mm_check_opt(cfg_get_value(cfg,"p_prod",mm_p_prod),mm_p_prod,wlog=mm_log)
822    IF (err/=0) RETURN
823    err = mm_check_opt(cfg_get_value(cfg,"tx_prod",mm_tx_prod),mm_tx_prod,wlog=mm_log)
824    IF (err/=0) RETURN
825    err = mm_check_opt(cfg_get_value(cfg,"rc_prod",mm_rc_prod),mm_rc_prod,wlog=mm_log)
826    IF (err/=0) RETURN
827    err = mm_check_opt(cfg_get_value(cfg,"planet_radius",mm_rpla),mm_rpla,wlog=mm_log)
828    IF (err/=0) RETURN
829    err = mm_check_opt(cfg_get_value(cfg,"g0",mm_g0),mm_g0,wlog=mm_log)
830    IF (err/=0) RETURN
831    err = mm_check_opt(cfg_get_value(cfg,"timestep",mm_dt),mm_dt,wlog=mm_log)
832    IF (err/=0) RETURN
833    err = mm_check_opt(cfg_get_value(cfg,"air_radius",mm_air_rad),mm_air_rad,wlog=mm_log)
834    IF (err/=0) RETURN
835    err = mm_check_opt(cfg_get_value(cfg,"air_molarmass",mm_air_mmol),mm_air_mmol,wlog=mm_log)
836    IF (err/=0) RETURN
837    err = mm_check_opt(cfg_get_value(cfg,"haze_coag_interactions",mm_coag_choice),mm_coag_choice,wlog=mm_log)
838    IF (err/=0) RETURN
839    err = mm_check_opt(cfg_get_value(cfg,"clouds_microphysics",mm_w_clouds),mm_w_clouds,wlog=mm_log)
840    IF (err/=0) RETURN
841
842    ! computes the conversion factor for bulk -> fractal radius
843    mm_rb2ra = mm_rm**((mm_df-3._mm_wp)/mm_df)
844
845    ! Check coagulation interactions choice
846    IF (mm_coag_choice < 0 .OR. mm_coag_choice > 7) THEN
847      err = error("mm_global_init: Invalid choice for coagulation interactions activation",-1)
848      RETURN
849    ENDIF
850
851    ! Check clouds microphysics input
852    ! it is read only if clouds is activated. We must to check if it is self-consistent...
853    IF (mm_w_clouds) THEN
854      ! Gets species property file path
855      err = cfg_get_value(cfg,'specie_cfg',spcpath) ; IF (err /= 0) RETURN
856      ! Reads species properties configuration file 
857      err = cfg_read_config(spccfg,trim(spcpath)) ; IF (err /= 0) RETURN
858      err = cfg_get_value(spccfg,"used_species",species)
859      IF (err /= 0) THEN
860        err = error("mm_global_init: cannot retrieve 'used_species' values",-1)
861        RETURN
862      ENDIF
863      ! Now attempts to find specides properties !!!
864      mm_nesp = SIZE(species)
865      ALLOCATE(mm_spcname(mm_nesp),mm_xESPS(mm_nesp))
866      !mm_spcname(1:mm_nesp) = species(:)
867      DO i=1,mm_nesp
868        mm_spcname(i) = str_to_lower(species(i))
869        IF (.NOT.cfg_has_section(spccfg,TRIM(mm_spcname(i)))) THEN
870          err = error("mm_global_init: Cannot find "//TRIM(mm_spcname(i))//" properties",-1)
871          RETURN
872        ELSE
873          err = read_esp(spccfg,TRIM(mm_spcname(i)),mm_xESPS(i))
874          ! compute conversion factor: mol.mol-1 => kg.kg-1
875          mm_xESPS(i)%fmol2fmas = mm_xESPS(i)%masmol / mm_air_mmol
876          IF (err/=0) THEN
877            err = error(TRIM(mm_spcname(i))//": "//TRIM(err%msg),-2)
878            RETURN
879          ENDIF
880        ENDIF
881      ENDDO
882    ENDIF
883
884    zwhp = .true. ; zwhs = .true. ; zwhc = .true.
885    zwcs = mm_w_clouds ; zwcn = mm_w_clouds
886    znofia = .false. ; zfiamin = 0.1_mm_wp ; zfiamax = 10._mm_wp
887    zwstom0 = .true. ; zwstom3 = .false.
888
889    ! MP2M Optional parameters
890    err = mm_check_opt(cfg_get_value(cfg,"haze_production",mm_w_haze_prod),mm_w_haze_prod,zwhp,wlog=mm_log)
891    err = mm_check_opt(cfg_get_value(cfg,"haze_sedimentation",mm_w_haze_sed),mm_w_haze_sed,zwhs,wlog=mm_log)
892    err = mm_check_opt(cfg_get_value(cfg,"haze_coagulation",mm_w_haze_coag),mm_w_haze_coag,zwhc,wlog=mm_log)
893    err = mm_check_opt(cfg_get_value(cfg,"clouds_sedimentation",mm_w_cloud_sed),mm_w_cloud_sed,zwcs,wlog=mm_log)
894    err = mm_check_opt(cfg_get_value(cfg,"clouds_nucl_cond",mm_w_cloud_nucond),mm_w_cloud_nucond,zwcn,wlog=mm_log)
895    err = mm_check_opt(cfg_get_value(cfg,"wsed_m0",mm_wsed_m0),mm_wsed_m0,zwstom0,wlog=mm_log)
896    err = mm_check_opt(cfg_get_value(cfg,"wsed_m3",mm_wsed_m3),mm_wsed_m3,zwstom3,wlog=mm_log)
897    err = mm_check_opt(cfg_get_value(cfg,"no_fiadero",mm_no_fiadero_w),mm_no_fiadero_w,znofia,wlog=mm_log)
898    err = mm_check_opt(cfg_get_value(cfg,"fiadero_min_ratio",mm_fiadero_min),mm_fiadero_min,zfiamin,wlog=mm_log)
899    err = mm_check_opt(cfg_get_value(cfg,"fiadero_max_ratio",mm_fiadero_max),mm_fiadero_max,zfiamax,wlog=mm_log)
900
901    err = noerror
902    ! special check for settling velocity
903    IF (mm_wsed_m0 .AND. mm_wsed_m3) THEN
904      err = error("'wsed_m0' and 'wsed_m3' options are mutually exclusive",-1)
905    ENDIF
906    mm_ini = err == noerror
907  END FUNCTION mm_global_init_1
908
909  FUNCTION mm_column_init(plev,zlev,play,zlay,temp) RESULT(err)
910    !! Initialize vertical atmospheric fields.
911    !!
912    !! This subroutine initializes vertical fields needed by the microphysics:
913    !!
914    !! 1. Save reversed input field into "local" array
915    !! 2. Compute thicknesses layers and levels
916    !! 3. Interpolate temperature at levels
917    !!
918    !! The method should be called whenever the vertical structure of the atmosphere changes.
919    !!
920    !! @attention
921    !! All the input vectors should be defined from __GROUND__ to __TOP__ of the atmosphere,
922    !! otherwise nasty things will occur in computations.
923    REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: plev !! Pressure levels (Pa).
924    REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: zlev !! Altitude levels (m).
925    REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: play !! Pressure layers (Pa).
926    REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: zlay !! Altitude at the center of each layer (m).
927    REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: temp !! Temperature at the center of each layer (K).
928    TYPE(error) :: err                                 !! Error status of the function.
929    INTEGER :: i
930    mm_ini_col = .false.                         
931    err = noerror
932    IF (.NOT.mm_ini) THEN
933      err = error("mm_column_init: Global initialization not done yet",-1)
934      RETURN
935    ENDIF
936    IF (mm_nla < 0) THEN
937      mm_nla = SIZE(play)
938    ELSE
939      IF (mm_nla /= SIZE(play)) THEN
940        err = error("mm_column_init: mm_nla cannot be modified dynamically within the run",-1)
941        RETURN
942      ENDIF
943    ENDIF
944    IF (mm_nle < 0) THEN
945      mm_nle = SIZE(plev)
946    ELSE
947      IF (mm_nle /= SIZE(plev)) THEN
948        err = error("mm_column_init: mm_nle cannot be modified dynamically within the run",-1)
949        RETURN
950      ENDIF
951    ENDIF
952    ! should be trashed soon or later
953    IF (mm_nla+1 /= mm_nle) THEN
954      err = error("mm_column_init: Inconsistent number of layers/levels",-1)
955      RETURN
956    ENDIF
957    ! Allocates if required
958    IF (.NOT.ALLOCATED(mm_plev))   ALLOCATE(mm_plev(mm_nle))
959    IF (.NOT.ALLOCATED(mm_zlev))   ALLOCATE(mm_zlev(mm_nle))
960    IF (.NOT.ALLOCATED(mm_play))   ALLOCATE(mm_play(mm_nla))
961    IF (.NOT.ALLOCATED(mm_zlay))   ALLOCATE(mm_zlay(mm_nla))
962    IF (.NOT.ALLOCATED(mm_temp))   ALLOCATE(mm_temp(mm_nla))
963    IF (.NOT.ALLOCATED(mm_btemp))  ALLOCATE(mm_btemp(mm_nle))
964    IF (.NOT.ALLOCATED(mm_dzlev))  ALLOCATE(mm_dzlev(mm_nla))
965    IF (.NOT.ALLOCATED(mm_dzlay))  ALLOCATE(mm_dzlay(mm_nla))
966    IF (.NOT.ALLOCATED(mm_rhoair)) ALLOCATE(mm_rhoair(mm_nla))
967    ! Saves reversed input vectors
968    mm_zlay = zlay(mm_nla:1:-1) ; mm_zlev = zlev(mm_nle:1:-1)
969    mm_play = play(mm_nla:1:-1) ; mm_plev = plev(mm_nle:1:-1)
970    mm_temp = temp(mm_nla:1:-1)
971    ! Computes others vectors
972    mm_dzlay(1:mm_nla-1) = mm_zlay(1:mm_nla-1)-mm_zlay(2:mm_nla)
973    mm_dzlay(mm_nla)     = mm_dzlay(mm_nla-1) ! actually arbitrary
974    mm_dzlev(1:mm_nla)   = mm_zlev(1:mm_nle-1)-mm_zlev(2:mm_nle)
975    mm_btemp(2:mm_nla)   = (mm_temp(1:mm_nla-1)+mm_temp(2:mm_nla))/2._mm_wp
976    mm_btemp(1)          = mm_temp(1)
977    mm_btemp(mm_nle)     = mm_temp(mm_nla)+(mm_temp(mm_nla)-mm_temp(mm_nla-1))/2._mm_wp
978    ! Hydrostatic equilibrium
979    mm_rhoair(1:mm_nla) = (mm_plev(2:mm_nle)-mm_plev(1:mm_nla)) / &
980                          (mm_effg(mm_zlay)*mm_dzlev)
981    mm_ini_col = .true.                         
982    ! write out profiles (only if BOTH debug and log are enabled).
983    IF (mm_log.AND.mm_debug) THEN
984      WRITE(*,'(a)') '# TEMP             PLAY             ZLAY             DZLAY            RHOAIR'
985      DO i=1,mm_nla
986        WRITE(*,'(5(ES15.7,2X))') mm_temp(i),mm_play(i),mm_zlay(i),mm_dzlay(i), mm_rhoair(i)
987      ENDDO
988      WRITE(*,'(a)') '# TEMP             PLEV             ZLEV             DZLEV'
989      DO i=1,mm_nle
990        IF (i /= mm_nle) THEN
991          WRITE(*,'(4(ES15.7,2X))') mm_btemp(i),mm_plev(i),mm_zlev(i),mm_dzlev(i)
992        ELSE
993          WRITE(*,'(3(ES15.7,2X))') mm_btemp(i),mm_plev(i),mm_zlev(i)
994        ENDIF
995      ENDDO
996    ENDIF
997
998    RETURN
999  END FUNCTION mm_column_init
1000
1001  FUNCTION mm_aerosols_init(m0aer_s,m3aer_s,m0aer_f,m3aer_f) RESULT(err)
1002    !! Initialize clouds tracers vertical grid.
1003    !!
1004    !! The subroutine initializes aerosols microphysics tracers columns. It allocates variables if
1005    !! required and stores input vectors in reversed order. It also computes the characteristic radii
1006    !! of each mode.
1007    !! @note
1008    !! All the input arguments should be defined from ground to top.
1009    !!
1010    !! @attention
1011    !! [[mm_globals(module):mm_global_init(interface)]] and [[mm_globals(module):mm_column_init(function)]]
1012    !! must have been called at least once before this method is called. Moreover, this method should be
1013    !! after each call of [[mm_globals(module):mm_column_init(function)]] to reflect changes in the
1014    !! vertical atmospheric structure.
1015    REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m0aer_s !! \(0^{th}\) order moment of the spherical mode (\(m^{-2}\)).
1016    REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m3aer_s !! \(3^{rd}\) order moment of the spherical mode (\(m^{3}.m^{-2}\)).
1017    REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m0aer_f !! \(0^{th}\) order moment of the fractal mode (\(m^{-2}\)).
1018    REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m3aer_f !! \(3^{rd}\) order moment of the fractal mode (\(m^{3}.m^{-2}\)).
1019    TYPE(error) :: err                                    !! Error status of the function.
1020    err = noerror
1021    IF (.NOT.mm_ini) THEN
1022      err = error("mm_aerosols_init: Global initialization not done yet",-1) ; RETURN
1023    ENDIF
1024    IF (.NOT.mm_ini_col) THEN
1025      err = error("mm_aerosols_init: Column initialization not done yet",-1) ; RETURN
1026    ENDIF
1027    ! Check input size ???
1028    IF (SIZE(m0aer_s) /= mm_nla) THEN
1029      err = error("mm_aerosols_init: Invalid size for input arrays",-1) ; RETURN
1030    ENDIF
1031
1032    ! Allocate variable if required
1033    IF (.NOT.ALLOCATED(mm_m0aer_s)) ALLOCATE(mm_m0aer_s(mm_nla))
1034    IF (.NOT.ALLOCATED(mm_m3aer_s)) ALLOCATE(mm_m3aer_s(mm_nla))
1035    IF (.NOT.ALLOCATED(mm_m0aer_f)) ALLOCATE(mm_m0aer_f(mm_nla))
1036    IF (.NOT.ALLOCATED(mm_m3aer_f)) ALLOCATE(mm_m3aer_f(mm_nla))
1037    IF (.NOT.ALLOCATED(mm_rcs))     ALLOCATE(mm_rcs(mm_nla))
1038    IF (.NOT.ALLOCATED(mm_rcf))     ALLOCATE(mm_rcf(mm_nla))
1039    ! Allocate memory for diagnostics
1040    IF (.NOT.ALLOCATED(mm_aer_s_flux)) THEN
1041      ALLOCATE(mm_aer_s_flux(mm_nla)) ; mm_aer_s_flux(:) = 0._mm_wp
1042    ENDIF
1043    IF (.NOT.ALLOCATED(mm_aer_f_flux)) THEN
1044      ALLOCATE(mm_aer_f_flux(mm_nla)) ; mm_aer_f_flux(:) = 0._mm_wp
1045    ENDIF
1046    IF (.NOT.ALLOCATED(mm_m0as_vsed)) THEN
1047      ALLOCATE(mm_m0as_vsed(mm_nla)) ; mm_m0as_vsed(:) = 0._mm_wp
1048    ENDIF
1049    IF (.NOT.ALLOCATED(mm_m3as_vsed)) THEN
1050      ALLOCATE(mm_m3as_vsed(mm_nla)) ; mm_m3as_vsed(:) = 0._mm_wp
1051    ENDIF
1052    IF (.NOT.ALLOCATED(mm_m0af_vsed)) THEN
1053      ALLOCATE(mm_m0af_vsed(mm_nla)) ; mm_m0af_vsed(:) = 0._mm_wp
1054    ENDIF
1055    IF (.NOT.ALLOCATED(mm_m3af_vsed)) THEN
1056      ALLOCATE(mm_m3af_vsed(mm_nla)) ; mm_m3af_vsed(:) = 0._mm_wp
1057    ENDIF
1058    ! note : mm_dzlev is already from top to ground
1059    mm_m0aer_s = m0aer_s(mm_nla:1:-1)/mm_dzlev(:)
1060    mm_m3aer_s = m3aer_s(mm_nla:1:-1)/mm_dzlev(:)
1061    mm_m0aer_f = m0aer_f(mm_nla:1:-1)/mm_dzlev(:)
1062    mm_m3aer_f = m3aer_f(mm_nla:1:-1)/mm_dzlev(:)
1063    ! aerosols characteristic radii
1064    ! il faudrait peut etre revoir la gestion des zeros ici et la
1065    ! remplacer par une valeur seuil des moments.
1066    WHERE(mm_m3aer_s > 0._mm_wp .AND. mm_m0aer_s > 0._mm_wp)
1067      mm_rcs = mm_get_rcs(mm_m0aer_s,mm_m3aer_s)
1068    ELSEWHERE
1069      mm_rcs = 0._mm_wp
1070    ENDWHERE
1071    WHERE(mm_m3aer_f > 0._mm_wp .AND. mm_m0aer_f > 0._mm_wp)
1072      mm_rcf = mm_get_rcf(mm_m0aer_f,mm_m3aer_f)
1073    ELSEWHERE
1074      mm_rcf = 0._mm_wp
1075    ENDWHERE
1076    mm_ini_aer = .true.
1077  END FUNCTION mm_aerosols_init
1078
1079  FUNCTION mm_clouds_init(m0ccn,m3ccn,m3ice,gazs) RESULT(err)
1080    !! Initialize clouds tracers vertical grid.
1081    !!
1082    !! The subroutine initializes cloud microphysics tracers columns. It allocates variables if
1083    !! required and stores input vectors in reversed order. It also computes the mean drop radius
1084    !! and density and allocates diagnostic vectors.
1085    !! @note
1086    !! All the input arguments should be defined from ground to top.
1087    !!
1088    !! @attention
1089    !! [[mm_globals(module):mm_global_init(interface)]] and [[mm_globals(module):mm_column_init(function)]]
1090    !! must have been called at least once before this method is called. Moreover, this method should be
1091    !! after each call of [[mm_globals(module):mm_column_init(function)]] to reflect changes in the
1092    !! vertical atmospheric structure.
1093    REAL(kind=mm_wp), DIMENSION(:), INTENT(in)   :: m0ccn !! 0th order moment of the CCN distribution (\(m^{-2}\)).
1094    REAL(kind=mm_wp), DIMENSION(:), INTENT(in)   :: m3ccn !! 3rd order moment of the CCN distribution (\(m^{3}.m^{-2}\)).
1095    REAL(kind=mm_wp), DIMENSION(:,:), INTENT(in) :: m3ice !! 3rd order moments of the ice components (\(m^{3}.m^{-2}\)).
1096    REAL(kind=mm_wp), DIMENSION(:,:), INTENT(in) :: gazs  !! Condensible species gazs molar fraction (\(mol.mol^{-1}\)).
1097    TYPE(error) :: err                                    !! Error status of the function.
1098    INTEGER :: i
1099    err = noerror
1100    IF (.NOT.mm_ini) THEN
1101      err = error("Global initialization not done yet",-8)
1102      RETURN
1103    ENDIF
1104
1105    IF (.NOT.mm_w_clouds) THEN
1106      IF (mm_debug) WRITE(*,'(a)') "WARNING: Cloud microphysic is not enabled..."
1107      RETURN
1108    ENDIF
1109
1110    ! Note:
1111    !  Here we could check that mm_nla is the same size of gazs(DIM=1)
1112    !  Actually, mm_nla should always initialized the first time mm_column_init is called, NOT HERE.
1113    IF (mm_nla < 0)  mm_nla  = SIZE(gazs,DIM=1)
1114    ! Note:
1115    !   here we could check that mm_nesp is the same size of gazs(DIM=2)
1116    !   Actually, mm_nesp should be always initialized in mm_global_init, NOT HERE.
1117    IF (mm_nesp < 0) mm_nesp = SIZE(gazs,DIM=2)
1118
1119    ! Allocate variable if required
1120    IF (.NOT.ALLOCATED(mm_m0ccn))   ALLOCATE(mm_m0ccn(mm_nla))
1121    IF (.NOT.ALLOCATED(mm_m3ccn))   ALLOCATE(mm_m3ccn(mm_nla))
1122    IF (.NOT.ALLOCATED(mm_m3ice))   ALLOCATE(mm_m3ice(mm_nla,mm_nesp))
1123    IF (.NOT.ALLOCATED(mm_gazs))    ALLOCATE(mm_gazs(mm_nla,mm_nesp))
1124    IF (.NOT.ALLOCATED(mm_drad))    ALLOCATE(mm_drad(mm_nla))
1125    IF (.NOT.ALLOCATED(mm_drho))    ALLOCATE(mm_drho(mm_nla))
1126    ! Allocate memory for diagnostics
1127    IF (.NOT.ALLOCATED(mm_ccn_flux)) THEN
1128      ALLOCATE(mm_ccn_flux(mm_nla)) ; mm_ccn_flux(:) = 0._mm_wp
1129    ENDIF
1130    IF (.NOT.ALLOCATED(mm_ice_prec))   THEN
1131      ALLOCATE(mm_ice_prec(mm_nesp)) ; mm_ice_prec(:) = 0._mm_wp
1132    ENDIF
1133    IF (.NOT.ALLOCATED(mm_ice_fluxes)) THEN
1134      ALLOCATE(mm_ice_fluxes(mm_nla,mm_nesp)) ; mm_ice_fluxes(:,:) = 0._mm_wp
1135    ENDIF
1136    IF (.NOT.ALLOCATED(mm_gazs_sat)) THEN
1137      ALLOCATE(mm_gazs_sat(mm_nla,mm_nesp)) ; mm_gazs_sat(:,:) = 0._mm_wp
1138    ENDIF
1139
1140    ! note mm_dzlev already from top to ground
1141    mm_m0ccn = m0ccn(mm_nla:1:-1)/mm_dzlev(:)
1142    mm_m3ccn = m3ccn(mm_nla:1:-1)/mm_dzlev(:)
1143    DO i=1,mm_nesp
1144      mm_m3ice(:,i) = m3ice(mm_nla:1:-1,i)/mm_dzlev(:)
1145      mm_gazs(:,i)  = gazs(mm_nla:1:-1,i)
1146    ENDDO
1147    ! drop mean radius
1148    call mm_cloud_properties(mm_m0ccn,mm_m3ccn,mm_m3ice,mm_drad,mm_drho)
1149    mm_ini_cld = .true.
1150  END FUNCTION mm_clouds_init
1151
1152  SUBROUTINE mm_dump_parameters()
1153    !! Dump model global parameters on stdout.
1154    WRITE(*,'(a)')         "========= YAMMS PARAMETERS ============"
1155    WRITE(*,'(a,L2)')      "mm_w_haze_prod         : ", mm_w_haze_prod
1156    WRITE(*,'(a,ES14.7)')  "   mm_p_prod           : ", mm_p_prod
1157    WRITE(*,'(a,ES14.7)')  "   mm_tx_prod          : ", mm_tx_prod
1158    WRITE(*,'(a,ES14.7)')  "   mm_rc_prod          : ", mm_rc_prod
1159    WRITE(*,'(a,L2)')      "mm_w_haze_coag         : ", mm_w_haze_coag
1160    WRITE(*,'(a,I2.2)')    "   mm_coag_interactions: ", mm_coag_choice
1161    WRITE(*,'(a,L2)')      "mm_w_haze_sed          : ", mm_w_haze_sed
1162    WRITE(*,'(a,L2)')      "   mm_wsed_m0          : ", mm_wsed_m0
1163    WRITE(*,'(a,L2)')      "   mm_wsed_m3          : ", mm_wsed_m3
1164    WRITE(*,'(a,L2)')      "   mm_no_fiadero_w     : ", mm_no_fiadero_w
1165    WRITE(*,'(a,ES14.7)')  "   mm_fiadero_min      : ", mm_fiadero_min
1166    WRITE(*,'(a,ES14.7)')  "   mm_fiadero_max      : ", mm_fiadero_max
1167    WRITE(*,'(a,L2)')      "mm_w_clouds            : ", mm_w_clouds
1168    WRITE(*,'(a,L2)')      "   mm_w_cloud_sed      : ", mm_w_cloud_sed
1169    WRITE(*,'(a,L2)')      "   mm_w_cloud_nucond   : ", mm_w_cloud_nucond
1170    WRITE(*,'(a)')         "---------------------------------------"
1171    WRITE(*,'(a,ES14.7)')  "mm_dt                  : ", mm_dt
1172    IF (mm_nla > -1) THEN
1173      WRITE(*,'(a,I3.3)')    "mm_nla                 : ", mm_nla
1174    ELSE
1175      WRITE(*,'(a)')         "mm_nla                 : not initialized yet"
1176    ENDIF
1177    WRITE(*,'(a,ES14.7)')  "mm_df                  : ", mm_df
1178    WRITE(*,'(a,ES14.7)')  "mm_rm                  : ", mm_rm
1179    WRITE(*,'(a,ES14.7)')  "mm_rpla                : ", mm_rpla
1180    WRITE(*,'(a,ES14.7)')  "mm_g0                  : ", mm_g0
1181    WRITE(*,'(a)')         "======================================="
1182  END SUBROUTINE mm_dump_parameters
1183
1184  ELEMENTAL FUNCTION mm_get_rcs(m0,m3) RESULT(res)
1185    !! Get the characteristic radius for the spherical aerosols size distribution.
1186    !!
1187    !! The method computes the characteristic radius of the size distribution law
1188    !! of the spherical aerosols mode according to its moments and its inter-moments
1189    !! relation.
1190    REAL(kind=mm_wp), INTENT(in) :: m0 !! \(0^{th}\) order moment
1191    REAL(kind=mm_wp), INTENT(in) :: m3 !! \(3^{rd}\) order moment
1192    REAL(kind=mm_wp) :: res            !! Radius
1193    ! arbitrary: if there is no way to compute radius
1194    IF (m3 <= 0._mm_wp .OR. m0 <= 0._mm_wp) res = 1._mm_wp
1195    res = (m3/m0/mm_alpha_s(3._mm_wp))**(1._mm_wp/3._mm_wp)
1196  END FUNCTION mm_get_rcs
1197
1198  ELEMENTAL FUNCTION mm_get_rcf(m0,m3) RESULT(res)
1199    !! Get the characteristic radius for the fractal aerosols size distribution.
1200    !!
1201    !! The method computes the characteristic radius of the size distribution law
1202    !! of the fractal aerosols mode according to its moments and its inter-moments
1203    !! relation.
1204    REAL(kind=mm_wp), INTENT(in) :: m0 !! \(0^{th}\) order moment
1205    REAL(kind=mm_wp), INTENT(in) :: m3 !! \(3^{rd}\) order moment
1206    REAL(kind=mm_wp) :: res            !! Radius
1207    ! arbitrary: if there is no way to compute radius
1208    IF (m3 <= 0._mm_wp .OR. m0 <= 0._mm_wp) res = 1._mm_wp
1209    res = (m3/m0/mm_alpha_f(3._mm_wp))**(1._mm_wp/3._mm_wp)
1210  END FUNCTION mm_get_rcf
1211
1212  ELEMENTAL FUNCTION mm_effg(z) RESULT(effg)
1213    !! Compute effective gravitational acceleration.
1214    REAL(kind=mm_wp), INTENT(in) :: z !! Altitude in meters
1215    REAL(kind=mm_wp) :: effg          !! Effective gravitational acceleration in \(m.s^{-2}\)
1216    effg = mm_g0 * (mm_rpla/(mm_rpla+z))**2
1217    RETURN
1218  END FUNCTION mm_effg
1219
1220  !==================================
1221  ! --- private methods -------------
1222  !==================================
1223
1224  SUBROUTINE cldprop_sc(m0ccn,m3ccn,m3ice,drad,drho)
1225    !! Get cloud drop properties (scalar).
1226    !!
1227    !! The method computes the mean radius and mean density of cloud drops.
1228    !!
1229    !! @bug
1230    !! A possible bug can happen because of threshold snippet. If __drad__ is greater than
1231    !! __drmax__ (== 100 microns) it is automatically set to __drmax__, but computation of
1232    !! __drho__ remains unmodified. So __drho__ is not correct in that case.
1233    !!
1234    !! @todo
1235    !! Fix the bug of the subroutine, but it is rather minor, since theoretically we do not
1236    !! need the density of the drop.
1237    !!
1238    !! @todo
1239    !! Think about a better implementation of thresholds, and get sure of their consequences in
1240    !! the other parts of the model.
1241    REAL(kind=mm_wp), INTENT(in)               :: m0ccn !! \(0^{th}\) order moment of the ccn
1242    REAL(kind=mm_wp), INTENT(in)               :: m3ccn !! \(3^{rd}\) order moment of the ccn
1243    REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: m3ice !! \(3^{rd}\) order moments of each ice component
1244    REAL(kind=mm_wp), INTENT(out)              :: drad  !! Output mean drop radius
1245    REAL(kind=mm_wp), INTENT(out), OPTIONAL    :: drho  !! Optional output mean drop density
1246    REAL(kind=mm_wp)            :: vtot, wtot, ntot
1247    REAL(kind=mm_wp), PARAMETER :: threshold = 1.e-25_mm_wp,         & 
1248                                       drmin = 1.e-10_mm_wp,         &
1249                                       drmax = 1.e-4_mm_wp,          &
1250                                      athird = 1._mm_wp/3._mm_wp,    &
1251                                       pifac = 4._mm_wp/3._mm_wp*mm_pi
1252    drad = 0._mm_wp
1253    ntot = m0ccn
1254    vtot = pifac*m3ccn+SUM(m3ice)
1255    wtot = pifac*m3ccn*mm_rhoaer+SUM(m3ice*mm_xESPS(:)%rho)
1256    IF (ntot <= threshold .OR. vtot <= 0._mm_wp) THEN
1257      drad  = drmin
1258      IF (PRESENT(drho)) drho  = mm_rhoaer
1259    ELSE
1260      drad = (vtot/ntot/pifac)**athird
1261      drad = MAX(MIN(drad,drmax),drmin)
1262      IF (PRESENT(drho)) drho = wtot/vtot
1263    ENDIF
1264    RETURN
1265  END SUBROUTINE cldprop_sc
1266
1267  SUBROUTINE cldprop_ve(m0ccn,m3ccn,m3ice,drad,drho)
1268    !! Get cloud drop properties (vector).
1269    !!
1270    !! The method performs the same computations than [[mm_globals(module):cldprop_sc(subroutine)]]
1271    !! but for the entire vertical atmospheric structure.
1272    !! Same remarks apply here.
1273    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)            :: m0ccn !! 0th order moment of the ccn.
1274    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)            :: m3ccn !! 3rd order moment of the ccn.
1275    REAL(kind=mm_wp), INTENT(in), DIMENSION(:,:)          :: m3ice !! 3rd order moments of each ice component.
1276    REAL(kind=mm_wp), INTENT(out), DIMENSION(:)           :: drad  !! Output mean drop radius.
1277    REAL(kind=mm_wp), INTENT(out), DIMENSION(:), OPTIONAL :: drho  !! Optional output mean drop density.
1278    INTEGER                                     :: i,ns
1279    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: vtot,wtot,ntot,rho
1280    REAL(kind=mm_wp), PARAMETER                 :: threshold = 1.e-25_mm_wp,         & 
1281                                                       drmin = 1.e-10_mm_wp,         &
1282                                                       drmax = 1.e-4_mm_wp,          &
1283                                                      athird = 1._mm_wp/3._mm_wp,    &
1284                                                       pifac = 4._mm_wp/3._mm_wp*mm_pi
1285                                                     
1286    ns = SIZE(m0ccn) ; ALLOCATE(vtot(ns),wtot(ns),ntot(ns),rho(ns))
1287    drad = 0._mm_wp
1288    ntot = m0ccn
1289    vtot = pifac*m3ccn+SUM(m3ice,DIM=2)
1290    wtot = pifac*m3ccn*mm_rhoaer
1291    DO i=1,SIZE(m3ice,DIM=2)
1292      wtot = wtot+m3ice(:,i)*mm_xESPS(i)%rho
1293    ENDDO
1294    WHERE(ntot <= threshold .OR. vtot <= 0._mm_wp)
1295      drad  = drmin
1296      rho = mm_rhoaer
1297    ELSEWHERE
1298      drad = (vtot/ntot/pifac)**athird
1299      drad = MAX(MIN(drad,drmax),drmin)
1300      rho = wtot/vtot
1301    END WHERE
1302    IF (PRESENT(drho)) drho  = rho
1303    RETURN
1304  END SUBROUTINE cldprop_ve
1305
1306  ! For configuration file (requires fccp library).
1307
1308  FUNCTION read_esp(parser,sec,pp) RESULT (err)
1309    !! Read and store [[mm_globals(module):mm_esp(type)]] parameters.
1310    TYPE(cfgparser), INTENT(in)   :: parser !! Configuration parser.
1311    CHARACTER(len=*), INTENT(in)  :: sec    !! Name of the specie (should match a section of the configuration.
1312    TYPE(mm_esp), INTENT(out)     :: pp     !! [[mm_globals(module):mm_esp(type)]] object that stores the parameters.
1313    TYPE(error)                   :: err    !! Error status of the function.
1314    err = cfg_get_value(parser,TRIM(sec)//'name',pp%name)       ; IF (err /= 0) RETURN
1315    err = cfg_get_value(parser,TRIM(sec)//'mas',pp%mas)         ; IF (err /= 0) RETURN
1316    err = cfg_get_value(parser,TRIM(sec)//'vol',pp%vol)         ; IF (err /= 0) RETURN
1317    err = cfg_get_value(parser,TRIM(sec)//'ray',pp%ray)         ; IF (err /= 0) RETURN
1318    err = cfg_get_value(parser,TRIM(sec)//'mas',pp%mas)         ; IF (err /= 0) RETURN
1319    err = cfg_get_value(parser,TRIM(sec)//'vol',pp%vol)         ; IF (err /= 0) RETURN
1320    err = cfg_get_value(parser,TRIM(sec)//'ray',pp%ray)         ; IF (err /= 0) RETURN
1321    err = cfg_get_value(parser,TRIM(sec)//'masmol',pp%masmol)   ; IF (err /= 0) RETURN
1322    err = cfg_get_value(parser,TRIM(sec)//'rho',pp%rho)         ; IF (err /= 0) RETURN
1323    err = cfg_get_value(parser,TRIM(sec)//'tc',pp%tc)           ; IF (err /= 0) RETURN
1324    err = cfg_get_value(parser,TRIM(sec)//'pc',pp%pc)           ; IF (err /= 0) RETURN
1325    err = cfg_get_value(parser,TRIM(sec)//'tb',pp%tb)           ; IF (err /= 0) RETURN
1326    err = cfg_get_value(parser,TRIM(sec)//'w',pp%w)             ; IF (err /= 0) RETURN
1327    err = cfg_get_value(parser,TRIM(sec)//'a_sat',pp%a_sat)     ; IF (err /= 0) RETURN
1328    err = cfg_get_value(parser,TRIM(sec)//'b_sat',pp%b_sat)     ; IF (err /= 0) RETURN
1329    err = cfg_get_value(parser,TRIM(sec)//'c_sat',pp%c_sat)     ; IF (err /= 0) RETURN
1330    err = cfg_get_value(parser,TRIM(sec)//'d_sat',pp%d_sat)     ; IF (err /= 0) RETURN
1331    err = cfg_get_value(parser,TRIM(sec)//'mteta',pp%mteta)     ; IF (err /= 0) RETURN
1332    err = cfg_get_value(parser,TRIM(sec)//'tx_prod',pp%tx_prod) ; IF (err /= 0) RETURN
1333    RETURN
1334  END FUNCTION read_esp
1335
1336  ! =========================================================================
1337  ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1338  !                CONFIGURATION PARSER checking methods
1339  ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1340  ! =========================================================================
1341
1342  FUNCTION check_r1(err,var,def,wlog) RESULT(ret)
1343    !! Check an option value (float).
1344    !!
1345    !! The method checks an option value and optionally set a default value, __def__ to initialize
1346    !! __var__ on error if given.
1347    TYPE(error), INTENT(in)                :: err  !! Error object from value getter.
1348    REAL(kind=mm_wp), INTENT(inout)        :: var  !! Input/output option value.
1349    REAL(kind=mm_wp), INTENT(in), OPTIONAL :: def  !! Default value to set.
1350    LOGICAL, INTENT(in), OPTIONAL          :: wlog !! .true. to print warning/error message.
1351    TYPE(error) :: ret                             !! Input error.
1352    CHARACTER(len=*), PARAMETER :: defmsg = '... Using default value: '
1353    LOGICAL                     :: zlog
1354    ret = err
1355    zlog = .false. ; IF (PRESENT(wlog)) zlog = wlog
1356    IF (err == 0) RETURN
1357    IF (PRESENT(def)) THEN
1358      var = def
1359      IF (zlog) WRITE(*,'(a,a,a)') error_to_string(err,'',.true.),defmsg,to_string(var)
1360      ret = noerror
1361    ELSE
1362      IF (zlog) WRITE(*,'(a)') error_to_string(err,'',.true.)
1363    ENDIF
1364  END FUNCTION check_r1
1365
1366  FUNCTION check_l1(err,var,def,wlog) RESULT(ret)
1367    !! Check an option value (logical).
1368    !!
1369    !! The method checks an option value and optionally set a default value, __def__ to initialize
1370    !! __var__ on error if given.
1371    TYPE(error), INTENT(in)       :: err  !! Error object from value getter.
1372    LOGICAL, INTENT(inout)        :: var  !! Input/output option value.
1373    LOGICAL, INTENT(in), OPTIONAL :: def  !! Default value to set.
1374    LOGICAL, INTENT(in), OPTIONAL :: wlog !! .true. to print warning/error message.
1375    TYPE(error) :: ret                    !! Input error.
1376    CHARACTER(len=*), PARAMETER :: defmsg = '... Using default value: '
1377    LOGICAL                     :: zlog
1378    ret = err
1379     zlog = .false. ; IF (PRESENT(wlog)) zlog = wlog
1380    IF (err == 0) RETURN
1381    IF (PRESENT(def)) THEN
1382      var = def
1383      IF (zlog) WRITE(*,'(a,a,a)') error_to_string(err,'',.true.),defmsg,to_string(var)
1384      ret = noerror
1385    ELSE
1386      IF (zlog) WRITE(*,'(a)') error_to_string(err,'',.true.)
1387    ENDIF
1388  END FUNCTION check_l1
1389
1390  FUNCTION check_i1(err,var,def,wlog) RESULT(ret)
1391    !! Check an option value (integer).
1392    !!
1393    !! The method checks an option value and optionally set a default value, __def__ to initialize
1394    !! __var__ on error if given.
1395    TYPE(error), INTENT(in)       :: err  !! Error object from value getter.
1396    INTEGER, INTENT(inout)        :: var  !! Input/output option value.
1397    INTEGER, INTENT(in), OPTIONAL :: def  !! Default value to set.
1398    LOGICAL, INTENT(in), OPTIONAL :: wlog !! .true. to print warning/error message.
1399    TYPE(error) :: ret                    !! Input error.
1400    CHARACTER(len=*), PARAMETER :: defmsg = '... Using default value: '
1401    LOGICAL                     :: zlog
1402    ret = err
1403     zlog = .false. ; IF (PRESENT(wlog)) zlog = wlog
1404    IF (err == 0) RETURN
1405    IF (PRESENT(def)) THEN
1406      var = def
1407      IF (zlog) WRITE(*,'(a,a,a)') error_to_string(err,'',.true.),defmsg,to_string(var)
1408      ret = noerror
1409    ELSE
1410      IF (zlog) WRITE(*,'(a)') error_to_string(err,'',.true.)
1411    ENDIF
1412  END FUNCTION check_i1
1413
1414  FUNCTION check_s1(err,var,def,wlog) RESULT(ret)
1415    !! Check an option value (string).
1416    !!
1417    !! The method checks an option value and optionally set a default value, __def__ to initialize
1418    !! __var__ on error if given.
1419    TYPE(error), INTENT(in)                      :: err  !! Error object from value getter.
1420    CHARACTER(len=*), INTENT(inout)              :: var  !! Input/output option value.
1421    CHARACTER(len=*), INTENT(in), OPTIONAL       :: def  !! Default value to set.
1422    LOGICAL, INTENT(in), OPTIONAL                :: wlog !! .true. to print warning/error message.
1423    TYPE(error) :: ret                                   !! Input error.
1424    CHARACTER(len=*), PARAMETER :: defmsg = '... Using default value: '
1425    LOGICAL                     :: zlog
1426    ret = err
1427    zlog = .false. ; IF (PRESENT(wlog)) zlog = wlog
1428    IF (err == 0) RETURN
1429    IF (PRESENT(def)) THEN
1430      var = TRIM(def)
1431      IF (zlog) WRITE(*,'(a,a,a)') error_to_string(err,'',.true.),defmsg,var
1432      ret = noerror
1433    ELSE
1434      IF (zlog) WRITE(*,'(a)') error_to_string(err,'')
1435    ENDIF
1436    RETURN
1437  END FUNCTION check_s1
1438
1439END MODULE MM_GLOBALS
Note: See TracBrowser for help on using the repository browser.