Merging in Emily Ball's Isca-Mars developments

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: <https://github.com/emilyrball/Isca-Mars/tree/bp1_lh_dev/exp/socrates_mars/input>.
+
+## 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.

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)