diff --git a/applications/adjoint_tests/example/configuration.nml b/applications/adjoint_tests/example/configuration.nml index c358d14e69..eb76180860 100644 --- a/applications/adjoint_tests/example/configuration.nml +++ b/applications/adjoint_tests/example/configuration.nml @@ -34,10 +34,8 @@ stochastic_physics='none', surface='none', / &convection -dx_ref=50000.0, l_cvdiag_ctop_qmax=.false., qlmin=4.0e-4, -resdep_precipramp=.false., / &cosp l_cosp=.false., diff --git a/applications/jedi_lfric_tests/example/configuration.nml b/applications/jedi_lfric_tests/example/configuration.nml index 4e3c939fb0..717bbe2c01 100644 --- a/applications/jedi_lfric_tests/example/configuration.nml +++ b/applications/jedi_lfric_tests/example/configuration.nml @@ -33,10 +33,8 @@ stochastic_physics='none', surface='none', / &convection -dx_ref=50000.0, l_cvdiag_ctop_qmax=.false., qlmin=4.0e-4, -resdep_precipramp=.false., / &cosp l_cosp=.false., diff --git a/applications/jedi_lfric_tests/example_id_tlm_tests/configuration.nml b/applications/jedi_lfric_tests/example_id_tlm_tests/configuration.nml index 0744592c20..2f60bc925a 100644 --- a/applications/jedi_lfric_tests/example_id_tlm_tests/configuration.nml +++ b/applications/jedi_lfric_tests/example_id_tlm_tests/configuration.nml @@ -70,10 +70,8 @@ stochastic_physics='none', surface='none', / &convection -dx_ref=50000.0, l_cvdiag_ctop_qmax=.false., qlmin=4.0e-4, -resdep_precipramp=.false., / &cosp l_cosp=.false., diff --git a/applications/jedi_lfric_tests/example_tlm_forecast_tl/configuration.nml b/applications/jedi_lfric_tests/example_tlm_forecast_tl/configuration.nml index 4ddd6423f3..079def8985 100644 --- a/applications/jedi_lfric_tests/example_tlm_forecast_tl/configuration.nml +++ b/applications/jedi_lfric_tests/example_tlm_forecast_tl/configuration.nml @@ -69,10 +69,8 @@ stochastic_physics='none', surface='none', / &convection -dx_ref=50000.0, l_cvdiag_ctop_qmax=.false., qlmin=4.0e-4, -resdep_precipramp=.false., / &cosp l_cosp=.false., diff --git a/applications/jedi_lfric_tests/example_tlm_tests/configuration.nml b/applications/jedi_lfric_tests/example_tlm_tests/configuration.nml index ccef538775..e86dccd0fc 100644 --- a/applications/jedi_lfric_tests/example_tlm_tests/configuration.nml +++ b/applications/jedi_lfric_tests/example_tlm_tests/configuration.nml @@ -69,10 +69,8 @@ stochastic_physics='none', surface='none', / &convection -dx_ref=50000.0, l_cvdiag_ctop_qmax=.false., qlmin=4.0e-4, -resdep_precipramp=.false., / &cosp l_cosp=.false., diff --git a/applications/ngarch/example/configuration_bl.nml b/applications/ngarch/example/configuration_bl.nml index c132b13b31..8f71dc1efb 100644 --- a/applications/ngarch/example/configuration_bl.nml +++ b/applications/ngarch/example/configuration_bl.nml @@ -60,10 +60,8 @@ stochastic_physics='none', surface='jules', / &convection -dx_ref=50000.0, l_cvdiag_ctop_qmax=.false., qlmin=4.0e-4, -resdep_precipramp=.false., / &cosp l_cosp=.false., diff --git a/applications/ngarch/example/configuration_casim.nml b/applications/ngarch/example/configuration_casim.nml index e2ed420e9e..9c385834f9 100644 --- a/applications/ngarch/example/configuration_casim.nml +++ b/applications/ngarch/example/configuration_casim.nml @@ -100,14 +100,12 @@ chem_scheme='none', &convection cape_timescale=1800.0, cv_scheme='gregory_rowntree', -dx_ref=50000.0, efrac=1.0, l_cvdiag_ctop_qmax=.false., number_of_convection_substeps=2, orig_mdet_fac=1.0, prog_ent_min=0.5, qlmin=4.0e-4, -resdep_precipramp=.false., / &cosp l_cosp=.false., diff --git a/interfaces/physics_schemes_interface/rose-meta/um-convection/HEAD/rose-meta.conf b/interfaces/physics_schemes_interface/rose-meta/um-convection/HEAD/rose-meta.conf index 283c4d6b66..72db2adf37 100644 --- a/interfaces/physics_schemes_interface/rose-meta/um-convection/HEAD/rose-meta.conf +++ b/interfaces/physics_schemes_interface/rose-meta/um-convection/HEAD/rose-meta.conf @@ -2,6 +2,412 @@ # CONVECTION: #============================================================================== +[namelist:comorph] +compulsory=true +description=CoMorph convection scheme +ns=namelist/Science/UM Convection/CoMorph +sort-key=Section-A06b +title=CoMorph convection scheme + +[namelist:comorph=ass_min_radius] +compulsory=true +description=Asymptotic minimum allowed initial parcel radius / m +help=This prevents the turbulence-based parcel radius from becoming + =too small in elevated or nocturnal convection. +!kind=default +ns=namelist/Science/UM Convection/CoMorph/genesis +range=0.0:10000.0 +sort-key=gen-04a +type=real + +[namelist:comorph=autoc_opt] +compulsory=true +description=Switch for autoconversion method +!enumeration=true +help=Autoconversion in CoMorph can be linear or quadratic in parcel + =liquid water content. +!kind=default +ns=namelist/Science/UM Convection/CoMorph/micro +sort-key=micro-01a +value-titles=Linear autoconversion (rate = coef_auto * (q_cl - q_cl_auto)), + =Quadratic autoconversion (rate = coef_auto * q_cl^2) +values='linear','quadratic' + +[namelist:comorph=cf_area_coef] +compulsory=true +description=Ice crystal non-spherical area factor +help=This scales up the assumed surface area of ice crystals relative to + =what they would be if they were spheres. +!kind=default +ns=namelist/Science/UM Convection/CoMorph/micro +range=1.0:1000.0 +sort-key=micro-02b +type=real + +[namelist:comorph=cf_conv_fac] +compulsory=true +description=Scaling factor for convective cloud fraction +help=Controls the effective convective cloud fraction in Comorph. +!kind=default +ns=namelist/Science/UM Convection/CoMorph +sort-key=02a +type=real + +[namelist:comorph=coef_auto] +compulsory=true +description=In-parcel cloud-to-rain autoconversion rate coefficient / s-1 +help=Coefficient in the CoMorph autoconversion rate expression. +!kind=default +ns=namelist/Science/UM Convection/CoMorph/micro +sort-key=micro-01b +type=real + +[namelist:comorph=col_eff_coef] +compulsory=true +description=Coefficient in the deflection term affecting hydrometeor collisions +help=This reduces collection efficiency and therefore tends to reduce riming + =and accretion. +!kind=default +ns=namelist/Science/UM Convection/CoMorph/micro +range=0.0:10.0 +sort-key=micro-04a +type=real + +[namelist:comorph=core_ent_cmr] +compulsory=true +description=Include the factor of 1/core_mean_ratio in core dilution +help=Switches the core dilution correction on or off. +!kind=default +ns=namelist/Science/UM Convection/CoMorph/plume +sort-key=plume-01b +type=logical + +[namelist:comorph=core_ent_fac] +compulsory=true +description=Factor for parcel core entrainment (dimensionless) +help=core_ent = core_ent_fac * mean_ent. +!kind=default +ns=namelist/Science/UM Convection/CoMorph/plume +sort-key=plume-01c +type=real + +[namelist:comorph=drag_coef_cond] +compulsory=true +description=Asymptotic drag coefficient for a sphere at high Reynolds number +help=Mainly affects fall-speed for rain and graupel. +!kind=default +ns=namelist/Science/UM Convection/CoMorph/micro +range=0.2:1.0 +sort-key=micro-04b +type=real + +[namelist:comorph=drag_coef_par] +compulsory=true +description=Drag coefficient applied to parcel winds in CoMorph CMT +help=Scales the drag term used to relax parcel winds toward the environment. +!kind=default +ns=namelist/Science/UM Convection/CoMorph/plume +sort-key=plume-03a +type=real + +[namelist:comorph=dx_ref] +compulsory=true +description=Reference grid-length for resolution-dependence / m +help=When using the resolution-dependent link between parcel radius + =and precipitation rate, dx_ref is the grid-length at which the + =namelist input value of par_radius_ppn_max applies. + =Typically set to 50000.0m (ie. 50km). +!kind=default +ns=namelist/Science/UM Convection/CoMorph/genesis +range=100.0:1000000.0 +sort-key=gen-03f +type=real + +[namelist:comorph=ent_coef] +compulsory=true +description=Fractional entrainment coefficient for a given parcel radius +help=Scales the entrainment rate for the parcel radius formulation. +!kind=default +ns=namelist/Science/UM Convection/CoMorph/plume +range=0.1:0.4 +sort-key=plume-01a +type=real + +[namelist:comorph=hetnuc_temp] +compulsory=true +description=Heterogeneous nucleation temperature / K +help=Gradual freezing starts below this temperature. +!kind=default +ns=namelist/Science/UM Convection/CoMorph/micro +range=230.0:273.0 +sort-key=micro-02a +type=real + +[namelist:comorph=max_cmr] +compulsory=true +description=Maximum limit on the core-mean-ratio +help=Constrains the assumed PDF shape used in the detrainment calculation. +!kind=default +ns=namelist/Science/UM Convection/CoMorph/plume +range=3.0:10.0 +sort-key=plume-02b +type=real + +[namelist:comorph=min_cmr] +compulsory=true +description=Minimum limit on the core-mean-ratio +help=Constrains the assumed PDF shape used in the detrainment calculation. +!kind=default +ns=namelist/Science/UM Convection/CoMorph/plume +range=1.0:3.0 +sort-key=plume-02a +type=real + +[namelist:comorph=n_dndraft_types] +compulsory=true +description=Number of downdraught types +help=Controls the number of downdraughts in CoMorph. +!kind=default +ns=namelist/Science/UM Convection/CoMorph +range=0:1 +sort-key=01a +type=integer + +[namelist:comorph=nconc_cf] +compulsory=true +description=Number concentration for ice-cloud / kg-1 +help=Prescribed fixed number concentration for parcel ice cloud. +!kind=default +ns=namelist/Science/UM Convection/CoMorph/micro +sort-key=micro-05c +type=real + +[namelist:comorph=nconc_cl] +compulsory=true +description=Number concentration for liquid-cloud / kg-1 +help=Prescribed fixed number concentration for parcel liquid cloud. +!kind=default +ns=namelist/Science/UM Convection/CoMorph/micro +sort-key=micro-05a +type=real + +[namelist:comorph=nconc_graup] +compulsory=true +description=Number concentration for graupel / kg-1 +help=Prescribed fixed number concentration for parcel graupel. +!kind=default +ns=namelist/Science/UM Convection/CoMorph/micro +sort-key=micro-05e +type=real + +[namelist:comorph=nconc_rain] +compulsory=true +description=Number concentration for rain / kg-1 +help=Prescribed fixed number concentration for parcel rain. +!kind=default +ns=namelist/Science/UM Convection/CoMorph/micro +sort-key=micro-05b +type=real + +[namelist:comorph=nconc_snow] +compulsory=true +description=Number concentration for snow / kg-1 +help=Prescribed fixed number concentration for parcel snow. +!kind=default +ns=namelist/Science/UM Convection/CoMorph/micro +sort-key=micro-05d +type=real + +[namelist:comorph=overlap_power] +compulsory=true +description=Controls overlap between detrained liquid and ice cloud +help=Applies in mixed-phase parcels. 0 is total overlap, 1 is no overlap. +!kind=default +ns=namelist/Science/UM Convection/CoMorph +range=1.0e-6:1.0 +sort-key=03a +type=real + +[namelist:comorph=par_gen_core_fac] +compulsory=true +description=Ratio of parcel core initial perturbations over parcel mean +help=Controls how extreme the parcel core is compared to the parcel mean. +!kind=default +ns=namelist/Science/UM Convection/CoMorph/genesis +range=2.0:6.0 +sort-key=gen-02b +type=real + +[namelist:comorph=par_gen_mass_fac] +compulsory=true +description=Coefficient scaling the initiating mass sources (dimensionless) +help=Higher values could mean that the layer-cloud in unstable environments + =convects away more rapidly, leading to reduced layer-cloud in such + =environments. + =Recommended value 0.25 +!kind=default +ns=namelist/Science/UM Convection/CoMorph/genesis +range=0.01:1.0 +sort-key=gen-01a +type=real + +[namelist:comorph=par_gen_pert_fac] +compulsory=true +description=Scaling factor for CoMorph turbulent parcel initial perturbations +help=Sets the magnitude of turbulent contributions to parcel initial + =perturbations to T, q, u and v. +!kind=default +ns=namelist/Science/UM Convection/CoMorph/genesis +sort-key=gen-02a +type=real + +[namelist:comorph=par_gen_rhpert] +compulsory=true +description=Comorph neutrally-buoyant parcel initial RH perturbation (fractional) +help=Sets the magnitude of the "background" parcel initial moisture excess + = in Relative Humidity (decimal fraction). +!kind=default +ns=namelist/Science/UM Convection/CoMorph/genesis +sort-key=gen-02c +type=real + +[namelist:comorph=par_radius_evol_method] +compulsory=true +description=Method for height-evolution of the parcel radius +!enumeration=true +help=Controls how the parcel radius changes with height during ascent. +!kind=default +ns=namelist/Science/UM Convection/CoMorph/plume +sort-key=plume-04a +value-titles=Keep parcel radius constant with height, + =Parcel radius scales with volume change, + =Scales with volume change but radius cannot decrease, + =Scales with volume change but neglect the detrainment term +values='constant','volume','no_decrease','no_detrain' + +[namelist:comorph=par_radius_init_method] +compulsory=true +description=Method to use for CoMorph parcel initial radius scaling +!enumeration=true +help=Use this to apply dependence of the parcel initial radius on precip, + =boundary-layer moisture and vertical velocity. +!kind=default +ns=namelist/Science/UM Convection/CoMorph/genesis +sort-key=gen-03a +value-titles=Constant scaling factor, + =Linear dependence on precip rate, + =Linear dependence on precip/q, + =Linear dependence on precip/q and vertical velocity, + =Linear dependence on both precip rate and q +values='constant','linear_p','linear_prq','linear_prq_w','linear_p_q' + +[namelist:comorph=par_radius_knob] +compulsory=true +description=Tuning knob for Comorph parcel initial radius +help=Scales the parcel initial radius, which determines the entrainment rate. +!kind=double +ns=namelist/Science/UM Convection/CoMorph/genesis +sort-key=gen-03b +type=real + +[namelist:comorph=par_radius_knob_max] +compulsory=true +description=Maximum parcel radius scaling from precip dependence +help=Sets the upper bound for the precip-dependent parcel radius scaling. +!kind=default +ns=namelist/Science/UM Convection/CoMorph/genesis +sort-key=gen-03c +type=real + +[namelist:comorph=par_radius_ppn_max] +compulsory=true +description=Comorph precip rate at which max parcel radius scaling occurs. +help=Scales the parcel initial radius, which determines the entrainment rate. +!kind=default +ns=namelist/Science/UM Convection/CoMorph/genesis +sort-key=gen-03d +type=real + +[namelist:comorph=rain_area_min] +compulsory=true +description=Minimum allowed area fraction for convective precip +help=Sets the minimum area fraction for CoMorph precipitation. +!kind=default +ns=namelist/Science/UM Convection/CoMorph +sort-key=04a +type=real + +[namelist:comorph=resdep_precipramp] +compulsory=true +description=Include grid-length dependence in Comorph parcel radius precip ramp +help=Option to make the CoMorph parcel radius dependence on the + =precipitation rate a function of the grid size, and so allow for + =the ability to resolve higher precip rates at higher resolution. +!kind=default +ns=namelist/Science/UM Convection/CoMorph/genesis +sort-key=gen-03e +trigger=namelist:comorph=dx_ref: .true.; +type=logical + +[namelist:comorph=rho_rim] +compulsory=true +description=Density of rimed ice used for graupel / kg m-3 +help=Recommended value is 600.0. +!kind=default +ns=namelist/Science/UM Convection/CoMorph/micro +range=200.0:1000.0 +sort-key=micro-04d +type=real + +[namelist:comorph=tdep_n_cf] +compulsory=true +description=Temperature-dependent ice number concentration slope / K +help=Sets the temperature dependence of the in-plume ice number concentration. +!kind=default +ns=namelist/Science/UM Convection/CoMorph/micro +sort-key=micro-06b +type=real + +[namelist:comorph=tdep_n_cl] +compulsory=true +description=Temperature-dependent liquid-cloud number concentration slope / K +help=Sets the temperature dependence of the in-plume liquid-cloud number. +!kind=default +ns=namelist/Science/UM Convection/CoMorph/micro +sort-key=micro-06a +type=real + +[namelist:comorph=vent_factor] +compulsory=true +description=Coefficient scaling ventilation-driven hydrometeor exchange +help=Scales the extra turbulent exchange due to hydrometeor fall-speed + =ventilation. +!kind=default +ns=namelist/Science/UM Convection/CoMorph/micro +range=0.0:1.0 +sort-key=micro-04c +type=real + +[namelist:comorph=wind_w_buoy_fac] +compulsory=true +description=Tuning constant for buoyancy-dependent convective fraction +help=Assuming w = fac * sqrt(buoyancy * radius). +!kind=default +ns=namelist/Science/UM Convection/CoMorph +range=0.5:2.0 +sort-key=02b +type=real + +[namelist:comorph=wind_w_fac] +compulsory=true +description=Prescribed parcel ascent rate profile / m s-1 +help=Used when the vertical momentum equation is disabled. +!kind=default +ns=namelist/Science/UM Convection/CoMorph/micro +range=0.1:10.0 +sort-key=micro-03a +type=real + [namelist:convection] compulsory=true description=Convection parametrization @@ -65,7 +471,7 @@ description=Convection scheme help=Convection scheme to use, options are: =GregoryRowntree - UM 6A mass flux scheme =LambertLewis - UM other convection scheme 10 - =Comorph - ComorphA scheme + =Comorph - CoMorph scheme sort-key=Panel-01 trigger=namelist:convection=number_of_convection_substeps: 'gregory_rowntree'; =namelist:convection=cape_timescale: 'gregory_rowntree'; @@ -83,25 +489,10 @@ trigger=namelist:convection=number_of_convection_substeps: 'gregory_rowntree'; =namelist:convection=c_mass_sh: 'gregory_rowntree'; =namelist:convection=l_conv_prog_dq: 'gregory_rowntree', 'comorph'; =namelist:convection=l_conv_prog_dtheta: 'gregory_rowntree', 'comorph'; - =namelist:convection=par_gen_mass_fac: 'comorph'; - =namelist:convection=par_gen_rhpert: 'comorph'; - =namelist:convection=par_radius_ppn_max: 'comorph'; - =namelist:convection=resdep_precipramp: 'comorph'; + =namelist:comorph: 'comorph'; value-titles=GregoryRowntree, LambertLewis, Comorph values='gregory_rowntree','lambert_lewis','comorph' -[namelist:convection=dx_ref] -compulsory=true -description=Reference grid-length for Comorph A resolution-dependence / m -help=When using the resolution-dependent link between parcel radius - =and precipitation rate, dx_ref is the grid-length at which the - =namelist input value of par_radius_ppn_max applies. - =Typically set to 50000.0m (ie. 50km). -!kind=default -range=100.0:1000000.0 -sort-key=Panel-08 -type=real - [namelist:convection=efrac] compulsory=true description=Evaporative fraction of fluxes at mid-level cloud base @@ -218,35 +609,6 @@ range=0.0:5.0 sort-key=Panel-03d type=real -[namelist:convection=par_gen_mass_fac] -compulsory=true -description=Coefficient scaling the initiating mass sources (dimensionless) -help=Higher values could mean that the layer-cloud in unstable environments - =convects away more rapidly, leading to reduced layer-cloud in such - =environments. - =Recommended value 0.25 -!kind=default -range=0.01:1.0 -sort-key=Panel-04 -type=real - -[namelist:convection=par_gen_rhpert] -compulsory=true -description=Comorph neutrally-buoyant parcel initial RH perturbation (fractional) -help=Sets the magnitude of the "background" parcel initial moisture excess - = in Relative Humidity (decimal fraction). -!kind=default -sort-key=Panel-05 -type=real - -[namelist:convection=par_radius_ppn_max] -compulsory=true -description=Comorph precip rate at which max parcel radius scaling occurs. -help=Scales the parcel initial radius, which determines the entrainment rate. -!kind=default -sort-key=Panel-06 -type=real - [namelist:convection=prog_ent_grad] compulsory=true description=Gradient used to calculate entrainment scaling from prognostic field @@ -318,17 +680,6 @@ range=0:1 sort-key=Panel-03e type=real -[namelist:convection=resdep_precipramp] -compulsory=true -description=Include grid-length dependence in Comorph A parcel radius precip ramp -help=Option to make the CoMorph A parcel radius dependence on the - =precipitation rate a function of the grid size, and so allow for - =the ability to resolve higher precip rates at higher resolution. -!kind=default -sort-key=Panel-07 -trigger=namelist:convection=dx_ref: .true.; -type=logical - [namelist:convection=thpixs_mid] compulsory=true description=The initial mid-level parcel potential temperature perturbation (K) @@ -343,3 +694,19 @@ help=The initial mid-level temperature perturbation prior to any scaling. range=0.0:2.0 sort-key=Panel-03h type=real + +[ns=namelist/Science/UM Convection/CoMorph] +sort-key=Section-A06b +title=CoMorph Convection Scheme + +[ns=namelist/Science/UM Convection/CoMorph/genesis] +sort-key=Section-A06b-1 +title=Convective triggering and parcel initialisation + +[ns=namelist/Science/UM Convection/CoMorph/micro] +sort-key=Section-A06b-3 +title=Plume microphysics + +[ns=namelist/Science/UM Convection/CoMorph/plume] +sort-key=Section-A06b-2 +title=Plume model diff --git a/interfaces/physics_schemes_interface/rose-meta/um-convection/versions.py b/interfaces/physics_schemes_interface/rose-meta/um-convection/versions.py index 01e65a2b91..5304819ca1 100644 --- a/interfaces/physics_schemes_interface/rose-meta/um-convection/versions.py +++ b/interfaces/physics_schemes_interface/rose-meta/um-convection/versions.py @@ -31,3 +31,85 @@ def upgrade(self, config, meta_config=None): # Add settings return config, self.reports """ + +class vn32_t717(MacroUpgrade): + # Upgrade macro for PR#717 by Mike Whitall + + BEFORE_TAG = "vn3.2" + AFTER_TAG = "vn3.2_t717" + + def upgrade(self, config, meta_config=None): + # Add settings + + # Add new comorph namelist + # Append after 'convection' in configuration.nml + source = self.get_setting_value( + config, ["file:configuration.nml", "source"] + ) + source = re.sub( + r"(\(?)namelist:convection(\)?)(\n)", + r"\1namelist:convection\2\3 (namelist:comorph)\n", + source + ) + self.change_setting_value( + config, ["file:configuration.nml", "source"], source + ) + + # Move existing comorph namelist entries from the "convection" + # namelist to the new "comorph" namelist. + nml1 = "namelist:convection" + nml2 = "namelist:comorph" + for entry in ["par_gen_mass_fac", "par_gen_rhpert", + "par_radius_ppn_max", "resdep_precipramp", "dx_ref"]: + source = self.get_setting_value(config, [nml1, entry]) + self.remove_setting(config, [nml1, entry]) + self.add_setting(config, [nml2, entry], source) + + # Add new namelist entries with previous hardwired default values + nml = "namelist:comorph" + + # Top-level settings + self.add_setting(config, [nml, "n_dndraft_types"], "1") + self.add_setting(config, [nml, "cf_conv_fac"], "2.0") + self.add_setting(config, [nml, "wind_w_buoy_fac"], "1.0") + self.add_setting(config, [nml, "overlap_power"], "0.5") + self.add_setting(config, [nml, "rain_area_min"], "0.05") + + # Conv triggering and parcel initialisation + self.add_setting(config, [nml, "par_gen_pert_fac"], "0.333") + self.add_setting(config, [nml, "par_gen_core_fac"], "3.0") + self.add_setting(config, [nml, "par_radius_init_method"], + "'linear_p_q'") + self.add_setting(config, [nml, "par_radius_knob"], "0.45") + self.add_setting(config, [nml, "par_radius_knob_max"], "2.0") + self.add_setting(config, [nml, "ass_min_radius"], "500.0") + + # Plume model + self.add_setting(config, [nml, "ent_coef"], "0.2") + self.add_setting(config, [nml, "core_ent_cmr"], ".true.") + self.add_setting(config, [nml, "core_ent_fac"], "1.0") + self.add_setting(config, [nml, "min_cmr"], "2.0") + self.add_setting(config, [nml, "max_cmr"], "6.0") + self.add_setting(config, [nml, "drag_coef_par"], "0.5") + self.add_setting(config, [nml, "par_radius_evol_method"], + "'no_detrain'") + + # In-plume microphysics + self.add_setting(config, [nml, "autoc_opt"], "'quadratic'") + self.add_setting(config, [nml, "coef_auto"], "0.025") + self.add_setting(config, [nml, "hetnuc_temp"], "263.0") + self.add_setting(config, [nml, "cf_area_coef"], "10.0") + self.add_setting(config, [nml, "wind_w_fac"], "1.0") + self.add_setting(config, [nml, "col_eff_coef"], "1.0") + self.add_setting(config, [nml, "drag_coef_cond"], "0.5") + self.add_setting(config, [nml, "vent_factor"], "0.25") + self.add_setting(config, [nml, "rho_rim"], "600.0") + self.add_setting(config, [nml, "nconc_cl"], "1.0E8") + self.add_setting(config, [nml, "nconc_rain"], "1000.0") + self.add_setting(config, [nml, "nconc_cf"], "300.0") + self.add_setting(config, [nml, "nconc_snow"], "300.0") + self.add_setting(config, [nml, "nconc_graup"], "100.0") + self.add_setting(config, [nml, "tdep_n_cl"], "0.0") + self.add_setting(config, [nml, "tdep_n_cf"], "8.18") + + return config, self.reports 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 b3e2e73252..41d9962d17 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 @@ -151,13 +151,61 @@ module um_physics_init_mod c_mass_sh_in => c_mass_sh, & l_conv_prog_dtheta_in => l_conv_prog_dtheta, & l_conv_prog_dq_in => l_conv_prog_dq, & - par_gen_mass_fac_in => par_gen_mass_fac, & - par_gen_rhpert_in => par_gen_rhpert, & - par_radius_ppn_max_in => par_radius_ppn_max, & - resdep_precipramp, dx_ref_in => dx_ref, & l_cvdiag_ctop_qmax_in => l_cvdiag_ctop_qmax, & llcs_first_outer + use comorph_config_mod, only : & + ass_min_radius_in => ass_min_radius, & + autoc_opt_in => autoc_opt, & + autoc_opt_linear, & + autoc_opt_quadratic, & + cf_area_coef_in => cf_area_coef, & + cf_conv_fac_in => cf_conv_fac, & + coef_auto_in => coef_auto, & + col_eff_coef_in => col_eff_coef, & + core_ent_cmr_in => core_ent_cmr, & + core_ent_fac_in => core_ent_fac, & + drag_coef_cond_in => drag_coef_cond, & + drag_coef_par_in => drag_coef_par, & + dx_ref_in => dx_ref, & + ent_coef_in => ent_coef, & + hetnuc_temp_in => hetnuc_temp, & + max_cmr_in => max_cmr, & + min_cmr_in => min_cmr, & + n_dndraft_types_in => n_dndraft_types, & + nconc_cf_in => nconc_cf, & + nconc_cl_in => nconc_cl, & + nconc_graup_in => nconc_graup, & + nconc_rain_in => nconc_rain, & + nconc_snow_in => nconc_snow, & + overlap_power_in => overlap_power, & + par_gen_core_fac_in => par_gen_core_fac, & + par_gen_mass_fac_in => par_gen_mass_fac, & + par_gen_pert_fac_in => par_gen_pert_fac, & + par_gen_rhpert_in => par_gen_rhpert, & + par_radius_evol_method_in => par_radius_evol_method, & + par_radius_evol_method_constant, & + par_radius_evol_method_volume, & + par_radius_evol_method_no_decrease, & + par_radius_evol_method_no_detrain, & + par_radius_init_method_in => par_radius_init_method, & + par_radius_init_method_constant, & + par_radius_init_method_linear_p, & + par_radius_init_method_linear_prq, & + par_radius_init_method_linear_prq_w, & + par_radius_init_method_linear_p_q, & + par_radius_knob_in => par_radius_knob, & + par_radius_knob_max_in => par_radius_knob_max, & + par_radius_ppn_max_in => par_radius_ppn_max, & + rain_area_min_in => rain_area_min, & + rho_rim_in => rho_rim, & + resdep_precipramp_in => resdep_precipramp, & + tdep_n_cf_in => tdep_n_cf, & + tdep_n_cl_in => tdep_n_cl, & + vent_factor_in => vent_factor, & + wind_w_buoy_fac_in => wind_w_buoy_fac, & + wind_w_fac_in => wind_w_fac + use extrusion_config_mod, only : domain_height, number_of_layers use formulation_config_mod, only : moisture_formulation, & @@ -386,15 +434,27 @@ subroutine um_physics_init() l_pc2_homog_conv_pressure, l_cloud_call_b4_conv, & i_bm_ez_orig, i_bm_ez_subcrit, i_bm_ez_entpar use cloud_config_mod, only: cld_fsd_hill - use comorph_um_namelist_mod, only: ass_min_radius, autoc_opt, & - cf_conv_fac, coef_auto, col_eff_coef, core_ent_fac, drag_coef_cond, & - drag_coef_par, ent_coef, hetnuc_temp, l_core_ent_cmr, & - n_dndraft_types, overlap_power, par_gen_core_fac, par_gen_mass_fac, & - par_gen_pert_fac, par_gen_rhpert, par_radius_evol_method, & - par_radius_init_method, par_radius_knob, par_radius_knob_max, & - par_radius_ppn_max, r_fac_tdep_n, rain_area_min, rho_rim, & - vent_factor, wind_w_buoy_fac, wind_w_fac, check_run_comorph, & - l_resdep_precipramp, dx_ref + use comorph_um_namelist_mod, only: & + ! UM namelist entries + ass_min_radius, autoc_opt, cf_area_coef, cf_conv_fac, coef_auto, & + col_eff_coef, core_ent_fac, drag_coef_cond, drag_coef_par, dx_ref, & + ent_coef, hetnuc_temp, l_core_ent_cmr, l_resdep_precipramp, & + max_cmr, min_cmr, n_dndraft_types, & + nconc_cf, nconc_cl, nconc_graup, nconc_rain, nconc_snow, & + overlap_power, par_gen_core_fac, par_gen_mass_fac, & + par_gen_pert_fac, par_gen_rhpert, par_radius_evol_method, & + par_radius_init_method, par_radius_knob, par_radius_knob_max, & + par_radius_ppn_max, rain_area_min, rho_rim, tdep_n_cf, tdep_n_cl, & + vent_factor, wind_w_buoy_fac, wind_w_fac, & + ! Namelis checking subroutine + check_run_comorph, & + ! Allowed values for multi-option switches + no_dependence, rain_dependence, qfacrain_dependence, w_dependence, & + linear_qfacrain_dep + use comorph_constants_mod, only: & + autoc_linear, autoc_quadratic, & + par_radius_evol_const, par_radius_evol_volume, & + par_radius_evol_no_decrease, par_radius_evol_no_detrain use cv_run_mod, only: icvdiag, cvdiag_inv, cvdiag_sh_wtest, & limit_pert_opt, tv1_sd_opt, iconv_congestus, iconv_deep, & ent_fac_dp, cldbase_opt_dp, cldbase_opt_sh, w_cape_limit, & @@ -875,48 +935,90 @@ subroutine um_physics_init() tau_conv_prog_dtheta = 2700.0_r_um tau_conv_prog_dq = 2700.0_r_um - ! main Comorph options - ass_min_radius = 500.0_r_um - autoc_opt = 2 - cf_conv_fac = 2.0_r_um - coef_auto = 0.025_r_um - col_eff_coef = 1.0_r_um - core_ent_fac = 1.0_r_um - drag_coef_cond = 0.5_r_um - drag_coef_par = 0.5_r_um - dx_ref = dx_ref_in - ent_coef = 0.2_r_um - hetnuc_temp = 263.0_r_um - l_core_ent_cmr = .true. - l_resdep_precipramp = resdep_precipramp - n_dndraft_types = 1 - overlap_power = 0.5_r_um - par_gen_core_fac = 3.0_r_um - par_gen_mass_fac = par_gen_mass_fac_in - par_gen_pert_fac = 0.333_r_um - par_gen_rhpert = par_gen_rhpert_in - par_radius_evol_method = 3 - par_radius_init_method = 4 - par_radius_knob = 0.45_r_um - par_radius_knob_max = 2.0_r_um - par_radius_ppn_max = par_radius_ppn_max_in - r_fac_tdep_n = 8.18_r_um - rain_area_min = 0.05_r_um - rho_rim = 600.0_r_um - vent_factor = 0.25_r_um - wind_w_buoy_fac = 1.0_r_um - wind_w_fac = 1.0_r_um + ! CoMorph namelist settings... + + ! Top-level settings: + n_dndraft_types = n_dndraft_types_in + cf_conv_fac = real( cf_conv_fac_in, r_um ) + wind_w_buoy_fac = real( wind_w_buoy_fac_in, r_um ) + overlap_power = real( overlap_power_in, r_um ) + rain_area_min = real( rain_area_min_in, r_um ) + + ! Conv triggering and parcel initialisation + par_gen_mass_fac = real( par_gen_mass_fac_in, r_um ) + par_gen_pert_fac = real( par_gen_pert_fac_in, r_um ) + par_gen_core_fac = real( par_gen_core_fac_in, r_um ) + par_gen_rhpert = real( par_gen_rhpert_in, r_um ) + select case ( par_radius_init_method_in ) + case ( par_radius_init_method_constant ) + par_radius_init_method = no_dependence + case ( par_radius_init_method_linear_p ) + par_radius_init_method = rain_dependence + case ( par_radius_init_method_linear_prq ) + par_radius_init_method = qfacrain_dependence + case ( par_radius_init_method_linear_prq_w ) + par_radius_init_method = w_dependence + case ( par_radius_init_method_linear_p_q ) + par_radius_init_method = linear_qfacrain_dep + end select + par_radius_knob = real( par_radius_knob_in, r_um ) + par_radius_knob_max = real( par_radius_knob_max_in, r_um ) + par_radius_ppn_max = real( par_radius_ppn_max_in, r_um ) + l_resdep_precipramp = resdep_precipramp_in + dx_ref = real( dx_ref_in, r_um ) + ass_min_radius = real( ass_min_radius_in, r_um ) + + ! Plume model + ent_coef = real( ent_coef_in, r_um ) + l_core_ent_cmr = core_ent_cmr_in + core_ent_fac = real( core_ent_fac_in, r_um ) + min_cmr = real( min_cmr_in, r_um ) + max_cmr = real( max_cmr_in, r_um ) + drag_coef_par = real( drag_coef_par_in, r_um ) + select case ( par_radius_evol_method_in ) + case ( par_radius_evol_method_constant ) + par_radius_evol_method = par_radius_evol_const + case ( par_radius_evol_method_volume ) + par_radius_evol_method = par_radius_evol_volume + case ( par_radius_evol_method_no_decrease ) + par_radius_evol_method = par_radius_evol_no_decrease + case ( par_radius_evol_method_no_detrain ) + par_radius_evol_method = par_radius_evol_no_detrain + end select + + ! In-plume microphysics + select case ( autoc_opt_in ) + case ( autoc_opt_linear ) + autoc_opt = autoc_linear + case (autoc_opt_quadratic ) + autoc_opt = autoc_quadratic + end select + coef_auto = real( coef_auto_in, r_um ) + hetnuc_temp = real( hetnuc_temp_in, r_um ) + cf_area_coef = real( cf_area_coef_in, r_um ) + wind_w_fac = real( wind_w_fac_in, r_um ) + col_eff_coef = real( col_eff_coef_in, r_um ) + drag_coef_cond = real( drag_coef_cond_in, r_um ) + vent_factor = real( vent_factor_in, r_um ) + rho_rim = real( rho_rim_in, r_um ) + nconc_cl = real( nconc_cl_in, r_um ) + nconc_rain = real( nconc_rain_in, r_um ) + nconc_cf = real( nconc_cf_in, r_um ) + nconc_snow = real( nconc_snow_in, r_um ) + nconc_graup = real( nconc_graup_in, r_um ) + tdep_n_cl = real( tdep_n_cl_in, r_um ) + tdep_n_cf = real( tdep_n_cf_in, r_um ) ! check the namelist call check_run_comorph() case(cv_scheme_gregory_rowntree) - if ( boundary_layer /= boundary_layer_um ) then - write( log_scratch_space, '(A)' ) & - 'UM boundary layer is required for GR convection - please switch on' - call log_event( log_scratch_space, LOG_LEVEL_ERROR ) - end if + if ( boundary_layer /= boundary_layer_um ) then + write( log_scratch_space, '(A)' ) & + 'UM boundary layer is required for GR convection - please switch on' + call log_event( log_scratch_space, LOG_LEVEL_ERROR ) + end if i_convection_vn = i_convection_vn_6a adapt = 8 diff --git a/science/linear/integration-test/nwp_gal9/resources/nwp_gal9_configuration.nml b/science/linear/integration-test/nwp_gal9/resources/nwp_gal9_configuration.nml index 456d2e84c2..06e0f1c287 100644 --- a/science/linear/integration-test/nwp_gal9/resources/nwp_gal9_configuration.nml +++ b/science/linear/integration-test/nwp_gal9/resources/nwp_gal9_configuration.nml @@ -30,10 +30,8 @@ stochastic_physics='none', surface='none', / &convection -dx_ref=50000.0, l_cvdiag_ctop_qmax=.false., qlmin=4.0e-4, -resdep_precipramp=.false., / &cosp l_cosp=.false., diff --git a/science/physics_schemes/source/convection/comorph/control/cloudfracs_type_mod.F90 b/science/physics_schemes/source/convection/comorph/control/cloudfracs_type_mod.F90 index 0972fdf894..37f59ecb69 100644 --- a/science/physics_schemes/source/convection/comorph/control/cloudfracs_type_mod.F90 +++ b/science/physics_schemes/source/convection/comorph/control/cloudfracs_type_mod.F90 @@ -9,7 +9,7 @@ module cloudfracs_type_mod -use comorph_constants_mod, only: real_hmprec, name_length +use comorph_constants_mod, only: real_cvprec, real_hmprec, name_length use fields_type_mod, only: fields_list_type implicit none @@ -86,6 +86,9 @@ module cloudfracs_type_mod ! (used for labelling diagnostics and error print-outs) character(len=name_length), allocatable :: convcloud_names(:) +! Max allowed value for each concloud field (all must be > 0) +real(kind=real_cvprec), allocatable :: convcloud_max(:) + contains @@ -94,8 +97,9 @@ module cloudfracs_type_mod ! depending on switches subroutine cloudfracs_set_addresses() -use comorph_constants_mod, only: i_convcloud, i_convcloud_bulkonly, & - i_convcloud_liqonly, i_convcloud_mph +use comorph_constants_mod, only: real_cvprec, one, i_convcloud, & + i_convcloud_bulkonly, i_convcloud_liqonly, & + i_convcloud_mph implicit none @@ -111,9 +115,6 @@ subroutine cloudfracs_set_addresses() n_convcloud = 2 i_frac_bulk_conv = 1 i_q_c_conv = 2 - allocate( convcloud_names(n_convcloud) ) - convcloud_names(i_frac_bulk_conv) = "frac_bulk_conv" - convcloud_names(i_q_c_conv) = "q_c_conv" case ( i_convcloud_liqonly ) ! Only using liquid convective cloud ! (but we also still output the bulk convective cloud fraction, @@ -123,10 +124,6 @@ subroutine cloudfracs_set_addresses() i_frac_liq_conv = 1 i_frac_bulk_conv = 2 i_q_cl_conv = 3 - allocate( convcloud_names(n_convcloud) ) - convcloud_names(i_frac_liq_conv) = "frac_liq_conv" - convcloud_names(i_frac_bulk_conv) = "frac_bulk_conv" - convcloud_names(i_q_cl_conv) = "q_cl_conv" case ( i_convcloud_mph ) ! Using separate liquid and ice cloud, with variable overlap n_convcloud = 5 @@ -135,14 +132,26 @@ subroutine cloudfracs_set_addresses() i_frac_bulk_conv = 3 i_q_cl_conv = 4 i_q_cf_conv = 5 - allocate( convcloud_names(n_convcloud) ) - convcloud_names(i_frac_liq_conv) = "frac_liq_conv" - convcloud_names(i_frac_ice_conv) = "frac_ice_conv" - convcloud_names(i_frac_bulk_conv) = "frac_bulk_conv" - convcloud_names(i_q_cl_conv) = "q_cl_conv" - convcloud_names(i_q_cf_conv) = "q_cf_conv" end select +! Set names of the fields +allocate( convcloud_names(n_convcloud) ) +if (i_frac_liq_conv>0) convcloud_names(i_frac_liq_conv) = "frac_liq_conv" +if (i_frac_ice_conv>0) convcloud_names(i_frac_ice_conv) = "frac_ice_conv" +if (i_frac_bulk_conv>0) convcloud_names(i_frac_bulk_conv) = "frac_bulk_conv" +if (i_q_cl_conv>0) convcloud_names(i_q_cl_conv) = "q_cl_conv" +if (i_q_cf_conv>0) convcloud_names(i_q_cf_conv) = "q_cf_conv" +if (i_q_c_conv>0) convcloud_names(i_q_c_conv) = "q_c_conv" + +! Set max allowed values used in run-time checks +allocate( convcloud_max(n_convcloud) ) +if (i_frac_liq_conv>0) convcloud_max(i_frac_liq_conv) = one +if (i_frac_ice_conv>0) convcloud_max(i_frac_ice_conv) = one +if (i_frac_bulk_conv>0) convcloud_max(i_frac_bulk_conv) = one +if (i_q_cl_conv>0) convcloud_max(i_q_cl_conv) = 0.1_real_cvprec +if (i_q_cf_conv>0) convcloud_max(i_q_cf_conv) = 0.1_real_cvprec +if (i_q_c_conv>0) convcloud_max(i_q_c_conv) = 0.1_real_cvprec + return end subroutine cloudfracs_set_addresses @@ -282,7 +291,7 @@ end subroutine cloudfracs_list_clear subroutine cloudfracs_check_bad_values( cloudfracs, & where_string ) -use comorph_constants_mod, only: name_length, l_cv_cloudfrac +use comorph_constants_mod, only: zero, one, name_length, l_cv_cloudfrac use check_bad_values_mod, only: check_bad_values_3d implicit none @@ -300,10 +309,6 @@ subroutine cloudfracs_check_bad_values( cloudfracs, & ! Lower and upper bounds of array integer :: lb(3), ub(3) -! Flag passed into check_bad_values; -! all fields checked are positive-only so hardwire to true -logical, parameter :: l_positive = .true. - ! Loop counter integer :: i_field @@ -317,21 +322,21 @@ subroutine cloudfracs_check_bad_values( cloudfracs, & ub = ubound( cloudfracs % frac_liq ) call check_bad_values_3d( lb, ub, cloudfracs % frac_liq, & where_string, field_name, & - l_positive ) + field_min=zero, field_max=one ) field_name = "frac_ice" lb = lbound( cloudfracs % frac_ice ) ub = ubound( cloudfracs % frac_ice ) call check_bad_values_3d( lb, ub, cloudfracs % frac_ice, & where_string, field_name, & - l_positive ) + field_min=zero, field_max=one ) field_name = "frac_bulk" lb = lbound( cloudfracs % frac_bulk ) ub = ubound( cloudfracs % frac_bulk ) call check_bad_values_3d( lb, ub, cloudfracs % frac_bulk, & where_string, field_name, & - l_positive ) + field_min=zero, field_max=one ) end if @@ -341,7 +346,7 @@ subroutine cloudfracs_check_bad_values( cloudfracs, & ub = ubound( cloudfracs % frac_precip ) call check_bad_values_3d( lb, ub, cloudfracs % frac_precip, & where_string, field_name, & - l_positive ) + field_min=zero, field_max=one ) if ( n_convcloud > 0 ) then ! Check convective cloud fields if they are used @@ -350,7 +355,8 @@ subroutine cloudfracs_check_bad_values( cloudfracs, & ub = ubound( cloudfracs % convcloud_list(i_field)%pt ) call check_bad_values_3d( lb, ub, cloudfracs % convcloud_list(i_field)%pt, & where_string, convcloud_names(i_field), & - l_positive ) + field_min=zero, & + field_max=convcloud_max(i_field) ) end do end if diff --git a/science/physics_schemes/source/convection/comorph/control/comorph_constants_mod.F90 b/science/physics_schemes/source/convection/comorph/control/comorph_constants_mod.F90 index 59f7315cd4..5088713096 100644 --- a/science/physics_schemes/source/convection/comorph/control/comorph_constants_mod.F90 +++ b/science/physics_schemes/source/convection/comorph/control/comorph_constants_mod.F90 @@ -96,6 +96,9 @@ module comorph_constants_mod ! Square-root of the above real(kind=real_cvprec), parameter :: sqrt_min_delta = sqrt(min_delta) +! Largest possible floating point number +real(kind=real_cvprec), parameter :: max_float = huge(zero) + !--------------------------------------------------------------- ! Things set at run-time by the host-model @@ -233,7 +236,7 @@ module comorph_constants_mod ! GENERAL CONTROL OPTIONS... ! Number of updraft types for cloud spectrum model option -integer, parameter :: n_updraft_types = 1 +integer :: n_updraft_types = 1 ! Number of independent downdraft types ! (distinct from "fall-back downdrafts", which are each tied @@ -451,11 +454,13 @@ module comorph_constants_mod ! Number concentration of hydrometeor per unit dry-mass / kg-1 real(kind=real_cvprec) :: n - ! Flag for temperature-dependent number concentration. - ! If true, n above is taken to be the number concentration - ! at the melting point, with a parameterised increase - ! at colder temperatures. - logical :: l_tdep_n + ! Temperature-dependence of the number concentration. + ! If non-zero, n above is taken to be the number concentration + ! at the melting point, with a parameterised increase or decrease + ! at colder temperatures. The number concentration n(T) will be given by: + ! n(T) = n0 exp( fac_tdep_n ( T - Tmelt ) ) + ! (but limited above Tmelt and below T_homnuc) + real(kind=real_cvprec) :: fac_tdep_n ! Effective area coefficient for hydrometeors; ratio of actual ! area to the area you get by assuming a sphere @@ -499,7 +504,7 @@ module comorph_constants_mod cond_name = "cl", & l_ice = .false., & n = 1.0E8_real_cvprec, & ! Liq cloud number conc. - l_tdep_n = .false., & ! Use T-dependent number + fac_tdep_n= zero, & ! Not using T-dependent number area_coef = 1.0_real_cvprec, & ! 1.0 for spheres r_min = 1.0e-6_real_cvprec, & ! 1 micron CCN size i_sg = i_sg_frac_liq, & ! Lives in the liquid cloud fraction @@ -510,7 +515,7 @@ module comorph_constants_mod cond_name = "rain", & l_ice = .false., & n = 1000.0_real_cvprec, & ! Rain number conc. ~ 1 per l - l_tdep_n = .false., & ! Use T-dependent number + fac_tdep_n= zero, & ! Not using T-dependent number area_coef = 1.0_real_cvprec, & ! 1.0 for spheres r_min = 0.0_real_cvprec, & ! No CCN for rain i_sg = i_sg_frac_prec, & ! Lives in the precip fraction @@ -521,7 +526,7 @@ module comorph_constants_mod cond_name = "cf", & l_ice = .true., & n = 300.0_real_cvprec, & ! Ice cloud number at 0oC - l_tdep_n = .true., & ! Use T-dependent number + fac_tdep_n= -one/8.18_real_cvprec, & ! Use T-dependent number area_coef = 10.0_real_cvprec, & ! 10.0 for crystals r_min = 5.0e-6_real_cvprec, & ! 5 micron CCN size i_sg = i_sg_frac_ice, & ! Lives in the ice cloud fraction @@ -532,7 +537,7 @@ module comorph_constants_mod cond_name = "snow", & l_ice = .true., & n = 300.0_real_cvprec, & ! Snow number conc. ~0.3 per l - l_tdep_n = .false., & ! T-depentent number conc. off. + fac_tdep_n= zero, & ! T-dependent number conc. off. area_coef = 10.0_real_cvprec, & ! 10.0 for aggregates r_min = 0.0_real_cvprec, & ! No CCN for snow i_sg = i_sg_frac_ice, & ! Lives in the ice cloud fraction @@ -543,22 +548,15 @@ module comorph_constants_mod cond_name = "graup", & l_ice = .true., & n = 100.0_real_cvprec, & ! Graupel number conc. ~0.1 per l - l_tdep_n = .false., & ! T-depentent number conc. off. + fac_tdep_n= zero, & ! T-dependent number conc. off. area_coef = 1.0_real_cvprec, & ! 1.0 for spheres r_min = 0.0_real_cvprec, & ! No CCN for graupel - i_sg = i_sg_homog, & ! Assumed homogeneous across grid-box + i_sg = i_sg_frac_prec, & ! Lives in the precip fraction i_frzmlt = i_rain ) ! Melt into rain ! Density of rimed ice (used for graupel) real(kind=real_cvprec) :: rho_rim = 600.0_real_cvprec -! Temperature-dependent ice number concentration slope -! The number concentration n(T) will be given by: -! n(T) = n0 exp( fac_tdep_n ( T - Tmelt ) ) -! ( but limited above Tmelt and below T_homnuc) -real(kind=real_cvprec) :: fac_tdep_n = -one/8.18_real_cvprec -! set to 1/8.18 K-1, consistent with the Wilson-Ballard microphysics. - ! Ice nucleation ! Homogeneous freezing temperature / K ! All liquid is instantly frozen below this @@ -703,8 +701,6 @@ module comorph_constants_mod ! a) Those used if NOT using turbulence-based parcel properties ! (l_turb_par_gen = .FALSE.) -! Prescribed initial parcel radius -real(kind=real_cvprec), parameter :: par_gen_radius = 500.0_real_cvprec ! Fractional moisture perturbation to apply to the initial parcel real(kind=real_cvprec), parameter :: par_gen_qpert = 0.05_real_cvprec @@ -745,8 +741,8 @@ module comorph_constants_mod ! Minimum and maximum allowed values of the core/mean ratio of parcel ! buoyancies, which sets the power of the assumed power-law PDF of ! in-parcel buoyancy used for detrainment -real(kind=real_cvprec), parameter :: min_cmr = 2.0_real_cvprec -real(kind=real_cvprec), parameter :: max_cmr = 6.0_real_cvprec +real(kind=real_cvprec) :: min_cmr = 2.0_real_cvprec +real(kind=real_cvprec) :: max_cmr = 6.0_real_cvprec ! Maximum allowed convective area fraction real(kind=real_cvprec), parameter :: max_sigma = 0.5_real_cvprec diff --git a/science/physics_schemes/source/convection/comorph/control/comorph_ctl.F90 b/science/physics_schemes/source/convection/comorph/control/comorph_ctl.F90 index cfbd71e51d..6ffa116695 100644 --- a/science/physics_schemes/source/convection/comorph/control/comorph_ctl.F90 +++ b/science/physics_schemes/source/convection/comorph/control/comorph_ctl.F90 @@ -31,7 +31,7 @@ subroutine comorph_ctl( l_tracer, n_segments, & use fields_type_mod, only: fields_type, l_init_fields_type_mod, & fields_set_addresses, & fields_list_make, fields_list_clear, & - field_names, field_positive, n_fields + field_names, field_min, field_max, n_fields use grid_type_mod, only: grid_type, grid_check_bad_values use turb_type_mod, only: turb_type, & turb_list_make, turb_list_clear, & @@ -242,7 +242,8 @@ subroutine comorph_ctl( l_tracer, n_segments, & ub_1 = ubound( fields_n % list(i_field)%pt ) call check_bad_values_3d( lb_1, ub_1, fields_n%list(i_field)%pt, & where_string, field_names(i_field), & - field_positive(i_field) ) + field_min=field_min(i_field), & + field_max=field_max(i_field) ) end do ! Check latest fields: where_string = "On input to CoMorph: latest fields:" @@ -251,7 +252,8 @@ subroutine comorph_ctl( l_tracer, n_segments, & ub_1 = ubound( fields_np1 % list(i_field)%pt ) call check_bad_values_3d( lb_1, ub_1, fields_np1%list(i_field)%pt, & where_string, field_names(i_field), & - field_positive(i_field) ) + field_min=field_min(i_field), & + field_max=field_max(i_field) ) end do ! Check input grid fields @@ -514,7 +516,8 @@ subroutine comorph_ctl( l_tracer, n_segments, & ub_1 = ubound( fields_np1 % list(i_field)%pt ) call check_bad_values_3d( lb_1, ub_1, fields_np1%list(i_field)%pt, & where_string, field_names(i_field), & - field_positive(i_field) ) + field_min=field_min(i_field), & + field_max=field_max(i_field) ) end do ! Check cloud and rain fractions (the convective cloud fields diff --git a/science/physics_schemes/source/convection/comorph/control/conv_genesis_ctl.F90 b/science/physics_schemes/source/convection/comorph/control/conv_genesis_ctl.F90 index 06f9673d43..ddaab63f4e 100644 --- a/science/physics_schemes/source/convection/comorph/control/conv_genesis_ctl.F90 +++ b/science/physics_schemes/source/convection/comorph/control/conv_genesis_ctl.F90 @@ -315,7 +315,6 @@ subroutine conv_genesis_ctl( max_points, ij_first, ij_last, & k_half=k2, grid_half_super=grid_half(:,:,k2_c) ) end do - ! Compress grid-mean primary fields from current level ! and the levels above and below do i_field = 1, n_fields_tot diff --git a/science/physics_schemes/source/convection/comorph/control/core_combine.F90 b/science/physics_schemes/source/convection/comorph/control/core_combine.F90 new file mode 100644 index 0000000000..667a7e4a22 --- /dev/null +++ b/science/physics_schemes/source/convection/comorph/control/core_combine.F90 @@ -0,0 +1,157 @@ +! *****************************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******************************* + +! Code Owner: Please refer to the UM file CodeOwners.txt +! This file belongs in section: convection_comorph + +module core_combine_mod + +implicit none + +contains + +! Subroutine to calculate combined parcel core properties of 2 different +! parcels; there are various different options for how to do this. +! This needs to be done when: +! a) Combining initiating parcels from multiple different sub-grid +! regions on a given model-level. +! b) Combining initiating parcels from subsequent model-levels higher-up +! into the rising plume. +subroutine core_combine( n_points_a, n_points_m, index_ic, & + n_points_super_a, n_points_super_m, & + i_field_first, i_field_last, l_down, & + core_a_fields, core_m_fields, & + edge_a_virt_temp, edge_m_virt_temp ) + +use comorph_constants_mod, only: real_cvprec, zero, one +use fields_type_mod, only: i_temperature, i_q_vap, i_qc_first, i_qc_last +use calc_virt_temp_mod, only: calc_virt_temp + +implicit none + +! Number of points in the parcel to be added "_a" +integer, intent(in) :: n_points_a + +! Number of points in the existing parcel arrays "_m" +integer, intent(in) :: n_points_m + +! Indices for referencing the "_m" arrays from the "_a" compression list +integer, intent(in) :: index_ic(n_points_a) + +! Sizes of the fields super-arrays; maybe > n_points due to reusing arrays +integer, intent(in) :: n_points_super_a +integer, intent(in) :: n_points_super_m + +! First and last primary fields to set +integer, intent(in) :: i_field_first +integer, intent(in) :: i_field_last + +! Flag for downdraft versus updraft +logical, intent(in) :: l_down + +! Parcel core fields to be combined: +! "_a": Input properties of the "added" parcel (one of the parcels to combine) +real(kind=real_cvprec), intent(in) :: core_a_fields & + ( n_points_super_a, i_field_first:i_field_last ) +! "_m": IN: properties of the other parcel to combine "_a" into +! OUT: properties of the "merged" parcel after combining with "_a" +real(kind=real_cvprec), intent(in out) :: core_m_fields & + ( n_points_super_m, i_field_first:i_field_last ) + +! Parcel edge virtual temperatures to be combined +real(kind=real_cvprec), intent(in) :: edge_a_virt_temp(n_points_a) +real(kind=real_cvprec), intent(in out) :: edge_m_virt_temp(n_points_m) + +! Weights for combining the core properties +real(kind=real_cvprec) :: weight_core_a(n_points_a) +real(kind=real_cvprec) :: weight_core_m(n_points_a) + +! Virtual temperatures of the core properties +real(kind=real_cvprec) :: core_a_virt_temp(n_points_a) +real(kind=real_cvprec) :: core_m_virt_temp(n_points_m) + +! Loop counters +integer :: ic, ic2, i_field + + +! Compute the core virtual temperature of the two parcels +! NOTE: this calculation relies on the fact that the parcel +! core properties are NOT in conserved variable form. +call calc_virt_temp( n_points_a, n_points_super_a, & + core_a_fields(:,i_temperature), & + core_a_fields(:,i_q_vap), & + core_a_fields(:,i_qc_first:i_qc_last), & + core_a_virt_temp ) +call calc_virt_temp( n_points_m, n_points_super_m, & + core_m_fields(:,i_temperature), & + core_m_fields(:,i_q_vap), & + core_m_fields(:,i_qc_first:i_qc_last), & + core_m_virt_temp ) + +! Choose properties from the parcel with the more buoyant core. + +! Reset the weights so that the core fields will inherit only +! the values from the most buoyant of the two. +! Combine the edge virtual temperatures by choosing the least buoyant edge +! (i.e. we always try to make the PDF of Tv as wide as possible) +if ( l_down ) then + do ic = 1, n_points_a + ic2 = index_ic(ic) + ! Choose most negatively buoyant core for downdrafts + ! TEMPORARY CODE TO PRESERVE KGO: should really test on mass-fluxes > 0 + ! (this code can use core properties from parcel m with zero mass-flux, + ! which is wrong; fix this soon...) + if ( core_a_virt_temp(ic) <= core_m_virt_temp(ic2) .or. & + ( .not. core_m_virt_temp(ic2) > zero ) ) then + weight_core_a(ic) = one + weight_core_m(ic) = zero + else + weight_core_a(ic) = zero + weight_core_m(ic) = one + end if + ! Choose least negatively buoyant edge for downdrafts + if ( edge_a_virt_temp(ic) > edge_m_virt_temp(ic2) .or. & + ( .not. edge_m_virt_temp(ic2) > zero ) ) then + edge_m_virt_temp(ic2) = edge_a_virt_temp(ic) + end if + end do +else + do ic = 1, n_points_a + ic2 = index_ic(ic) + ! Choose most positively buoyant core for updrafts + ! TEMPORARY CODE TO PRESERVE KGO: should really test on mass-fluxes > 0 + ! (this code can use core properties from parcel m with zero mass-flux, + ! which is wrong; fix this soon...) + if ( core_a_virt_temp(ic) >= core_m_virt_temp(ic2) .or. & + ( .not. core_m_virt_temp(ic2) > zero ) ) then + weight_core_a(ic) = one + weight_core_m(ic) = zero + else + weight_core_a(ic) = zero + weight_core_m(ic) = one + end if + ! Choose least positively buoyant edge for downdrafts + if ( edge_a_virt_temp(ic) < edge_m_virt_temp(ic2) .or. & + ( .not. edge_m_virt_temp(ic2) > zero ) ) then + edge_m_virt_temp(ic2) = edge_a_virt_temp(ic) + end if + end do +end if + +! Compute combined parcel core properties using the weights set above +do i_field = i_field_first, i_field_last + do ic = 1, n_points_a + ic2 = index_ic(ic) + core_m_fields(ic2,i_field) = weight_core_m(ic) * core_m_fields(ic2,i_field)& + + weight_core_a(ic) * core_a_fields(ic,i_field) + end do +end do + + +return +end subroutine core_combine + +end module core_combine_mod diff --git a/science/physics_schemes/source/convection/comorph/control/diags_super_type_mod.F90 b/science/physics_schemes/source/convection/comorph/control/diags_super_type_mod.F90 index 789116a5bd..a4fe1a9d48 100644 --- a/science/physics_schemes/source/convection/comorph/control/diags_super_type_mod.F90 +++ b/science/physics_schemes/source/convection/comorph/control/diags_super_type_mod.F90 @@ -756,7 +756,7 @@ end subroutine diags_super_expand subroutine diags_super_combine( n_diags_super, l_weight, index_ic, & diags_super_a, diags_super_m ) -use comorph_constants_mod, only: real_cvprec, one, min_float +use comorph_constants_mod, only: real_cvprec, zero, one, min_float implicit none @@ -791,18 +791,17 @@ subroutine diags_super_combine( n_diags_super, l_weight, index_ic, & do ic = 1, diags_super_a % cmpr % n_points ic2 = index_ic(ic) - ! Add contribution from diags_super_a to the combined total mass - diags_super_m % super(ic2,i_ds) = diags_super_m % super(ic2,i_ds) & - + diags_super_a % super(ic,i_ds) - ! Calculate mass-fraction weight to be applied to subsequent diagnostics - weight_a(ic) = diags_super_a % super(ic,i_ds) & - / max( diags_super_m % super(ic2,i_ds), min_float ) - end do - do ic = 1, diags_super_a % cmpr % n_points - ic2 = index_ic(ic) - ! Force weight to be exactly 1.0 at points where the mass-fluxes are - ! equal (needed for reproducibility on different decompositions). - if ( diags_super_a%super(ic,i_ds) == diags_super_m%super(ic2,i_ds) ) then + if ( diags_super_m % super(ic2,i_ds) > zero ) then + ! Add contribution from diags_super_a to the combined total mass + diags_super_m % super(ic2,i_ds) = diags_super_m % super(ic2,i_ds) & + + diags_super_a % super(ic,i_ds) + ! Calculate mass-fraction weight to be applied to subsequent diagnostics + weight_a(ic) = diags_super_a % super(ic,i_ds) & + / max( diags_super_m % super(ic2,i_ds), min_float ) + else + ! Existing value in "_m" is zero; force weight to be exactly 1.0 + ! (needed for reproducibility on different decompositions). + diags_super_m % super(ic2,i_ds) = diags_super_a % super(ic,i_ds) weight_a(ic) = one end if end do diff --git a/science/physics_schemes/source/convection/comorph/control/fields_type_mod.F90 b/science/physics_schemes/source/convection/comorph/control/fields_type_mod.F90 index 913c7c3fdc..562e245d8e 100644 --- a/science/physics_schemes/source/convection/comorph/control/fields_type_mod.F90 +++ b/science/physics_schemes/source/convection/comorph/control/fields_type_mod.F90 @@ -150,6 +150,13 @@ module fields_type_mod ! (e.g. mixing-ratios are not allowed to go negative) logical, allocatable :: field_positive(:) +! Min and max plausible values for each field, used in run-time checks +real(kind=real_cvprec), allocatable :: field_min(:) +real(kind=real_cvprec), allocatable :: field_max(:) +! Min and max values when in conserved-variable form +real(kind=real_cvprec), allocatable :: field_min_cons(:) +real(kind=real_cvprec), allocatable :: field_max_cons(:) + contains @@ -165,7 +172,8 @@ subroutine fields_set_addresses() cond_params, k_bot_conv, k_top_conv, & i_cond_cl, i_cond_rain, & i_cond_cf, i_cond_snow, i_cond_graup, & - tracer_positive + tracer_positive, & + max_float, zero, one, cp_dry implicit none @@ -243,6 +251,10 @@ subroutine fields_set_addresses() ! Allocate lists allocate( field_names(n_fields+n_tracers) ) allocate( field_positive(n_fields+n_tracers) ) +allocate( field_min(n_fields+n_tracers) ) +allocate( field_max(n_fields+n_tracers) ) +allocate( field_min_cons(n_fields+n_tracers) ) +allocate( field_max_cons(n_fields+n_tracers) ) ! Set the field names... field_names(i_wind_u) = "wind_u" @@ -289,6 +301,62 @@ subroutine fields_set_addresses() end if end if +! Set min and max plausible values for each field +! Winds +do i_field = i_wind_u, i_wind_w + field_min(i_field) = -100.0_real_cvprec + field_max(i_field) = 100.0_real_cvprec +end do +! Temperature: use bounds of the qsat look-up table +field_min(i_temperature) = zero ! 183.15_real_cvprec +field_max(i_temperature) = 338.15_real_cvprec +! q_vap: max 50 g kg-1? +field_min(i_q_vap) = zero +field_max(i_q_vap) = 0.05_real_cvprec +! Condensate: max 50 g kg-1? +do i_field = i_qc_first, i_qc_last + field_min(i_field) = zero + field_max(i_field) = 0.05_real_cvprec +end do +! Cloud fractions: Should be between 0 and 1. +if ( l_cv_cloudfrac ) then + do i_field = i_cf_first, i_cf_last + field_min(i_field) = zero + field_max(i_field) = one + end do +end if +! Tracers: don't know so just set to +/- max possible float +if ( n_tracers > 0 ) then + do i_field = 1, n_tracers + field_min(i_tracers(i_field)) = -max_float + field_max(i_tracers(i_field)) = max_float + end do + if ( allocated(tracer_positive) ) then + ! If host-model has specified tracers positive-only, reset min value + ! to zero for those tracers + do i_field = 1, n_tracers + if ( tracer_positive(i_field) ) field_min(i_tracers(i_field)) = zero + end do + end if +end if + +! Max and min limits for fields in conserved-variable form: +! Initialise the same as above +do i_field = 1, n_fields + n_tracers + field_min_cons(i_field) = field_min(i_field) + field_max_cons(i_field) = field_max(i_field) +end do +! Winds are scaled by 1+qt but that doesn't change it much. +! Temperature is scaled by cp: +field_min_cons(i_temperature) = field_min(i_temperature) * cp_dry +field_max_cons(i_temperature) = field_max(i_temperature) * cp_dry +! Cloud fractions are scaled by Tv: +if ( l_cv_cloudfrac ) then + do i_field = i_cf_first, i_cf_last + field_max_cons(i_field) = field_max(i_field) * field_max(i_temperature) + end do +end if + return end subroutine fields_set_addresses diff --git a/science/physics_schemes/source/convection/comorph/control/grid_type_mod.F90 b/science/physics_schemes/source/convection/comorph/control/grid_type_mod.F90 index 78caccba5c..313d0b74a4 100644 --- a/science/physics_schemes/source/convection/comorph/control/grid_type_mod.F90 +++ b/science/physics_schemes/source/convection/comorph/control/grid_type_mod.F90 @@ -157,7 +157,7 @@ end subroutine grid_compress !---------------------------------------------------------------- subroutine grid_check_bad_values( grid, where_string ) -use comorph_constants_mod, only: name_length +use comorph_constants_mod, only: name_length, zero use check_bad_values_mod, only: check_bad_values_3d implicit none @@ -167,42 +167,41 @@ subroutine grid_check_bad_values( grid, where_string ) character(len=name_length) :: field_name integer :: lb(3), ub(3) -logical, parameter :: l_positive = .true. field_name = "height_full" lb = lbound( grid % height_full ) ub = ubound( grid % height_full ) call check_bad_values_3d( lb, ub, grid % height_full, & where_string, field_name, & - l_positive ) + field_min=zero ) field_name = "height_half" lb = lbound( grid % height_half ) ub = ubound( grid % height_half ) call check_bad_values_3d( lb, ub, grid % height_half, & where_string, field_name, & - l_positive, l_half=.true. ) + field_min=zero, l_half=.true. ) field_name = "pressure_full" lb = lbound( grid % pressure_full ) ub = ubound( grid % pressure_full ) call check_bad_values_3d( lb, ub, grid % pressure_full, & where_string, field_name, & - l_positive ) + field_min=zero ) field_name = "pressure_half" lb = lbound( grid % pressure_half ) ub = ubound( grid % pressure_half ) call check_bad_values_3d( lb, ub, grid % pressure_half, & where_string, field_name, & - l_positive, l_half=.true. ) + field_min=zero, l_half=.true. ) field_name = "rho_dry" lb = lbound( grid % rho_dry ) ub = ubound( grid % rho_dry ) call check_bad_values_3d( lb, ub, grid % rho_dry, & where_string, field_name, & - l_positive ) + field_min=zero ) ! Note: no check on grid % r_surf yet; need to make a 2-D ! version of check_bad_values. diff --git a/science/physics_schemes/source/convection/comorph/control/parcel_type_mod.F90 b/science/physics_schemes/source/convection/comorph/control/parcel_type_mod.F90 index e5abaaa67c..8c56607a6f 100644 --- a/science/physics_schemes/source/convection/comorph/control/parcel_type_mod.F90 +++ b/science/physics_schemes/source/convection/comorph/control/parcel_type_mod.F90 @@ -10,7 +10,7 @@ module parcel_type_mod use cmpr_type_mod, only: cmpr_type -use comorph_constants_mod, only: real_cvprec +use comorph_constants_mod, only: real_cvprec, name_length implicit none @@ -83,6 +83,13 @@ module parcel_type_mod ! Parcel edge virtual temperature (for constructing assumed PDF) integer :: i_edge_virt_temp = 0 +! Name of each parcel field +character(len=name_length), allocatable :: par_names(:) + +! Min and max plausible values for parcel fields, used in bad-value checking +real(kind=real_cvprec), allocatable :: par_min(:) +real(kind=real_cvprec), allocatable :: par_max(:) + contains @@ -92,7 +99,8 @@ module parcel_type_mod !---------------------------------------------------------------- subroutine parcel_set_addresses() -use comorph_constants_mod, only: l_par_core +use comorph_constants_mod, only: real_cvprec, zero, l_par_core +use fields_type_mod, only: field_min, field_max, i_temperature implicit none @@ -111,6 +119,27 @@ subroutine parcel_set_addresses() n_par = n_par + 1 end if +! Set name of each field in the parcel super-array +allocate( par_names(n_par) ) +par_names(i_massflux_d) = "massflux_d" +par_names(i_radius) = "radius" +if ( l_par_core ) par_names(i_edge_virt_temp) = "edge_virt_temp" + +! Set min and max plausible values of parcel fields, used in run-time checking +allocate( par_min(n_par) ) +allocate( par_max(n_par) ) +! Mass-flux +par_min(i_massflux_d) = zero +par_max(i_massflux_d) = 100.0_real_cvprec ! 100 kg m-2 s-1 +! Updraft radius +par_min(i_radius) = zero +par_max(i_radius) = 5.0E4_real_cvprec ! 50km +! Edge virtual temperature +if ( l_par_core ) then + par_min(i_edge_virt_temp) = field_min(i_temperature) + par_max(i_edge_virt_temp) = field_max(i_temperature) +end if + return end subroutine parcel_set_addresses @@ -610,9 +639,8 @@ subroutine parcel_combine( l_tracer, l_down, index_ic, & parcel_a, parcel_m ) use comorph_constants_mod, only: real_cvprec, zero, one, n_tracers, l_par_core -use fields_type_mod, only: n_fields, i_temperature, i_q_vap, & - i_qc_first, i_qc_last -use calc_virt_temp_mod, only: calc_virt_temp +use fields_type_mod, only: n_fields +use core_combine_mod, only: core_combine implicit none @@ -622,11 +650,11 @@ subroutine parcel_combine( l_tracer, l_down, index_ic, & ! Flag for downdraft versus updraft logical, intent(in) :: l_down -! Input properties of one of the parcels to combine +! "_a": Input properties of the "added" parcel (one of the parcels to combine) type(parcel_type), intent(in) :: parcel_a -! IN: properties of the other parcel to combine -! OUT: combined merged parcel properties +! "_m": IN: properties of the other parcel to combine "_a" into +! OUT: properties of the "merged" parcel after combining with "_a" type(parcel_type), intent(in out) :: parcel_m ! Index list for referencing the parcel_m compression list @@ -640,18 +668,7 @@ subroutine parcel_combine( l_tracer, l_down, index_ic, & real(kind=real_cvprec) :: weight_a( parcel_a % cmpr % n_points ) ! Weight to apply to the existing properties of parcel m real(kind=real_cvprec) :: weight_m( parcel_a % cmpr % n_points ) - -! Weights for computing parcel core properties, if used -real(kind=real_cvprec) :: weight_core_a( parcel_a % cmpr % n_points ) -real(kind=real_cvprec) :: weight_core_m( parcel_a % cmpr % n_points ) - -! Virtual temperature of the cores of parcels a and m -real(kind=real_cvprec) :: core_a_virt_temp & - ( parcel_a % cmpr % n_points ) -real(kind=real_cvprec) :: core_m_virt_temp & - ( parcel_m % cmpr % n_points ) - -! Normalisation for weights +! Sum of mass-fluxes used to compute the above real(kind=real_cvprec) :: norm ! Loop counters @@ -724,88 +741,12 @@ subroutine parcel_combine( l_tracer, l_down, index_ic, & if ( l_par_core ) then ! Set parcel core properties... - - ! Compute the core virtual temperature of the two parcels - ! NOTE: this calculation relies on the fact that the parcel - ! core properties are NOT in conserved variable form at this - ! point, whereas the parcel mean properties are. - call calc_virt_temp( parcel_a % cmpr % n_points, & - size(parcel_a % cmpr % index_i), & - parcel_a % core_super(:,i_temperature), & - parcel_a % core_super(:,i_q_vap), & - parcel_a % core_super(:,i_qc_first:i_qc_last), & - core_a_virt_temp ) - call calc_virt_temp( parcel_m % cmpr % n_points, & - size(parcel_m % cmpr % index_i), & - parcel_m % core_super(:,i_temperature), & - parcel_m % core_super(:,i_q_vap), & - parcel_m % core_super(:,i_qc_first:i_qc_last), & - core_m_virt_temp ) - - ! Choose properties from the parcel with the more buoyant core. - - ! Reset the weights so that the core fields will inherit only - ! the values from the most buoyant of the two. - ! Combine the edge virtual temperatures by choosing the least buoyant edge - ! (i.e. we always try to make the PDF of Tv as wide as possible) - if ( l_down ) then - do ic = 1, parcel_a % cmpr % n_points - ic2 = index_ic(ic) - ! Choose most negatively buoyant core for downdrafts - ! TEMPORARY CODE TO PRESERVE KGO: should really test on mass-fluxes > 0 - ! (this code can use core properties from parcel m with zero mass-flux, - ! which is wrong; fix this soon...) - if ( core_a_virt_temp(ic) <= core_m_virt_temp(ic2) .or. & - ( .not. core_m_virt_temp(ic2) > zero ) ) then - weight_core_a(ic) = one - weight_core_m(ic) = zero - else - weight_core_a(ic) = zero - weight_core_m(ic) = one - end if - ! Choose least negatively buoyant edge for downdrafts - if ( parcel_a%par_super(ic,i_edge_virt_temp) > & - parcel_m%par_super(ic2,i_edge_virt_temp) .or. & - ( .not. parcel_m%par_super(ic2,i_edge_virt_temp) > zero ) ) then - parcel_m % par_super(ic2,i_edge_virt_temp) & - = parcel_a % par_super(ic,i_edge_virt_temp) - end if - end do - else - do ic = 1, parcel_a % cmpr % n_points - ic2 = index_ic(ic) - ! Choose most positively buoyant core for updrafts - ! TEMPORARY CODE TO PRESERVE KGO: should really test on mass-fluxes > 0 - ! (this code can use core properties from parcel m with zero mass-flux, - ! which is wrong; fix this soon...) - if ( core_a_virt_temp(ic) >= core_m_virt_temp(ic2) .or. & - ( .not. core_m_virt_temp(ic2) > zero ) ) then - weight_core_a(ic) = one - weight_core_m(ic) = zero - else - weight_core_a(ic) = zero - weight_core_m(ic) = one - end if - ! Choose least positively buoyant edge for downdrafts - if ( parcel_a%par_super(ic,i_edge_virt_temp) < & - parcel_m%par_super(ic2,i_edge_virt_temp) .or. & - ( .not. parcel_m%par_super(ic2,i_edge_virt_temp) > zero ) ) then - parcel_m % par_super(ic2,i_edge_virt_temp) & - = parcel_a % par_super(ic,i_edge_virt_temp) - end if - end do - end if - - ! Compute combined parcel core properties using the weights set above - do i_field = 1, n_fields_tot - do ic = 1, parcel_a % cmpr % n_points - ic2 = index_ic(ic) - parcel_m % core_super(ic2,i_field) & - = weight_core_m(ic) * parcel_m % core_super(ic2,i_field) & - + weight_core_a(ic) * parcel_a % core_super(ic,i_field) - end do - end do - + call core_combine( parcel_a%cmpr%n_points, parcel_m%cmpr%n_points, index_ic, & + size(parcel_a%cmpr%index_i), size(parcel_m%cmpr%index_i), & + 1, n_fields_tot, l_down, & + parcel_a%core_super, parcel_m%core_super, & + parcel_a%par_super(:,i_edge_virt_temp), & + parcel_m%par_super(:,i_edge_virt_temp) ) end if ! ( l_par_core ) @@ -820,7 +761,7 @@ subroutine parcel_check_bad_values( parcel, n_fields_tot, k, & where_string ) use comorph_constants_mod, only: name_length, l_par_core -use fields_type_mod, only: field_names, field_positive +use fields_type_mod, only: field_names, field_min, field_max use check_bad_values_mod, only: check_bad_values_cmpr implicit none @@ -840,27 +781,19 @@ subroutine parcel_check_bad_values( parcel, n_fields_tot, k, & ! Name of individual field character(len=name_length) :: field_name -! Flag for whether field is positive-only -logical :: l_positive ! Loop counter integer :: i_field -! Check mass-flux, parcel radius, turb_len and environment virtual temperature -l_positive = .true. -field_name = "massflux_d" -call check_bad_values_cmpr( parcel % cmpr, k, & - parcel % par_super(:,i_massflux_d), & - where_string, field_name, l_positive) -field_name = "radius" -call check_bad_values_cmpr( parcel % cmpr, k, & - parcel % par_super(:,i_radius), & - where_string, field_name, l_positive) -field_name = "edge_virt_temp" -call check_bad_values_cmpr( parcel % cmpr, k, & - parcel % par_super(:,i_edge_virt_temp), & - where_string, field_name, l_positive) +! Check parcel super-array fields +do i_field = 1, n_par + call check_bad_values_cmpr( parcel % cmpr, k, & + parcel % par_super(:,i_field), & + where_string, par_names(i_field), & + field_min=par_min(i_field), & + field_max=par_max(i_field) ) +end do ! Check mean primary fields do i_field = 1, n_fields_tot @@ -868,7 +801,8 @@ subroutine parcel_check_bad_values( parcel, n_fields_tot, k, & call check_bad_values_cmpr( parcel % cmpr, k, & parcel % mean_super(:,i_field), & where_string, field_name, & - field_positive(i_field) ) + field_min=field_min(i_field), & + field_max=field_max(i_field) ) end do ! Check parcel core fields if used @@ -878,7 +812,8 @@ subroutine parcel_check_bad_values( parcel, n_fields_tot, k, & call check_bad_values_cmpr( parcel % cmpr, k, & parcel % core_super(:,i_field), & where_string, field_name, & - field_positive(i_field) ) + field_min=field_min(i_field), & + field_max=field_max(i_field) ) end do end if diff --git a/science/physics_schemes/source/convection/comorph/control/res_source_type_mod.F90 b/science/physics_schemes/source/convection/comorph/control/res_source_type_mod.F90 index 62f293638b..f735342092 100644 --- a/science/physics_schemes/source/convection/comorph/control/res_source_type_mod.F90 +++ b/science/physics_schemes/source/convection/comorph/control/res_source_type_mod.F90 @@ -516,9 +516,9 @@ end subroutine res_source_combine subroutine res_source_check_bad_values( res_source, n_fields_tot, & k, where_string ) -use comorph_constants_mod, only: name_length +use comorph_constants_mod, only: zero, name_length use fields_type_mod, only: field_names -use cloudfracs_type_mod, only: n_convcloud, convcloud_names +use cloudfracs_type_mod, only: n_convcloud, convcloud_names, convcloud_max use check_bad_values_mod, only: check_bad_values_cmpr implicit none @@ -536,44 +536,34 @@ subroutine res_source_check_bad_values( res_source, n_fields_tot, & ! we are, for constructing error message if bad value found. character(len=name_length), intent(in) :: where_string -! Flag for whether field is positive-only -logical :: l_positive - ! Loop counter integer :: i_field ! Check fields in the res_super array -l_positive = .true. do i_field = 1, n_res call check_bad_values_cmpr( res_source % cmpr, k, & res_source % res_super(:,i_field), & where_string, res_source_names(i_field), & - l_positive) + field_min=zero ) end do ! Check source terms for primary fields -l_positive = .false. do i_field = 1, n_fields_tot call check_bad_values_cmpr( res_source % cmpr, k, & res_source % fields_super(:,i_field), & - where_string, field_names(i_field), & - l_positive ) + where_string, field_names(i_field) ) end do if ( n_convcloud > 0 ) then ! Check convective cloud fields... - - ! These must all be positive - l_positive = .true. - do i_field = 1, n_convcloud call check_bad_values_cmpr( res_source % cmpr, k, & res_source % convcloud_super(:,i_field), & where_string, convcloud_names(i_field), & - l_positive ) + field_min=zero, & + field_max=convcloud_max(i_field) ) end do - end if diff --git a/science/physics_schemes/source/convection/comorph/control/set_cloudfracs_k.F90 b/science/physics_schemes/source/convection/comorph/control/set_cloudfracs_k.F90 index b77c86eb6b..2f3ceb4ac5 100644 --- a/science/physics_schemes/source/convection/comorph/control/set_cloudfracs_k.F90 +++ b/science/physics_schemes/source/convection/comorph/control/set_cloudfracs_k.F90 @@ -13,7 +13,7 @@ module set_cloudfracs_k_mod contains -! Subroutine to set the super-array containing environment cloud +! Subroutine to set the super-array containing environment cloud and precip ! fractions at the current full model-level. Either needs to copy them ! from the primary fields array (if they are treated as prognostics), ! or compress from separate diagnostic cloud arrays otherwise. diff --git a/science/physics_schemes/source/convection/comorph/control/set_dependent_constants.F90 b/science/physics_schemes/source/convection/comorph/control/set_dependent_constants.F90 index ffd21f0b03..035df91da1 100644 --- a/science/physics_schemes/source/convection/comorph/control/set_dependent_constants.F90 +++ b/science/physics_schemes/source/convection/comorph/control/set_dependent_constants.F90 @@ -23,22 +23,20 @@ module set_dependent_constants_mod subroutine set_dependent_constants() use comorph_constants_mod, only: newline, l_init_constants, & - L_con_ref, L_sub_ref, L_fus_ref, & - L_con_0, L_sub_0, L_fus_0, & - melt_temp, & - l_cv_rain, l_cv_cf, l_cv_snow, l_cv_graup, & - n_cond_species, & - n_cond_species_liq, n_cond_species_ice, & - i_cond_cl, i_cond_rain, & - i_cond_cf, i_cond_snow, i_cond_graup, & - k_bot_conv, k_top_conv, & - params_cl, params_rain, & - params_cf, params_snow, params_graup, & - cond_params, & - i_cl, i_rain, i_cf, i_snow, i_graup, & - cp_vap, cp_liq, cp_ice, & - rho_liq, rho_ice, rho_rim, & - nx_full, ny_full, k_bot_conv, k_top_conv, k_top_init + L_con_ref, L_sub_ref, L_fus_ref, & + L_con_0, L_sub_0, L_fus_0, melt_temp, & + l_cv_rain, l_cv_cf, l_cv_snow, l_cv_graup, & + n_cond_species, & + n_cond_species_liq, n_cond_species_ice, & + i_cond_cl, i_cond_rain, & + i_cond_cf, i_cond_snow, i_cond_graup, & + cond_params, params_cl, params_rain, & + params_cf, params_snow, params_graup, & + i_cl, i_rain, i_cf, i_snow, i_graup, & + cp_vap, cp_liq, cp_ice, & + rho_liq, rho_ice, rho_rim, & + nx_full, ny_full, k_bot_conv, k_top_conv, & + k_top_init use raise_error_mod, only: raise_fatal diff --git a/science/physics_schemes/source/convection/comorph/control/turb_type_mod.F90 b/science/physics_schemes/source/convection/comorph/control/turb_type_mod.F90 index 7de762a9cd..e7b4570b4e 100644 --- a/science/physics_schemes/source/convection/comorph/control/turb_type_mod.F90 +++ b/science/physics_schemes/source/convection/comorph/control/turb_type_mod.F90 @@ -82,12 +82,6 @@ module turb_type_mod character(len=name_length), parameter :: turb_names(n_turb) & = [ "w_var ", "f_templ ", "f_q_tot ", "f_wind_u", "f_wind_v" ] -! Flag for whether each turbulence field should only be positive -logical, parameter :: turb_positive(n_turb) & - = [ .true., .false., .false., .false., .false. ] -! TKE and diffusivities must be positive, -! but fluxes could go either way. - ! Max allowed ratio of heat-flux over sqrt(TKE) / K real(kind=real_hmprec), parameter :: max_templ = 20.0 ! Max allowed ratio of moisture-flux over sqrt(TKE) / kg kg-1 @@ -276,7 +270,7 @@ end subroutine turb_list_clear !---------------------------------------------------------------- subroutine turb_check_bad_values( turb, where_string ) -use comorph_constants_mod, only: name_length +use comorph_constants_mod, only: name_length, zero use check_bad_values_mod, only: check_bad_values_3d implicit none @@ -293,9 +287,6 @@ subroutine turb_check_bad_values( turb, where_string ) ! Lower and upper bounds of array integer :: lb(3), ub(3) -! Flag passed into check_bad_values; -logical, parameter :: l_positive_true = .true. - ! Loop counter integer :: i_turb @@ -305,7 +296,6 @@ subroutine turb_check_bad_values( turb, where_string ) ub = ubound( turb % list(i_turb)%pt ) call check_bad_values_3d( lb, ub, turb%list(i_turb)%pt, & where_string, turb_names(i_turb), & - turb_positive(i_turb), & l_half=.true., l_init=.true. ) end do @@ -316,7 +306,7 @@ subroutine turb_check_bad_values( turb, where_string ) ub = ubound( turb % lengthscale ) call check_bad_values_3d( lb, ub, turb % lengthscale, & where_string, field_name, & - l_positive_true, l_init=.true. ) + field_min=zero, l_init=.true. ) return diff --git a/science/physics_schemes/source/convection/comorph/interface/um/comorph_um_namelist_mod.F90 b/science/physics_schemes/source/convection/comorph/interface/um/comorph_um_namelist_mod.F90 index 45ec7e44e9..a98da61cd4 100644 --- a/science/physics_schemes/source/convection/comorph/interface/um/comorph_um_namelist_mod.F90 +++ b/science/physics_schemes/source/convection/comorph/interface/um/comorph_um_namelist_mod.F90 @@ -56,13 +56,13 @@ module comorph_um_namelist_mod integer, parameter :: w_dependence = 3 ! add further dependence on max w integer, parameter :: linear_qfacrain_dep = 4 ! scale precip dependence by q +! Comorph internal switches +! (allowed values are stored in comorph_constants_mod). + integer :: par_radius_evol_method = imdi ! Switch for how parcel radius evolves ! with height in the plume -! (allowed values are stored in comorph_constants_mod). integer :: autoc_opt = imdi ! Switch for autoconversion option -! (allowed values are stored in comorph_constants_mod). - integer :: n_dndraft_types = imdi ! Number of independent downdraft ! types @@ -134,6 +134,11 @@ module comorph_um_namelist_mod real(kind=real_umphys) :: overlap_power = rmdi ! Note: do not set to exactly zero as this causes a singularity! +! Max and min limits on the core-mean-ratio (determines assumed in-plume +! PDF shape used for detraiment) +real(kind=real_umphys) :: min_cmr = rmdi +real(kind=real_umphys) :: max_cmr = rmdi + ! Plume microphysics parameters @@ -146,17 +151,25 @@ module comorph_um_namelist_mod ! Density of rimed ice (used for graupel) real(kind=real_umphys) :: rho_rim = rmdi -! Reciprocal of fac_tdep_n +! Temperature-dependent liquid-cloud number concentration slope +! Liquid-cloud number n(T) = n0 exp( ( T - Tmelt ) / tdep_n_cl ) +! ( but limited above Tmelt and below T_homnuc) +real(kind=real_umphys) :: tdep_n_cl = rmdi + ! Temperature-dependent ice number concentration slope -! The number concentration n(T) will be given by: -! n(T) = n0 exp( fac_tdep_n ( T - Tmelt ) ) +! Ice number n(T) = n0 exp( -( T - Tmelt ) / tdep_n_cf ) ! ( but limited above Tmelt and below T_homnuc) -real(kind=real_umphys) :: r_fac_tdep_n = rmdi +real(kind=real_umphys) :: tdep_n_cf = rmdi ! Heterogeneous nucleation temeprature / K ! Gradual freezing starts below this real(kind=real_umphys) :: hetnuc_temp = rmdi +! Ice crystal non-spherical area factor. This scales up the assumed surface +! area of ice crystals relative to what it would be if they were spheres +! (affects vapour deposition, fall-speed, riming...) +real(kind=real_umphys) :: cf_area_coef = rmdi + ! Asymptotic drag coefficient for a sphere at high Reynolds ! number limit real(kind=real_umphys) :: drag_coef_cond = rmdi @@ -169,6 +182,13 @@ module comorph_um_namelist_mod ! deflection flow around hydrometeors real(kind=real_umphys) :: col_eff_coef = rmdi +! Number concentrations for liquid-cloud, ice-cloud, rain, snow, graupel +real(kind=real_umphys) :: nconc_cl = rmdi +real(kind=real_umphys) :: nconc_cf = rmdi +real(kind=real_umphys) :: nconc_rain = rmdi +real(kind=real_umphys) :: nconc_snow = rmdi +real(kind=real_umphys) :: nconc_graup = rmdi + !------------------------------------------------------------------------------ ! Define namelist &Run_Comorph read in from CNTLATM control file. !------------------------------------------------------------------------------ @@ -184,10 +204,12 @@ module comorph_um_namelist_mod core_ent_fac, rain_area_min, cf_conv_fac, drag_coef_par, par_gen_rhpert, & par_gen_mass_fac, wind_w_fac, wind_w_buoy_fac, par_gen_pert_fac, & ass_min_radius, par_gen_core_fac, overlap_power, ent_coef, & +min_cmr, max_cmr, & ! Plume microphysics parameters -rho_rim, r_fac_tdep_n, hetnuc_temp, drag_coef_cond, & +rho_rim, tdep_n_cl, tdep_n_cf, hetnuc_temp, cf_area_coef, drag_coef_cond, & vent_factor, col_eff_coef, q_cl_auto, coef_auto, & +nconc_cl, nconc_cf, nconc_rain, nconc_snow, nconc_graup, & ! Logical switches l_core_ent_cmr, l_resdep_precipramp @@ -243,6 +265,8 @@ subroutine check_run_comorph() ! Checking reals within allowed range - ranges as in meta-data used for GUI +if (l_resdep_precipramp) call chk_var(dx_ref,'dx_ref','[100.0:1000000.0]') + call chk_var(par_gen_mass_fac,'par_gen_mass_fac','[0.01:1.0]') call chk_var(drag_coef_cond,'drag_coef_cond','[0.2:1.0]') @@ -253,6 +277,8 @@ subroutine check_run_comorph() call chk_var(hetnuc_temp,'hetnuc_temp','[230.0:273.0]') +call chk_var(cf_area_coef,'cf_area_coef','[1.0:1000.0]') + call chk_var(wind_w_fac,'wind_w_fac','[0.1:10.0]') call chk_var(wind_w_buoy_fac,'wind_w_buoy_fac','[0.5:2.0]') @@ -265,7 +291,8 @@ subroutine check_run_comorph() call chk_var(ent_coef,'ent_coef','[0.1:0.4]') -if (l_resdep_precipramp) call chk_var(dx_ref,'dx_ref','[100.0:1000000.0]') +call chk_var(min_cmr,'min_cmr','[1.0:3.0]') +call chk_var(max_cmr,'max_cmr','[3.0:10.0]') !--------------------------------------------------------------------------- if (lhook) call dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) @@ -324,9 +351,15 @@ subroutine print_nlist_run_comorph() call umPrint(lineBuffer,src=ModuleName) write(lineBuffer,"(A,ES14.6)")' ent_coef = ',ent_coef call umPrint(lineBuffer,src=ModuleName) +write(lineBuffer,"(A,ES14.6)")' min_cmr = ',min_cmr +call umPrint(lineBuffer,src=ModuleName) +write(lineBuffer,"(A,ES14.6)")' max_cmr = ',max_cmr +call umPrint(lineBuffer,src=ModuleName) write(lineBuffer,"(A,ES14.6)")' rho_rim = ',rho_rim call umPrint(lineBuffer,src=ModuleName) -write(lineBuffer,"(A,ES14.6)")' r_fac_tdep_n = ',r_fac_tdep_n +write(lineBuffer,"(A,ES14.6)")' tdep_n_cl = ',tdep_n_cl +call umPrint(lineBuffer,src=ModuleName) +write(lineBuffer,"(A,ES14.6)")' tdep_n_cf = ',tdep_n_cf call umPrint(lineBuffer,src=ModuleName) write(lineBuffer,"(A,ES14.6)")' par_radius_knob = ',par_radius_knob @@ -348,9 +381,10 @@ subroutine print_nlist_run_comorph() write(lineBuffer,"(A,ES14.6)")' par_gen_rhpert = ',par_gen_rhpert call umPrint(lineBuffer,src=ModuleName) - write(lineBuffer,"(A,ES14.6)")' hetnuc_temp = ',hetnuc_temp call umPrint(lineBuffer,src=ModuleName) +write(lineBuffer,"(A,ES14.6)")' cf_area_coef = ',cf_area_coef +call umPrint(lineBuffer,src=ModuleName) write(lineBuffer,"(A,ES14.6)")' drag_coef_cond = ',drag_coef_cond call umPrint(lineBuffer,src=ModuleName) write(lineBuffer,"(A,ES14.6)")' vent_factor = ',vent_factor @@ -361,7 +395,16 @@ subroutine print_nlist_run_comorph() call umPrint(lineBuffer,src=ModuleName) write(lineBuffer,"(A,ES14.6)")' coef_auto = ',coef_auto call umPrint(lineBuffer,src=ModuleName) - +write(lineBuffer,"(A,ES14.6)")' nconc_cl = ',nconc_cl +call umPrint(lineBuffer,src=ModuleName) +write(lineBuffer,"(A,ES14.6)")' nconc_cf = ',nconc_cf +call umPrint(lineBuffer,src=ModuleName) +write(lineBuffer,"(A,ES14.6)")' nconc_rain = ',nconc_rain +call umPrint(lineBuffer,src=ModuleName) +write(lineBuffer,"(A,ES14.6)")' nconc_snow = ',nconc_snow +call umPrint(lineBuffer,src=ModuleName) +write(lineBuffer,"(A,ES14.6)")' nconc_graup = ',nconc_graup +call umPrint(lineBuffer,src=ModuleName) ! Logicals @@ -404,7 +447,7 @@ subroutine read_nml_run_comorph(unit_in) ! set number of each type of variable in my_namelist type integer, parameter :: no_of_types = 3 integer, parameter :: n_int = 4 -integer, parameter :: n_real = 25 +integer, parameter :: n_real = 34 integer, parameter :: n_log = 2 type :: my_namelist @@ -430,14 +473,23 @@ subroutine read_nml_run_comorph(unit_in) real(kind=real_umphys) :: par_gen_core_fac real(kind=real_umphys) :: overlap_power real(kind=real_umphys) :: ent_coef + real(kind=real_umphys) :: min_cmr + real(kind=real_umphys) :: max_cmr real(kind=real_umphys) :: rho_rim - real(kind=real_umphys) :: r_fac_tdep_n + real(kind=real_umphys) :: tdep_n_cl + real(kind=real_umphys) :: tdep_n_cf real(kind=real_umphys) :: hetnuc_temp + real(kind=real_umphys) :: cf_area_coef real(kind=real_umphys) :: drag_coef_cond real(kind=real_umphys) :: vent_factor real(kind=real_umphys) :: col_eff_coef real(kind=real_umphys) :: q_cl_auto real(kind=real_umphys) :: coef_auto + real(kind=real_umphys) :: nconc_cl + real(kind=real_umphys) :: nconc_cf + real(kind=real_umphys) :: nconc_rain + real(kind=real_umphys) :: nconc_snow + real(kind=real_umphys) :: nconc_graup logical :: l_core_ent_cmr logical :: l_resdep_precipramp end type my_namelist @@ -482,14 +534,23 @@ subroutine read_nml_run_comorph(unit_in) my_nml % par_gen_core_fac = par_gen_core_fac my_nml % overlap_power = overlap_power my_nml % ent_coef = ent_coef + my_nml % min_cmr = min_cmr + my_nml % max_cmr = max_cmr my_nml % rho_rim = rho_rim - my_nml % r_fac_tdep_n = r_fac_tdep_n + my_nml % tdep_n_cl = tdep_n_cl + my_nml % tdep_n_cf = tdep_n_cf my_nml % hetnuc_temp = hetnuc_temp + my_nml % cf_area_coef = cf_area_coef my_nml % drag_coef_cond = drag_coef_cond my_nml % vent_factor = vent_factor my_nml % col_eff_coef = col_eff_coef my_nml % q_cl_auto = q_cl_auto my_nml % coef_auto = coef_auto + my_nml % nconc_cl = nconc_cl + my_nml % nconc_cf = nconc_cf + my_nml % nconc_rain = nconc_rain + my_nml % nconc_snow = nconc_snow + my_nml % nconc_graup = nconc_graup ! end of reals ! logicals my_nml % l_core_ent_cmr = l_core_ent_cmr @@ -523,14 +584,23 @@ subroutine read_nml_run_comorph(unit_in) par_gen_core_fac = my_nml % par_gen_core_fac overlap_power = my_nml % overlap_power ent_coef = my_nml % ent_coef + min_cmr = my_nml % min_cmr + max_cmr = my_nml % max_cmr rho_rim = my_nml % rho_rim - r_fac_tdep_n = my_nml % r_fac_tdep_n + tdep_n_cl = my_nml % tdep_n_cl + tdep_n_cf = my_nml % tdep_n_cf hetnuc_temp = my_nml % hetnuc_temp + cf_area_coef = my_nml % cf_area_coef drag_coef_cond = my_nml % drag_coef_cond vent_factor = my_nml % vent_factor col_eff_coef = my_nml % col_eff_coef q_cl_auto = my_nml % q_cl_auto coef_auto = my_nml % coef_auto + nconc_cl = my_nml % nconc_cl + nconc_cf = my_nml % nconc_cf + nconc_rain = my_nml % nconc_rain + nconc_snow = my_nml % nconc_snow + nconc_graup = my_nml % nconc_graup ! end of reals l_core_ent_cmr = my_nml % l_core_ent_cmr l_resdep_precipramp = my_nml % l_resdep_precipramp diff --git a/science/physics_schemes/source/convection/comorph/interface/um/set_constants_from_um.F90 b/science/physics_schemes/source/convection/comorph/interface/um/set_constants_from_um.F90 index cac56e1c38..ec0ca9effe 100644 --- a/science/physics_schemes/source/convection/comorph/interface/um/set_constants_from_um.F90 +++ b/science/physics_schemes/source/convection/comorph/interface/um/set_constants_from_um.F90 @@ -57,12 +57,16 @@ subroutine set_constants_from_um( n_conv_levels, ntra_fld, i_tr_vars ) par_gen_core_fac_um => par_gen_core_fac, & ent_coef_um => ent_coef, & overlap_power_um => overlap_power, & + min_cmr_um => min_cmr, & + max_cmr_um => max_cmr, & rho_rim_um => rho_rim, & hetnuc_temp_um => hetnuc_temp, & + cf_area_coef_um => cf_area_coef, & drag_coef_cond_um => drag_coef_cond, & vent_factor_um => vent_factor, & col_eff_coef_um => col_eff_coef, & - r_fac_tdep_n + tdep_n_cl, tdep_n_cf, & + nconc_cl, nconc_cf, nconc_rain, nconc_snow, nconc_graup ! comorph settings and constants set by this routine use comorph_constants_mod, only: real_cvprec, nx_full, ny_full, & @@ -87,7 +91,8 @@ subroutine set_constants_from_um( n_conv_levels, ntra_fld, i_tr_vars ) drag_coef_par, rho_rim, & par_gen_mass_fac, wind_w_fac, wind_w_buoy_fac,& ass_min_radius, par_gen_core_fac, ent_coef, & - overlap_power, fac_tdep_n, hetnuc_temp, & + overlap_power, min_cmr, max_cmr, & + hetnuc_temp, & drag_coef_cond, vent_factor, col_eff_coef implicit none @@ -294,18 +299,34 @@ subroutine set_constants_from_um( n_conv_levels, ntra_fld, i_tr_vars ) ! inside the parcel overlap_power = real(overlap_power_um, real_cvprec ) +! Max and min limits on core-mean-ratio +min_cmr = real(min_cmr_um, real_cvprec ) +max_cmr = real(max_cmr_um, real_cvprec ) + ! Density of rimed ice (used for graupel) rho_rim = real(rho_rim_um, real_cvprec ) +! Temperature-dependent liquid-cloud number concentration slope +if ( abs(tdep_n_cl) > 0.0 ) then + params_cl % fac_tdep_n = one/real( tdep_n_cl, real_cvprec ) +else ! Disable T-dependence if T-scale is zero + params_cl % fac_tdep_n = zero +end if + ! Temperature-dependent ice number concentration slope -! The number concentration n(T) will be given by: -! n(T) = n0 exp( fac_tdep_n ( T - Tmelt ) ) -! ( but limited above Tmelt and below T_homnuc) -fac_tdep_n = -one/real( r_fac_tdep_n, real_cvprec ) +if ( abs(tdep_n_cf) > 0.0 ) then + params_cf % fac_tdep_n = -one/real( tdep_n_cf, real_cvprec ) +else ! Disable T-dependence if T-scale is zero + params_cf % fac_tdep_n = zero +end if ! Heterogeneous nucleation temeprature / K hetnuc_temp = real(hetnuc_temp_um, real_cvprec ) +! Ice crystal non-spherical area factor. +params_cf % area_coef = real(cf_area_coef_um, real_cvprec ) +params_snow % area_coef = real(cf_area_coef_um, real_cvprec ) + ! Asymptotic drag coefficient for a sphere at high Reynolds ! number limit drag_coef_cond = real(drag_coef_cond_um, real_cvprec ) @@ -318,6 +339,13 @@ subroutine set_constants_from_um( n_conv_levels, ntra_fld, i_tr_vars ) ! deflection flow around hydrometeors col_eff_coef = real(col_eff_coef_um, real_cvprec ) +! Set prescribed number concentrations from namelist +params_cl % n = real( nconc_cl, real_cvprec ) +params_cf % n = real( nconc_cf, real_cvprec ) +params_rain % n = real( nconc_rain, real_cvprec ) +params_snow % n = real( nconc_snow, real_cvprec ) +params_graup % n = real( nconc_graup, real_cvprec ) + ! Set threshold for using indirect indexing versus straight do-loops ! over all points in various calculations inside comorph. ! A value somewhere between 0 and 1 should optimise performance, diff --git a/science/physics_schemes/source/convection/comorph/microphysics/autoconversion.F90 b/science/physics_schemes/source/convection/comorph/microphysics/autoconversion.F90 new file mode 100644 index 0000000000..d779b7b134 --- /dev/null +++ b/science/physics_schemes/source/convection/comorph/microphysics/autoconversion.F90 @@ -0,0 +1,267 @@ +! *****************************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******************************* + +! Code Owner: Please refer to the UM file CodeOwners.txt +! This file belongs in section: convection_comorph + +module autoconversion_mod + +implicit none + +contains + +! Subroutine to do autoconversion of liquid-cloud to rain or +! ice-cloud to snow in BOTEMS (Back-Of-The-Envelope Microphysical Scheme) +subroutine autoconversion( n_points, nc, index_ic, & + delta_t, vert_len, wf_cond1, & + q_cond1, q_cond2, & + l_diags, diags_cond1, diags_cond2, & + n_points_diag, n_diags, diags_super ) + +use comorph_constants_mod, only: real_cvprec, zero, one, two, four, & + indi_thresh, autoc_opt, autoc_linear, & + autoc_quadratic, q_cl_auto, coef_auto +use moist_proc_diags_type_mod, only: diags_cond_type + +implicit none + +! Number of points +integer, intent(in) :: n_points + +! Number of points where autoconverted species is non-zero +integer, intent(in) :: nc +! Indices of those points +integer, intent(in) :: index_ic(n_points) + +! Time interval for converting process rates to increments. +real(kind=real_cvprec), intent(in) :: delta_t(n_points) + +! Vertical length-scale of the parcel. +real(kind=real_cvprec), intent(in) :: vert_len(n_points) + +! Fall-speed of cloud-water species +real(kind=real_cvprec), intent(in) :: wf_cond1(n_points) + +! Mixing ratio of cloud-water species +real(kind=real_cvprec), intent(in out) :: q_cond1(n_points) +! Mixing ratio of precipitation-water species which q_cond1 autoconverts to +real(kind=real_cvprec), intent(in out) :: q_cond2(n_points) + +! Master switch for whether or not to calculate any diagnostics +logical, intent(in) :: l_diags +! Structures containing diagnostic flags etc. +type(diags_cond_type), intent(in) :: diags_cond1 +type(diags_cond_type), intent(in) :: diags_cond2 +! Number of points in the diagnostics super-array +integer, intent(in) :: n_points_diag +! Total number of diagnostics in the super-array +integer, intent(in) :: n_diags +! Super-array to store all the diagnostics +real(kind=real_cvprec), intent(in out) :: diags_super & + ( n_points_diag, n_diags ) + + +! Increment to q_rain due to autoconversion +real(kind=real_cvprec) :: dq_auto(n_points) + +! Temporary store for ratio fall-speed / height-interval +real(kind=real_cvprec) :: wf_over_lz + +! Coefficients of quadratic equation +real(kind=real_cvprec) :: a_quad, b_quad, c_quad + +! Solution to quadratic equation for final q_cl +real(kind=real_cvprec) :: q_cl_out + +! Loop counters +integer :: ic, ic2, i_super + + +! If calculation required at majority of points +if ( real(nc,real_cvprec) > indi_thresh * real(n_points,real_cvprec) ) then + ! Full-field calculation + + ! Implicit discretisation of autoconversion and fall-out: + ! Two choices, linear or quadratic + if ( autoc_opt==autoc_linear ) then + ! ---------------------------------------------------------------- + ! Linear, with threshold + ! + ! d_q_cl_auto = delta_t coef_auto ( q_cl_out - q_cl_auto ) + ! d_q_cl_fall = delta_t q_cl_out wf/lz + ! + ! Therefore + ! + ! q_cl_out = q_cl_in - delta_t ( + ! coef_auto ( q_cl_out - q_cl_auto ) + ! + q_cl_out wf/lz ) + ! + ! Rearranging: + ! + ! q_cl_out - q_cl_auto = q_cl_in - q_cl_auto - delta_t ( + ! coef_auto ( q_cl_out - q_cl_auto ) + ! + ( q_cl_out - q_cl_auto + q_cl_auto ) wf/lz ) + ! + ! => ( q_cl_out - q_cl_auto ) ( 1 + delta_t ( coef_auto + wf/lz ) ) + ! = ( q_cl_in - q_cl_auto ) - delta_t wf/lz q_cl_auto + ! + ! => ( q_cl_out - q_cl_auto ) + ! = ( ( q_cl_in - q_cl_auto ) - delta_t wf/lz q_cl_auto ) + ! / ( 1 + delta_t ( coef_auto + wf/lz ) ) + ! + ! Then multiply this by delta_t * coef_auto to get + ! the autoconversion increment + ! + ! Note: the solution only makes sense where the numerator + ! ( q_cl_in - q_cl_auto ) - delta_t wf/lz q_cl_auto + ! is positive (otherwise, the solution with no autoconversion + ! yields q_cl < q_cl_auto, and we should have no autoconversion + ! as we're below the threshold). + + do ic = 1, n_points + ! Store ratio fall-speed / height-interval + wf_over_lz = wf_cond1(ic) / vert_len(ic) + ! Compute autoconversion increment, using the above formula + dq_auto(ic) = delta_t(ic) * coef_auto & + * max( ( q_cond1(ic) - q_cl_auto ) & + - delta_t(ic) * wf_over_lz * q_cl_auto, zero ) & + / ( one + delta_t(ic) * ( coef_auto + wf_over_lz ) ) + end do + + else if ( autoc_opt==autoc_quadratic ) then + ! ----------------------------------------------------------- + ! Quadratic, no threshold + ! + ! d_q_cl_auto = delta_t * coef_auto * q_cl_out^2 + ! d_q_cl_fall = delta_t * q_cl_out * wf/lz + ! + ! Therefore + ! + ! q_cl_out = q_cl_in - delta_t ( coef_auto q_cl_out^2 + ! + q_cl_out wf/lz ) + ! + ! Rearranging: + ! + ! delta_t*coef_auto * q_cl_out^2 + (1+delta_t*wf/lz) * q_cl_out + ! - q_cl_in = 0 + ! + ! => quadratic in q_cl_out = (-b +/- sqrt(b^2-4ac) )/2a + ! Only +ve root will be +ve + ! => q_cl_out = ( -(1+delta_t*wf/lz) + + ! sqrt( (1+delta_t*wf/lz)^2 + 4*delta_t*coef_auto*q_cl_in ) ) / + ! (2*delta_t*coef_auto) + ! + ! Then autoconversion increment is delta_t coef_auto * q_cl_out^2 + ! + ! Note: the solution only makes sense where the numerator + ! is positive. + + do ic = 1, n_points + ! Store coefficients in the quadratic formula + a_quad = delta_t(ic) * coef_auto + b_quad = one + delta_t(ic) * wf_cond1(ic) / vert_len(ic) + c_quad = - q_cond1(ic) + ! Compute the solution to the quadratic + q_cl_out = max( zero, & + ( -b_quad + sqrt( b_quad*b_quad - four*a_quad*c_quad ) ) & + /( two*a_quad ) ) + ! Find increment consistent with this + dq_auto(ic) = delta_t(ic) * coef_auto * q_cl_out * q_cl_out + end do + + end if ! ( autoc_opt ) + + do ic = 1, n_points + ! Don't allow increment to reduce q_cl below min allowed value + dq_auto(ic) = min( dq_auto(ic), q_cond1(ic) ) + ! Add increment to q_cl and q_rain + q_cond1(ic) = q_cond1(ic) - dq_auto(ic) + q_cond2(ic) = q_cond2(ic) + dq_auto(ic) + end do + + if ( l_diags ) then + ! Save autoconversion increment diagnostics, if requested + if ( diags_cond1 % dq_aut % flag ) then + ! Extract super-array address + i_super = diags_cond1 % dq_aut % i_super + ! Copy increment + do ic = 1, n_points + diags_super(ic,i_super) = diags_super(ic,i_super) - dq_auto(ic) + end do + end if + if ( diags_cond2 % dq_aut % flag ) then + ! Extract super-array address + i_super = diags_cond2 % dq_aut % i_super + ! Copy increment + do ic = 1, n_points + diags_super(ic,i_super) = diags_super(ic,i_super) + dq_auto(ic) + end do + end if + end if ! ( l_diags ) + + ! Liquid-cloud present at only a small fraction of points +else + ! Compressed version of exactly the same calculation + + if ( autoc_opt==autoc_linear ) then + ! linear, with threshold q_cl_auto + + do ic2 = 1, nc + ic = index_ic(ic2) + wf_over_lz = wf_cond1(ic) / vert_len(ic) + dq_auto(ic) = delta_t(ic) * coef_auto & + * max( ( q_cond1(ic) - q_cl_auto ) & + - delta_t(ic) * wf_over_lz * q_cl_auto, zero ) & + / ( one + delta_t(ic) * ( coef_auto + wf_over_lz ) ) + end do + + else if ( autoc_opt==autoc_quadratic ) then + ! quadratic, with no threshold q_cl + + do ic2 = 1, nc + ic = index_ic(ic2) + a_quad = delta_t(ic) * coef_auto + b_quad = one + delta_t(ic) * wf_cond1(ic) / vert_len(ic) + c_quad = - q_cond1(ic) + q_cl_out = max( zero, & + ( -b_quad + sqrt( b_quad*b_quad - four*a_quad*c_quad ) ) & + /( two*a_quad ) ) + dq_auto(ic) = delta_t(ic) * coef_auto * q_cl_out * q_cl_out + end do + + end if ! ( autoc_opt ) + + do ic2 = 1, nc + ic = index_ic(ic2) + dq_auto(ic) = min( dq_auto(ic), q_cond1(ic) ) + q_cond1(ic) = q_cond1(ic) - dq_auto(ic) + q_cond2(ic) = q_cond2(ic) + dq_auto(ic) + end do + + if ( l_diags ) then + if ( diags_cond1 % dq_aut % flag ) then + i_super = diags_cond1 % dq_aut % i_super + do ic2 = 1, nc + ic = index_ic(ic2) + diags_super(ic,i_super) = diags_super(ic,i_super) - dq_auto(ic) + end do + end if + if ( diags_cond2 % dq_aut % flag ) then + i_super = diags_cond2 % dq_aut % i_super + do ic2 = 1, nc + ic = index_ic(ic2) + diags_super(ic,i_super) = diags_super(ic,i_super) + dq_auto(ic) + end do + end if + end if ! ( l_diags ) + +end if ! Test whether to use compressed vs full-field calculation + + +return +end subroutine autoconversion + +end module autoconversion_mod diff --git a/science/physics_schemes/source/convection/comorph/microphysics/calc_cond_properties.F90 b/science/physics_schemes/source/convection/comorph/microphysics/calc_cond_properties.F90 index 9374a0d823..7e449b4b75 100644 --- a/science/physics_schemes/source/convection/comorph/microphysics/calc_cond_properties.F90 +++ b/science/physics_schemes/source/convection/comorph/microphysics/calc_cond_properties.F90 @@ -34,9 +34,8 @@ subroutine calc_cond_properties( n_points, & r_cond, wf_cond, & kq_cond, kt_cond ) -use comorph_constants_mod, only: real_cvprec, cond_params_type, & - melt_temp, homnuc_temp, fac_tdep_n, & - solve_wf_tolerance +use comorph_constants_mod, only: real_cvprec, zero, cond_params_type, & + melt_temp, homnuc_temp, solve_wf_tolerance use set_cond_radius_mod, only: set_cond_radius use fall_speed_mod, only: fall_speed use fall_out_mod, only: fall_out @@ -101,10 +100,10 @@ subroutine calc_cond_properties( n_points, & end do ! If this species uses temperature-dependent number concentration -if ( cond_params % l_tdep_n ) then +if ( abs( cond_params % fac_tdep_n ) > zero ) then do ic = 1, n_points ! Scale n by temperature-dependent factor - n_cond(ic) = n_cond(ic) * exp( fac_tdep_n & + n_cond(ic) = n_cond(ic) * exp( cond_params % fac_tdep_n & * ( max(min( ref_temp(ic), melt_temp ),homnuc_temp) & - melt_temp ) ) end do @@ -117,8 +116,8 @@ subroutine calc_cond_properties( n_points, & ! Set first-guess fall-speed call fall_speed( n_points, & - cond_params % area_coef, cond_params % rho, & - rho_wet, r_cond, wf_cond ) + cond_params % area_coef, cond_params % rho, & + rho_wet, r_cond, wf_cond ) ! Find the fall-speed that you'd get by using the current-guess @@ -208,14 +207,12 @@ subroutine calc_cond_properties_cmpr( n_points, nc, index_ic, & real(kind=real_cvprec), intent(in) :: rho_dry(n_points) real(kind=real_cvprec), intent(in) :: rho_wet(n_points) real(kind=real_cvprec), intent(in) :: dt_over_lz(n_points) -real(kind=real_cvprec), intent(in out) :: q_loc_cond(n_points) -real(kind=real_cvprec), intent(in out) :: n_cond(n_points) -real(kind=real_cvprec), intent(in out) :: r_cond(n_points) -real(kind=real_cvprec), intent(in out) :: wf_cond(n_points) -real(kind=real_cvprec), intent(in out) :: kq_cond(n_points) -real(kind=real_cvprec), intent(in out) :: kt_cond(n_points) -! Outputs need intent inout so that values initialised to zero -! are preserved at points that aren't in the compression list. +real(kind=real_cvprec), intent(out) :: q_loc_cond(n_points) +real(kind=real_cvprec), intent(out) :: n_cond(n_points) +real(kind=real_cvprec), intent(out) :: r_cond(n_points) +real(kind=real_cvprec), intent(out) :: wf_cond(n_points) +real(kind=real_cvprec), intent(out) :: kq_cond(n_points) +real(kind=real_cvprec), intent(out) :: kt_cond(n_points) ! Super-array to store compressed copies of the array arguments real(kind=real_cvprec), allocatable :: args_cmpr(:,:) @@ -279,46 +276,61 @@ subroutine calc_cond_properties_cmpr( n_points, nc, index_ic, & else ! ( REAL(nc,real_cvprec) <= cmpr_thresh * REAL(n_points,real_cvprec) ) - ! If only a minority of points have non-zero mixing ratio... + ! If only a minority (or none) of the points have non-zero mixing ratio... ! Do a compressed call... - ! Allocate compression array - allocate( args_cmpr ( nc, n_args ) ) - - ! Compress inputs - do ic2 = 1, nc - ic = index_ic(ic2) - args_cmpr(ic2,i_ref_temp) = ref_temp(ic) - args_cmpr(ic2,i_q_cond) = q_cond(ic) - args_cmpr(ic2,i_rho_dry) = rho_dry(ic) - args_cmpr(ic2,i_rho_wet) = rho_wet(ic) - args_cmpr(ic2,i_dt_over_lz) = dt_over_lz(ic) + ! Initialise outputs to zero + do ic = 1, n_points + q_loc_cond(ic) = zero + n_cond(ic) = zero + r_cond(ic) = zero + wf_cond(ic) = zero + kq_cond(ic) = zero + kt_cond(ic) = zero end do - ! Call the main routine on the compressed arrays - call calc_cond_properties( nc, & - cond_params, & - args_cmpr(:,i_ref_temp), args_cmpr(:,i_q_cond), & - args_cmpr(:,i_rho_dry), args_cmpr(:,i_rho_wet), & - args_cmpr(:,i_dt_over_lz), & - args_cmpr(:,i_q_loc_cond), args_cmpr(:,i_n_cond), & - args_cmpr(:,i_r_cond), args_cmpr(:,i_wf_cond), & - args_cmpr(:,i_kq_cond), args_cmpr(:,i_kt_cond) ) - - ! Decompress outputs - do ic2 = 1, nc - ic = index_ic(ic2) - q_loc_cond(ic) = args_cmpr(ic2,i_q_loc_cond) - n_cond(ic) = args_cmpr(ic2,i_n_cond) - r_cond(ic) = args_cmpr(ic2,i_r_cond) - wf_cond(ic) = args_cmpr(ic2,i_wf_cond) - kq_cond(ic) = args_cmpr(ic2,i_kq_cond) - kt_cond(ic) = args_cmpr(ic2,i_kt_cond) - end do + if ( nc > 0 ) then + ! If any points have nonzero mixing-ratio (nothing else to do otherwise) + + ! Allocate compression array + allocate( args_cmpr ( nc, n_args ) ) + + ! Compress inputs + do ic2 = 1, nc + ic = index_ic(ic2) + args_cmpr(ic2,i_ref_temp) = ref_temp(ic) + args_cmpr(ic2,i_q_cond) = q_cond(ic) + args_cmpr(ic2,i_rho_dry) = rho_dry(ic) + args_cmpr(ic2,i_rho_wet) = rho_wet(ic) + args_cmpr(ic2,i_dt_over_lz) = dt_over_lz(ic) + end do + + ! Call the main routine on the compressed arrays + call calc_cond_properties( nc, & + cond_params, & + args_cmpr(:,i_ref_temp), args_cmpr(:,i_q_cond), & + args_cmpr(:,i_rho_dry), args_cmpr(:,i_rho_wet), & + args_cmpr(:,i_dt_over_lz), & + args_cmpr(:,i_q_loc_cond), args_cmpr(:,i_n_cond), & + args_cmpr(:,i_r_cond), args_cmpr(:,i_wf_cond), & + args_cmpr(:,i_kq_cond), args_cmpr(:,i_kt_cond) ) + + ! Decompress outputs + do ic2 = 1, nc + ic = index_ic(ic2) + q_loc_cond(ic) = args_cmpr(ic2,i_q_loc_cond) + n_cond(ic) = args_cmpr(ic2,i_n_cond) + r_cond(ic) = args_cmpr(ic2,i_r_cond) + wf_cond(ic) = args_cmpr(ic2,i_wf_cond) + kq_cond(ic) = args_cmpr(ic2,i_kq_cond) + kt_cond(ic) = args_cmpr(ic2,i_kt_cond) + end do + + ! Deallocate + deallocate( args_cmpr ) - ! Deallocate - deallocate( args_cmpr ) + end if ! ( nc > 0 ) end if ! ( REAL(nc,real_cvprec) <= cmpr_thresh * REAL(n_points,real_cvprec) ) diff --git a/science/physics_schemes/source/convection/comorph/microphysics/collision_ctl.F90 b/science/physics_schemes/source/convection/comorph/microphysics/collision_ctl.F90 index ab3558aedf..68c7e18190 100644 --- a/science/physics_schemes/source/convection/comorph/microphysics/collision_ctl.F90 +++ b/science/physics_schemes/source/convection/comorph/microphysics/collision_ctl.F90 @@ -24,9 +24,9 @@ subroutine collision_ctl( n_points, n_points_super, nc, index_ic, & diags_super ) use comorph_constants_mod, only: real_cvprec, zero, & - n_cond_species, n_cond_species_liq, & - cond_params, l_cv_cf, l_cv_graup, & - i_cond_rain, i_cond_cf, i_cond_graup + n_cond_species, n_cond_species_liq, & + cond_params, l_cv_cf, l_cv_graup, & + i_cond_rain, i_cond_cf, i_cond_graup use moist_proc_diags_type_mod, only: moist_proc_diags_type use collision_rate_mod, only: collision_rate_cmpr use lat_heat_mod, only: lat_heat_incr, i_phase_change_frz @@ -138,7 +138,7 @@ subroutine collision_ctl( n_points, n_points_super, nc, index_ic, & logical :: l_collect ! Loop counters -integer :: ic, ic2, i_liq, i_cond, i_super +integer :: ic, ic2, i_cond, i_cond1, i_cond2, i_super if ( l_diags ) then @@ -162,173 +162,167 @@ subroutine collision_ctl( n_points, n_points_super, nc, index_ic, & end if ! ( l_diags ) -! Loop over all species which have non-zero mixing ratio +! Loop over all collected species which have non-zero mixing ratio ! (currently only doing collection of liquids though) -do i_liq = 1, n_cond_species_liq - if ( nc(i_liq) > 0 ) then +do i_cond1 = 1, n_cond_species_liq + if ( nc(i_cond1) > 0 ) then - ! Initialise total collection increment to species i_liq + ! Initialise total collection increment to species i_cond1 ! to zero - do ic2 = 1, nc(i_liq) - ic = index_ic(ic2,i_liq) - dq_col_cond(ic,i_liq) = zero + do ic2 = 1, nc(i_cond1) + ic = index_ic(ic2,i_cond1) + dq_col_cond(ic,i_cond1) = zero end do - ! Initialise flag for whether or not species i_liq has + ! Initialise flag for whether or not species i_cond1 has ! experienced any collection. l_collect = .false. - ! Loop over all species after the current one. - ! Note loop from i_liq+1 ensures all pairs of species + ! Loop over all collecting species after the current one. + ! Note loop from i_cond1+1 ensures all pairs of species ! are compared once but not twice (except we don't ! currently bother with ice-ice collisions). - do i_cond = i_liq+1, n_cond_species + do i_cond2 = i_cond1+1, n_cond_species ! Initialise number of points where collection occurs - nc_col(i_cond) = 0 - if ( nc(i_cond) > 0 ) then + nc_col(i_cond2) = 0 + if ( nc(i_cond2) > 0 ) then ! Find points where both species coincide - if ( nc(i_liq) == n_points ) then - ! Species i_liq present at all points; - ! => coinciding points are just i_cond points - nc_col(i_cond) = nc(i_cond) - index_ic_col(1:nc_col(i_cond),i_cond) & - = index_ic(1:nc_col(i_cond),i_cond) - else if ( nc(i_cond) == n_points ) then - ! Species i_cond present at all points; - ! => coinciding points are just i_liq points - nc_col(i_cond) = nc(i_liq) - index_ic_col(1:nc_col(i_cond),i_cond) & - = index_ic(1:nc_col(i_cond),i_liq) + if ( nc(i_cond1) == n_points ) then + ! Species i_cond1 present at all points; + ! => coinciding points are just i_cond2 points + nc_col(i_cond2) = nc(i_cond2) + index_ic_col(1:nc(i_cond2),i_cond2) = index_ic(1:nc(i_cond2),i_cond2) + else if ( nc(i_cond2) == n_points ) then + ! Species i_cond2 present at all points; + ! => coinciding points are just i_cond1 points + nc_col(i_cond2) = nc(i_cond1) + index_ic_col(1:nc(i_cond1),i_cond2) = index_ic(1:nc(i_cond1),i_cond1) else ! Both species have partial coverage; find overlap - do ic2 = 1, nc(i_liq) - ic = index_ic(ic2,i_liq) - if ( q_cond(ic,i_cond) > zero ) then - nc_col(i_cond) = nc_col(i_cond) + 1 - index_ic_col(nc_col(i_cond),i_cond) = ic + do ic2 = 1, nc(i_cond2) + ic = index_ic(ic2,i_cond2) + if ( q_cond(ic,i_cond1) > zero ) then + nc_col(i_cond2) = nc_col(i_cond2) + 1 + index_ic_col(nc_col(i_cond2),i_cond2) = ic end if end do end if ! If any points: - if ( nc_col(i_cond) > 0 ) then + if ( nc_col(i_cond2) > 0 ) then ! Call routine to calculate increment to species - ! i_liq due to collection by species i_cond + ! i_cond1 due to collection by species i_cond2 call collision_rate_cmpr( n_points, & - nc_col(i_cond), index_ic_col(:,i_cond), & - q_loc_cond(:,i_liq), n_cond(:,i_cond), & - r_cond(:,i_liq), r_cond(:,i_cond), & - wf_cond(:,i_liq), wf_cond(:,i_cond), & + nc_col(i_cond2), index_ic_col(:,i_cond2), & + q_loc_cond(:,i_cond1), n_cond(:,i_cond2), & + r_cond(:,i_cond1), r_cond(:,i_cond2), & + wf_cond(:,i_cond1), wf_cond(:,i_cond2), & delta_t, rho_dry, rho_wet, & - cond_params(i_liq)%pt % rho, & - cond_params(i_cond)%pt % rho, & - cond_params(i_liq)%pt % area_coef, & - cond_params(i_cond)%pt % area_coef, & - dq_col_cond(:,i_cond) ) + cond_params(i_cond1)%pt % rho, & + cond_params(i_cond2)%pt % rho, & + cond_params(i_cond1)%pt % area_coef, & + cond_params(i_cond2)%pt % area_coef, & + dq_col_cond(:,i_cond2) ) ! Add on contribution to total increment to species - ! i_liq due to collection - do ic2 = 1, nc_col(i_cond) - ic = index_ic_col(ic2,i_cond) - dq_col_cond(ic,i_liq) = dq_col_cond(ic,i_liq) & - + dq_col_cond(ic,i_cond) + ! i_cond1 due to collection + do ic2 = 1, nc_col(i_cond2) + ic = index_ic_col(ic2,i_cond2) + dq_col_cond(ic,i_cond1) = dq_col_cond(ic,i_cond1) & + + dq_col_cond(ic,i_cond2) end do ! Set flag indicating some collection has been done l_collect = .true. - end if ! ( nc_col(i_cond) > 0 ) + end if ! ( nc_col(i_cond2) > 0 ) - end if ! ( nc(i_cond) > 0 ) - end do ! i_cond = i_liq+1, n_cond_species + end if ! ( nc(i_cond2) > 0 ) + end do ! i_cond2 = i_cond1+1, n_cond_species ! If any collection was done... if ( l_collect ) then - ! Check to avoid collecting more of species i_liq's + ! Check to avoid collecting more of species i_cond1's ! mixing ratio than actually exists; rescale the increments ! down if this occurs nc_tmp = 0 - do ic2 = 1, nc(i_liq) - ic = index_ic(ic2,i_liq) - ! If the increment to q_cond(i_liq) is greater than its value - if ( dq_col_cond(ic,i_liq) > q_cond(ic,i_liq) ) then + do ic2 = 1, nc(i_cond1) + ic = index_ic(ic2,i_cond1) + ! If the increment to q_cond(i_cond1) is greater than its value + if ( dq_col_cond(ic,i_cond1) > q_cond(ic,i_cond1) ) then ! Save indices of points nc_tmp = nc_tmp + 1 index_ic_tmp(nc_tmp) = ic end if end do - ! If any points were found to have negative q_cond(i_liq): + ! If any points were found to have negative q_cond(i_cond1): if ( nc_tmp > 0 ) then do ic2 = 1, nc_tmp ic = index_ic_tmp(ic2) ! Calculate a scaling factor by which the increments ! need to be reduced to avoid creating negative q_cond ! = pre-existing mixing-ratio q / increment (dq) - factor(ic) = q_cond(ic,i_liq) / dq_col_cond(ic,i_liq) + factor(ic) = q_cond(ic,i_cond1) / dq_col_cond(ic,i_cond1) ! Reset total collection increment to current value - dq_col_cond(ic,i_liq) = q_cond(ic,i_liq) + dq_col_cond(ic,i_cond1) = q_cond(ic,i_cond1) end do ! Reduce the collection increments accordingly at - ! points where we removed too much of species i_liq - do i_cond = i_liq+1, n_cond_species - if ( nc_col(i_cond) > 0 ) then + ! points where we removed too much of species i_cond1 + do i_cond2 = i_cond1+1, n_cond_species + if ( nc_col(i_cond2) > 0 ) then do ic2 = 1, nc_tmp ic = index_ic_tmp(ic2) - dq_col_cond(ic,i_cond) = dq_col_cond(ic,i_cond) & - * factor(ic) + dq_col_cond(ic,i_cond2) = dq_col_cond(ic,i_cond2) * factor(ic) end do end if end do end if ! Decrement the collected species' mixing ratio - do ic2 = 1, nc(i_liq) - ic = index_ic(ic2,i_liq) - q_cond(ic,i_liq) = q_cond(ic,i_liq) & - - dq_col_cond(ic,i_liq) - ! dq_col(i_liq) stores sum of increments from collection + do ic2 = 1, nc(i_cond1) + ic = index_ic(ic2,i_cond1) + q_cond(ic,i_cond1) = q_cond(ic,i_cond1) - dq_col_cond(ic,i_cond1) + ! dq_col(i_cond1) stores sum of increments from collection ! by all other species end do ! Add on increments for each collecting species... - do i_cond = i_liq+1, n_cond_species - if ( nc_col(i_cond) > 0 ) then + do i_cond2 = i_cond1+1, n_cond_species + if ( nc_col(i_cond2) > 0 ) then ! Increment the collecting species' mixing ratio - do ic2 = 1, nc_col(i_cond) - ic = index_ic_col(ic2,i_cond) - q_cond(ic,i_cond) = q_cond(ic,i_cond) & - + dq_col_cond(ic,i_cond) + do ic2 = 1, nc_col(i_cond2) + ic = index_ic_col(ic2,i_cond2) + q_cond(ic,i_cond2) = q_cond(ic,i_cond2) + dq_col_cond(ic,i_cond2) end do ! If this collection process entails a phase-change ! (i.e. riming) - if ( cond_params(i_cond)%pt % l_ice .and. & - (.not. cond_params(i_liq)%pt % l_ice) ) then + if ( cond_params(i_cond2)%pt % l_ice .and. & + (.not. cond_params(i_cond1)%pt % l_ice) ) then ! Add on latent heating - call lat_heat_incr( n_points, nc_col(i_cond), & + call lat_heat_incr( n_points, nc_col(i_cond2), & i_phase_change_frz, & cp_tot, temperature, & - index_ic=index_ic_col(:,i_cond), & - dq=dq_col_cond(:,i_cond) ) - ! Increment the total rate of freezing onto - ! species i_cond - do ic2 = 1, nc_col(i_cond) - ic = index_ic_col(ic2,i_cond) - dq_frz_cond(ic,i_cond) = dq_frz_cond(ic,i_cond) & - + dq_col_cond(ic,i_cond) + index_ic=index_ic_col(:,i_cond2), & + dq=dq_col_cond(:,i_cond2) ) + ! Increment the total rate of freezing onto species i_cond2 + do ic2 = 1, nc_col(i_cond2) + ic = index_ic_col(ic2,i_cond2) + dq_frz_cond(ic,i_cond2) = dq_frz_cond(ic,i_cond2) & + + dq_col_cond(ic,i_cond2) end do end if end if end do - ! If i_liq is rain, and ice & graupel are both on - if ( i_liq == i_cond_rain .and. l_cv_cf & + ! If i_cond1 is rain, and ice & graupel are both on + if ( i_cond1 == i_cond_rain .and. l_cv_cf & .and. l_cv_graup ) then ! If any collisions between rain and ice if ( nc_col(i_cond_cf) > 0 ) then @@ -339,24 +333,22 @@ subroutine collision_ctl( n_points, n_points_super, nc, index_ic, & nc_col(i_cond_cf), index_ic_col(:,i_cond_cf), & nc(i_cond_graup), index_ic(:,i_cond_graup), & dq_col_cond(:,i_cond_cf), & - dq_col_cond(:,i_liq), q_cond(:,i_liq), & - kq_cond(:,i_liq), kt_cond(:,i_liq), & - n_cond(:,i_cond_cf), n_cond(:,i_liq), & - q_loc_cond(:,i_cond_cf), q_loc_cond(:,i_liq), & + dq_col_cond(:,i_cond1), q_cond(:,i_cond1), & + kq_cond(:,i_cond1), kt_cond(:,i_cond1), & + n_cond(:,i_cond_cf), n_cond(:,i_cond1), & + q_loc_cond(:,i_cond_cf), q_loc_cond(:,i_cond1), & q_cond(:,i_cond_cf), q_cond(:,i_cond_graup), & - dq_frz_cond(:,i_cond_cf), & - dq_frz_cond(:,i_cond_graup), & - kq_cond(:,i_cond_graup), & - kt_cond(:,i_cond_graup) ) + dq_frz_cond(:,i_cond_cf), dq_frz_cond(:,i_cond_graup), & + kq_cond(:,i_cond_graup), kt_cond(:,i_cond_graup) ) end if ! ( nc_col(i_cond_cf) > 0 ) - end if ! ( i_liq == i_cond_rain .AND. l_cv_cf + end if ! ( i_cond1 == i_cond_rain .AND. l_cv_cf ! .AND. l_cv_graup ) end if ! ( l_collect ) - end if ! ( nc(i_liq) > 0 ) -end do ! i_liq = 1, n_cond_species_liq + end if ! ( nc(i_cond1) > 0 ) +end do ! i_cond1 = 1, n_cond_species_liq if ( l_diags ) then diff --git a/science/physics_schemes/source/convection/comorph/microphysics/collision_rate.F90 b/science/physics_schemes/source/convection/comorph/microphysics/collision_rate.F90 index 83a58bc599..ad8eb55b11 100644 --- a/science/physics_schemes/source/convection/comorph/microphysics/collision_rate.F90 +++ b/science/physics_schemes/source/convection/comorph/microphysics/collision_rate.F90 @@ -25,10 +25,10 @@ subroutine collision_rate( n_points, & area_coef_1, area_coef_2, & dq_col ) -use comorph_constants_mod, only: pi, coef_wf_spread, col_eff_coef, & - kin_visc, drag_coef_cond, & - real_cvprec, six, half, four_thirds, min_float, & - sqrt_min_float +use comorph_constants_mod, only: real_cvprec, min_float, sqrt_min_float, & + six, half, four_thirds, pi, & + coef_wf_spread, col_eff_coef, & + kin_visc, drag_coef_cond implicit none diff --git a/science/physics_schemes/source/convection/comorph/microphysics/ice_nucleation.F90 b/science/physics_schemes/source/convection/comorph/microphysics/ice_nucleation.F90 index a403a34410..4db8bd97a8 100644 --- a/science/physics_schemes/source/convection/comorph/microphysics/ice_nucleation.F90 +++ b/science/physics_schemes/source/convection/comorph/microphysics/ice_nucleation.F90 @@ -31,7 +31,7 @@ subroutine ice_nucleation( n_points, & n_points_diag, n_diags, diags_super ) use comorph_constants_mod, only: real_cvprec, zero, & - homnuc_temp, hetnuc_temp, coef_hetnuc + homnuc_temp, hetnuc_temp, coef_hetnuc use moist_proc_diags_type_mod, only: moist_proc_diags_type use lat_heat_mod, only: lat_heat_incr, i_phase_change_frz @@ -44,7 +44,7 @@ subroutine ice_nucleation( n_points, & ! Points where the liquid species exists integer, intent(in out) :: nc_liq integer, intent(in out) :: index_ic_liq(n_points) -! (this list get altered if liquid is completely removed +! (this list gets altered if liquid is completely removed ! from any points by freezing) ! Points where the ice species exists @@ -190,7 +190,7 @@ subroutine ice_nucleation( n_points, & if ( l_added_where_none ) then nc_ice = 0 do ic = 1, n_points - if ( q_ice(ic) > 0 ) then + if ( q_ice(ic) > zero ) then nc_ice = nc_ice + 1 index_ic_ice(nc_ice) = ic end if diff --git a/science/physics_schemes/source/convection/comorph/microphysics/ice_rain_to_graupel.F90 b/science/physics_schemes/source/convection/comorph/microphysics/ice_rain_to_graupel.F90 index 15ee26d785..8c2dd2a8fc 100644 --- a/science/physics_schemes/source/convection/comorph/microphysics/ice_rain_to_graupel.F90 +++ b/science/physics_schemes/source/convection/comorph/microphysics/ice_rain_to_graupel.F90 @@ -124,7 +124,7 @@ subroutine ice_rain_to_graupel( n_points, & if ( l_added_where_none ) then nc_graup = 0 do ic = 1, n_points - if ( q_graup(ic) > 0 ) then + if ( q_graup(ic) > zero ) then nc_graup = nc_graup + 1 index_ic_graup(nc_graup) = ic end if diff --git a/science/physics_schemes/source/convection/comorph/microphysics/microphysics_1.F90 b/science/physics_schemes/source/convection/comorph/microphysics/microphysics_1.F90 index be4fd88b48..121d9a9436 100644 --- a/science/physics_schemes/source/convection/comorph/microphysics/microphysics_1.F90 +++ b/science/physics_schemes/source/convection/comorph/microphysics/microphysics_1.F90 @@ -50,8 +50,8 @@ subroutine microphysics_1( n_points, n_points_super, nc, index_ic, & n_points_diag, n_diags, diags_super ) use comorph_constants_mod, only: real_cvprec, cond_params, zero, & - n_cond_species, & - n_cond_species_liq, n_cond_species_ice + n_cond_species, & + n_cond_species_liq, n_cond_species_ice use moist_proc_diags_type_mod, only: moist_proc_diags_type use activate_cond_mod, only: activate_cond @@ -253,32 +253,21 @@ subroutine microphysics_1( n_points, n_points_super, nc, index_ic, & ! Loop over all condensed water species do i_cond = 1, n_cond_species + ! Note we still need to call calc_cond_properties_cmpr even if + ! no points have nonzero mass of this species, just to set the + ! required fields to zero. + + ! Routine implicitly solves the fall-speed / fall-out + ! and number concentration / particle radius relationship, + ! and calculates moisture and heat exchange coefficients + call calc_cond_properties_cmpr( & + n_points, nc(i_cond), index_ic(:,i_cond), & + cond_params(i_cond)%pt, ref_temp, q_cond(:,i_cond), & + rho_dry, rho_wet, dt_over_lz, & + q_loc_cond(:,i_cond), n_cond(:,i_cond), & + r_cond(:,i_cond), wf_cond(:,i_cond), & + kq_cond(:,i_cond), kt_cond(:,i_cond) ) - ! Initialise outputs to zero - do ic = 1, n_points - q_loc_cond(ic,i_cond) = zero - n_cond(ic,i_cond) = zero - r_cond(ic,i_cond) = zero - wf_cond(ic,i_cond) = zero - kq_cond(ic,i_cond) = zero - kt_cond(ic,i_cond) = zero - end do - - ! If any points - if ( nc(i_cond) > 0 ) then - - ! Routine implicitly solves the fall-speed / fall-out - ! and number concentration / particle radius relationship, - ! and calculates moisture and heat exchange coefficients - call calc_cond_properties_cmpr( & - n_points, nc(i_cond), index_ic(:,i_cond), & - cond_params(i_cond)%pt, ref_temp, q_cond(:,i_cond), & - rho_dry, rho_wet, dt_over_lz, & - q_loc_cond(:,i_cond), n_cond(:,i_cond), & - r_cond(:,i_cond), wf_cond(:,i_cond), & - kq_cond(:,i_cond), kt_cond(:,i_cond) ) - - end if end do diff --git a/science/physics_schemes/source/convection/comorph/microphysics/microphysics_2.F90 b/science/physics_schemes/source/convection/comorph/microphysics/microphysics_2.F90 index 9b47a8c4a9..8311694dbd 100644 --- a/science/physics_schemes/source/convection/comorph/microphysics/microphysics_2.F90 +++ b/science/physics_schemes/source/convection/comorph/microphysics/microphysics_2.F90 @@ -18,17 +18,16 @@ module microphysics_2_mod ! ! This includes any microphysical processes that need to be done ! after the implicit solution of phase-changes. -! Current only does autoconversion of liquid cloud to rain +! Currently only does autoconversion of liquid-cloud to rain subroutine microphysics_2( n_points, n_points_super, nc, index_ic, & delta_t, vert_len, wf_cond, q_cond, & l_diags, moist_proc_diags, & n_points_diag, n_diags, diags_super ) -use comorph_constants_mod, only: real_cvprec, zero, one, two, four, & - n_cond_species, i_cond_cl, i_cond_rain, & - indi_thresh, autoc_opt, autoc_linear, & - autoc_quadratic, q_cl_auto, coef_auto +use comorph_constants_mod, only: real_cvprec, n_cond_species, i_cond_cl, & + i_cond_rain, l_cv_rain use moist_proc_diags_type_mod, only: moist_proc_diags_type +use autoconversion_mod, only: autoconversion implicit none @@ -71,222 +70,16 @@ subroutine microphysics_2( n_points, n_points_super, nc, index_ic, & ( n_points_diag, n_diags ) -! Increment to q_rain due to autoconversion -real(kind=real_cvprec) :: dq_auto(n_points) - -! Temporary store for ratio fall-speed / height-interval -real(kind=real_cvprec) :: wf_over_lz - -! Coefficients of quadratic equation -real(kind=real_cvprec) :: a_quad, b_quad, c_quad - -! Solution to quadratic equation for final q_cl -real(kind=real_cvprec) :: q_cl_out - -! Loop counters -integer :: ic, ic2, i_super - - -! Autoconversion of liquid cloud to rain - -! If any liquid cloud present -if ( nc(i_cond_cl) > 0 ) then - - ! If calculation required at majority of points - if ( real(nc(i_cond_cl),real_cvprec) > indi_thresh & - * real(n_points,real_cvprec) ) then - ! Full-field calculation - - ! Implicit discretisation of autoconversion and fall-out: - ! Two choices, linear or quadratic - if ( autoc_opt==autoc_linear ) then - ! ---------------------------------------------------------------- - ! Linear, with threshold - ! - ! d_q_cl_auto = delta_t coef_auto ( q_cl_out - q_cl_auto ) - ! d_q_cl_fall = delta_t q_cl_out wf/lz - ! - ! Therefore - ! - ! q_cl_out = q_cl_in - delta_t ( - ! coef_auto ( q_cl_out - q_cl_auto ) - ! + q_cl_out wf/lz ) - ! - ! Rearranging: - ! - ! q_cl_out - q_cl_auto = q_cl_in - q_cl_auto - delta_t ( - ! coef_auto ( q_cl_out - q_cl_auto ) - ! + ( q_cl_out - q_cl_auto + q_cl_auto ) wf/lz ) - ! - ! => ( q_cl_out - q_cl_auto ) ( 1 + delta_t ( coef_auto + wf/lz ) ) - ! = ( q_cl_in - q_cl_auto ) - delta_t wf/lz q_cl_auto - ! - ! => ( q_cl_out - q_cl_auto ) - ! = ( ( q_cl_in - q_cl_auto ) - delta_t wf/lz q_cl_auto ) - ! / ( 1 + delta_t ( coef_auto + wf/lz ) ) - ! - ! Then multiply this by delta_t * coef_auto to get - ! the autoconversion increment - ! - ! Note: the solution only makes sense where the numerator - ! ( q_cl_in - q_cl_auto ) - delta_t wf/lz q_cl_auto - ! is positive (otherwise, the solution with no autoconversion - ! yields q_cl < q_cl_auto, and we should have no autoconversion - ! as we're below the threshold). - - do ic = 1, n_points - ! Store ratio fall-speed / height-interval - wf_over_lz = wf_cond(ic,i_cond_cl) / vert_len(ic) - ! Compute autoconversion increment, using the above formula - dq_auto(ic) = delta_t(ic) * coef_auto & - * max( ( q_cond(ic,i_cond_cl) - q_cl_auto ) & - - delta_t(ic) * wf_over_lz * q_cl_auto, zero ) & - / ( one + delta_t(ic) * ( coef_auto + wf_over_lz ) ) - end do - - else if ( autoc_opt==autoc_quadratic ) then - ! ----------------------------------------------------------- - ! Quadratic, no threshold - ! - ! d_q_cl_auto = delta_t * coef_auto * q_cl_out^2 - ! d_q_cl_fall = delta_t * q_cl_out * wf/lz - ! - ! Therefore - ! - ! q_cl_out = q_cl_in - delta_t ( coef_auto q_cl_out^2 - ! + q_cl_out wf/lz ) - ! - ! Rearranging: - ! - ! delta_t*coef_auto * q_cl_out^2 + (1+delta_t*wf/lz) * q_cl_out - ! - q_cl_in = 0 - ! - ! => quadratic in q_cl_out = (-b +/- sqrt(b^2-4ac) )/2a - ! Only +ve root will be +ve - ! => q_cl_out = ( -(1+delta_t*wf/lz) + - ! sqrt( (1+delta_t*wf/lz)^2 + 4*delta_t*coef_auto*q_cl_in ) ) / - ! (2*delta_t*coef_auto) - ! - ! Then autoconversion increment is delta_t coef_auto * q_cl_out^2 - ! - ! Note: the solution only makes sense where the numerator - ! is positive. - - do ic = 1, n_points - ! Store coefficients in the quadratic formula - a_quad = delta_t(ic) * coef_auto - b_quad = one + delta_t(ic) * wf_cond(ic,i_cond_cl) / vert_len(ic) - c_quad = - q_cond(ic,i_cond_cl) - ! Compute the solution to the quadratic - q_cl_out = max( zero, & - ( -b_quad + sqrt( b_quad*b_quad - four*a_quad*c_quad ) ) & - /( two*a_quad ) ) - ! Find increment consistent with this - dq_auto(ic) = delta_t(ic) * coef_auto * q_cl_out * q_cl_out - end do - - end if ! ( autoc_opt ) - - ! Add increment to q_cl and q_rain - do ic = 1, n_points - q_cond(ic,i_cond_cl) = q_cond(ic,i_cond_cl) & - - dq_auto(ic) - q_cond(ic,i_cond_rain) = q_cond(ic,i_cond_rain) & - + dq_auto(ic) - end do - - if ( l_diags ) then - ! Save autoconversion increment diagnostics, if requested - if ( moist_proc_diags % diags_cond(i_cond_cl)%pt & - % dq_aut % flag ) then - ! Extract super-array address - i_super = moist_proc_diags % diags_cond(i_cond_cl)%pt & - % dq_aut % i_super - ! Copy increment - do ic = 1, n_points - diags_super(ic,i_super) = diags_super(ic,i_super) & - - dq_auto(ic) - end do - end if - if ( moist_proc_diags % diags_cond(i_cond_rain)%pt & - % dq_aut % flag ) then - ! Extract super-array address - i_super = moist_proc_diags % diags_cond(i_cond_rain)%pt & - % dq_aut % i_super - ! Copy increment - do ic = 1, n_points - diags_super(ic,i_super) = diags_super(ic,i_super) & - + dq_auto(ic) - end do - end if - end if ! ( l_diags ) - - ! Liquid-cloud present at only a small fraction of points - else - ! Compressed version of exactly the same calculation - - if ( autoc_opt==autoc_linear ) then - ! linear, with threshold q_cl_auto - - do ic2 = 1, nc(i_cond_cl) - ic = index_ic(ic2,i_cond_cl) - wf_over_lz = wf_cond(ic,i_cond_cl) / vert_len(ic) - dq_auto(ic) = delta_t(ic) * coef_auto & - * max( ( q_cond(ic,i_cond_cl) - q_cl_auto ) & - - delta_t(ic) * wf_over_lz * q_cl_auto, zero ) & - / ( one + delta_t(ic) * ( coef_auto + wf_over_lz ) ) - end do - - else if ( autoc_opt==autoc_quadratic ) then - ! quadratic, with no threshold q_cl - - do ic2 = 1, nc(i_cond_cl) - ic = index_ic(ic2,i_cond_cl) - a_quad = delta_t(ic) * coef_auto - b_quad = one + delta_t(ic) * wf_cond(ic,i_cond_cl) / vert_len(ic) - c_quad = - q_cond(ic,i_cond_cl) - q_cl_out = max( zero, & - ( -b_quad + sqrt( b_quad*b_quad - four*a_quad*c_quad ) ) & - /( two*a_quad ) ) - dq_auto(ic) = delta_t(ic) * coef_auto * q_cl_out * q_cl_out - end do - - end if ! ( autoc_opt ) - - do ic2 = 1, nc(i_cond_cl) - ic = index_ic(ic2,i_cond_cl) - q_cond(ic,i_cond_cl) = q_cond(ic,i_cond_cl) & - - dq_auto(ic) - q_cond(ic,i_cond_rain) = q_cond(ic,i_cond_rain) & - + dq_auto(ic) - end do - - if ( l_diags ) then - if ( moist_proc_diags % diags_cond(i_cond_cl)%pt & - % dq_aut % flag ) then - i_super = moist_proc_diags % diags_cond(i_cond_cl)%pt & - % dq_aut % i_super - do ic2 = 1, nc(i_cond_cl) - ic = index_ic(ic2,i_cond_cl) - diags_super(ic,i_super) = diags_super(ic,i_super) & - - dq_auto(ic) - end do - end if - if ( moist_proc_diags % diags_cond(i_cond_rain)%pt & - % dq_aut % flag ) then - i_super = moist_proc_diags % diags_cond(i_cond_rain)%pt & - % dq_aut % i_super - do ic2 = 1, nc(i_cond_cl) - ic = index_ic(ic2,i_cond_cl) - diags_super(ic,i_super) = diags_super(ic,i_super) & - + dq_auto(ic) - end do - end if - end if ! ( l_diags ) - - end if ! Test whether to use compressed vs full-field calculation - -end if ! ( nc(i_cond_cl) > 0 ) +if ( nc(i_cond_cl) > 0 .and. l_cv_rain) then + ! If any liquid cloud present + ! Autoconversion of liquid cloud to rain + call autoconversion( n_points, nc(i_cond_cl), index_ic(:,i_cond_cl), & + delta_t, vert_len, wf_cond(:,i_cond_cl), & + q_cond(:,i_cond_cl), q_cond(:,i_cond_rain), & + l_diags, moist_proc_diags % diags_cl, & + moist_proc_diags % diags_rain, & + n_points_diag, n_diags, diags_super ) +end if return diff --git a/science/physics_schemes/source/convection/comorph/moist_proc/fall_in.F90 b/science/physics_schemes/source/convection/comorph/moist_proc/fall_in.F90 index 1cb55be988..545534edb2 100644 --- a/science/physics_schemes/source/convection/comorph/moist_proc/fall_in.F90 +++ b/science/physics_schemes/source/convection/comorph/moist_proc/fall_in.F90 @@ -31,7 +31,7 @@ subroutine fall_in( n_points, nc, index_ic, & l_diags, i_cond, moist_proc_diags, & n_points_diag, n_diags, diags_super ) -use comorph_constants_mod, only: real_cvprec, indi_thresh, one +use comorph_constants_mod, only: real_cvprec, indi_thresh, zero, one use moist_proc_diags_type_mod, only: moist_proc_diags_type implicit none @@ -55,7 +55,8 @@ subroutine fall_in( n_points, nc, index_ic, & ! Inward flux of the current hydrometeor species / kg m-2 s-1 real(kind=real_cvprec), intent(in out) :: flux_cond(n_points) ! (this is an input, but needs intent(inout) as we convert it -! to an increment) +! to an increment, and then reset it to zero ready for +! reusing the array for the fall-out flux) ! Winds and temperature of the air from which the falling-in ! hydrometeros have fallen @@ -176,6 +177,10 @@ subroutine fall_in( n_points, nc, index_ic, & end if end if + ! Zero the flux ready for reusing the array as the fall-out flux + do ic = 1, n_points + flux_cond(ic) = zero + end do ! If non-zero fall-in flux at minority of points, use the ! stored indices to do calculations at only those points @@ -218,6 +223,11 @@ subroutine fall_in( n_points, nc, index_ic, & end if end if + do ic2 = 1, nc + ic = index_ic(ic2) + flux_cond(ic) = zero + end do + end if !( REAL(nc,real_cvprec) > indi_thresh * REAL(n_points,real_cvprec) ) diff --git a/science/physics_schemes/source/convection/comorph/moist_proc/moist_proc.F90 b/science/physics_schemes/source/convection/comorph/moist_proc/moist_proc.F90 index 1d10adec8d..eaf2612fbf 100644 --- a/science/physics_schemes/source/convection/comorph/moist_proc/moist_proc.F90 +++ b/science/physics_schemes/source/convection/comorph/moist_proc/moist_proc.F90 @@ -33,6 +33,8 @@ subroutine moist_proc( n_points, n_points_super, linear_qs, & use linear_qs_mod, only: n_linear_qs_fields, i_ref_temp, & i_qsat_liq_ref, i_dqsatdT_liq use cmpr_type_mod, only: cmpr_type +use fields_type_mod, only: field_min, field_max, & + i_temperature, i_q_vap, i_qc_first use moist_proc_conservation_mod, only: moist_proc_conservation use set_cp_tot_mod, only: set_cp_tot @@ -212,11 +214,13 @@ subroutine moist_proc( n_points, n_points_super, linear_qs, & ! this gets incremented by phase-changes real(kind=real_cvprec) :: cp_tot(n_points) +! Total water mixing ratio q_vap + sum q_cond +real(kind=real_cvprec) :: q_tot(n_points) + ! Dry-density real(kind=real_cvprec) :: rho_dry(n_points) ! Wet density -! (also used to store total-water q_vap + sum q_cond) real(kind=real_cvprec) :: rho_wet(n_points) ! Work array used for conservation checks @@ -230,8 +234,6 @@ subroutine moist_proc( n_points, n_points_super, linear_qs, & ! Name of a field (for error message) character(len=name_length) :: field_name -! Flag for whether field is positive-only -logical :: l_positive ! Description of where we are in the code, for error messages character(len=name_length) :: where_string @@ -241,20 +243,21 @@ subroutine moist_proc( n_points, n_points_super, linear_qs, & ! Check inputs for bad values (NaN, Inf, etc) if ( i_check_bad_values_cmpr > i_check_bad_none ) then - where_string = "Start of moist_proc call for " // & - trim(adjustl(call_string)) - l_positive = .true. + where_string = "Start of moist_proc call for " // trim(adjustl(call_string)) field_name = "temperature" - call check_bad_values_cmpr( cmpr, k, temperature, & - where_string, field_name, l_positive ) + call check_bad_values_cmpr( cmpr, k, temperature, where_string, field_name, & + field_min=field_min(i_temperature), & + field_max=field_max(i_temperature) ) field_name = "q_vap" - call check_bad_values_cmpr( cmpr, k, q_vap, & - where_string, field_name, l_positive ) + call check_bad_values_cmpr( cmpr, k, q_vap, where_string, field_name, & + field_min=field_min(i_q_vap), & + field_max=field_max(i_q_vap) ) do i_cond = 1, n_cond_species - field_name = "q_" // & - trim(adjustl( cond_params(i_cond)%pt % cond_name )) + field_name = "q_" // trim(adjustl( cond_params(i_cond)%pt % cond_name )) call check_bad_values_cmpr( cmpr, k, q_cond(:,i_cond), & - where_string, field_name, l_positive ) + where_string, field_name, & + field_min=field_min(i_qc_first-1+i_cond), & + field_max=field_max(i_qc_first-1+i_cond) ) end do end if @@ -277,9 +280,18 @@ subroutine moist_proc( n_points, n_points_super, linear_qs, & call set_cp_tot( n_points, n_points_super, & q_vap, q_cond, cp_tot ) -! Set total-water mixing-ratio (stored in rho_wet) +! Set total-water mixing-ratio call calc_q_tot( n_points, n_points_super, & - q_vap, q_cond, rho_wet ) + q_vap, q_cond, q_tot ) + +! Calculate air dry-density +call calc_rho_dry( n_points, temperature, q_vap, pressure, & + rho_dry ) + +! Calculate rho_wet +do ic = 1, n_points + rho_wet(ic) = rho_dry(ic) * ( one + q_tot(ic) ) +end do !---------------------------------------------------------------- @@ -308,10 +320,11 @@ subroutine moist_proc( n_points, n_points_super, linear_qs, & dt_over_rhod_lz, flux_cond(:,i_cond), & fallin_wind_u, fallin_wind_v, fallin_wind_w, & fallin_temp, & - cp_tot, rho_wet, wind_u, wind_v, wind_w, & + cp_tot, q_tot, wind_u, wind_v, wind_w, & temperature, q_cond(:,i_cond), & l_diags, i_cond, moist_proc_diags, & n_points_diag, n_diags, diags_super ) + ! This routine resets flux_cond to zero ready for storing the fall-out flux end if @@ -335,16 +348,6 @@ subroutine moist_proc( n_points, n_points_super, linear_qs, & end do end do -! Calculate air dry-density -call calc_rho_dry( n_points, temperature, q_vap, pressure, & - rho_dry ) - -! Complete calculation of rho_wet (currently stores q_tot; -! dry-mass to wet-mass conversion factor = 1 + q_tot) -do ic = 1, n_points - rho_wet(ic) = rho_dry(ic) * ( one + rho_wet(ic) ) -end do - ! Call microphysics routine: ! Back-Of-The-Envelope Microphysics Scheme (BOTEMS) call microphysics_1( n_points, n_points_super, nc, index_ic, & @@ -359,7 +362,7 @@ subroutine moist_proc( n_points, n_points_super, linear_qs, & ! Note: the above call to microphysics_1 needs to be done even ! if all the condensed water species are currently zero -! everywhere, as the test for activation of new condesation +! everywhere, as the test for activation of new condensation ! is done inside microphysics_1. ! If any condensed water species are non-zero anywhere @@ -481,20 +484,21 @@ subroutine moist_proc( n_points, n_points_super, linear_qs, & ! Check outputs for bad values (NaN, Inf, etc) if ( i_check_bad_values_cmpr > i_check_bad_none ) then - where_string = "End of moist_proc call for " // & - trim(adjustl(call_string)) - l_positive = .true. + where_string = "End of moist_proc call for " // trim(adjustl(call_string)) field_name = "temperature" - call check_bad_values_cmpr( cmpr, k, temperature, & - where_string, field_name, l_positive ) + call check_bad_values_cmpr( cmpr, k, temperature, where_string, field_name, & + field_min=field_min(i_temperature), & + field_max=field_max(i_temperature) ) field_name = "q_vap" - call check_bad_values_cmpr( cmpr, k, q_vap, & - where_string, field_name, l_positive ) + call check_bad_values_cmpr( cmpr, k, q_vap, where_string, field_name, & + field_min=field_min(i_q_vap), & + field_max=field_max(i_q_vap) ) do i_cond = 1, n_cond_species - field_name = "q_" // & - trim(adjustl( cond_params(i_cond)%pt % cond_name )) + field_name = "q_" // trim(adjustl( cond_params(i_cond)%pt % cond_name )) call check_bad_values_cmpr( cmpr, k, q_cond(:,i_cond), & - where_string, field_name, l_positive ) + where_string, field_name, & + field_min=field_min(i_qc_first-1+i_cond), & + field_max=field_max(i_qc_first-1+i_cond) ) end do end if diff --git a/science/physics_schemes/source/convection/comorph/moist_proc/phase_change_solve.F90 b/science/physics_schemes/source/convection/comorph/moist_proc/phase_change_solve.F90 index 30eb3db75b..a896c0d279 100644 --- a/science/physics_schemes/source/convection/comorph/moist_proc/phase_change_solve.F90 +++ b/science/physics_schemes/source/convection/comorph/moist_proc/phase_change_solve.F90 @@ -47,6 +47,8 @@ subroutine phase_change_solve( n_points, n_points_super, & i_dqsatdT_liq, i_dqsatdT_ice use phase_change_coefs_mod, only: n_coefs use moist_proc_diags_type_mod, only: moist_proc_diags_type +use fields_type_mod, only: field_min, field_max, & + i_temperature, i_q_vap, i_qc_first use cmpr_type_mod, only: cmpr_type use lat_heat_mod, only: set_l_con, set_l_sub, set_l_fus, & @@ -227,8 +229,6 @@ subroutine phase_change_solve( n_points, n_points_super, & ! Name of a field (for error message) character(len=name_length) :: field_name -! Flag for whether field is positive-only -logical :: l_positive ! Description of where we are in the code, for error messages character(len=name_length) :: where_string @@ -238,50 +238,50 @@ subroutine phase_change_solve( n_points, n_points_super, & ! Check inputs for bad values (NaN, Inf, etc) if ( i_check_bad_values_cmpr > i_check_bad_none ) then - where_string = "Start of phase_change_solve call for " // & + where_string = "Start of phase_change_solve call for " // & trim(adjustl(call_string)) - l_positive = .true. field_name = "temperature" - call check_bad_values_cmpr( cmpr, k, temperature, & - where_string, field_name, l_positive ) + call check_bad_values_cmpr( cmpr, k, temperature, where_string, field_name, & + field_min=field_min(i_temperature), & + field_max=field_max(i_temperature) ) field_name = "q_vap" - call check_bad_values_cmpr( cmpr, k, q_vap, & - where_string, field_name, l_positive ) + call check_bad_values_cmpr( cmpr, k, q_vap, where_string, field_name, & + field_min=field_min(i_q_vap), & + field_max=field_max(i_q_vap) ) do i_cond = 1, n_cond_species - field_name = "q_" // & - trim(adjustl( cond_params(i_cond)%pt % cond_name )) + field_name = "q_" // trim(adjustl( cond_params(i_cond)%pt % cond_name )) call check_bad_values_cmpr( cmpr, k, q_cond(:,i_cond), & - where_string, field_name, l_positive ) + where_string, field_name, & + field_min=field_min(i_qc_first-1+i_cond), & + field_max=field_max(i_qc_first-1+i_cond) & + *100.0_real_cvprec ) end do do i_cond = 1, n_cond_species - field_name = "kq_" // & - trim(adjustl( cond_params(i_cond)%pt % cond_name )) + field_name = "kq_" // trim(adjustl( cond_params(i_cond)%pt % cond_name )) call check_bad_values_cmpr( cmpr, k, kq_cond(:,i_cond), & - where_string, field_name, l_positive ) + where_string, field_name, field_min=zero ) end do do i_cond = 1, n_cond_species - field_name = "kt_" // & - trim(adjustl( cond_params(i_cond)%pt % cond_name )) + field_name = "kt_" // trim(adjustl( cond_params(i_cond)%pt % cond_name )) call check_bad_values_cmpr( cmpr, k, kt_cond(:,i_cond), & - where_string, field_name, l_positive ) + where_string, field_name, field_min=zero ) end do do i_cond = 1, n_cond_species - field_name = "wf_" // & - trim(adjustl( cond_params(i_cond)%pt % cond_name )) + field_name = "wf_" // trim(adjustl( cond_params(i_cond)%pt % cond_name )) call check_bad_values_cmpr( cmpr, k, wf_cond(:,i_cond), & - where_string, field_name, l_positive ) + where_string, field_name, field_min=zero ) end do do i_cond = 1, n_cond_species - field_name = "q_loc_" // & - trim(adjustl( cond_params(i_cond)%pt % cond_name )) + field_name = "q_loc_" // trim(adjustl( cond_params(i_cond)%pt%cond_name )) call check_bad_values_cmpr( cmpr, k, q_loc_cond(:,i_cond), & - where_string, field_name, l_positive ) + where_string, field_name, & + field_min=field_min(i_qc_first-1+i_cond), & + field_max=field_max(i_qc_first-1+i_cond) ) end do do i_ice = n_cond_species_liq+1, n_cond_species - field_name = "dq_frz_" // & - trim(adjustl( cond_params(i_ice)%pt % cond_name )) + field_name = "dq_frz_" // trim(adjustl( cond_params(i_ice)%pt%cond_name )) call check_bad_values_cmpr( cmpr, k, dq_frz_cond(:,i_ice), & - where_string, field_name, l_positive ) + where_string, field_name, field_min=zero ) end do end if @@ -394,28 +394,25 @@ subroutine phase_change_solve( n_points, n_points_super, & where_string = "phase_change_solve, call for " // & trim(adjustl(call_string)) // & "; after calc_phase_change_coefs" - l_positive = .false. ! Check condensation and melting coefficients do i_cond = 1, n_cond_species do i_coef = 1, n_coefs write(field_name,"(I1)") i_coef field_name = "coefs_cond_" // & - trim(adjustl( cond_params(i_cond)%pt % cond_name )) // & - "_c" // trim(adjustl( field_name )) - call check_bad_values_cmpr( cmpr, k, & - coefs_cond(:,i_coef,i_cond), & - where_string, field_name, l_positive ) + trim(adjustl( cond_params(i_cond)%pt % cond_name )) // & + "_c" // trim(adjustl( field_name )) + call check_bad_values_cmpr( cmpr, k, coefs_cond(:,i_coef,i_cond), & + where_string, field_name ) end do end do do i_ice = n_cond_species_liq+1, n_cond_species do i_coef = 1, n_coefs write(field_name,"(I1)") i_coef field_name = "coefs_melt_" // & - trim(adjustl( cond_params(i_ice)%pt % cond_name )) // & - "_c" // trim(adjustl( field_name )) - call check_bad_values_cmpr( cmpr, k, & - coefs_melt(:,i_coef,i_ice), & - where_string, field_name, l_positive ) + trim(adjustl( cond_params(i_ice)%pt % cond_name )) // & + "_c" // trim(adjustl( field_name )) + call check_bad_values_cmpr( cmpr, k, coefs_melt(:,i_coef,i_ice), & + where_string, field_name ) end do end do @@ -467,28 +464,25 @@ subroutine phase_change_solve( n_points, n_points_super, & where_string = "phase_change_solve, call for " // & trim(adjustl(call_string)) // & "; after melt_ctl" - l_positive = .false. ! Check condensation and melting coefficients do i_cond = 1, n_cond_species do i_coef = 1, n_coefs write(field_name,"(I1)") i_coef field_name = "coefs_cond_" // & - trim(adjustl( cond_params(i_cond)%pt % cond_name )) // & - "_c" // trim(adjustl( field_name )) - call check_bad_values_cmpr( cmpr, k, & - coefs_cond(:,i_coef,i_cond), & - where_string, field_name, l_positive ) + trim(adjustl( cond_params(i_cond)%pt % cond_name )) // & + "_c" // trim(adjustl( field_name )) + call check_bad_values_cmpr( cmpr, k, coefs_cond(:,i_coef,i_cond), & + where_string, field_name ) end do end do do i_ice = n_cond_species_liq+1, n_cond_species do i_coef = 1, n_coefs write(field_name,"(I1)") i_coef field_name = "coefs_melt_" // & - trim(adjustl( cond_params(i_ice)%pt % cond_name )) // & - "_c" // trim(adjustl( field_name )) - call check_bad_values_cmpr( cmpr, k, & - coefs_melt(:,i_coef,i_ice), & - where_string, field_name, l_positive ) + trim(adjustl( cond_params(i_ice)%pt % cond_name )) // & + "_c" // trim(adjustl( field_name )) + call check_bad_values_cmpr( cmpr, k, coefs_melt(:,i_coef,i_ice), & + where_string, field_name ) end do end do @@ -793,20 +787,23 @@ subroutine phase_change_solve( n_points, n_points_super, & if ( i_check_bad_values_cmpr > i_check_bad_none ) then where_string = "End of phase_change_solve call for " // & trim(adjustl(call_string)) - l_positive = .true. ! Check temperature, vapour and condensed water mixing-ratios field_name = "temperature" - call check_bad_values_cmpr( cmpr, k, temperature, & - where_string, field_name, l_positive ) + call check_bad_values_cmpr( cmpr, k, temperature, where_string, field_name, & + field_min=field_min(i_temperature), & + field_max=field_max(i_temperature) ) field_name = "q_vap" - call check_bad_values_cmpr( cmpr, k, q_vap, & - where_string, field_name, l_positive ) + call check_bad_values_cmpr( cmpr, k, q_vap, where_string, field_name, & + field_min=field_min(i_q_vap), & + field_max=field_max(i_q_vap) ) do i_cond = 1, n_cond_species - field_name = "q_" // & - trim(adjustl( cond_params(i_cond)%pt % cond_name )) + field_name = "q_" // trim(adjustl( cond_params(i_cond)%pt % cond_name )) call check_bad_values_cmpr( cmpr, k, q_cond(:,i_cond), & - where_string, field_name, l_positive ) + where_string, field_name, & + field_min=field_min(i_qc_first-1+i_cond), & + field_max=field_max(i_qc_first-1+i_cond) & + *100.0_real_cvprec ) end do ! Check condensation and melting increments @@ -814,14 +811,13 @@ subroutine phase_change_solve( n_points, n_points_super, & field_name = "dq_melt_" // & trim(adjustl( cond_params(i_ice)%pt % cond_name )) call check_bad_values_cmpr( cmpr, k, dq_melt(:,i_ice), & - where_string, field_name, l_positive ) + where_string, field_name, field_min=zero ) end do - l_positive = .false. do i_cond = 1, n_cond_species field_name = "dq_cond_" // & trim(adjustl( cond_params(i_cond)%pt % cond_name )) call check_bad_values_cmpr( cmpr, k, dq_cond(:,i_cond), & - where_string, field_name, l_positive ) + where_string, field_name ) end do ! Check condensation and melting coefficeints @@ -829,22 +825,20 @@ subroutine phase_change_solve( n_points, n_points_super, & do i_coef = 1, n_coefs write(field_name,"(I1)") i_coef field_name = "coefs_cond_" // & - trim(adjustl( cond_params(i_cond)%pt % cond_name )) // & - "_c" // trim(adjustl( field_name )) - call check_bad_values_cmpr( cmpr, k, & - coefs_cond(:,i_coef,i_cond), & - where_string, field_name, l_positive ) + trim(adjustl( cond_params(i_cond)%pt % cond_name )) // & + "_c" // trim(adjustl( field_name )) + call check_bad_values_cmpr( cmpr, k, coefs_cond(:,i_coef,i_cond), & + where_string, field_name ) end do end do do i_ice = n_cond_species_liq+1, n_cond_species do i_coef = 1, n_coefs write(field_name,"(I1)") i_coef field_name = "coefs_melt_" // & - trim(adjustl( cond_params(i_ice)%pt % cond_name )) // & - "_c" // trim(adjustl( field_name )) - call check_bad_values_cmpr( cmpr, k, & - coefs_melt(:,i_coef,i_ice), & - where_string, field_name, l_positive ) + trim(adjustl( cond_params(i_ice)%pt % cond_name )) // & + "_c" // trim(adjustl( field_name )) + call check_bad_values_cmpr( cmpr, k, coefs_melt(:,i_coef,i_ice), & + where_string, field_name ) end do end do diff --git a/science/physics_schemes/source/convection/comorph/moist_thermo/calc_q_tot.F90 b/science/physics_schemes/source/convection/comorph/moist_thermo/calc_q_tot.F90 index a256ec3030..c55b78d328 100644 --- a/science/physics_schemes/source/convection/comorph/moist_thermo/calc_q_tot.F90 +++ b/science/physics_schemes/source/convection/comorph/moist_thermo/calc_q_tot.F90 @@ -76,16 +76,16 @@ end subroutine calc_q_tot !---------------------------------------------------------------- -! Version for full 3-D arrays +! Version for full 2-D arrays !---------------------------------------------------------------- -subroutine calc_q_tot_3d( lb_v, ub_v, q_vap, lb_l, ub_l, q_cl, & +subroutine calc_q_tot_2d( lb_v, ub_v, q_vap, lb_l, ub_l, q_cl, & lb_r, ub_r, q_rain, lb_f, ub_f, q_cf, & lb_s, ub_s, q_snow, lb_g, ub_g, q_graup, & q_tot ) use comorph_constants_mod, only: real_hmprec, & - l_cv_rain, l_cv_cf, l_cv_snow, l_cv_graup, & - nx_full, ny_full, k_bot_conv, k_top_conv + l_cv_rain, l_cv_cf, l_cv_snow, l_cv_graup, & + nx_full, ny_full implicit none @@ -96,49 +96,47 @@ subroutine calc_q_tot_3d( lb_v, ub_v, q_vap, lb_l, ub_l, q_cl, & ! indices of the arrays regardless of whether or not they ! have halos. The lower and upper bounds are passed in through the ! argument list in the lb_* and ub_* integer arrays. Those storing the -! bounds for 3D arrays must have 3 elements; one for each dimension +! bounds for 2D arrays must have 2 elements; one for each dimension ! of the array. ! Water vapour mixing ratio -integer, intent(in) :: lb_v(3), ub_v(3) +integer, intent(in) :: lb_v(2), ub_v(2) real(kind=real_hmprec), intent(in) :: q_vap & - ( lb_v(1):ub_v(1), lb_v(2):ub_v(2), lb_v(3):ub_v(3) ) + ( lb_v(1):ub_v(1), lb_v(2):ub_v(2) ) ! Condensed water species mixing ratios -integer, intent(in) :: lb_l(3), ub_l(3) +integer, intent(in) :: lb_l(2), ub_l(2) real(kind=real_hmprec), intent(in) :: q_cl & - ( lb_l(1):ub_l(1), lb_l(2):ub_l(2), lb_l(3):ub_l(3) ) -integer, intent(in) :: lb_r(3), ub_r(3) + ( lb_l(1):ub_l(1), lb_l(2):ub_l(2) ) +integer, intent(in) :: lb_r(2), ub_r(2) real(kind=real_hmprec), intent(in) :: q_rain & - ( lb_r(1):ub_r(1), lb_r(2):ub_r(2), lb_r(3):ub_r(3) ) -integer, intent(in) :: lb_f(3), ub_f(3) + ( lb_r(1):ub_r(1), lb_r(2):ub_r(2) ) +integer, intent(in) :: lb_f(2), ub_f(2) real(kind=real_hmprec), intent(in) :: q_cf & - ( lb_f(1):ub_f(1), lb_f(2):ub_f(2), lb_f(3):ub_f(3) ) -integer, intent(in) :: lb_s(3), ub_s(3) + ( lb_f(1):ub_f(1), lb_f(2):ub_f(2) ) +integer, intent(in) :: lb_s(2), ub_s(2) real(kind=real_hmprec), intent(in) :: q_snow & - ( lb_s(1):ub_s(1), lb_s(2):ub_s(2), lb_s(3):ub_s(3) ) -integer, intent(in) :: lb_g(3), ub_g(3) + ( lb_s(1):ub_s(1), lb_s(2):ub_s(2) ) +integer, intent(in) :: lb_g(2), ub_g(2) real(kind=real_hmprec), intent(in) :: q_graup & - ( lb_g(1):ub_g(1), lb_g(2):ub_g(2), lb_g(3):ub_g(3) ) + ( lb_g(1):ub_g(1), lb_g(2):ub_g(2) ) ! Output dry-mass to wet-mass conversion factor real(kind=real_hmprec), intent(out) :: q_tot & - ( nx_full, ny_full, k_bot_conv:k_top_conv ) + ( nx_full, ny_full ) ! Loop counters -integer :: i, j, k +integer :: i, j -!$OMP PARALLEL DEFAULT(NONE) PRIVATE( i, j, k ) & -!$OMP SHARED( nx_full, ny_full, l_cv_snow, k_bot_conv, k_top_conv, & +!$OMP PARALLEL DEFAULT(NONE) PRIVATE( i, j ) & +!$OMP SHARED( nx_full, ny_full, l_cv_snow, & !$OMP q_tot, q_cl, q_rain, q_cf, q_snow, q_graup, q_vap ) ! Initialise using liquid cloud !$OMP DO SCHEDULE(STATIC) -do k = k_bot_conv, k_top_conv - do j = 1, ny_full - do i = 1, nx_full - q_tot(i,j,k) = q_cl(i,j,k) - end do +do j = 1, ny_full + do i = 1, nx_full + q_tot(i,j) = q_cl(i,j) end do end do !$OMP END DO @@ -147,11 +145,9 @@ subroutine calc_q_tot_3d( lb_v, ub_v, q_vap, lb_l, ub_l, q_cl, & if ( l_cv_rain ) then !$OMP DO SCHEDULE(STATIC) - do k = k_bot_conv, k_top_conv - do j = 1, ny_full - do i = 1, nx_full - q_tot(i,j,k) = q_tot(i,j,k) + q_rain(i,j,k) - end do + do j = 1, ny_full + do i = 1, nx_full + q_tot(i,j) = q_tot(i,j) + q_rain(i,j) end do end do !$OMP END DO @@ -159,11 +155,9 @@ subroutine calc_q_tot_3d( lb_v, ub_v, q_vap, lb_l, ub_l, q_cl, & if ( l_cv_cf ) then !$OMP DO SCHEDULE(STATIC) - do k = k_bot_conv, k_top_conv - do j = 1, ny_full - do i = 1, nx_full - q_tot(i,j,k) = q_tot(i,j,k) + q_cf(i,j,k) - end do + do j = 1, ny_full + do i = 1, nx_full + q_tot(i,j) = q_tot(i,j) + q_cf(i,j) end do end do !$OMP END DO @@ -171,11 +165,9 @@ subroutine calc_q_tot_3d( lb_v, ub_v, q_vap, lb_l, ub_l, q_cl, & if ( l_cv_snow ) then !$OMP DO SCHEDULE(STATIC) - do k = k_bot_conv, k_top_conv - do j = 1, ny_full - do i = 1, nx_full - q_tot(i,j,k) = q_tot(i,j,k) + q_snow(i,j,k) - end do + do j = 1, ny_full + do i = 1, nx_full + q_tot(i,j) = q_tot(i,j) + q_snow(i,j) end do end do !$OMP END DO @@ -183,11 +175,9 @@ subroutine calc_q_tot_3d( lb_v, ub_v, q_vap, lb_l, ub_l, q_cl, & if ( l_cv_graup ) then !$OMP DO SCHEDULE(STATIC) - do k = k_bot_conv, k_top_conv - do j = 1, ny_full - do i = 1, nx_full - q_tot(i,j,k) = q_tot(i,j,k) + q_graup(i,j,k) - end do + do j = 1, ny_full + do i = 1, nx_full + q_tot(i,j) = q_tot(i,j) + q_graup(i,j) end do end do !$OMP END DO @@ -195,11 +185,9 @@ subroutine calc_q_tot_3d( lb_v, ub_v, q_vap, lb_l, ub_l, q_cl, & ! Add on water-vapour !$OMP DO SCHEDULE(STATIC) -do k = k_bot_conv, k_top_conv - do j = 1, ny_full - do i = 1, nx_full - q_tot(i,j,k) = q_tot(i,j,k) + q_vap(i,j,k) - end do +do j = 1, ny_full + do i = 1, nx_full + q_tot(i,j) = q_tot(i,j) + q_vap(i,j) end do end do !$OMP END DO @@ -207,7 +195,7 @@ subroutine calc_q_tot_3d( lb_v, ub_v, q_vap, lb_l, ub_l, q_cl, & !$OMP END PARALLEL return -end subroutine calc_q_tot_3d +end subroutine calc_q_tot_2d end module calc_q_tot_mod diff --git a/science/physics_schemes/source/convection/comorph/moist_thermo/calc_virt_temp.F90 b/science/physics_schemes/source/convection/comorph/moist_thermo/calc_virt_temp.F90 index 4bbd9192a0..713bd9f9df 100644 --- a/science/physics_schemes/source/convection/comorph/moist_thermo/calc_virt_temp.F90 +++ b/science/physics_schemes/source/convection/comorph/moist_thermo/calc_virt_temp.F90 @@ -54,7 +54,7 @@ subroutine calc_virt_temp( n_points, n_points_super, & real(kind=real_cvprec) :: q_tot(n_points) ! Loop counters -integer :: i +integer :: ic ! In terms of mixing ratios qv, qc, virtual temperature is: ! Tv = T ( 1 + Rv/Rd qv ) / ( 1 + qv + qc ) @@ -69,9 +69,9 @@ subroutine calc_virt_temp( n_points, n_points_super, & call calc_q_tot( n_points, n_points_super, & q_vap, q_cond_super, q_tot ) -! Complete the formula for virtual temperature -do i = 1, n_points - virt_temp(i) = virt_temp(i) / ( one + q_tot(i) ) +! Complete the formula for virtual temperature including water loading +do ic = 1, n_points + virt_temp(ic) = virt_temp(ic) / ( one + q_tot(ic) ) end do return @@ -88,10 +88,10 @@ subroutine calc_virt_temp_3d( lb_T, ub_T, temperature, & virt_temp ) use comorph_constants_mod, only: real_hmprec, & - nx_full, ny_full, k_bot_conv, k_top_conv + nx_full, ny_full, k_bot_conv, k_top_conv use calc_virt_temp_dry_mod, only: calc_virt_temp_dry_3d -use calc_q_tot_mod, only: calc_q_tot_3d +use calc_q_tot_mod, only: calc_q_tot_2d implicit none @@ -138,7 +138,7 @@ subroutine calc_virt_temp_3d( lb_T, ub_T, temperature, & ! Work variable: total-water mixing-ratio real(kind=real_hmprec) :: q_tot & - ( nx_full, ny_full, k_bot_conv:k_top_conv ) + ( nx_full, ny_full ) real(kind=real_hmprec), parameter :: one = 1.0_real_hmprec @@ -154,19 +154,23 @@ subroutine calc_virt_temp_3d( lb_T, ub_T, temperature, & call calc_virt_temp_dry_3d( lb_T, ub_T, temperature, lb_v, ub_v, q_vap, & virt_temp ) -! Calculate the total-water qv + qc -call calc_q_tot_3d( lb_v, ub_v, q_vap, lb_l, ub_l, q_cl, & - lb_r, ub_r, q_rain, lb_f, ub_f, q_cf, & - lb_s, ub_s, q_snow, lb_g, ub_g, q_graup, & - q_tot ) - -! Complete the formula for virtual temperature -!$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) PRIVATE( i, j, k ) & -!$OMP SHARED( nx_full, ny_full, k_bot_conv, k_top_conv, virt_temp, q_tot ) +!$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) PRIVATE( i, j, k, q_tot ) & +!$OMP SHARED( nx_full, ny_full, k_bot_conv, k_top_conv, virt_temp, & +!$OMP lb_v, ub_v, q_vap, lb_l, ub_l, q_cl, lb_r, ub_r, q_rain, & +!$OMP lb_f, ub_f, q_cf, lb_s, ub_s, q_snow, lb_g, ub_g, q_graup ) do k = k_bot_conv, k_top_conv + ! Calculate the total-water qv + qc + call calc_q_tot_2d( lb_v(1:2), ub_v(1:2), q_vap(:,:,k), & + lb_l(1:2), ub_l(1:2), q_cl(:,:,k), & + lb_r(1:2), ub_r(1:2), q_rain(:,:,k), & + lb_f(1:2), ub_f(1:2), q_cf(:,:,k), & + lb_s(1:2), ub_s(1:2), q_snow(:,:,k), & + lb_g(1:2), ub_g(1:2), q_graup(:,:,k), & + q_tot ) + ! Complete the formula for virtual temperature including water loading do j = 1, ny_full do i = 1, nx_full - virt_temp(i,j,k) = virt_temp(i,j,k) / (one + q_tot(i,j,k)) + virt_temp(i,j,k) = virt_temp(i,j,k) / (one + q_tot(i,j)) end do end do end do diff --git a/science/physics_schemes/source/convection/comorph/moist_thermo/set_dqsatdt.F90 b/science/physics_schemes/source/convection/comorph/moist_thermo/set_dqsatdt.F90 index c1621602c9..49659ee336 100644 --- a/science/physics_schemes/source/convection/comorph/moist_thermo/set_dqsatdt.F90 +++ b/science/physics_schemes/source/convection/comorph/moist_thermo/set_dqsatdt.F90 @@ -16,7 +16,28 @@ module set_dqsatdt_mod ! Estimates the gradient of the saturation water vapour mixing ! ratio curve as a function of temperature T, ! by assuming d/dT of the saturation vapour pressure follows -! the Claussius-Clapeyron equation +! the Claussius-Clapeyron equation: +! +! des/dT = Lc es / (Rv T^2) +! +! p = pd + ev +! qv = rhov/rhod +! pd = rhod Rd T +! ev = rhov Rv T +! => ev/pd = Rv/Rd qv +! => qv = Rd/Rv ev / (p - ev) +! +! Assuming constant total-pressure p, +! +! dqs/dT = dqs/des des/dT +! = Rd/Rv ( 1/(p - es) + es/(p - es)^2 ) des/dT +! = Rd/Rv p/(p - es)^2 des/dT +! = Rd/Rv p/(p - es)^2 Lc es / (Rv T^2) +! = qs p/(p - es) Lc / (Rv T^2) +! = qs ( rhod Rd + rhov Rv ) / ( rhod Rd ) Lc / (Rv T^2) +! = qs ( 1 + Rv/Rd qs ) Lc / (Rv T^2) +! +! This formula is used below... !---------------------------------------------------------------- ! Routine for liquid at all temperatures diff --git a/science/physics_schemes/source/convection/comorph/plume_genesis/calc_env_region_tq_nb.F90 b/science/physics_schemes/source/convection/comorph/plume_genesis/calc_env_region_tq_nb.F90 index a3dc6250ab..5422baaec7 100644 --- a/science/physics_schemes/source/convection/comorph/plume_genesis/calc_env_region_tq_nb.F90 +++ b/science/physics_schemes/source/convection/comorph/plume_genesis/calc_env_region_tq_nb.F90 @@ -108,6 +108,8 @@ subroutine calc_env_region_tq_nb( n_points, n_points_super, & ! Local total condensed water within each region real(kind=real_cvprec) :: qc_tot_loc(n_points,n_regions) +! Condensate water-loading excess of each region relative to the grid-mean +real(kind=real_cvprec) :: qc_excess(n_points,n_regions) ! Local total condensed water outside the liquid-cloud real(kind=real_cvprec) :: qc_tot_noliq @@ -188,6 +190,13 @@ subroutine calc_env_region_tq_nb( n_points, n_points_super, & end select end do ! i_cond = 1, n_cond_species +! Set water-loading excess of each region used in buoyancy calculations +do i_region = 1, n_regions + do ic = 1, n_points + qc_excess(ic,i_region) = qc_tot_loc(ic,i_region) - qc_tot(ic) + end do +end do + ! Set frequently used ratio dTv/dT / dqsat/dT do ic = 1, n_points dtv_dqsat_liq(ic) = dtv_dt(ic) / dqsatdt_liq(ic) @@ -252,7 +261,7 @@ subroutine calc_env_region_tq_nb( n_points, n_points_super, & q_vap_r(ic,i_liq) = & ( ( qsat_liq(ic) + supersat(ic) ) * dtv_dqsat_liq(ic) & + q_vap(ic) * dtv_dqv(ic) & - - ( qc_tot_loc(ic,i_liq) - qc_tot(ic) ) * dtv_dqc(ic) ) & + - qc_excess(ic,i_liq) * dtv_dqc(ic) ) & / ( dtv_dqsat_liq(ic) + dtv_dqv(ic) ) end do @@ -282,7 +291,7 @@ subroutine calc_env_region_tq_nb( n_points, n_points_super, & q_vap_r(ic,i_mph) = & ( ( qsat_liq(ic) + supersat(ic) ) * dtv_dqsat_liq(ic) & + q_vap(ic) * dtv_dqv(ic) & - - ( qc_tot_loc(ic,i_mph) - qc_tot(ic) ) * dtv_dqc(ic) ) & + - qc_excess(ic,i_mph) * dtv_dqc(ic) ) & / ( dtv_dqsat_liq(ic) + dtv_dqv(ic) ) end do @@ -368,8 +377,8 @@ subroutine calc_env_region_tq_nb( n_points, n_points_super, & ic = index_ic(ic2) ! Find local total-condensed-water in the no-liquid cloud region qc_tot_noliq = qc_tot(ic) & - - ( frac_r(ic,i_liq) * ( qc_tot_loc(ic,i_liq) - qc_tot(ic) ) & - + frac_r(ic,i_mph) * ( qc_tot_loc(ic,i_mph) - qc_tot(ic) ) ) & + - ( frac_r(ic,i_liq) * qc_excess(ic,i_liq) & + + frac_r(ic,i_mph) * qc_excess(ic,i_mph) ) & / ( one - cloudfracs(ic,i_frac_liq) ) ! Set T for specified Tv temperature_noliq(ic) = temperature(ic) & @@ -408,7 +417,7 @@ subroutine calc_env_region_tq_nb( n_points, n_points_super, & q_vap_r(ic,i_icr) = & ( qsat_liq(ic) * dtv_dqsat_liq(ic) & + q_vap(ic) * dtv_dqv(ic) & - - ( qc_tot_loc(ic,i_icr) - qc_tot(ic) ) * dtv_dqc(ic) ) & + - qc_excess(ic,i_icr) * dtv_dqc(ic) ) & / ( dtv_dqsat_liq(ic) + dtv_dqv(ic) ) ! If below freezing and ice is present: @@ -426,7 +435,7 @@ subroutine calc_env_region_tq_nb( n_points, n_points_super, & + ice_frac_icr * ( & ( qsat_ice(ic) * dtv_dqsat_ice(ic) & + q_vap(ic) * dtv_dqv(ic) & - - ( qc_tot_loc(ic,i_icr) - qc_tot(ic) ) * dtv_dqc(ic) ) & + - qc_excess(ic,i_icr) * dtv_dqc(ic) ) & / ( dtv_dqsat_ice(ic) + dtv_dqv(ic) ) ) end if @@ -481,14 +490,14 @@ subroutine calc_env_region_tq_nb( n_points, n_points_super, & ! Upper bound is liquid saturation q_lim = ( ( qsat_liq(ic) + supersat(ic) ) * dtv_dqsat_liq(ic) & + q_vap(ic) * dtv_dqv(ic) & - - ( qc_tot_loc(ic,i_icr) - qc_tot(ic) ) * dtv_dqc(ic) ) & + - qc_excess(ic,i_icr) * dtv_dqc(ic) ) & / ( dtv_dqsat_liq(ic) + dtv_dqv(ic) ) q_vap_r(ic,i_icr) = min( q_vap_r(ic,i_icr), q_lim ) ! Upper bound for dry region is also liquid saturation q_lim = ( ( qsat_liq(ic) + supersat(ic) ) * dtv_dqsat_liq(ic) & + q_vap(ic) * dtv_dqv(ic) & - - ( qc_tot_loc(ic,i_dry) - qc_tot(ic) ) * dtv_dqc(ic) ) & + - qc_excess(ic,i_dry) * dtv_dqc(ic) ) & / ( dtv_dqsat_liq(ic) + dtv_dqv(ic) ) ! This translates into a lower bound for the icr region: ! f_dry qv_dry + f_icr qv_icr = (f_dry + f_icr) qv_noliq @@ -569,7 +578,7 @@ subroutine calc_env_region_tq_nb( n_points, n_points_super, & ! Calculate T so-as to yield the correct buoyancy temperature_r(ic,i_region) = temperature(ic) & - ( ( q_vap_r(ic,i_region) - q_vap(ic) ) * dtv_dqv(ic) & - + ( qc_tot_loc(ic,i_region) - qc_tot(ic) ) * dtv_dqc(ic) ) / dtv_dt(ic) + + qc_excess(ic,i_region) * dtv_dqc(ic) ) / dtv_dt(ic) end do end if diff --git a/science/physics_schemes/source/convection/comorph/plume_genesis/calc_env_regions.F90 b/science/physics_schemes/source/convection/comorph/plume_genesis/calc_env_regions.F90 index 1a7aa0537d..f0fad599e7 100644 --- a/science/physics_schemes/source/convection/comorph/plume_genesis/calc_env_regions.F90 +++ b/science/physics_schemes/source/convection/comorph/plume_genesis/calc_env_regions.F90 @@ -56,7 +56,7 @@ subroutine calc_env_regions( n_points, n_points_super, & i_frac_bulk, i_frac_precip use cmpr_type_mod, only: cmpr_type use check_bad_values_mod, only: check_bad_values_cmpr -use fields_type_mod, only: field_names, n_fields, & +use fields_type_mod, only: field_names, field_min, field_max, n_fields, & i_temperature, i_q_vap, i_qc_first use genesis_diags_type_mod, only: genesis_diags_type @@ -165,7 +165,6 @@ subroutine calc_env_regions( n_points, n_points_super, & ! Description of where we are in the code, for error messages character(len=name_length) :: where_string character(len=name_length) :: field_name -logical :: field_positive ! Loop counter integer :: ic, ic2, i_cond, i_region, i_field, i_super @@ -446,36 +445,34 @@ subroutine calc_env_regions( n_points, n_points_super, & if ( i_check_bad_values_cmpr > i_check_bad_none ) then where_string = "calc_env_regions call for " // trim(adjustl(call_string)) - field_positive = .true. do i_region = 1, n_regions ! Check area fractions field_name = trim(adjustl(region_names(i_region))) // "_frac" - call check_bad_values_cmpr( cmpr, k, & - frac_r(:,i_region), & + call check_bad_values_cmpr( cmpr, k, frac_r(:,i_region), & where_string, field_name, & - field_positive ) + field_min=zero, field_max=one ) ! Check temperatures field_name = trim(adjustl(region_names(i_region))) // "_temperature" - call check_bad_values_cmpr( cmpr, k, & - temperature_r(:,i_region), & + call check_bad_values_cmpr( cmpr, k, temperature_r(:,i_region), & where_string, field_name, & - field_positive ) + field_min=field_min(i_temperature), & + field_max=field_max(i_temperature) ) ! Check water-vapour mixing-ratio field_name = trim(adjustl(region_names(i_region))) // "_q_vap" - call check_bad_values_cmpr( cmpr, k, & - q_vap_r(:,i_region), & + call check_bad_values_cmpr( cmpr, k, q_vap_r(:,i_region), & where_string, field_name, & - field_positive ) + field_min=field_min(i_q_vap), & + field_max=field_max(i_q_vap) ) end do ! Check condensed water species do i_cond = 1, n_cond_species field_name = "loc_" // trim(adjustl(field_names(i_qc_first-1+i_cond))) - call check_bad_values_cmpr( cmpr, k, & - q_cond_loc(:,i_cond), & + call check_bad_values_cmpr( cmpr, k, q_cond_loc(:,i_cond), & where_string, field_name, & - field_positive ) + field_min=field_min(i_qc_first-1+i_cond), & + field_max=field_max(i_qc_first-1+i_cond) ) end do end if ! ( i_check_bad_values_cmpr > i_check_bad_none ) diff --git a/science/physics_schemes/source/convection/comorph/plume_genesis/calc_turb_parcel.F90 b/science/physics_schemes/source/convection/comorph/plume_genesis/calc_turb_parcel.F90 index c6b8807923..e46a60ca37 100644 --- a/science/physics_schemes/source/convection/comorph/plume_genesis/calc_turb_parcel.F90 +++ b/science/physics_schemes/source/convection/comorph/plume_genesis/calc_turb_parcel.F90 @@ -113,7 +113,6 @@ subroutine calc_turb_parcel( n_points, n_points_super, cmpr_init, k, & ! String containing info to print in error messages character(len=name_length) :: call_string character(len=name_length) :: field_name -logical :: l_positive ! Loop counters integer :: ic, i_field @@ -253,15 +252,9 @@ subroutine calc_turb_parcel( n_points, n_points_super, cmpr_init, k, & ! Parcel initial perturbations do i_field = i_wind_u, i_q_vap - ! Perturbations can be +ive or -ive, except w' which must be positive - if ( i_field==i_wind_w) then - l_positive = .true. - else - l_positive = .false. - end if field_name = "turb_pert_k_" // trim(adjustl(field_names(i_field))) call check_bad_values_cmpr( cmpr_init, k, turb_pert_k(:,i_field), & - call_string, field_name, l_positive ) + call_string, field_name ) end do end if ! ( i_check_bad_values_cmpr > i_check_bad_none ) diff --git a/science/physics_schemes/source/convection/comorph/plume_genesis/init_mass_moist_frac.F90 b/science/physics_schemes/source/convection/comorph/plume_genesis/init_mass_moist_frac.F90 index 20acede8b4..21c82d7c0f 100644 --- a/science/physics_schemes/source/convection/comorph/plume_genesis/init_mass_moist_frac.F90 +++ b/science/physics_schemes/source/convection/comorph/plume_genesis/init_mass_moist_frac.F90 @@ -41,10 +41,9 @@ subroutine init_mass_moist_frac( n_points, n_points_super, & turb_kmh, turb_kph, & grid_km1, grid_kmh, grid_k, & grid_kph, grid_kp1, & - fields_km1, fields_k, & - fields_kp1, cloudfracs_k, & - virt_temp_km1, virt_temp_k, & - virt_temp_kp1, & + fields_km1, fields_k, fields_kp1, & + cloudfracs_k, & + virt_temp_km1, virt_temp_k, virt_temp_kp1, & updraft_par_gen, dndraft_par_gen, & genesis_diags, diags_super ) @@ -59,10 +58,10 @@ subroutine init_mass_moist_frac( n_points, n_points_super, & use grid_type_mod, only: n_grid, i_pressure use parcel_type_mod, only: parcel_type, i_massflux_d, i_radius, & - parcel_check_bad_values + par_min, par_max, parcel_check_bad_values use fields_type_mod, only: n_fields, i_wind_u, i_temperature, & i_q_vap, i_qc_first, i_qc_last, i_q_cl, & - field_names, field_positive + field_names, field_min, field_max use turb_type_mod, only: n_turb use cloudfracs_type_mod, only: n_cloudfracs use cmpr_type_mod, only: cmpr_type, cmpr_alloc, cmpr_dealloc @@ -265,7 +264,6 @@ subroutine init_mass_moist_frac( n_points, n_points_super, & ! Character string for error messages character(len=name_length) :: call_string type(cmpr_type) :: cmpr_check -logical, parameter :: l_positive = .true. character(len=name_length) :: field_name ! Loop counters @@ -281,15 +279,18 @@ subroutine init_mass_moist_frac( n_points, n_points_super, & call_string = "Start of init_mass_moist_frac, fields_km1" call check_bad_values_cmpr( cmpr_init, k, fields_km1(:,i_field), & call_string, field_names(i_field), & - field_positive(i_field) ) + field_min=field_min(i_field), & + field_max=field_max(i_field) ) call_string = "Start of init_mass_moist_frac, fields_k" call check_bad_values_cmpr( cmpr_init, k, fields_k(:,i_field), & call_string, field_names(i_field), & - field_positive(i_field) ) + field_min=field_min(i_field), & + field_max=field_max(i_field) ) call_string = "Start of init_mass_moist_frac, fields_kp1" call check_bad_values_cmpr( cmpr_init, k, fields_kp1(:,i_field), & call_string, field_names(i_field), & - field_positive(i_field) ) + field_min=field_min(i_field), & + field_max=field_max(i_field) ) end do end if @@ -572,12 +573,15 @@ subroutine init_mass_moist_frac( n_points, n_points_super, & do i_type = 1, n_updraft_types write(field_name,"(A,I3)") "massflux_d type", i_type call check_bad_values_cmpr( cmpr_check, k, init_mass_up_t(:,i_type), & - call_string, field_name, l_positive ) + call_string, field_name, & + field_min=par_min(i_massflux_d), & + field_max=par_max(i_massflux_d) ) end do do i_field = i_temperature, n_fields call check_bad_values_cmpr( cmpr_check, k, fields_par_up(:,i_field), & call_string, field_names(i_field), & - field_positive(i_field) ) + field_min=field_min(i_field), & + field_max=field_max(i_field) ) end do do i_type = 1, n_updraft_types do ic2 = 1, nc_up @@ -588,10 +592,14 @@ subroutine init_mass_moist_frac( n_points, n_points_super, & end do write(field_name,"(A,I3)") "temperature + tl_pert conv type ", i_type call check_bad_values_cmpr( cmpr_check, k, work1, & - call_string, field_name, l_positive ) + call_string, field_name, & + field_min=field_min(i_temperature), & + field_max=field_max(i_temperature) ) write(field_name,"(A,I3)") "q_vap + qt_pert conv type ", i_type call check_bad_values_cmpr( cmpr_check, k, work2, & - call_string, field_name, l_positive ) + call_string, field_name, & + field_min=field_min(i_q_vap), & + field_max=field_max(i_q_vap) ) end do call cmpr_dealloc( cmpr_check ) end if ! ( nc_up > 0 ) @@ -609,12 +617,15 @@ subroutine init_mass_moist_frac( n_points, n_points_super, & do i_type = 1, n_dndraft_types write(field_name,"(A,I3)") "massflux_d type", i_type call check_bad_values_cmpr( cmpr_check, k, init_mass_dn_t(:,i_type), & - call_string, field_name, l_positive ) + call_string, field_name, & + field_min=par_min(i_massflux_d), & + field_max=par_max(i_massflux_d) ) end do do i_field = i_temperature, n_fields call check_bad_values_cmpr( cmpr_check, k, fields_par_dn(:,i_field), & call_string, field_names(i_field), & - field_positive(i_field) ) + field_min=field_min(i_field), & + field_max=field_max(i_field) ) end do do i_type = 1, n_dndraft_types do ic2 = 1, nc_dn @@ -625,10 +636,14 @@ subroutine init_mass_moist_frac( n_points, n_points_super, & end do write(field_name,"(A,I3)") "temperature + tl_pert conv type ", i_type call check_bad_values_cmpr( cmpr_check, k, work1, & - call_string, field_name, l_positive ) + call_string, field_name, & + field_min=field_min(i_temperature), & + field_max=field_max(i_temperature) ) write(field_name,"(A,I3)") "q_vap + qt_pert conv type ", i_type call check_bad_values_cmpr( cmpr_check, k, work2, & - call_string, field_name, l_positive ) + call_string, field_name, & + field_min=field_min(i_q_vap), & + field_max=field_max(i_q_vap) ) end do call cmpr_dealloc( cmpr_check ) end if ! ( nc_dn > 0 ) diff --git a/science/physics_schemes/source/convection/comorph/plume_genesis/region_parcel_calcs.F90 b/science/physics_schemes/source/convection/comorph/plume_genesis/region_parcel_calcs.F90 index b9bd71b3e0..b650f72ad4 100644 --- a/science/physics_schemes/source/convection/comorph/plume_genesis/region_parcel_calcs.F90 +++ b/science/physics_schemes/source/convection/comorph/plume_genesis/region_parcel_calcs.F90 @@ -43,7 +43,9 @@ subroutine region_parcel_calcs( n_points, n_points_super, n_conv_types, & use cmpr_type_mod, only: cmpr_type, cmpr_alloc, cmpr_dealloc use fields_type_mod, only: n_fields, i_wind_u, i_wind_w, & i_temperature, i_q_vap, i_q_cl, & - i_qc_first, i_qc_last, field_names, field_positive + i_qc_first, i_qc_last, field_names, & + field_min, field_max +use parcel_type_mod, only: par_min, par_max, i_massflux_d use grid_type_mod, only: n_grid, i_height, i_pressure use subregion_mod, only: n_regions, region_names use cloudfracs_type_mod, only: n_cloudfracs, i_frac_liq, i_frac_ice @@ -221,7 +223,6 @@ subroutine region_parcel_calcs( n_points, n_points_super, n_conv_types, & ! Structure storing i,j indices of points currently being worked ! on (for error reporting) type(cmpr_type) :: cmpr -logical :: l_positive character(len=name_length) :: field_name ! Strings used for error reporting @@ -525,25 +526,26 @@ subroutine region_parcel_calcs( n_points, n_points_super, n_conv_types, & call_string = "On output from region_parcel_calcs, region: " // & trim(adjustl(region_names(i_region))) // " " // & trim(adjustl(draft_string)) - l_positive = .true. do i_type = 1, n_conv_types write(field_name,"(A,I3)") "massflux_d conv type ", i_type call check_bad_values_cmpr( cmpr, k, init_mass_t(:,i_type), & - call_string, field_name, l_positive ) + call_string, field_name, & + field_min=par_min(i_massflux_d), & + field_max=par_max(i_massflux_d) ) end do do i_field = i_temperature, n_fields call check_bad_values_cmpr( cmpr, k, fields_par(:,i_field), & call_string, field_names(i_field), & - field_positive(i_field) ) + field_min=field_min(i_field), & + field_max=field_max(i_field) ) end do - l_positive = .false. do i_type = 1, n_conv_types write(field_name,"(A,I3)") "pert_tl conv type ", i_type call check_bad_values_cmpr( cmpr, k, pert_tl_t(:,i_type), & - call_string, field_name, l_positive ) + call_string, field_name ) write(field_name,"(A,I3)") "pert_qt conv type ", i_type call check_bad_values_cmpr( cmpr, k, pert_qt_t(:,i_type), & - call_string, field_name, l_positive ) + call_string, field_name ) end do end if ! ( i_check_bad_values_cmpr > i_check_bad_none ) diff --git a/science/physics_schemes/source/convection/comorph/plume_genesis/set_par_fields.F90 b/science/physics_schemes/source/convection/comorph/plume_genesis/set_par_fields.F90 index 2e0063ef85..6ec4a689f2 100644 --- a/science/physics_schemes/source/convection/comorph/plume_genesis/set_par_fields.F90 +++ b/science/physics_schemes/source/convection/comorph/plume_genesis/set_par_fields.F90 @@ -25,8 +25,8 @@ subroutine set_par_fields( n_points, n_points_super, n_fields_tot, & par_gen_par, par_gen_mean, par_gen_core, & rhpert_t, frac_r_t ) -use comorph_constants_mod, only: real_cvprec, zero, one, l_turb_par_gen, & - par_gen_radius, par_gen_radius_fac, & +use comorph_constants_mod, only: real_cvprec, zero, one, & + par_gen_radius_fac, & ass_min_radius, min_radius_fac, & par_gen_core_fac, & n_tracers, name_length, l_par_core, & @@ -35,7 +35,7 @@ subroutine set_par_fields( n_points, n_points_super, n_fields_tot, & use grid_type_mod, only: n_grid, i_height use fields_type_mod, only: i_wind_u, i_wind_w, i_tracers, & i_temperature, i_q_vap, i_qc_first, i_qc_last, & - field_positive, field_names + field_min, field_max, field_names use subregion_mod, only: n_regions use parcel_type_mod, only: n_par, i_radius, i_edge_virt_temp use cmpr_type_mod, only: cmpr_type @@ -118,29 +118,18 @@ subroutine set_par_fields( n_points, n_points_super, n_fields_tot, & integer :: ic, i_field, i_region -if ( l_turb_par_gen ) then - - ! Use max of turbulence-based radius and an arbitrary linear - ! ramp from the surface - do ic = 1, n_points - par_gen_par(ic,i_radius) = max( par_gen_radius_fac * turb_len_k(ic), & - min( min_radius_fac * grid_k(ic,i_height), & - ass_min_radius ) ) - end do - - ! Amplify the parcel radius using input variable scaling factor... - do ic = 1, n_points - par_gen_par(ic,i_radius) = par_gen_par(ic,i_radius) * par_radius_amp(ic) - end do - -else - - ! Set arbitrary fixed parcel radius - do ic = 1, n_points - par_gen_par(ic,i_radius) = par_gen_radius - end do +! Use max of turbulence-based radius and an arbitrary linear +! ramp from the surface +do ic = 1, n_points + par_gen_par(ic,i_radius) = max( par_gen_radius_fac * turb_len_k(ic), & + min( min_radius_fac * grid_k(ic,i_height), & + ass_min_radius ) ) +end do -end if +! Amplify the parcel radius using input variable scaling factor... +do ic = 1, n_points + par_gen_par(ic,i_radius) = par_gen_par(ic,i_radius) * par_radius_amp(ic) +end do ! Set environment virtual temperature stored in the parcel call calc_virt_temp( n_points, n_points_super, & @@ -232,12 +221,14 @@ subroutine set_par_fields( n_points, n_points_super, n_fields_tot, & field_name = "par_gen_mean_" // trim(adjustl(field_names(i_field))) call check_bad_values_cmpr( cmpr_init, k, par_gen_mean(:,i_field), & call_string, field_name, & - field_positive(i_field) ) + field_min=field_min(i_field), & + field_max=field_max(i_field) ) if ( l_par_core ) then field_name = "par_gen_core_" // trim(adjustl(field_names(i_field))) call check_bad_values_cmpr( cmpr_init, k, par_gen_core(:,i_field), & call_string, field_name, & - field_positive(i_field) ) + field_min=field_min(i_field), & + field_max=field_max(i_field) ) end if end do if ( l_tracer .and. n_tracers > 0 ) then @@ -245,12 +236,14 @@ subroutine set_par_fields( n_points, n_points_super, n_fields_tot, & field_name = "par_gen_mean_" // trim(adjustl(field_names(i_field))) call check_bad_values_cmpr( cmpr_init, k, par_gen_mean(:,i_field), & call_string, field_name, & - field_positive(i_field) ) + field_min=field_min(i_field), & + field_max=field_max(i_field) ) if ( l_par_core ) then field_name = "par_gen_core_" // trim(adjustl(field_names(i_field))) call check_bad_values_cmpr( cmpr_init, k, par_gen_core(:,i_field), & call_string, field_name, & - field_positive(i_field) ) + field_min=field_min(i_field), & + field_max=field_max(i_field) ) end if end do end if diff --git a/science/physics_schemes/source/convection/comorph/plume_genesis/test_unstable.F90 b/science/physics_schemes/source/convection/comorph/plume_genesis/test_unstable.F90 index 903a322e42..8f0688c0a0 100644 --- a/science/physics_schemes/source/convection/comorph/plume_genesis/test_unstable.F90 +++ b/science/physics_schemes/source/convection/comorph/plume_genesis/test_unstable.F90 @@ -80,7 +80,6 @@ subroutine test_unstable( virt_temp_1, virt_temp_2, & real(kind=real_hmprec), parameter :: zero_p = 0.0_real_hmprec real(kind=real_hmprec) :: gravity_p - ! Loop counters integer :: i, j diff --git a/science/physics_schemes/source/convection/comorph/plume_model/calc_cape.F90 b/science/physics_schemes/source/convection/comorph/plume_model/calc_cape.F90 index d97464b03a..146d53d483 100644 --- a/science/physics_schemes/source/convection/comorph/plume_model/calc_cape.F90 +++ b/science/physics_schemes/source/convection/comorph/plume_model/calc_cape.F90 @@ -124,9 +124,8 @@ subroutine calc_cape( n_points, sublevs, i_next, l_within_bl, & + sublevs(ic,j_massflux_d,i_prev+1) ) ! Contribution M^2 dz - mm_dz(ic) = half * dz & - *( sublevs(ic,j_massflux_d,i_prev) * sublevs(ic,j_massflux_d,i_prev) & - + sublevs(ic,j_massflux_d,i_prev+1)*sublevs(ic,j_massflux_d,i_prev+1) ) + mm_dz(ic) = half * dz * ( sublevs(ic,j_massflux_d,i_prev)**2 & + + sublevs(ic,j_massflux_d,i_prev+1)**2 ) end do @@ -226,8 +225,7 @@ subroutine calc_cape( n_points, sublevs, i_next, l_within_bl, & m_dz(ic) = m_dz(ic) + dz * sublevs(ic,j_massflux_d,i_lev) ! Contribution M^2 dz - mm_dz(ic) = mm_dz(ic) + dz & - * sublevs(ic,j_massflux_d,i_lev) * sublevs(ic,j_massflux_d,i_lev) + mm_dz(ic) = mm_dz(ic) + dz * sublevs(ic,j_massflux_d,i_lev)**2 end do ! i_lev = i_prev, i_next(ic) diff --git a/science/physics_schemes/source/convection/comorph/plume_model/calc_mean_q_cl_with_core.F90 b/science/physics_schemes/source/convection/comorph/plume_model/calc_mean_q_cl_with_core.F90 index 56bbed631f..451fd4c9a7 100644 --- a/science/physics_schemes/source/convection/comorph/plume_model/calc_mean_q_cl_with_core.F90 +++ b/science/physics_schemes/source/convection/comorph/plume_model/calc_mean_q_cl_with_core.F90 @@ -28,7 +28,7 @@ subroutine calc_mean_q_cl_with_core( n_points, n_points_super, n_fields_tot, & use comorph_constants_mod, only: real_cvprec, one, sqrt_min_delta, & l_cv_cloudfrac, i_mean_q_cl, & i_mean_q_cl_full -use fields_type_mod, only: i_temperature, i_q_vap, i_q_cl, i_q_cf, & +use fields_type_mod, only: i_temperature, i_q_vap, i_q_cl, & i_qc_first, i_qc_last, i_cf_liq, i_cf_bulk use set_cp_tot_mod, only: set_cp_tot @@ -125,9 +125,8 @@ subroutine calc_mean_q_cl_with_core( n_points, n_points_super, n_fields_tot, & if ( l_cv_cloudfrac ) then ! Reset the in-parcel cloud-fractions using the modified q_cl - call set_par_cloudfrac( n_points, n_points_super, & - par_mean_fields(:,i_q_cl), & - par_mean_fields(:,i_q_cf), & + call set_par_cloudfrac( n_points, n_points_super, n_points_super, & + par_mean_fields(:,i_qc_first:i_qc_last), & par_mean_fields(:,i_cf_liq:i_cf_bulk) ) end if ! ( l_cv_cloudfrac ) diff --git a/science/physics_schemes/source/convection/comorph/plume_model/conv_level_step.F90 b/science/physics_schemes/source/convection/comorph/plume_model/conv_level_step.F90 index 4bd19c9537..faab44c9ef 100644 --- a/science/physics_schemes/source/convection/comorph/plume_model/conv_level_step.F90 +++ b/science/physics_schemes/source/convection/comorph/plume_model/conv_level_step.F90 @@ -56,7 +56,6 @@ subroutine conv_level_step( & use plume_model_diags_type_mod, only: plume_model_diags_type use calc_virt_temp_mod, only: calc_virt_temp -use dry_adiabat_mod, only: dry_adiabat use calc_core_mean_ratio_mod, only: calc_core_mean_ratio use calc_env_nsq_mod, only: calc_env_nsq use set_ent_mod, only: set_ent @@ -65,7 +64,6 @@ subroutine conv_level_step( & use entdet_res_source_mod, only: entdet_res_source use init_sublevs_mod, only: init_sublevs use parcel_dyn_mod, only: parcel_dyn, i_call_mean, i_call_core, i_call_det -use set_par_cloudfrac_mod, only: set_par_cloudfrac use calc_rho_dry_mod, only: calc_rho_dry use update_edge_virt_temp_mod, only: update_edge_virt_temp use update_par_radius_mod, only: update_par_radius @@ -393,44 +391,18 @@ subroutine conv_level_step( & ! into the parcel !---------------------------------------------------------------- -! Set entrained air properties the same as the mean environment at level k -do i_field = 1, n_fields_tot - do ic = 1, n_points - ent_fields(ic,i_field) = env_k_fields(ic,i_field) - end do -end do - -if ( l_to_full_level ) then - ! If this is the first of the two half-level-steps - ! (from previous half-level to level k), - ! then the environment fields to entrain are at level k, but the - ! parcel is at a different pressure, at the previous half-level. - ! Adjust the environment temperature to what it would be at the - ! start of the level-step, so that we entrain it into the parcel - ! consistently... - do ic = 1, n_points - exner_ratio(ic) = one - end do - call dry_adiabat( n_points, n_points, & - grid_next_super(:,i_pressure), grid_prev_super(:,i_pressure),& - ent_fields(:,i_q_vap), & - ent_fields(:,i_qc_first:i_qc_last), & - exner_ratio ) - do ic = 1, n_points - ent_fields(ic,i_temperature) = ent_fields(ic,i_temperature) & - * exner_ratio(ic) - end do -end if - -! Set amount of entrained dry-mass over the current half-level-step +! Set amount of entrained dry-mass over the current half-level-step, +! and properties of the entrained air. call set_ent( n_points, n_fields_tot, max_points, & - max_ent_frac, & - par_conv_mean_fields, ent_fields, & + n_points_diag, n_diags_super, & + l_to_full_level, max_ent_frac, & + par_conv_mean_fields, env_k_fields, & grid_prev_super, grid_next_super, & par_conv_super, & l_within_bl, core_mean_ratio, & layer_mass_step, sum_massflux, & - ent_mass_d, core_ent_ratio ) + ent_fields, exner_ratio, ent_mass_d, core_ent_ratio, & + plume_model_diags, diags_super ) ! Add the entrained mass onto the mass-flux do ic = 1, n_points @@ -524,14 +496,6 @@ subroutine conv_level_step( & ent_mass_d, ent_fields, & res_source_super, res_source_fields ) -! Core environment entrainment ratio diagnostic -if ( plume_model_diags % core_ent_ratio % flag ) then - i_diag = plume_model_diags % core_ent_ratio % i_super - do ic = 1, n_points - diags_super(ic,i_diag) = core_ent_ratio(ic) - end do -end if - ! Save diagnostics of entrained mass and air properties ! (in conserved variable form ready for finding mean over types) if ( plume_model_diags % ent_mass_d % flag ) then diff --git a/science/physics_schemes/source/convection/comorph/plume_model/init_sublevs.F90 b/science/physics_schemes/source/convection/comorph/plume_model/init_sublevs.F90 index e8e31a8ffb..cd037312a0 100644 --- a/science/physics_schemes/source/convection/comorph/plume_model/init_sublevs.F90 +++ b/science/physics_schemes/source/convection/comorph/plume_model/init_sublevs.F90 @@ -161,7 +161,7 @@ subroutine init_sublevs( n_points, n_points_super, l_down, & end do -! For safety, reset delta_tv_next to zero in statically-unstable +! For safety, reset delta_tv to zero in statically-unstable ! layers so we just use explicit value of env Tv for detrainment. ! Expected sign depends on whether this is updraft or downdraft if ( l_down ) then diff --git a/science/physics_schemes/source/convection/comorph/plume_model/parcel_dyn.F90 b/science/physics_schemes/source/convection/comorph/plume_model/parcel_dyn.F90 index b80acdc7df..e100b96bf0 100644 --- a/science/physics_schemes/source/convection/comorph/plume_model/parcel_dyn.F90 +++ b/science/physics_schemes/source/convection/comorph/plume_model/parcel_dyn.F90 @@ -66,7 +66,7 @@ subroutine parcel_dyn( n_points, n_points_prev, n_points_next, & i_qc_first, i_qc_last, & i_cf_liq, i_cf_bulk, & i_wind_u, i_wind_v, i_wind_w, i_tracers, & - n_fields, field_names, field_positive + n_fields, field_names, field_min, field_max use cmpr_type_mod, only: cmpr_type use linear_qs_mod, only: linear_qs_set_ref, & @@ -290,9 +290,6 @@ subroutine parcel_dyn( n_points, n_points_prev, n_points_next, & ! q_vap and q_cl when calculating a saturated reference T logical, parameter :: l_update_q_false = .false. -! Flag passed into check_bad_values to indicate whether fields must be positive -logical :: l_positive - ! Flag for parcel mean ascent with an accompanying core logical :: l_mean_with_core @@ -364,24 +361,19 @@ subroutine parcel_dyn( n_points, n_points_prev, n_points_next, & trim(adjustl(call_string)) // "; " // & "par_next_fields" do i_field = 1, n_fields_tot - call check_bad_values_cmpr( cmpr, k, & - par_next_fields(:,i_field), & - where_string, & - field_names(i_field), & - field_positive(i_field) ) + call check_bad_values_cmpr( cmpr, k, par_next_fields(:,i_field), & + where_string, field_names(i_field), & + field_min = field_min(i_field), & + field_max = field_max(i_field) ) end do if ( present( res_source_fields ) ) then ! Check resolved-scale source-terms where_string = "Start of parcel_dyn call for " // & trim(adjustl(call_string)) // "; " // & "res_source_fields" - l_positive = .false. ! Source terms may be positive or negative do i_field = 1, n_fields_tot - call check_bad_values_cmpr( cmpr, k, & - res_source_fields(:,i_field), & - where_string, & - field_names(i_field), & - l_positive ) + call check_bad_values_cmpr( cmpr, k, res_source_fields(:,i_field), & + where_string, field_names(i_field) ) end do end if end if @@ -618,9 +610,8 @@ subroutine parcel_dyn( n_points, n_points_prev, n_points_next, & ! Update in-parcel cloud-fractions and precipitation fraction if ( l_cv_cloudfrac ) then - call set_par_cloudfrac( n_points, n_points_next, & - par_next_fields(:,i_q_cl), & - par_next_fields(:,i_q_cf), & + call set_par_cloudfrac( n_points, n_points_next, n_points_next, & + par_next_fields(:,i_qc_first:i_qc_last), & par_next_fields(:,i_cf_liq:i_cf_bulk) ) end if @@ -773,24 +764,19 @@ subroutine parcel_dyn( n_points, n_points_prev, n_points_next, & trim(adjustl(call_string)) // "; " // & "par_next_fields" do i_field = 1, n_fields_tot - call check_bad_values_cmpr( cmpr, k, & - par_next_fields(:,i_field), & - where_string, & - field_names(i_field), & - field_positive(i_field) ) + call check_bad_values_cmpr( cmpr, k, par_next_fields(:,i_field), & + where_string, field_names(i_field), & + field_min = field_min(i_field), & + field_max = field_max(i_field) ) end do if ( present( res_source_fields ) ) then ! Check resolved-scale source-terms where_string = "End of parcel_dyn call for " // & trim(adjustl(call_string)) // "; " // & "res_source_fields" - l_positive = .false. ! Source terms may be positive or negative do i_field = 1, n_fields_tot - call check_bad_values_cmpr( cmpr, k, & - res_source_fields(:,i_field), & - where_string, & - field_names(i_field), & - l_positive ) + call check_bad_values_cmpr( cmpr, k, res_source_fields(:,i_field), & + where_string, field_names(i_field) ) end do end if end if diff --git a/science/physics_schemes/source/convection/comorph/plume_model/set_det.F90 b/science/physics_schemes/source/convection/comorph/plume_model/set_det.F90 index 64b30605b8..2d36dc624d 100644 --- a/science/physics_schemes/source/convection/comorph/plume_model/set_det.F90 +++ b/science/physics_schemes/source/convection/comorph/plume_model/set_det.F90 @@ -36,12 +36,13 @@ subroutine set_det( n_points, max_points, n_points_res, n_fields_tot, & use cmpr_type_mod, only: cmpr_type use fields_type_mod, only: field_positive, field_names, i_wind_w, & + field_min_cons, field_max_cons, & i_qc_first, i_qc_last use sublevs_mod, only: max_sublevs, n_sublev_vars, i_prev, & j_height, j_mean_buoy, j_core_buoy, j_delta_tv, & j_massflux_d, j_env_tv, j_env_w, & j_mean_wex, j_core_wex -use parcel_type_mod, only: n_par, i_massflux_d +use parcel_type_mod, only: n_par, i_massflux_d, par_min, par_max use solve_detrainment_mod, only: solve_detrainment use wind_w_eqn_mod, only: wind_w_eqn @@ -186,8 +187,6 @@ subroutine set_det( n_points, max_points, n_points_res, n_fields_tot, & character(len=name_length) :: where_string ! Field name for error messages character(len=name_length) :: field_name -! Flag for checking positivity (input to bad value checker) -logical :: l_positive ! Numerical tolerance used in safety-check to avoid negative values real(kind=real_cvprec) :: tolerance @@ -211,7 +210,8 @@ subroutine set_det( n_points, max_points, n_points_res, n_fields_tot, & do i_field = 1, n_fields_tot call check_bad_values_cmpr( cmpr, k, par_next_mean_fields(:,i_field), & where_string, field_names(i_field), & - field_positive(i_field) ) + field_min = field_min_cons(i_field), & + field_max = field_max_cons(i_field) ) end do end if @@ -259,17 +259,16 @@ subroutine set_det( n_points, max_points, n_points_res, n_fields_tot, & do i_field = 1, n_fields_tot call check_bad_values_cmpr( cmpr, k, par_next_core_fields(:,i_field), & where_string, field_names(i_field), & - field_positive(i_field) ) + field_min = field_min_cons(i_field), & + field_max = field_max_cons(i_field) ) end do where_string = "On input to set_det call for " // & trim(adjustl(draft_string)) field_name = "core_mean_ratio" - l_positive = .true. call check_bad_values_cmpr( cmpr, k, core_mean_ratio, & - where_string, & - field_name, l_positive ) + where_string, field_name, field_min=zero ) end if @@ -480,14 +479,12 @@ subroutine set_det( n_points, max_points, n_points_res, n_fields_tot, & trim(adjustl(draft_string)) field_name = "x_edge" - l_positive = .true. - call check_bad_values_cmpr( cmpr, k, x_edge, where_string, & - field_name, l_positive ) + call check_bad_values_cmpr( cmpr, k, x_edge, where_string, field_name, & + field_min=zero, field_max=one ) field_name = "frac (non-detrained fraction)" - l_positive = .true. - call check_bad_values_cmpr( cmpr, k, frac, where_string, & - field_name, l_positive ) + call check_bad_values_cmpr( cmpr, k, frac, where_string, field_name, & + field_min=zero, field_max=one ) end if @@ -899,7 +896,8 @@ subroutine set_det( n_points, max_points, n_points_res, n_fields_tot, & do i_field = 1, n_fields_tot call check_bad_values_cmpr( cmpr, k, par_next_mean_fields(:,i_field), & where_string, field_names(i_field), & - field_positive(i_field) ) + field_min = field_min_cons(i_field), & + field_max = field_max_cons(i_field) ) end do ! Check detrained properties @@ -909,20 +907,22 @@ subroutine set_det( n_points, max_points, n_points_res, n_fields_tot, & do i_field = 1, n_fields_tot call check_bad_values_cmpr( cmpr, k, det_fields(:,i_field), & where_string, field_names(i_field), & - field_positive(i_field) ) + field_min = field_min_cons(i_field), & + field_max = field_max_cons(i_field) ) end do ! Check detrained mass field_name = "det_mass_d" - l_positive = .true. - call check_bad_values_cmpr( cmpr, k, det_mass_d, where_string, & - field_name, l_positive ) + call check_bad_values_cmpr( cmpr, k, det_mass_d, & + where_string, field_name, & + field_min = par_min(i_massflux_d), & + field_max = par_max(i_massflux_d) ) ! Check mass-flux field_name = "next_massflux_d" - l_positive = .true. call check_bad_values_cmpr( cmpr, k, par_conv_super(:,i_massflux_d), & - where_string, & - field_name, l_positive ) + where_string, field_name, & + field_min = par_min(i_massflux_d), & + field_max = par_max(i_massflux_d) ) end if diff --git a/science/physics_schemes/source/convection/comorph/plume_model/set_diag_conv_cloud_a.F90 b/science/physics_schemes/source/convection/comorph/plume_model/set_diag_conv_cloud_a.F90 index 612e54cba1..fbbdfde5b4 100644 --- a/science/physics_schemes/source/convection/comorph/plume_model/set_diag_conv_cloud_a.F90 +++ b/science/physics_schemes/source/convection/comorph/plume_model/set_diag_conv_cloud_a.F90 @@ -28,9 +28,11 @@ subroutine set_diag_conv_cloud_a( n_points, n_points_env, n_points_res, & use comorph_constants_mod, only: real_cvprec, zero, half, & gravity, cf_conv_fac, & + n_cond_species, cond_params, & + i_sg_frac_liq, i_sg_frac_ice, & wind_w_buoy_fac, w_min, l_cv_cloudfrac -use fields_type_mod, only: n_fields, i_temperature, & - i_q_vap, i_q_cl, i_q_cf, i_cf_liq, i_cf_bulk +use fields_type_mod, only: n_fields, i_temperature, i_q_vap, & + i_qc_first, i_qc_last, i_cf_liq, i_cf_bulk use cloudfracs_type_mod, only: n_convcloud, i_frac_liq, i_frac_bulk use sublevs_mod, only: max_sublevs, n_sublev_vars, i_prev, & j_height, j_mean_buoy @@ -130,7 +132,7 @@ subroutine set_diag_conv_cloud_a( n_points, n_points_env, n_points_res, & real(kind=real_cvprec) :: convcloud_step ( n_points, n_convcloud ) ! Loop counter -integer :: ic, i_lev, i_field +integer :: ic, i_lev, i_field, i_cond ! Calculate the convective fraction irrespective of whether @@ -195,12 +197,34 @@ subroutine set_diag_conv_cloud_a( n_points, n_points_env, n_points_res, & * next_radius(ic), w_min*w_min ) ) ) end do -! Compute grid-mean convective liquid and ice mixing ratios +! Sum total parcel liquid and ice cloud mixing-ratios over species do ic = 1, n_points - prev_q_cl_conv(ic) = prev_cf_conv(ic)*par_prev_mean(ic,i_q_cl) - prev_q_cf_conv(ic) = prev_cf_conv(ic)*par_prev_mean(ic,i_q_cf) - next_q_cl_conv(ic) = next_cf_conv(ic)*par_next_mean(ic,i_q_cl) - next_q_cf_conv(ic) = next_cf_conv(ic)*par_next_mean(ic,i_q_cf) + prev_q_cl_conv(ic) = zero + prev_q_cf_conv(ic) = zero + next_q_cl_conv(ic) = zero + next_q_cf_conv(ic) = zero +end do +do i_cond = 1, n_cond_species + i_field = i_qc_first - 1 + i_cond + if ( cond_params(i_cond)%pt % i_sg == i_sg_frac_liq ) then + do ic = 1, n_points + prev_q_cl_conv(ic) = prev_q_cl_conv(ic) + par_prev_mean(ic,i_field) + next_q_cl_conv(ic) = next_q_cl_conv(ic) + par_next_mean(ic,i_field) + end do + end if + if ( cond_params(i_cond)%pt % i_sg == i_sg_frac_ice ) then + do ic = 1, n_points + prev_q_cf_conv(ic) = prev_q_cf_conv(ic) + par_prev_mean(ic,i_field) + next_q_cf_conv(ic) = next_q_cf_conv(ic) + par_next_mean(ic,i_field) + end do + end if +end do +! Scale by convective fraction to convert to grid-means +do ic = 1, n_points + prev_q_cl_conv(ic) = prev_q_cl_conv(ic) * prev_cf_conv(ic) + prev_q_cf_conv(ic) = prev_q_cf_conv(ic) * prev_cf_conv(ic) + next_q_cl_conv(ic) = next_q_cl_conv(ic) * next_cf_conv(ic) + next_q_cf_conv(ic) = next_q_cf_conv(ic) * next_cf_conv(ic) end do ! Inflate the fraction without inflating the grid-mean @@ -256,11 +280,11 @@ subroutine set_diag_conv_cloud_a( n_points, n_points_env, n_points_res, & allocate( par_next_cloudfracs( n_points, i_frac_liq:i_frac_bulk ) ) ! Call routine to diagnose in-parcel cloud-fractions based on q_cl, q_cf - call set_par_cloudfrac( n_points, n_points, & - par_prev_mean(:,i_q_cl), par_prev_mean(:,i_q_cf), & + call set_par_cloudfrac( n_points, n_points, n_points, & + par_prev_mean(:,i_qc_first:i_qc_last), & par_prev_cloudfracs ) - call set_par_cloudfrac( n_points, n_points, & - par_next_mean(:,i_q_cl), par_next_mean(:,i_q_cf), & + call set_par_cloudfrac( n_points, n_points, n_points, & + par_next_mean(:,i_qc_first:i_qc_last), & par_next_cloudfracs ) ! Use diagnosed in-parcel cloud-fractions to interpolate to full-level diff --git a/science/physics_schemes/source/convection/comorph/plume_model/set_ent.F90 b/science/physics_schemes/source/convection/comorph/plume_model/set_ent.F90 index c8f7225ae7..42f2c38b3b 100644 --- a/science/physics_schemes/source/convection/comorph/plume_model/set_ent.F90 +++ b/science/physics_schemes/source/convection/comorph/plume_model/set_ent.F90 @@ -16,22 +16,26 @@ module set_ent_mod ! Subroutine sets the entrained mass from the current layer, ! and sets the properties of the entrained air subroutine set_ent( n_points, n_fields_tot, max_points, & - max_ent_frac, & - par_conv_mean_fields, ent_fields, & + n_points_diag, n_diags_super, & + l_to_full_level, max_ent_frac, & + par_conv_mean_fields, env_k_fields, & grid_prev_super, grid_next_super, & par_conv_super, & l_within_bl, core_mean_ratio, & layer_mass_step, sum_massflux, & - ent_mass_d, core_ent_ratio ) + ent_fields, exner_ratio, ent_mass_d, core_ent_ratio, & + plume_model_diags, diags_super ) use comorph_constants_mod, only: real_cvprec, min_float, one, & ent_coef, comorph_timestep, & core_ent_fac, l_core_ent_cmr, & i_cfl_local, i_cfl_local_all, & i_cfl_local_nobl -use fields_type_mod, only: i_temperature, i_q_vap +use fields_type_mod, only: i_temperature, i_q_vap, i_qc_first, i_qc_last use grid_type_mod, only: n_grid, i_height, i_pressure use parcel_type_mod, only: n_par, i_massflux_d, i_radius +use plume_model_diags_type_mod, only: plume_model_diags_type +use dry_adiabat_mod, only: dry_adiabat use calc_rho_dry_mod, only: calc_rho_dry implicit none @@ -48,15 +52,23 @@ subroutine set_ent( n_points, n_fields_tot, max_points, & ! which will often be bigger than the number of points here) integer, intent(in) :: max_points +! Dimensions of the diagnostics super-array +integer, intent(in) :: n_points_diag +integer, intent(in) :: n_diags_super + +! Flag for first half-level ascent from half-level to full-level +logical, intent(in) :: l_to_full_level + ! Maximum allowed entrained fraction of layer mass for current draft real(kind=real_cvprec), intent(in) :: max_ent_frac ! Super-array containing the parcel mean primary fields real(kind=real_cvprec), intent(in) :: par_conv_mean_fields & ( max_points, n_fields_tot ) -! Entrained primary field values -real(kind=real_cvprec), intent(in) :: ent_fields & - ( n_points, n_fields_tot ) +! Super-array containing the environment primary fields +! at the current thermodynamic level k; +real(kind=real_cvprec), intent(in) :: env_k_fields & + ( max_points, n_fields_tot ) ! Height and pressure of parcel before and after the level-step real(kind=real_cvprec), intent(in) :: grid_prev_super ( max_points, n_grid ) @@ -77,12 +89,25 @@ subroutine set_ent( n_points, n_fields_tot, max_points, & ! Sum of previous level-interface mass-fluxes over all convection types/layers real(kind=real_cvprec), intent(in) :: sum_massflux(n_points) +! Entrained primary field values +real(kind=real_cvprec), intent(out) :: ent_fields & + ( n_points, n_fields_tot ) + +! Exner pressure factor for dry-adiabatic adjustment to prev +real(kind=real_cvprec), intent(out) :: exner_ratio(n_points) + ! Rate of entrainment of dry-mass from current level / kg m-2 s-1 real(kind=real_cvprec), intent(out) :: ent_mass_d(n_points) ! Weight used to calculate properties of air entrained into the core real(kind=real_cvprec), intent(out) :: core_ent_ratio(n_points) +! Structure storing flags and super-array addresses for +! various diagnostics +type(plume_model_diags_type), intent(in) :: plume_model_diags +! Super-array to contain output diagnostics +real(kind=real_cvprec), intent(in out) :: diags_super & + ( n_points_diag, n_diags_super ) ! Dry-density of the entrained air and the parcel real(kind=real_cvprec) :: ent_rho_dry(n_points) @@ -92,11 +117,45 @@ subroutine set_ent( n_points, n_fields_tot, max_points, & real(kind=real_cvprec) :: max_ent(n_points) ! Loop counters -integer :: ic +integer :: ic, i_field, i_diag !------------------------------------------------------------------------------ -! 1) Calculate 1/R "mixing" entrainment rate +! 1) Set properties of the entrained air +!------------------------------------------------------------------------------ + +! Set entrained air properties the same as the mean environment at level k +do i_field = 1, n_fields_tot + do ic = 1, n_points + ent_fields(ic,i_field) = env_k_fields(ic,i_field) + end do +end do + +if ( l_to_full_level ) then + ! If this is the first of the two half-level-steps + ! (from previous half-level to level k), + ! then the environment fields to entrain are at level k, but the + ! parcel is at a different pressure, at the previous half-level. + ! Adjust the temperature of the entrained air to what it would be at the + ! start of the level-step, so that we entrain it into the parcel + ! consistently... + do ic = 1, n_points + exner_ratio(ic) = one + end do + call dry_adiabat( n_points, n_points, & + grid_next_super(:,i_pressure), grid_prev_super(:,i_pressure),& + ent_fields(:,i_q_vap), & + ent_fields(:,i_qc_first:i_qc_last), & + exner_ratio ) + do ic = 1, n_points + ent_fields(ic,i_temperature) = ent_fields(ic,i_temperature) & + * exner_ratio(ic) + end do +end if + + +!------------------------------------------------------------------------------ +! 2) Calculate 1/R "mixing" entrainment rate !------------------------------------------------------------------------------ ! Set the fractional volume entrainment rate in m-1 @@ -152,7 +211,7 @@ subroutine set_ent( n_points, n_fields_tot, max_points, & !------------------------------------------------------------------------------ -! 2) Apply CFL limit to entrainment for numerical stability +! 3) Apply CFL limit to entrainment for numerical stability !------------------------------------------------------------------------------ ! Compute maximum allowed total entrainment rate @@ -193,6 +252,19 @@ subroutine set_ent( n_points, n_fields_tot, max_points, & end select +!------------------------------------------------------------------------------ +! 4) Copy diagnostics into a super-array for output +!------------------------------------------------------------------------------ + +! Core environment entrainment ratio diagnostic +if ( plume_model_diags % core_ent_ratio % flag ) then + i_diag = plume_model_diags % core_ent_ratio % i_super + do ic = 1, n_points + diags_super(ic,i_diag) = core_ent_ratio(ic) + end do +end if + + return end subroutine set_ent diff --git a/science/physics_schemes/source/convection/comorph/plume_model/set_par_cloudfrac.F90 b/science/physics_schemes/source/convection/comorph/plume_model/set_par_cloudfrac.F90 index 5e83978a2c..61e44c6e41 100644 --- a/science/physics_schemes/source/convection/comorph/plume_model/set_par_cloudfrac.F90 +++ b/science/physics_schemes/source/convection/comorph/plume_model/set_par_cloudfrac.F90 @@ -17,11 +17,13 @@ module set_par_cloudfrac_mod ! Subroutine to set the in-parcel cloud fraction fields. ! Just sets them based on the presence of non-zero liquid and ! ice cloud mixing-ratios. -subroutine set_par_cloudfrac( n_points, n_points_super, & - q_cl, q_cf, cloudfracs ) +subroutine set_par_cloudfrac( n_points, n_points_cond, n_points_cf, & + q_cond_super, cloudfracs ) use comorph_constants_mod, only: real_cvprec, zero, one, & - l_cv_cf, i_par_cloudfrac, & + n_cond_species, cond_params, & + i_sg_frac_liq, i_sg_frac_ice, & + i_par_cloudfrac, & i_par_cloudfrac_hom, i_par_cloudfrac_mph, & overlap_power use cloudfracs_type_mod, only: i_frac_liq, i_frac_ice, i_frac_bulk @@ -31,23 +33,24 @@ subroutine set_par_cloudfrac( n_points, n_points_super, & ! Number of points integer, intent(in) :: n_points -! Number of points in the cloud-fractions super-array +! Number of points in the condensate and cloud-fractions super-arrays ! (maybe larger than needed here, to save having to reallocate) -integer, intent(in) :: n_points_super +integer, intent(in) :: n_points_cond +integer, intent(in) :: n_points_cf -! Liquid and ice cloud mixing-ratios -real(kind=real_cvprec), intent(in) :: q_cl(n_points) -real(kind=real_cvprec), intent(in) :: q_cf(n_points) +! Super-array containing all condensed water mixing ratios +real(kind=real_cvprec), intent(in) :: q_cond_super & + ( n_points_cond, n_cond_species ) ! Cloud-fractions super-array real(kind=real_cvprec), intent(out) :: cloudfracs & - ( n_points_super, i_frac_liq:i_frac_bulk ) + ( n_points_cf, i_frac_liq:i_frac_bulk ) ! Sum of liquid and ice cloud mixing-ratio real(kind=real_cvprec) :: qc -! Loop counter -integer :: ic +! Loop counters +integer :: ic, i_cond do ic = 1, n_points ! Initialise cloud fractions to zero @@ -63,32 +66,50 @@ subroutine set_par_cloudfrac( n_points, n_points_super, & ! associated with each is just 0 if none present, ! 1 if any is present. - ! Set liquid cloud fraction to 1 if any liquid cloud present - do ic = 1, n_points - if ( q_cl(ic) > zero ) then - cloudfracs(ic,i_frac_liq) = one - cloudfracs(ic,i_frac_bulk) = one + do i_cond = 1, n_cond_species + ! Set liquid cloud fraction to 1 if any liquid cloud present + if ( cond_params(i_cond)%pt % i_sg == i_sg_frac_liq ) then + do ic = 1, n_points + if ( q_cond_super(ic,i_cond) > zero ) then + cloudfracs(ic,i_frac_liq) = one + cloudfracs(ic,i_frac_bulk) = one + end if + end do end if - end do - - if ( l_cv_cf ) then ! Set ice cloud fraction to 1 if any ice cloud present - do ic = 1, n_points - if ( q_cf(ic) > zero ) then - cloudfracs(ic,i_frac_ice) = one - cloudfracs(ic,i_frac_bulk) = one - end if - end do - end if + if ( cond_params(i_cond)%pt % i_sg == i_sg_frac_ice ) then + do ic = 1, n_points + if ( q_cond_super(ic,i_cond) > zero ) then + cloudfracs(ic,i_frac_ice) = one + cloudfracs(ic,i_frac_bulk) = one + end if + end do + end if + end do case (i_par_cloudfrac_mph) ! Alternative option; when liquid and ice both present, ! they may not be fully overlapped. - ! Note: not allowed to use this option if l_cv_cf is false - ! (ice-cloud mass not in use). + + ! Accumulate total liquid and ice cloud mixing-ratios + ! (temporarily storing in the cloud-fraction fields) + do i_cond = 1, n_cond_species + if ( cond_params(i_cond)%pt % i_sg == i_sg_frac_liq ) then + do ic = 1, n_points + cloudfracs(ic,i_frac_liq) = cloudfracs(ic,i_frac_liq) & + + q_cond_super(ic,i_cond) + end do + end if + if ( cond_params(i_cond)%pt % i_sg == i_sg_frac_ice ) then + do ic = 1, n_points + cloudfracs(ic,i_frac_ice) = cloudfracs(ic,i_frac_ice) & + + q_cond_super(ic,i_cond) + end do + end if + end do do ic = 1, n_points - qc = q_cl(ic) + q_cf(ic) + qc = cloudfracs(ic,i_frac_liq) + cloudfracs(ic,i_frac_ice) ! If any condensed water if ( qc > zero ) then @@ -97,8 +118,8 @@ subroutine set_par_cloudfrac( n_points, n_points_super, & ! Set liquid and ice cloud fractions to their respective ! fractions of total condensate - cloudfracs(ic,i_frac_liq) = q_cl(ic) / qc - cloudfracs(ic,i_frac_ice) = q_cf(ic) / qc + cloudfracs(ic,i_frac_liq) = cloudfracs(ic,i_frac_liq) / qc + cloudfracs(ic,i_frac_ice) = cloudfracs(ic,i_frac_ice) / qc ! This yields zero overlap between liquid and ice. ! Parameterise some overlap by raising both fractions diff --git a/science/physics_schemes/source/convection/comorph/plume_model/update_par_radius.F90 b/science/physics_schemes/source/convection/comorph/plume_model/update_par_radius.F90 index efd3b02fbe..58148c1b90 100644 --- a/science/physics_schemes/source/convection/comorph/plume_model/update_par_radius.F90 +++ b/science/physics_schemes/source/convection/comorph/plume_model/update_par_radius.F90 @@ -122,9 +122,9 @@ subroutine update_par_radius( n_points, det_mass_d, & * ( prev_rho_dry(ic) / next_rho_dry(ic) ) )**third end do - end select ! ( par_radius_vol_method ) + end select ! ( par_radius_evol_method ) -end if ! ( par_radius_evol_method==par_radius_evol_const ) +end if ! ( .NOT. par_radius_evol_method==par_radius_evol_const ) ! Limit the parcel radius to not fall below a hardwired ! safety-limit (some calculations that use it run into diff --git a/science/physics_schemes/source/convection/comorph/unit_tests/build_test_comorph.sh b/science/physics_schemes/source/convection/comorph/unit_tests/build_test_comorph.sh index ed4d872983..9282e0a060 100755 --- a/science/physics_schemes/source/convection/comorph/unit_tests/build_test_comorph.sh +++ b/science/physics_schemes/source/convection/comorph/unit_tests/build_test_comorph.sh @@ -60,6 +60,7 @@ gfortran -O0 -g -Wall -Wextra \ "$comorph/control/env_half_mod.F90" \ "$comorph/control/cloudfracs_type_mod.F90" \ "$comorph/control/turb_type_mod.F90" \ + "$comorph/control/core_combine.F90" \ "$comorph/control/parcel_type_mod.F90" \ "$comorph/control/res_source_type_mod.F90" \ "$comorph/control/subregion_mod.F90" \ @@ -101,6 +102,7 @@ gfortran -O0 -g -Wall -Wextra \ "$comorph/microphysics/collision_rate.F90" \ "$comorph/microphysics/ice_rain_to_graupel.F90" \ "$comorph/microphysics/collision_ctl.F90" \ + "$comorph/microphysics/autoconversion.F90" \ "$comorph/microphysics/microphysics_1.F90" \ "$comorph/microphysics/microphysics_2.F90" \ "$comorph/moist_proc/moist_proc.F90" \ diff --git a/science/physics_schemes/source/convection/comorph/unit_tests/build_test_moist_proc.sh b/science/physics_schemes/source/convection/comorph/unit_tests/build_test_moist_proc.sh index 8e006357fb..7c542471b8 100755 --- a/science/physics_schemes/source/convection/comorph/unit_tests/build_test_moist_proc.sh +++ b/science/physics_schemes/source/convection/comorph/unit_tests/build_test_moist_proc.sh @@ -41,6 +41,9 @@ gfortran -O0 -g -Wall -ffpe-trap=invalid,zero -fbounds-check -Warray-bounds \ "$comorph/moist_thermo/sat_adjust.F90" \ "$comorph/moist_thermo/linear_qs_mod.F90" \ \ + "$comorph/control/grid_type_mod.F90" \ + "$comorph/control/fields_type_mod.F90" \ + \ "$comorph/moist_proc/moist_proc_diags_type_mod.F90" \ "$comorph/moist_proc/phase_change_coefs_mod.F90" \ "$comorph/moist_proc/calc_phase_change_coefs.F90" \ @@ -68,6 +71,7 @@ gfortran -O0 -g -Wall -ffpe-trap=invalid,zero -fbounds-check -Warray-bounds \ "$comorph/microphysics/collision_rate.F90" \ "$comorph/microphysics/ice_rain_to_graupel.F90" \ "$comorph/microphysics/collision_ctl.F90" \ + "$comorph/microphysics/autoconversion.F90" \ "$comorph/microphysics/microphysics_1.F90" \ "$comorph/microphysics/microphysics_2.F90" \ \ diff --git a/science/physics_schemes/source/convection/comorph/unit_tests/build_test_phase_change_solve.sh b/science/physics_schemes/source/convection/comorph/unit_tests/build_test_phase_change_solve.sh new file mode 100755 index 0000000000..6d587fd050 --- /dev/null +++ b/science/physics_schemes/source/convection/comorph/unit_tests/build_test_phase_change_solve.sh @@ -0,0 +1,61 @@ +#!/bin/bash + +# Code Owner: Please refer to the UM file CodeOwners.txt +# This file belongs in section: convection_comorph + +# 1st argument $1 is path to compile directory + +# Find full path to comorph directory containing this script +comorph=$(readlink -f "$(dirname "$(readlink -f "$0")")""/..") + +echo "Compiling unit test at: $comorph/unit_tests/test_comorph.F90" + +echo "Compile directory: $1" + +cd "$1" || exit + +# ifort -O0 -g -debug all -check all -warn interface -traceback \ +gfortran -O0 -g -Wall -ffpe-trap=invalid,zero -fbounds-check -Warray-bounds \ + -fcheck-array-temporaries -finit-real=nan -fimplicit-none \ + -std=f2008ts -Wtabs -fbacktrace \ + -o test_phase_change_solve.exe \ + \ + "$comorph/control/comorph_constants_mod.F90" \ + "$comorph/interface/standalone/raise_error.F90" \ + "$comorph/interface/standalone/qsat_data.F90" \ + "$comorph/interface/standalone/set_qsat.F90" \ + "$comorph/control/set_dependent_constants.F90" \ + "$comorph/control/cmpr_type_mod.F90" \ + "$comorph/util/compress.F90" \ + "$comorph/util/check_bad_values.F90" \ + "$comorph/moist_thermo/lat_heat_mod.F90" \ + "$comorph/moist_thermo/set_dqsatdt.F90" \ + "$comorph/moist_thermo/linear_qs_mod.F90" \ + "$comorph/moist_thermo/calc_q_tot.F90" \ + "$comorph/moist_thermo/set_cp_tot.F90" \ + "$comorph/moist_thermo/calc_virt_temp_dry.F90" \ + "$comorph/moist_thermo/dry_adiabat.F90" \ + "$comorph/control/grid_type_mod.F90" \ + "$comorph/control/fields_type_mod.F90" \ + "$comorph/control/diag_type_mod.F90" \ + \ + "$comorph/moist_proc/moist_proc_diags_type_mod.F90" \ + "$comorph/moist_proc/phase_change_coefs_mod.F90" \ + "$comorph/moist_proc/calc_phase_change_coefs.F90" \ + "$comorph/moist_proc/proc_incr.F90" \ + "$comorph/moist_proc/solve_tq.F90" \ + "$comorph/moist_proc/calc_cond_temp.F90" \ + "$comorph/moist_proc/toggle_melt.F90" \ + "$comorph/moist_proc/melt_ctl.F90" \ + "$comorph/moist_proc/modify_coefs_liq.F90" \ + "$comorph/moist_proc/modify_coefs_ice.F90" \ + "$comorph/moist_proc/check_negatives.F90" \ + "$comorph/moist_proc/moist_proc_consistency.F90" \ + "$comorph/moist_proc/phase_change_solve.F90" \ + \ + "$comorph/unit_tests/test_phase_change_solve.F90" \ +1> std_out.txt 2> std_err.txt + +echo "Done. Compile output is in:" +echo "$1/std_out.txt" +echo "$1/std_err.txt" diff --git a/science/physics_schemes/source/convection/comorph/unit_tests/test_check_bad_values.F90 b/science/physics_schemes/source/convection/comorph/unit_tests/test_check_bad_values.F90 index 72b063558f..d7e8a21403 100644 --- a/science/physics_schemes/source/convection/comorph/unit_tests/test_check_bad_values.F90 +++ b/science/physics_schemes/source/convection/comorph/unit_tests/test_check_bad_values.F90 @@ -31,8 +31,6 @@ program test_check_bad_values ! (used to make error messages more informative) character(len=name_length) :: where_string character(len=name_length) :: field_name -! Flag to test for negative values. -logical :: l_positive ! Full 3-D field to test real(kind=real_hmprec) :: field(nx,ny,nz) @@ -84,7 +82,6 @@ program test_check_bad_values ! Set other inputs to check_bad_values where_string = "Unit test for check_bad_values" -l_positive = .true. ! Test whether routine can detect negative values @@ -93,14 +90,14 @@ program test_check_bad_values end do field_name = "field_with_negative" call check_bad_values_3d( lb, ub, field, where_string, & - field_name, l_positive ) + field_name, field_min=0.0_real_cvprec ) k = 1 do ic = 1, cmpr%n_points field_cmpr(ic) = real( & field( cmpr%index_i(ic), cmpr%index_j(ic), k ), real_cvprec ) end do call check_bad_values_cmpr( cmpr, k, field_cmpr, where_string, & - field_name, l_positive ) + field_name, field_min=0.0_real_cvprec ) ! Test whether routine can detect div-by-zero tmp = 0.0_real_hmprec @@ -109,14 +106,14 @@ program test_check_bad_values end do field_name = "field_with_div_by_zero" call check_bad_values_3d( lb, ub, field, where_string, & - field_name, l_positive ) + field_name ) k = 1 do ic = 1, cmpr%n_points field_cmpr(ic) = real( & field( cmpr%index_i(ic), cmpr%index_j(ic), k ), real_cvprec ) end do call check_bad_values_cmpr( cmpr, k, field_cmpr, where_string, & - field_name, l_positive ) + field_name ) ! Test whether routine can detect SQRT(negative) @@ -126,14 +123,14 @@ program test_check_bad_values end do field_name = "field_with_sqrt_negative" call check_bad_values_3d( lb, ub, field, where_string, & - field_name, l_positive ) + field_name ) k = 1 do ic = 1, cmpr%n_points field_cmpr(ic) = real( & field( cmpr%index_i(ic), cmpr%index_j(ic), k ), real_cvprec ) end do call check_bad_values_cmpr( cmpr, k, field_cmpr, where_string, & - field_name, l_positive ) + field_name ) end program test_check_bad_values diff --git a/science/physics_schemes/source/convection/comorph/unit_tests/test_phase_change_solve.F90 b/science/physics_schemes/source/convection/comorph/unit_tests/test_phase_change_solve.F90 new file mode 100644 index 0000000000..1bb67040d5 --- /dev/null +++ b/science/physics_schemes/source/convection/comorph/unit_tests/test_phase_change_solve.F90 @@ -0,0 +1,180 @@ + +#if !defined(LFRIC) +program test_phase_change_solve + +use cmpr_type_mod, only: cmpr_type +use comorph_constants_mod, only: real_cvprec, zero, name_length, & + n_cond_species, n_cond_species_liq, & + nx_full, ny_full, k_bot_conv, k_top_conv, & + k_top_init, & + melt_temp, R_dry, R_vap, & + cp_dry, cp_vap, cp_liq, cp_ice, & + L_con_ref, L_fus_ref, rho_liq, rho_ice +use set_dependent_constants_mod, only: set_dependent_constants +use linear_qs_mod, only: n_linear_qs_fields +use moist_proc_diags_type_mod, only: moist_proc_diags_type +use phase_change_solve_mod, only: phase_change_solve + +implicit none + +integer, parameter :: n_points = 1 + + +! Number of points where each condensed water species is non-zero +integer, allocatable :: nc(:) +! Indices of those points +integer, allocatable :: index_ic(:,:) + +type(cmpr_type) :: cmpr +integer :: k +character(len=name_length) :: call_string + +! Super-array containing qsat at a reference temperature, +! and dqsat/dT, for linearised qsat calculations +real(kind=real_cvprec) :: linear_qs ( n_points, n_linear_qs_fields ) + +! Height interval for this step; = the vertical distance +! between the current point and the previous point where +! prev_temp is defined. +! If integrating downwards (as is conventional for +! Eulerian calculations), it should be negative. +real(kind=real_cvprec) :: delta_z(n_points) + +! Time interval for converting process rates to increments. +! For Eulerian calculations, this is the model timestep length, +! but for Lagrangian ascents, it is the time taken for the +! parcel to rise over the height interval delta_z, +! so delta_t = delta_z/wind_w +real(kind=real_cvprec) :: delta_t(n_points) + +! Vertical wind velocity +! For Euelerian calculations, this is the vertical wind-speed. +! For Lagrangian ascents, it is the vertical velocity of the +! parcel relative to the environment, such that +! wind_w = delta_z/delta_t +real(kind=real_cvprec) :: wind_w(n_points) + +! Parcel temperature at the previous model-level +real(kind=real_cvprec) :: prev_temp(n_points) + +! Vapour exchange coefficient for each hydrometeor species +real(kind=real_cvprec), allocatable :: kq_cond(:,:) +! Heat exchange coefficient for each hydrometeor species +real(kind=real_cvprec), allocatable :: kt_cond(:,:) +! These have intent inout because they occasionally need to be +! limited for numerical safety reasons. Also we scale them +! by the timestep delta_t in this routine. + +! Fall-speed of each hydrometeor species +real(kind=real_cvprec), allocatable :: wf_cond(:,:) + +! Total freezing increment onto each ice hydrometeor species +! (includes homogeneous and heterogeneous freezing and riming) +! Needed for the hydrometeor surface heat budget, important for +! determining the melting rate +real(kind=real_cvprec), allocatable :: dq_frz_cond(:,:) + +! Local mixing ratio of each condensed water species, +! implicitly accounting for fall-out from current level / parcel +real(kind=real_cvprec), allocatable :: q_loc_cond(:,:) + +! Total available mixing ratio of each condensed water species +! (includes amount that falls through during this step, which +! maybe considerably larger than the amount actually present +! at a given instant). +! These are the values updated here; fall-out is calculated +! after this routine. +real(kind=real_cvprec), allocatable :: q_cond(:,:) + +! Total heat capacity incremented by phase-changes +real(kind=real_cvprec) :: cp_tot(n_points) + +! Parcel air temperature and water vapour mixing ratio +real(kind=real_cvprec) :: temperature(n_points) +real(kind=real_cvprec) :: q_vap(n_points) + +! Master switch for diagnostics +logical, parameter :: l_diags = .false. +! Structure storing diagnostics switches and meta-data +type(moist_proc_diags_type) :: moist_proc_diags +integer, parameter :: n_diags = 1 +real(kind=real_cvprec) :: diags_super( n_points, n_diags ) + +integer :: ic, i_cond + + + +! Not used here, but need setting to avoid error trap +nx_full = n_points +ny_full = 1 +k_bot_conv = 1 +k_top_conv = 10 +k_top_init = 9 + +! Set thermodynamics constants +melt_temp = real( 273.15, real_cvprec ) +R_dry = real( 287.05, real_cvprec ) +R_vap = real( 287.05/0.62198, real_cvprec ) +cp_dry = real( 1005.0, real_cvprec ) +cp_vap = zero +cp_liq = zero +cp_ice = zero +L_con_ref = real( 2.501e6, real_cvprec ) +L_fus_ref = real( 0.334e6, real_cvprec ) +rho_liq = real( 1000.0, real_cvprec ) +rho_ice = real( 917.0, real_cvprec ) + +! Setup constants +call set_dependent_constants() + +call_string = "test_phase_change_solve" + +allocate( nc ( n_cond_species ) ) +allocate( index_ic ( n_points, n_cond_species ) ) +allocate( cmpr%index_i(n_points) ) +allocate( cmpr%index_j(n_points) ) +cmpr%n_points = n_points +allocate( kq_cond ( n_points, n_cond_species ) ) +allocate( kt_cond ( n_points, n_cond_species ) ) +allocate( wf_cond ( n_points, n_cond_species ) ) +allocate( dq_frz_cond ( n_points, n_cond_species_liq+1 : n_cond_species ) ) +allocate( q_loc_cond ( n_points, n_cond_species ) ) +allocate( q_cond ( n_points, n_cond_species ) ) + +cmpr%index_i(1) = 1 +cmpr%index_j(1) = 1 +linear_qs(1,:) = [273.14948, 4.49374365E-3, 4.49329196E-3, & + 3.28750612E-4, 3.726164E-4] +delta_z(1) = 86.81665 +delta_t(1) = 347.2666 +wind_w(1) = 0.25 +prev_temp(1) = 273.67719 +kq_cond(1,:) = [0.17745258, 7.7014578E-5, 1.61309999E-5, 5.90306045E-5] +kt_cond(1,:) = [0.1204145, 5.54143517E-5, 1.09460952E-5, 5.06376382E-5] +wf_cond(1,:) = [2.1684817E-3, 1.8552991, 1.82522414E-3, 4.1219993] +dq_frz_cond(1,:) = [0., 1.01640635E-5] +q_loc_cond(1,:) = [4.22776575E-5, 2.29610287E-5, 3.23526983E-9, 8.08190598E-5] +q_cond(1,:) = [3.59028309E-5, 2.79465403E-5, 3.23648353E-9, 1.59455813E-4] +cp_tot(1) = 1005.0 +temperature(1) = 272.83563 +q_vap(1) = 4.62010596E-3 + +do i_cond = 1, n_cond_species + nc(i_cond) = 0.0 + do ic = 1, n_points + if ( kq_cond(ic,i_cond) > zero ) then + nc(i_cond) = nc(i_cond) + 1 + index_ic(nc(i_cond),i_cond) = ic + end if + end do +end do + +call phase_change_solve( n_points, n_points, & + nc, index_ic, cmpr, k, call_string, linear_qs, & + delta_z, delta_t, wind_w, prev_temp, & + kq_cond, kt_cond, wf_cond, dq_frz_cond, & + q_loc_cond, q_cond, cp_tot, temperature, q_vap, & + l_diags, moist_proc_diags, n_points, n_diags, diags_super ) + +end program test_phase_change_solve +#endif diff --git a/science/physics_schemes/source/convection/comorph/util/check_bad_values.F90 b/science/physics_schemes/source/convection/comorph/util/check_bad_values.F90 index 841c82434c..69aa19ec54 100644 --- a/science/physics_schemes/source/convection/comorph/util/check_bad_values.F90 +++ b/science/physics_schemes/source/convection/comorph/util/check_bad_values.F90 @@ -20,8 +20,8 @@ module check_bad_values_mod !---------------------------------------------------------------- ! Bad value check for full 3-D fields !---------------------------------------------------------------- -subroutine check_bad_values_3d( lb, ub, field, where_string, & - field_name, l_positive, l_half, l_init ) +subroutine check_bad_values_3d( lb, ub, field, where_string, field_name, & + field_min, field_max, l_half, l_init ) use comorph_constants_mod, only: real_hmprec, real_cvprec, name_length, & nx_full, ny_full, k_bot_conv, k_top_conv, & @@ -44,9 +44,9 @@ subroutine check_bad_values_3d( lb, ub, field, where_string, & ! Name of the field (for error message) character(len=name_length), intent(in) :: field_name -! Flag for whether the field is positive-only -! (in which case we check for negative values) -logical, intent(in) :: l_positive +! Optionally specify the min and max allowed values for the current field: +real(kind=real_cvprec), optional, intent(in) :: field_min +real(kind=real_cvprec), optional, intent(in) :: field_max ! Flag for the field being checked is on half-levels; in this case, ! the uppermost model-level used is k_top_conv+1, since this is the @@ -103,8 +103,8 @@ subroutine check_bad_values_3d( lb, ub, field, where_string, & ! Note: crucially this will catch instances where the value in field ! was representable in the host-model at 64-bit, but becomes garbage ! when converted to comorph native precision at 32-bit. - call check_bad_values_cmpr( cmpr, k, field_cmpr, where_string, & - field_name, l_positive, & + call check_bad_values_cmpr( cmpr, k, field_cmpr, where_string, field_name, & + field_min, field_max, & i_check_bad=i_check_bad_values_3d ) end do @@ -120,11 +120,11 @@ end subroutine check_bad_values_3d !---------------------------------------------------------------- ! Bad value check for compressed 1-D fields !---------------------------------------------------------------- -subroutine check_bad_values_cmpr( cmpr, k, field, where_string, & - field_name, l_positive, i_check_bad ) +subroutine check_bad_values_cmpr( cmpr, k, field, where_string, field_name, & + field_min, field_max, i_check_bad ) -use comorph_constants_mod, only: real_cvprec, zero, name_length, newline, & - i_check_bad_values_cmpr +use comorph_constants_mod, only: real_cvprec, max_float, & + name_length, newline, i_check_bad_values_cmpr use cmpr_type_mod, only: cmpr_type implicit none @@ -144,9 +144,9 @@ subroutine check_bad_values_cmpr( cmpr, k, field, where_string, & ! Name of the field (for error message) character(len=name_length), intent(in) :: field_name -! Flag for whether the field is positive-only -! (in which case we check for negative values) -logical, intent(in) :: l_positive +! Optionally specify the min and max allowed values for the current field: +real(kind=real_cvprec), optional, intent(in) :: field_min +real(kind=real_cvprec), optional, intent(in) :: field_max ! Optionally override default setting for whether to do warning or fatal error integer, optional, intent(in) :: i_check_bad @@ -158,8 +158,6 @@ subroutine check_bad_values_cmpr( cmpr, k, field, where_string, & real(kind=real_cvprec) :: min_val real(kind=real_cvprec) :: max_val -real(kind=real_cvprec), parameter :: large_number = huge(field) - ! Warning vs fatal error switch to use integer :: i_check_bad_use @@ -167,17 +165,14 @@ subroutine check_bad_values_cmpr( cmpr, k, field, where_string, & integer :: ic -! Set max allowed value to the max possible floating point real -! (so really just checking for Infinities) -max_val = large_number -if ( l_positive ) then - ! If the input is supposed to be positive-only, set minimum - ! allowed value to zero - min_val = zero -else - ! Otherwise set to largest possible negative to just check for -Inf - min_val = -large_number -end if +! Set default allowed bounds to the min and max possible floating-point +! values (so just checking for NaN or Inf): +min_val = -max_float +max_val = max_float + +! Override with more stringent limits if these are input: +if ( present(field_min) ) min_val = field_min +if ( present(field_max) ) max_val = field_max ! Loop over points do ic = 1, cmpr%n_points