diff --git a/exp/test_cases/socrates_mars_dust/input/README.md b/exp/test_cases/socrates_mars_dust/input/README.md new file mode 100644 index 000000000..d8d9ee41c --- /dev/null +++ b/exp/test_cases/socrates_mars_dust/input/README.md @@ -0,0 +1,45 @@ +# Mars dust + latent-heating input files + +`socrates_mars_dust_test_case.py` requires the following files to be placed in this +folder before it can be run: + +* `sp_lw_17_dsa_mars_dust` +* `sp_lw_17_dsa_mars_dust_k` +* `sp_sw_42_dsa_mars_sun_dust` +* `sp_sw_42_dsa_mars_sun_dust_k` +* `t42_mola_mars.nc` +* `cdod_clim_MY24.nc` .. `cdod_clim_MY34.nc` (one file per Mars Year, 24-34) +* `cdod_clim_scenario.nc` +* `cdod_cold.nc` +* `cdod_warm_25.nc` +* `cdod_all_years.nc` +* `cdod_all_years_long.nc` + +These are not included in this repository (the dust-climatology files alone are around +100MB). They can be found in Emily Ball's own fork, which this test case is ported +from: . + +## What these files are + +* The `sp_*_dsa_mars_dust*` files are Mars-specific Socrates spectral files with dust + optical properties included, distinct from the plain (non-dust) Mars spectral files + used by the `socrates_mars` test case. +* `t42_mola_mars.nc` is MOLA-derived Mars topography, also used by `socrates_mars`. +* The `cdod_*.nc` files are column dust optical depth climatologies, used to set the + reference dust mass mixing ratio via `do_read_cdod`/`cdod_file_name` in + `socrates_rad_nml`. `cdod_field_name` selects which variable within the chosen file + to read. + +## Physics this exercises + +This test case switches on two capabilities from Emily Ball's PhD work (Ball et al. +2021, "The roles of latent heating and dust in the structure and variability of the +northern Martian polar vortex"), both gated behind their own namelist flags and off by +default everywhere else in Isca: + +* `do_lscale_cond_lh` (`idealized_moist_phys_nml`) - a CO2-condensation latent-heating + condensation scheme, an alternative to the standard `do_lscale_cond`. +* `do_dust_forcing` (`socrates_rad_nml`) - builds a Conrath-type Mars dust vertical + profile from the dust optical depth climatology and feeds it to Socrates as a + radiatively-active aerosol. `dust_scale` controls the optical-depth-to-mass-mixing- + ratio conversion. diff --git a/exp/test_cases/socrates_mars_dust/socrates_mars_dust_test_case.py b/exp/test_cases/socrates_mars_dust/socrates_mars_dust_test_case.py new file mode 100644 index 000000000..20d725734 --- /dev/null +++ b/exp/test_cases/socrates_mars_dust/socrates_mars_dust_test_case.py @@ -0,0 +1,281 @@ +import os + +import numpy as np + +from isca import SocratesCodeBase, DiagTable, Experiment, Namelist, GFDL_BASE + +NCORES = 8 +base_dir = os.path.dirname(os.path.realpath(__file__)) + +# a CodeBase can be a directory on the computer, +# useful for iterative development +cb = SocratesCodeBase.from_directory(GFDL_BASE) + +# or it can point to a specific git repo and commit id. +# This method should ensure future, independent, reproducibility of results. +# cb = SocratesCodeBase.from_repo(repo='https://github.com/isca/isca', commit='isca1.1') + +# compilation depends on computer specific settings. The $GFDL_ENV +# environment variable is used to determine which `$GFDL_BASE/src/extra/env` file +# is used to load the correct compilers. The env file is always loaded from +# $GFDL_BASE and not the checked out git repo. + +# NOTE: this test case requires the Mars-dust-specific Socrates spectral files and dust +# climatology files below to be present in the `input` folder next to this script: +# sp_lw_17_dsa_mars_dust, sp_lw_17_dsa_mars_dust_k, +# sp_sw_42_dsa_mars_sun_dust, sp_sw_42_dsa_mars_sun_dust_k, +# t42_mola_mars.nc, cdod_clim_MY24.nc .. cdod_clim_MY34.nc, cdod_clim_scenario.nc, +# cdod_cold.nc, cdod_warm_25.nc, cdod_all_years.nc, cdod_all_years_long.nc +# These are not currently included in the repository - see input/README.md. +dust_clim = 'cdod_clim_MY24' + +inputfiles = [os.path.join(base_dir,'input/sp_lw_17_dsa_mars_dust'), + os.path.join(base_dir,'input/sp_sw_42_dsa_mars_sun_dust'), + os.path.join(base_dir,'input/sp_lw_17_dsa_mars_dust_k'), + os.path.join(base_dir,'input/sp_sw_42_dsa_mars_sun_dust_k'), + os.path.join(base_dir,'input/t42_mola_mars.nc'), + os.path.join(base_dir,'input/cdod_clim_MY24.nc'), + os.path.join(base_dir,'input/cdod_clim_MY25.nc'), + os.path.join(base_dir,'input/cdod_clim_MY26.nc'), + os.path.join(base_dir,'input/cdod_clim_MY27.nc'), + os.path.join(base_dir,'input/cdod_clim_MY28.nc'), + os.path.join(base_dir,'input/cdod_clim_MY29.nc'), + os.path.join(base_dir,'input/cdod_clim_MY30.nc'), + os.path.join(base_dir,'input/cdod_clim_MY31.nc'), + os.path.join(base_dir,'input/cdod_clim_MY32.nc'), + os.path.join(base_dir,'input/cdod_clim_MY33.nc'), + os.path.join(base_dir,'input/cdod_clim_MY34.nc'), + os.path.join(base_dir,'input/cdod_clim_scenario.nc'), + os.path.join(base_dir,'input/cdod_cold.nc'), + os.path.join(base_dir,'input/cdod_warm_25.nc'), + os.path.join(base_dir,'input/cdod_all_years_long.nc'), + os.path.join(base_dir,'input/cdod_all_years.nc')] + +# create an Experiment object to handle the configuration of model parameters +# and output diagnostics +exp = Experiment('socrates_mars_dust_test_experiment', codebase=cb) + +#Tell model how to write diagnostics +diag = DiagTable() +diag.add_file('atmos_daily', 88440, 'seconds', time_units='days') + +#Tell model which diagnostics to write +diag.add_field('dynamics', 'ps', time_avg=True) +diag.add_field('dynamics', 'bk', time_avg=True) +diag.add_field('dynamics', 'pk', time_avg=True) +diag.add_field('dynamics', 'zsurf', time_avg=True) +diag.add_field('dynamics', 'ucomp', time_avg=True) +diag.add_field('dynamics', 'vcomp', time_avg=True) +diag.add_field('dynamics', 'temp', time_avg=True) +diag.add_field('dynamics', 'omega', time_avg=True) +diag.add_field('dynamics', 'height', time_avg=True) +diag.add_field('dynamics', 'sphum', time_avg=True) + +diag.add_field('atmosphere', 'precipitation', time_avg=True) +diag.add_field('atmosphere', 'dt_tg_convection', time_avg=True) +diag.add_field('atmosphere', 'lh_rel', time_avg=True) +diag.add_field('atmosphere', 'dt_tg_lh_condensation', time_avg=True) + +diag.add_field('mixed_layer', 't_surf', time_avg=True) +diag.add_field('mixed_layer', 'flux_lhe', time_avg=True) +diag.add_field('mixed_layer', 'flux_t', time_avg=True) + +diag.add_field('socrates', 'mars_solar_long', time_avg=True) +diag.add_field('socrates', 'soc_coszen', time_avg=True) +diag.add_field('socrates', 'rrsun', time_avg=True) +diag.add_field('socrates', 'soc_toa_sw_down', time_avg=True) +diag.add_field('socrates', 'time_since_ae', time_avg=True) +diag.add_field('socrates', 'true_anomaly', time_avg=True) +diag.add_field('socrates', 'dec', time_avg=True) +diag.add_field('socrates', 'ang', time_avg=True) +diag.add_field('socrates', 'soc_dust', time_avg=True) +diag.add_field('socrates', 'soc_cdod', time_avg=True) + +exp.diag_table = diag +exp.inputfiles = inputfiles + +#Empty the run directory ready to run +exp.clear_rundir() + +#Define values for the 'core' namelist. Physical constants (grav, pstd, radius, etc.) are Mars' +#values, taken from Emily Ball's Mars dust/latent-heating setup (Ball et al. 2021, +#github.com/emilyrball/Isca-Mars), reused unchanged here. +exp.namelist = namelist = Namelist({ + 'main_nml': { + 'dt_atmos': 110, + 'days': 0, + 'seconds': 30*88440, + 'calendar': 'no_calendar' + }, + + 'idealized_moist_phys_nml': { + 'do_damping': True, + 'turb':True, + 'mixed_layer_bc':True, + 'do_virtual' :False, + 'do_simple': True, + 'roughness_mom':3.21e-05, + 'roughness_heat':3.21e-05, + 'roughness_moist':0., + 'two_stream_gray': False, + 'do_lscale_cond': False, + 'do_lscale_cond_lh': True, + 'do_socrates_radiation': True, + 'convection_scheme': 'none', + }, + + 'vert_turb_driver_nml': { + 'do_mellor_yamada': False, # default: True + 'do_diffusivity': True, # default: False + 'do_simple': True, # default: False + 'constant_gust': 0.0, # default: 1.0 + 'use_tau': False + }, + + 'diffusivity_nml': { + 'do_entrain':False, + 'do_simple': True, + }, + + 'surface_flux_nml': { + 'use_virtual_temp': False, + 'do_simple': True, + 'old_dtaudv': True, + 'use_actual_surface_temperatures':False, + }, + + 'atmosphere_nml': { + 'idealized_moist_model': True + }, + + 'mixed_layer_nml': { + 'tconst' : 285., + 'prescribe_initial_dist':True, + 'evaporation':False, + 'albedo_value': 0.3, + 'depth': 2.0, + }, + + 'qe_moist_convection_nml': { + 'rhbm':0.0, + 'tau_bm':3600., + }, + + 'dry_convection_nml': { + 'tau':7200, + }, + + 'betts_miller_nml': { + 'rhbm': .7, + 'do_simp': False, + 'do_shallower': True, + }, + + 'lscale_cond_nml': { + 'do_simple':True, + 'do_evap':True + }, + + 'sat_vapor_pres_nml': { + 'do_simple':True, + 'tcmin': -223, #Make sure low temperature limit of saturation vapour pressure is low enough that it doesn't cause an error (this planet has no moisture anyway, so doesn't directly affect the calculation). + 'tcmax': 350, # NOTE: tcmin/tcmax are Fortran INTEGERs (sat_vapor_pres.F90) - a trailing '.' here makes this a Python float, which gfortran correctly rejects as a namelist type mismatch (FATAL "Unknown namelist, or mistyped namelist variable in namelist sat_vapor_pres_nml"). Keep this as an int. + }, + + 'damping_driver_nml': { + 'do_rayleigh': True, + 'trayfric': -0.125, # neg. value: time in *days* + 'sponge_pbottom': 0.5, #Bottom of the model's sponge down to 50hPa (units are Pa) + 'do_conserve_energy': True, + }, + + 'spectral_dynamics_nml': { + 'num_levels': 25, + 'exponent': 2.5, + 'scale_heights': 4, + 'surf_res': 0.1, + 'robert_coeff': 4e-2, + 'do_water_correction': False, + 'vert_coord_option': 'input', + 'initial_sphum': 0., + 'valid_range_T': [0, 700], + 'ocean_topog_smoothing': 0.8, + 'reference_sea_level_press': 610.0, + }, + + 'vert_coordinate_nml': { + 'bk': [ 0.00000000e+00, 1.53008955e-04, 4.63790800e-04, 1.10977640e-03, 2.37044150e-03, 4.74479200e-03, 9.12245300e-03, 1.70677050e-02, 3.12516100e-02, 5.59939500e-02, 9.76165000e-02, 1.63754900e-01, 2.60315150e-01, 3.85974250e-01, 5.28084800e-01, 6.65956600e-01, 7.81088000e-01, 8.65400050e-01, 9.21109250e-01, 9.55343500e-01, 9.75416950e-01, 9.86856800e-01, 9.93269300e-01, 9.96830200e-01, 9.98799150e-01, 1.00000000e+00], + 'pk': [0.]*26, + }, + + 'socrates_rad_nml': { + 'stellar_constant':589., + 'lw_spectral_filename':'INPUT/sp_lw_17_dsa_mars_dust', + 'sw_spectral_filename':'INPUT/sp_sw_42_dsa_mars_sun_dust', + 'do_read_ozone': False, + 'dt_rad':1320, + 'store_intermediate_rad':True, + 'chunk_size': 16, + 'use_pressure_interp_for_half_levels':False, + 'tidally_locked':False, + 'inc_h2o':False, + 'inc_o3':False, + 'inc_co2':True, + 'inc_n2':True, + 'equinox_day': 0.0, + 'co2_ppmv': 0.949*1e6, + 'n2_mix_ratio': 0.026*(28.0134)/(1000.*8.314/192.0), + 'do_read_cdod': True, + 'cdod_file_name': dust_clim, + 'cdod_field_name': 'cdod', + # do_dust_forcing and account_for_effect_of_dust both need to be True to get + # the full Mars dust scheme (Ball et al. 2021) as intended here: + # - do_dust_forcing is the master switch for whether the Conrath-type dust + # vertical profile is computed at all. If False, there's no dust profile + # regardless of account_for_effect_of_dust below. + # - account_for_effect_of_dust then decides whether that computed profile is + # actually passed into the radiative transfer calculation, or discarded at + # that last step. Set it False (with do_dust_forcing still True) to run with + # the dust profile computed and diagnosed (soc_dust, dust_mmr_ref) but not + # radiatively active - e.g. to isolate dust's radiative effect from runs + # that only differ by this flag. + # See the matching comments in socrates_config_mod.f90/socrates_interface.F90 + # for the Fortran-side detail. + 'account_for_effect_of_dust': True, + 'do_dust_forcing': True, + 'dust_scale': 7.4e-5, + }, + + 'astronomy_nml': { + 'ecc':0.0935, + 'obliq':25.19, + 'per':70.85, + 'use_mean_anom_in_rrsun_calc':True, + 'use_old_r_inv_squared':False + }, + + 'constants_nml': { + 'orbital_period': 59166360, + 'solar_const':589.0, + 'radius':3396.0e3, + 'rdgas':192.0, + 'kappa':0.22727, + 'rotation_period':88308, + 'grav': 3.71, + 'pstd': 6100.0, + 'pstd_mks': 610.0, + }, + + 'spectral_init_cond_nml': { + 'topog_file_name': 't42_mola_mars.nc', + 'topography_option': 'input', + }, +}) + +#Lets do a run! +if __name__=="__main__": + cb.compile() # compile the source code to working directory $GFDL_WORK/codebase + + exp.run(1, use_restart=False, num_cores=NCORES) + for i in range(2, 65): + exp.run(i, num_cores=NCORES) diff --git a/exp/test_cases/trip_test/trip_test_functions.py b/exp/test_cases/trip_test/trip_test_functions.py index 2afbee836..37df1c1d8 100644 --- a/exp/test_cases/trip_test/trip_test_functions.py +++ b/exp/test_cases/trip_test/trip_test_functions.py @@ -15,6 +15,18 @@ import sys import f90nml +# Mars test cases (grey_mars, radiative_eq_mars, socrates_mars, socrates_mars_dust) use dt_atmos=110s, +# chosen to divide evenly into a Mars "day" of MARS_DAY_LENGTH_SECONDS=88440s (804 steps) - a calendar +# convention (see e.g. grey_mars_test_case.py) rather than the true Martian sol (rotation_period=88308s +# in those same namelists), chosen so an integer number of days fits into a Mars year. main_nml's +# days/seconds units and DiagTable's 'days' time_units always mean Earth's fixed 86400s/day regardless +# of planet, and 86400 is not a multiple of 110 - so the generic short-test-run overrides below (which +# assume a run length in whole Earth days divides evenly by dt_atmos) need a Mars-specific equivalent. +MARS_DAY_LENGTH_SECONDS = 88440 + +def is_mars_test_case(test_case_name): + return any(name in test_case_name for name in ('grey_mars', 'radiative_eq_mars', 'socrates_mars')) + def get_nml_diag(test_case_name): """Gets the appropriate namelist and input files from each of the test case scripts in the test_cases folder """ @@ -72,6 +84,15 @@ def get_nml_diag(test_case_name): nml_out = exp_temp.namelist codebase_to_use = SocratesCodeBase + if 'socrates_mars_dust' in test_case_name: + # Note: 'socrates_mars_dust' also matches the 'socrates_mars' check above - this + # block runs afterwards and overwrites nml_out/input_files/codebase_to_use, so it wins. + sys.path.insert(0, os.path.join(GFDL_BASE, 'exp/test_cases/socrates_mars_dust/')) + from socrates_mars_dust_test_case import exp as exp_temp + input_files = exp_temp.inputfiles + nml_out = exp_temp.namelist + codebase_to_use = SocratesCodeBase + if 'column_test' in test_case_name: sys.path.insert(0, os.path.join(GFDL_BASE, 'exp/test_cases/column_test_case/')) from column_test_case import exp as exp_temp @@ -212,6 +233,7 @@ def list_all_test_cases_implemented_in_trip_test(): 'grey_mars', 'radiative_eq_mars', #'socrates_mars', # requires Mars-specific Socrates spectral files not yet included in the repo - see exp/test_cases/socrates_mars/input/README.md + #'socrates_mars_dust', # requires Mars-dust-specific Socrates spectral files and dust climatology data not yet included in the repo - see exp/test_cases/socrates_mars_dust/input/README.md #'frierson_dry_heating', # exercises the new local_heating feature, which the current ExeClim master doesn't have - included for opt-in testing, not run by default ] @@ -235,6 +257,26 @@ def define_simple_diag_table(): return diag +def define_simple_diag_table_mars(): + """Defines a simple diag table for the Mars test cases, with the output cadence set to + MARS_DAY_LENGTH_SECONDS (one Mars day) rather than the generic 1 Earth day used elsewhere - + see the MARS_DAY_LENGTH_SECONDS comment above for why.""" + + diag = DiagTable() + diag.add_file('atmos_daily', MARS_DAY_LENGTH_SECONDS, 'seconds', time_units='days') + + #Tell model which diagnostics to write + diag.add_field('dynamics', 'ps', time_avg=True) + diag.add_field('dynamics', 'bk') + diag.add_field('dynamics', 'pk') + diag.add_field('dynamics', 'ucomp', time_avg=True) + diag.add_field('dynamics', 'vcomp', time_avg=True) + diag.add_field('dynamics', 'temp', time_avg=True) + diag.add_field('dynamics', 'vor', time_avg=True) + diag.add_field('dynamics', 'div', time_avg=True) + + return diag + def define_simple_diag_table_2d(shallow_or_baro): """Defines a simple diag table for the shallow water and barotropic vorticity test cases.""" @@ -304,6 +346,8 @@ def conduct_comparison_on_test_case(base_commit, later_commit, test_case_name, r diag_use = define_simple_diag_table_2d('barotropic') elif 'column_test' in test_case_name: diag_use = define_simple_diag_table_column() + elif is_mars_test_case(test_case_name): + diag_use = define_simple_diag_table_mars() else: diag_use = define_simple_diag_table() @@ -322,11 +366,21 @@ def conduct_comparison_on_test_case(base_commit, later_commit, test_case_name, r exp.diag_table = diag_use exp.inputfiles = input_files_use - #Only run for 3 days to keep things short. - exp.update_namelist({ - 'main_nml': { - 'days': 3, - }}) + #Only run for a short time to keep things short. + if is_mars_test_case(test_case_name): + #Override in whole Mars days (see MARS_DAY_LENGTH_SECONDS above) rather than + #'days', since main_nml's 'days' always means Earth's fixed 86400s regardless + #of planet, and 86400 is not a multiple of these test cases' dt_atmos=110s. + exp.update_namelist({ + 'main_nml': { + 'days': 0, + 'seconds': 3*MARS_DAY_LENGTH_SECONDS, + }}) + else: + exp.update_namelist({ + 'main_nml': { + 'days': 3, + }}) except: run_complete = False test_pass = False diff --git a/src/atmos_param/lscale_cond/lscale_cond.F90 b/src/atmos_param/lscale_cond/lscale_cond.F90 index 231c03803..1341f615d 100644 --- a/src/atmos_param/lscale_cond/lscale_cond.F90 +++ b/src/atmos_param/lscale_cond/lscale_cond.F90 @@ -34,7 +34,7 @@ module lscale_cond_mod !----------------------------------------------------------------------- ! ---- public interfaces ---- - public lscale_cond, lscale_cond_init, lscale_cond_end + public lscale_cond, lscale_cond_co2_mars, lscale_cond_init, lscale_cond_end !----------------------------------------------------------------------- ! ---- version number ---- @@ -254,6 +254,139 @@ subroutine precip_evap (pmass, tin, qin, qsat, dqsat, hlcp, & end subroutine precip_evap +!####################################################################### + + subroutine lscale_cond_co2_mars (tin, qin, pfull, phalf, lh_rel, coldT, & + rain, snow, tdel, qdel, mask, conv) + +!----------------------------------------------------------------------- +! +! Mars-specific CO2 condensation (Way et al. 2017), ported from +! Emily Ball's lscale_cond_lh (github.com/emilyrball/Isca-Mars). +! +! Unlike lscale_cond above, which condenses water vapor wherever +! specific humidity exceeds a saturation curve, this adjusts +! temperature toward a pressure-dependent CO2 frost point and +! diagnoses the associated latent heat release. Not applicable to +! Earth or other non-CO2-dominated atmospheres - only wire this up +! behind its own namelist option, never as a default code path. +! +!----------------------------------------------------------------------- +! +! input: tin temperature at full model levels +! qin specific humidity of water vapor at full +! model levels (unused - retained for interface +! compatibility with lscale_cond) +! pfull pressure at full model levels +! phalf pressure at half (interface) model levels +! coldT should precipitation be snow at this point? +! optional: +! mask optional mask (0 or 1.) (unused) +! conv logical flag (unused) +! +! output: rain liquid precipitation (kg/m2) - always zero; see +! qdel note below +! snow frozen precipitation (kg/m2) - always zero; see +! qdel note below +! tdel temperature tendency at full model levels +! qdel specific humidity tendency - always zero; see note +! below +! lh_rel latent heat released by CO2 condensation (J/kg) +! +!----------------------------------------------------------------------- +!--------------------- interface arguments ----------------------------- + + real , intent(in) , dimension(:,:,:) :: tin, qin, pfull, phalf + logical , intent(in) , dimension(:,:):: coldT + real , intent(out), dimension(:,:) :: rain,snow + real , intent(out), dimension(:,:,:) :: tdel, qdel, lh_rel + real , intent(in) , dimension(:,:,:), optional :: mask + logical, intent(in) , dimension(:,:,:), optional :: conv +!----------------------------------------------------------------------- +!---------------------- local data ------------------------------------- + + real,dimension(size(tin,1),size(tin,2),size(tin,3)) :: pmass, tcond, mdel + integer k, kx, j, jx, i, ix + + ! Reference latent heat of CO2 condensation at T0=150K (J/kg), the + ! anchor point for the Kirchhoff-type expansion of lh_rel below. + real, parameter :: L_c0 = 5.902e05 + ! Specific heat capacity of gaseous CO2 (J/kg/K), evaluated near the + ! ~150K polar-cap condensation temperature - NOT the same value as + ! CP_AIR (constants_mod), which is derived from the model's configured + ! rdgas/kappa and represents Mars' bulk/mean atmospheric temperature + ! (~845 J/kg/K for the existing Mars test cases, vs 770.2 here). + ! two_stream_gray_rad.F90, socrates_interface.F90, spectral_dynamics.F90, + ! mixed_layer.F90 and hs_forcing.F90 all use CP_AIR for the bulk + ! atmosphere's energy budget - this hasn't been reconciled with c_pg/c_p + ! below, so a Mars run currently has two different atmospheric specific + ! heats in play simultaneously. Check before trusting lh_rel/mdel + ! quantitatively. + real, parameter :: c_pg = 770.2 + ! Specific heat capacity of solid CO2 ice (J/kg/K). + real, parameter :: c_pc = 1070.7 + ! Specific heat capacity used to convert the diagnosed temperature + ! adjustment into an implied condensed CO2 mass (see mdel below). + ! Same CP_AIR consistency caveat as c_pg above. + real, parameter :: c_p = 735.9 + +!----------------------------------------------------------------------- +! computation of latent heating from CO2 condensation +!----------------------------------------------------------------------- + + if (.not. module_is_initialized) call error_mesg ('lscale_cond_co2_mars', & + 'lscale_cond_init has not been called.', FATAL) + + kx=size(tin,3) + jx=size(tin,2) + ix=size(tin,1) + + tcond(:,:,:) = 149.2+6.48*LOG(0.00135*pfull(:,:,:)) ! CO2 condensation temperature, Way 2017 + + do i=1,ix + do j=1,jx + do k=1,kx + if (tin(i,j,k) < tcond(i,j,k)) then + tdel(i,j,k) = tcond(i,j,k)-tin(i,j,k) + lh_rel(i,j,k) = L_c0 + c_pg*(tin(i,j,k)-150.) & + - c_pc*(tcond(i,j,k)-150.) + else + tdel(i,j,k)=0.0 + lh_rel(i,j,k)=0.0 + endif + enddo + enddo + enddo + +!------------ pressure mass of each layer ------------------------------ + + do k=1,kx + pmass(:,:,k)=(phalf(:,:,k+1)-phalf(:,:,k))/Grav + do j=1,jx + do i=1,ix + mdel(i,j,k)=c_p*pmass(i,j,k)/L_c0 * tdel(i,j,k) + lh_rel(i,j,k) = lh_rel(i,j,k)*mdel(i,j,k) + enddo + enddo + enddo + +!----------------------------------------------------------------------- +! NOTE: this scheme does not yet track the CO2 mass condensing as a +! tracer tendency - it only adjusts temperature (tdel) and diagnoses the +! associated latent heat (lh_rel). qdel/rain/snow are set to zero (they +! were previously left uninitialized - do_evap/precip_evap depended on +! qsat/dqsat/hlcp, which this formulation never computes, so that path +! has been dropped rather than run against undefined data). +!----------------------------------------------------------------------- + + qdel(:,:,:) = 0.0 + rain(:,:) = 0.0 + snow(:,:) = 0.0 + +!----------------------------------------------------------------------- + + end subroutine lscale_cond_co2_mars + !####################################################################### subroutine lscale_cond_init () diff --git a/src/atmos_param/socrates/interface/read_control.F90 b/src/atmos_param/socrates/interface/read_control.F90 index c8d3785e5..fb8ac422b 100644 --- a/src/atmos_param/socrates/interface/read_control.F90 +++ b/src/atmos_param/socrates/interface/read_control.F90 @@ -15,9 +15,10 @@ SUBROUTINE read_control(control, spectrum, do_clouds) USE def_control, ONLY: StrCtrl, allocate_control USE def_spectrum, ONLY: StrSpecData -USE socrates_config_mod, ONLY: l_planet_grey_surface, inc_h2o, inc_co2, inc_co, & +USE socrates_config_mod, ONLY: l_planet_grey_surface, inc_h2o, inc_co2, inc_co, & inc_o3, inc_n2o, inc_ch4, inc_o2, inc_so2, inc_cfc11, & - inc_cfc12, inc_cfc113, inc_hcfc22, inc_hfc134a, inc_n2 + inc_cfc12, inc_cfc113, inc_hcfc22, inc_hfc134a, inc_n2, & + do_dust_forcing IMPLICIT NONE @@ -178,7 +179,19 @@ SUBROUTINE read_control(control, spectrum, do_clouds) END IF ! Aerosols -control%l_aerosol = .FALSE. +! Only switch on aerosol radiative effects when the Mars dust forcing scheme +! (Ball et al. 2021) is actually requested via do_dust_forcing. Note this +! can NOT be gated on spectrum%aerosol%n_aerosol > 0: the standard ga7 +! spectral files used by every non-Mars Socrates run also carry aerosol +! species blocks (inherited from the Met Office UM), but with no +! i_aerosol_parametrization wired up for Isca's use of them - switching +! l_aerosol on for those crashes Socrates' opt_prop_aerosol with "illegal +! value" for those species. do_dust_forcing is the actual user-facing +! namelist flag for whether this capability is wanted - deliberately NOT +! account_for_effect_of_dust (socrates_config_mod.f90), which only decides +! whether an already-computed dust profile is fed to the radiation call and +! is meaningless when do_dust_forcing is off in the first place. +control%l_aerosol = do_dust_forcing control%l_aerosol_mode = .FALSE. control%l_aerosol_ccn = .FALSE. diff --git a/src/atmos_param/socrates/interface/socrates_calc.F90 b/src/atmos_param/socrates/interface/socrates_calc.F90 index bc40be70f..47e109623 100644 --- a/src/atmos_param/socrates/interface/socrates_calc.F90 +++ b/src/atmos_param/socrates/interface/socrates_calc.F90 @@ -26,7 +26,7 @@ subroutine socrates_calc(Time_diag,control, spectrum, & n_profile, n_layer, n_cloud_layer, n_aer_mode, & cld_subcol_gen, cld_subcol_req, & p_layer, t_layer, t_layer_boundaries, d_mass, density, & - h2o, o3, co2, & + h2o, o3, co2, dust, & t_rad_surf, cos_zenith_angle, solar_irrad, orog_corr, & l_planet_grey_surface, planet_albedo, planet_emissivity, & layer_heat_capacity, & @@ -50,7 +50,7 @@ subroutine socrates_calc(Time_diag,control, spectrum, & use set_atm_mod, only: set_atm use set_bound_mod, only: set_bound use socrates_set_cld, only: set_cld -use set_aer_mod, only: set_aer +use socrates_set_aer, only: set_aer use soc_constants_mod, only: i_def, r_def @@ -95,6 +95,9 @@ subroutine socrates_calc(Time_diag,control, spectrum, & real(r_def), intent(in) :: co2(n_profile, n_layer) ! Mass mixing ratio of carbon dioxide +real(r_def), intent(in) :: dust(n_profile, n_layer) +! Mass mixing ratio of dust + real(r_def), intent(in) :: t_rad_surf(n_profile) ! Effective radiative temperature over whole grid-box real(r_def), intent(in) :: cos_zenith_angle(n_profile) @@ -195,7 +198,7 @@ subroutine socrates_calc(Time_diag,control, spectrum, & liq_dim = reff_rad, & ice_dim = zeros_cld ) -call set_aer(control, dimen, spectrum, aer, n_profile) +call set_aer(control, dimen, spectrum, aer, n_profile, n_layer, dust) ! DEPENDS ON: radiance_calc call radiance_calc(control, dimen, spectrum, atm, cld, aer, bound, radout) diff --git a/src/atmos_param/socrates/interface/socrates_config_mod.f90 b/src/atmos_param/socrates/interface/socrates_config_mod.f90 index f7588cf0b..e7f582427 100644 --- a/src/atmos_param/socrates/interface/socrates_config_mod.f90 +++ b/src/atmos_param/socrates/interface/socrates_config_mod.f90 @@ -29,13 +29,37 @@ module socrates_config_mod character(len=256) :: sw_hires_spectral_filename='unset' logical :: account_for_effect_of_water=.TRUE. !if False then radiation is fed water mixing ratios = 0. If true it's fed mixing ratios based on model specific humidity. logical :: account_for_effect_of_ozone=.TRUE. !if False then radiation is fed ozone mixing ratios = 0. If true it's fed mixing ratios based on model ozone field. + ! account_for_effect_of_dust vs do_dust_forcing (below): two different gates on the + ! same field, applied at two different stages. do_dust_forcing controls whether a + ! Mars dust profile is COMPUTED at all (run_socrates); account_for_effect_of_dust + ! controls whether that computed profile is then actually PASSED to the radiative + ! transfer call (socrates_interface), vs. zeroed at that last step. So: + ! do_dust_forcing=F -> no dust, regardless of this flag (both default here) + ! do_dust_forcing=T, account_for_...=T -> dust computed AND radiatively active (the normal Mars dust case - + ! both must be set explicitly, e.g. socrates_mars_dust_test_case.py) + ! do_dust_forcing=T, account_for_...=F -> dust computed and diagnosed (soc_dust, + ! dust_mmr_ref diagnostics) but NOT fed to radiation - + ! useful for isolating dust's radiative effect from its + ! mere presence in diagnostics/other schemes. + ! Mirrors the existing account_for_effect_of_water/_ozone pattern above, which exists + ! because water/ozone are always available as fields; do_dust_forcing exists as well + ! because, unlike water/ozone, the dust field has no other reason to be computed. + ! Defaults to .FALSE. (unlike water/ozone above) so that setting do_dust_forcing alone + ! is never enough to silently turn on dust radiative effects - both flags must be set. + logical :: account_for_effect_of_dust=.FALSE. !if False then radiation is fed dust mixing ratios = 0. If true it's fed mixing ratios based on model dust field (see comment above). logical :: do_read_ozone = .FALSE. ! read ozone from an external file? character(len=256) :: ozone_file_name='ozone' !Name of file containing ozone field - n.b. don't need to include '.nc' character(len=256) :: ozone_field_name='ozone' !Name of ozone variable in ozone file logical :: input_o3_file_is_mmr=.true. ! Does the ozone input file contain values as a mass mixing ratio (set to true) or a volume mixing ratio (set to false)? logical :: do_read_co2 = .FALSE. ! read ozone from an external file? character(len=256) :: co2_file_name='co2' !Name of file containing co2 field - n.b. don't need to include '.nc' - character(len=256) :: co2_field_name='co2' !Name of co2 variable in co2 file + character(len=256) :: co2_field_name='co2' !Name of co2 variable in co2 file + logical :: do_read_cdod = .FALSE. ! read dust optical depth from external file? + character(len=256) :: cdod_file_name='dust' + character(len=256) :: cdod_field_name='cdod' !Name of dust optical depth variable in cdod file + logical :: do_dust_forcing=.FALSE. ! Master switch: compute a Mars dust vertical profile (Ball et al. 2021)? + ! See account_for_effect_of_dust above for how this interacts with whether + ! that profile actually reaches the radiation calculation. real(r_def) :: input_planet_emissivity = 1.0 !Emissivity of surface. Defined as constant all over surface. real :: co2_ppmv = 300. !Default CO2 concentration in PPMV logical :: input_co2_mmr=.false. !Socrates wants input concentrations as mmr not vmr, so need to make sure input data supplied is converted if necessary @@ -43,8 +67,8 @@ module socrates_config_mod logical :: do_scm_ozone=.false. ! read single column ozone from namelist? note: ONLY when using SCM real(r_def), dimension(100) :: scm_ozone = -1 ! input array for single column ozone. max number of levels = 100 - logical :: use_pressure_interp_for_half_levels = .False. !By default (.False.) does linear interpolation in height for half-level temperatures. True does linear interp using pressure. - + logical :: use_pressure_interp_for_half_levels = .False. !By default (.False.) does linear interpolation in height for half-level temperatures. True does linear interp using pressure. + real(r_def) :: dust_scale = 2.6e-4 ! Incoming radiation options for namelist integer :: solday=0 ! if >0, do perpetual run corresponding to day of the year = solday \in [0,days per year] @@ -71,6 +95,9 @@ module socrates_config_mod REAL(r_def) :: cfc113_mix_ratio = 4.826e-10 REAL(r_def) :: hcfc22_mix_ratio = 6.866e-10 REAL(r_def) :: hfc134a_mix_ratio = 2.536e-10 + + REAL(r_def) :: dust_mix_ratio = 1.0e-06 + REAL(r_def) :: nu_dust = 0.007 ! Whether to include radiative effects of particular gases @@ -121,9 +148,10 @@ module socrates_config_mod NAMELIST/socrates_rad_nml/ stellar_constant, tidally_locked, lw_spectral_filename, lw_hires_spectral_filename, & sw_spectral_filename, sw_hires_spectral_filename, socrates_hires_mode, & input_planet_emissivity, co2_ppmv, & - account_for_effect_of_water, account_for_effect_of_ozone, & + account_for_effect_of_water, account_for_effect_of_ozone, account_for_effect_of_dust, & do_read_ozone, ozone_file_name, ozone_field_name, input_o3_file_is_mmr, & - do_read_co2, co2_file_name, co2_field_name, input_co2_mmr, & + do_read_co2, co2_file_name, co2_field_name, input_co2_mmr, & + do_read_cdod, cdod_file_name, cdod_field_name, do_dust_forcing, & solday, do_rad_time_avg, equinox_day, & store_intermediate_rad, dt_rad_avg, dt_rad, & chunk_size, & @@ -134,6 +162,7 @@ module socrates_config_mod inc_h2o, inc_co2, inc_co, inc_o3, inc_n2o, inc_n2, inc_ch4, inc_o2, & inc_so2, inc_cfc11, inc_cfc12, inc_cfc113, inc_hcfc22, inc_hfc134a, & use_pressure_interp_for_half_levels, & - frierson_solar_rad, del_sol, del_sw, do_scm_ozone, scm_ozone + frierson_solar_rad, del_sol, del_sw, do_scm_ozone, scm_ozone, & + dust_mix_ratio, nu_dust, dust_scale end module socrates_config_mod diff --git a/src/atmos_param/socrates/interface/socrates_interface.F90 b/src/atmos_param/socrates/interface/socrates_interface.F90 index 0127d4fc6..88381fd80 100644 --- a/src/atmos_param/socrates/interface/socrates_interface.F90 +++ b/src/atmos_param/socrates/interface/socrates_interface.F90 @@ -63,7 +63,8 @@ MODULE socrates_interface_mod INTEGER :: id_soc_olr_clr, id_soc_toa_sw_clr, id_soc_toa_sw_up_clr ! clear-sky case INTEGER :: id_soc_surf_flux_sw_clr, id_soc_surf_flux_sw_down_clr, & id_soc_surf_flux_lw_clr, id_soc_surf_flux_lw_down_clr - INTEGER :: id_soc_ozone, id_soc_co2, id_soc_coszen + INTEGER :: id_soc_ozone, id_soc_co2, id_soc_cdod, id_soc_dust_mmr_ref, & + id_soc_dust, id_soc_coszen INTEGER :: n_soc_bands_lw, n_soc_bands_sw INTEGER :: n_soc_bands_lw_hires, n_soc_bands_sw_hires INTEGER :: id_soc_bins_lw, id_soc_bins_sw @@ -75,19 +76,21 @@ MODULE socrates_interface_mod CHARACTER(len=10), PARAMETER :: soc_mod_name = 'socrates' REAL :: missing_value = -999 - type(interpolate_type),save :: o3_interp, co2_interp ! use external file for ozone and co2 + type(interpolate_type),save :: o3_interp, co2_interp + type(interpolate_type),save :: cdod_interp ! use external file for dust optical depth REAL :: dt_last !Time of last radiation calculation - used to tell whether it is time to recompute radiation or not REAL(r_def), allocatable, dimension(:,:,:) :: tdt_soc_sw_store, tdt_soc_lw_store REAL(r_def), allocatable, dimension(:,:,:) :: thd_sw_flux_net_store, thd_lw_flux_net_store, & thd_sw_flux_clr_net_store, thd_lw_flux_clr_net_store - REAL(r_def), allocatable, dimension(:,:,:) :: thd_co2_store, thd_ozone_store + REAL(r_def), allocatable, dimension(:,:,:) :: thd_co2_store, thd_ozone_store, thd_dust_store REAL(r_def), allocatable, dimension(:,:) :: net_surf_sw_down_store, surf_lw_down_store, surf_lw_net_store, & surf_sw_down_store, toa_sw_down_store, & toa_sw_store, olr_store, coszen_store, & toa_sw_clr_store, olr_clr_store, toa_sw_up_store, toa_sw_up_clr_store, & net_surf_sw_down_clr_store, surf_lw_down_clr_store, & - surf_lw_net_clr_store, surf_sw_down_clr_store + surf_lw_net_clr_store, surf_sw_down_clr_store, & + thd_cdod_store, thd_dust_mmr_ref_store REAL(r_def), allocatable, dimension(:,:,:) :: outputted_soc_spectral_olr, spectral_olr_store, outputted_soc_spectral_olr_clr REAL(r_def), allocatable, dimension(:) :: soc_bins_lw, soc_bins_sw @@ -385,6 +388,21 @@ SUBROUTINE socrates_init(is, ie, js, je, num_levels, axes, Time, lat, lonb, latb 'socrates Co2', & 'mmr', missing_value=missing_value ) + id_soc_dust_mmr_ref = & + register_diag_field ( soc_mod_name, 'dust_mmr_ref', axes(1:2), Time, & + 'socrates dust_mmr_ref', & + 'mmr', missing_value=missing_value ) + + id_soc_cdod = & + register_diag_field ( soc_mod_name, 'soc_cdod', axes(1:2), Time, & + 'socrates CDOD', & + 'none', missing_value=missing_value ) + + id_soc_dust = & + register_diag_field ( soc_mod_name, 'soc_dust', axes(1:3), Time, & + 'socrates Dust', & + 'mmr', missing_value=missing_value ) + id_mars_solar_long = register_diag_field ( soc_mod_name, 'mars_solar_long', & Time, 'Martian solar longitude', 'deg') @@ -421,6 +439,11 @@ SUBROUTINE socrates_init(is, ie, js, je, num_levels, axes, Time, lat, lonb, latb WARNING) endif + if(do_read_cdod)then + call interpolator_init (cdod_interp, trim(cdod_file_name)//'.nc', lonb, latb, data_out_of_bounds=(/ZERO/)) + print*, 'CDOD interpolator initialized' + endif + if (mod((size(lonb,1)-1)*(size(latb,1)-1), chunk_size) .ne. 0) then call error_mesg( 'socrates_init', & 'chunk_size must equally divide number of points per processor, which it currently does not.', FATAL) @@ -528,6 +551,18 @@ SUBROUTINE socrates_init(is, ie, js, je, num_levels, axes, Time, lat, lonb, latb allocate(thd_co2_store(size(lonb,1)-1, size(latb,2)-1, num_levels)) endif + if (id_soc_cdod > 0 ) then + allocate(thd_cdod_store(size(lonb,1)-1, size(latb,2)-1)) + endif + + if (id_soc_dust_mmr_ref > 0) then + allocate(thd_dust_mmr_ref_store(size(lonb,1)-1, size(latb,2)-1)) + endif + + if (id_soc_dust > 0) then + allocate(thd_dust_store(size(lonb,1)-1, size(latb,2)-1, num_levels)) + endif + ! spectral output currently not available as required axis not present in diag file if (id_soc_spectral_olr > 0) then if (socrates_hires_mode .eqv. .True.) then @@ -590,7 +625,7 @@ end subroutine socrates_init ! Set up the call to the Socrates radiation scheme ! ----------------------------------------------------------------------------- subroutine socrates_interface(Time_diag, rlat, rlon, soc_lw_mode, & - fms_temp, fms_spec_hum, fms_ozone, fms_co2, fms_t_surf, & + fms_temp, fms_spec_hum, fms_ozone, fms_co2, fms_dust, fms_t_surf, & fms_p_full, fms_p_half, fms_z_full, fms_z_half, fms_albedo, & fms_coszen, fms_rrsun, n_profile, n_layer, & fms_cld_frac, fms_reff_rad, fms_mmr_cl_rad, & @@ -628,7 +663,7 @@ subroutine socrates_interface(Time_diag, rlat, rlon, soc_lw_mode, & INTEGER(i_def) :: nlat ! Input arrays - real(r_def), intent(in) :: fms_temp(:,:,:), fms_spec_hum(:,:,:), fms_ozone(:,:,:), fms_co2(:,:,:) + real(r_def), intent(in) :: fms_temp(:,:,:), fms_spec_hum(:,:,:), fms_ozone(:,:,:), fms_co2(:,:,:), fms_dust(:,:,:) real(r_def), intent(in) :: fms_p_full(:,:,:) real(r_def), intent(in) :: fms_p_half(:,:,:) real(r_def), intent(in) :: fms_t_surf(:,:), fms_albedo(:,:) @@ -666,6 +701,8 @@ subroutine socrates_interface(Time_diag, rlat, rlon, soc_lw_mode, & soc_heating_rate, input_o3_mixing_ratio, & input_co2_mixing_ratio,z_full_reshaped, input_cld_frac, input_reff_rad, input_mmr_cl_rad + real(r_def), dimension(n_profile, n_layer) :: input_dust_mixing_ratio + real(r_def), dimension(n_profile, 0:n_layer) :: input_p_level, input_t_level, soc_flux_direct, & soc_flux_down, soc_flux_up, soc_flux_direct_clr, soc_flux_down_clr, soc_flux_up_clr, z_half_reshaped @@ -728,6 +765,17 @@ subroutine socrates_interface(Time_diag, rlat, rlon, soc_lw_mode, & input_co2_mixing_ratio = reshape(fms_co2(:,:,:),(/si*sj,sk /)) + ! Second (radiative-transfer) gate on dust: fms_dust here is whatever run_socrates + ! computed, which is already zero if do_dust_forcing=.FALSE. (see its Conrath-profile + ! block below). account_for_effect_of_dust decides whether a non-zero profile is + ! actually passed on to Socrates, or zeroed right here - see the comment by both + ! flags' declarations in socrates_config_mod.f90 for the full picture. + if (account_for_effect_of_dust .eqv. .true.) then + input_dust_mixing_ratio = reshape(fms_dust(:,:,:),(/si*sj,sk /)) + else + input_dust_mixing_ratio = 0.0 + endif + !------------- !Default parameters @@ -830,6 +878,7 @@ subroutine socrates_interface(Time_diag, rlat, rlon, soc_lw_mode, & input_mixing_ratio(idx_chunk_start:idx_chunk_end,:), & input_o3_mixing_ratio(idx_chunk_start:idx_chunk_end,:), & input_co2_mixing_ratio(idx_chunk_start:idx_chunk_end,:), & + input_dust_mixing_ratio(idx_chunk_start:idx_chunk_end,:), & input_t_surf(idx_chunk_start:idx_chunk_end), & input_cos_zenith_angle(idx_chunk_start:idx_chunk_end), & input_solar_irrad(idx_chunk_start:idx_chunk_end), & @@ -863,6 +912,7 @@ subroutine socrates_interface(Time_diag, rlat, rlon, soc_lw_mode, & input_mixing_ratio(idx_chunk_start:idx_chunk_end,:), & input_o3_mixing_ratio(idx_chunk_start:idx_chunk_end,:), & input_co2_mixing_ratio(idx_chunk_start:idx_chunk_end,:), & + input_dust_mixing_ratio(idx_chunk_start:idx_chunk_end,:), & input_t_surf(idx_chunk_start:idx_chunk_end), & input_cos_zenith_angle(idx_chunk_start:idx_chunk_end), & input_solar_irrad(idx_chunk_start:idx_chunk_end), & @@ -933,10 +983,13 @@ subroutine run_socrates(Time, Time_diag, rad_lat, rad_lon, temp_in, q_in, t_surf real, intent(in), dimension(:,:,:) :: cf_rad, reff_rad, qcl_rad integer(i_def) :: n_profile, n_layer + ! loop variables + integer :: i, j, k - real(r_def), dimension(size(temp_in,1), size(temp_in,2)) :: t_surf_for_soc, rad_lat_soc, rad_lon_soc, albedo_soc + real(r_def), dimension(size(temp_in,1), size(temp_in,2)) :: t_surf_for_soc, rad_lat_soc, rad_lon_soc, albedo_soc, sin_lat, zmax + real(r_def), dimension(size(temp_in,1), size(temp_in,2)) :: dust_mmr_ref, cdod_in real(r_def), dimension(size(temp_in,1), & - size(temp_in,2), size(temp_in,3)) :: tg_tmp_soc, q_soc, ozone_soc, co2_soc, p_full_soc, & + size(temp_in,2), size(temp_in,3)) :: tg_tmp_soc, q_soc, ozone_soc, co2_soc, dust_soc, p_full_soc, & output_heating_rate_sw, output_heating_rate_lw, output_heating_rate_total, & output_heating_rate_sw_clr, output_heating_rate_lw_clr, & z_full_soc, cld_frac_soc, reff_rad_soc, mmr_cl_rad_soc, qcl_rad_soc, & @@ -957,7 +1010,7 @@ subroutine run_socrates(Time, Time_diag, rad_lat, rad_lon, temp_in, q_in, t_surf p2, toa_sw_down, surf_sw_down, & olr_clr, toa_sw_clr, toa_sw_up, toa_sw_up_clr, & net_surf_sw_down_clr, surf_sw_down_clr, surf_lw_net_clr, surf_lw_down_clr - real, dimension(size(temp_in,1), size(temp_in,2), size(temp_in,3)) :: ozone_in, co2_in + real, dimension(size(temp_in,1), size(temp_in,2), size(temp_in,3)) :: ozone_in, co2_in, dust_in real, dimension(size(temp_in,1), size(temp_in,2), size(temp_in,3)+1) :: thd_sw_flux_net, thd_lw_flux_net, thd_sw_flux_clr_net, thd_lw_flux_clr_net type(time_type) :: Time_loc @@ -1059,6 +1112,18 @@ subroutine run_socrates(Time, Time_diag, rad_lat, rad_lon, temp_in, q_in, t_surf co2_in = thd_co2_store endif + if (id_soc_cdod > 0) then + cdod_in = thd_cdod_store + endif + + if (id_soc_dust_mmr_ref > 0) then + dust_mmr_ref = thd_dust_mmr_ref_store + endif + + if (id_soc_dust > 0) then + dust_in = thd_dust_store + endif + if (id_soc_spectral_olr > 0) then outputted_soc_spectral_olr = spectral_olr_store endif @@ -1068,19 +1133,19 @@ subroutine run_socrates(Time, Time_diag, rad_lat, rad_lon, temp_in, q_in, t_surf endif if (id_time_since_ae > 0) then time_since_ae = time_since_ae_store - endif + endif if (id_dec > 0) then dec = dec_store - endif + endif if (id_ang > 0) then - ang_out= ang_out_store - endif + ang_out= ang_out_store + endif if (id_true_anom > 0) then true_anomaly = true_anomaly_store - endif + endif if (id_rrsun > 0) then rrsun = rrsun_store - endif + endif else output_heating_rate_sw = 0. @@ -1107,13 +1172,16 @@ subroutine run_socrates(Time, Time_diag, rad_lat, rad_lon, temp_in, q_in, t_surf coszen = 0. ozone_in = 0. co2_in = 0. + cdod_in = 0. + dust_mmr_ref = 0. + dust_in = 0. outputted_soc_spectral_olr = 0. mars_solar_long = 0. time_since_ae = 0. dec = 0. ang_out = 0. true_anomaly = 0. - rrsun = 0. + rrsun = 0. endif temp_tend(:,:,:) = temp_tend(:,:,:) + real(output_heating_rate_sw)+real(output_heating_rate_lw) @@ -1198,28 +1266,37 @@ subroutine run_socrates(Time, Time_diag, rad_lat, rad_lon, temp_in, q_in, t_surf if(id_soc_ozone > 0) then used = send_data ( id_soc_ozone, ozone_in, Time_diag) endif + if(id_soc_cdod > 0) then + used = send_data ( id_soc_cdod, cdod_in, Time_diag) + endif + if(id_soc_dust_mmr_ref > 0) then + used = send_data ( id_soc_dust_mmr_ref, dust_mmr_ref, Time_diag) + endif + if(id_soc_dust > 0) then + used = send_data ( id_soc_dust, dust_in, Time_diag) + endif if(id_soc_spectral_olr > 0) then used = send_data ( id_soc_spectral_olr, outputted_soc_spectral_olr, Time_diag) - endif + endif if (id_mars_solar_long > 0) then used = send_data ( id_mars_solar_long, mars_solar_long, Time_diag) endif if (id_time_since_ae > 0) then used = send_data ( id_time_since_ae, time_since_ae, Time_diag) - endif + endif if (id_dec > 0) then used = send_data ( id_dec, dec, Time_diag) - endif + endif if (id_ang > 0) then used = send_data ( id_ang, ang_out, Time_diag) - endif + endif if (id_true_anom > 0) then used = send_data ( id_true_anom, true_anomaly, Time_diag) - endif + endif if (id_rrsun > 0) then - used = send_data ( id_rrsun, rrsun, Time_diag) - endif - ! Diagnostics sent + used = send_data ( id_rrsun, rrsun, Time_diag) + endif + ! Diagnostics sent return !not time yet @@ -1282,33 +1359,33 @@ subroutine run_socrates(Time, Time_diag, rad_lat, rad_lon, temp_in, q_in, t_surf ozone_in = 0.0 - !get ozone + !get ozone if(do_read_ozone)then call interpolator( o3_interp, Time_diag, p_half_in, ozone_in, trim(ozone_field_name)) - endif - if(do_scm_ozone)then ! Allows for option to specify ozone vertical profile in namelist for SCM. - if(do_read_ozone)then + endif + if(do_scm_ozone)then ! Allows for option to specify ozone vertical profile in namelist for SCM. + if(do_read_ozone)then call error_mesg('socrates_interface', 'Cannot set do_scm_ozone and do_read_ozone = .true.', FATAL) - endif - if((size(temp_in,1)>1).or.(size(temp_in,2)>1))then + endif + if((size(temp_in,1)>1).or.(size(temp_in,2)>1))then call error_mesg('socrates_interface', 'Cannot set do_scm_ozone if simulating more than one column, use do_read_ozone instead', FATAL) - endif - if(scm_ozone(size(temp_in,3)).eq.-1)then + endif + if(scm_ozone(size(temp_in,3)).eq.-1)then call error_mesg('socrates_interface', 'Input o3 must be specified on model pressure levels but not enough levels specified', FATAL) - endif - if(scm_ozone(size(temp_in,3)+1).ne.-1)then + endif + if(scm_ozone(size(temp_in,3)+1).ne.-1)then call error_mesg('socrates_interface', 'Input o3 must be specified on model pressure levels but too many levels specified', FATAL) - endif + endif ozone_in(1,1,:) = scm_ozone(1:size(temp_in,3)) - !PUT THIS WARNING SOMEWHERE ELSE + !PUT THIS WARNING SOMEWHERE ELSE endif if (do_read_ozone .or. do_scm_ozone) then if (input_o3_file_is_mmr.eqv..false.) then ozone_in = ozone_in * wtmozone / (1000. * gas_constant / rdgas ) !Socrates expects all abundances to be mass mixing ratio. So if input file is volume mixing ratio, it must be converted to mass mixing ratio using the molar masses of dry air and ozone ! Molar mass of dry air calculated from gas_constant / rdgas, and converted into g/mol from kg/mol by multiplying by 1000. This conversion is necessary because wtmozone is in g/mol. - - endif + + endif endif if (input_co2_mmr .eqv. .false.) then @@ -1330,6 +1407,35 @@ subroutine run_socrates(Time, Time_diag, rad_lat, rad_lon, temp_in, q_in, t_surf endif endif + !get dust optical depth / mass mixing ratio reference value + if(do_read_cdod)then + call interpolator( cdod_interp, Time_diag, cdod_in, trim(cdod_field_name) ) + dust_mmr_ref = dust_scale*cdod_in / (16.4 - dust_scale*cdod_in) ! Converts dust optical depth at 610Pa to dust mass mixing ratio at 610Pa. + else + dust_mmr_ref = dust_mix_ratio + endif + + ! Build the dust vertical profile (Conrath-type distribution) from the reference + ! mixing ratio computed above, following Ball et al. 2021's Mars dust representation. + ! This is the first (computation) gate on dust: dust_in ends up all zero here if + ! do_dust_forcing=.FALSE., regardless of account_for_effect_of_dust. If this profile + ! IS computed, whether it's then actually radiatively active is decided separately, + ! further down the call chain in socrates_interface (see the comment there, and by + ! both flags' declarations in socrates_config_mod.f90). + if (do_dust_forcing .eqv. .true.) then + sin_lat(:,:) = sin(rad_lat(:,:)) + zmax(:,:) = 60 + 18*sin((mars_solar_long-158.)*pi/180.) & + -(32+18*sin((mars_solar_long-158.)*pi/180.))*(sin_lat(:,:))**4 & + -8*sin((mars_solar_long-158.)*pi/180.)*(sin_lat(:,:))**5 + do i=1, size(temp_in,1) + do j=1, size(temp_in,2) + dust_in(i, j, :) = dust_mmr_ref(i, j)*exp(nu_dust*(1-MAX((700/p_full_in(i,j,:))**(70./zmax(i,j)),1.))) + end do + end do + else + dust_in = 0. + endif + if(do_cloud_simple .or. do_cloud_spookie) then cld_frac_soc = REAL(cf_rad, kind(r_def)) reff_rad_soc = REAL(reff_rad, kind(r_def)) @@ -1359,6 +1465,7 @@ subroutine run_socrates(Time, Time_diag, rad_lat, rad_lon, temp_in, q_in, t_surf q_soc = REAL(q_in, kind(r_def)) ozone_soc = REAL(ozone_in, kind(r_def)) co2_soc = REAL(co2_in, kind(r_def)) + dust_soc = REAL(dust_in, kind(r_def)) p_full_soc = REAL(p_full_in, kind(r_def)) p_half_soc = REAL(p_half_in, kind(r_def)) albedo_soc = REAL(albedo_in, kind(r_def)) @@ -1366,7 +1473,7 @@ subroutine run_socrates(Time, Time_diag, rad_lat, rad_lon, temp_in, q_in, t_surf z_half_soc = REAL(z_half_in, kind(r_def)) CALL socrates_interface(Time, rad_lat_soc, rad_lon_soc, soc_lw_mode, & - tg_tmp_soc, q_soc, ozone_soc, co2_soc, t_surf_for_soc, p_full_soc, & + tg_tmp_soc, q_soc, ozone_soc, co2_soc, dust_soc, t_surf_for_soc, p_full_soc, & p_half_soc, z_full_soc, z_half_soc, albedo_soc, coszen, rrsun, & n_profile, n_layer, cld_frac_soc, reff_rad_soc, mmr_cl_rad_soc, & output_heating_rate_lw, output_soc_flux_lw_down, output_soc_flux_lw_up, & @@ -1394,7 +1501,7 @@ subroutine run_socrates(Time, Time_diag, rad_lat, rad_lon, temp_in, q_in, t_surf ! Retrieve output_heating_rate, and downward surface SW and LW fluxes soc_lw_mode = .FALSE. CALL socrates_interface(Time, rad_lat_soc, rad_lon_soc, soc_lw_mode, & - tg_tmp_soc, q_soc, ozone_soc, co2_soc, t_surf_for_soc, p_full_soc, & + tg_tmp_soc, q_soc, ozone_soc, co2_soc, dust_soc, t_surf_for_soc, p_full_soc, & p_half_soc, z_full_soc, z_half_soc, albedo_soc, coszen, rrsun, & n_profile, n_layer, cld_frac_soc, reff_rad_soc, mmr_cl_rad_soc, & output_heating_rate_sw, output_soc_flux_sw_down, output_soc_flux_sw_up, & @@ -1506,6 +1613,18 @@ subroutine run_socrates(Time, Time_diag, rad_lat, rad_lon, temp_in, q_in, t_surf thd_co2_store = co2_in endif + if (id_soc_cdod > 0) then + thd_cdod_store = cdod_in + endif + + if (id_soc_dust_mmr_ref > 0) then + thd_dust_mmr_ref_store = dust_mmr_ref + endif + + if (id_soc_dust > 0) then + thd_dust_store = dust_in + endif + if (id_soc_spectral_olr > 0) then spectral_olr_store = outputted_soc_spectral_olr endif @@ -1515,19 +1634,19 @@ subroutine run_socrates(Time, Time_diag, rad_lat, rad_lon, temp_in, q_in, t_surf endif if (id_time_since_ae > 0) then time_since_ae_store = time_since_ae - endif + endif if (id_dec > 0) then dec_store = dec - endif + endif if (id_ang > 0) then ang_out_store = ang_out - endif + endif if (id_true_anom > 0) then true_anomaly_store = true_anomaly - endif + endif if (id_rrsun > 0) then rrsun_store = rrsun - endif + endif endif @@ -1608,6 +1727,15 @@ subroutine run_socrates(Time, Time_diag, rad_lat, rad_lon, temp_in, q_in, t_surf if(id_soc_ozone > 0) then used = send_data ( id_soc_ozone, ozone_in, Time_diag) endif + if(id_soc_cdod > 0) then + used = send_data ( id_soc_cdod, cdod_in, Time_diag) + endif + if(id_soc_dust_mmr_ref > 0) then + used = send_data ( id_soc_dust_mmr_ref, dust_mmr_ref, Time_diag) + endif + if(id_soc_dust > 0) then + used = send_data ( id_soc_dust, dust_in, Time_diag) + endif if(id_soc_spectral_olr > 0) then used = send_data ( id_soc_spectral_olr, outputted_soc_spectral_olr, Time_diag) endif @@ -1619,19 +1747,19 @@ subroutine run_socrates(Time, Time_diag, rad_lat, rad_lon, temp_in, q_in, t_surf endif if (id_time_since_ae > 0) then used = send_data ( id_time_since_ae, time_since_ae, Time_diag) - endif + endif if (id_dec > 0) then used = send_data ( id_dec, dec, Time_diag) - endif + endif if (id_ang > 0) then used = send_data ( id_ang, ang_out, Time_diag) - endif + endif if (id_true_anom > 0) then used = send_data ( id_true_anom, true_anomaly, Time_diag) - endif + endif if (id_rrsun > 0) then - used = send_data ( id_rrsun, rrsun, Time_diag) - endif + used = send_data ( id_rrsun, rrsun, Time_diag) + endif ! Diagnostics sent @@ -1644,6 +1772,7 @@ subroutine run_socrates_end if(do_read_ozone) call interpolator_end(o3_interp) if(do_read_co2) call interpolator_end(co2_interp) + if(do_read_cdod) call interpolator_end(cdod_interp) end subroutine run_socrates_end diff --git a/src/atmos_param/socrates/interface/socrates_set_aer.F90 b/src/atmos_param/socrates/interface/socrates_set_aer.F90 new file mode 100644 index 000000000..66403a01a --- /dev/null +++ b/src/atmos_param/socrates/interface/socrates_set_aer.F90 @@ -0,0 +1,71 @@ +! *****************************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******************************* +MODULE socrates_set_aer +IMPLICIT NONE +CONTAINS + +! Subroutine to set the aerosol fields for the core radiation code. +!------------------------------------------------------------------------------ +SUBROUTINE set_aer(control, dimen, spectrum, aer, n_profile, n_layer, dust) + +USE rad_pcf, ONLY: ip_aersrc_classic_ron +USE def_spectrum, ONLY: StrSpecData +USE def_dimen, ONLY: StrDim +USE def_control, ONLY: StrCtrl +USE def_aer, ONLY: StrAer, allocate_aer, allocate_aer_prsc +use soc_constants_mod, only: i_def, r_def + +IMPLICIT NONE + +integer :: i_aer, i, l + +! Control options: +TYPE(StrCtrl), INTENT(IN) :: control + +! Dimensions: +TYPE(StrDim), INTENT(IN) :: dimen + +! Spectral data: +TYPE (StrSpecData), INTENT(IN) :: spectrum + +! Aerosol properties: +TYPE(StrAer), INTENT(OUT) :: aer + +INTEGER(i_def), INTENT(IN) :: n_profile +! Number of atmospheric profiles for radiation calculations +INTEGER(i_def), INTENT(IN) :: n_layer +! Number of atmospheric layers for radiation calculations + +REAL(r_def), INTENT(IN) :: dust(n_profile, n_layer) +! Dust mass mixing ratio + +! Allocate structure for the core radiation code interface +CALL allocate_aer(aer, dimen, spectrum) +CALL allocate_aer_prsc(aer, dimen, spectrum) + +! Only touch mr_source/mix_ratio when aerosol radiative effects are actually +! switched on (control%l_aerosol, set from do_dust_forcing in read_control.F90) +! - for ordinary runs aer%mr_source must be left exactly as +! allocate_aer/allocate_aer_prsc set it, matching the upstream set_aer_mod's +! behaviour. This can NOT be gated on spectrum%aerosol%n_aerosol > 0: the +! standard ga7 spectral files also define aerosol species blocks (inherited +! from the Met Office UM) despite Isca never wiring up their parametrizations, +! so that check is true even for non-Mars, non-dust runs. +IF (control%l_aerosol) THEN + aer%mr_source = ip_aersrc_classic_ron + + DO i_aer=1, spectrum%aerosol%n_aerosol + aer%mr_type_index(i_aer)=i_aer + DO i=1, n_layer + DO l=1, n_profile + aer%mix_ratio(l,i,i_aer) = dust(l,i) + END DO + END DO + END DO +END IF + +END SUBROUTINE set_aer +END MODULE socrates_set_aer diff --git a/src/atmos_shared/interpolator/interpolator.F90 b/src/atmos_shared/interpolator/interpolator.F90 index de9a2e9db..d40d46467 100644 --- a/src/atmos_shared/interpolator/interpolator.F90 +++ b/src/atmos_shared/interpolator/interpolator.F90 @@ -78,6 +78,7 @@ module interpolator_mod get_calendar_type, & JULIAN, NOLEAP, & THIRTY_DAY_MONTHS, & !mj + NO_CALENDAR, & get_date_julian, set_date_no_leap, & set_date_julian, get_date_no_leap, & print_date, & @@ -147,6 +148,7 @@ module interpolator_mod integer :: vertical_indices ! direction of vertical ! data axis logical :: climatological_year ! Is data for year = 0000? +logical :: mars !Field specific data for nfields type(fieldtype), pointer :: field_type(:) =>NULL() ! NetCDF field type @@ -261,6 +263,7 @@ subroutine interpolate_type_eq (Out, In) Out%je = In%je Out%vertical_indices = In%vertical_indices Out%climatological_year = In%climatological_year + Out%mars = In%mars Out%field_type => In%field_type if (associated(In%field_name )) Out%field_name => In%field_name if (associated(In%time_init )) Out%time_init => In%time_init @@ -542,8 +545,10 @@ subroutine interpolator_init( clim_type, file_name, lonb_mod, latb_mod, & call mpp_error(FATAL,'Interpolator_init : Time units not recognised in file '//file_name) end select - clim_type%climatological_year = (fileyr == 0) - if (.not. clim_type%climatological_year) then + clim_type%climatological_year = (fileyr == 0 .and. model_calendar/=NO_CALENDAR) + clim_type%mars = (fileyr == 0 .and. model_calendar==NO_CALENDAR) + + if (.not. clim_type%climatological_year .and. .not. clim_type%mars) then !---------------------------------------------------------------------- ! if file date has a non-zero year in the base time, determine that @@ -637,7 +642,10 @@ subroutine interpolator_init( clim_type, file_name, lonb_mod, latb_mod, & !! year, not the displacement from a base_time. clim_type%time_slice(n) = & set_time(INT( ( time_in(n) - INT(time_in(n)) ) * 86400 ),INT(time_in(n))) - else + else if (clim_type%mars) then + clim_type%time_slice(n) = & + set_time(INT( ( time_in(n) - INT(time_in(n)) ) * 88440 ),INT(time_in(n))) + else !-------------------------------------------------------------------- ! if fileyr /= 0 (i.e., climatological_year=F), diff --git a/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 b/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 index cfcbb4606..3a10c655a 100644 --- a/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 +++ b/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 @@ -26,7 +26,7 @@ module idealized_moist_phys_mod use mixed_layer_mod, only: mixed_layer_init, mixed_layer, mixed_layer_end, albedo_calc -use lscale_cond_mod, only: lscale_cond_init, lscale_cond, lscale_cond_end +use lscale_cond_mod, only: lscale_cond_init, lscale_cond, lscale_cond_co2_mars, lscale_cond_end use qe_moist_convection_mod, only: qe_moist_convection_init, qe_moist_convection, qe_moist_convection_end @@ -119,6 +119,7 @@ module idealized_moist_phys_mod logical :: do_ras = .false. logical :: do_lscale_cond = .true. +logical :: do_lscale_cond_lh = .false. ! Cloud options logical :: do_cloud_simple = .false. ! SimCloud cloud scheme logical :: do_cloud_spookie = .false. ! SPOOKIE protocol cloud scheme @@ -175,7 +176,7 @@ module idealized_moist_phys_mod max_bucket_depth_land, robert_bucket, raw_bucket, & do_lscale_cond, do_socrates_radiation, do_lcl_diffusivity_depth, damping_coeff_bucket, & finite_bucket_depth_over_land, & - do_local_heating + do_local_heating, do_lscale_cond_lh integer, parameter :: num_time_levels = 2 ! Add bucket - number of time levels added to allow timestepping in this module @@ -245,7 +246,10 @@ module idealized_moist_phys_mod conv_dt_tg, & ! temperature tendency from convection conv_dt_qg, & ! moisture tendency from convection cond_dt_tg, & ! temperature tendency from condensation - cond_dt_qg ! moisture tendency from condensation + cond_dt_qg, & ! moisture tendency from condensation + cond_lh_dt_tg, & + cond_lh_dt_qg, & + lh_rel logical, allocatable, dimension(:,:) :: & @@ -293,11 +297,14 @@ module idealized_moist_phys_mod id_conv_dt_qg, & ! temperature tendency from convection id_cond_dt_tg, & ! temperature tendency from condensation id_cond_dt_qg, & ! temperature tendency from condensation + id_cond_lh_dt_tg, & + id_cond_lh_dt_qg, & + id_lh_rel, & id_bucket_depth, & ! bucket depth variable for output id_bucket_depth_conv, & ! bucket depth variation induced by convection id_bucket_depth_cond, & ! bucket depth variation induced by condensation id_bucket_depth_lh, & ! bucket depth variation induced by LH - id_bucket_diffusion, & ! diffused surface water depth + id_bucket_diffusion, & ! diffused surface water depth id_rh, & ! Relative humidity id_diss_heat_ray,& ! Heat dissipated by rayleigh bottom drag if gp_surface=.True. id_z_tg, & ! Relative humidity @@ -557,6 +564,9 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l allocate(conv_dt_qg (is:ie, js:je, num_levels)) allocate(cond_dt_tg (is:ie, js:je, num_levels)) allocate(cond_dt_qg (is:ie, js:je, num_levels)) +allocate(cond_lh_dt_tg (is:ie, js:je, num_levels)) +allocate(cond_lh_dt_qg (is:ie, js:je, num_levels)) +allocate(lh_rel (is:ie, js:je, num_levels)) allocate(coldT (is:ie, js:je)); coldT = .false. allocate(klzbs (is:ie, js:je)); klzbs = 0 @@ -693,6 +703,12 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l axes(1:3), Time, 'Moisture tendency from condensation','kg/kg/s') id_cond_dt_tg = register_diag_field(mod_name, 'dt_tg_condensation', & axes(1:3), Time, 'Temperature tendency from condensation','K/s') +id_cond_lh_dt_qg = register_diag_field(mod_name, 'dt_qg_lh_condensation', & + axes(1:3), Time, 'Moisture tendency from condensation','kg/kg/s') +id_cond_lh_dt_tg = register_diag_field(mod_name, 'dt_tg_lh_condensation', & + axes(1:3), Time, 'Temperature tendency from condensation','K/s') +id_lh_rel = register_diag_field(mod_name, 'lh_rel', & + axes(1:3), Time, 'Latent heat released from condensation', 'K') id_cond_rain = register_diag_field(mod_name, 'condensation_rain', & axes(1:2), Time, 'Rain from condensation','kg/m/m/s') id_precip = register_diag_field(mod_name, 'precipitation', & @@ -821,6 +837,9 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l if (r_conv_scheme .eq. DRY_CONV .and. do_lscale_cond .eqv. .true.) then call error_mesg('idealized_moist_phys','do_lscale_cond is .true. but r_conv_scheme is dry. These options may not be consistent.', WARNING) endif +if (r_conv_scheme .eq. DRY_CONV .and. do_lscale_cond_lh .eqv. .true.) then + call error_mesg('idealized_moist_phys','do_lscale_cond_lh is .true. but r_conv_scheme is dry. These options may not be consistent.', WARNING) +endif if(two_stream_gray) call two_stream_gray_rad_init(is, ie, js, je, num_levels, get_axis_id(), Time, rad_lonb_2d, rad_latb_2d, dt_real) @@ -1033,6 +1052,35 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, ps dt_tg = dt_tg + conv_dt_tg dt_tracers(:,:,:,nsphum) = dt_tracers(:,:,:,nsphum) + conv_dt_qg +! Perform CO2 condensation with latent heating (Mars) +if ( do_lscale_cond_lh .eqv. .true.) then + ! Large scale convection is a function of humidity only. This is + ! inconsistent with the dry convection scheme, don't run it! + rain = 0.0; snow = 0.0; lh_rel = 0.0 + call lscale_cond_co2_mars (tg_tmp, qg_tmp, & + p_full(:,:,:,previous), p_half(:,:,:,previous), & + lh_rel, coldT, rain, & + snow, cond_lh_dt_tg, & + cond_lh_dt_qg ) + + cond_lh_dt_tg = cond_lh_dt_tg/delta_t + cond_lh_dt_qg = cond_lh_dt_qg/delta_t + depth_change_cond = rain/dens_vapor + rain = rain/delta_t + snow = snow/delta_t + precip = precip + rain + snow + + dt_tg = dt_tg + cond_lh_dt_tg + dt_tracers(:,:,:,nsphum) = dt_tracers(:,:,:,nsphum) + cond_lh_dt_qg + + if(id_cond_lh_dt_qg > 0) used = send_data(id_cond_lh_dt_qg, cond_lh_dt_qg, Time) + if(id_cond_lh_dt_tg > 0) used = send_data(id_cond_lh_dt_tg, cond_lh_dt_tg, Time) + if(id_cond_rain > 0) used = send_data(id_cond_rain, rain, Time) + if(id_precip > 0) used = send_data(id_precip, precip, Time) + if(id_lh_rel > 0) used = send_data(id_lh_rel, lh_rel, Time) + +endif + convective_rain = precip ! Perform large scale convection diff --git a/src/extra/model/socrates/path_names b/src/extra/model/socrates/path_names index cab81e74d..c860ab9f0 100644 --- a/src/extra/model/socrates/path_names +++ b/src/extra/model/socrates/path_names @@ -509,3 +509,4 @@ atmos_param/socrates/interface/set_control.F90 atmos_param/socrates/interface/socrates_calc.F90 atmos_param/socrates/interface/socrates_interface.F90 atmos_param/socrates/interface/soc_constants.f90 +atmos_param/socrates/interface/socrates_set_aer.F90 diff --git a/src/extra/python/isca/util.py b/src/extra/python/isca/util.py index 68696da2a..ecae6611a 100644 --- a/src/extra/python/isca/util.py +++ b/src/extra/python/isca/util.py @@ -162,6 +162,8 @@ def interpolate_output(infile, outfile, all_fields=True, var_names=[], p_levs = elif p_levs.upper() == "EVEN": #plev = ' -p "100000 95000 90000 85000 80000 75000 70000 65000 60000 55000 50000 45000 40000 35000 30000 25000 20000 15000 10000 5000" ' levels = [100000, 95000, 90000, 85000, 80000, 75000, 70000, 65000, 60000, 55000, 50000, 45000, 40000, 35000, 30000, 25000, 20000, 15000, 10000, 5000] + elif p_levs.upper() == "OTHER": + levels = [625, 615, 605, 595, 585, 575, 565, 555, 545, 535, 525, 515, 505, 495, 485, 475, 465, 455, 445, 435, 425, 420, 415, 410, 405, 400, 395, 390, 385, 380, 375, 350, 345, 340, 335, 330, 325, 320, 315, 310, 305, 300, 295, 290, 285, 280, 275, 250, 225, 200, 175, 150, 125, 100, 75, 50, 30, 20, 10, 1] else: raise ValueError("Unknown p_levs type '{}'".format(p_levs)) else: @@ -261,4 +263,4 @@ def run_cli(exp, fail_if_underconditioned=True): runs = [config['run']] for i in runs: with context(exp): - exp.run(i,**config['run_config']) \ No newline at end of file + exp.run(i,**config['run_config']) diff --git a/src/extra/python/requirements.txt b/src/extra/python/requirements.txt index 722b8d75b..4ceab4c42 100644 --- a/src/extra/python/requirements.txt +++ b/src/extra/python/requirements.txt @@ -5,3 +5,4 @@ numpy pandas xarray tqdm +netcdf4