diff --git a/dependencies.yaml b/dependencies.yaml index c1313e4d75..3a68a0d26d 100644 --- a/dependencies.yaml +++ b/dependencies.yaml @@ -20,6 +20,8 @@ casim: source: git@github.com:MetOffice/casim.git ref: 2026.07.1 + source: git@github.com:iboutle/casim.git + ref: activation_fix jules: source: git@github.com:MetOffice/jules.git diff --git a/interfaces/physics_schemes_interface/rose-meta/um-microphysics/HEAD/rose-meta.conf b/interfaces/physics_schemes_interface/rose-meta/um-microphysics/HEAD/rose-meta.conf index 7c3945f903..bb88181c0a 100644 --- a/interfaces/physics_schemes_interface/rose-meta/um-microphysics/HEAD/rose-meta.conf +++ b/interfaces/physics_schemes_interface/rose-meta/um-microphysics/HEAD/rose-meta.conf @@ -80,6 +80,25 @@ range=-1:1 sort-key=Panel-A05 type=real +[namelist:microphysics=casim_activation] +compulsory=true +description=CASIM droplet activation option +!enumeration=true +help=Fixed number uses a fixed cloud droplet number within liquid cloud and + =makes no use of any aerosol, which is the option used by RAL3. + = + =Abdul-Razzak and Ghan uses the mechanistic activation scheme of + =Abdul-Razzak and Ghan (2000, https://doi.org/10.1029/1999JD901161) to + =work the cloud droplet number out from the fraction activated in each + =aerosol mode. The aerosol is taken from the GLOMAP modal aerosol, which + =is passed into CASIM but is not modified by it, so UKCA must be running. +!kind=default +ns=namelist/Science/UM Microphysics +sort-key=Panel-A02c +trigger=namelist:microphysics=casim_cdnc_opt: 'fixed'; +value-titles=Fixed number, Abdul-Razzak and Ghan +values='fixed','arg' + [namelist:microphysics=casim_cdnc_opt] compulsory=true description=How to calculate cloud number in Casim @@ -89,9 +108,13 @@ help=External will take the cloud number calculated by an external aerosol = with the Jones relation = =Fixed number will use a fixed cloud drop number, as in RA3 + = + =This only applies to the fixed number droplet activation. Abdul-Razzak + =and Ghan activation always works the cloud number out from the aerosol + =supplied to CASIM, so there is nothing external left for it to take. !kind=default ns=namelist/Science/UM Microphysics -sort-key=Panel-A02 +sort-key=Panel-A02d value-titles=External code, Fixed number values='external','fixed' @@ -295,7 +318,7 @@ help=Setting to True will switch from the Wilson and Ballard (1999) !kind=default ns=namelist/Science/UM Microphysics sort-key=Panel-A00 -trigger=namelist:microphysics=casim_cdnc_opt: .true.; +trigger=namelist:microphysics=casim_activation: .true.; =namelist:microphysics=graupel_scheme: .false.; =namelist:microphysics=shape_rime: .false.; =namelist:microphysics=orog_rain: .false.; diff --git a/interfaces/physics_schemes_interface/rose-meta/um-microphysics/versions.py b/interfaces/physics_schemes_interface/rose-meta/um-microphysics/versions.py index 01e65a2b91..5c6fa8229b 100644 --- a/interfaces/physics_schemes_interface/rose-meta/um-microphysics/versions.py +++ b/interfaces/physics_schemes_interface/rose-meta/um-microphysics/versions.py @@ -31,3 +31,17 @@ def upgrade(self, config, meta_config=None): # Add settings return config, self.reports """ + + +class vn32_t709(MacroUpgrade): + """Upgrade macro for #709 by Ian Boutle.""" + + BEFORE_TAG = "vn3.2" + AFTER_TAG = "vn3.2_t709" + + def upgrade(self, config, meta_config=None): + # Commands From: rose-meta/um-microphysics + self.add_setting( + config, ["namelist:microphysics", "casim_activation"], "'fixed'" + ) + return config, self.reports diff --git a/interfaces/physics_schemes_interface/source/algorithm/casim_activate_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/casim_activate_alg_mod.x90 index b20232fa10..7fd32c9d82 100644 --- a/interfaces/physics_schemes_interface/source/algorithm/casim_activate_alg_mod.x90 +++ b/interfaces/physics_schemes_interface/source/algorithm/casim_activate_alg_mod.x90 @@ -3,6 +3,9 @@ ! The file LICENCE, distributed with this code, contains details of the terms ! under which the code may be used. !------------------------------------------------------------------------------- +! Some of the content of this file has been produced with the assistance of +! Anthropic Claude Opus 5 (Claude Code). +!------------------------------------------------------------------------------- !> @brief Interface to the CASIM activation scheme module casim_activate_alg_mod @@ -17,6 +20,7 @@ module casim_activate_alg_mod use mesh_mod, only: mesh_type use microphysics_config_mod, only: casim_cdnc_opt, & casim_cdnc_opt_fixed + use mphys_inputs_mod, only: casim_iopt_act, fixed_number use mr_indices_mod, only: nummr, imr_cl, imr_ci use sci_geometric_constants_mod, only: get_height_fv use section_choice_config_mod, only: convection, convection_um @@ -37,18 +41,25 @@ contains !> documentation. !>@param[in] config Application namelist configuration object !>@param[in] mr (in) Mixing ratios, in theta space + !>@param[in] m_cl_pre_fast (in) Cloud liq after advection !>@param[in] cloud_fields (in) Fields for cloud fractions !>@param[in,out] microphysics_fields (in,out) Fields for mphys scheme + !>@param[in] convection_fields (in) Fields for convection scheme + !>@param[in] aerosol_fields (in) Fields for the aerosol schemes + !>@param[in] turbulence_fields (in) Fields for turbulence scheme !>@param[in] initialise Flag to use routine to initialise fields subroutine casim_activate_alg(config, & - theta, mr, & + theta, mr, m_cl_pre_fast, & derived_fields, & cloud_fields, & microphysics_fields, & - convection_fields, initialise) + convection_fields, & + aerosol_fields, turbulence_fields, & + initialise) use casim_activate_kernel_mod, only: casim_activate_kernel_type use casim_ice_act_kernel_mod, only: casim_ice_act_kernel_type + use casim_arg_act_kernel_mod, only: casim_arg_act_kernel_type implicit none @@ -57,10 +68,12 @@ contains type( field_collection_type ), intent(in) :: derived_fields type( field_collection_type ), intent(in) :: cloud_fields type( field_collection_type ), intent(in) :: convection_fields + type( field_collection_type ), intent(in) :: aerosol_fields + type( field_collection_type ), intent(in) :: turbulence_fields type( field_collection_type ), intent(inout) :: microphysics_fields type( field_type ), intent( in ), target :: mr(nummr) - type( field_type ), intent( in ) :: theta + type( field_type ), intent( in ) :: theta, m_cl_pre_fast logical, intent(in) :: initialise @@ -73,13 +86,42 @@ contains type( field_type ), pointer :: rho_wth type( field_type ), pointer :: mci_init + ! Fields needed by the mechanistic activation + type( field_type ), pointer :: cf_liq_pre_fast + type( field_type ), pointer :: w_in_wth + type( field_type ), pointer :: wvar + + ! GLOMAP modal aerosol + type( field_type ), pointer :: n_ait_sol + type( field_type ), pointer :: ait_sol_su + type( field_type ), pointer :: ait_sol_bc + type( field_type ), pointer :: ait_sol_om + type( field_type ), pointer :: n_acc_sol + type( field_type ), pointer :: acc_sol_su + type( field_type ), pointer :: acc_sol_bc + type( field_type ), pointer :: acc_sol_om + type( field_type ), pointer :: acc_sol_ss + type( field_type ), pointer :: n_cor_sol + type( field_type ), pointer :: cor_sol_su + type( field_type ), pointer :: cor_sol_bc + type( field_type ), pointer :: cor_sol_om + type( field_type ), pointer :: cor_sol_ss + type( field_type ), pointer :: n_ait_ins + type( field_type ), pointer :: ait_ins_bc + type( field_type ), pointer :: ait_ins_om + type( field_type ), pointer :: n_acc_ins + type( field_type ), pointer :: acc_ins_du + type( field_type ), pointer :: n_cor_ins + type( field_type ), pointer :: cor_ins_du + type( mesh_type ), pointer :: mesh integer( tik ) :: id if ( LPROF ) call start_timing( id, 'microphysics.casim_activate' ) ! For fixed CDNC or if we want to initialise the CDNC - if (casim_cdnc_opt == casim_cdnc_opt_fixed .or. initialise) then + if ( ( casim_cdnc_opt == casim_cdnc_opt_fixed .and. & + casim_iopt_act == fixed_number ) .or. initialise ) then call cloud_fields%get_field('liquid_fraction', cf_liq) call microphysics_fields%get_field('nl_mphys', nl_mphys) @@ -97,6 +139,64 @@ contains end if + ! Mechanistic activation of any liquid cloud which the physics outside + ! CASIM has created since the start of the fast physics + if (casim_iopt_act > fixed_number .and. .not. initialise) then + + call cloud_fields%get_field('liquid_fraction', cf_liq) + call microphysics_fields%get_field('nl_mphys', nl_mphys) + call microphysics_fields%get_field('cf_liq_pre_fast', cf_liq_pre_fast) + call derived_fields%get_field('exner_in_wth', exner_wth) + call derived_fields%get_field('rho_in_wth', rho_wth) + call derived_fields%get_field('w_in_wth', w_in_wth) + call turbulence_fields%get_field('wvar', wvar) + + call aerosol_fields%get_field('n_ait_sol', n_ait_sol) + call aerosol_fields%get_field('ait_sol_su', ait_sol_su) + call aerosol_fields%get_field('ait_sol_bc', ait_sol_bc) + call aerosol_fields%get_field('ait_sol_om', ait_sol_om) + call aerosol_fields%get_field('n_acc_sol', n_acc_sol) + call aerosol_fields%get_field('acc_sol_su', acc_sol_su) + call aerosol_fields%get_field('acc_sol_bc', acc_sol_bc) + call aerosol_fields%get_field('acc_sol_om', acc_sol_om) + call aerosol_fields%get_field('acc_sol_ss', acc_sol_ss) + call aerosol_fields%get_field('n_cor_sol', n_cor_sol) + call aerosol_fields%get_field('cor_sol_su', cor_sol_su) + call aerosol_fields%get_field('cor_sol_bc', cor_sol_bc) + call aerosol_fields%get_field('cor_sol_om', cor_sol_om) + call aerosol_fields%get_field('cor_sol_ss', cor_sol_ss) + call aerosol_fields%get_field('n_ait_ins', n_ait_ins) + call aerosol_fields%get_field('ait_ins_bc', ait_ins_bc) + call aerosol_fields%get_field('ait_ins_om', ait_ins_om) + call aerosol_fields%get_field('n_acc_ins', n_acc_ins) + call aerosol_fields%get_field('acc_ins_du', acc_ins_du) + call aerosol_fields%get_field('n_cor_ins', n_cor_ins) + call aerosol_fields%get_field('cor_ins_du', cor_ins_du) + + call invoke( casim_arg_act_kernel_type( mr(imr_cl), m_cl_pre_fast, & + rho_wth, theta, exner_wth, & + cf_liq, cf_liq_pre_fast, & + w_in_wth, wvar, & + nl_mphys, & + n_ait_sol, ait_sol_su, & + ait_sol_bc, ait_sol_om, & + n_acc_sol, acc_sol_su, & + acc_sol_bc, acc_sol_om, & + acc_sol_ss, & + n_cor_sol, cor_sol_su, & + cor_sol_bc, cor_sol_om, & + cor_sol_ss, & + n_ait_ins, ait_ins_bc, & + ait_ins_om, & + n_acc_ins, acc_ins_du, & + n_cor_ins, cor_ins_du ) ) + + if (write_diag .and. use_xios_io) then + call nl_mphys%write_field('casim__nl_mphys_act') + end if + + end if + ! If convection has provided an increment to ice or we want to ! initialise the ice if (convection == convection_um .or. initialise) then diff --git a/interfaces/physics_schemes_interface/source/algorithm/casim_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/casim_alg_mod.x90 index 51e1ec17ba..afe7009887 100644 --- a/interfaces/physics_schemes_interface/source/algorithm/casim_alg_mod.x90 +++ b/interfaces/physics_schemes_interface/source/algorithm/casim_alg_mod.x90 @@ -3,6 +3,9 @@ ! The file LICENCE, distributed with this code, contains details of the terms ! under which the code may be used. !------------------------------------------------------------------------------- +! Some of the content of this file has been produced with the assistance of +! Anthropic Claude Opus 5 (Claude Code). +!------------------------------------------------------------------------------- !> @brief Interface to the CASIM Microphysics scheme module casim_alg_mod @@ -115,6 +118,29 @@ contains type( field_type ), pointer :: wvar + ! GLOMAP modal aerosol, used by the CASIM mechanistic activation + type( field_type ), pointer :: n_ait_sol + type( field_type ), pointer :: ait_sol_su + type( field_type ), pointer :: ait_sol_bc + type( field_type ), pointer :: ait_sol_om + type( field_type ), pointer :: n_acc_sol + type( field_type ), pointer :: acc_sol_su + type( field_type ), pointer :: acc_sol_bc + type( field_type ), pointer :: acc_sol_om + type( field_type ), pointer :: acc_sol_ss + type( field_type ), pointer :: n_cor_sol + type( field_type ), pointer :: cor_sol_su + type( field_type ), pointer :: cor_sol_bc + type( field_type ), pointer :: cor_sol_om + type( field_type ), pointer :: cor_sol_ss + type( field_type ), pointer :: n_ait_ins + type( field_type ), pointer :: ait_ins_bc + type( field_type ), pointer :: ait_ins_om + type( field_type ), pointer :: n_acc_ins + type( field_type ), pointer :: acc_ins_du + type( field_type ), pointer :: n_cor_ins + type( field_type ), pointer :: cor_ins_du + type( field_type ), pointer :: height_w3 => null() type( field_type ), pointer :: height_wth => null() type( field_type ), pointer :: dz_in_wth @@ -170,6 +196,28 @@ contains call aerosol_fields%get_field('cloud_drop_no_conc',cloud_drop_no_conc) call aerosol_fields%get_field('murk', murk) + call aerosol_fields%get_field('n_ait_sol', n_ait_sol) + call aerosol_fields%get_field('ait_sol_su', ait_sol_su) + call aerosol_fields%get_field('ait_sol_bc', ait_sol_bc) + call aerosol_fields%get_field('ait_sol_om', ait_sol_om) + call aerosol_fields%get_field('n_acc_sol', n_acc_sol) + call aerosol_fields%get_field('acc_sol_su', acc_sol_su) + call aerosol_fields%get_field('acc_sol_bc', acc_sol_bc) + call aerosol_fields%get_field('acc_sol_om', acc_sol_om) + call aerosol_fields%get_field('acc_sol_ss', acc_sol_ss) + call aerosol_fields%get_field('n_cor_sol', n_cor_sol) + call aerosol_fields%get_field('cor_sol_su', cor_sol_su) + call aerosol_fields%get_field('cor_sol_bc', cor_sol_bc) + call aerosol_fields%get_field('cor_sol_om', cor_sol_om) + call aerosol_fields%get_field('cor_sol_ss', cor_sol_ss) + call aerosol_fields%get_field('n_ait_ins', n_ait_ins) + call aerosol_fields%get_field('ait_ins_bc', ait_ins_bc) + call aerosol_fields%get_field('ait_ins_om', ait_ins_om) + call aerosol_fields%get_field('n_acc_ins', n_acc_ins) + call aerosol_fields%get_field('acc_ins_du', acc_ins_du) + call aerosol_fields%get_field('n_cor_ins', n_cor_ins) + call aerosol_fields%get_field('cor_ins_du', cor_ins_du) + call turbulence_fields%get_field('wvar', wvar) height_w3 => get_height_fv(config, theta%get_mesh(), W3) @@ -195,7 +243,16 @@ contains ls_rain_3d, ls_snow_3d, ls_graup_3d, & theta_inc, cloud_drop_no_conc, murk, & refl_tot, refl_1km, & - superc_liq, superc_rain ) ) + superc_liq, superc_rain, & + n_ait_sol, ait_sol_su, & + ait_sol_bc, ait_sol_om, & + n_acc_sol, acc_sol_su, & + acc_sol_bc, acc_sol_om, acc_sol_ss, & + n_cor_sol, cor_sol_su, & + cor_sol_bc, cor_sol_om, cor_sol_ss, & + n_ait_ins, ait_ins_bc, ait_ins_om, & + n_acc_ins, acc_ins_du, & + n_cor_ins, cor_ins_du ) ) if (turb_gen_mixph) then dz_in_wth => get_dz_at_wtheta(config, theta%get_mesh()) diff --git a/interfaces/physics_schemes_interface/source/kernel/casim_arg_act_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/casim_arg_act_kernel_mod.F90 new file mode 100644 index 0000000000..f92b6810fe --- /dev/null +++ b/interfaces/physics_schemes_interface/source/kernel/casim_arg_act_kernel_mod.F90 @@ -0,0 +1,273 @@ +!----------------------------------------------------------------------------- +! (c) Crown copyright 2026 Met Office. All rights reserved. +! The file LICENCE, distributed with this code, contains details of the terms +! under which the code may be used. +!----------------------------------------------------------------------------- +! Some of the content of this file has been produced with the assistance of +! Anthropic Claude Opus 5 (Claude Code). +!----------------------------------------------------------------------------- +!> @brief Mechanistic activation of the CASIM cloud droplet number. + +module casim_arg_act_kernel_mod + +use argument_mod, only: arg_type, & + GH_FIELD, GH_REAL, & + GH_READ, & + GH_READWRITE, & + CELL_COLUMN +use fs_continuity_mod, only: WTHETA +use kernel_mod, only: kernel_type + +implicit none + +private + +!------------------------------------------------------------------------------- +! Public types +!------------------------------------------------------------------------------- +!> The type declaration for the kernel. +!> Contains the metadata needed by the Psy layer + +type, public, extends(kernel_type) :: casim_arg_act_kernel_type + private + type(arg_type) :: meta_args(31) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! m_cl + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! m_cl_pre_fast + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! rho_in_wth + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! theta + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! exner_in_wth + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! cf_liq + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! cf_liq_pre_fast + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! w_in_wth + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! wvar + arg_type(GH_FIELD, GH_REAL, GH_READWRITE, WTHETA), & ! nl_mphys + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! n_ait_sol + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! ait_sol_su + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! ait_sol_bc + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! ait_sol_om + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! n_acc_sol + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! acc_sol_su + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! acc_sol_bc + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! acc_sol_om + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! acc_sol_ss + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! n_cor_sol + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! cor_sol_su + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! cor_sol_bc + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! cor_sol_om + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! cor_sol_ss + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! n_ait_ins + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! ait_ins_bc + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! ait_ins_om + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! n_acc_ins + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! acc_ins_du + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! n_cor_ins + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA) & ! cor_ins_du + /) + integer :: operates_on = CELL_COLUMN +contains + procedure, nopass :: casim_arg_act_code +end type + +public :: casim_arg_act_code + +contains + +!> @brief Mechanistic activation of the CASIM cloud droplet number +!> @details Liquid cloud can be created or removed by processes outside the +!! CASIM microphysics. For the mechanistic activation options, the +!! cloud droplet number for any newly created cloud is worked out +!! from the GLOMAP modal aerosol using the CASIM activation scheme, +!! and the number is removed again where the cloud has gone. +!> @param[in] nlayers Number of layers +!> @param[in] m_cl Liquid cloud mass mixing ratio +!> @param[in] m_cl_pre_fast Liquid cloud mass mixing ratio as it was on +!! entry to the fast physics +!> @param[in] rho_in_wth Dry density in potential temperature space +!> @param[in] theta Potential temperature field +!> @param[in] exner_in_wth Exner pressure in potential temperature space +!> @param[in] cf_liq Liquid cloud fraction +!> @param[in] cf_liq_pre_fast Liquid cloud fraction as it was on entry to +!! the fast physics +!> @param[in] w_in_wth 'Vertical' wind in theta space +!> @param[in] wvar Vertical velocity variance +!> @param[in,out] nl_mphys CASIM cloud-droplet number concentration +!> @param[in] n_ait_sol Soluble Aitken mode number mixing ratio +!> @param[in] ait_sol_su Soluble Aitken mode H2SO4 mass mixing ratio +!> @param[in] ait_sol_bc Soluble Aitken mode black carbon m.m.r. +!> @param[in] ait_sol_om Soluble Aitken mode organic m.m.r. +!> @param[in] n_acc_sol Soluble accumulation mode number m.r. +!> @param[in] acc_sol_su Soluble accumulation mode H2SO4 m.m.r. +!> @param[in] acc_sol_bc Soluble accumulation mode black carbon m.m.r. +!> @param[in] acc_sol_om Soluble accumulation mode organic m.m.r. +!> @param[in] acc_sol_ss Soluble accumulation mode sea salt m.m.r. +!> @param[in] n_cor_sol Soluble coarse mode number mixing ratio +!> @param[in] cor_sol_su Soluble coarse mode H2SO4 mass mixing ratio +!> @param[in] cor_sol_bc Soluble coarse mode black carbon m.m.r. +!> @param[in] cor_sol_om Soluble coarse mode organic m.m.r. +!> @param[in] cor_sol_ss Soluble coarse mode sea salt m.m.r. +!> @param[in] n_ait_ins Insoluble Aitken mode number mixing ratio +!> @param[in] ait_ins_bc Insoluble Aitken mode black carbon m.m.r. +!> @param[in] ait_ins_om Insoluble Aitken mode organic m.m.r. +!> @param[in] n_acc_ins Insoluble accumulation mode number m.r. +!> @param[in] acc_ins_du Insoluble accumulation mode dust m.m.r. +!> @param[in] n_cor_ins Insoluble coarse mode number mixing ratio +!> @param[in] cor_ins_du Insoluble coarse mode dust m.m.r. +!> @param[in] ndf_wth Number of degrees of freedom per cell for +!! potential temperature space +!> @param[in] undf_wth Number unique of degrees of freedom for +!! potential temperature space +!> @param[in] map_wth Dofmap for the cell at the base of the +!! column for potential temperature space + +subroutine casim_arg_act_code( nlayers, & + m_cl, m_cl_pre_fast, & + rho_in_wth, theta, exner_in_wth, & + cf_liq, cf_liq_pre_fast, & + w_in_wth, wvar, & + nl_mphys, & + n_ait_sol, ait_sol_su, & + ait_sol_bc, ait_sol_om, & + n_acc_sol, acc_sol_su, & + acc_sol_bc, acc_sol_om, & + acc_sol_ss, & + n_cor_sol, cor_sol_su, & + cor_sol_bc, cor_sol_om, & + cor_sol_ss, & + n_ait_ins, ait_ins_bc, & + ait_ins_om, & + n_acc_ins, acc_ins_du, & + n_cor_ins, cor_ins_du, & + ndf_wth, undf_wth, map_wth ) + + use constants_mod, only: r_def, i_def, r_um, i_um + + !--------------------------------------- + ! UM modules + !--------------------------------------- + + use planet_constants_mod, only: p_zero, kappa + use nlsizes_namelist_mod, only: tr_ukca + use mphys_inputs_mod, only: wvarfac + use casim_activation_in_um_mod, only: activate_column_ukca + use casim_ukca_tracer_mod, only: casim_ukca_tracer_column + + implicit none + + ! Arguments + integer(kind=i_def), intent(in) :: nlayers + integer(kind=i_def), intent(in) :: ndf_wth + integer(kind=i_def), intent(in) :: undf_wth + + real(kind=r_def), intent(in), dimension(undf_wth) :: m_cl + real(kind=r_def), intent(in), dimension(undf_wth) :: m_cl_pre_fast + real(kind=r_def), intent(in), dimension(undf_wth) :: rho_in_wth + real(kind=r_def), intent(in), dimension(undf_wth) :: theta + real(kind=r_def), intent(in), dimension(undf_wth) :: exner_in_wth + real(kind=r_def), intent(in), dimension(undf_wth) :: cf_liq + real(kind=r_def), intent(in), dimension(undf_wth) :: cf_liq_pre_fast + real(kind=r_def), intent(in), dimension(undf_wth) :: w_in_wth + real(kind=r_def), intent(in), dimension(undf_wth) :: wvar + + real(kind=r_def), intent(inout), dimension(undf_wth) :: nl_mphys + + real(kind=r_def), intent(in), dimension(undf_wth) :: n_ait_sol + real(kind=r_def), intent(in), dimension(undf_wth) :: ait_sol_su + real(kind=r_def), intent(in), dimension(undf_wth) :: ait_sol_bc + real(kind=r_def), intent(in), dimension(undf_wth) :: ait_sol_om + real(kind=r_def), intent(in), dimension(undf_wth) :: n_acc_sol + real(kind=r_def), intent(in), dimension(undf_wth) :: acc_sol_su + real(kind=r_def), intent(in), dimension(undf_wth) :: acc_sol_bc + real(kind=r_def), intent(in), dimension(undf_wth) :: acc_sol_om + real(kind=r_def), intent(in), dimension(undf_wth) :: acc_sol_ss + real(kind=r_def), intent(in), dimension(undf_wth) :: n_cor_sol + real(kind=r_def), intent(in), dimension(undf_wth) :: cor_sol_su + real(kind=r_def), intent(in), dimension(undf_wth) :: cor_sol_bc + real(kind=r_def), intent(in), dimension(undf_wth) :: cor_sol_om + real(kind=r_def), intent(in), dimension(undf_wth) :: cor_sol_ss + real(kind=r_def), intent(in), dimension(undf_wth) :: n_ait_ins + real(kind=r_def), intent(in), dimension(undf_wth) :: ait_ins_bc + real(kind=r_def), intent(in), dimension(undf_wth) :: ait_ins_om + real(kind=r_def), intent(in), dimension(undf_wth) :: n_acc_ins + real(kind=r_def), intent(in), dimension(undf_wth) :: acc_ins_du + real(kind=r_def), intent(in), dimension(undf_wth) :: n_cor_ins + real(kind=r_def), intent(in), dimension(undf_wth) :: cor_ins_du + + integer(kind=i_def), intent(in), dimension(ndf_wth) :: map_wth + + ! Local variables for the kernel + ! Bounds on the activating vertical velocity, as used by the UM + real(r_um), parameter :: min_velocity = 0.01_r_um ! ms-1 + real(r_um), parameter :: max_velocity = 4.00_r_um ! ms-1 + real(r_um), parameter :: eps = epsilon(1.0_r_um) + + real(r_um), dimension(nlayers) :: & + cloud_mass_post_qtbal, cloud_mass_pre_ap2, rho_col, t_col, p_col, & + cf_liquid, cf_liquid_pre_ap2, w_tke, cloudnumber_col + + ! UKCA tracer array holding the GLOMAP modes for this column + real(r_um), allocatable :: tracer_ukca(:,:,:,:) + + integer(i_um) :: k + + !------------------------------------------------------------------------- + ! End of Declarations + !------------------------------------------------------------------------- + + allocate( tracer_ukca(1,1,0:nlayers,tr_ukca) ) + + call casim_ukca_tracer_column( nlayers, undf_wth, map_wth(1), & + n_ait_sol, ait_sol_su, ait_sol_bc, & + ait_sol_om, & + n_acc_sol, acc_sol_su, acc_sol_bc, & + acc_sol_om, acc_sol_ss, & + n_cor_sol, cor_sol_su, cor_sol_bc, & + cor_sol_om, cor_sol_ss, & + n_ait_ins, ait_ins_bc, ait_ins_om, & + n_acc_ins, acc_ins_du, & + n_cor_ins, cor_ins_du, & + tracer_ukca ) + + do k = 1, nlayers + + cloud_mass_post_qtbal(k) = m_cl(map_wth(1) + k) + cloud_mass_pre_ap2(k) = m_cl_pre_fast(map_wth(1) + k) + rho_col(k) = rho_in_wth(map_wth(1) + k) + t_col(k) = exner_in_wth(map_wth(1) + k) * & + theta(map_wth(1) + k) + p_col(k) = p_zero * & + ( exner_in_wth(map_wth(1) + k) )**(1.0_r_um / kappa) + cf_liquid(k) = cf_liq(map_wth(1) + k) + cf_liquid_pre_ap2(k) = cf_liq_pre_fast(map_wth(1) + k) + cloudnumber_col(k) = nl_mphys(map_wth(1) + k) + + ! Activating vertical velocity, following the UM + if ( wvar(map_wth(1) + k) > eps ) then + w_tke(k) = w_in_wth(map_wth(1) + k) + & + wvarfac * sqrt( wvar(map_wth(1) + k) ) + else + w_tke(k) = w_in_wth(map_wth(1) + k) + end if + + w_tke(k) = min( max_velocity, max( min_velocity, w_tke(k) ) ) + + end do + + call activate_column_ukca( cloud_mass_post_qtbal, & + cloud_mass_pre_ap2, & + rho_col, t_col, p_col, & + cf_liquid, cf_liquid_pre_ap2, w_tke, & + cloudnumber_col, & + tracer_ukca ) + + do k = 1, nlayers + nl_mphys(map_wth(1) + k) = cloudnumber_col(k) + end do + + ! Set level 0 the same as level 1 (as done in the UM) + nl_mphys(map_wth(1)) = nl_mphys(map_wth(1) + 1) + + deallocate( tracer_ukca ) + +end subroutine casim_arg_act_code + +end module casim_arg_act_kernel_mod diff --git a/interfaces/physics_schemes_interface/source/kernel/casim_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/casim_kernel_mod.F90 index 9171f19ee5..aa9f45d6af 100644 --- a/interfaces/physics_schemes_interface/source/kernel/casim_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/casim_kernel_mod.F90 @@ -3,6 +3,9 @@ ! The file LICENCE, distributed with this code, contains details of the terms ! under which the code may be used. !----------------------------------------------------------------------------- +! Some of the content of this file has been produced with the assistance of +! Anthropic Claude Opus 5 (Claude Code). +!----------------------------------------------------------------------------- !> @brief Interface to CASIM microphysics scheme. module casim_kernel_mod @@ -33,7 +36,7 @@ module casim_kernel_mod type, public, extends(kernel_type) :: casim_kernel_type private - type(arg_type) :: meta_args(40) = (/ & + type(arg_type) :: meta_args(61) = (/ & arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! mv_wth arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! ml_wth arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! mi_wth @@ -73,7 +76,28 @@ module casim_kernel_mod arg_type(GH_FIELD, GH_REAL, GH_WRITE, WTHETA), & ! refl_tot arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_1), & ! refl_1km arg_type(GH_FIELD, GH_REAL, GH_WRITE, WTHETA), & ! superc_liq - arg_type(GH_FIELD, GH_REAL, GH_WRITE, WTHETA) & ! superc_rain + arg_type(GH_FIELD, GH_REAL, GH_WRITE, WTHETA), & ! superc_rain + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! n_ait_sol + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! ait_sol_su + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! ait_sol_bc + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! ait_sol_om + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! n_acc_sol + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! acc_sol_su + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! acc_sol_bc + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! acc_sol_om + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! acc_sol_ss + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! n_cor_sol + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! cor_sol_su + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! cor_sol_bc + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! cor_sol_om + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! cor_sol_ss + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! n_ait_ins + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! ait_ins_bc + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! ait_ins_om + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! n_acc_ins + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! acc_ins_du + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & ! n_cor_ins + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA) & ! cor_ins_du /) integer :: operates_on = CELL_COLUMN contains @@ -134,6 +158,27 @@ module casim_kernel_mod !! surface !> @param[in,out] superc_liq Supercooled liquid cloud mass mixing ratio !> @param[in,out] superc_rain Supercooled rain mass mixing ratio +!> @param[in] n_ait_sol Soluble Aitken mode number mixing ratio +!> @param[in] ait_sol_su Soluble Aitken mode H2SO4 mass mixing ratio +!> @param[in] ait_sol_bc Soluble Aitken mode black carbon m.m.r. +!> @param[in] ait_sol_om Soluble Aitken mode organic m.m.r. +!> @param[in] n_acc_sol Soluble accumulation mode number m.r. +!> @param[in] acc_sol_su Soluble accumulation mode H2SO4 m.m.r. +!> @param[in] acc_sol_bc Soluble accumulation mode black carbon m.m.r. +!> @param[in] acc_sol_om Soluble accumulation mode organic m.m.r. +!> @param[in] acc_sol_ss Soluble accumulation mode sea salt m.m.r. +!> @param[in] n_cor_sol Soluble coarse mode number mixing ratio +!> @param[in] cor_sol_su Soluble coarse mode H2SO4 mass mixing ratio +!> @param[in] cor_sol_bc Soluble coarse mode black carbon m.m.r. +!> @param[in] cor_sol_om Soluble coarse mode organic m.m.r. +!> @param[in] cor_sol_ss Soluble coarse mode sea salt m.m.r. +!> @param[in] n_ait_ins Insoluble Aitken mode number mixing ratio +!> @param[in] ait_ins_bc Insoluble Aitken mode black carbon m.m.r. +!> @param[in] ait_ins_om Insoluble Aitken mode organic m.m.r. +!> @param[in] n_acc_ins Insoluble accumulation mode number m.r. +!> @param[in] acc_ins_du Insoluble accumulation mode dust m.m.r. +!> @param[in] n_cor_ins Insoluble coarse mode number mixing ratio +!> @param[in] cor_ins_du Insoluble coarse mode dust m.m.r. !> @param[in] ndf_wth Number of degrees of freedom per cell for !! potential temperature space !> @param[in] undf_wth Number unique of degrees of freedom for @@ -174,6 +219,18 @@ subroutine casim_code( nlayers, & cloud_drop_no_conc, murk, & refl_tot, refl_1km, & superc_liq, superc_rain, & + n_ait_sol, ait_sol_su, & + ait_sol_bc, ait_sol_om, & + n_acc_sol, acc_sol_su, & + acc_sol_bc, acc_sol_om, & + acc_sol_ss, & + n_cor_sol, cor_sol_su, & + cor_sol_bc, cor_sol_om, & + cor_sol_ss, & + n_ait_ins, ait_ins_bc, & + ait_ins_om, & + n_acc_ins, acc_ins_du, & + n_cor_ins, cor_ins_du, & ndf_wth, undf_wth, map_wth, & ndf_w3, undf_w3, map_w3, & ndf_2d, undf_2d, map_2d ) @@ -189,12 +246,19 @@ subroutine casim_code( nlayers, & use atm_fields_bounds_mod, only: pdims + use nlsizes_namelist_mod, only: tr_ukca + use mphys_inputs_mod, only: casim_aerosol_option, & + casim_iopt_act, fixed_number + use aerosol_extract_convert_mod, only: aerosol_extract_convert + use casim_ukca_tracer_mod, only: casim_ukca_tracer_column + use planet_constants_mod, only: p_zero, kappa, planet_radius use water_constants_mod, only: tm use micro_main, only: shipway_microphysics use casim_switches, only: its, ite, jts, jte, kts, kte, & - ils, ile, jls, jle + ils, ile, jls, jle, & + l_ukca_aerosol use generic_diagnostic_variables, & only: allocate_diagnostic_space, & deallocate_diagnostic_space, & @@ -228,6 +292,28 @@ subroutine casim_code( nlayers, & real(kind=r_def), intent(in), dimension(undf_w3) :: dry_rho_in_w3 real(kind=r_def), intent(in), dimension(undf_w3) :: height_w3 + real(kind=r_def), intent(in), dimension(undf_wth) :: n_ait_sol + real(kind=r_def), intent(in), dimension(undf_wth) :: ait_sol_su + real(kind=r_def), intent(in), dimension(undf_wth) :: ait_sol_bc + real(kind=r_def), intent(in), dimension(undf_wth) :: ait_sol_om + real(kind=r_def), intent(in), dimension(undf_wth) :: n_acc_sol + real(kind=r_def), intent(in), dimension(undf_wth) :: acc_sol_su + real(kind=r_def), intent(in), dimension(undf_wth) :: acc_sol_bc + real(kind=r_def), intent(in), dimension(undf_wth) :: acc_sol_om + real(kind=r_def), intent(in), dimension(undf_wth) :: acc_sol_ss + real(kind=r_def), intent(in), dimension(undf_wth) :: n_cor_sol + real(kind=r_def), intent(in), dimension(undf_wth) :: cor_sol_su + real(kind=r_def), intent(in), dimension(undf_wth) :: cor_sol_bc + real(kind=r_def), intent(in), dimension(undf_wth) :: cor_sol_om + real(kind=r_def), intent(in), dimension(undf_wth) :: cor_sol_ss + real(kind=r_def), intent(in), dimension(undf_wth) :: n_ait_ins + real(kind=r_def), intent(in), dimension(undf_wth) :: ait_ins_bc + real(kind=r_def), intent(in), dimension(undf_wth) :: ait_ins_om + real(kind=r_def), intent(in), dimension(undf_wth) :: n_acc_ins + real(kind=r_def), intent(in), dimension(undf_wth) :: acc_ins_du + real(kind=r_def), intent(in), dimension(undf_wth) :: n_cor_ins + real(kind=r_def), intent(in), dimension(undf_wth) :: cor_ins_du + real(kind=r_def), intent(inout), dimension(undf_wth) :: nl_mphys real(kind=r_def), intent(inout), dimension(undf_wth) :: nr_mphys real(kind=r_def), intent(inout), dimension(undf_wth) :: ni_mphys @@ -309,6 +395,14 @@ subroutine casim_code( nlayers, & deltaz, rhodz_dry, rhodz_moist, rho_r2, dry_rho, r_rho_levels real(r_um), dimension(1,1,0:nlayers) :: r_theta_levels + ! Pressure and temperature on theta levels, shaped as the UM aerosol + ! extraction routine expects them + real(r_um), dimension(1,1,0:nlayers) :: p_layer_centres + real(r_um), dimension(1,1,nlayers) :: t_layer_centres + + ! UKCA tracer array holding the GLOMAP modes for this column + real(r_um), allocatable :: tracer_ukca(:,:,:,:) + integer(i_um) :: k logical :: supercooled_layer(nlayers) @@ -320,8 +414,11 @@ subroutine casim_code( nlayers, & ! Configure optional diagnostics casdiags % l_graupfall_3d = ls_graup_3d_flag - ! Set CDNC for radiation here as we need the start of timestep value - if (casim_cdnc_opt == casim_cdnc_opt_fixed) then + ! Set CDNC for radiation here as we need the start of timestep value. + ! The mechanistic activation options always work the cloud number out from + ! the aerosol, so the in-cloud number for radiation comes from CASIM too. + if (casim_iopt_act > fixed_number .or. & + casim_cdnc_opt == casim_cdnc_opt_fixed) then do k = 0, nlayers if (cfl_wth(map_wth(1) + k) > 0.001_r_def) then cloud_drop_no_conc(map_wth(1) + k) = max(nl_mphys(map_wth(1) + k) / & @@ -454,6 +551,49 @@ subroutine casim_code( nlayers, & dact_insol_number_casim(k,1,1) = 0.0_wp end do ! k + !----------------------------------------------------------------------- + ! Overwrite the zeroed aerosol above with the GLOMAP modal aerosol + !----------------------------------------------------------------------- + if ( l_ukca_aerosol .and. casim_aerosol_option > 0 ) then + + allocate( tracer_ukca(1,1,0:nlayers,tr_ukca) ) + + call casim_ukca_tracer_column( nlayers, undf_wth, map_wth(1), & + n_ait_sol, ait_sol_su, ait_sol_bc, & + ait_sol_om, & + n_acc_sol, acc_sol_su, acc_sol_bc, & + acc_sol_om, acc_sol_ss, & + n_cor_sol, cor_sol_su, cor_sol_bc, & + cor_sol_om, cor_sol_ss, & + n_ait_ins, ait_ins_bc, ait_ins_om, & + n_acc_ins, acc_ins_du, & + n_cor_ins, cor_ins_du, & + tracer_ukca ) + + do k = 0, nlayers + p_layer_centres(1,1,k) = p_zero * & + ( exner_in_wth(map_wth(1) + k) )**(1.0_r_um / kappa) + end do + + do k = 1, nlayers + t_layer_centres(1,1,k) = exner_in_wth(map_wth(1) + k) * & + theta_in_wth(map_wth(1) + k) + end do + + call aerosol_extract_convert( p_layer_centres, t_layer_centres, & + rho_casim, tracer_ukca, & + aitken_sol_mass, aitken_sol_number, & + accum_sol_mass, accum_sol_number, & + coarse_sol_mass, coarse_sol_number, & + accum_dust_mass, accum_dust_number, & + coarse_dust_mass, coarse_dust_number, & + aitken_sol_bk, accum_sol_bk, & + coarse_sol_bk ) + + deallocate( tracer_ukca ) + + end if + cfrain_casim(nlayers,:,:)=0.0_wp cfgr_casim(nlayers,:,:)=0.0_wp do k = nlayers-1, 1, -1 diff --git a/interfaces/physics_schemes_interface/source/legacy/grids/nlsizes_namelist_mod.F90 b/interfaces/physics_schemes_interface/source/legacy/grids/nlsizes_namelist_mod.F90 index ead10a514c..da88f29fc0 100644 --- a/interfaces/physics_schemes_interface/source/legacy/grids/nlsizes_namelist_mod.F90 +++ b/interfaces/physics_schemes_interface/source/legacy/grids/nlsizes_namelist_mod.F90 @@ -3,6 +3,10 @@ ! For further details please refer to the file COPYRIGHT.txt ! which you should have received as part of this distribution. ! *****************************COPYRIGHT******************************* +!------------------------------------------------------------------------------- +! Some of the content of this file has been produced with the assistance of +! Anthropic Claude Opus 5 (Claude Code). +!------------------------------------------------------------------------------- ! ! Sizes for many of the UM's main, dynamic data arrays ! @@ -46,4 +50,6 @@ module nlsizes_namelist_mod integer :: n_cca_lev ! Number of CCA levels +integer :: tr_ukca = 0 ! Number of UKCA tracers + end module nlsizes_namelist_mod diff --git a/interfaces/physics_schemes_interface/source/support/casim_ukca_tracer_mod.F90 b/interfaces/physics_schemes_interface/source/support/casim_ukca_tracer_mod.F90 new file mode 100644 index 0000000000..00d0685a44 --- /dev/null +++ b/interfaces/physics_schemes_interface/source/support/casim_ukca_tracer_mod.F90 @@ -0,0 +1,230 @@ +!------------------------------------------------------------------------------- +! (c) Crown copyright 2026 Met Office. All rights reserved. +! The file LICENCE, distributed with this code, contains details of the terms +! under which the code may be used. +!------------------------------------------------------------------------------- +! Some of the content of this file has been produced with the assistance of +! Anthropic Claude Opus 5 (Claude Code). +!------------------------------------------------------------------------------- +!> @brief Gather the GLOMAP aerosol fields into a UKCA tracer array column +!> @details The UM CASIM aerosol routines locate the GLOMAP mode numbers and +!! component mass mixing ratios within the UKCA tracer array, using +!! the indexing set up by ukca_set_conv_indices. This routine builds +!! that array for a single LFRic column from the individual GLOMAP +!! prognostic fields. Only the GLOMAP entries are filled, as the +!! chemistry tracers are never looked at by CASIM. + +module casim_ukca_tracer_mod + + use constants_mod, only: i_def, r_def, r_um + + implicit none + + private + public :: casim_ukca_tracer_column + +contains + + !> @brief Fill a column of the UKCA tracer array from the GLOMAP fields + !> @param[in] nlayers Number of layers + !> @param[in] undf_wth Number of unique degrees of freedom for potential + !! temperature space + !> @param[in] base_wth Dofmap entry for the bottom of the column in + !! potential temperature space + !> @param[in] n_ait_sol Soluble Aitken mode number mixing ratio + !> @param[in] ait_sol_su Soluble Aitken mode H2SO4 mass mixing ratio + !> @param[in] ait_sol_bc Soluble Aitken mode black carbon mass mixing ratio + !> @param[in] ait_sol_om Soluble Aitken mode organic mass mixing ratio + !> @param[in] n_acc_sol Soluble accumulation mode number mixing ratio + !> @param[in] acc_sol_su Soluble accumulation mode H2SO4 mass mixing ratio + !> @param[in] acc_sol_bc Soluble accumulation mode black carbon m.m.r. + !> @param[in] acc_sol_om Soluble accumulation mode organic m.m.r. + !> @param[in] acc_sol_ss Soluble accumulation mode sea salt m.m.r. + !> @param[in] n_cor_sol Soluble coarse mode number mixing ratio + !> @param[in] cor_sol_su Soluble coarse mode H2SO4 mass mixing ratio + !> @param[in] cor_sol_bc Soluble coarse mode black carbon m.m.r. + !> @param[in] cor_sol_om Soluble coarse mode organic m.m.r. + !> @param[in] cor_sol_ss Soluble coarse mode sea salt m.m.r. + !> @param[in] n_ait_ins Insoluble Aitken mode number mixing ratio + !> @param[in] ait_ins_bc Insoluble Aitken mode black carbon m.m.r. + !> @param[in] ait_ins_om Insoluble Aitken mode organic m.m.r. + !> @param[in] n_acc_ins Insoluble accumulation mode number mixing ratio + !> @param[in] acc_ins_du Insoluble accumulation mode dust m.m.r. + !> @param[in] n_cor_ins Insoluble coarse mode number mixing ratio + !> @param[in] cor_ins_du Insoluble coarse mode dust m.m.r. + !> @param[out] tracer_ukca_col UKCA tracer array for this column + subroutine casim_ukca_tracer_column( nlayers, undf_wth, base_wth, & + n_ait_sol, ait_sol_su, ait_sol_bc, & + ait_sol_om, & + n_acc_sol, acc_sol_su, acc_sol_bc, & + acc_sol_om, acc_sol_ss, & + n_cor_sol, cor_sol_su, cor_sol_bc, & + cor_sol_om, cor_sol_ss, & + n_ait_ins, ait_ins_bc, ait_ins_om, & + n_acc_ins, acc_ins_du, & + n_cor_ins, cor_ins_du, & + tracer_ukca_col ) + + use nlsizes_namelist_mod, only: tr_ukca + use um_ukca_init_mod, only: tracer_names, & + fldname_n_ait_sol, & + fldname_ait_sol_su, & + fldname_ait_sol_bc, & + fldname_ait_sol_om, & + fldname_n_acc_sol, & + fldname_acc_sol_su, & + fldname_acc_sol_bc, & + fldname_acc_sol_om, & + fldname_acc_sol_ss, & + fldname_n_cor_sol, & + fldname_cor_sol_su, & + fldname_cor_sol_bc, & + fldname_cor_sol_om, & + fldname_cor_sol_ss, & + fldname_n_ait_ins, & + fldname_ait_ins_bc, & + fldname_ait_ins_om, & + fldname_n_acc_ins, & + fldname_acc_ins_du, & + fldname_n_cor_ins, & + fldname_cor_ins_du + + implicit none + + integer(kind=i_def), intent(in) :: nlayers + integer(kind=i_def), intent(in) :: undf_wth + integer(kind=i_def), intent(in) :: base_wth + + real(kind=r_def), intent(in), dimension(undf_wth) :: n_ait_sol + real(kind=r_def), intent(in), dimension(undf_wth) :: ait_sol_su + real(kind=r_def), intent(in), dimension(undf_wth) :: ait_sol_bc + real(kind=r_def), intent(in), dimension(undf_wth) :: ait_sol_om + real(kind=r_def), intent(in), dimension(undf_wth) :: n_acc_sol + real(kind=r_def), intent(in), dimension(undf_wth) :: acc_sol_su + real(kind=r_def), intent(in), dimension(undf_wth) :: acc_sol_bc + real(kind=r_def), intent(in), dimension(undf_wth) :: acc_sol_om + real(kind=r_def), intent(in), dimension(undf_wth) :: acc_sol_ss + real(kind=r_def), intent(in), dimension(undf_wth) :: n_cor_sol + real(kind=r_def), intent(in), dimension(undf_wth) :: cor_sol_su + real(kind=r_def), intent(in), dimension(undf_wth) :: cor_sol_bc + real(kind=r_def), intent(in), dimension(undf_wth) :: cor_sol_om + real(kind=r_def), intent(in), dimension(undf_wth) :: cor_sol_ss + real(kind=r_def), intent(in), dimension(undf_wth) :: n_ait_ins + real(kind=r_def), intent(in), dimension(undf_wth) :: ait_ins_bc + real(kind=r_def), intent(in), dimension(undf_wth) :: ait_ins_om + real(kind=r_def), intent(in), dimension(undf_wth) :: n_acc_ins + real(kind=r_def), intent(in), dimension(undf_wth) :: acc_ins_du + real(kind=r_def), intent(in), dimension(undf_wth) :: n_cor_ins + real(kind=r_def), intent(in), dimension(undf_wth) :: cor_ins_du + + ! Shaped as the UM tracer array is, for a single column, so that it can be + ! handed straight to the UM CASIM aerosol routines + real(kind=r_um), intent(out) :: tracer_ukca_col(1, 1, 0:nlayers, tr_ukca) + + integer(kind=i_def) :: k, n + + ! The chemistry tracers are not used by CASIM, and neither are the + ! soluble accumulation and coarse mode dust masses, for which LFRic has + ! no prognostic. Everything therefore starts as zero and only the fields + ! which LFRic holds are copied in. + tracer_ukca_col(:,:,:,:) = 0.0_r_um + + do n = 1, tr_ukca + select case (tracer_names(n)) + + case (fldname_n_ait_sol) + do k = 0, nlayers + tracer_ukca_col(1,1,k,n) = real( n_ait_sol(base_wth + k), r_um ) + end do + case (fldname_ait_sol_su) + do k = 0, nlayers + tracer_ukca_col(1,1,k,n) = real( ait_sol_su(base_wth + k), r_um ) + end do + case (fldname_ait_sol_bc) + do k = 0, nlayers + tracer_ukca_col(1,1,k,n) = real( ait_sol_bc(base_wth + k), r_um ) + end do + case (fldname_ait_sol_om) + do k = 0, nlayers + tracer_ukca_col(1,1,k,n) = real( ait_sol_om(base_wth + k), r_um ) + end do + + case (fldname_n_acc_sol) + do k = 0, nlayers + tracer_ukca_col(1,1,k,n) = real( n_acc_sol(base_wth + k), r_um ) + end do + case (fldname_acc_sol_su) + do k = 0, nlayers + tracer_ukca_col(1,1,k,n) = real( acc_sol_su(base_wth + k), r_um ) + end do + case (fldname_acc_sol_bc) + do k = 0, nlayers + tracer_ukca_col(1,1,k,n) = real( acc_sol_bc(base_wth + k), r_um ) + end do + case (fldname_acc_sol_om) + do k = 0, nlayers + tracer_ukca_col(1,1,k,n) = real( acc_sol_om(base_wth + k), r_um ) + end do + case (fldname_acc_sol_ss) + do k = 0, nlayers + tracer_ukca_col(1,1,k,n) = real( acc_sol_ss(base_wth + k), r_um ) + end do + + case (fldname_n_cor_sol) + do k = 0, nlayers + tracer_ukca_col(1,1,k,n) = real( n_cor_sol(base_wth + k), r_um ) + end do + case (fldname_cor_sol_su) + do k = 0, nlayers + tracer_ukca_col(1,1,k,n) = real( cor_sol_su(base_wth + k), r_um ) + end do + case (fldname_cor_sol_bc) + do k = 0, nlayers + tracer_ukca_col(1,1,k,n) = real( cor_sol_bc(base_wth + k), r_um ) + end do + case (fldname_cor_sol_om) + do k = 0, nlayers + tracer_ukca_col(1,1,k,n) = real( cor_sol_om(base_wth + k), r_um ) + end do + case (fldname_cor_sol_ss) + do k = 0, nlayers + tracer_ukca_col(1,1,k,n) = real( cor_sol_ss(base_wth + k), r_um ) + end do + + case (fldname_n_ait_ins) + do k = 0, nlayers + tracer_ukca_col(1,1,k,n) = real( n_ait_ins(base_wth + k), r_um ) + end do + case (fldname_ait_ins_bc) + do k = 0, nlayers + tracer_ukca_col(1,1,k,n) = real( ait_ins_bc(base_wth + k), r_um ) + end do + case (fldname_ait_ins_om) + do k = 0, nlayers + tracer_ukca_col(1,1,k,n) = real( ait_ins_om(base_wth + k), r_um ) + end do + + case (fldname_n_acc_ins) + do k = 0, nlayers + tracer_ukca_col(1,1,k,n) = real( n_acc_ins(base_wth + k), r_um ) + end do + case (fldname_acc_ins_du) + do k = 0, nlayers + tracer_ukca_col(1,1,k,n) = real( acc_ins_du(base_wth + k), r_um ) + end do + + case (fldname_n_cor_ins) + do k = 0, nlayers + tracer_ukca_col(1,1,k,n) = real( n_cor_ins(base_wth + k), r_um ) + end do + case (fldname_cor_ins_du) + do k = 0, nlayers + tracer_ukca_col(1,1,k,n) = real( cor_ins_du(base_wth + k), r_um ) + end do + + end select + end do + + end subroutine casim_ukca_tracer_column + +end module casim_ukca_tracer_mod diff --git a/interfaces/physics_schemes_interface/source/support/um_physics_init_mod.f90 b/interfaces/physics_schemes_interface/source/support/um_physics_init_mod.f90 index d1ae29c5f0..a9613dbfe9 100644 --- a/interfaces/physics_schemes_interface/source/support/um_physics_init_mod.f90 +++ b/interfaces/physics_schemes_interface/source/support/um_physics_init_mod.f90 @@ -3,6 +3,9 @@ ! The file LICENCE, distributed with this code, contains details of the terms ! under which the code may be used. !---------------------------------------------------------------------------- +! Some of the content of this file has been produced with the assistance of +! Anthropic Claude Opus 5 (Claude Code). +!---------------------------------------------------------------------------- !> @brief Controls the setting of variables for UM physics schemes, which !> are either fixed in LFRic or derived from LFRic inputs @@ -179,6 +182,9 @@ module um_physics_init_mod nsigmasf_in => nsigmasf, & nscalesf_in => nscalesf, & microphysics_casim, & + casim_activation, & + casim_activation_fixed, & + casim_activation_arg, & ci_input_in => ci_input, & cic_input_in => cic_input, & c_r_correl_in => c_r_correl, & @@ -456,6 +462,7 @@ subroutine um_physics_init() a_ratio_fac, l_droplet_tpr, qclrime, l_shape_rime, ndrop_surf, & z_surf, l_fsd_generator, mp_dz_scal, l_subgrid_qcl_mp, aut_qc, & l_mphys_nonshallow, casim_iopt_act, casim_iopt_inuc, & + abdul_razzak_ghan, & casim_aerosol_couple_choice,l_casim, & casim_aerosol_process_level,casim_moments_choice, & casim_aerosol_option, & @@ -478,7 +485,8 @@ subroutine um_physics_init() irs, ire, jrs, jre, krs, kre, & casim_moments_option, n_casim_tracers, & l_casim_warm_only, & - l_ukca_aerosol, no_aerosol_modes + l_ukca_aerosol, no_aerosol_modes, & + soluble_insoluble_modes use casim_stph, only: l_rp2_casim use casim_set_dependent_switches_mod, only: & casim_set_dependent_switches, & @@ -1374,15 +1382,46 @@ subroutine um_physics_init() casim_moments_option = 22222 ! all double moment - casim_iopt_act = 0_i_um ! 'fixed number' casim_aerosol_option = 0_i_um ! no soluble or insoluble aerosol modes casim_aerosol_process_level = 0_i_um casim_aerosol_couple_choice = 0_i_um l_ukca_aerosol = .false. + select case (casim_activation) + + case (casim_activation_fixed) + casim_iopt_act = fixed_number + + case (casim_activation_arg) + ! Mechanistic activation following Abdul-Razzak and Ghan, driven by + ! the GLOMAP modal aerosol, which must therefore be available. + if ( glomap_mode /= glomap_mode_ukca ) then + write( log_scratch_space, '(A)' ) & + 'CASIM Abdul-Razzak and Ghan activation requires the UKCA '// & + 'GLOMAP aerosol, stopping' + call log_event( log_scratch_space, LOG_LEVEL_ERROR ) + end if + casim_iopt_act = abdul_razzak_ghan + + case default + write( log_scratch_space, '(A,I0)' ) & + 'Invalid CASIM activation option, stopping', casim_activation + call log_event( log_scratch_space, LOG_LEVEL_ERROR ) + + end select + casim_moments_choice = 1_i_um CALL casim_set_dependent_switches + if (casim_iopt_act == abdul_razzak_ghan) then + ! casim_set_dependent_switches only asks for the insoluble modes when + ! the ice nucleation scheme needs them, but LFRic always supplies the + ! full set of soluble and insoluble GLOMAP modes to CASIM. This is + ! safe here because the aerosol option is not passed on to CASIM until + ! the call to set_mphys_switches below. + casim_aerosol_option = soluble_insoluble_modes + end if + ! Tell CASIM that its parent model is the UM. This allows for any UM-specific ! operations to take place within CASIM. casim_parent = parent_um diff --git a/interfaces/physics_schemes_interface/source/support/um_ukca_init_mod.f90 b/interfaces/physics_schemes_interface/source/support/um_ukca_init_mod.f90 index f233d3085b..d240428473 100644 --- a/interfaces/physics_schemes_interface/source/support/um_ukca_init_mod.f90 +++ b/interfaces/physics_schemes_interface/source/support/um_ukca_init_mod.f90 @@ -3,6 +3,10 @@ ! The file LICENCE, distributed with this code, contains details of the terms ! under which the code may be used. !---------------------------------------------------------------------------- +!---------------------------------------------------------------------------- +! Some of the content of this file has been produced with the assistance of +! Anthropic Claude Opus 5 (Claude Code). +!---------------------------------------------------------------------------- !> @brief UKCA initialisation subroutine for UM science configuration module um_ukca_init_mod @@ -63,6 +67,12 @@ module um_ukca_init_mod use ukca_option_mod, only: i_mode_nucscav, l_ukca_plume_scav use ukca_scavenging_mod, only: ukca_set_conv_indices, tracer_info + ! CASIM options, to determine whether the CASIM aerosol interface needs the + ! UKCA tracer indexing + use microphysics_config_mod, only: microphysics_casim, & + casim_activation, & + casim_activation_arg + ! Other UM modules used use dms_flux_mod_4a, only: i_liss_merlivat @@ -916,9 +926,10 @@ subroutine ukca_init( row_length, rows, model_levels, bl_levels, & l_ukca_mode = .true. end if - ! If the Easy Aerosol climatology is being used to set CDNC values + ! If the Easy Aerosol climatology or Casim is being used to set CDNC values ! then there is no need to calculate CDNCs via an activation scheme in UKCA - if (easyaerosol_cdnc) then + if (easyaerosol_cdnc.or. ( microphysics_casim .and. & + casim_activation == casim_activation_arg )) then i_tmp_ukca_activation_scheme = ukca_activation_off else i_tmp_ukca_activation_scheme = ukca_activation_arg @@ -1081,9 +1092,13 @@ subroutine ukca_init( row_length, rows, model_levels, bl_levels, & ! Retrieve the lists of required fields for the configuration call set_ukca_field_lists() - ! Set up indexing data needed for plume scavenging of UKCA tracers in the - ! GR convection scheme - if (l_ukca_plume_scav) then + ! Set up the indexing data that locates the GLOMAP mode numbers and + ! component masses within the UKCA tracer array. This is needed both for + ! plume scavenging of UKCA tracers in the GR convection scheme and by the + ! CASIM aerosol interface, which reads the same index arrays to build the + ! aerosol distribution seen by Abdul-Razzak and Ghan activation. + if ( l_ukca_plume_scav .or. ( microphysics_casim .and. & + casim_activation == casim_activation_arg ) ) then n = size(tracer_names) allocate(nm_spec_active(n)) do i = 1, n @@ -1513,6 +1528,8 @@ end subroutine ukca_emiss_init !>@brief Set the lists of fields required for the UKCA configuration subroutine set_ukca_field_lists() + use nlsizes_namelist_mod, only: tr_ukca + implicit none ! Local variables @@ -1536,6 +1553,9 @@ subroutine set_ukca_field_lists() ukca_errmsg call log_event( log_scratch_space, LOG_LEVEL_ERROR ) end if + ! Record the number of UKCA tracers where the UM code expects to find it. + tr_ukca = size(tracer_names) + write( log_scratch_space, '(A,I0,A)' ) & 'Tracers required (', size(tracer_names), '):' call log_event( log_scratch_space, LOG_LEVEL_INFO ) @@ -1897,7 +1917,9 @@ subroutine aerosol_ukca_dust_only_init( row_length, rows, model_levels, & call set_ukca_field_lists() ! Set up indexing data needed for plume scavenging of UKCA tracers in the - ! GR convection scheme + ! GR convection scheme. The CASIM aerosol interface also reads this + ! indexing, but only under Abdul-Razzak and Ghan activation, which requires + ! the full GLOMAP aerosol and so never reaches this dust-only path. if (l_ukca_plume_scav) then n = size(tracer_names) allocate(nm_spec_active(n)) diff --git a/rose-stem/app/lfric_atm/opt/rose-app-casim_ukca.conf b/rose-stem/app/lfric_atm/opt/rose-app-casim_ukca.conf new file mode 100644 index 0000000000..dae1a0e67e --- /dev/null +++ b/rose-stem/app/lfric_atm/opt/rose-app-casim_ukca.conf @@ -0,0 +1,27 @@ +[namelist:blayer] +noice_in_turb=.true. + +[namelist:microphysics] +!!aut_qc=2.47 +!!c_r_correl=0.9 +casim_activation='arg' +!!casim_cdnc_opt='external' +!!ci_input=14.3 +!!cic_input=1024.0 +!!fcrit=1.0 +!!graupel_scheme='modified' +!!heavy_rain_evap_fac=0.0 +!!l_mcr_precfrac=.false. +!!l_proc_fluxes=.false. +microphysics_casim=.true. +!!nscalesf=1.0 +!!nsigmasf=2.82843 +!!orog_block=.false. +!!orog_rain=.false. +!!orog_rime=.false. +!!prog_tnuc=.false. +!!shape_rime=.true. + +[namelist:radiation] +i_cloud_ice_type_lw=13 +i_cloud_ice_type_sw=13 diff --git a/rose-stem/site/common/lfric_atm/tasks_lfric_atm.cylc b/rose-stem/site/common/lfric_atm/tasks_lfric_atm.cylc index a37b78aceb..85ae638219 100644 --- a/rose-stem/site/common/lfric_atm/tasks_lfric_atm.cylc +++ b/rose-stem/site/common/lfric_atm/tasks_lfric_atm.cylc @@ -525,6 +525,20 @@ "plot_str": "plot_map.py -encorr $NODAL_DATA_DIR/lfric_diagnostics.nc $PLOT_DIR", }) %} +{% elif task_ns.conf_name == "clim_casim_ukca-C12" %} + + {% do task_dict.update({ + "opt_confs": ["climate","lowres_stp","casim_ukca","um_dump"], + "resolution": "C12", + "DT": 1800, + "tsteps": 144, + "crun": 2, + "crun_compare": false, + "mpi_parts": 6, + "kgo_checks": ["checksum"], + "plot_str": "plot_map.py -encorr $NODAL_DATA_DIR/lfric_diagnostics.nc $PLOT_DIR", + }) %} + {% elif task_ns.conf_name == "clim_gal9-C48_MG" %} {% do task_dict.update({ diff --git a/rose-stem/site/meto/groups/groups_lfric_atm.cylc b/rose-stem/site/meto/groups/groups_lfric_atm.cylc index 17717320f0..240213ac80 100644 --- a/rose-stem/site/meto/groups/groups_lfric_atm.cylc +++ b/rose-stem/site/meto/groups/groups_lfric_atm.cylc @@ -26,6 +26,7 @@ "lfric_atm_nwp_gal9_ls_and_jedi-C12_azspice_gnu_fast-debug-32bit", "lfric_atm_nwp_gal9-pert-C12_azspice_gnu_fast-debug-32bit", "lfric_atm_nwp_casim-C12_azspice_gnu_fast-debug-32bit", + "lfric_atm_clim_casim_ukca-C12_azspice_gnu_fast-debug-32bit", "lfric_atm_nwp_coma9-C12_azspice_gnu_fast-debug-32bit", "lfric_atm_nwp_comorph_dev-C12_azspice_gnu_fast-debug-32bit", "lfric_atm_nwp_gal9_mol-C12_azspice_gnu_fast-debug-32bit", @@ -128,6 +129,7 @@ "lfric_atm_nwp_gal9_ls_and_jedi-C12_ex1a_cce_fast-debug-32bit", "lfric_atm_nwp_gal9-pert-C12_ex1a_cce_fast-debug-32bit", "lfric_atm_nwp_casim-C12_ex1a_cce_fast-debug-32bit", + "lfric_atm_clim_casim_ukca-C12_ex1a_cce_fast-debug-32bit", "lfric_atm_nwp_coma9-C12_ex1a_cce_fast-debug-32bit", "lfric_atm_nwp_comorph_tb-C12_ex1a_cce_fast-debug-32bit", "lfric_atm_nwp_comorph_dev-C12_ex1a_cce_fast-debug-32bit", diff --git a/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_clim_casim_ukca-C12_azspice_gnu_fast-debug-32bit.txt b/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_clim_casim_ukca-C12_azspice_gnu_fast-debug-32bit.txt new file mode 100644 index 0000000000..2e857224bd --- /dev/null +++ b/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_clim_casim_ukca-C12_azspice_gnu_fast-debug-32bit.txt @@ -0,0 +1,9 @@ +Inner product checksum rho = 46D5584E +Inner product checksum theta = 518ECF8B +Inner product checksum u = 6AF7E2A6 +Inner product checksum mr1 = 3FDDA206 +Inner product checksum mr2 = 3715935C +Inner product checksum mr3 = 360A34F2 +Inner product checksum mr4 = 3672250C +Inner product checksum mr5 = 310CE5EA +Inner product checksum mr6 = 3515D610 diff --git a/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_clim_casim_ukca-C12_ex1a_cce_fast-debug-32bit.txt b/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_clim_casim_ukca-C12_ex1a_cce_fast-debug-32bit.txt new file mode 100644 index 0000000000..8bc815f654 --- /dev/null +++ b/rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_clim_casim_ukca-C12_ex1a_cce_fast-debug-32bit.txt @@ -0,0 +1,9 @@ +Inner product checksum rho = 46D556F2 +Inner product checksum theta = 518EC977 +Inner product checksum u = 6AF88A4F +Inner product checksum mr1 = 3FDDD8DC +Inner product checksum mr2 = 3717A5E0 +Inner product checksum mr3 = 3623871F +Inner product checksum mr4 = 368D771E +Inner product checksum mr5 = 30C75F82 +Inner product checksum mr6 = 352CAAAA diff --git a/rose-stem/site/meto/lfric_atm/tasks_lfric_atm_ex1a.cylc b/rose-stem/site/meto/lfric_atm/tasks_lfric_atm_ex1a.cylc index 5238590dfb..16d4950828 100644 --- a/rose-stem/site/meto/lfric_atm/tasks_lfric_atm_ex1a.cylc +++ b/rose-stem/site/meto/lfric_atm/tasks_lfric_atm_ex1a.cylc @@ -34,6 +34,12 @@ "memory": [24, "GB"], }) %} +{% elif task_ns.conf_name == "clim_casim_ukca-C12" %} + + {% do task_dict.update({ + "memory": [24, "GB"], + }) %} + {% elif task_ns.conf_name == "clim_gal9_short-C12" %} {% do task_dict.update({ diff --git a/science/gungho/source/algorithm/timestepping/semi_implicit_timestep_alg_mod.X90 b/science/gungho/source/algorithm/timestepping/semi_implicit_timestep_alg_mod.X90 index 298decc8ba..fa57ebc5cf 100644 --- a/science/gungho/source/algorithm/timestepping/semi_implicit_timestep_alg_mod.X90 +++ b/science/gungho/source/algorithm/timestepping/semi_implicit_timestep_alg_mod.X90 @@ -3,6 +3,9 @@ ! The file LICENCE, distributed with this code, contains details of the terms ! under which the code may be used. !----------------------------------------------------------------------------- +! Some of the content of this file has been produced with the assistance of +! Anthropic Claude Opus 5 (Claude Code). +!----------------------------------------------------------------------------- !> @brief A two time-level iterative time-discretisation of the nonlinear !! 3D equations. @@ -150,6 +153,7 @@ module semi_implicit_timestep_alg_mod use cld_alg_mod, only: cld_alg use aerosol_ukca_alg_mod, only: aerosol_ukca_alg use casim_activate_alg_mod, only: casim_activate_alg + use mphys_inputs_mod, only: casim_iopt_act, fixed_number use lsp_precfrac_checks_kernel_mod, only: lsp_precfrac_checks_kernel_type use stph_main_alg_mod, only : spt_array_names, & spt_array_count, & @@ -712,6 +716,8 @@ contains #ifdef UM_PHYSICS type(field_type), pointer :: precfrac + type(field_type), pointer :: liquid_fraction + type(field_type), pointer :: cf_liq_pre_fast #endif type(operator_type), pointer :: mm_wt @@ -738,7 +744,7 @@ contains character(str_def) :: prime_mesh_name integer(i_def) :: lbc_option logical(l_def) :: l_mcr_precfrac - logical(l_def) :: microphysics_casim + logical(l_def) :: microphysics_casim = .false. logical(l_def) :: murk_lbc real(r_def) :: tau_r @@ -990,6 +996,17 @@ contains call cld_incs_output(cloud_fields, dcfl_adv, dcff_adv, dbcf_adv, & sec_adv, suffix_adv) end if +#ifdef UM_PHYSICS + ! The CASIM mechanistic droplet activation works out how much liquid + ! cloud the fast physics has created or removed, so keep a copy of the + ! liquid cloud as it is on entry to the fast physics. + if (microphysics_casim .and. casim_iopt_act > fixed_number .and. & + outer == outer_iterations) then + call cloud_fields%get_field('liquid_fraction', liquid_fraction) + call microphysics_fields%get_field('cf_liq_pre_fast', cf_liq_pre_fast) + call invoke( setval_X(cf_liq_pre_fast, liquid_fraction) ) + end if +#endif !-------------------------------------------------------------------- ! Call the fast physics terms !-------------------------------------------------------------------- @@ -1229,8 +1246,10 @@ contains if (microphysics_casim) then call casim_activate_alg( self%config, & self%state(igh_t), mr, & + self%mr_after_adv(imr_cl), & derived_fields, cloud_fields, & microphysics_fields, convection_fields, & + aerosol_fields, turbulence_fields, & initialise=.false. ) end if ! microphysics_casim diff --git a/science/gungho/source/driver/create_physics_prognostics_mod.F90 b/science/gungho/source/driver/create_physics_prognostics_mod.F90 index 6c2ffd64dc..235ae86b95 100644 --- a/science/gungho/source/driver/create_physics_prognostics_mod.F90 +++ b/science/gungho/source/driver/create_physics_prognostics_mod.F90 @@ -3,6 +3,9 @@ ! The file LICENCE, distributed with this code, contains details of the terms ! under which the code may be used. !------------------------------------------------------------------------------- +! Some of the content of this file has been produced with the assistance of +! Anthropic Claude Opus 5 (Claude Code). +!------------------------------------------------------------------------------- !> @brief create physics prognostics !> @details Creates the physics prognostic fields module create_physics_prognostics_mod @@ -134,7 +137,7 @@ module create_physics_prognostics_mod l_conv_prog_dq, & adv_conv_prog_dtheta, & adv_conv_prog_dq - use mphys_inputs_mod, only: casim_iopt_act, l_mcr_precfrac + use mphys_inputs_mod, only: casim_iopt_act, l_mcr_precfrac, fixed_number use bl_option_mod, only: l_calc_tau_at_p use cloud_inputs_mod, only: l_pc2_homog_conv_pressure use io_config_mod, only : checkpoint_read, checkpoint_write @@ -172,6 +175,7 @@ subroutine process_physics_prognostics(processor) logical(l_def) :: is_rad ! Flag for chemistry fields ! that are radiatively active logical(l_def) :: sst_pert_flag + logical(l_def) :: casim_mechanistic_activation character(len=str_def) :: mesh_name #endif @@ -550,7 +554,7 @@ subroutine process_physics_prognostics(processor) advection_flag = microphysics_casim call processor%apply(make_spec('nl_mphys', main%microphysics, & - adv_coll=if_adv((advection_flag .and. casim_iopt_act /= 0_i_def), & + adv_coll=if_adv((advection_flag .and. casim_iopt_act > fixed_number), & adv%last_adv), ckp=checkpoint_flag, empty = (.not. microphysics_casim))) call processor%apply(make_spec('nr_mphys', main%microphysics, & adv_coll=if_adv(advection_flag, adv%last_adv), ckp=checkpoint_flag, & @@ -565,6 +569,16 @@ subroutine process_physics_prognostics(processor) adv_coll=if_adv(advection_flag, adv%last_adv), ckp=checkpoint_flag, & empty = (.not. microphysics_casim) )) + ! Copy of the liquid cloud taken immediately before the fast physics is + ! called. The CASIM mechanistic activation uses this to work out how much + ! liquid cloud the fast physics has created or removed. It is rewritten + ! every timestep before it is read, so doesn't need checkpointing or + ! advecting. + casim_mechanistic_activation = ( microphysics_casim .and. & + casim_iopt_act > fixed_number ) + call processor%apply(make_spec('cf_liq_pre_fast', main%microphysics, & + Wtheta, empty = (.not. casim_mechanistic_activation) )) + ! 2D fields, don't need checkpointing call processor%apply(make_spec('ls_rain', main%microphysics, W3, twod=.true.)) call processor%apply(make_spec('ls_snow', main%microphysics, W3, twod=.true.)) diff --git a/science/gungho/source/driver/gungho_init_fields_mod.X90 b/science/gungho/source/driver/gungho_init_fields_mod.X90 index b00b7c2258..d85e1196b9 100644 --- a/science/gungho/source/driver/gungho_init_fields_mod.X90 +++ b/science/gungho/source/driver/gungho_init_fields_mod.X90 @@ -10,7 +10,7 @@ !> module gungho_init_fields_mod - use mr_indices_mod, only : nummr, mr_names + use mr_indices_mod, only : nummr, mr_names, imr_cl use moist_dyn_mod, only : num_moist_factors, moist_dyn_names use field_array_mod, only : field_array_type use field_mod, only : field_type @@ -956,10 +956,12 @@ subroutine create_model_data( modeldb, & ( prognostic_init_choice /= init_option_checkpoint_dump ) ) then call casim_activate_alg( modeldb%config, & theta, mr_array%bundle, & + mr_array%bundle(imr_cl), & derived_fields, & cloud_fields, & microphysics_fields, & convection_fields, & + aerosol_fields, turbulence_fields, & initialise=.true.) end if #endif diff --git a/science/physics_schemes/source/large_scale_precipitation/CASIM/aerosol_extract_convert_mod.F90 b/science/physics_schemes/source/large_scale_precipitation/CASIM/aerosol_extract_convert_mod.F90 new file mode 100644 index 0000000000..98f99382cd --- /dev/null +++ b/science/physics_schemes/source/large_scale_precipitation/CASIM/aerosol_extract_convert_mod.F90 @@ -0,0 +1,475 @@ +! *****************************COPYRIGHT******************************* +! (C) Crown copyright Met Office. All rights reserved. +! For further details please refer to the file COPYRIGHT.txt +! which you should have received as part of this distribution. +! *****************************COPYRIGHT******************************* +!------------------------------------------------------------------------------- +! Some of the content of this file has been produced with the assistance of +! Anthropic Claude Opus 5 (Claude Code). +!------------------------------------------------------------------------------- + +! Code Owner: Please refer to the UM file CodeOwners.txt +! This file belongs in section: large_scale_precipitation + +! This module contains the declarations and subroutines for the conversion +! and incorporation of aerosol from UKCA into the CASIM microphysics. +! +! +module aerosol_extract_convert_mod + +use um_types, only: real_umphys + +implicit none + +character(len=*), parameter, private :: ModuleName='AEROSOL_EXTRACT_CONVERT_MOD' + +contains + +subroutine aerosol_extract_convert( p_layer_centres, t_work, rhoCasim, & + tracer_ukca, & + AitkenSolMass, AitkenSolNumber, & + AccumSolMass, AccumSolNumber, & + CoarseSolMass, CoarseSolNumber, & + AccumDustMass, AccumDustNumber, & + CoarseDustMass, CoarseDustNumber, & + AitkenSolBk, AccumSolBk, CoarseSolBk ) + +! Subroutine to read aerosol fields from UKCA +! This routine returns the soluble aitken, accumulation, coarse, +! and insoluble accumulation, coarse mass and number. + +use atm_fields_bounds_mod, only: tdims, tdims_s +use nlsizes_namelist_mod, only: tr_ukca + +use ukca_config_specification_mod, only: glomap_variables + +use ukca_mode_setup, only: mode_ait_sol, mode_acc_sol, & + mode_cor_sol, mode_acc_insol, & + mode_cor_insol + +use ukca_scavenging_mod, only: nmr_index_um,mmr_index_um +use ukca_constants, only: mmw +use water_constants_mod, only: rho_water +use chemistry_constants_mod, only: boltzmann +use missing_data_mod, only: rmdi + +! Dr Hook Modules +use yomhook, only: lhook, dr_hook +use parkind1, only: jprb, jpim + +implicit none + +!------------------------------------------------------------------------------- +! Subroutine arguments +!------------------------------------------------------------------------------- + +! Pressure at layer centres [Pa] +real(kind=real_umphys), intent(in) :: & + p_layer_centres( tdims%i_start : tdims%i_end, & + tdims%j_start : tdims%j_end, & + 0 : tdims%k_end ) + +! Local Working temperature [K] +real(kind=real_umphys), intent(in) :: & + t_work( tdims%i_start : tdims%i_end, & + tdims%j_start : tdims%j_end, & + 1 : tdims%k_end ) + +!Air density for CASIM [kg m-3] +real(kind=real_umphys), intent(in) :: & + rhoCasim( 1 : tdims%k_end, & + tdims%i_start : tdims%i_end, & + tdims%j_start : tdims%j_end ) + +! UKCA tracers [] +real(kind=real_umphys), intent(in) :: & + tracer_ukca( tdims_s%i_start : tdims_s%i_end, & + tdims_s%j_start : tdims_s%j_end, & + tdims_s%k_start : tdims_s%k_end, tr_ukca ) + +! Accumulation Model Soluble Mass [kg kg-1] +real(kind=real_umphys), intent(out) :: & + AccumSolMass( 1 : tdims%k_end, & + tdims%i_start : tdims%i_end, & + tdims%j_start : tdims%j_end ) + +! Accumulation Model Soluble Number [kg-1] +real(kind=real_umphys), intent(out) :: & + AccumSolNumber( 1 : tdims%k_end, & + tdims%i_start : tdims%i_end, & + tdims%j_start : tdims%j_end ) + +! Aitken mode Soluble Mass [kg kg-1] +real(kind=real_umphys), intent(out) :: & + AitkenSolMass( 1 : tdims%k_end, & + tdims%i_start : tdims%i_end, & + tdims%j_start : tdims%j_end ) + +! Accumulation Model Soluble Number [kg-1] +real(kind=real_umphys), intent(out) :: & + AitkenSolNumber( 1 : tdims%k_end, & + tdims%i_start : tdims%i_end, & + tdims%j_start : tdims%j_end ) + +! Coarse Model Soluble Mass [kg kg-1] +real(kind=real_umphys), intent(out) :: & + CoarseSolMass( 1 : tdims%k_end, & + tdims%i_start : tdims%i_end, & + tdims%j_start : tdims%j_end ) + +! Coarse Model Soluble Number [kg-1] +real(kind=real_umphys), intent(out) :: & + CoarseSolNumber( 1 : tdims%k_end, & + tdims%i_start : tdims%i_end, & + tdims%j_start : tdims%j_end ) + +! Coarse Dust Mass [kg kg-1] +real(kind=real_umphys), intent(out) :: & + CoarseDustMass( 1 : tdims%k_end, & + tdims%i_start : tdims%i_end, & + tdims%j_start : tdims%j_end ) + +! Coarse Dust Number [kg-1] +real(kind=real_umphys), intent(out) :: & + CoarseDustNumber( 1 : tdims%k_end, & + tdims%i_start : tdims%i_end, & + tdims%j_start : tdims%j_end ) + +! Accumulation Mode Dust Mass [kg kg-1] +real(kind=real_umphys), intent(out) :: & + AccumDustMass( 1 : tdims%k_end, & + tdims%i_start : tdims%i_end, & + tdims%j_start : tdims%j_end ) + +! Accumulation Mode Dust Number [kg-1] +real(kind=real_umphys), intent(out) :: & + AccumDustNumber( 1 : tdims%k_end, & + tdims%i_start : tdims%i_end, & + tdims%j_start : tdims%j_end ) + +! Abdul-Razzak-Ghan parameters (volume weighted) +real(kind=real_umphys), intent(out) :: & + AccumSolBk( 1 : tdims%k_end, & + tdims%i_start : tdims%i_end, & + tdims%j_start : tdims%j_end) +real(kind=real_umphys), intent(out) :: & + AitkenSolBk( 1 : tdims%k_end, & + tdims%i_start : tdims%i_end, & + tdims%j_start : tdims%j_end) +real(kind=real_umphys), intent(out) :: & + CoarseSolBk( 1 : tdims%k_end, & + tdims%i_start : tdims%i_end, & + tdims%j_start : tdims%j_end) +!------------------------------------------------------------------------------- +! Local variables +!------------------------------------------------------------------------------- + +! Caution - pointers to type glomap_variables% +! have been included here to make the code easier to read +! take care when making changes involving pointers +logical, pointer :: component(:,:) +real, pointer :: mm (:) +logical, pointer :: mode (:) +integer, pointer :: ncp +real, pointer :: no_ions (:) +real, pointer :: rhocomp (:) + +! Local working density [kg m-3] + +real(kind=real_umphys) :: & + AccumSolVolume( 1 : tdims%k_end, & + tdims%i_start : tdims%i_end, & + tdims%j_start : tdims%j_end) +real(kind=real_umphys) :: & + AitkenSolVolume( 1 : tdims%k_end, & + tdims%i_start : tdims%i_end, & + tdims%j_start : tdims%j_end) +real(kind=real_umphys) :: & + CoarseSolVolume( 1 : tdims%k_end, & + tdims%i_start : tdims%i_end, & + tdims%j_start : tdims%j_end) +real(kind=real_umphys) :: aird +integer :: i, j, k ! Loop counters + +integer :: imode,icp,i_cpt ! UKCA counters + +real(kind=real_umphys), parameter :: real_eps = epsilon(1.0_real_umphys) + +character(len=*), parameter :: RoutineName='AEROSOL_EXTRACT_CONVERT' + +! Declarations for Dr Hook +integer(kind=jpim), parameter :: zhook_in = 0 +integer(kind=jpim), parameter :: zhook_out = 1 +real(kind=jprb) :: zhook_handle + +!------------------------------------------------------------------------------- +! End of declarations and start of subroutine +!------------------------------------------------------------------------------- +if (lhook) call dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle) + +! Caution - pointers to type glomap_variables% +! have been included here to make the code easier to read +! take care when making changes involving pointers +component => glomap_variables%component +mm => glomap_variables%mm +mode => glomap_variables%mode +ncp => glomap_variables%ncp +no_ions => glomap_variables%no_ions +rhocomp => glomap_variables%rhocomp + +! sum up component masses in each mode for total mass +! and copy number to inputs for CASIM. + +! initialise modal Mass arrays to zero (will contain sum of mass over all +! components) + +!$OMP PARALLEL DEFAULT(none) & +!$OMP SHARED( tdims, AitkenSolMass, AitkenSolNumber, CoarseSolMass, & +!$OMP CoarseSolNumber, AccumDustMass, AccumDustNumber, & +!$OMP CoarseDustMass, CoarseDustNumber, AccumSolMass, & +!$OMP AccumSolNumber, AitkenSolBk, AccumSolBk, CoarseSolBk, & +!$OMP AitkenSolVolume, AccumSolVolume, CoarseSolVolume, & +!$OMP tracer_ukca, p_layer_centres, no_ions, mm, rhocomp, & +!$OMP rhoCASIM, T_work, mode, component, ncp, mmr_index_um, & +!$OMP nmr_index_um) & +!$OMP private( i, j, k, aird, imode, icp, i_cpt ) +!$OMP do SCHEDULE(STATIC) +do j = tdims%j_start, tdims%j_end + do i = tdims%i_start, tdims%i_end + do k = 1, tdims%k_end + AitkenSolMass (k,i,j)=0.0 + AccumSolMass (k,i,j)=0.0 + CoarseSolMass (k,i,j)=0.0 + AccumDustMass (k,i,j)=0.0 + CoarseDustMass(k,i,j)=0.0 + AitkenSolBk (k,i,j)=0.0 + AccumSolBk (k,i,j)=0.0 + CoarseSolBk (k,i,j)=0.0 + AitkenSolVolume(k,i,j)=0.0 + AccumSolVolume(k,i,j)=0.0 + CoarseSolVolume(k,i,j)=0.0 + end do ! loop over k + end do ! loop over i +end do ! loop over j +!$OMP end do +! +! include modes 2 to 6 (Aitsol,accsol,corsol,Aitins,accins,corins) +do imode=mode_ait_sol,mode_cor_insol + if (mode(imode)) then + do icp=1,ncp + i_cpt=0 + if (component(imode,icp)) then + i_cpt=mmr_index_um(imode,icp) + !The code in the CASIM activation scheme that is replaced by this interface: + !start: + ! Bk=chem%vantHoff(i)*Mw*chem%density(i)/(chem%massMole(i)*rho_water) + ! no_ions,rho_comp from ukca_mode_setup,massMole is 0.132 + !which is ammonium sulphate (18*2+96). Density is set to 1777 in + ! mphys_constants.F90. + !end + ! The Bk value here does not assume ammonium sulphate but is a volume + ! weighted average of the Bk values that would come from the UKCA components + + ! For each active mode calculate mass to transfer to CASIM, and Bk numerator + ! and denominator + select case (imode) + + case (mode_ait_sol) +!$OMP do SCHEDULE(STATIC) + do j = tdims%j_start, tdims%j_end + do i = tdims%i_start, tdims%i_end + do k = 1, tdims%k_end + AitkenSolMass(k,i,j) = AitkenSolMass(k,i,j) + & + tracer_ukca(i,j,k,i_cpt) + AitkenSolBk(k,i,j) = AitkenSolBk(k,i,j) +no_ions(icp) * & + mmw * tracer_ukca(i,j,k,i_cpt) / & + (rho_water*mm(icp)) + AitkenSolVolume(k,i,j) = AitkenSolVolume(k,i,j) + & + tracer_ukca(i,j,k,i_cpt) / rhocomp(icp) + end do ! loop over k + end do ! loop over i + end do ! loop over j +!$OMP end do + case (mode_acc_sol) +!$OMP do SCHEDULE(STATIC) + do j = tdims%j_start, tdims%j_end + do i = tdims%i_start, tdims%i_end + do k = 1, tdims%k_end + AccumSolMass(k,i,j) = AccumSolMass(k,i,j) + & + tracer_ukca(i,j,k,i_cpt) + AccumSolBk(k,i,j) = AccumSolBk(k,i,j) + no_ions(icp) * & + mmw * tracer_ukca(i,j,k,i_cpt) / & + (rho_water*mm(icp)) + AccumSolVolume(k,i,j) = AccumSolVolume(k,i,j) + & + tracer_ukca(i,j,k,i_cpt)/rhocomp(icp) + end do ! loop over k + end do ! loop over i + end do ! loop over j +!$OMP end do + case (mode_cor_sol) +!$OMP do SCHEDULE(STATIC) + do j = tdims%j_start, tdims%j_end + do i = tdims%i_start, tdims%i_end + do k = 1, tdims%k_end + CoarseSolMass(k,i,j) = CoarseSolMass(k,i,j) + & + tracer_ukca(i,j,k,i_cpt) + CoarseSolBk(k,i,j) = CoarseSolBk(k,i,j) + no_ions(icp) * & + mmw * tracer_ukca(i,j,k,i_cpt) / & + (rho_water*mm(icp)) + CoarseSolVolume(k,i,j) = CoarseSolVolume(k,i,j) + & + tracer_ukca(i,j,k,i_cpt) / rhocomp(icp) + end do ! loop over k + end do ! loop over i + end do ! loop over j +!$OMP end do + case (mode_acc_insol) +!$OMP do SCHEDULE(STATIC) + do j = tdims%j_start, tdims%j_end + do i = tdims%i_start, tdims%i_end + do k = 1, tdims%k_end + AccumDustMass (k,i,j) = AccumDustMass(k,i,j) + & + tracer_ukca(i,j,k,i_cpt) + end do ! loop over k + end do ! loop over i + end do ! loop over j +!$OMP end do + case (mode_cor_insol) +!$OMP do SCHEDULE(STATIC) + do j = tdims%j_start, tdims%j_end + do i = tdims%i_start, tdims%i_end + do k = 1, tdims%k_end + CoarseDustMass(k,i,j) = CoarseDustMass(k,i,j) + & + tracer_ukca(i,j,k,i_cpt) + end do ! loop over k + end do ! loop over i + end do ! loop over j +!$OMP end do + end select + end if ! if(component(imode,icp) + end do ! loop over icp + + ! For each active mode calculate number to transfer to CASIM + i_cpt = nmr_index_um(imode) + select case (imode) + + case (mode_ait_sol) +!$OMP do SCHEDULE(STATIC) + do k=1,tdims%k_end + do j = tdims%j_start, tdims%j_end + do i = tdims%i_start, tdims%i_end + aird = p_layer_centres(i,j,k)/ & + ( T_work(i,j,k)*boltzmann ) + AitkenSolNumber(k,i,j) = & + tracer_ukca(i,j,k,i_cpt)*aird/rhoCASIM(k,i,j) + end do ! loop over i + end do ! loop over j + end do ! loop over k +!$OMP end do + case (mode_acc_sol) +!$OMP do SCHEDULE(STATIC) + do k=1,tdims%k_end + do j = tdims%j_start, tdims%j_end + do i = tdims%i_start, tdims%i_end + aird = p_layer_centres(i,j,k)/ & + ( T_work(i,j,k)*boltzmann ) + AccumSolNumber(k,i,j) = & + tracer_ukca(i,j,k,i_cpt)*aird/rhoCASIM(k,i,j) + end do ! loop over i + end do ! loop over j + end do ! loop over k +!$OMP end do + case (mode_cor_sol) +!$OMP do SCHEDULE(STATIC) + do k=1,tdims%k_end + do j = tdims%j_start, tdims%j_end + do i = tdims%i_start, tdims%i_end + aird = p_layer_centres(i,j,k)/ & + ( T_work(i,j,k)*boltzmann ) + CoarseSolNumber(k,i,j) = & + tracer_ukca(i,j,k,i_cpt)*aird/rhoCASIM(k,i,j) + end do ! loop over i + end do ! loop over j + end do ! loop over k +!$OMP end do + case (mode_acc_insol) +!$OMP do SCHEDULE(STATIC) + do k=1,tdims%k_end + do j = tdims%j_start, tdims%j_end + do i = tdims%i_start, tdims%i_end + aird = p_layer_centres(i,j,k)/ & + ( T_work(i,j,k)*boltzmann ) + AccumDustNumber(k,i,j) = & + tracer_ukca(i,j,k,i_cpt)*aird/rhoCASIM(k,i,j) + end do ! loop over i + end do ! loop over j + end do ! loop over k +!$OMP end do + case (mode_cor_insol) +!$OMP do SCHEDULE(STATIC) + do k=1,tdims%k_end + do j = tdims%j_start, tdims%j_end + do i = tdims%i_start, tdims%i_end + aird = p_layer_centres(i,j,k)/ & + ( T_work(i,j,k)*boltzmann ) + CoarseDustNumber(k,i,j) = & + tracer_ukca(i,j,k,i_cpt)*aird/rhoCASIM(k,i,j) + end do ! loop over i + end do ! loop over j + end do ! loop over k +!$OMP end do + end select + end if ! if(mode(imode)) +end do ! loop over imode + +! calculate activation parameters Bk for soluble Aitken, accumulation, +! coarse modes. These modes are always used in UKCA, even if not always +! used in CASIM +! If there is no aerosol in a mode then CASIM will not carry out activation +! in that mode. In order to avoid carrying out spurious divides by zero here, +! we set the Bk to missing data if the volume is zero. +!$OMP do SCHEDULE(STATIC) +do j = tdims%j_start, tdims%j_end + do i = tdims%i_start, tdims%i_end + do k = 1, tdims%k_end + if (abs(AitkenSolVolume(k,i,j)) < real_eps) then + AitkenSolBk(k,i,j) = rmdi + else + AitkenSolBk(k,i,j) = AitkenSolBk(k,i,j)/AitkenSolVolume(k,i,j) + end if + end do + end do +end do +!$OMP end do NOWAIT +!$OMP do SCHEDULE(STATIC) +do j = tdims%j_start, tdims%j_end + do i = tdims%i_start, tdims%i_end + do k = 1, tdims%k_end + if (abs(AccumSolBk(k,i,j)) < real_eps) then + AccumSolBk(k,i,j) = rmdi + else + AccumSolBk(k,i,j) = AccumSolBk(k,i,j)/AccumSolVolume(k,i,j) + end if + end do + end do +end do +!$OMP end do NOWAIT +!$OMP do SCHEDULE(STATIC) +do j = tdims%j_start, tdims%j_end + do i = tdims%i_start, tdims%i_end + do k = 1, tdims%k_end + if (abs(CoarseSolBk(k,i,j)) < real_eps) then + CoarseSolBk(k,i,j) = rmdi + else + CoarseSolBk(k,i,j) = CoarseSolBk(k,i,j)/CoarseSolVolume(k,i,j) + end if + end do + end do +end do +!$OMP end do NOWAIT +!$OMP end PARALLEL + +if (lhook) call dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) + +end subroutine aerosol_extract_convert + +end module aerosol_extract_convert_mod diff --git a/science/physics_schemes/source/large_scale_precipitation/CASIM/casim_activation_in_um_mod.F90 b/science/physics_schemes/source/large_scale_precipitation/CASIM/casim_activation_in_um_mod.F90 new file mode 100644 index 0000000000..07408fa534 --- /dev/null +++ b/science/physics_schemes/source/large_scale_precipitation/CASIM/casim_activation_in_um_mod.F90 @@ -0,0 +1,473 @@ +! *****************************COPYRIGHT******************************* +! (C) Crown copyright Met Office. All rights reserved. +! For further details please refer to the file COPYRIGHT.txt +! which you should have received as part of this distribution. +! *****************************COPYRIGHT******************************* +!------------------------------------------------------------------------------- +! Some of the content of this file has been produced with the assistance of +! Anthropic Claude Opus 5 (Claude Code). +!------------------------------------------------------------------------------- + +! Code Owner: Please refer to the UM file CodeOwners.txt +! This file belongs in section: large_scale_precipitation + +module casim_activation_in_um_mod + +use micro_main, only: aerophys, aeroact, aerochem, dustact, dustphys, & + dustchem, aeroice, dustliq + + + +implicit none + +private + +public :: examine_ukca_aerosol_column, activate_column_ukca + +character(len=*), parameter, private :: & + ModuleName = 'CASIM_ACTIVATION_IN_UM_MOD' + +contains + +subroutine examine_ukca_aerosol_column(tracer_ukca, t_col, p_col, rho_col) + +!========================================================================== +! Description: +! This subroutine takes the ukca tracers (code variable tracer_ukca) and +! Sets up the CASIM aerosol species ready for an activation call outside +! of the main CASIM call. +!========================================================================== + +use um_types, only: real_umphys +use atm_fields_bounds_mod, only: tdims +use nlsizes_namelist_mod, only: tr_ukca +use casim_switches, only: l_casim_warm_only + +use mphys_switches, only: aero_index +use lognormal_funcs, only: MNtoRm +use mphys_parameters, only: nz +use thresholds, only: ccn_tidy, aeromass_small, aeronumber_small + +use ukca_config_specification_mod, only: glomap_variables + +use ukca_mode_setup, only: mode_ait_sol, mode_acc_sol, & + mode_cor_sol, mode_acc_insol, & + mode_cor_insol + +use ukca_scavenging_mod, only: nmr_index_um,mmr_index_um +use ukca_constants, only: mmw +use water_constants_mod, only: rho_water +use chemistry_constants_mod, only: boltzmann +use yomhook, only: lhook, dr_hook +use parkind1, only: jprb, jpim + +implicit none + +! Subroutine arguments (all intent(in) for now). +real(kind=real_umphys), intent(in) :: tracer_ukca(0:tdims%k_end, tr_ukca) +real(kind=real_umphys), intent(in) :: t_col(1:tdims%k_end) +real(kind=real_umphys), intent(in) :: p_col(1:tdims%k_end) +real(kind=real_umphys), intent(in) :: rho_col(1:tdims%k_end) + +! Local variables + +! Caution - pointers to type glomap_variables% +! have been included here to make the code easier to read +! take care when making changes involving pointers +logical, pointer :: component(:,:) +real, pointer :: mm (:) +logical, pointer :: mode (:) +integer, pointer :: ncp +real, pointer :: no_ions (:) +real, pointer :: rhocomp (:) + +! Aerosol density and air density +real(kind=real_umphys) :: density, aird(nz) + +integer :: k ! loop counter + +integer :: top_level ! top model level to work with + +integer :: imode, jmode, i_cpt, icp ! Aerosol modes and counters + +! Aerosol real variables (number, mass, numerator and volume). +real(kind=real_umphys) :: mode_N, mode_M, bk_numerator, volume, mode_Bk + +character(len=*), parameter :: RoutineName='EXAMINE_UKCA_AEROSOL_COLUMN' + +real(kind=real_umphys), parameter :: zero = 0.0_real_umphys +real(kind=real_umphys), parameter :: rcrit_high = 999.0_real_umphys + +integer(kind=jpim), parameter :: zhook_in = 0 +integer(kind=jpim), parameter :: zhook_out = 1 +real(kind=jprb) :: zhook_handle + +real(kind=real_umphys) :: eps_1 ! epsilon of 1.0 + +if (lhook) call dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle) + +! Caution - pointers to type glomap_variables% +! have been included here to make the code easier to read +! take care when making changes involving pointers +component => glomap_variables%component +mm => glomap_variables%mm +mode => glomap_variables%mode +ncp => glomap_variables%ncp +no_ions => glomap_variables%no_ions +rhocomp => glomap_variables%rhocomp + +eps_1 = epsilon(1.0_real_umphys) +top_level = min(tdims%k_end, nz) + +!-------------------------------------------------- +! Set up Soluble Aerosol Species +!-------------------------------------------------- + +do k = 1, top_level + density = aerochem(k) % density(1) + aird(k) = p_col(k) / (t_col(k)*boltzmann) + + do imode = 1, aero_index%nccn + + mode_N = zero + mode_M = zero + bk_numerator = zero + volume = zero + mode_Bk = zero + + if (imode == aero_index%i_aitken) then + jmode = mode_ait_sol + else if (imode == aero_index%i_accum) then + jmode = mode_acc_sol + else if (imode == aero_index%i_coarse) then + jmode = mode_cor_sol + end if + + ! Loop over aerosol modes to set up the aerosol physics + if (mode(jmode)) then + i_cpt = nmr_index_um(jmode) + mode_N = tracer_ukca(k, i_cpt) * aird(k) / rho_col(k) + + do icp = 1, ncp + if ( component(jmode, icp) ) then + + i_cpt = mmr_index_um(jmode,icp) + mode_M = mode_M + tracer_ukca(k, i_cpt) + + bk_numerator = bk_numerator + no_ions(icp) * mmw * & + tracer_ukca(k, i_cpt) / (rho_water*mm(icp)) + + volume = volume + tracer_ukca(k, i_cpt) / rhocomp(icp) + end if ! component(jmode, icp) + end do ! icp + + if (volume > zero) then + mode_Bk = bk_numerator / volume + else + mode_Bk = zero + end if + + end if ! mode(jmode) + + if (mode_N > ccn_tidy .and. mode_M > ccn_tidy * eps_1 .and. & + mode_Bk > zero) then + aerophys(k)%n(imode) = mode_N + aerophys(k)%m(imode) = mode_M + aerophys(k)%rd(imode) = MNtoRm( mode_M, mode_N, density, & + aerophys(k)%sigma(imode) ) + aerochem(k)%bk(imode) = mode_Bk + else + ! Zero the processes. + aerophys(k)%n(imode) = zero + aerophys(k)%m(imode) = zero + aerophys(k)%rd(imode) = zero + aerochem(k)%bk(imode) = zero + + end if ! mode_N > ccn_tidy etc + + + end do ! imode + + ! Set other aerosol properties for this level + aeroact(k) % nact = zero + aeroact(k) % mact = zero + aeroact(k) % rcrit = rcrit_high + aeroact(k) % mact_mean = zero + aeroact(k) % nact2 = zero + aeroact(k) % nact1 = zero + aeroact(k) % mact1 = zero + aeroact(k) % rcrit1 = rcrit_high + aeroact(k) % nact2 = zero + aeroact(k) % rcrit2 = rcrit_high + aeroact(k) % mact2 = zero + aeroact(k) % mact2_mean = zero + aeroact(k) % mact1_mean = zero + aeroact(k) % rcrit1 = rcrit_high + aeroact(k) % nact2 = zero + aeroact(k) % rcrit2 = rcrit_high + aeroact(k) % mact2 = zero + aeroact(k) % mact2_mean = zero + aeroact(k) % mact1_mean = zero + aeroact(k) % nratio1 = zero + aeroact(k) % nratio2 = zero +end do ! loop over levels (k) + +!-------------------------------------------------- +! Set up Insoluble Aerosol Species +!-------------------------------------------------- + +if (.not. l_casim_warm_only ) then + ! Activated dust + do k = 1, nz + density = dustchem(k) % density(1) + ! Examine interstitial dust + do imode = 1, aero_index % nin + + mode_N = zero + mode_M = zero + + if (imode == aero_index % i_coarse_dust) then + jmode = mode_cor_insol + end if + if (imode == aero_index % i_accum_dust) then + jmode = mode_acc_insol + end if + if (mode(jmode)) then + i_cpt = nmr_index_um(jmode) + mode_N = tracer_ukca(k, i_cpt) * aird(k) / rho_col(k) + + do icp = 1, ncp + if ( component(jmode, icp) ) then + i_cpt = mmr_index_um(jmode, icp) + mode_M = mode_M + tracer_ukca(k, i_cpt) + end if + end do + + end if + if (mode_m < zero .or. mode_n < zero ) then + if (mode_m < zero) mode_m = aeromass_small + if (mode_n < zero) mode_n = aeronumber_small + end if + + dustphys(k)%n(imode) = mode_N + dustphys(k)%m(imode) = mode_M + + dustphys(k)%rd(imode) = MNtoRm(mode_M, mode_N, density, & + dustphys(k)%sigma(imode) ) + end do + + ! Examine activated dust + ! Initialize to zero/defaults + dustact(k) % nact = zero + dustact(k) % mact = zero + dustact(k) % rcrit = rcrit_high + dustact(k) % mact_mean = zero + dustact(k) % nact1 = zero + dustact(k) % nratio1 = zero + dustact(k) % mact1 = zero + dustact(k) % rcrit1 = rcrit_high + dustact(k) % mact1_mean = zero + dustact(k) % mact2 = zero + dustact(k) % nact2 = zero + dustact(k) % nratio2 = zero + dustact(k) % rcrit2 = rcrit_high + dustact(k) % mact2_mean = zero + dustact(k) % mact3 = zero + dustact(k) % nact3 = zero + dustact(k) % nratio3 = zero + dustact(k) % rcrit3 = rcrit_high + dustact(k) % mact3_mean = zero + + end do +end if + +if (lhook) call dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) + +end subroutine examine_ukca_aerosol_column + +!============================================================================== + +subroutine activate_column_ukca( cloud_mass_post_qtbal, & + cloud_mass_pre_ap2, & + rho_col, t_col, p_col, & + cf_liquid, cf_liquid_pre_ap2, w_tke, & + cloudnumber_col, & + tracer_ukca_col ) + +!---------------------------------------------------------------------------- +! Code Description: +! Changes to the liquid and ice cloud mass can occur due to processes outside +! of the CASIM microphysics. These are often in the slow physics terms +! such as convection and boundary layer. For double-moment CASIM simulations, +! the cloud number and ice number have to be set in order to ensure that +! CASIM does not reject this new liquid and ice cloud. For UKCA simulations, +! this also means checking and activating the aerosol within the model column +! before calling the 'activate' routine from the CASIM repository to +! calculate the cloud number in the same way it would normally be done in +! the main CASIM call. +!---------------------------------------------------------------------------- + +use atm_fields_bounds_mod, only: tdims +use timestep_mod, only: timestep, recip_timestep +use um_types, only: real_umphys +use nlsizes_namelist_mod, only: tr_ukca +use yomhook, only: lhook, dr_hook +use parkind1, only: jprb, jpim +use mphys_parameters, only: nz +use mphys_switches, only: aero_index +use thresholds, only: ql_tidy, cfliq_small +use condensation, only: dnccn_all, dmac_all, dnccnd_all, dmad_all +use activation, only: activate + +implicit none + +! Subroutine arguments +real(kind=real_umphys), intent(in) :: cloud_mass_post_qtbal(1:tdims%k_end) + +real(kind=real_umphys), intent(in) :: cloud_mass_pre_ap2(1:tdims%k_end) + +real(kind=real_umphys), intent(in) :: rho_col(1:tdims%k_end) + +real(kind=real_umphys), intent(in) :: t_col(1:tdims%k_end) + +real(kind=real_umphys), intent(in) :: p_col(1:tdims%k_end) + +real(kind=real_umphys), intent(in) :: cf_liquid(1:tdims%k_end) + +real(kind=real_umphys), intent(in) :: cf_liquid_pre_ap2(1:tdims%k_end) + +real(kind=real_umphys), intent(in) :: w_tke(1:tdims%k_end) + +real(kind=real_umphys), intent(in) :: tracer_ukca_col(0:tdims%k_end, tr_ukca) + +real(kind=real_umphys), intent(in out) :: cloudnumber_col(1:tdims%k_end) + +! Local variables + +real(kind=real_umphys) :: w_cloud_number + +real(kind=real_umphys) :: delta_mass +real(kind=real_umphys) :: delta_cfliq +real(kind=real_umphys) :: dnumber, dnumber_d, dnumber_a +real(kind=real_umphys) :: dmac +real(kind=real_umphys) :: dmad + +real(kind=real_umphys) :: smax,ait_cdnc,accum_cdnc,tot_cdnc, & + activated_cloud, activated_arg + + + +real(kind=real_umphys), parameter :: zero = 0.0_real_umphys + +integer :: k ! loop counter +integer :: top_level ! highest level to use. + +integer(kind=jpim), parameter :: zhook_in = 0 +integer(kind=jpim), parameter :: zhook_out = 1 +real(kind=jprb) :: zhook_handle +character(len=*), parameter :: RoutineName='ACTIVATE_COLUMN_UKCA' + +if (lhook) call dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle) + +!----------------------------------- +! Initial set up of constants +!----------------------------------- +! Check for bounds - the top level cannot be greater than nz or +! tdims%k_end. +top_level = min(nz, tdims%k_end) + +!------------------------------------------------------------------------------ +! Call examine_ukca_aerosol_column. This uses the UKCA tracers to set up the +! aerosol amounts within CASIM for a single column. The internal aerosol +! amounts within CASIM are required for the activate routine to be able to +! set the cloud number based on the aerosol in each mode. +!------------------------------------------------------------------------------ +call examine_ukca_aerosol_column( tracer_ukca_col, t_col, p_col, rho_col ) + +!------------------------------------------------------------------------ +! Main loop over levels to call the activate routine from within CASIM. +!------------------------------------------------------------------------ +do k = 1, top_level + + ! Work out if cloud mass or liquid cloud fraction has changed. + delta_mass = cloud_mass_post_qtbal(k) - cloud_mass_pre_ap2(k) + delta_cfliq = cf_liquid(k) - cf_liquid_pre_ap2(k) + + ! Set local working cloud number for this level. + w_cloud_number = cloudnumber_col(k) + + if ( delta_mass > zero .and. cloud_mass_post_qtbal(k) > ql_tidy .and. & + cf_liquid(k) > cfliq_small ) then + + !----------------------------------------------------------------- + ! There has been an increase in liquid cloud due to non-CASIM + ! physical processes in the model. + !----------------------------------------------------------------- + ! Call CASIM routine activate (operates over one grid box) in + ! order to work out what the new cloud number is based on the + ! aerosol and the cloud mass added by the other processes. + !----------------------------------------------------------------- + + call activate(timestep, cloud_mass_post_qtbal(k), w_cloud_number, & + w_tke(k), rho_col(k), dnumber, dmac, t_col(k), p_col(k), & + cf_liquid(k), cf_liquid_pre_ap2(k), & + aerophys(k), aerochem(k), aeroact(k), & + dustphys(k), dustchem(k), dustliq(k), & + dnccn_all, dmac_all, dnumber_d, dmad, & + dnccnd_all, dmad_all, & + smax,ait_cdnc,accum_cdnc, tot_cdnc,activated_arg, & + activated_cloud) + + dnumber_a = dnumber + + else if ( cloud_mass_post_qtbal(k) < ql_tidy .or. & + cf_liquid(k) < cfliq_small ) then + + !----------------------------------------------------------------- + ! Non-CASIM processes have removed all of the liquid cloud or the + ! liquid cloud fraction. Therefore, we need to remove the cloud + ! number + !----------------------------------------------------------------- + + dnumber = -w_cloud_number * recip_timestep + + w_cloud_number = zero + + else + + if (delta_cfliq < zero) then + dnumber = w_cloud_number * delta_cfliq / cf_liquid_pre_ap2(k) * & + recip_timestep + dmac = zero + dnumber_a = dnumber + dmac_all( aero_index % i_accum) = dmac + ! put back into accum + dnccn_all( aero_index % i_accum) = dnumber_a + dnumber_d = zero + dnccnd_all = zero + dmad_all = zero + else + dnumber = zero + dnumber_a = zero ! No aerosol processing required + + dnccn_all = zero ! we assume no change in number during evap + dmac = zero ! No aerosol processing required + dmac_all = zero ! No aerosol processing required + dnumber_d = zero ! No aerosol processing required + dnccnd_all = zero + dmad_all = zero + end if ! delta_cfliq < zero + + end if ! delta_mass > zero etc. + + ! Update cloud number in the column + w_cloud_number = max(w_cloud_number + (dnumber * timestep), zero) + cloudnumber_col(k) = w_cloud_number + +end do + +if (lhook) call dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) + +end subroutine activate_column_ukca + +end module casim_activation_in_um_mod