mainpage.F90

!> \mainpage NMMB-MONARCH Code Documentation
!! The Multiscale Online Nonhydrostatic AtmospheRe CHemistry model
!! (NMMB-MONARCH), a chemical weather prediction system formerly known as
!! NMMB/BSC-CTM, that can be run either globally or regionally. The
!! NMMB-MONARCH, developed at the Barcelona Supercomputing Center, is based on
!! the online coupling of the meteorological Nonhydrostatic Multiscale Model on
!! the B-grid \cite Janjic2012 with a full aerosol-chemistry module.
!! This documentation is under development and in its current state is
!! not comprehensive and may describe developments not completely included in 
!! the present release or miss some of the ones included.
!!
!! \htmlonly
!! <div class="row"><div class="col-sm-4">
!! \endhtmlonly
!!     \anchor main_what_is_monarch
!!     ### What is MONARCH ###
!!
!!     - \subpage what_is_monarch "Meteorology"
!!     - \ref what_is_monarch_chemistry "Chemistry"
!!     - \ref what_is_monarch_feedbacks "Feedbacks"
!!
!!     \anchor main_how_monarch_works
!!     ### How MONARCH works ###
!!
!!     - \subpage how_monarch_works "Preprocessing"
!!     - \ref how_monarch_works_model "Model"
!!
!! \htmlonly
!! </div><div class="col-sm-4">
!! \endhtmlonly
!!
!!     \anchor main_how_tos
!!     ### How To's ###
!!
!!     - \ref running_monarch "Running MONARCH"
!!
!!     See the full list \ref howto_user "here".
!!
!!     \anchor main_config_formats
!!     ### Configuration Formats ###
!!
!!       See \ref config_formats for a descripton of the configuration
!!       options for MONARCH
!!
!! \htmlonly
!! </div><div class="col-sm-4">
!! \endhtmlonly
!!
!!     \anchor main_publications
!!     ### Publications ###
!!
!!     - Badia and Jorba 2015 \cite Badia2015
!!     - Badia et al. 2017 \cite Badia2017
!!     - Badia 2014 \cite Badia2014
!!     - Basart et al. 2016 \cite Basart2016
!!     - DiTomaso et al. 2017 \cite DiTomaso2017
!!     - Gama et al. 2015 \cite Gama2015
!!     - Gkikas et al. 2018 \cite Gkikas2018
!!     - Haustein et al. 2015 \cite Haustein2012
!!     - Jorba et al. 2012 \cite Jorba2012
!!     - Marti et al. 2017 \cite Marti2017
!!     - P&eacute;rez et al. 2011 \cite Perez2011
!!     - P&eacute;rez Garc&iacute;a-Pando et al. 2014 \cite PerezGarcia-Pando2014
!!     - Spada et al. 2013 \cite Spada2013
!!     - Spada 2015 \cite Spada2015a
!!
!! \htmlonly
!! </div></div>
!! <div class="header"><div class="headertitle">
!! <div class="title">For Developers</div>
!! </div></div>
!! <div class="row"><div class="col-sm-4">
!! \endhtmlonly
!!
!!     \anchor main_dev_how_tos
!!     ### How To's ###
!!
!!     - \ref coding_style "Coding style"
!!     - \ref adding_documentation "Documenting the code"
!!
!!       See the full list \ref howto_developer "here".
!!
!! \htmlonly
!! </div><div class="col-sm-4">
!! <h3>Ongoing Work</h3>
!! <p>
!!   <a class="el" href="todo.html">Todo List</a>
!! </p><p>
!!   <a class="el" href="bug.html">Bug List</a>
!! </p>
!! </div><div class="col-sm-4">
!! <h3>License</h3>
!! \endhtmlonly
!!
!!     GNU General Public License version 3.
!!
!! \htmlonly
!! </div></div>
!! \endhtmlonly

!******************************************************************************

!> \page what_is_monarch What is NMMB-MONARCH
!!
!! The Multiscale Online Nonhydrostatic AtmospheRe CHemistry model (NMMB-MONARCH) 
!! is a fully online multiscale chemical weather prediction system for regional 
!! and global-scale applications developed at the Barcelona Supercomputing Center 
!! (Pérez et al., 2011; Jorba et al., 2012). The model couples online 
!! a meteorological driver with the gas-phase and aerosol continuity equations 
!! to solve the atmospheric chemistry processes in detail. The model is designed 
!! to account for the feedbacks among gases, aerosol particles and meteorology.
!!
!! A description of the major components of NMMB-MONARCH follows.
!!
!!
!! ## Meteorology ##
!!
!! The meteorological driver of the NMMB-MONARCH is the Nonhydrostatic Multiscale Model 
!! on the B-grid (NMMB; Janjic and Black, 2005; Janjic and Gall, 2012) developed at NCEP. 
!! The NMMB is conceived for short- and medium-range forecasting over a wide range of spatial
!! and temporal scales, from large eddy simulations (LES) to global simulations. 
!! Its unified nonhydrostatic dynamical core allows for running either regional or global simulations,
!! both including embedded regional nests. The NMMB has been developed within the Earth 
!! System Modeling Framework (ESMF).
!!
!! ### Dynamics ###
!!
!! The Dynamical core of the NMMB is described in Janjic and Gall (2012). 
!! The NMMB has been developed following the general modeling philosophy of 
!! the NCEP regional Weather Research and Forecasting (WRF) Nonhydrostatic Mesoscale Model 
!! (NMM; Janjic et al., 2001; Janjic, 2003). The numerical schemes used in the model were designed
!! following the principles presented in Janjic (1977, 1979, 1984, 2003). Isotropic horizontal finite volume differencing
!! is employed so a variety of basic and derived dynamical and quadratic quantities are conserved. Among these, the conservation
!! of energy and enstrophy (Arakawa, 1966) improves the accuracy of the nonlinear dynamics. The hybrid pressure-sigma
!! coordinate is used in the vertical direction and the Arakawa B-grid is applied in the horizontal direction. The global model on
!! the latitude-longitude grid with polar filtering was developed as the reference version, and other geometries such the cubedsphere
!! are currently being tested. The regional model is formulated on a rotated longitude-latitude grid, with the Equator of
!! the rotated system running through the middle of the integration domain resulting in more uniform grid distances. The nonhydrostatic
!! component of the model dynamics is introduced through an add-on module that can be turned on or off, depending
!! on the resolution.
!!
!! ### Physics ###
!!
!! The NMMB has several different physical packages, the operational configuration includes: (1) the Mellor-Yamada-Janjic (MYJ) level 2.5 turbulence closure
!! for the treatment of turbulence in the Planetary Boundary Layer (PBL) and in the free atmosphere (Janjic et al., 2001), (2)
!! the surface layer scheme based on the Monin-Obukhov similarity theory (Monin and Obukhov, 1954) with introduced viscous
!! sublayer over land and water (Zilitinkevich, 1965; Janjic, 1994), (3) the NCEP NOAH (Ek et al., 2003) or the LISS land surface
!! model (Vukovic et al., 2010) for the computation of the heat and moisture surface fluxes, (4) the GFDL or RRTMG long-wave
!! and shortwave radiation package (Fels and Schwarzkopf, 1975; Mlawer et al., 1997), (5) the Ferrier gridscale clouds and
!! microphysics (Ferrier et al., 2002), and (6) the Betts-Miller-Janjic convective parametrization (Betts, 1986; Betts and Miller,
!! 1986; Janjic, 1994, 2000). Vertical diffusion is handled by the surface layer scheme and by the PBL scheme. Lateral diffusion
!! is formulated following the Smagorinsky non-linear approach (Janjic, 1990).
!!
!! \todo Further details on each physical package should be found in the following links (currently not working):
!!
!!     - \ref RADIATION "Radiation"
!!     - \ref TURBULENCE "Turbulence"
!!     - \ref MICROPHYSICS "Microphysics grid scale clouds"
!!     - \ref CONVECTION "Convective sub-grid clouds" 
!!
!! \anchor what_is_monarch_chemistry
!! ## Chemistry ##
!!
!! ### Emissions ###
!!
!! The High-Elective Resolution Modelling Emission System version 3 (HERMESv3; Guevara et al., 2019) is used to pre-process
!! the anthropogenic, biomass burning, soil, and ocean emissions for the NMMB-MONARCH model. HERMESv3 is an open source, parallel 
!! and stand-alone multiscale atmospheric emission modelling framework that processes gaseous and aerosol emissions for use in 
!! atmospheric chemistry models. The user can flexibly define combinations of existing up-to-date global and regional emission 
!! inventories and apply country specific scaling factors and masks. Each emission inventory is individually remapped onto the desired
!! destination grid and processed using user-defined vertical, temporal and speciation profiles that allow obtaining emission outputs 
!! compatible with multiple chemical mechanisms. The selection and combination of emission inventories and databases is done through 
!! detailed configuration files.
!!
!! Further details on HERMESv3 can be found in https://earth.bsc.es/gitlab/es/hermesv3_gr/
!!
!! The biogenic emissions are computed online within NMMB-MONARCH by using the Model of Emissions of Gases and Aerosols from Nature
!! (MEGAN; Guenther et al. 2006) version 2.04. MEGAN estimates the net emission rate of gases (more than 130 Non-Methane Volatile Organic 
!! Compounds (NMVOCs) and aerosols from terrestrial ecosystems into the above-canopy atmosphere. The model uses the 2 meter surface temperature
!! and shortwave incoming radiation from the meteorological driver.
!!
!! ### Photolysis ###
!!
!! To compute the photolysis rates, the NMMB-MONARCH uses online the Fast-J (Wild et al., 2000) photolysis scheme. Fast-J has been coupled 
!! with the physics of each model layer (e.g., clouds and absorbers such as O3). The optical depths of grid-scale clouds from the atmospheric 
!! driver are considered by using the fractional cloudiness based on relative humidity. 
!! The Fast-J scheme has been upgraded with CB05 photolytic reactions. The quantum yields and cross section for the
!! CB05 photolysis reactions have been revised and updated following the recommendations of Atkinson et al. (2004) and Sander
!! et al. (2006). 
!!
!! ### Aerosol sedimentation ###
!!
!! Sedimentation or gravitational settling is the most efficient removal process for large aerosols. The NMMB-MONARCH solves implicitly the
!! sedimentation for the aerosol mixing ratio from column top to bottom. The gravitational settling velocity is calculated following 
!! the Stokes-Cunningham approximation considering the size, diameter, and density of each aerosol component. The Cunningham correction factor 
!! accounts for the reduced resistance of viscosity as particle size approaches the mean free path of air molecules.
!!
!! ### Dry deposition ###
!! #### Gas ####
!!
!! The dry deposition scheme follows the classical deposition velocity analogy, enabling the calculation of deposition fluxes from airborne
!! concentrations. The aerodynamic resistance (Ra; depends only on atmospheric conditions) and the quasilaminar sublayer resistance (Rb; 
!! depends on friction velocity and modelcular characteristics of gases) are computed following their common definition (Seinfeld and Pandis, 1998),
!! while the canopy or surface resistance (Rc) is simulated following Wesely (1989), where Rc is derived from the resistances of the surfaces
!! of the soil and the plants. The properties of the plants are determined using land-use data (from the land use of the meteorological driver) 
!! and depend on the season. The surface resistance also depends on the diffusion coefficient, the reactivity, and the water solubility of the reactive trace gases.
!!
!! #### Aerosols ####
!!
!! The aerosol dry deposition is based on Zhang et al. (2001) which includes simplified empirical parameterizations for the deposition processes of
!! Brownian diffusion, impaction, interception and gravitational settling detailed in Slinn (1982). Aerosol rebound at the surface is not taken into 
!! account due to limited knowledge of this process. 
!!
!! ### Wet deposition ###
!! #### Gas ####
!!
!! The gas-phase cloud-chemistry processes are included in the system considering both sub-grid and grid-scale processes following Byun and Ching (1999)
!! and Foley et al. (2010). The processes included are scavenging, vertical mixing and wet-deposition. Only in-cloud scavenging is considered currently.
!!
!! #### Aerosols ####
!!
!! Wet scavenging of aerosols by precipitation is computed separately for convective and grid-scale (stratiform) precipitation. It represents the most 
!! efficient process for the deposition of the smallest particles. The model includes parameterizations for in-cloud scavenging, and for below cloud scavenging.
!! Detailed description of the schemes can be found in Pérez et al. (2011). 
!!
!! ### Chemistry mechanisms ###
!!
!! Two different different solutions are implementated in NMMB-MONARCH to solve a chemistry mechanism: (1) Fix code, and (2) Run-time configured mechanism.
!!
!! The Fix implementation uses the Carbon Bond 2005 chemical mechanism (CB05; Yarwood, 2005) extended with Toluene and Chlorine chemistry. Two options for 
!! solving the gas-phase chemistry are implemented: (1) an efficient and fast solver based on the Euler-Backward-Iterative shceme, and (2) a chemical mechanism
!! and chemistry solver based on the Kinetic PreProcessor KPP package with the main purpose of maintaining a wide flexibility when configuring the model. 
!! The CB05 is well formulated for urban to remote tropospheric conditions and it considers 51 chemical species and solves 156 reactions.
!!
!! With the aim of extending the flexibility of the model, the NMMB-MONARCH includes a novel option of adding run-time configured mechanisms for gas- and aerosol-phase
!! chemistry and mass transfer using the \ref module_phlex_chem "Chemistry Across Multiple Phases (CAMP)" library.
!!
!! ### Aerosol representation ###
!!
!! The aerosol scheme of the NMMB-MONARCH model follows a hybrid sectional-bulk multicomponent approach.
!! It describes the life cycle of mineral dust, sea-salt, black carbon, organic matter (both primary and secondary), sulfate, nitrate, and ammonium aerosols.
!! While a sectional approach is used for mineral dust and sea-salt, a bulk description of the other aerosol species is adopted. A simplified gas-aqueous-aerosol mechanism
!! has been introduced in the module to account for the sulfur chemistry, and a two-product scheme or a Simple non-volatile scheme is used for the formation of secondary 
!! organic aerosols. The EQSAM thermodynamic equilibrium model (Metzger et al., 2002) is used to solve the production of secondary inorganic nitrate and ammonium.
!! Mineral dust is treated as an hydrophobic particle, sea salt as hydrophilic, both black carbon and primary organic aerosol are emitted as hydrophobic particles and
!! an aging process transforms them to hydrophilic, secondary organic aerosol (SOA) is considered hydrophilic, and sulfate/nitrate/ammonium are treated as a hydrophilic particles.
!!
!! To extend the flexibility of the model when solving different types of problems, the \ref module_phlex_chem "CAMP" development allows an abstract non-fixed 
!! representation of the aerosols that can be configured during run-time.
!!
!!
!! \anchor what_is_monarch_feedbacks
!! ## Feedbacks ##
!! ### Aerosol&ndash;Radiation Interactions ###
!!
!! Documentation on Aerosol&ndash;Radiation Interactions (ARI) in MONARCH can be found \ref ARI_vobiso "here".
!!
!! ### Aerosol&ndash;Cloud Interactions ###
!!
!! No aerosol-cloud-interaction is implemented in the NMMB-MONARCH, yet.
!!

!*******************************************************************************

!> \page how_monarch_works How NMMB-MONARCH works
!!
!! \todo Finish the technical description for MONARCH
!!
!!
!! ## Preprocessing ##
!!
!!   <I>Document me!</I>
!!
!! \anchor how_monarch_works_model
!! ## Model ##
!!
!!   <I>Document me!</I>
!!

!*******************************************************************************

!> \page how_tos How To MONARCH
!!
!! \anchor howto_user
!! ## Users ##
!!
!! - \subpage running_monarch "Running MONARCH"
!!
!! \anchor howto_developer
!! ## Developers ##
!!
!! - \subpage coding_style "Coding style"
!! - \subpage adding_documentation "Documenting the code"
!! - \subpage using_diagnostics "Using diagnostics"
!! - \subpage using_column_diagnostics "Using column diagnostics"
!! - \subpage adding_diagnostics "Adding new diagnostics"
!! - \subpage adding_column_diagnostics "Adding new column diagnostics"
!! - \subpage setup_output_diagnostics "Setting up a new diagnostic for output"
!!

!*******************************************************************************

!> \page config_formats Configuration File Formats
!!
!! \todo Develop a structure for documenting configuration files for MONARCH
!! and update existing configuration format pages
!! (see issue:
!! https://earth.bsc.es/gitlab/ac/nmmb-monarch/issues/101 )
!!
!! - \subpage config_output_diagnostics "Output diagnostics"
!! - \subpage config_dry_deposition "Dry deposition"
!! - \subpage config_wet_deposition "Wet deposition"

!*******************************************************************************

!> \page ARI_vobiso ARI implementation in MONARCH
!! 
!! \htmlonly
!! <div class="row"><div class="col-sm-10">
!! \endhtmlonly
!! 
!! ### MODIFIED FILES ###
!! 
!! - nmmb-monarch/MODEL/SRC/nmmbphys/
!!   + radiation_aerosols_nmmb.f
!!   + radsw_main_nmmb.f
!!   + grrad_nmmb.f
!!   + rad_initialize_nmmb.f
!!   + module_RA_RRTM.F90
!! 
!! - nmmb-monarch/MODEL/SRC/
!!   + module_RADIATION.F90
!!   + module_SOLVER_GRID_COMP.F90
!!   + module_SOLVER_INTERNAL_STATE.F90
!!   + module_GET_CONFIG.F90
!!   + module_WRITE_GRID_COMP.F90
!!   + module_WRITE_ROUTINES.F90
!!   + module_WRITE_INTERNAL_STATE.F90
!! 
!! - nmmb-monarch/MODEL/NAMELIST/
!!   + solver_state.*.txt
!! 
!! - nmmb-monarch/JOB/TEMPLATE/
!!   + configfile_rrtm_chem.*.tmp
!!   + NMMBrrtm_RUN_TEMPLATE.sh
!! 
!! - auto-monarch/templates/conf/
!!   + proj_*.conf
!! 
!! - auto-monarch/templates/
!!   + 05b_sim.sh
!! 
!! 
!! <br/> 
!! ### MAIN CODE FEATURES ###
!! 
!! The ARI implementation in MONARCH is based on two main routine fluxes: (I) initialization and
!! (II) calculation of aerosol layer optical properties (optical depth: TAU; single scattering
!! albedo: SSA; asymmetry factor: ASY). The routines involved in these two processes are:
!! 
!! <pre>
!! (I) : SOLVER_INITIALIZE (module_SOLVER_GRID_COMP.F90)
!!         -> PHYSICS_INITIALIZE (module_SOLVER_GRID_COMP.F90)
!!              -> rad_initialize_nmmb (rad_initialize_nmmb.f)
!!                   -> radinit_nmmb (grrad_nmmb.f)
!!                        -> aer_init (radiation_aerosols_nmmb.f)
!!                             -> mnch_aerinit (radiation_aerosols_nmmb.f)
!! 
!! (II): SOLVER_RUN (module_SOLVER_GRID_COMP.F90)
!!         -> RADIATION (module_RADIATION.F90)
!!              -> RRTM (module_RA_RRTM.F90)
!!                   -> grrad_nmmb (grrad_nmmb.f)
!!                        -> setaer (radiation_aerosols_nmmb.f)
!!                             -> mnch_optprops (radiation_aerosols_nmmb.f) </pre>
!! 
!! Flux (I) only occurs at the beginning of the simulation (before the first time step) while flux (II)
!! at each time step.
!! 
!! In order to assure flexibility, this implementation is based on single aerosol bins/modes and not on
!! aerosol species (as it was before). Note that in this document each aerosol size bin (for dust and sea
!! salt) and mode (for the other species) is generally indicated as aerosol "type".
!! 
!! 
!! <br/>
!! #### module_radiation_aerosols_nmmb ####
!! 
!! The <I>module_radiation_aerosols_nmmb</I> module (<I>radiation_aerosols_nmmb.f</I>) contains the 2 main routines for ARI:
!! <I>mnch_aerinit</I> and <I>mnch_optprops</I>, in which initialization and calculation, respectively, of MONARCH
!! aerosol optical properties are performed. These 2 routines are called ONLY if the ARI module flag <I>iaermdl == 3</I>
!! (<I>mdlari</I> in configfile) and are ALTERNATIVE to OPAC routines (which work for <I>iaermdl == 0</I>). This means
!! that (at the moment) it is not possible to mix MONARCH and OPAC aerosols (not trivial work, given the OPAC structure)
!! and so if the MONARCH ARI is selected no climatological background is considered to affect the radiative fluxes.
!! Note that a specific instruction in <I>PHYSICS_INITIALIZE</I> routine imposes <I>iaermdl == 0</I> (OPAC) if there
!! are not running aerosols from MONARCH (<I>int_state\%NUM_AERO == 0</I>). 
!! 
!! Important module variables are:
!! 
!! <pre>
!! setari (allocatable, save) -> re-mapping array (to read optical table according to ARI configuration)
!! intopt1 (allocatable, save) -> intensive optical properties mapped onto RRTMG SW/LW spectral bands
!! optfls(1:4) -> string array containing names of main optical tables available </pre>
!! 
!! 
!! <br/>
!! #### mnch_aerinit ####
!! 
!! The <I>mnch_aerinit</I> routine initializes the optical properties for MONARCH aerosols.
!! The relevant inputs are:
!! 
!! <pre>
!! NSPARI -> number of available aerosol species (nspari in configfile)
!! NBNARI -> number of available aerosol types (nbnari in configfile)
!! numari(1:NSPARI) -> number of running types for each aerosol species (manually filled with num_* in PHYSICS_INITIALIZE routine)
!! nxbari(1:NSPARI) -> number of available types for each aerosol species (nxbari in configfile)
!! fhbari(1:NBNARI) -> ARI configuration for each type (fhbari in configfile)
!! iaertbl -> flag for optical table (tblari in configfile) </pre>
!! 
!! As outputs this routine fills the module arrays <I>setari</I> and <I>intopt1</I> for subsequent runtime calls. 
!! 
!! At first, re-mapping for ARI is performed. The array <I>setari</I> is allocated (<I>setari(1:sum(numari),1:2)</I>,
!! that is as first dimension it takes the same dynamic size of the aerosol tracer array) and filled:
!! 
!! <pre>
!! setari(:,1) -> id numbers of running types indicating their position in optical table
!! setari(:,2) -> fhbari values for running types (0 -> no ARI, 1 -> hydrophobic, 2 -> hydrophilic) </pre>
!! 
!! Note that, since the aerosol types do not have yet an own id number, the re-mapping is done by following
!! the order in the tracer array. In other words, the id number that each type gets here is the position it
!! would have in the tracer array if all the types were running (the optical table follows the same order).
!! 
!! Then, the actual optical table (<I>opttbl</I>) to be read is selected depending on <I>iaertbl</I>:
!! 
!! <pre>
!! opttbl = optfls(iaertbl) </pre>
!! 
!! Finally, two internal subroutines are called:
!! 
!! <pre>
!! call intprops -> reads optical table and calculates flux weights at each OPAC wavenumber
!! call maprrtmg -> maps intensive properties onto RRTMG SW/LW spectral bands and stores them in intopt1 </pre>
!! 
!! Note that mapping the intensive optical properties means calculating averages for each band, weighted
!! by the solar flux. The array <I>intopt1</I> is previously allocated as:
!! 
!! <pre>
!! intopt1(1:3,1:NRHARI,1:NSWLWBD,1:NBNARI) </pre>
!! 
!! where the first dimension is for the 3 optical properties (TAU, SSA, ASY), <I>NRHARI</I> is the number of
!! relative humidity (RH) levels for optical properties and <I>NSWLWBD</I> the number of SW/LW spectral bands.
!! 
!! 
!! <br/>
!! #### mnch_optprops ####
!! 
!! The <I>mnch_optprops</I> routine computes layer optical properties for MONARCH aerosols in SW/LW bands.
!! 
!! The relevant inputs are:
!! 
!! <pre>
!! rhlay(:,:) -> layer mean RH
!! dz(:,:) -> layer thickness (km)
!! RTPARI -> number of running aerosol types (int_state\%NUM_AERO)
!! trcari(:,:,:,1:RTPARI) -> MONARCH aerosol mass concentrations (kg/m³) (int_state\%TRACERS(:,:,:,INDX_AERO_START:INDX_AERO_END)) </pre>
!! 
!! The outputs are the aerosol optical properties (<I>aerosw</I> and <I>aerolw</I>) for each layer and spectral band:
!! 
!! <pre>
!! aerosw/aerolw(:,:,:,1) -> TAU
!! aerosw/aerolw(:,:,:,2) -> SSA
!! aerosw/aerolw(:,:,:,3) -> ASY </pre>
!! 
!! The RH levels containing <I>rhlay</I> (<I>ih1</I>, <I>ih2</I>) and their weight (<I>rdrh</I>)
!! are firstly computed for each layer (to interpolate the intensive optical properties).
!! 
!! The main loop over all the running types then starts, within layer and band loops, and the optical
!! properties of the aerosol mixture are calculated:
!! 
!! <pre>
!! do t = 1, max(1,RTPARI)
!!   b = setari(t,1)
!!   if ( setari(t,2) == 1 ) then   ! hydrophobic type
!!     do p = 1, 3
!!       iopts(p) = intopt1(p,1,nb,b)
!!     enddo
!!     mssari = trcari(i,k,t)
!!   elseif ( setari(t,2) == 2 ) then   ! hydrophilic type
!!     do p = 1, 3
!!       wopt1 = intopt1(p,ih1,nb,b)
!!       wopt2 = intopt1(p,ih2,nb,b)
!!       iopts(p) = wopt1 + rdrh*(wopt2 - wopt1)
!!     enddo
!!     mssari = trcari(i,k,t)
!!   ...
!!   else   ! setari(t,2) == 0 -> no ARI for this type
!!     mssari = f_zero
!!   endif
!!   mdz1 = mssari*dz1(k)*1000.0*1000.0   ! g/m²
!!   aod1 = aod1 + mdz1*iopts(1)   ! iopts(1) in m²/g
!!   ssa1 = ssa1 + mdz1*iopts(1)*iopts(2)
!!   asy1 = asy1 + mdz1*iopts(1)*iopts(2)*iopts(3)
!! enddo </pre>
!! 
!! Note that this calculation loop is without "cases". Thanks to <I>setari(t,1)</I> and <I>setari(t,2)</I> it is possible to know
!! "who is" each running type.
!! 
!! Finally, the arrays <I>aerosw</I> and <I>aerolw</I> are filled and returned to <I>grrad_nmmb</I>. These optical properties are
!! also stored as model outputs, column-integrated and at 0.550µm (actually at the RRTMG band containing 0.550µm):
!! 
!! <pre>
!! AERO_OPT_R(:,:,1) -> TAU_550
!! AERO_OPT_R(:,:,2) -> SSA_550
!! AERO_OPT_R(:,:,3) -> ASY_550 </pre>
!! 
!! 
!! <br/> 
!! #### mineralogy ####
!! 
!! The implementation of ARI with mineralogy in MONARCH has advanced but still requires to be refined and concluded
!! (because it depends on other ongoing implementations). The ARI code for 8 dust bins has been extended to 64 dust
!! types (8 components * 8 bins), 2 additional optical tables (<I>tblari: 3, 4</I>) have been generated for minerals
!! and a dynamic internal mixing mechanism (between hematite and clays: illite, kaolinite, montmorillonite) has been
!! implemented in both <I>mnch_aerinit</I> (initialization) and <I>mnch_optprops</I> (dynamic calculations) routines.
!! 
!! In case of external mixing (<I>tblari: 3</I>) between minerals, the code would work as for homogeneous dust
!! components. The table contains optical properties for the different minerals (in a fixed order as explained
!! later) and so in this case the 64 dust types would just contribute separately to the optical properties.
!! 
!! In case of internal mixing (<I>tblari: 4</I>), the enabled dynamic mechanism would allow selecting suitable
!! optical properties for the hematite-clays mixture depending on the actual volume ratio of hematite at each
!! time step (note that the other minerals are considered externally mixed in this case).
!! 
!! IMPORTANT: Due to the current lack of own id numbers for the aerosol types, the code is optimized to work with
!! ALL MINERALS and ALL BINS in the right order! If only some minerals or some bins for each mineral are activated
!! the model would not crash, but the optical properties would not be the right ones. If for example calcite and
!! quartz are activated (with 8 bins each one), the re-mapping would assign to these minerals just the first 16
!! positions available for dust and the optical table would be read accordingly (in case of external mixing the
!! first 16 dust positions of the table are for illite and kaolinite). For internal mixing the situation is even
!! more complex. This implementation will be improved when type id numbers will be available and num_* will disappear.
!! 
!! In <I>mnch_aerinit</I>, the input flag for mineralogy is:
!! 
!! <pre>
!! iaermin -> flag for dust mineralogy (now forced to 0, in future mine in configfile) </pre>
!! 
!! A specific instruction does not allow the tables for mineralogy being accidentally selected if minerals
!! are not running (in this case <I>iaertbl = 2</I> is set). Moreover, currently options with more than
!! 8 minerals are not allowed.
!! 
!! The dynamic internal mixing mechanism (enabled only if <I>iaertbl == 4</I>) requires at first 2 module
!! variables to be filled (for later use in <I>mnch_optprops</I>):
!! 
!! <pre>
!! tbl (save) -> flag for optical table (tbl = iaertbl)
!! hem(1:8), ilt(1:8), kao(1:8), mtm(1:8) (save) -> dynamic positions of hematite and clays in setari(:,1) </pre>
!! 
!! Once calculated <I>hem</I>, <I>ilt</I>, <I>kao</I> and <I>mtm</I>, the ARI configuration for internally mixed
!! dust types is OVERWRITTEN: from <I>setari(:,2) == 1</I> (hydrophobic) to <I>setari(:,2) == 3</I> (internally
!! mixed), in order to assure a special treatment for these types during the calculation of the optical properties.
!! 
!! The optical table for internal mixing contains a higher number of columns (currently 34) than the other tables
!! because it includes optical properties for 11 hematite/clays volume ratios: from 0.00 (base case) to 0.20 with
!! 0.02 step (I fixed the upper limit to 0.20 for the moment because in the emission dataset I found volume ratios
!! always < 0.16, but this has to be tested). So, in routines <I>intprops</I> and <I>maprrtmg</I> a specific code
!! is implemented to read and map onto RRTMG bands this special table.
!! 
!! In <I>mnch_optprops</I> (called for each time step and each grid cell), firstly, just before the start of the main
!! loop for the calculation of the optical properties, the actual volume ratios of hematite to clays are calculated,
!! separately for each bin, starting from tracers (dynamic) and mass densities (module parameters) of the involved
!! minerals, and stored in the auxiliary local variable <I>hemcly(1:RTPARI)</I>.
!! 
!! Then, within the calculation main loop, the option for internally mixed dust types has been inserted:
!! 
!! <pre>
!! ...
!! elseif ( setari(t,2) == 3 ) then   ! hydrophobic internally mixed type
!!   hc1 = int(hemcly(t)/HEMSTP)
!!   hc2 = hc1 + 1
!!   if ( hc2 > 10 ) then
!!     hc1 = 10
!!     hc2 = 10
!!     rdvr = 0.0
!!   else
!!     rdvr = (hemcly(t)-hc1*HEMSTP)/(hc2*HEMSTP-hc1*HEMSTP)
!!   endif
!!   do p = 1, 3
!!     xopt1 = intopt1(p+3*hc1,1,nb,b)
!!     xopt2 = intopt1(p+3*hc2,1,nb,b)
!!     iopts(p) = xopt1 + rdvr*xopt2
!!   enddo
!!   mssari = trcari(i,k,t)
!! ... </pre>
!! 
!! Substantially, the integer hematite levels (<I>hc1</I>, <I>hc2</I>) are firstly calculated (<I>HEMSTP</I> is the
!! hematite/clays ratio step, set to 0.02 as already reported) and then used to pick up the right columns in <I>intopt1</I>.
!! The actual optical properties are thus calculated by linearly interpolating these columns based on the level weight (<I>rdvr</I>).
!! 
!! 
!! <br/>
!! #### radiation multiple call ####
!! 
!! It is possible to activate the radiation multiple call to compute instantaneous Direct Radiative Effect (DRE)
!! of the running aerosols (at each time step and with fixed meteorological conditions). Note that this option
!! does not affect the aerosol impact on radiation, it just provides aerosol DRE in 4 model output variables:
!! 
!! <pre>
!! DRE_SW_TOA(:,:,1:(RTPARI+1)) -> aerosol direct radiative effect for SW at TOA (w/m²)
!! DRE_SW_SFC(:,:,1:(RTPARI+1)) ->               "                 for SW at SFC (w/m²)
!! DRE_LW_TOA(:,:,1:(RTPARI+1)) ->               "                 for LW at TOA (w/m²)
!! DRE_LW_SFC(:,:,1:(RTPARI+1)) ->               "                 for LW at SFC (w/m²) </pre>
!! 
!! The calculation is performed in the <I>RRTM</I> routine (<I>module_RA_RRTM.F90</I>) through single/double/multiple calls of
!! <I>grrad_nmmb</I>, depending on the flag <I>iaermlc</I> (<I>mlcari</I> in configfile). Currently, 4 options are available:
!! 
!! <pre>
!! iaermlc == 0 -> 1 call: all tracers
!! iaermlc == 1 -> 2 calls: no tracers + all tracers
!! iaermlc == 2 -> (RTPARI+2) calls: no tracers + tracer of each running type + all tracers
!! iaermlc == 3 -> (RTPARI+2) calls: no tracers + all tracers except each running type + all tracers </pre>
!! 
!! When <I>iaermlc > 0</I>, the SW/LW DREs at the top-of-atmosphere (<I>toa</I>) and surface (<I>sfc</I>) are calculated
!! as flux (<I>F</I>) differences with (<I>a</I>) and without (<I>0</I>) aerosols, just after the <I>grrad_nmmb</I> calls:
!! 
!! \f[
!! \begin{align}
!! DRE_{toa} &=  F_{toa}^{\uparrow 0} - F_{toa}^{\uparrow a} \\
!! \\
!! DRE_{sfc} &= (F_{sfc}^{\downarrow a} - F_{sfc}^{\uparrow a}) - (F_{sfc}^{\downarrow 0} - F_{sfc}^{\uparrow 0})
!! \end{align}
!! \f]
!! 
!! Note that in this case (<I>iaermlc > 0</I>), always the first call is without aerosols and the last call with all types.
!! This implies two things:
!! - the output radiative fluxes are not affected by <I>iaermlc</I>;
!! - the DRE of all aerosols is always stored in the last position of the DRE arrays (<I>RTPARI+1</I>).
!! 
!! The options <I>iaermlc == 2</I> and <I>iaermlc == 3</I> are quite different. In the first case, the DRE of each type
!! (alone) is calculated. In the second case, instead, the impact of each type (missing) on the total DRE is calculated.
!! 
!! Other important notes are:
!! - the DRE of a single running type is non-zero only if <I>fhbari > 0</I> for that type;
!! - if <I>iaermdl == 0</I> (OPAC routines), the single call of radiation is automatically set.
!! 
!! 
!! <br/>
!! #### DNI-DHI-GHI ####
!! 
!! For solar energy studies, 3 additional output fluxes are stored:
!! 
!! <pre>
!! DNI -> Direct Normalized Irradiance (SW): direct flux (from the Sun) collected by an unit normal surface
!! DHI -> Diffuse Horizontal Irradiance (SW): diffuse flux (from all directions) collected by an unit horizontal surface
!! GHI -> Global Horizontal Irradiance (SW): direct + diffuse fluxes collected by an unit horizontal surface </pre>
!! 
!! From these definitions it follows that: <I>GHI = DNI*cos(zen) + DHI</I> (where <I>zen</I> is the zenith angle).
!! These fluxes are computed in the <I>swrad</I> routine (<I>radsw_main_nmmb.f</I>) and stored in the output variable:
!! 
!! <pre>
!! ENG_SW_SFC(:,:,1) -> DNI (w/m²)
!! ENG_SW_SFC(:,:,2) -> DHI (w/m²)
!! ENG_SW_SFC(:,:,3) -> GHI (w/m²) </pre>
!! 
!! Currently only total-sky fluxes are stored but the clear-sky fluxes are also computed.
!! If the radiation multiple call is activated, the differences between these 3 fluxes
!! with and without aerosols are also computed and stored in the same variable:
!! 
!! <pre>
!! ENG_SW_SFC(:,:,4) -> DNI - DNI_0 (w/m²)
!! ENG_SW_SFC(:,:,5) -> DHI - DHI_0 (w/m²)
!! ENG_SW_SFC(:,:,6) -> GHI - GHI_0 (w/m²) </pre>
!! 
!! Finally, just for evaluation purpose, <I>cos(zen)</I> is stored in the last position of the array:
!! 
!! <pre>
!! ENG_SW_SFC(:,:,7) -> cos(zen) </pre>
!! 
!! 
!! <br/>
!! ### ARI CONFIGURATION ###
!! 
!! The ARI configuration requires the setup of the configfile (mainly via <I>proj_*.conf</I>) and the optical table.
!! 
!! 
!! <br/>
!! #### configfile_rrtm_chem.*.tmp ####
!! 
!! Only 4 new flags control the ARI setup in the configfile:
!! 
!! <pre>
!! slfari -> flag for ARI LW/SW calculations (old IAER)
!!   : 0 -> no ARI
!!   : 1 -> only SW
!!   : 10 -> only LW
!!   : 11 -> LW + SW
!!   : 1** -> as before + volcanic ash (always from OPAC)
!! 
!! mdlari -> flag for ARI module (XMDL in proj_*.conf)
!!   : 0 -> OPAC climatology
!!   : 3 -> MONARCH dynamic coupling
!! 
!! mlcari -> flag for multiple call (XMLC in proj_*.conf)
!!   : 0 -> no DRE
!!   : 1 -> DRE of all types
!!   : 2 -> DRE of single types (alone) + all types
!!   : 3 -> DRE of single types (missing) + all types
!! 
!! tblari -> flag for optical table (XTBL in proj_*.conf)
!!   : 1 -> OPAC refractive indices (homogeneous dust components)
!!   : 2 -> OPAC refractive indices + emission p50 for dust imaginary index in (0.25-2.00)µm
!!   : 3 -> mineral refractive indices: external mixing
!!   : 4 -> mineral refractive indices: internal mixing for hematite-clays (Maxwell-Garnett mixing rule) </pre>
!! 
!! In <I>proj_*.conf</I> another flag (<I>XKOK</I>) has been added, which allows controlling the shape for dust
!! (only for <I>tblari: 1, 2</I>):
!! 
!! <pre>
!! XKOK -> flag for shape (only for XTBL=1 or XTBL=2)
!!   : 0 -> spheres
!!   : 1 -> tri-axial spheroids for dust in (0.25-2.00)µm </pre>
!! 
!! Other parameters are required but fixed until the general aerosol parameterization is not changed:
!! 
!! <pre>
!! nspari: 10 -> number of available aerosol species (for allocation of nxbari)
!! nbnari: 90 -> number of available aerosol types (for allocation of fhbari)
!! nxbari: 1   64     8                 7               2     1   3       1   2     1 -> number of available types for each aerosol species
!!        |-| |----| |---------------| |-------------| |---| |-| |-----| |-| |---| |-|
!!        ash  dust   salt              om              bc   so4  no3    nh4 unspc aero_uci
!!        |-| |----| |---------------| |-------------| |---| |-| |-----| |-| |---| |-|
!! fhbari: 0  1(x64)  2 2 2 2 2 2 2 2   1 2 1 2 2 2 2   1 2   2   0 0 0   0   0 0   0 -> ARI configuration for each type
!!                                                                                         : 0 -> no ARI
!!                                                                                         : 1 -> hydrophobic
!!                                                                                         : 2 -> hydrophilic </pre>
!! 
!! Important notes:
!! - the order of species in <I>nxbari</I> must be the same as in the aerosol tracer array;
!! - the order of types in <I>fhbari</I> (based on <I>nxbari</I>) must be the same as in the optical tables.
!! 
!! 
!! <br/>
!! #### optical tables ####
!! 
!! Currently 4 main optical tables are available (+ 2 accessory tables with tri-axial spheroid optical properties for dust),
!! prepared for 5 global aerosols: dust, salt, om, bc and so4 from Michele's thesis (the other MONARCH species are not optically
!! parameterized yet and so are included in the tables as "phantom" types):
!! 
!! <pre>
!! Main tables:
!!   mnch_tables_opc.txt (tblari: 1) -> opc stands for OPAC
!!   mnch_tables_clq.txt (tblari: 2) -> clq stands for Claquin
!!   mnch_tables_ext.txt (tblari: 3) -> ext stands for external mixing
!!   mnch_tables_int.txt (tblari: 4) -> int stands for internal mixing
!! 
!! Accessory tables:
!!   mnch_tables_opc_kok.txt (tblari: 1 + XKOK=1)
!!   mnch_tables_clq_kok.txt (tblari: 2 + XKOK=1) </pre>
!! 
!! In the main tables, only the refractive indices for dust change from one table to another one, while the accessory tables
!! only affect the shape for dust (the parameterization for the other species is unvaried and their optical properties are
!! repeated in all tables). All the tables have the following structure (<I>b</I> -> size bin; <I>m</I> -> mode):
!! 
!! <pre>
!!     1      m1        ash  ("phantom" type)
!!   2-9   b1-b8       dust     (tblari: 1,2)      or illite            (tblari: 3)      or hem/clays mixture (tblari: 4)
!! 10-17   b1-b8       dust     (tblari: 1,2)      or kaolinite         (tblari: 3)      or hem/clays mixture (tblari: 4)
!! 18-25   b1-b8       dust     (tblari: 1,2)      or montmorillonite   (tblari: 3)      or hem/clays mixture (tblari: 4)
!! 26-33   b1-b8       dust     (tblari: 1,2)      or calcite         (tblari: 3,4)
!! 34-41   b1-b8       dust     (tblari: 1,2)      or quartz          (tblari: 3,4)
!! 42-49   b1-b8       dust     (tblari: 1,2)      or feldspar        (tblari: 3,4)
!! 50-57   b1-b8       dust     (tblari: 1,2)      or hematite          (tblari: 3)      or hem/clays mixture (tblari: 4)
!! 58-65   b1-b8       dust     (tblari: 1,2)      or gypsum          (tblari: 3,4)
!! 66-73   b1-b8       salt (tblari: 1,2,3,4)
!! 74-80   m1-m7         om (tblari: 1,2,3,4)
!! 81-82   m1-m2         bc (tblari: 1,2,3,4)
!!    83      m1        so4 (tblari: 1,2,3,4)
!! 84-86   m1-m3        no3  ("phantom" type)
!!    87      m1        nh4  ("phantom" type)
!! 88-89   m1-m2      unspc  ("phantom" type)
!!    90      m1   aero_uci  ("phantom" type) </pre>
!! 
!! The position in the tables of each type (first column) is the id number contained in <I>setari(:,1)</I> and identifies
!! the subtable (61 wavelengths: 7 RH levels: 3 optical properties) relative to that type.
!! 
!! 
!! <br/>
!! #### optical table calculation ####
!! 
!! The scripts for the calculation of the optical tables are in the Gitlab project:
!! 
!! <pre>
!! https://earth.bsc.es/gitlab/ac/monarch-radiative-tables </pre>
!! 
!! The uppercase directories contain input data:
!! 
!! <pre>
!! OPAC/ (used for all tables) -> refractive indices from OPAC database:
!!   miam00 <-> dust
!!   ssam00 <-> salt
!!   waso00 <-> om
!!   soot00 <-> bc
!!   suso00 <-> so4
!!   inso00 <-> "phantom" type
!! 
!! MCPH/ (used for all tables) -> microphysical properties for all types (line <-> type):
!!   growth_factors.txt
!!   size_distributions.txt
!!   opac_files.txt
!!   mass_densities.txt
!!   segelstein.txt -> water refractive indices [Segelstein, 1981]
!! 
!! SCNZ/ (used for mnch_tables_clq.txt) -> mineral refractive indices [Scanza et al., 2015]:
!!   illite.txt
!!   kaolinite.txt
!!   montmorillonite.txt
!!   calcite.txt
!!   quartz.txt
!!   feldspar.txt
!!   hematite.txt
!!   gypsum.txt
!! 
!! OXFD/ (used for mnch_tables_ext/int.txt) -> modified mineral refractive indices (Oxford University ARIA):
!!   calcite.txt -> (3.0-15.0)µm (E-ray + O-ray) [Posch et al., 2007]
!!   quartz.txt -> (7.2-15.0)µm [Kitamura et al., 2007]
!! 
!! VLFR/ (used for mnch_tables_clq.txt) -> mineral volume fractions for 28 soil types and 8 bins:
!!   volume_fractions_* -> from emission dataset [Claquin et al., 1999] </pre>
!! 
!! The lowercase directories contain the calculation scripts for the tables:
!! 
!! <pre>
!! opc/ clq/ ext/ int/:
!!   input_*.ipy -> creates input text file (input_*.txt) for spher_sed (Mie code)
!!   table_*.bash -> reads input_*.txt, runs spher_sed and generates mnch_tables_*.txt
!!   sbatch_*.cmd -> launches table_*.bash in Marenostrum4
!! 
!! kok/:
!!   substitute.py -> substitutes optical properties for tri-axial spheroids [Kok et al., 2017] in mnch_tables_opc/clq.txt </pre>
!! 
!! For parallel computing in Marenostrum4 (for main tables) it is necessary to use <I>input_*.txt</I>, <I>table_*.bash</I>,
!! <I>sbatch_*.cmd</I> and <I>spher_sed</I> in different directories. In case of internal mixing, this must be done separately
!! for each hematite volume ratio (<I>v00</I>, <I>v02</I>, ..., <I>v20</I>) and then running <I>table_int.ipy</I> generates
!! <I>mnch_tables_int.txt</I>.
!! 
!! Once calculated <I>mnch_tables_*.txt</I>, these must be copied within the <I>tables2018</I> directory in the 2 HPC
!! repositories for MONARCH:
!! 
!! <pre>
!! /gpfs/archive/bsc32/esarchive/exp/monarch/static/fix/fix_aeropt_bsc/
!! /gpfs/projects/bsc32/models/monarch/data/static/fix/fix_aeropt_bsc/ </pre>
!! 
!! A specific instruction in <I>05b_sim.sh</I> creates virtual links to all these tables in <I>CURRENT_RUN</I>
!! (<I>XKOK</I> controls the virtual link creation for mnch_tables_opc/clq.txt).
!! 
!! Note that in the HPC repositories, other directories (<I>tables_oldfmt_chinhg</I>, <I>tables_oldfmt_opachg</I>
!! and <I>tables_newfmt_basepk</I>) contain no longer used optical tables (described in a local <I>README</I>).

!*******************************************************************************