! Copyright (2013-2015,2017,2022-2023) Université de Reims Champagne-Ardenne ! Contributors : J. Burgalat (GSMA, URCA), B. de Batz de Trenquelléon (GSMA, URCA) ! email of the author : jeremie.burgalat@univ-reims.fr ! ! This software is a computer program whose purpose is to compute ! microphysics processes using a two-moments scheme. ! ! This library is governed by the CeCILL-B license under French law and ! abiding by the rules of distribution of free software. You can use, ! modify and/ or redistribute the software under the terms of the CeCILL-B ! license as circulated by CEA, CNRS and INRIA at the following URL ! "http://www.cecill.info". ! ! As a counterpart to the access to the source code and rights to copy, ! modify and redistribute granted by the license, users are provided only ! with a limited warranty and the software's author, the holder of the ! economic rights, and the successive licensors have only limited ! liability. ! ! In this respect, the user's attention is drawn to the risks associated ! with loading, using, modifying and/or developing or reproducing the ! software by the user in light of its specific status of free software, ! that may mean that it is complicated to manipulate, and that also ! therefore means that it is reserved for developers and experienced ! professionals having in-depth computer knowledge. Users are therefore ! encouraged to load and test the software's suitability as regards their ! requirements in conditions enabling the security of their systems and/or ! data to be ensured and, more generally, to use and operate it in the ! same conditions as regards security. ! ! The fact that you are presently reading this means that you have had ! knowledge of the CeCILL-B license and that you accept its terms. !! file: mm_globals.f90 !! summary: Parameters and global variables module. !! author: J. Burgalat !! date: 2013-2015,2017,2022-2023 !! modifications: B. de Batz de Trenquelléon MODULE MM_GLOBALS !! Parameters and global variables module. !! !! # Module overview !! !! The module defines all the parameters and global variables that are common !! to all other modules of the library. !! !! It is separated in two parts : !! !! - Main parameters and global saved variables. Most of these variables should !! be initialized once and should hold the same value during run-time. These !! variables are completly public and initialized by [[mm_globals(module):mm_global_init(interface)]] !! method. !! - The second part defines a set of vectors that defines the vertical structure of the atmosphere. !! Each time a new atmospheric column has to be computed (either on a new timestep or on a new couple !! of longitude/latitude), these vectors should be intialized with new values by calling !! [[mm_globals(module):mm_column_init(function)]] method. !! This part is separated in two sets : !! !! - The atmospheric structure with temperature, pressure levels and altitude definitions. !! - The vertical profiles of tracers with the moments of the two aerosols modes (both \(M_{0}\) !! and \(M_{3}\) for a total of 4 vectors), the _clouds_ microphysics moments tracers (i.e. !! \(M_{0}\) and \(M_{3}\) for the ccn and \(M_{3}\) for the ice components). !! Additionally, the module also stores intermediates variables of interest such as the !! characteristic radii of the aerosols modes, the mean drop radius and the drop density, !! the molar fraction of each condensible species (related to ice components) and some !! scalar variables that holds arrays sizes. !! !! @note !! All the vectors that represent the vertical structure of the atmosphere (altitude, pressure and !! temperature...) are oriented from the __TOP__ of the atmosphere to the __GROUND__. !! !! @note !! The module also imports errors module from __SWIFT__ library to get definitions of the error object !! everywhere in the library ([[mm_globals(module)]] is always imported, except in [[mm_mprec(module)]]). !! !! # Global variables !! !! [[mm_globals(module)]] module contains the declaration of all global/common variable that are shared !! by all other modules of the model. Except for few physical constant which are declared as parameters, !! these variables are onlu SAVEd. They are initialized by [[mm_globals(module):mm_global_init(interface)]] !! methods. !! the following sections list all the global variables by category. !! !! ## Control flags !! !! | Name | Description !! | :----------------- | :----------------- !! | mm_log | Enable log mode (verbose) !! | mm_w_haze_prod | Enable/Disable haze production !! | mm_w_haze_sed | Enable/Disable haze sedimentation !! | mm_w_haze_coag | Enable/Disable haze coagulation !! | mm_w_clouds   | Enable/Disable clouds microphysics !! | mm_w_clouds_sed  | Enable/Disable clouds microphysics sedimentation !! | mm_w_clouds_nucond | Enable/Disable clouds microphysics nucleation/condensation !! | mm_wsed_m0 | Force all aerosols moments to fall at M0 settling velocity !! | mm_wsed_m3 | Force all aerosols moments to fall at M3 settling velocity !! | mm_no_fiadero_w | Enable/Disable __Fiadero__ correction !! !! ### Related free parameters: !! !! | Name | Description !! | :-------------- | :----------------- !! | mm_fiadero_min | Minimum ratio for __Fiadero__'s correction !! | mm_fiadero_max | Maximum ratio for __Fiadero__'s correction !! | 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)]]. !! !! ## Physical constants !! !! | Name | Description !! | :-------- | :----------------- !! | mm_pi | Pi number !! | mm_navo | Avogadro number !! | mm_kboltz | Boltzmann constant (\(J.K^{-1}\)) !! | mm_rgas | Perfect gas constant (\(J.mol^{-1}.K^{-1}\)) !! | mm_fdes | Desorption energy (\(J\)) (nucleation) !! | mm_fdif | Surface diffusion energy (\(J\)) (nucleation) !! | mm_fnus | Jump frequency (\(s^{-1}\)) (nucleation) !! | mm_akn | Approximated slip-flow correction coefficient ( !! !! ## Free parameters !! !! | Name | Description !! | :---------- | :----------------- !! | mm_rhoaer | Aerosol density (in \(kg.m^{-3}\)) !! | mm_df | Fractal dimension !! | mm_rm | Monomer radius (in m) !! | mm_p_prod | Spherical aerosols production pressure level (Pa) !! | mm_p_rcprod | Spherical aerosols equivalent radius production (m) !! | mm_tx_prod | Production rate of spherical aerosols (\(kg.m^{-2}.s^{-1}\)) !! | mm_d_prod | Time-dependent sine wve pre-factor. !! | mm_w_prod | Angular frequency of the time-dependent production rate. !! | mm_ne | Electric charging of aerosols (\(e^{-}.m^{-1}\)) (unused) !! | mm_rb2ra  | Bulk to apparent radius conversion pre-factor (\(m^X\)) !! | mm_rpla | Planet radius (m) !! | mm_g0 | Planet acceleration due to gravity constant (ground) (\(m.s^{-2}\)) !! | mm_air_rad | Air molecules mean radius (m) !! | mm_air_mmol | Air molecules molar mass (\(kg.mol^{-1}\)) !! | mm_dt | Microphysic time step (s) USE MM_MPREC USE MM_INTERFACES ! from swift USE CFGPARSE USE STRING_OP USE ERRORS IMPLICIT NONE PUBLIC PRIVATE :: cldprop_sc,cldprop_ve,read_esp,check_r1,check_i1,check_l1,check_s1 ! Protected variables ! the following variables are read-only outside this module. ! One must call the afferent subroutine to update them. ! initialization control flags (cannot be updated) PROTECTED :: mm_ini,mm_ini_col,mm_ini_aer,mm_ini_cld ! model parameters (mm_global_init) PROTECTED :: mm_dt,mm_rhoaer,mm_df,mm_rm,mm_p_prod,mm_rc_prod,mm_tx_prod,mm_rpla,mm_g0,mm_rb2ra ! atmospheric vertical structure (mm_column_init) PROTECTED :: mm_nla,mm_nle,mm_zlay,mm_zlev,mm_play,mm_plev,mm_temp,mm_rhoair,mm_btemp,mm_dzlev,mm_dzlay ! Condensible species parameters (mm_global_init) PROTECTED :: mm_nesp,mm_spcname,mm_xESPS ! Moments parameters (mm_aerosols_init / mm_clouds_init) PROTECTED :: mm_m0aer_s, mm_m3aer_s, mm_m0aer_f, mm_m3aer_f, mm_m0ccn, mm_m3ccn, mm_m3ice ! Moments parameters (derived, are updated with moments parameters) PROTECTED :: mm_rcs, mm_rcf, mm_drad, mm_drho ! Thresholds parameters PROTECTED :: mm_m0as_min, mm_m3as_min, mm_rcs_min, mm_m0af_min, mm_m3af_min, mm_rcf_min, mm_m0n_min, mm_m3cld_min LOGICAL, SAVE :: mm_debug = .false. !! Enable QnD debug mode (can be used for devel). LOGICAL, SAVE :: mm_log = .false. !! Enable log mode (for configuration only). LOGICAL, SAVE :: mm_w_haze_prod = .true. !! Enable/Disable haze production. LOGICAL, SAVE :: mm_w_haze_sed = .true. !! Enable/Disable haze sedimentation. LOGICAL, SAVE :: mm_w_haze_coag = .true. !! Activate haze coagulation. LOGICAL, SAVE :: mm_wsed_m0 = .false. !! Force all aerosols moments to fall at M0 settling velocity. LOGICAL, SAVE :: mm_wsed_m3 = .false. !! Force all aerosols moments to fall at M3 settling velocity. LOGICAL, SAVE :: mm_var_prod = .false. !! Time variation of production rate control flag. LOGICAL, SAVE :: mm_use_effg = .true. !! Enable/Disable effective G for computations. !> Enable/Disable __Fiadero__'s correction. !! !! This flag enables/disables the __Fiadero__ correction alogrithm for fractal mode settling velocity !! computation. !! !! @bug !! Currently, the Fiadero correction creates instatibilities on the vertical structure. It seems to be !! related to the coupling between the two moments. In order to reduce the instabilities, settling !! velocity of moments are forced to be the same, see [[mm_globals(module):mm_wsed_m0(variable)]] and !! [[mm_globals(module):mm_wsed_m3(variable)]]). LOGICAL, SAVE :: mm_no_fiadero_w = .false. !> Minimum ratio for __Fiadero__ correction. !! !! When [[mm_globals(module):mm_no_fiadero_w(variable)]] is disabled, this variable defines the minimum !! value of the moment's ratio between two adjacents vertical cells to be used within the correction. REAL(kind=mm_wp), SAVE :: mm_fiadero_min = 0.1_mm_wp !> Maximum ratio for __Fiadero__ correction. !! !! When [[mm_globals(module):mm_no_fiadero_w(variable)]] is disabled, this variable defines the maximum !! value of the moment's ratio between two adjacents vertical cells to be used within the correction. REAL(kind=mm_wp), SAVE :: mm_fiadero_max = 10._mm_wp LOGICAL, SAVE :: mm_w_clouds = .true. !! Enable/Disable clouds microphysics. LOGICAL, SAVE :: mm_w_cloud_sed = .true. !! Enable/Disable cloud sedimentation. LOGICAL, SAVE :: mm_w_cloud_nucond = .true. !! Activate cloud nucleation/condensation. INTEGER, PARAMETER :: mm_coag_no = 0 !! no mode interaction for coagulation (i.e. no coagulation at all). INTEGER, PARAMETER :: mm_coag_ss = 1 !! SS mode interaction for coagulation. INTEGER, PARAMETER :: mm_coag_sf = 2 !! SF mode interaction for coagulation. INTEGER, PARAMETER :: mm_coag_ff = 4 !! FF mode interaction for coagulation. !> Default interactions to activate (all by default). INTEGER, SAVE :: mm_coag_choice = mm_coag_ss+mm_coag_sf+mm_coag_ff !> Pi number. REAL(kind=mm_wp), PARAMETER :: mm_pi = 4._mm_wp*atan(1._mm_wp) !> Avogadro number. REAL(kind=mm_wp), PARAMETER :: mm_navo = 6.0221367e23_mm_wp !> Boltzmann constant (\(J.K^{-1}\)). REAL(kind=mm_wp), PARAMETER :: mm_kboltz = 1.3806488e-23_mm_wp !> Perfect gas constant (\(J.mol^{-1}.K^{-1}\)). REAL(kind=mm_wp), PARAMETER :: mm_rgas = mm_kboltz * mm_navo !> Desorption energy (\(J\)) (nucleation). REAL(kind=mm_wp), PARAMETER :: mm_fdes = 0.288e-19_mm_wp !> Surface diffusion energy (\(J\)) (nucleation). REAL(kind=mm_wp), PARAMETER :: mm_fdif = 0.288e-20_mm_wp !> Jump frequency (\(s^{-1}\)) (nucleation). REAL(kind=mm_wp), PARAMETER :: mm_nus = 1.e+13_mm_wp !> Approximated slip-flow correction coefficient. REAL(kind=mm_wp), PARAMETER :: mm_akn = 1.591_mm_wp !> Aerosols density (\(kg.m^{-3}\)). REAL(kind=mm_wp), SAVE :: mm_rhoaer = 1.e3_mm_wp !> Fractal dimension of fractal aerosols. REAL(kind=mm_wp), SAVE :: mm_df = 3._mm_wp !> Monomer radius (m). REAL(kind=mm_wp), SAVE :: mm_rm = 6.66e-8_mm_wp !> Spherical aerosols production pressure level (Pa). REAL(kind=mm_wp), SAVE :: mm_p_prod = 1._mm_wp !> Spherical aerosols equivalent radius production (m) REAL(kind=mm_wp), SAVE :: mm_rc_prod = 1.3101721857598102e-9_mm_wp !> Production rate of spherical aerosols (\(kg.m^{-2}.s^{-1}\)). REAL(kind=mm_wp), SAVE :: mm_tx_prod = 3.5e-13_mm_wp !> Aerosol production delta if time variations is enabled (fraction). REAL(kind=mm_wp), SAVE :: mm_d_prod = 0.25_mm_wp !> Aerosol production variations angular frequency if time variations is enabled (\(rad.s^{-1}\)). REAL(kind=mm_wp), SAVE :: mm_w_prod = 2.*mm_pi / (86400.*16.) !> Electric charging of aerosols (\(e^{-}.m^{-1}\)). REAL(kind=mm_wp), SAVE :: mm_ne = -15.e6_mm_wp !> Bulk to apparent radius conversion pre-factor (\(m^{X}\)). !! !! It is initialized using [[mm_globals(module):mm_rm(variable)]] in !! [[mm_globals(module):mm_global_init(interface)]] from the following equation: !! !! $$ r_{a} = r_{b}^{3/D_{f}}\times r_{m}^{\frac{D_{f}-3}{D_{f}}} $$ !! !! Where \(r_{a}\) is the apparent radius, \(r_{b}\) the bulk radius and !! \(rb2ra = r_{m}^{\frac{D_{f}-3}{D_{f}}}\) is the returned pre-factor REAL(kind=mm_wp), SAVE :: mm_rb2ra = 1._mm_wp ! Thresholds ! !> (min.) Total number of aerosols minimum threshold for the spherical mode. REAL(kind=mm_wp), SAVE :: mm_m0as_min = 1.e-10_mm_wp !> (min.) Total volume of aerosols minimum threshold for the spherical mode. REAL(kind=mm_wp), SAVE :: mm_m3as_min = 1.e-40_mm_wp !> Characteristic radius minimum threshold for the spherical mode. REAL(kind=mm_wp), SAVE :: mm_rcs_min = 1.e-9_mm_wp !> (min.) Total number of aerosols minimum threshold for the fractal mode. REAL(kind=mm_wp), SAVE :: mm_m0af_min = 1.e-10_mm_wp !> (min.) Total volume of aerosols minimum threshold for the fractal mode. REAL(kind=mm_wp), SAVE :: mm_m3af_min = 1.e-40_mm_wp !> Characteristic radius minimum threshold for the fractal mode. REAL(kind=mm_wp), SAVE :: mm_rcf_min = 1.e-9_mm_wp !> (min.) Total number of cloud drop minimum threshold. REAL(kind=mm_wp), SAVE :: mm_m0n_min = 1.e-10_mm_wp !> (min.) Total volume of cloud drop minimum threshold. REAL(kind=mm_wp), SAVE :: mm_m3cld_min = 1.e-40_mm_wp !> Minimum cloud drop radius REAL(kind=mm_wp), SAVE :: mm_drad_min = 1.e-9_mm_wp !> Maximum cloud drop radius REAL(kind=mm_wp), SAVE :: mm_drad_max = 1.e-3_mm_wp !> Characteristic radius threshold. REAL(kind=mm_wp), SAVE :: mm_rc_min = 1.e-200_mm_wp !> Name of condensible species. CHARACTER(len=30), DIMENSION(:), ALLOCATABLE, SAVE :: mm_spcname TYPE, PUBLIC :: mm_esp !! Cloud related chemical specie properties. !! !! This derived type is used in thermodynamic methods related to cloud microphysics. !! Most of its fields represent parameters of equations from \cite{reid1986}. CHARACTER(LEN=10) :: name !! Specie name. REAL(kind=mm_wp) :: mas !! Molecular weight (kg). REAL(kind=mm_wp) :: vol !! Molecular volume (\(m^{3}\)). REAL(kind=mm_wp) :: ray !! Molecular radius (m). REAL(kind=mm_wp) :: masmol !! Molar mass (\(kg.mol^{-1}\)). REAL(kind=mm_wp) :: rho !! density (liquid) (\(kg.m^{-3}\)). REAL(kind=mm_wp) :: tc !! Critical temperature (K). REAL(kind=mm_wp) :: pc !! Critical pressure (Bar). REAL(kind=mm_wp) :: tb !! Boiling point temperature (K). REAL(kind=mm_wp) :: w !! Acentric factor (--). REAL(kind=mm_wp) :: a_sat !! Saturation equation A coefficient. REAL(kind=mm_wp) :: b_sat !! Saturation equation B coefficient. REAL(kind=mm_wp) :: c_sat !! saturation equation C coefficient. REAL(kind=mm_wp) :: d_sat !! Saturation equation D coefficient. REAL(kind=mm_wp) :: mteta !! Wettability. REAL(kind=mm_wp) :: tx_prod !! Production rate. REAL(kind=mm_wp) :: fmol2fmas !! molar fraction to mass fraction coefficient. ! = masmol(X)/masmol(AIR) END TYPE mm_esp !> Planet radius (m). REAL(kind=mm_wp), SAVE :: mm_rpla = 2575000._mm_wp !> Planet acceleration due to gravity constant (ground) (\(m.s^{-2}\)). REAL(kind=mm_wp), SAVE :: mm_g0 = 1.35_mm_wp !> Air molecules mean radius (m). REAL(kind=mm_wp), SAVE :: mm_air_rad = 1.75e-10_mm_wp !> Air molecules molar mass (\(kg.mol^{-1}\)). REAL(kind=mm_wp), SAVE :: mm_air_mmol = 28e-3_mm_wp !> Microphysic time step (s). REAL(kind=mm_wp), SAVE :: mm_dt = 5529.6_mm_wp !> Model current time tracer (s). REAL(kind=mm_wp), SAVE :: mm_ct = 0.0 !> Total number of clouds condensible species. INTEGER, SAVE :: mm_nesp = -1 !> Clouds chemical species properties. TYPE(mm_esp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_xESPS !------------------------ ! Vertical structure part !------------------------ !> Number of vertical layers. INTEGER, SAVE :: mm_nla = -1 !> Number of vertical levels. INTEGER, SAVE :: mm_nle = -1 !> Altitude layers (m). REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_zlay !> Altitude levels (m). REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_zlev !> Pressure layers (Pa). REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_play !> Pressure levels (Pa). REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_plev !> Temperature vertical profile (K). REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_temp !> Air density vertical profile (\(kg.m^{-3}\)). REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_rhoair !> Temperature vertical profil at interfaces (K). REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_btemp !> Atmospheric levels thickness (m). !! !! Atmospheric thickness between two adjacent levels (\(m\)) from the !! __TOP__ to the __GROUND__. !! @note __mm_dzlev__ is defined on the total number of layers and actually !! corresponds to the thickness of a given layer. REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_dzlev !> Atmospheric layers "thickness" (m). !! !! Atmospheric thickness between the center of two adjacent layers (\(m\)) !! from the __TOP__ to the __GROUND__. !! @note !! __mm_dzlay__ is defined on the total number of layers. The last !! value of __mm_dzlay__ is set to twice the altitude of the ground layer. !! @note This value corresponds to the thickness between the center of the !! __GROUND__ layer and below the surface. It is arbitrary and not used. REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_dzlay !> Spherical mode \(0^{th}\) order moment (\(m^{-3}\)). REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m0aer_s !> Spherical mode \(3^{rd}\) order moment (\(m^{3}.m^{-3}\)). REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m3aer_s !> Fractal mode \(0^{th}\) order moment (\(m^{-3}\)). REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m0aer_f !> Fractal mode \(3^{rd}\) order moment (\(m^{3}.m^{-3}\)). REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m3aer_f !> CCN \(0^{th}\) order moment (\(m^{-3}\)). REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m0ccn !> CCN \(3^{rd}\) order moment (\(m^{3}.m^{-3}\)). REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m3ccn !> Ice components 3rd order moments (\(m^{3}.m^{-3}\)). !! !! It is a 2D array with the vertical layers in first dimension, and the number of ice !! components in the second. !! @note !! Both [[mm_globals(module):mm_m3ice(variable)]] and [[mm_globals(module):mm_gazs(variable)]] !! share the same indexing (related to species order). REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: mm_m3ice !> Condensible species molar fraction (\(mol.mol^{-1}\)). !! !! It is a 2D array with the vertical layers in first dimension, and !! the number of condensible species in the second. !! @note !! Both [[mm_globals(module):mm_m3ice(variable)]] and [[mm_globals(module):mm_gazs(variable)]] !! share the same indexing (related to species order). REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: mm_gazs !> Spherical mode characteristic radius (m). REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_rcs !> Fractal mode characteristic radius (m). REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_rcf !> Mean Drop radius (m). REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_drad !> Mean Drop density (\(kg.m^{-3}\)). REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_drho !> Aerosols precipitations (m). !! !! Aerosols precipitations take into account both spherical and fractal modes. !! It is updated in [[mm_haze(module):mm_haze_microphysics(subroutine)]]. REAL(kind=mm_wp), SAVE :: mm_aer_prec = 0._mm_wp !> CCN precipitations (m). !! It is updated in [[mm_clouds(module):mm_cloud_microphysics(subroutine)]]. REAL(kind=mm_wp), SAVE :: mm_ccn_prec = 0._mm_wp !> Spherical mode \(M_{0}\) settling velocity (\(m.s^{-1}\)). !! !! It is a vector with the vertical layers that contains the settling velocity for !! the \(0^{th}\) order moment of the spherical mode. !! It is updated in [[mm_haze(module):mm_haze_sedimentation(subroutine)]]. !! @note !! This variable is always negative. REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m0as_vsed !> Spherical mode \(M_{3}\) settling velocity (\(m.s^{-1}\)). !! !! It is a vector with the vertical layers that contains the settling velocity for the !! \(3^{rd}\) order moment of the spherical mode. !! It is updated in [[mm_haze(module):mm_haze_sedimentation(subroutine)]]. !! @note !! This variable is always negative. REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m3as_vsed !> Fractal mode \(M_{0}\) settling velocity (\(m.s^{-1}\)). !! !! It is a vector with the vertical layers that contains the settling velocity for the !! \(0^{th}\) order moment of the fractal mode. !! It is updated in [[mm_haze(module):mm_haze_sedimentation(subroutine)]]. !! @note !! This variable is always negative. REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m0af_vsed !> Fractal mode \(M_{3}\) settling velocity (\(m.s^{-1}\)). !! !! It is a vector with the vertical layers that contains the settling velocity for the !! \(3^{rd}\) order moment of the fractal mode. !! It is updated in [[mm_haze(module):mm_haze_sedimentation(subroutine)]]. !! @note !! This variable is always negative. REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m3af_vsed !> CCN settling velocity (\(m.s^{-1}\)). !! !! It is a vector with the vertical layers that contains the !! settling velocity for CCN (and ices). !! It is updated in [[mm_clouds(module):mm_cloud_microphysics(subroutine)]]. REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_ccn_vsed !> Spherical aerosol mass fluxes (\(kg.m^{-2}.s^{-1}\)). !! !! It is a vector with the vertical layers that contains the mass fluxes for spherical aerosols. !! It is updated in [[mm_haze(module):mm_haze_sedimentation(subroutine)]]. !! @note !! This variable is always negative. REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_aer_s_flux !> Fractal aerosol mass fluxes (\(kg.m^{-2}.s^{-1}\)). !! !! It is a vector with the vertical layers that contains the mass fluxes for fractal aerosols !! It is updated in [[mm_haze(module):mm_haze_sedimentation(subroutine)]]. !! @note !! This variable is always negative. REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_aer_f_flux !> CCN mass fluxes (\(kg.m^{-2}.s^{-1}\)). !! !! It is a vector with the vertical layers that contains the !! mass fluxes for CCN. !! It is updated in [[mm_clouds(module):mm_cloud_microphysics(subroutine)]]. REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_ccn_flux !> Ice components precipitations (m). !! !! It is a vector of [[mm_globals(module):mm_nesp(variable)]] values which share the same indexing !! than [[mm_globals(module):mm_m3ice(variable)]] and [[mm_globals(module):mm_gazs(variable)]]. !! It is updated in [[mm_clouds(module):mm_cloud_microphysics(subroutine)]]. !! @note !! This variable is always negative. REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_ice_prec !> Ice components sedimentation fluxes (\(kg.m^{-2}.s-1\)). !! !! It is a 2D-array with the vertical layers in first dimension and the number of ice components !! in the second. It is updated in [[mm_clouds(module):mm_cloud_microphysics(subroutine)]]. !! @note !! This variable is always negative. REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: mm_ice_fluxes !> Condensible species saturation ratio (--). !! !! It is a 2D-array with the vertical layers in first dimension and the number of condensible !! species in the second. !! It is updated in [[mm_clouds(module):mm_cloud_microphysics(subroutine)]]. REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: mm_gazs_sat !> [[mm_globals(module):mm_global_init(interface)]] initialization control flag. LOGICAL, PUBLIC, SAVE :: mm_ini = .false. !> [[mm_globals(module):mm_column_init(function)]] initialization control flag. LOGICAL, PUBLIC, SAVE :: mm_ini_col = .false. !> [[mm_globals(module):mm_aerosols_init(function)]] initialization control flag. LOGICAL, PUBLIC, SAVE :: mm_ini_aer = .false. !> [[mm_globals(module):mm_clouds_init(function)]] initialization control flag. LOGICAL, PUBLIC, SAVE :: mm_ini_cld = .false. !> Interface to cloud properties methods. !! !! The method computes clouds properties (mean drop radius and denstity) from their afferent !! moments. It is overloaded to compute properties at a single level or over all the vertical !! atmospheric structure. INTERFACE mm_cloud_properties MODULE PROCEDURE cldprop_sc,cldprop_ve END INTERFACE mm_cloud_properties !> Interface to global initialization. !! !! The method performs the global initialization of the model. !! @warning !! If OpenMP is activated, this subroutine must be called in an $OMP SINGLE statement as it !! initializes global variable that are not thread private. !! !! !$OMP SINGLE !! call mm_global_init(...) !! !$OMP END SINGLE INTERFACE mm_global_init MODULE PROCEDURE mm_global_init_0,mm_global_init_1 END INTERFACE mm_global_init !> Check an option from the configuration system. !! !! The method checks for an option in the configuration system and optionally !! set a default value if the option is not found. This is an overloaded method !! that can take in input either a floating point, integer, logical or string !! option value. INTERFACE mm_check_opt MODULE PROCEDURE check_r1,check_i1,check_l1,check_s1 END INTERFACE mm_check_opt ! --- OPENMP --------------- ! All variables related to column computations should be private to each thread ! !$OMP THREADPRIVATE(mm_ini_col,mm_ini_aer,mm_ini_cld) !$OMP THREADPRIVATE(mm_zlay,mm_zlev,mm_play,mm_plev,mm_temp,mm_rhoair,mm_btemp,mm_dzlev,mm_dzlay) !$OMP THREADPRIVATE(mm_m0aer_s,mm_m3aer_s,mm_m0aer_f,mm_m3aer_f) !$OMP THREADPRIVATE(mm_m0ccn,mm_m3ccn,mm_m3ice,mm_gazs) !$OMP THREADPRIVATE(mm_rcs,mm_rcf,mm_drad,mm_drho) !$OMP THREADPRIVATE(mm_m0as_vsed,mm_m3as_vsed,mm_m0af_vsed,mm_m3af_vsed) !$OMP THREADPRIVATE(mm_aer_s_flux,mm_aer_f_flux,mm_ccn_vsed,mm_ccn_flux,mm_ice_prec,mm_ice_fluxes,mm_gazs_sat) !$OMP THREADPRIVATE(mm_m0as_min,mm_m3as_min,mm_rcs_min,mm_m0af_min,mm_m3af_min,mm_rcf_min,mm_m0n_min,mm_m3cld_min) !$OMP THREADPRIVATE(mm_nla,mm_nle) ! -------------------------- CONTAINS FUNCTION mm_global_init_0(dt,df,rm,rho_aer,p_prod,tx_prod,rc_prod,rplanet,g0, & air_rad,air_mmol,coag_interactions,clouds,spcfile, & w_haze_prod,w_haze_sed,w_haze_coag,w_cloud_nucond, & w_cloud_sed,force_wsed_to_m0,force_wsed_to_m3, & no_fiadero,fiadero_min,fiadero_max, & m0as_min,rcs_min,m0af_min,rcf_min,m0n_min,debug) RESULT(err) !! Initialize global parameters of the model. !! !! The function initializes all the global parameters of the model from direct input. !! Boolean, Fiadero and thresholds parameters are optional as they are rather testing parameters. !! Their default values are suitable for production runs. !! @note !! If the method fails to initialize parameters (i.e. returned error is not 0). Then the model !! should probably be aborted as the global variables of the model will not be correctly setup. !! @warning !! If OpenMP is activated, this subroutine must be called in an $OMP SINGLE statement as it !! initializes (only) global variables that are not thread private. !! !! !$OMP SINGLE !! call mm_global_init_0(...) !! !$OMP END SINGLE REAL(kind=mm_wp), INTENT(in) :: dt !! Microphysics timestep in seconds. REAL(kind=mm_wp), INTENT(in) :: df !! Fractal dimension of fractal aerosol. REAL(kind=mm_wp), INTENT(in) :: rm !! Monomer radius in meter. REAL(kind=mm_wp), INTENT(in) :: rho_aer !! Aerosol density in \(kg.m^{-3}\). REAL(kind=mm_wp), INTENT(in) :: p_prod !! Aerosol production pressure level in Pa. REAL(kind=mm_wp), INTENT(in) :: tx_prod !! Spherical aerosol mode production rate in \(kg.m^{-2}.s^{-1}\). REAL(kind=mm_wp), INTENT(in) :: rc_prod !! Spherical mode characteristic radius for production in meter. REAL(kind=mm_wp), INTENT(in) :: rplanet !! Planet radius in meter REAL(kind=mm_wp), INTENT(in) :: g0 !! Planet gravity acceleration at ground level in \(m.s^{-2}\). REAL(kind=mm_wp), INTENT(in) :: air_rad !! Air molecules mean radius in meter. REAL(kind=mm_wp), INTENT(in) :: air_mmol !! Air molecules mean molar mass in \(kg.mol^{-1}\). INTEGER, INTENT(in) :: coag_interactions !! Coagulation interactions process control flag. LOGICAL, INTENT(in) :: clouds !! Clouds microphysics control flag. CHARACTER(len=*), INTENT(in) :: spcfile !! Clouds microphysics condensible species properties file. REAL(kind=mm_wp), INTENT(in), OPTIONAL :: fiadero_max !! Maximum moment ratio threshold for Fiadero correction (default: 10.) . REAL(kind=mm_wp), INTENT(in), OPTIONAL :: fiadero_min !! Minimum moment ratio threshold for Fiadero correction (default: 0.1). LOGICAL, INTENT(in), OPTIONAL :: w_haze_prod !! Haze microphysics production process control flag (default: T). LOGICAL, INTENT(in), OPTIONAL :: w_haze_sed !! Haze microphysics sedimentation process control flag (default: T). LOGICAL, INTENT(in), OPTIONAL :: w_haze_coag !! Haze microphysics coagulation process control flag (default: T). LOGICAL, INTENT(in), OPTIONAL :: w_cloud_sed !! Cloud microphysics nucleation/conensation process control flag (default: __clouds__ value). LOGICAL, INTENT(in), OPTIONAL :: w_cloud_nucond !! Cloud microphysics production process control flag (default: __clouds__ value). LOGICAL, INTENT(in), OPTIONAL :: no_fiadero !! Disable Fiadero correction for haze sedimentation process (default: F). LOGICAL, INTENT(in), OPTIONAL :: force_wsed_to_m0 !! force __all__ aerosols moments to fall at M0 settling velocity (default: T). LOGICAL, INTENT(in), OPTIONAL :: force_wsed_to_m3 !! Force __all__ aerosols moments to fall at M3 settling velocity (default: F). REAL(kind=mm_wp), INTENT(in), OPTIONAL :: m0as_min !! Minimum threshold for M0 of the spherical mode (default: 1e-10). REAL(kind=mm_wp), INTENT(in), OPTIONAL :: rcs_min !! Minimum threshold for the characteristic radius of the spherical mode in meter (default: 1e-9). REAL(kind=mm_wp), INTENT(in), OPTIONAL :: m0af_min !! Minimum threshold for M0 of the factal mode (default: 1e-10). REAL(kind=mm_wp), INTENT(in), OPTIONAL :: rcf_min !! Minimum threshold for the characteristic radius of the fractal mode in meter (default: _monomer radius_). REAL(kind=mm_wp), INTENT(in), OPTIONAL :: m0n_min !! Minimum threshold for M0 of cloud drop (default: 1e-10). LOGICAL, INTENT(in), OPTIONAL :: debug !! Debug mode control flag (may print lot of stuff if enabled) TYPE(error) :: err !! Error status of the function. INTEGER :: i TYPE(cfgparser) :: cp CHARACTER(len=st_slen), DIMENSION(:), ALLOCATABLE :: species REAL(kind=mm_wp) :: zfiamin,zfiamax LOGICAL :: zwhp,zwhs,zwhc,zwcs,zwcn,znofia, & zwstom0,zwstom3 zwhp = .true. ; zwhs = .true. ; zwhc = .true. zwcs = clouds ; zwcn = clouds znofia = .false. ; zfiamin = 0.1_mm_wp ; zfiamax = 10._mm_wp zwstom0 = .true. ; zwstom3 = .false. err = noerror IF (mm_ini) THEN err = error("mm_global_init: YAMMS global initialization already performed !",-1) RETURN ENDIF ! Store options values in global variables... mm_df = df mm_rm = rm mm_rb2ra = mm_rm**((mm_df-3._mm_wp)/mm_df) ! conversion factor for bulk -> fractal radius mm_rhoaer = rho_aer mm_p_prod = p_prod mm_tx_prod = tx_prod mm_rc_prod = rc_prod mm_rpla = rplanet mm_g0 = g0 mm_dt = dt mm_air_rad = air_rad mm_air_mmol = air_mmol mm_coag_choice = coag_interactions ! check coagulation interactions choice IF (mm_coag_choice < 0 .OR. mm_coag_choice > 7) THEN err = error("mm_global_init: Invalid choice for coagulation interactions activation",-1) RETURN ENDIF ! force fractal radius minimum threshold to monomer radius ^^ mm_rcf_min = mm_rm mm_w_clouds = clouds ! Check clouds microphysics species file ! (only if clouds is activated) IF (mm_w_clouds) THEN IF (LEN_TRIM(spcfile) == 0) THEN err = error("mm_global_init: No species properties file given",-1) RETURN ENDIF ! Reads species properties configuration file err = cfg_read_config(cp,TRIM(spcfile)) ; IF (err /= 0) RETURN err = cfg_get_value(cp,"used_species",species) IF (err /= 0) THEN err = error("mm_global_init: cannot retrieve 'used_species' values",-1) RETURN ENDIF ! Now attempts to find species properties !!! mm_nesp = SIZE(species) ALLOCATE(mm_spcname(mm_nesp),mm_xESPS(mm_nesp)) DO i=1,mm_nesp mm_spcname(i) = TRIM(species(i)) IF(.NOT.cfg_has_section(cp,TRIM(mm_spcname(i)))) THEN err = error("mm_global_init: Cannot find "//TRIM(mm_spcname(i))//" properties",-1) RETURN ELSE err = read_esp(cp,TRIM(mm_spcname(i)),mm_xESPS(i)) ! compute conversion factor: mol.mol-1 => kg.kg-1 mm_xESPS(i)%fmol2fmas = mm_xESPS(i)%masmol / mm_air_mmol IF (err/=0) THEN err = error("mm_global_init: "//TRIM(mm_spcname(i))//": "//TRIM(err%msg),-1) RETURN ENDIF ENDIF ENDDO ENDIF ! optional flags ! debug mode IF (PRESENT(debug)) THEN mm_debug = debug ELSE mm_debug = .false. call printw("mm_debug",to_string(mm_debug)) ENDIF ! haze control flags IF (PRESENT(w_haze_prod)) THEN mm_w_haze_prod = w_haze_prod ELSE mm_w_haze_prod = zwhp call printw("mm_haze_production",to_string(mm_w_haze_prod)) ENDIF IF (PRESENT(w_haze_sed)) THEN mm_w_haze_sed = w_haze_sed ELSE mm_w_haze_sed = zwhs call printw("mm_haze_sedimentation",to_string(mm_w_haze_sed)) ENDIF IF (PRESENT(w_haze_coag)) THEN mm_w_haze_coag = w_haze_coag ELSE mm_w_haze_coag = zwhc call printw("mm_haze_coagulation",to_string(mm_w_haze_coag)) ENDIF IF (PRESENT(force_wsed_to_m0)) THEN mm_wsed_m0 = force_wsed_to_m0 ELSE mm_wsed_m0 = zwstom0 call printw("mm_wsed_m0",to_string(mm_wsed_m0)) ENDIF IF (PRESENT(force_wsed_to_m3)) THEN mm_wsed_m3 = force_wsed_to_m3 ELSE mm_wsed_m3 = zwstom3 call printw("mm_wsed_m3",to_string(mm_wsed_m3)) ENDIF IF (PRESENT(no_fiadero)) THEN mm_no_fiadero_w = no_fiadero ELSE mm_no_fiadero_w = znofia call printw("mm_no_fiadero",to_string(mm_no_fiadero_w)) ENDIF IF (PRESENT(fiadero_min)) THEN mm_fiadero_min = fiadero_min ELSE mm_fiadero_min = zfiamin call printw("mm_fiadero_min",to_string(mm_fiadero_min)) ENDIF IF (PRESENT(fiadero_max)) THEN mm_fiadero_max = fiadero_max ELSE mm_fiadero_max = zfiamax call printw("mm_fiadero_max",to_string(mm_fiadero_max)) ENDIF ! moments threshold flags IF (PRESENT(m0as_min)) THEN mm_m0as_min = MAX(0._mm_wp,m0as_min) ELSE call printw("mm_m0as_min",to_string(mm_m0as_min)) ENDIF IF (PRESENT(rcs_min)) THEN mm_rcs_min = MAX(1.e-9_mm_wp,rcs_min) ELSE call printw("mm_rcs_min",to_string(mm_rcs_min)) ENDIF IF (PRESENT(m0af_min)) THEN mm_m0af_min = MAX(0._mm_wp,m0af_min) ELSE call printw("mm_m0af_min",to_string(mm_m0af_min)) ENDIF IF (PRESENT(rcf_min)) THEN mm_rcf_min = MAX(rcf_min,mm_rm) ELSE mm_rcf_min = mm_rm call printw("mm_rcf_min",to_string(mm_rcf_min)) ENDIF IF (PRESENT(m0n_min)) THEN mm_m0n_min = MAX(0._mm_wp,m0n_min) ELSE call printw("mm_m0n_min",to_string(mm_m0n_min)) ENDIF ! compute m3 thresholds from user-defined thresholds. mm_m3as_min = mm_m0as_min*mm_alpha_s(3._mm_wp) * mm_rcs_min**3._mm_wp mm_m3af_min = mm_m0af_min*mm_alpha_f(3._mm_wp) * mm_rcf_min**3._mm_wp mm_m3cld_min = mm_m0n_min * (4._mm_wp * mm_pi / 3._mm_wp) * mm_drad_min**3._mm_wp ! clouds control flags IF (mm_w_clouds) THEN IF (PRESENT(w_cloud_sed)) THEN mm_w_cloud_sed = w_cloud_sed ELSE mm_w_cloud_sed = zwcs call printw("mm_cloud_sed",to_string(mm_w_cloud_sed)) ENDIF IF (PRESENT(w_cloud_nucond)) THEN mm_w_cloud_nucond = w_cloud_nucond ELSE mm_w_cloud_nucond = zwcs call printw("mm_cloud_nucond",to_string(mm_w_cloud_nucond)) ENDIF ENDIF ! check w sed flags err = noerror ! special check for settling velocity IF (mm_wsed_m0 .AND. mm_wsed_m3) THEN err = error("'wsed_m0' and 'wsed_m3' options are mutually exclusive",-1) ENDIF mm_ini = err == noerror CONTAINS SUBROUTINE printw(string,value) !! Print a warning message. CHARACTER(len=*), INTENT(in) :: string !! Name of the option. CHARACTER(len=*), INTENT(in) :: value !! (string) Value of the option. IF (mm_log) & WRITE(*,'(a,a,a)') "warning: Parameter "//string//"not given... Using default value: "//value END SUBROUTINE printw END FUNCTION mm_global_init_0 FUNCTION mm_global_init_1(cfg) RESULT(err) !! Set global configuration from a configuration file. !! !! See [[mm_globals(module):mm_global_init_0(function)]]. TYPE(cfgparser), INTENT(in) :: cfg !! Configuration file path. TYPE(error) :: err !! Error status of the function. INTEGER :: i TYPE(cfgparser) :: spccfg CHARACTER(len=st_slen) :: spcpath CHARACTER(len=st_slen), DIMENSION(:), ALLOCATABLE :: species REAL(kind=mm_wp) :: zfiamin,zfiamax LOGICAL :: zwhp,zwhs,zwhc,zwcs,zwcn,znofia, & zwstom0,zwstom3 err = noerror IF (mm_ini) THEN err = error("mm_global_init: YAMMS global initialization already performed !",-1) RETURN ENDIF ! MP2M mandatory parameters err = mm_check_opt(cfg_get_value(cfg,"df",mm_df),mm_df,wlog=mm_log) IF (err/=0) RETURN err = mm_check_opt(cfg_get_value(cfg,"rm",mm_rm),mm_rm,wlog=mm_log) IF (err/=0) RETURN err = mm_check_opt(cfg_get_value(cfg,"rho_aer",mm_rhoaer),mm_rhoaer,wlog=mm_log) IF (err/=0) RETURN err = mm_check_opt(cfg_get_value(cfg,"p_prod",mm_p_prod),mm_p_prod,wlog=mm_log) IF (err/=0) RETURN err = mm_check_opt(cfg_get_value(cfg,"tx_prod",mm_tx_prod),mm_tx_prod,wlog=mm_log) IF (err/=0) RETURN err = mm_check_opt(cfg_get_value(cfg,"rc_prod",mm_rc_prod),mm_rc_prod,wlog=mm_log) IF (err/=0) RETURN err = mm_check_opt(cfg_get_value(cfg,"planet_radius",mm_rpla),mm_rpla,wlog=mm_log) IF (err/=0) RETURN err = mm_check_opt(cfg_get_value(cfg,"g0",mm_g0),mm_g0,wlog=mm_log) IF (err/=0) RETURN err = mm_check_opt(cfg_get_value(cfg,"timestep",mm_dt),mm_dt,wlog=mm_log) IF (err/=0) RETURN err = mm_check_opt(cfg_get_value(cfg,"air_radius",mm_air_rad),mm_air_rad,wlog=mm_log) IF (err/=0) RETURN err = mm_check_opt(cfg_get_value(cfg,"air_molarmass",mm_air_mmol),mm_air_mmol,wlog=mm_log) IF (err/=0) RETURN err = mm_check_opt(cfg_get_value(cfg,"haze_coag_interactions",mm_coag_choice),mm_coag_choice,wlog=mm_log) IF (err/=0) RETURN err = mm_check_opt(cfg_get_value(cfg,"clouds_microphysics",mm_w_clouds),mm_w_clouds,wlog=mm_log) IF (err/=0) RETURN ! computes the conversion factor for bulk -> fractal radius mm_rb2ra = mm_rm**((mm_df-3._mm_wp)/mm_df) ! Check coagulation interactions choice IF (mm_coag_choice < 0 .OR. mm_coag_choice > 7) THEN err = error("mm_global_init: Invalid choice for coagulation interactions activation",-1) RETURN ENDIF ! Check clouds microphysics input ! it is read only if clouds is activated. We must to check if it is self-consistent... IF (mm_w_clouds) THEN ! Gets species property file path err = cfg_get_value(cfg,'specie_cfg',spcpath) ; IF (err /= 0) RETURN ! Reads species properties configuration file err = cfg_read_config(spccfg,trim(spcpath)) ; IF (err /= 0) RETURN err = cfg_get_value(spccfg,"used_species",species) IF (err /= 0) THEN err = error("mm_global_init: cannot retrieve 'used_species' values",-1) RETURN ENDIF ! Now attempts to find specides properties !!! mm_nesp = SIZE(species) ALLOCATE(mm_spcname(mm_nesp),mm_xESPS(mm_nesp)) !mm_spcname(1:mm_nesp) = species(:) DO i=1,mm_nesp mm_spcname(i) = TRIM(species(i)) IF (.NOT.cfg_has_section(spccfg,TRIM(mm_spcname(i)))) THEN err = error("mm_global_init: Cannot find "//TRIM(mm_spcname(i))//" properties",-1) RETURN ELSE err = read_esp(spccfg,TRIM(mm_spcname(i)),mm_xESPS(i)) ! compute conversion factor: mol.mol-1 => kg.kg-1 mm_xESPS(i)%fmol2fmas = mm_xESPS(i)%masmol / mm_air_mmol IF (err/=0) THEN err = error(TRIM(mm_spcname(i))//": "//TRIM(err%msg),-2) RETURN ENDIF ENDIF ENDDO ENDIF zwhp = .true. ; zwhs = .true. ; zwhc = .true. zwcs = mm_w_clouds ; zwcn = mm_w_clouds znofia = .false. ; zfiamin = 0.1_mm_wp ; zfiamax = 10._mm_wp zwstom0 = .true. ; zwstom3 = .false. ! MP2M Optional parameters err = mm_check_opt(cfg_get_value(cfg,"debug",mm_debug),mm_debug,.false.,wlog=mm_log) err = mm_check_opt(cfg_get_value(cfg,"haze_production",mm_w_haze_prod),mm_w_haze_prod,zwhp,wlog=mm_log) err = mm_check_opt(cfg_get_value(cfg,"haze_sedimentation",mm_w_haze_sed),mm_w_haze_sed,zwhs,wlog=mm_log) err = mm_check_opt(cfg_get_value(cfg,"haze_coagulation",mm_w_haze_coag),mm_w_haze_coag,zwhc,wlog=mm_log) err = mm_check_opt(cfg_get_value(cfg,"clouds_sedimentation",mm_w_cloud_sed),mm_w_cloud_sed,zwcs,wlog=mm_log) err = mm_check_opt(cfg_get_value(cfg,"clouds_nucl_cond",mm_w_cloud_nucond),mm_w_cloud_nucond,zwcn,wlog=mm_log) err = mm_check_opt(cfg_get_value(cfg,"wsed_m0",mm_wsed_m0),mm_wsed_m0,zwstom0,wlog=mm_log) err = mm_check_opt(cfg_get_value(cfg,"wsed_m3",mm_wsed_m3),mm_wsed_m3,zwstom3,wlog=mm_log) err = mm_check_opt(cfg_get_value(cfg,"no_fiadero",mm_no_fiadero_w),mm_no_fiadero_w,znofia,wlog=mm_log) err = mm_check_opt(cfg_get_value(cfg,"fiadero_min_ratio",mm_fiadero_min),mm_fiadero_min,zfiamin,wlog=mm_log) err = mm_check_opt(cfg_get_value(cfg,"fiadero_max_ratio",mm_fiadero_max),mm_fiadero_max,zfiamax,wlog=mm_log) err = mm_check_opt(cfg_get_value(cfg,"m0as_min",mm_m0as_min),mm_m0as_min,1e-10_mm_wp,wlog=mm_log) err = mm_check_opt(cfg_get_value(cfg,"rcs_min",mm_rcs_min),mm_rcs_min,1e-9_mm_wp,wlog=mm_log) err = mm_check_opt(cfg_get_value(cfg,"m0af_min",mm_m0af_min),mm_m0af_min,1e-10_mm_wp,wlog=mm_log) err = mm_check_opt(cfg_get_value(cfg,"rcf_min",mm_rcf_min),mm_rcf_min,mm_rm,wlog=mm_log) err = mm_check_opt(cfg_get_value(cfg,"m0n_min",mm_m0n_min),mm_m0n_min,1e-10_mm_wp,wlog=mm_log) ! force fractal radius minimum threshold to monomer radius ^^ mm_rcf_min = MAX(mm_rm,mm_rcf_min) ! compute m3 thresholds from user-defined thresholds. mm_m3as_min = mm_m0as_min*mm_alpha_s(3._mm_wp) * mm_rcs_min**3._mm_wp mm_m3af_min = mm_m0af_min*mm_alpha_f(3._mm_wp) * mm_rcf_min**3._mm_wp mm_m3cld_min = mm_m0n_min * (4._mm_wp * mm_pi / 3._mm_wp) * mm_drad_min**3._mm_wp err = noerror ! special check for settling velocity IF (mm_wsed_m0 .AND. mm_wsed_m3) THEN err = error("'wsed_m0' and 'wsed_m3' options are mutually exclusive",-1) ENDIF mm_ini = err == noerror END FUNCTION mm_global_init_1 FUNCTION mm_column_init(plev,zlev,play,zlay,temp) RESULT(err) !! Initialize vertical atmospheric fields. !! !! This subroutine initializes vertical fields needed by the microphysics: !! !! 1. Save reversed input field into "local" array !! 2. Compute thicknesses layers and levels !! 3. Interpolate temperature at levels !! !! The method should be called whenever the vertical structure of the atmosphere changes. !! !! @attention !! All the input vectors should be defined from __GROUND__ to __TOP__ of the atmosphere, !! otherwise nasty things will occur in computations. REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: plev !! Pressure levels (Pa). REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: zlev !! Altitude levels (m). REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: play !! Pressure layers (Pa). REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: zlay !! Altitude at the center of each layer (m). REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: temp !! Temperature at the center of each layer (K). TYPE(error) :: err !! Error status of the function. INTEGER :: i mm_ini_col = .false. err = noerror IF (.NOT.mm_ini) THEN err = error("mm_column_init: Global initialization not done yet",-1) RETURN ENDIF IF (mm_nla < 0) THEN mm_nla = SIZE(play) ELSE IF (mm_nla /= SIZE(play)) THEN err = error("mm_column_init: mm_nla cannot be modified dynamically within the run",-1) RETURN ENDIF ENDIF IF (mm_nle < 0) THEN mm_nle = SIZE(plev) ELSE IF (mm_nle /= SIZE(plev)) THEN err = error("mm_column_init: mm_nle cannot be modified dynamically within the run",-1) RETURN ENDIF ENDIF ! should be trashed soon or later IF (mm_nla+1 /= mm_nle) THEN err = error("mm_column_init: Inconsistent number of layers/levels",-1) RETURN ENDIF ! Allocates if required IF (.NOT.ALLOCATED(mm_plev)) ALLOCATE(mm_plev(mm_nle)) IF (.NOT.ALLOCATED(mm_zlev)) ALLOCATE(mm_zlev(mm_nle)) IF (.NOT.ALLOCATED(mm_play)) ALLOCATE(mm_play(mm_nla)) IF (.NOT.ALLOCATED(mm_zlay)) ALLOCATE(mm_zlay(mm_nla)) IF (.NOT.ALLOCATED(mm_temp)) ALLOCATE(mm_temp(mm_nla)) IF (.NOT.ALLOCATED(mm_btemp)) ALLOCATE(mm_btemp(mm_nle)) IF (.NOT.ALLOCATED(mm_dzlev)) ALLOCATE(mm_dzlev(mm_nla)) IF (.NOT.ALLOCATED(mm_dzlay)) ALLOCATE(mm_dzlay(mm_nla)) IF (.NOT.ALLOCATED(mm_rhoair)) ALLOCATE(mm_rhoair(mm_nla)) ! Saves reversed input vectors mm_zlay = zlay(mm_nla:1:-1) ; mm_zlev = zlev(mm_nle:1:-1) mm_play = play(mm_nla:1:-1) ; mm_plev = plev(mm_nle:1:-1) mm_temp = temp(mm_nla:1:-1) ! Computes others vectors mm_dzlay(1:mm_nla-1) = mm_zlay(1:mm_nla-1)-mm_zlay(2:mm_nla) mm_dzlay(mm_nla) = mm_dzlay(mm_nla-1) ! actually arbitrary mm_dzlev(1:mm_nla) = mm_zlev(1:mm_nle-1)-mm_zlev(2:mm_nle) mm_btemp(2:mm_nla) = (mm_temp(1:mm_nla-1)+mm_temp(2:mm_nla))/2._mm_wp mm_btemp(1) = mm_temp(1) mm_btemp(mm_nle) = mm_temp(mm_nla)+(mm_temp(mm_nla)-mm_temp(mm_nla-1))/2._mm_wp ! Hydrostatic equilibrium mm_rhoair(1:mm_nla) = (mm_plev(2:mm_nle)-mm_plev(1:mm_nla)) / & (mm_effg(mm_zlay)*mm_dzlev) mm_ini_col = .true. ! write out profiles (only if BOTH debug and log are enabled). IF (mm_log.AND.mm_debug) THEN WRITE(*,'(a)') '# TEMP PLAY ZLAY DZLAY RHOAIR' DO i=1,mm_nla WRITE(*,'(5(ES15.7,2X))') mm_temp(i),mm_play(i),mm_zlay(i),mm_dzlay(i), mm_rhoair(i) ENDDO WRITE(*,'(a)') '# TEMP PLEV ZLEV DZLEV' DO i=1,mm_nle IF (i /= mm_nle) THEN WRITE(*,'(4(ES15.7,2X))') mm_btemp(i),mm_plev(i),mm_zlev(i),mm_dzlev(i) ELSE WRITE(*,'(3(ES15.7,2X))') mm_btemp(i),mm_plev(i),mm_zlev(i) ENDIF ENDDO ENDIF RETURN END FUNCTION mm_column_init FUNCTION mm_aerosols_init(m0aer_s,m3aer_s,m0aer_f,m3aer_f) RESULT(err) !! Initialize clouds tracers vertical grid. !! !! The subroutine initializes aerosols microphysics tracers columns. It allocates variables if !! required and stores input vectors in reversed order. It also computes the characteristic radii !! of each mode. !! @note !! All the input arguments should be defined from ground to top. !! !! @attention !! [[mm_globals(module):mm_global_init(interface)]] and [[mm_globals(module):mm_column_init(function)]] !! must have been called at least once before this method is called. Moreover, this method should be !! after each call of [[mm_globals(module):mm_column_init(function)]] to reflect changes in the !! vertical atmospheric structure. REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m0aer_s !! \(0^{th}\) order moment of the spherical mode (\(m^{-2}\)). REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m3aer_s !! \(3^{rd}\) order moment of the spherical mode (\(m^{3}.m^{-2}\)). REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m0aer_f !! \(0^{th}\) order moment of the fractal mode (\(m^{-2}\)). REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m3aer_f !! \(3^{rd}\) order moment of the fractal mode (\(m^{3}.m^{-2}\)). TYPE(error) :: err !! Error status of the function. err = noerror IF (.NOT.mm_ini) THEN err = error("mm_aerosols_init: Global initialization not done yet",-1) ; RETURN ENDIF IF (.NOT.mm_ini_col) THEN err = error("mm_aerosols_init: Column initialization not done yet",-1) ; RETURN ENDIF ! Check input size ??? IF (SIZE(m0aer_s) /= mm_nla) THEN err = error("mm_aerosols_init: Invalid size for input arrays",-1) ; RETURN ENDIF ! Allocate variable if required IF (.NOT.ALLOCATED(mm_m0aer_s)) ALLOCATE(mm_m0aer_s(mm_nla)) IF (.NOT.ALLOCATED(mm_m3aer_s)) ALLOCATE(mm_m3aer_s(mm_nla)) IF (.NOT.ALLOCATED(mm_m0aer_f)) ALLOCATE(mm_m0aer_f(mm_nla)) IF (.NOT.ALLOCATED(mm_m3aer_f)) ALLOCATE(mm_m3aer_f(mm_nla)) IF (.NOT.ALLOCATED(mm_rcs)) ALLOCATE(mm_rcs(mm_nla)) IF (.NOT.ALLOCATED(mm_rcf)) ALLOCATE(mm_rcf(mm_nla)) ! Allocate memory for diagnostics IF (.NOT.ALLOCATED(mm_m0as_vsed)) THEN ALLOCATE(mm_m0as_vsed(mm_nla)) ; mm_m0as_vsed(:) = 0._mm_wp ENDIF IF (.NOT.ALLOCATED(mm_m3as_vsed)) THEN ALLOCATE(mm_m3as_vsed(mm_nla)) ; mm_m3as_vsed(:) = 0._mm_wp ENDIF IF (.NOT.ALLOCATED(mm_m0af_vsed)) THEN ALLOCATE(mm_m0af_vsed(mm_nla)) ; mm_m0af_vsed(:) = 0._mm_wp ENDIF IF (.NOT.ALLOCATED(mm_m3af_vsed)) THEN ALLOCATE(mm_m3af_vsed(mm_nla)) ; mm_m3af_vsed(:) = 0._mm_wp ENDIF IF (.NOT.ALLOCATED(mm_aer_s_flux)) THEN ALLOCATE(mm_aer_s_flux(mm_nla)) ; mm_aer_s_flux(:) = 0._mm_wp ENDIF IF (.NOT.ALLOCATED(mm_aer_f_flux)) THEN ALLOCATE(mm_aer_f_flux(mm_nla)) ; mm_aer_f_flux(:) = 0._mm_wp ENDIF ! note : mm_dzlev is already from top to ground mm_m0aer_s = m0aer_s(mm_nla:1:-1)/mm_dzlev(:) mm_m3aer_s = m3aer_s(mm_nla:1:-1)/mm_dzlev(:) mm_m0aer_f = m0aer_f(mm_nla:1:-1)/mm_dzlev(:) mm_m3aer_f = m3aer_f(mm_nla:1:-1)/mm_dzlev(:) ! Setup threshold: call mm_set_moments_thresholds() ! aerosols characteristic radii WHERE(mm_m3aer_s > 0._mm_wp .AND. mm_m0aer_s > 0._mm_wp) mm_rcs = mm_get_rcs(mm_m0aer_s,mm_m3aer_s) ELSEWHERE mm_rcs = 0._mm_wp ENDWHERE WHERE(mm_m3aer_f > 0._mm_wp .AND. mm_m0aer_f > 0._mm_wp) mm_rcf = mm_get_rcf(mm_m0aer_f,mm_m3aer_f) ELSEWHERE mm_rcf = 0._mm_wp ENDWHERE mm_ini_aer = .true. END FUNCTION mm_aerosols_init FUNCTION mm_clouds_init(m0ccn,m3ccn,m3ice,gazs) RESULT(err) !! Initialize clouds tracers vertical grid. !! !! The subroutine initializes cloud microphysics tracers columns. It allocates variables if !! required and stores input vectors in reversed order. It also computes the mean drop radius !! and density and allocates diagnostic vectors. !! @note !! All the input arguments should be defined from ground to top. !! !! @attention !! [[mm_globals(module):mm_global_init(interface)]] and [[mm_globals(module):mm_column_init(function)]] !! must have been called at least once before this method is called. Moreover, this method should be !! after each call of [[mm_globals(module):mm_column_init(function)]] to reflect changes in the !! vertical atmospheric structure. REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m0ccn !! 0th order moment of the CCN distribution (\(m^{-2}\)). REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: m3ccn !! 3rd order moment of the CCN distribution (\(m^{3}.m^{-2}\)). REAL(kind=mm_wp), DIMENSION(:,:), INTENT(in) :: m3ice !! 3rd order moments of the ice components (\(m^{3}.m^{-2}\)). REAL(kind=mm_wp), DIMENSION(:,:), INTENT(in) :: gazs !! Condensible species gazs molar fraction (\(mol.mol^{-1}\)). TYPE(error) :: err !! Error status of the function. INTEGER :: i err = noerror IF (.NOT.mm_ini) THEN err = error("Global initialization not done yet",-8) RETURN ENDIF IF (.NOT.mm_w_clouds) THEN IF (mm_debug) WRITE(*,'(a)') "WARNING: Cloud microphysic is not enabled..." RETURN ENDIF ! Note: ! Here we could check that mm_nla is the same size of gazs(DIM=1) ! Actually, mm_nla should always initialized the first time mm_column_init is called, NOT HERE. IF (mm_nla < 0) mm_nla = SIZE(gazs,DIM=1) ! Note: ! here we could check that mm_nesp is the same size of gazs(DIM=2) ! Actually, mm_nesp should be always initialized in mm_global_init, NOT HERE. IF (mm_nesp < 0) mm_nesp = SIZE(gazs,DIM=2) ! Allocate variable if required IF (.NOT.ALLOCATED(mm_m0ccn)) ALLOCATE(mm_m0ccn(mm_nla)) IF (.NOT.ALLOCATED(mm_m3ccn)) ALLOCATE(mm_m3ccn(mm_nla)) IF (.NOT.ALLOCATED(mm_m3ice)) ALLOCATE(mm_m3ice(mm_nla,mm_nesp)) IF (.NOT.ALLOCATED(mm_gazs)) ALLOCATE(mm_gazs(mm_nla,mm_nesp)) IF (.NOT.ALLOCATED(mm_drad)) ALLOCATE(mm_drad(mm_nla)) IF (.NOT.ALLOCATED(mm_drho)) ALLOCATE(mm_drho(mm_nla)) ! Allocate memory for diagnostics IF (.NOT.ALLOCATED(mm_ccn_vsed)) THEN ALLOCATE(mm_ccn_vsed(mm_nla)) ; mm_ccn_vsed(:) = 0._mm_wp ENDIF IF (.NOT.ALLOCATED(mm_ccn_flux)) THEN ALLOCATE(mm_ccn_flux(mm_nla)) ; mm_ccn_flux(:) = 0._mm_wp ENDIF IF (.NOT.ALLOCATED(mm_ice_prec)) THEN ALLOCATE(mm_ice_prec(mm_nesp)) ; mm_ice_prec(:) = 0._mm_wp ENDIF IF (.NOT.ALLOCATED(mm_ice_fluxes)) THEN ALLOCATE(mm_ice_fluxes(mm_nla,mm_nesp)) ; mm_ice_fluxes(:,:) = 0._mm_wp ENDIF IF (.NOT.ALLOCATED(mm_gazs_sat)) THEN ALLOCATE(mm_gazs_sat(mm_nla,mm_nesp)) ; mm_gazs_sat(:,:) = 0._mm_wp ENDIF ! note mm_dzlev already from top to ground mm_m0ccn = m0ccn(mm_nla:1:-1)/mm_dzlev(:) mm_m3ccn = m3ccn(mm_nla:1:-1)/mm_dzlev(:) DO i=1,mm_nesp mm_m3ice(:,i) = m3ice(mm_nla:1:-1,i)/mm_dzlev(:) mm_gazs(:,i) = gazs(mm_nla:1:-1,i) ENDDO ! Setup threshold : call mm_set_moments_cld_thresholds() ! drop mean radius call mm_cloud_properties(mm_m0ccn,mm_m3ccn,mm_m3ice,mm_drad,mm_drho) mm_ini_cld = .true. END FUNCTION mm_clouds_init SUBROUTINE mm_dump_parameters() !! Dump model global parameters on stdout. WRITE(*,'(a)') "========= YAMMS PARAMETERS ============" WRITE(*,'(a,a)') "mm_fp_precision : ", mm_wp_s WRITE(*,'(a,L2)') "mm_debug : ", mm_debug WRITE(*,'(a,L2)') "mm_w_haze_prod : ", mm_w_haze_prod WRITE(*,'(a,ES14.7)') " mm_p_prod : ", mm_p_prod WRITE(*,'(a,ES14.7)') " mm_tx_prod : ", mm_tx_prod WRITE(*,'(a,ES14.7)') " mm_rc_prod : ", mm_rc_prod WRITE(*,'(a,L2)') "mm_w_haze_coag : ", mm_w_haze_coag WRITE(*,'(a,I2.2)') " mm_coag_interactions: ", mm_coag_choice WRITE(*,'(a,L2)') "mm_w_haze_sed : ", mm_w_haze_sed WRITE(*,'(a,L2)') " mm_wsed_m0 : ", mm_wsed_m0 WRITE(*,'(a,L2)') " mm_wsed_m3 : ", mm_wsed_m3 WRITE(*,'(a,L2)') " mm_no_fiadero_w : ", mm_no_fiadero_w WRITE(*,'(a,ES14.7)') " mm_fiadero_min : ", mm_fiadero_min WRITE(*,'(a,ES14.7)') " mm_fiadero_max : ", mm_fiadero_max WRITE(*,'(a,L2)') "mm_w_clouds : ", mm_w_clouds WRITE(*,'(a,L2)') " mm_w_cloud_sed : ", mm_w_cloud_sed WRITE(*,'(a,L2)') " mm_w_cloud_nucond : ", mm_w_cloud_nucond WRITE(*,'(a)') "---------------------------------------" WRITE(*,'(a)') "Thresholds spherical mode" WRITE(*,'(a,ES14.7)') " mm_m0as_min : ", mm_m0as_min WRITE(*,'(a,ES14.7)') " mm_rcs_min : ", mm_rcs_min WRITE(*,'(a)') "Thresholds fractal mode" WRITE(*,'(a,ES14.7)') " mm_m0af_min : ", mm_m0af_min WRITE(*,'(a,ES14.7)') " mm_rcf_min : ", mm_rcf_min WRITE(*,'(a)') "Thresholds clouds drop" WRITE(*,'(a,ES14.7)') " mm_m0n_min : ", mm_m0n_min WRITE(*,'(a,ES14.7)') " mm_drad_min : ", mm_drad_min WRITE(*,'(a,ES14.7)') " mm_drad_max : ", mm_drad_max WRITE(*,'(a)') "---------------------------------------" WRITE(*,'(a,ES14.7)') "mm_dt : ", mm_dt IF (mm_nla > -1) THEN WRITE(*,'(a,I3.3)') "mm_nla : ", mm_nla ELSE WRITE(*,'(a)') "mm_nla : not initialized yet" ENDIF WRITE(*,'(a,ES14.7)') "mm_df : ", mm_df WRITE(*,'(a,ES14.7)') "mm_rm : ", mm_rm WRITE(*,'(a,ES14.7)') "mm_rpla : ", mm_rpla WRITE(*,'(a,ES14.7)') "mm_g0 : ", mm_g0 WRITE(*,'(a)') "=======================================" END SUBROUTINE mm_dump_parameters SUBROUTINE mm_set_moments_thresholds() !! Apply minimum threshold for the aerosols moments. !! !! The method resets moments (for both modes and orders, 0 and 3) values to zero if !! their current value is below the minimum threholds. !! !! See also [[mm_globals(module):mm_m0as_min(variable)]], [[mm_globals(module):mm_rcs_min(variable)]], !! [[mm_globals(module):mm_rcf_min(variable)]] and [[mm_globals(module):mm_m0as_min(variable)]]. INTEGER :: i DO i=1,mm_nla IF ((mm_m0aer_s(i) < mm_m0as_min) .OR. (mm_m3aer_s(i) < mm_m3as_min)) THEN mm_m0aer_s(i) = 0._mm_wp ! mm_m0as_min mm_m3aer_s(i) = 0._mm_wp ! mm_m0as_min * mm_rcs_min**3._mm_wp * mm_alpha_s(3._mm_wp) ENDIF IF ((mm_m0aer_f(i) < mm_m0af_min) .OR. (mm_m3aer_f(i) < mm_m3af_min)) THEN mm_m0aer_f(i) = 0._mm_wp ! mm_m0af_min mm_m3aer_f(i) = 0._mm_wp ! mm_m0af_min * mm_rcf_min**3._mm_wp * mm_alpha_f(3._mm_wp) ENDIF ENDDO END SUBROUTINE mm_set_moments_thresholds SUBROUTINE mm_set_moments_cld_thresholds() !! Apply minimum threshold for the cloud drop moments. !! !! The method resets moments (for both modes and orders, 0 and 3) values to zero if !! their current value is below the minimum threholds. INTEGER :: i, j REAL(kind=mm_wp) :: m3cld DO i = 1, mm_nla m3cld = mm_m3ccn(i) DO j = 1, mm_nesp m3cld = m3cld + mm_m3ice(i,j) ENDDO IF ((mm_m0ccn(i) < mm_m0n_min) .OR. (m3cld < mm_m3cld_min)) THEN mm_m0ccn(i) = 0._mm_wp mm_m3ccn(i) = 0._mm_wp DO j = 1, mm_nesp mm_m3ice(i,j) = 0._mm_wp ENDDO ENDIF ENDDO END SUBROUTINE mm_set_moments_cld_thresholds ELEMENTAL SUBROUTINE mm_check_tendencies(v,dv) !! Check that tendencies is not greater than value. !! !! the purpose of the subroutine is to update dvalue so that v+dv is not negative. REAL(kind=mm_wp), INTENT(in) :: v !! Value to check. REAL(kind=mm_wp), INTENT(inout) :: dv !! Value tendencies to check and update consequently. REAL(kind=mm_wp), PARAMETER :: a = (epsilon(1._mm_wp)-1._mm_wp) IF (v+dv < 0._mm_wp) THEN dv = a*v ENDIF END SUBROUTINE mm_check_tendencies ELEMENTAL FUNCTION mm_get_rcs(m0,m3) RESULT(res) !! Get the characteristic radius for the spherical aerosols size distribution. !! !! The method computes the characteristic radius of the size distribution law !! of the spherical aerosols mode according to its moments and its inter-moments !! relation. REAL(kind=mm_wp), INTENT(in) :: m0 !! \(0^{th}\) order moment REAL(kind=mm_wp), INTENT(in) :: m3 !! \(3^{rd}\) order moment REAL(kind=mm_wp) :: res !! Radius res = (m3/m0/mm_alpha_s(3._mm_wp))**(1._mm_wp/3._mm_wp) END FUNCTION mm_get_rcs ELEMENTAL FUNCTION mm_get_rcf(m0,m3) RESULT(res) !! Get the characteristic radius for the fractal aerosols size distribution. !! !! The method computes the characteristic radius of the size distribution law !! of the fractal aerosols mode according to its moments and its inter-moments !! relation. REAL(kind=mm_wp), INTENT(in) :: m0 !! \(0^{th}\) order moment REAL(kind=mm_wp), INTENT(in) :: m3 !! \(3^{rd}\) order moment REAL(kind=mm_wp) :: res !! Radius res = (m3/m0/mm_alpha_f(3._mm_wp))**(1._mm_wp/3._mm_wp) END FUNCTION mm_get_rcf ELEMENTAL FUNCTION mm_effg(z) RESULT(effg) !! Compute effective gravitational acceleration. REAL(kind=mm_wp), INTENT(in) :: z !! Altitude in meters REAL(kind=mm_wp) :: effg !! Effective gravitational acceleration in \(m.s^{-2}\) effg = mm_g0 IF (mm_use_effg) effg = effg * (mm_rpla/(mm_rpla+z))**2 RETURN END FUNCTION mm_effg !================================== ! --- private methods ------------- !================================== SUBROUTINE cldprop_sc(m0ccn,m3ccn,m3ice,drad,drho) !! Get cloud drop properties (scalar). !! !! The method computes the mean radius and mean density of cloud drops. !! !! @bug !! A possible bug can happen because of threshold snippet. If __drad__ is greater than !! __drmax__ (== 1e3 microns) it is automatically set to __drmax__, but computation of !! __drho__ remains unmodified. So __drho__ is not correct in that case. !! !! @todo !! Fix the bug of the subroutine, but it is rather minor, since theoretically we do not !! need the density of the drop. !! !! @todo !! Think about a better implementation of thresholds, and get sure of their consequences in !! the other parts of the model. REAL(kind=mm_wp), INTENT(in) :: m0ccn !! \(0^{th}\) order moment of the ccn REAL(kind=mm_wp), INTENT(in) :: m3ccn !! \(3^{rd}\) order moment of the ccn REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: m3ice !! \(3^{rd}\) order moments of each ice component REAL(kind=mm_wp), INTENT(out) :: drad !! Output mean drop radius REAL(kind=mm_wp), INTENT(out), OPTIONAL :: drho !! Optional output mean drop density REAL(kind=mm_wp) :: Ntot, Vtot, Wtot REAL(kind=mm_wp), PARAMETER :: athird = 1._mm_wp / 3._mm_wp REAL(kind=mm_wp), PARAMETER :: pifac = (4._mm_wp * mm_pi) / 3._mm_wp ! Set to zero : drad = 0._mm_wp IF (PRESENT(drho)) drho = 0._mm_wp ! Initialization : Ntot = m0ccn Vtot = pifac*m3ccn + SUM(m3ice) Wtot = pifac*m3ccn*mm_rhoaer + SUM(m3ice*mm_xESPS(:)%rho) IF (Ntot <= mm_m0n_min .OR. Vtot <= mm_m3cld_min) THEN drad = mm_drad_min IF (PRESENT(drho)) drho = mm_rhoaer ELSE drad = (Vtot / (pifac*Ntot))**athird drad = MAX(MIN(drad,mm_drad_max),mm_drad_min) IF (PRESENT(drho)) drho = Wtot / Vtot ENDIF RETURN END SUBROUTINE cldprop_sc SUBROUTINE cldprop_ve(m0ccn,m3ccn,m3ice,drad,drho) !! Get cloud drop properties (vector). !! !! The method performs the same computations than [[mm_globals(module):cldprop_sc(subroutine)]] !! but for the entire vertical atmospheric structure. !! Same remarks apply here. REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: m0ccn !! 0th order moment of the ccn. REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: m3ccn !! 3rd order moment of the ccn. REAL(kind=mm_wp), INTENT(in), DIMENSION(:,:) :: m3ice !! 3rd order moments of each ice component. REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: drad !! Output mean drop radius. REAL(kind=mm_wp), INTENT(out), DIMENSION(:), OPTIONAL :: drho !! Optional output mean drop density. INTEGER :: i IF (PRESENT(drho)) THEN DO i = 1, SIZE(m0ccn) ; call cldprop_sc(m0ccn(i),m3ccn(i),m3ice(i,:),drad(i),drho(i)) ; ENDDO ELSE DO i = 1, SIZE(m0ccn) ; call cldprop_sc(m0ccn(i),m3ccn(i),m3ice(i,:),drad(i)) ; ENDDO ENDIF RETURN END SUBROUTINE cldprop_ve ! For configuration file (requires swift library). FUNCTION read_esp(parser,sec,pp) RESULT (err) !! Read and store [[mm_globals(module):mm_esp(type)]] parameters. TYPE(cfgparser), INTENT(in) :: parser !! Configuration parser. CHARACTER(len=*), INTENT(in) :: sec !! Name of the specie (should match a section of the configuration. TYPE(mm_esp), INTENT(out) :: pp !! [[mm_globals(module):mm_esp(type)]] object that stores the parameters. TYPE(error) :: err !! Error status of the function. err = cfg_get_value(parser,TRIM(sec)//'/name',pp%name) ; IF (err /= 0) RETURN err = cfg_get_value(parser,TRIM(sec)//'/mas',pp%mas) ; IF (err /= 0) RETURN err = cfg_get_value(parser,TRIM(sec)//'/vol',pp%vol) ; IF (err /= 0) RETURN err = cfg_get_value(parser,TRIM(sec)//'/ray',pp%ray) ; IF (err /= 0) RETURN err = cfg_get_value(parser,TRIM(sec)//'/mas',pp%mas) ; IF (err /= 0) RETURN err = cfg_get_value(parser,TRIM(sec)//'/vol',pp%vol) ; IF (err /= 0) RETURN err = cfg_get_value(parser,TRIM(sec)//'/ray',pp%ray) ; IF (err /= 0) RETURN err = cfg_get_value(parser,TRIM(sec)//'/masmol',pp%masmol) ; IF (err /= 0) RETURN err = cfg_get_value(parser,TRIM(sec)//'/rho',pp%rho) ; IF (err /= 0) RETURN err = cfg_get_value(parser,TRIM(sec)//'/tc',pp%tc) ; IF (err /= 0) RETURN err = cfg_get_value(parser,TRIM(sec)//'/pc',pp%pc) ; IF (err /= 0) RETURN err = cfg_get_value(parser,TRIM(sec)//'/tb',pp%tb) ; IF (err /= 0) RETURN err = cfg_get_value(parser,TRIM(sec)//'/w',pp%w) ; IF (err /= 0) RETURN err = cfg_get_value(parser,TRIM(sec)//'/a_sat',pp%a_sat) ; IF (err /= 0) RETURN err = cfg_get_value(parser,TRIM(sec)//'/b_sat',pp%b_sat) ; IF (err /= 0) RETURN err = cfg_get_value(parser,TRIM(sec)//'/c_sat',pp%c_sat) ; IF (err /= 0) RETURN err = cfg_get_value(parser,TRIM(sec)//'/d_sat',pp%d_sat) ; IF (err /= 0) RETURN err = cfg_get_value(parser,TRIM(sec)//'/mteta',pp%mteta) ; IF (err /= 0) RETURN err = cfg_get_value(parser,TRIM(sec)//'/tx_prod',pp%tx_prod) ; IF (err /= 0) RETURN RETURN END FUNCTION read_esp ! ========================================================================= ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! CONFIGURATION PARSER checking methods ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ========================================================================= FUNCTION check_r1(err,var,def,wlog) RESULT(ret) !! Check an option value (float). !! !! The method checks an option value and optionally set a default value, __def__ to initialize !! __var__ on error if given. TYPE(error), INTENT(in) :: err !! Error object from value getter. REAL(kind=mm_wp), INTENT(inout) :: var !! Input/output option value. REAL(kind=mm_wp), INTENT(in), OPTIONAL :: def !! Default value to set. LOGICAL, INTENT(in), OPTIONAL :: wlog !! .true. to print warning/error message. TYPE(error) :: ret !! Input error. CHARACTER(len=*), PARAMETER :: defmsg = '... Using default value: ' LOGICAL :: zlog ret = err zlog = .false. ; IF (PRESENT(wlog)) zlog = wlog IF (err == 0) RETURN IF (PRESENT(def)) THEN var = def IF (zlog) WRITE(*,'(a,a,a)') error_to_string(err,'',.true.),defmsg,to_string(var) ret = noerror ELSE IF (zlog) WRITE(*,'(a)') error_to_string(err,'',.true.) ENDIF END FUNCTION check_r1 FUNCTION check_l1(err,var,def,wlog) RESULT(ret) !! Check an option value (logical). !! !! The method checks an option value and optionally set a default value, __def__ to initialize !! __var__ on error if given. TYPE(error), INTENT(in) :: err !! Error object from value getter. LOGICAL, INTENT(inout) :: var !! Input/output option value. LOGICAL, INTENT(in), OPTIONAL :: def !! Default value to set. LOGICAL, INTENT(in), OPTIONAL :: wlog !! .true. to print warning/error message. TYPE(error) :: ret !! Input error. CHARACTER(len=*), PARAMETER :: defmsg = '... Using default value: ' LOGICAL :: zlog ret = err zlog = .false. ; IF (PRESENT(wlog)) zlog = wlog IF (err == 0) RETURN IF (PRESENT(def)) THEN var = def IF (zlog) WRITE(*,'(a,a,a)') error_to_string(err,'',.true.),defmsg,to_string(var) ret = noerror ELSE IF (zlog) WRITE(*,'(a)') error_to_string(err,'',.true.) ENDIF END FUNCTION check_l1 FUNCTION check_i1(err,var,def,wlog) RESULT(ret) !! Check an option value (integer). !! !! The method checks an option value and optionally set a default value, __def__ to initialize !! __var__ on error if given. TYPE(error), INTENT(in) :: err !! Error object from value getter. INTEGER, INTENT(inout) :: var !! Input/output option value. INTEGER, INTENT(in), OPTIONAL :: def !! Default value to set. LOGICAL, INTENT(in), OPTIONAL :: wlog !! .true. to print warning/error message. TYPE(error) :: ret !! Input error. CHARACTER(len=*), PARAMETER :: defmsg = '... Using default value: ' LOGICAL :: zlog ret = err zlog = .false. ; IF (PRESENT(wlog)) zlog = wlog IF (err == 0) RETURN IF (PRESENT(def)) THEN var = def IF (zlog) WRITE(*,'(a,a,a)') error_to_string(err,'',.true.),defmsg,to_string(var) ret = noerror ELSE IF (zlog) WRITE(*,'(a)') error_to_string(err,'',.true.) ENDIF END FUNCTION check_i1 FUNCTION check_s1(err,var,def,wlog) RESULT(ret) !! Check an option value (string). !! !! The method checks an option value and optionally set a default value, __def__ to initialize !! __var__ on error if given. TYPE(error), INTENT(in) :: err !! Error object from value getter. CHARACTER(len=*), INTENT(inout) :: var !! Input/output option value. CHARACTER(len=*), INTENT(in), OPTIONAL :: def !! Default value to set. LOGICAL, INTENT(in), OPTIONAL :: wlog !! .true. to print warning/error message. TYPE(error) :: ret !! Input error. CHARACTER(len=*), PARAMETER :: defmsg = '... Using default value: ' LOGICAL :: zlog ret = err zlog = .false. ; IF (PRESENT(wlog)) zlog = wlog IF (err == 0) RETURN IF (PRESENT(def)) THEN var = TRIM(def) IF (zlog) WRITE(*,'(a,a,a)') error_to_string(err,'',.true.),defmsg,var ret = noerror ELSE IF (zlog) WRITE(*,'(a)') error_to_string(err,'') ENDIF RETURN END FUNCTION check_s1 END MODULE MM_GLOBALS