diff --git a/exp/test_cases/cg_drag/cg_drag_qbo_test_case.py b/exp/test_cases/cg_drag/cg_drag_qbo_test_case.py new file mode 100644 index 000000000..fc9a93403 --- /dev/null +++ b/exp/test_cases/cg_drag/cg_drag_qbo_test_case.py @@ -0,0 +1,192 @@ +""" +Non-orographic (convective) gravity wave drag test case, using the cg_drag +scheme -- originally implemented by Chaim Garfinkel (based on Alexander & +Dunkerton, JAS 1999), with edits by Stephen Thomson and further development +by Ross Castle (see src/atmos_param/cg_drag/cg_drag.f90 for the parameter +history). + +Physics settings (do_rayleigh off, do_cg_drag on, tuned cg_drag_nml) are the +QBO-capable configuration: do_rayleigh is off so the two forcings don't fight +each other. NUM_MONTHS is a short comparison run -- long enough for the +cg_drag forcing to develop and show up in the zonal-mean wind, well short of +the multi-year integration a full QBO cycle would need. +""" +import os + +from isca import IscaCodeBase, DiagTable, Experiment, Namelist, GFDL_BASE + +NCORES = 16 +RESOLUTION = 'T42', 50 +NUM_MONTHS = 6 + +cb = IscaCodeBase.from_directory(GFDL_BASE) + +exp = Experiment('cg_drag_qbo_test', codebase=cb) +exp.clear_rundir() + +exp.inputfiles = [os.path.join(GFDL_BASE, 'input/rrtm_input_files/ozone_1990.nc')] + +diag = DiagTable() +diag.add_file('atmos_monthly', 30, 'days', time_units='days') +diag.add_field('dynamics', 'ps', time_avg=True) +diag.add_field('dynamics', 'bk') +diag.add_field('dynamics', 'pk') +diag.add_field('atmosphere', 'precipitation', time_avg=True) +diag.add_field('mixed_layer', 't_surf', time_avg=True) +diag.add_field('dynamics', 'sphum', 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('rrtm_radiation', 'co2', time_avg=True) +diag.add_field('damping', 'udt_cgwd', time_avg=True) # cg_drag zonal wind tendency +exp.diag_table = diag + +exp.namelist = namelist = Namelist({ + 'main_nml': { + 'days': 30, + 'hours': 0, + 'minutes': 0, + 'seconds': 0, + 'dt_atmos': 120, + 'current_date': [1, 1, 1, 0, 0, 0], + 'calendar': 'thirty_day' + }, + + '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': 3.21e-05, + 'two_stream_gray': False, # Use RRTM, not grey radiation + 'do_rrtm_radiation': True, + 'convection_scheme': 'FULL_BETTS_MILLER' + }, + + 'vert_turb_driver_nml': { + 'do_mellor_yamada': False, + 'do_diffusivity': True, + 'do_simple': True, + 'constant_gust': 0.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 + }, + + 'atmosphere_nml': { + 'idealized_moist_model': True + }, + + 'mixed_layer_nml': { + 'tconst': 285., + 'prescribe_initial_dist': True, + 'evaporation': True, + 'depth': 100., + 'albedo_value': 0.25, + 'do_qflux': False, + }, + + '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 + }, + + 'damping_driver_nml': { + 'do_rayleigh': False, # off, so cg_drag is the only source of GWD forcing + 'trayfric': -0.5, + 'sponge_pbottom': 50., + 'do_conserve_energy': True, + 'do_mg_drag': False, + 'do_cg_drag': True, + }, + + 'cg_drag_nml': { + 'Bt_0': 0.0043, + 'Bt_nh': 0.0014, + 'Bt_eq': 0.0043, + 'Bt_sh': 0.0014, + 'phi0n': 15., + 'phi0s': -15., + 'dphin': 10., + 'dphis': -10., + 'flag': 1, + 'Bw': 0.4, + 'Bn': 8.4, + 'cw': 35.0, + 'cwtropics': 35.0, + 'cn': 2.0, + 'kelvin_kludge': 1.0, + 'weighttop': 0.7, + 'weightminus1': 0.28, + 'weightminus2': 0.02, + 'source_level_pressure': 315.e+02, + 'damp_level_pressure': 0.85e+02, + 'cg_drag_freq': 21600 + }, + + 'rrtm_radiation_nml': { + 'do_read_ozone': True, + 'ozone_file': 'ozone_1990', + 'solr_cnst': 1360., + 'dt_rad': 3600, + 'do_read_co2': False, + 'co2ppmv': 300, + }, + + 'diag_manager_nml': { + 'mix_snapshot_average_fields': False + }, + + 'fms_nml': { + 'domains_stack_size': 600000 + }, + + 'fms_io_nml': { + 'threading_write': 'single', + 'fileset_write': 'single', + }, + + 'spectral_dynamics_nml': { + 'damping_order': 4, + 'water_correction_limit': 200.e2, + 'reference_sea_level_press': 1.0e5, + 'num_levels': 50, + 'valid_range_t': [100., 800.], + 'initial_sphum': [2.e-6], + 'vert_coord_option': 'uneven_sigma', + 'surf_res': 0.1, + 'scale_heights': 7.9, + 'exponent': 1.4, + 'robert_coeff': 0.03, + } +}) + +exp.set_resolution(*RESOLUTION) + +if __name__ == '__main__': + cb.compile() + exp.run(1, use_restart=False, num_cores=NCORES) + for i in range(2, NUM_MONTHS + 1): + exp.run(i, num_cores=NCORES) diff --git a/exp/test_cases/mg_drag/input/mg_drag.res.nc b/exp/test_cases/mg_drag/input/mg_drag.res.nc new file mode 100644 index 000000000..dbaf842f7 Binary files /dev/null and b/exp/test_cases/mg_drag/input/mg_drag.res.nc differ diff --git a/exp/test_cases/mg_drag/mg_drag_socrates_topo_test_case.py b/exp/test_cases/mg_drag/mg_drag_socrates_topo_test_case.py new file mode 100644 index 000000000..50a3774a7 --- /dev/null +++ b/exp/test_cases/mg_drag/mg_drag_socrates_topo_test_case.py @@ -0,0 +1,228 @@ +""" +Orographic (mountain) gravity wave drag test case, using the mg_drag scheme +from Ross Castle's realistic_mgdrag branch, built on the official Socrates +aquaplanet-with-topography test case +(exp/test_cases/socrates_test/socrates_aquaplanet_amip_with_topo.py). + +mg_drag_nml settings are Ross's own tuned values, except source_of_sgsmtn. +His own configuration used 'computed', which reads a raw high-resolution +elevation dataset via topography_nml and computes subgrid-topography +variance internally (see src/shared/topography/topography.F90); that raw +elevation dataset isn't included in this repo. + +This test case instead uses source_of_sgsmtn='input', reading +input/mg_drag.res.nc (variable 'ghprime') directly -- a real subgrid- +orography-standard-deviation field regridded from ERA5 (see +src/extra/python/scripts/regrid_era5_sdor_to_t42.py for the regridding +method and source dataset). +""" +import os + +from isca import SocratesCodeBase, DiagTable, Experiment, Namelist, GFDL_BASE + +NCORES = 16 +RESOLUTION = 'T42', 40 +NUM_MONTHS = 6 + +base_dir = os.path.dirname(os.path.realpath(__file__)) +mg_drag_ghprime_input = os.path.join(base_dir, 'input', 'mg_drag.res.nc') + +cb = SocratesCodeBase.from_directory(GFDL_BASE) + +exp = Experiment('mg_drag_socrates_topo_test', codebase=cb) +exp.clear_rundir() + +exp.inputfiles = [ + os.path.join(GFDL_BASE, 'input/rrtm_input_files/ozone_1990.nc'), + os.path.join(GFDL_BASE, 'exp/test_cases/realistic_continents/input/era-spectral7_T42_64x128.out.nc'), + os.path.join(GFDL_BASE, 'exp/test_cases/realistic_continents/input/sst_clim_amip.nc'), + os.path.join(GFDL_BASE, 'exp/test_cases/realistic_continents/input/siconc_clim_amip.nc'), + mg_drag_ghprime_input, +] + +diag = DiagTable() +diag.add_file('atmos_monthly', 30, 'days', time_units='days') +diag.add_field('dynamics', 'ps', time_avg=True) +diag.add_field('dynamics', 'bk') +diag.add_field('dynamics', 'pk') +diag.add_field('dynamics', 'zsurf') +diag.add_field('atmosphere', 'precipitation', time_avg=True) +diag.add_field('atmosphere', 'rh', time_avg=True) +diag.add_field('mixed_layer', 't_surf', time_avg=True) +diag.add_field('dynamics', 'sphum', 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('socrates', 'soc_olr', time_avg=True) +# mg_drag diagnostics +diag.add_field('damping', 'udt_gwd', time_avg=True) +diag.add_field('damping', 'vdt_gwd', time_avg=True) +diag.add_field('damping', 'taubx', time_avg=True) +diag.add_field('damping', 'tauby', time_avg=True) +diag.add_field('damping', 'taus', time_avg=True) +diag.add_field('damping', 'tdt_diss_gwd', time_avg=True) +diag.add_field('damping', 'sgsmtn', time_avg=True) +exp.diag_table = diag + +exp.namelist = namelist = Namelist({ + 'main_nml': { + 'days': 30, + 'hours': 0, + 'minutes': 0, + 'seconds': 0, + 'dt_atmos': 180, # finer than the plain topo test case's 600s, needed for mg_drag stability + 'current_date': [1, 1, 1, 0, 0, 0], + 'calendar': 'thirty_day' + }, + + 'socrates_rad_nml': { + 'stellar_constant': 1370., + 'lw_spectral_filename': os.path.join(GFDL_BASE, 'src/atmos_param/socrates/src/trunk/data/spectra/ga7/sp_lw_ga7'), + 'sw_spectral_filename': os.path.join(GFDL_BASE, 'src/atmos_param/socrates/src/trunk/data/spectra/ga7/sp_sw_ga7'), + 'do_read_ozone': True, + 'ozone_file_name': 'ozone_1990', + 'ozone_field_name': 'ozone_1990', + 'dt_rad': 3600, + 'store_intermediate_rad': True, + 'chunk_size': 16, + 'use_pressure_interp_for_half_levels': False, + 'tidally_locked': False, + 'solday': 90 + }, + + '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': 3.21e-05, + 'two_stream_gray': False, + 'do_socrates_radiation': True, + 'convection_scheme': 'SIMPLE_BETTS_MILLER', + 'do_cloud_simple': False, + 'land_option': 'input', + 'land_file_name': 'INPUT/era-spectral7_T42_64x128.out.nc', + 'land_roughness_prefactor': 10.0, + }, + + 'vert_turb_driver_nml': { + 'do_mellor_yamada': False, + 'do_diffusivity': True, + 'do_simple': True, + 'constant_gust': 0.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, + 'land_humidity_prefactor': 0.7, + }, + + 'atmosphere_nml': { + 'idealized_moist_model': True + }, + + 'mixed_layer_nml': { + 'tconst': 285., + 'prescribe_initial_dist': True, + 'evaporation': True, + 'land_option': 'input', + 'land_h_capacity_prefactor': 0.1, + 'albedo_value': 0.25, + 'land_albedo_prefactor': 1.3, + 'do_qflux': False, + 'do_read_sst': True, + 'do_sc_sst': True, + 'sst_file': 'sst_clim_amip', + 'specify_sst_over_ocean_only': True, + }, + + 'qe_moist_convection_nml': { + 'rhbm': 0.7, + 'Tmin': 160., + 'Tmax': 350. + }, + + 'lscale_cond_nml': { + 'do_simple': True, + 'do_evap': True + }, + + 'sat_vapor_pres_nml': { + 'do_simple': True, + 'construct_table_wrt_liq_and_ice': True + }, + + 'damping_driver_nml': { + 'do_rayleigh': True, + 'trayfric': -0.5, + 'sponge_pbottom': 150., + 'do_conserve_energy': True, + 'do_cg_drag': False, + 'do_mg_drag': True, + }, + + 'mg_drag_nml': { + 'do_netcdf_restart': True, + 'xl_mtn': 1.0e5, + 'gmax': 1.0, + 'acoef': 1.0, + 'rho': 1.13, + 'low_lev_frac': 0.23, + 'do_conserve_energy': True, + 'do_mcm_mg_drag': False, + 'source_of_sgsmtn': 'input', # reads INPUT/mg_drag.res.nc ('ghprime'), ERA5-derived -- see module docstring + 'flux_cut_level': 0.0 + }, + + 'diag_manager_nml': { + 'mix_snapshot_average_fields': False + }, + + 'fms_nml': { + 'domains_stack_size': 600000 + }, + + 'fms_io_nml': { + 'threading_write': 'single', + 'fileset_write': 'single', + }, + + 'spectral_dynamics_nml': { + 'damping_order': 4, + 'water_correction_limit': 200.e2, + 'reference_sea_level_press': 1.0e5, + 'num_levels': 40, + 'valid_range_t': [100., 800.], + 'initial_sphum': [2.e-6], + 'vert_coord_option': 'uneven_sigma', + 'surf_res': 0.2, + 'scale_heights': 11.0, + 'exponent': 7.0, + 'robert_coeff': 0.03, + 'ocean_topog_smoothing': 0.0 + }, + + 'spectral_init_cond_nml': { + 'topog_file_name': 'era-spectral7_T42_64x128.out.nc', + 'topography_option': 'input' + }, +}) + +exp.set_resolution(*RESOLUTION) + +if __name__ == '__main__': + cb.compile(debug=False) + exp.run(1, use_restart=False, num_cores=NCORES) + for i in range(2, NUM_MONTHS + 1): + 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..f7ecae1c0 100644 --- a/exp/test_cases/trip_test/trip_test_functions.py +++ b/exp/test_cases/trip_test/trip_test_functions.py @@ -80,6 +80,20 @@ def get_nml_diag(test_case_name): nml_out = exp_temp.namelist codebase_to_use = ColumnCodeBase + if 'cg_drag_qbo' in test_case_name: + sys.path.insert(0, os.path.join(GFDL_BASE, 'exp/test_cases/cg_drag/')) + from cg_drag_qbo_test_case import exp as exp_temp + input_files = exp_temp.inputfiles + nml_out = exp_temp.namelist + codebase_to_use = IscaCodeBase + + if 'mg_drag_socrates_topo' in test_case_name: + sys.path.insert(0, os.path.join(GFDL_BASE, 'exp/test_cases/mg_drag/')) + from mg_drag_socrates_topo_test_case import exp as exp_temp + input_files = exp_temp.inputfiles + nml_out = exp_temp.namelist + codebase_to_use = SocratesCodeBase + if 'giant_planet' in test_case_name: sys.path.insert(0, os.path.join(GFDL_BASE, 'exp/test_cases/giant_planet/')) from giant_planet_test_case import exp as exp_temp @@ -211,6 +225,8 @@ def list_all_test_cases_implemented_in_trip_test(): 'column_test', 'grey_mars', 'radiative_eq_mars', + 'cg_drag_qbo', + 'mg_drag_socrates_topo', #'socrates_mars', # requires Mars-specific Socrates spectral files not yet included in the repo - see exp/test_cases/socrates_mars/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 ] @@ -356,13 +372,21 @@ def conduct_comparison_on_test_case(base_commit, later_commit, test_case_name, r base_commit_dataset = xar.open_dataset(data_dir_dict[base_commit] +'/run0001/'+diag_file_entry+'.nc', decode_times=False) later_commit_dataset = xar.open_dataset(data_dir_dict[later_commit]+'/run0001/'+diag_file_entry+'.nc', decode_times=False) - diff = later_commit_dataset - base_commit_dataset - - #Check each of the output variables for differences - for var in diff.data_vars.keys(): - maxval = np.abs(diff[var]).max() + #Compare raw array values position-by-position rather than taking an + #xarray Dataset difference. Dataset subtraction aligns the two datasets + #by coordinate value first, which raises if a coordinate has duplicate + #values -- e.g. the single-column model's 'lonb' is two identical + #(near-zero placeholder) values, since a 1-column grid has no real + #longitude extent. Both datasets are the same test case run at the same + #resolution by construction (only the commit differs), so they already + #share the same grid and a positional comparison is exactly what's + #wanted here. + for var in base_commit_dataset.data_vars.keys(): + base_vals = base_commit_dataset[var].values + later_vals = later_commit_dataset[var].values + maxval = np.abs(later_vals - base_vals).max() if maxval !=0.: - print('Test failed for '+var+' max diff value = '+str(maxval.values)) + print('Test failed for '+var+' max diff value = '+str(maxval)) test_pass = False base_experiment_input_nml = f90nml.read(data_dir_dict[base_commit] +'/run0001/input.nml') diff --git a/src/atmos_param/cg_drag/cg_drag.f90 b/src/atmos_param/cg_drag/cg_drag.f90 index 41bee433c..879021181 100644 --- a/src/atmos_param/cg_drag/cg_drag.f90 +++ b/src/atmos_param/cg_drag/cg_drag.f90 @@ -1,3 +1,4 @@ +!Chaim cg_drag with Stepehens cg_drag edits module cg_drag_mod use fms_mod, only: fms_init, mpp_pe, mpp_root_pe, & @@ -29,11 +30,11 @@ module cg_drag_mod private !--------------------------------------------------------------------- -! cg_drag_mod computes the convective gravity wave forcing on -! the zonal flow. the parameterization is described in Alexander and -! Dunkerton [JAS, 15 December 1999]. +! cg_drag_mod computes the convective gravity wave forcing on +! the zonal flow. the parameterization is described in Alexander and +! Dunkerton [JAS, 15 December 1999]. !-------------------------------------------------------------------- - + !--------------------------------------------------------------------- !----------- ****** VERSION NUMBER ******* --------------------------- @@ -63,7 +64,7 @@ module cg_drag_mod !!$integer :: vers, old_time_step !wfc++ Addition for regular use - integer, allocatable, dimension(:,:) :: source_level + integer, allocatable, dimension(:,:) :: source_level, damp_level real, allocatable, dimension(:,:) :: source_amp real, allocatable, dimension(:,:,:) :: gwd_u, gwd_v @@ -79,68 +80,100 @@ module cg_drag_mod ! at 00Z and calculations are not done ! every time step -real :: source_level_pressure= 315.e+02 - ! highest model level with pressure +real :: source_level_pressure= 315.e+02 + ! highest model level with pressure ! greater than this value (or sigma ! greater than this value normalized ! by 1013.25 hPa) will be the gravity - ! wave source level at the equator + ! wave source level at the equator ! [ Pa ] -integer :: nk=1 ! number of wavelengths contained in +real :: damp_level_pressure=0.8e+02 + ! added by cig, feb 27, 2017. any waves reaching the top level will be deposited down to this level +integer :: nk=1 ! number of wavelengths contained in ! the gravity wave spectrum real :: cmax=99.6 ! maximum phase speed in gravity wave ! spectrum [ m/s ] -real :: dc=1.2 ! gravity wave spectral resolution +real :: dc=1.2 ! gravity wave spectral resolution ! [ m/s ] ! previous values: 0.6 -real :: Bt_0=.004 ! sum across the wave spectrum of +real :: Bt_0=.004 ! sum across the wave spectrum of ! the magnitude of total GW stress [Pa] - -real :: Bt_aug=.000 ! magnitude of momentum flux divided by density + +real :: Bt_aug=.000 ! magnitude of momentum flux divided by density real :: Bt_nh=.001 ! additional momentum stress for NH [Pa] real :: Bt_sh=-.001 ! additional momentum stress for SH [Pa] -real :: Bt_eq=.000 ! additional momentum stress at equator - CURRENTLY NOT USED! +! epg - 30.6.16 - I shifted these spectral parameters to the name list + +!--------------------------------------------------------------------- +!--------------------------------------------------------------------- +! wave spectrum parameters. +!--------------------------------------------------------------------- + +integer :: flag = 1 ! flag = 1 for peak flux at c = 0 + ! flag = 0 for peak flux at (c-u) = 0 +real :: Bw = 0.4 ! amplitude for the wide spectrum [ m^2/s^2 ] + ! ~ u'w' +real :: Bn = 0.0 ! amplitude for the narrow spectrum [ m^2/s^2 ] + ! ~ u'w'; previous values: 5.4 +real :: cw = 40.0 ! half-width for the wide c spectrum [ m/s ] + ! previous values: 50.0, 25.0 +real :: cwtropics = 40.0 ! half-width for the wide c spectrum [ m/s ] + ! previous values: 50.0, 25.0 +real :: cn = 2.0 ! half-width for the narrow c spectrum [ m/s ] + +real :: Bt_eq=.000 ! additional momentum stress at equator - CURRENTLY NOT USED! real :: Bt_eq_width=4.0 ! scaling for width of equtorial momentum flux (equator) CURRENTLY NOT USED! real :: phi0n = 30., phi0s = -30., dphin = 5., dphis = -5. -logical :: calculate_ked=.false. +!add by chaim jan 2017 +real :: weightminus2=0. + +real :: weightminus1=0. + +real :: weighttop=1. + +real :: kelvin_kludge=1. + +logical :: calculate_ked=.false. ! calculate ked diagnostic ? integer :: num_diag_pts_ij=0 ! number of diagnostic columns specif- ! ied by global (i,j) coordinates -integer :: num_diag_pts_latlon=0 +integer :: num_diag_pts_latlon=0 ! number of diagnostic columns ! specified by lat-lon coordinates integer, parameter :: MAX_PTS= 20 ! maximum number of diagnostic columns -integer, dimension(MAX_PTS) :: i_coords_gl=-100 - ! global i coordinates for ij - ! diagnostic columns -integer, dimension(MAX_PTS) :: j_coords_gl=-100 - ! global j coordinates for ij - ! diagnostic columns -real, dimension(MAX_PTS) :: lat_coords_gl=-999. - ! latitudes for latlon diagnostic +integer, dimension(MAX_PTS) :: i_coords_gl=-100 + ! global i coordinates for ij + ! diagnostic columns +integer, dimension(MAX_PTS) :: j_coords_gl=-100 + ! global j coordinates for ij + ! diagnostic columns +real, dimension(MAX_PTS) :: lat_coords_gl=-999. + ! latitudes for latlon diagnostic ! columns [degrees, -90. -> 90. ] -real, dimension(MAX_PTS) :: lon_coords_gl=-999. - ! longitudes for latlon diagnostic +real, dimension(MAX_PTS) :: lon_coords_gl=-999. + ! longitudes for latlon diagnostic ! columns [ degrees, 0. -> 360. ] namelist / cg_drag_nml / & cg_drag_freq, cg_drag_offset, & - source_level_pressure, & + source_level_pressure, damp_level_pressure, & nk, cmax, dc, Bt_0, Bt_aug, & Bt_sh, Bt_nh, Bt_eq, Bt_eq_width, & calculate_ked, & num_diag_pts_ij, num_diag_pts_latlon, & i_coords_gl, j_coords_gl, & lat_coords_gl, lon_coords_gl, & - phi0n,phi0s,dphin,dphis + phi0n,phi0s,dphin,dphis, Bw, Bn, cw, cwtropics, cn, flag, & + weightminus2, weightminus1, weighttop,kelvin_kludge + !-------------------------------------------------------------------- !-------- public data ----- @@ -154,16 +187,16 @@ module cg_drag_mod !-------------------------------------------------------------------- !mj remove restart stuff integer, dimension(3) :: restart_versions = (/ 1, 2, 3 /) ! v1 : -! v2 : +! v2 : ! v3 : Now use NetCDF for restart file. ! !-------------------------------------------------------------------- ! these arrays must be preserved across timesteps in case the ! parameterization is not called every timestep: ! -! gwd time tendency for u eqn due to gravity wave forcing +! gwd time tendency for u eqn due to gravity wave forcing ! [ m/s^2 ] -! ked effective eddy diffusion coefficient resulting from +! ked effective eddy diffusion coefficient resulting from ! gravity wave forcing [ m^2/s ] ! !-------------------------------------------------------------------- @@ -180,28 +213,18 @@ module cg_drag_mod !------------------------------------------------------------------- real, dimension(:), allocatable :: c0, kwv, k2 - -!--------------------------------------------------------------------- -! wave spectrum parameters. -!--------------------------------------------------------------------- integer :: nc ! number of wave speeds in spectrum ! (symmetric around c = 0) -integer :: flag = 1 ! flag = 1 for peak flux at c = 0 - ! flag = 0 for peak flux at (c-u) = 0 -real :: Bw = 0.4 ! amplitude for the wide spectrum [ m^2/s^2 ] - ! ~ u'w' -real :: Bn = 0.0 ! amplitude for the narrow spectrum [ m^2/s^2 ] - ! ~ u'w'; previous values: 5.4 -real :: cw = 40.0 ! half-width for the wide c spectrum [ m/s ] - ! previous values: 50.0, 25.0 -real :: cn = 2.0 ! half-width for the narrow c spectrum [ m/s ] -integer :: klevel_of_source +integer :: klevel_of_source, klevel_of_damp ! k index of the gravity wave source level at ! the equator in a standard atmosphere + ! also k index of level up to where mesosphere drag is dumped (cig, feb 27 2017) + + !--------------------------------------------------------------------- ! variables which control module calculations: -! +! ! cgdrag_alarm time remaining until next cg_drag calculation [ s ] ! !--------------------------------------------------------------------- @@ -211,29 +234,30 @@ module cg_drag_mod ! variables used with column diagnostics: ! ! diag_units output unit numbers -! num_diag_pts number of columns where diagnostics are desired +! num_diag_pts number of columns where diagnostics are desired ! column_diagnostics_desired ! column diagnostics are desired ? -! do_column_diagnostics -! a diagnostic column is in this jrow ? +! do_column_diagnostics +! a diagnostic column is in this jrow ? ! diag_lon longitude of diagnostic columns [ degrees ] ! diag_lat latiude of diagnostic columns [ degrees ] ! diag_i processor-based i index of diagnostic columns ! diag_j processor-based j index of diagnostic columns ! !-------------------------------------------------------------------- -integer :: num_diag_pts = 0 +integer :: num_diag_pts = 0 logical :: column_diagnostics_desired=.false. -integer, dimension(:), allocatable :: diag_units +integer, dimension(:), allocatable :: diag_units logical, dimension(:), allocatable :: do_column_diagnostics real, dimension(:), allocatable :: diag_lon, diag_lat -integer, dimension(:), allocatable :: diag_j, diag_i +integer, dimension(:), allocatable :: diag_j, diag_i !--------------------------------------------------------------------- ! variables for netcdf diagnostic fields. !--------------------------------------------------------------------- integer :: id_kedx_cgwd, id_kedy_cgwd, id_bf_cgwd, & - id_gwfx_cgwd, id_gwfy_cgwd +! id_gwfx_cgwd, id_gwfy_cgwd + id_gwfx_cgwd, id_gwfy_cgwd, id_source_level real :: missing_value = -999. character(len=7) :: mod_name = 'cg_drag' @@ -256,8 +280,8 @@ module cg_drag_mod !#################################################################### -subroutine cg_drag_init (lonb, latb, pref, Time, axes) - +!subroutine cg_drag_init (lonb, latb, pref, Time, axes) RC 5/12/19 +subroutine cg_drag_init (is, ie, js, je, lonb, latb, pref, Time, axes) !------------------------------------------------------------------- ! cg_drag_init is the constructor for cg_drag_mod. @@ -268,6 +292,7 @@ subroutine cg_drag_init (lonb, latb, pref, Time, axes) real, dimension(:), intent(in) :: lonb, latb, pref integer, dimension(4), intent(in) :: axes type(time_type), intent(in) :: Time +integer, intent(in) :: is, ie, js, je !RC 5/12/19 !------------------------------------------------------------------- !------------------------------------------------------------------- @@ -284,7 +309,7 @@ subroutine cg_drag_init (lonb, latb, pref, Time, axes) !------------------------------------------------------------------ !------------------------------------------------------------------- -! local variables: +! local variables: integer :: unit, ierr, io, logunit integer :: n, i, j, k @@ -297,11 +322,11 @@ subroutine cg_drag_init (lonb, latb, pref, Time, axes) ! real :: lat(size(lonb,1) - 1, size(latb,2) - 1) real :: lat(size(lonb) - 1, size(latb) - 1) !------------------------------------------------------------------- -! local variables: -! -! unit unit number for nml file -! ierr error return flag -! io error return code +! local variables: +! +! unit unit number for nml file +! ierr error return flag +! io error return code ! n loop index ! k loop index ! idf number of i points on this processor @@ -323,7 +348,7 @@ subroutine cg_drag_init (lonb, latb, pref, Time, axes) call diag_manager_init call constants_init #ifdef COL_DIAG - call column_diagnostics_init +! call column_diagnostics_init #endif SKIP !--------------------------------------------------------------------- ! read namelist. @@ -345,7 +370,7 @@ subroutine cg_drag_init (lonb, latb, pref, Time, axes) if (mpp_pe() == mpp_root_pe()) write (logunit, nml=cg_drag_nml) !------------------------------------------------------------------- -! define the grid dimensions. idf and jdf are the (i,j) dimensions of +! define the grid dimensions. idf and jdf are the (i,j) dimensions of ! domain on this processor, kmax is the number of model layers. !------------------------------------------------------------------- kmax = size(pref(:)) - 1 @@ -356,6 +381,7 @@ subroutine cg_drag_init (lonb, latb, pref, Time, axes) idf = size(lonb) - 1 allocate( source_level(idf,jdf) ) + allocate( damp_level(idf,jdf) ) allocate( source_amp(idf,jdf) ) ! allocate( lat(idf,jdf) ) @@ -365,12 +391,16 @@ subroutine cg_drag_init (lonb, latb, pref, Time, axes) ! ied as the source location via namelist input. !-------------------------------------------------------------------- do k=1,kmax + if (pref(k) < damp_level_pressure) then + klevel_of_damp = k + endif if (pref(k) > source_level_pressure) then klevel_of_source = k exit endif end do + do j=1,jdf !mj change of dimensions ! lat(:,j)= 0.5*( latb(:,j+1)+latb(:,j) ) @@ -378,12 +408,34 @@ subroutine cg_drag_init (lonb, latb, pref, Time, axes) lat(i,j)= 0.5*( latb(j+1)+latb(j) ) source_level(i,j) = (kmax + 1) - ((kmax + 1 - & klevel_of_source)*cos(lat(i,j)) + 0.5) - source_amp(i,j) = Bt_0 + & - Bt_nh*0.5*(1.+tanh((lat(i,j)/pif-phi0n)/dphin)) + & - Bt_sh*0.5*(1.+tanh((lat(i,j)/pif-phi0s)/dphis)) + + damp_level(i,j) = klevel_of_damp !cig + +!code added by ipw - nov 23, 2016 + if (lat(i,j) > phi0n) then + source_amp(i,j) = Bt_0 + Bt_nh*0.5*(1.+tanh((lat(i,j)/pif-phi0n)/dphin))+ & + Bt_sh*0.5*(1.+tanh((lat(i,j)/pif-phi0s)/dphis)); + elseif (lat(i,j) < phi0s) then + source_amp(i,j) = Bt_0 + Bt_nh*0.5*(1.+tanh((lat(i,j)/pif-phi0n)/dphin))+ & + Bt_sh*0.5*(1.+tanh((lat(i,j)/pif-phi0s)/dphis)); + elseif ((lat(i,j) <= dphin) .and. (lat(i,j) >= dphis)) then + source_amp(i,j) = Bt_eq + elseif ((lat(i,j) <= phi0n) .and. (lat(i,j) > dphin)) then + source_amp(i,j) = Bt_0 + (Bt_eq-Bt_0)/(phi0n-dphin)*(phi0n-lat(i,j)) + elseif ((lat(i,j) < dphis) .and. (lat(i,j) >= phi0s)) then + source_amp(i,j) = Bt_0 + (Bt_eq-Bt_0)/(phi0s-dphis)*(phi0s-lat(i,j)) + endif + +! source_amp(i,j) = Bt_0 + & +! Bt_nh*0.5*(1.+tanh((lat(i,j)/pif-phi0n)/dphin)) + & +! Bt_sh*0.5*(1.+tanh((lat(i,j)/pif-phi0s)/dphis)) end do end do source_level = MIN (source_level, kmax-1) + damp_level = MIN (damp_level, kmax) + +!cig: make sure everyhing is ok +! write (*,*) "damp",pref(klevel_of_damp), ' ', klevel_of_damp, ' ', damp_level_pressure, ' ', damp_level(2,2), ' ', damp_level(12,2) ! deallocate( lat ) @@ -398,7 +450,7 @@ subroutine cg_drag_init (lonb, latb, pref, Time, axes) !--------------------------------------------------------------------- ! if column diagnostics are desired, check that array dimensions are -! sufficiently large for the number of requests. +! sufficiently large for the number of requests. !--------------------------------------------------------------------- #ifdef COL_DIAG if (column_diagnostics_desired) then @@ -409,7 +461,7 @@ subroutine cg_drag_init (lonb, latb, pref, Time, axes) endif !--------------------------------------------------------------------- -! allocate arrays needed for column diagnostics. +! allocate arrays needed for column diagnostics. !--------------------------------------------------------------------- allocate (do_column_diagnostics (jdf) ) allocate (diag_units (num_diag_pts) ) @@ -419,15 +471,15 @@ subroutine cg_drag_init (lonb, latb, pref, Time, axes) allocate (diag_j (num_diag_pts) ) !--------------------------------------------------------------------- -! call initialize_diagnostic_columns to determine the locations -! (i, j, lat and lon) of any diagnostic columns in this processsor's +! call initialize_diagnostic_columns to determine the locations +! (i, j, lat and lon) of any diagnostic columns in this processsor's ! space and to open output files for the diagnostics. !--------------------------------------------------------------------- call initialize_diagnostic_columns & (mod_name, num_diag_pts_latlon, num_diag_pts_ij, & i_coords_gl, j_coords_gl, lat_coords_gl, & !mj dimensions are different -! lon_coords_gl, lonb(:,1), latb(1,:), do_column_diagnostics, & + lon_coords_gl, lonb(:,1), latb(1,:), do_column_diagnostics, & lon_coords_gl, lonb, latb, do_column_diagnostics, & diag_lon, diag_lat, diag_i, diag_j, diag_units) endif @@ -443,11 +495,11 @@ subroutine cg_drag_init (lonb, latb, pref, Time, axes) do n=1,nc c0(n) = (n-1)*dc - cmax end do - + !-------------------------------------------------------------------- -! define the wavenumber kwv and its square k2 for the gravity waves -! contained in the spectrum. currently nk = 1, which means that the -! wavelength of all gravity waves considered is 300 km. +! define the wavenumber kwv and its square k2 for the gravity waves +! contained in the spectrum. currently nk = 1, which means that the +! wavelength of all gravity waves considered is 300 km. !-------------------------------------------------------------------- allocate ( kwv(nk) ) allocate ( k2 (nk) ) @@ -479,13 +531,20 @@ subroutine cg_drag_init (lonb, latb, pref, Time, axes) register_diag_field (mod_name, 'kedy_cgwd', axes(1:3), Time, & 'effective eddy viscosity from cg_drag', 'm^2/s', & missing_value=missing_value) + id_source_level = & + register_diag_field (mod_name, 'source_level', axes(1:2), Time, & + 'effective eddy viscosity from cg_drag', 'm^2/s', & + missing_value=missing_value) !-------------------------------------------------------------------- -! allocate and define module variables to hold values across +! allocate and define module variables to hold values across ! timesteps, in the event that cg_drag is not called on every step. !-------------------------------------------------------------------- - allocate ( gwd_u(idf,jdf,kmax) ) - allocate ( gwd_v(idf,jdf,kmax) ) + +! allocate ( gwd_u(idf,jdf,kmax) ) RC 06/12/19 +! allocate ( gwd_v(idf,jdf,kmax) ) RC 06/12/19 + allocate ( gwd_u(is:ie,js:je,kmax) ) + allocate ( gwd_v(is:ie,js:je,kmax) ) !-------------------------------------------------------------------- ! if present, read the restart data file. @@ -516,12 +575,12 @@ subroutine cg_drag_init (lonb, latb, pref, Time, axes) gwd_v(:,:,:) = 0.0 if (cg_drag_offset > 0) then cgdrag_alarm = cg_drag_offset - else + else cgdrag_alarm = cg_drag_freq endif !!$ endif !!$ vers = restart_versions(size(restart_versions(:))) -!!$ old_time_step = cgdrag_alarm +!!$ old_time_step = cgdrag_alarm !--------------------------------------------------------------------- ! mark the module as initialized. !--------------------------------------------------------------------- @@ -535,7 +594,7 @@ end subroutine cg_drag_init !#################################################################### - + subroutine cg_drag_time_vary (delt) real , intent(in) :: delt @@ -546,17 +605,17 @@ subroutine cg_drag_time_vary (delt) cgdrag_alarm = cgdrag_alarm - delt !--------------------------------------------------------------------- - + end subroutine cg_drag_time_vary !#################################################################### - + subroutine cg_drag_endts - + !-------------------------------------------------------------------- -! if this was a calculation step, reset cgdrag_alarm to indicate -! the time remaining before the next calculation of gravity wave +! if this was a calculation step, reset cgdrag_alarm to indicate +! the time remaining before the next calculation of gravity wave ! forcing. !-------------------------------------------------------------------- if (cgdrag_alarm <= 0 ) then @@ -570,9 +629,9 @@ end subroutine cg_drag_endts subroutine cg_drag_calc (is, js, lat, pfull, zfull, temp, uuu, vvv, & Time, delt, gwfcng_x, gwfcng_y) -!-------------------------------------------------------------------- +!-------------------------------------------------------------------- ! cg_drag_calc defines the arrays needed to calculate the convective -! gravity wave forcing, calls gwfc to calculate the forcing, returns +! gravity wave forcing, calls gwfc to calculate the forcing, returns ! the desired output fields, and saves the values for later retrieval ! if they are not calculated on every timestep. ! @@ -589,7 +648,7 @@ subroutine cg_drag_calc (is, js, lat, pfull, zfull, temp, uuu, vvv, & !------------------------------------------------------------------- ! intent(in) variables: ! -! is,js starting subdomain i,j indices of data in +! is,js starting subdomain i,j indices of data in ! the physics_window being integrated ! lat array of model latitudes at cell boundaries [radians] ! pfull pressure at model full levels [ Pa ] @@ -632,24 +691,25 @@ subroutine cg_drag_calc (is, js, lat, pfull, zfull, temp, uuu, vvv, & ! local variables: ! ! dtdz temperature lapse rate [ deg K/m ] -! ked_gwfc effective diffusion coefficient from cg_drag_mod +! ked_gwfc effective diffusion coefficient from cg_drag_mod ! [ m^2/s ] ! zzchm heights at model levels [ m ] ! zu zonal velocity [ m/s ] ! zden atmospheric density [ kg/m^3 ] ! zbf buoyancy frequency [ /s ] -! gwd_xtnd zonal wind tendency resulting from cg_drag_mod +! gwd_xtnd zonal wind tendency resulting from cg_drag_mod ! [ m/s^2 ] -! ked_xtnd effective diffusion coefficient from cg_drag_mod +! ked_xtnd effective diffusion coefficient from cg_drag_mod ! [ m^2/s ] ! source_level k index of gravity wave source level ((i,j) array) +! damp_level k index of gravity wave mesospheric dumping level ((i,j) array) ! iz0 k index of gravity wave source level in a column ! used return code for netcdf diagnostics -! bflim minimum allowable value of squared buoyancy +! bflim minimum allowable value of squared buoyancy ! frequency [ /s^2 ] -! ie, je ending subdomain indices of data in the current +! ie, je ending subdomain indices of data in the current ! physics window being integrated -! imax, jmax, kmax +! imax, jmax, kmax ! physics window dimensions ! i, j, k, nn do loop indices ! @@ -665,24 +725,25 @@ subroutine cg_drag_calc (is, js, lat, pfull, zfull, temp, uuu, vvv, & je = js + jmax - 1 !--------------------------------------------------------------------- -! if the convective gravity wave forcing should be calculated on +! if the convective gravity wave forcing should be calculated on ! this timestep (i.e., the alarm has gone off), proceed with the ! calculation. !--------------------------------------------------------------------- if (cgdrag_alarm <= 0) then - + write(6,*) 'GWD Calculated', cgdrag_alarm !----------------------------------------------------------------------- -! calculate temperature lapse rate. do one-sided differences over -! delta z at upper boundary and centered differences over 2 delta z +! calculate temperature lapse rate. do one-sided differences over +! delta z at upper boundary and centered differences over 2 delta z ! in the interior. dtdz is not needed at the lower boundary, since ! the source level is constrained to be above level kmax. !---------------------------------------------------------------------- do j=1,jmax do i=1,imax ! The following index-offsets are needed in case a physics_window is being used. - iz0 = source_level(i +is-1,j+js-1) + !iz0 = source_level(i +is-1,j+js-1) + iz0 = source_level(i,j) dtdz(i,j,1) = (temp (i,j,1) - temp (i,j,2))/ & (zfull(i,j,1) - zfull(i,j,2)) do k=2,iz0 @@ -698,14 +759,14 @@ subroutine cg_drag_calc (is, js, lat, pfull, zfull, temp, uuu, vvv, & end do !---------------------------------------------------------------------- -! calculate buoyancy frequency. restrict the squared buoyancy +! calculate buoyancy frequency. restrict the squared buoyancy ! frequency to be no smaller than bflim. !---------------------------------------------------------------------- - do k=1,iz0 + do k=1,iz0 zbf(i,j,k) = (GRAV/temp(i,j,k))*(dtdz(i,j,k) + GRAV/CP_AIR) if (zbf(i,j,k) < bflim) then zbf(i,j,k) = sqrt(bflim) - else + else zbf(i,j,k) = sqrt(zbf(i,j,k)) endif end do @@ -725,7 +786,7 @@ subroutine cg_drag_calc (is, js, lat, pfull, zfull, temp, uuu, vvv, & do k=1,iz0+1 zzchm(i,j,k) = zfull(i,j,k) end do - do k=1,iz0 + do k=1,iz0 zu(i,j,k) = uuu(i,j,k) zv(i,j,k) = vvv(i,j,k) end do @@ -734,8 +795,8 @@ subroutine cg_drag_calc (is, js, lat, pfull, zfull, temp, uuu, vvv, & ! add an extra level above model top so that the gravity wave forcing ! occurring between the topmost model level and the upper boundary ! may be calculated. define variable values at the new top level as -! follows: z - use delta z of layer just below; u - extend vertical -! gradient occurring just below; density - geometric mean; buoyancy +! follows: z - use delta z of layer just below; u - extend vertical +! gradient occurring just below; density - geometric mean; buoyancy ! frequency - constant across model top. !---------------------------------------------------------------------- zzchm(i,j,0) = zzchm(i,j,1) + zzchm(i,j,1) - zzchm(i,j,2) @@ -745,25 +806,25 @@ subroutine cg_drag_calc (is, js, lat, pfull, zfull, temp, uuu, vvv, & zbf(i,j,0) = zbf(i,j,1) end do end do - + !--------------------------------------------------------------------- ! pass the vertically-extended input arrays to gwfc. gwfc will cal- -! culate the gravity-wave forcing and, if desired, an effective eddy +! culate the gravity-wave forcing and, if desired, an effective eddy ! diffusion coefficient at each level above the source level. output ! is returned in the vertically-extended arrays gwfcng and ked_gwfc. -! upon return move the output fields into model-sized arrays. +! upon return move the output fields into model-sized arrays. !--------------------------------------------------------------------- - call gwfc (is, ie, js, je, source_level, source_amp, & + call gwfc (is, ie, js, je, damp_level, source_level, source_amp, lat, & zden, zu, zbf,zzchm, gwd_xtnd, ked_xtnd) gwfcng_x (:,:,1:kmax) = gwd_xtnd(:,:,1:kmax ) ked_gwfc_x(:,:,1:kmax) = ked_xtnd(:,:,1:kmax ) - - call gwfc (is, ie, js, je, source_level, source_amp, & + + call gwfc (is, ie, js, je, damp_level, source_level, source_amp, lat, & zden, zv, zbf,zzchm, gwd_ytnd, ked_ytnd) gwfcng_y (:,:,1:kmax) = gwd_ytnd(:,:,1:kmax ) ked_gwfc_y(:,:,1:kmax) = ked_ytnd(:,:,1:kmax ) - + !-------------------------------------------------------------------- ! store the gravity wave forcing into a processor-global array. !------------------------------------------------------------------- @@ -771,91 +832,89 @@ subroutine cg_drag_calc (is, js, lat, pfull, zfull, temp, uuu, vvv, & gwd_v(is:ie,js:je,:) = gwfcng_y(:,:,:) -#ifdef COL_DIAG -!-------------------------------------------------------------------- -! if column diagnostics are desired, determine if any columns are on -! this processor. if so, call column_diagnostics_header to write -! out location and timestamp information. then output desired -! quantities to the diag_unit file. -!--------------------------------------------------------------------- - if (column_diagnostics_desired) then - do j=1,jmax - if (do_column_diagnostics(j+js-1)) then - do nn=1,num_diag_pts - if (js + j - 1 == diag_j(nn)) then - call column_diagnostics_header & - (mod_name, diag_units(nn), Time, nn, diag_lon, & - diag_lat, diag_i, diag_j) - iz0 = source_level (diag_i(nn), j) - write (diag_units(nn),'(a, i5)') & - ' source_level =', iz0 - write (diag_units(nn),'(a)') & - ' k u z density& - & bf gwforcing' - do k=0,iz0 - write (diag_units(nn), '(i5, 2x, 5e12.5)') & - k, & - zu (diag_i(nn),j,k), & - zzchm (diag_i(nn),j,k), & - zden (diag_i(nn),j,k), & - zbf (diag_i(nn),j,k), & - gwd_xtnd (diag_i(nn),j,k) - end do - write (diag_units(nn), '(i5, 14x, 2e12.5)') & - iz0+1, & - zzchm (diag_i(nn),j,iz0+1), & - zden (diag_i(nn),j,iz0+1) - endif - end do ! (nn loop) - endif ! (do_column_diagnostics) - end do ! (j loop) - endif ! (column_diagnostics_desired) -#endif +! #ifdef COL_DIAG +! !-------------------------------------------------------------------- +! ! if column diagnostics are desired, determine if any columns are on +! ! this processor. if so, call column_diagnostics_header to write +! ! out location and timestamp information. then output desired +! ! quantities to the diag_unit file. +! !--------------------------------------------------------------------- +! if (column_diagnostics_desired) then +! do j=1,jmax +! if (do_column_diagnostics(j+js-1)) then +! do nn=1,num_diag_pts +! if (js + j - 1 == diag_j(nn)) then +! call column_diagnostics_header & +! (mod_name, diag_units(nn), Time, nn, diag_lon, & +! diag_lat, diag_i, diag_j) +! iz0 = source_level (diag_i(nn), j) +! write (diag_units(nn),'(a, i5)') & +! ' source_level =', iz0 +! write (diag_units(nn),'(a)') & +! ' k u z density& +! & bf gwforcing' +! do k=0,iz0 +! write (diag_units(nn), '(i5, 2x, 5e12.5)') & +! k, & +! zu (diag_i(nn),j,k), & +! zzchm (diag_i(nn),j,k), & +! zden (diag_i(nn),j,k), & +! zbf (diag_i(nn),j,k), & +! gwd_xtnd (diag_i(nn),j,k) +! end do +! write (diag_units(nn), '(i5, 14x, 2e12.5)') & +! iz0+1, & +! zzchm (diag_i(nn),j,iz0+1), & +! zden (diag_i(nn),j,iz0+1) +! endif +! end do ! (nn loop) +! endif ! (do_column_diagnostics) +! end do ! (j loop) +! endif ! (column_diagnostics_desired) +! #endif !-------------------------------------------------------------------- -! if activated, store the effective eddy diffusivity into a -! processor-global array, and if desired as a netcdf diagnostic, +! if activated, store the effective eddy diffusivity into a +! processor-global array, and if desired as a netcdf diagnostic, ! send the data to diag_manager_mod. !------------------------------------------------------------------- if (id_kedx_cgwd > 0) then -! used = send_data (id_kedx_cgwd, ked_gwfc_x, Time, is, js, 1) - used = send_data (id_kedx_cgwd, ked_gwfc_x, Time) !st 2013 FMS seems to have paralelisation issues when called with ...Time, is, js, 1). + used = send_data (id_kedx_cgwd, ked_gwfc_x, Time)!, is, js, 1) endif if (id_kedy_cgwd > 0) then -! used = send_data (id_kedy_cgwd, ked_gwfc_y, Time, is, js, 1) - used = send_data (id_kedy_cgwd, ked_gwfc_y, Time) + used = send_data (id_kedy_cgwd, ked_gwfc_y, Time)!, is, js, 1) endif - + if (id_source_level > 0) then + used = send_data (id_source_level, real(source_level), Time) + endif !-------------------------------------------------------------------- ! save any other netcdf file diagnostics that are desired. !-------------------------------------------------------------------- if (id_bf_cgwd > 0) then -! used = send_data (id_bf_cgwd, zbf(:,:,1:), Time, is, js ) - used = send_data (id_bf_cgwd, zbf(:,:,1:), Time) + used = send_data (id_bf_cgwd, zbf(:,:,1:), Time)!, is, js ) endif if (id_gwfx_cgwd > 0) then -! used = send_data (id_gwfx_cgwd, gwfcng_x, Time, is, js, 1) - used = send_data (id_gwfx_cgwd, gwfcng_x, Time) + used = send_data (id_gwfx_cgwd, gwfcng_x, Time)!, is, js, 1) endif if (id_gwfy_cgwd > 0) then -! used = send_data (id_gwfy_cgwd, gwfcng_y, Time, is, js, 1) - used = send_data (id_gwfy_cgwd, gwfcng_y, Time) + used = send_data (id_gwfy_cgwd, gwfcng_y, Time)!, is, js, 1) endif !-------------------------------------------------------------------- -! if this is not a timestep on which gravity wave forcing is to be +! if this is not a timestep on which gravity wave forcing is to be ! calculated, retrieve the values calculated previously from storage ! and return to the calling subroutine. !-------------------------------------------------------------------- else ! (cgdrag_alarm <= 0) + write(6,*) 'GWD not recalculated',cgdrag_alarm gwfcng_x(:,:,:) = gwd_u(is:ie,js:je,:) gwfcng_y(:,:,:) = gwd_v(is:ie,js:je,:) endif ! (cgdrag_alarm <= 0) @@ -865,7 +924,7 @@ subroutine cg_drag_calc (is, js, lat, pfull, zfull, temp, uuu, vvv, & ! will recalculate the NOGWD tendencies call cg_drag_endts call cg_drag_time_vary(delt) - + end subroutine cg_drag_calc @@ -927,8 +986,8 @@ end subroutine cg_drag_end !!$ unit = open_restart_file ('RESTART/cg_drag.res', 'write') !!$ !!$!------------------------------------------------------------------- -!!$! the root pe writes out the restart version, the time remaining -!!$! before the next call to cg_drag_mod and the current cg_drag +!!$! the root pe writes out the restart version, the time remaining +!!$! before the next call to cg_drag_mod and the current cg_drag !!$! timestep. !!$!------------------------------------------------------------------- !!$ if (mpp_pe() == mpp_root_pe() ) then @@ -937,7 +996,7 @@ end subroutine cg_drag_end !!$ endif !!$ !!$!------------------------------------------------------------------- -!!$! each processor writes out its gravity wave forcing tendency +!!$! each processor writes out its gravity wave forcing tendency !!$! on the zonal flow. !!$!------------------------------------------------------------------- !!$ call write_data (unit, gwd_u) @@ -968,11 +1027,11 @@ end subroutine cg_drag_end !!$ real :: secs_per_day = SECONDS_PER_DAY !!$ !!$!------------------------------------------------------------------- -!!$! local variables: -!!$! -!!$! unit unit number for nml file -!!$! chvers character representation of restart version -!!$! vers restart version +!!$! local variables: +!!$! +!!$! unit unit number for nml file +!!$! chvers character representation of restart version +!!$! vers restart version !!$! dummy array to hold restart version 1 control variables !!$! old_time_step cg_drag timestep used in previous model run [ s ] !!$! secs_per_day seconds in a day [ s ] @@ -981,7 +1040,7 @@ end subroutine cg_drag_end !!$ !!$ !!$!-------------------------------------------------------------------- -!!$! open file to read restart data. +!!$! open file to read restart data. !!$!--------------------------------------------------------------------- !!$ unit = open_restart_file ('INPUT/cg_drag.res','read') !!$ @@ -997,36 +1056,36 @@ end subroutine cg_drag_end !!$ endif !!$ !!$!-------------------------------------------------------------------- -!!$! read control information from restart file. +!!$! read control information from restart file. !!$!-------------------------------------------------------------------- !!$ if (vers == 1) then !!$ !!$!-------------------------------------------------------------------- !!$! if reading restart version 1, use the contents of array dummy to -!!$! define the cg_drag timestep that was used in the run which wrote -!!$! the restart. define the time remaining before the next cg_drag +!!$! define the cg_drag timestep that was used in the run which wrote +!!$! the restart. define the time remaining before the next cg_drag !!$! calculation to either be the previous timestep or the current !!$! offset, if that is specified. this assumes that the restart was !!$! written at 00Z. !!$!-------------------------------------------------------------------- -!!$ read (unit) dummy +!!$ read (unit) dummy !!$ old_time_step = secs_per_day*dummy(4) + dummy(3) !!$ if (cg_drag_offset == 0) then !!$ cgdrag_alarm = old_time_step !!$ else -!!$ cgdrag_alarm = cg_drag_offset +!!$ cgdrag_alarm = cg_drag_offset !!$ endif -!!$ else +!!$ else !!$ !!$!-------------------------------------------------------------------- -!!$! for restart version 2, read the time remaining until the next +!!$! for restart version 2, read the time remaining until the next !!$! cg_drag calculation, and the previously used timestep. !!$!--------------------------------------------------------------------- !!$ read (unit) cgdrag_alarm, old_time_step !!$ endif !!$ !!$!------------------------------------------------------------------- -!!$! read restart data (gravity wave forcing tendency terms) and close +!!$! read restart data (gravity wave forcing tendency terms) and close !!$! unit. !!$!------------------------------------------------------------------- !!$ call read_data (unit, gwd_u) @@ -1034,8 +1093,8 @@ end subroutine cg_drag_end !!$ call close_file (unit) !!$ !!$!-------------------------------------------------------------------- -!!$! if current cg_drag calling frequency differs from that previously -!!$! used, adjust the time remaining before the next calculation. +!!$! if current cg_drag calling frequency differs from that previously +!!$! used, adjust the time remaining before the next calculation. !!$!-------------------------------------------------------------------- !!$ if (cg_drag_freq /= old_time_step) then !!$ cgdrag_alarm = cgdrag_alarm - old_time_step + cg_drag_freq @@ -1048,7 +1107,7 @@ end subroutine cg_drag_end !!$ !!$!-------------------------------------------------------------------- !!$! if cg_drag_offset is specified and is smaller than the time remain- -!!$! ing until the next calculation, modify the time remaining to be +!!$! ing until the next calculation, modify the time remaining to be !!$! that offset time. the assumption is made that the restart was !!$! written at 00Z. !!$!-------------------------------------------------------------------- @@ -1066,8 +1125,8 @@ end subroutine cg_drag_end !!$ !!$subroutine read_nc_restart_file !!$!----------------------------------------------------------------------- -!!$! subroutine read_restart_nc reads a netcdf restart file to obtain -!!$! the variables needed upon experiment restart. +!!$! subroutine read_restart_nc reads a netcdf restart file to obtain +!!$! the variables needed upon experiment restart. !!$!----------------------------------------------------------------------- !!$ !!$!--------------------------------------------------------------------- @@ -1111,8 +1170,8 @@ end subroutine cg_drag_end !!$ vers = restart_versions(size(restart_versions(:))) !!$ !!$!-------------------------------------------------------------------- -!!$! if current cg_drag calling frequency differs from that previously -!!$! used, adjust the time remaining before the next calculation. +!!$! if current cg_drag calling frequency differs from that previously +!!$! used, adjust the time remaining before the next calculation. !!$!-------------------------------------------------------------------- !!$ if (cg_drag_freq /= old_time_step) then !!$ cgdrag_alarm = cgdrag_alarm - old_time_step + cg_drag_freq @@ -1126,7 +1185,7 @@ end subroutine cg_drag_end !!$ !!$!-------------------------------------------------------------------- !!$! if cg_drag_offset is specified and is smaller than the time remain- -!!$! ing until the next calculation, modify the time remaining to be +!!$! ing until the next calculation, modify the time remaining to be !!$! that offset time. the assumption is made that the restart was !!$! written at 00Z. !!$!-------------------------------------------------------------------- @@ -1144,10 +1203,10 @@ end subroutine cg_drag_end !!$subroutine cg_drag_register_restart !!$ !!$ character(len=64) :: fname = 'cg_drag.res.nc' ! name of restart file -!!$ character(len=64) :: fname2 +!!$ character(len=64) :: fname2 !!$ integer :: id_restart !!$ -!!$ call get_mosaic_tile_file(fname, fname2, .false. ) +!!$ call get_mosaic_tile_file(fname, fname2, .false. ) !!$ allocate(Cg_restart) !!$ if(trim(fname2) == trim(fname)) then !!$ Til_restart => Cg_restart @@ -1172,10 +1231,10 @@ end subroutine cg_drag_end !!$! !!$! !!$! write out restart file. -!!$! Arguments: -!!$! timestamp (optional, intent(in)) : A character string that represents the model time, +!!$! Arguments: +!!$! timestamp (optional, intent(in)) : A character string that represents the model time, !!$! used for writing restart. timestamp will append to -!!$! the any restart file name as a prefix. +!!$! the any restart file name as a prefix. !!$! !!$! !!$subroutine cg_drag_restart(timestamp) @@ -1190,21 +1249,21 @@ end subroutine cg_drag_end !#################################################################### -subroutine gwfc (is, ie, js, je, source_level, source_amp, rho, u, & +subroutine gwfc (is, ie, js, je, damp_level, source_level, source_amp, lat, rho, u, & bf, z, gwf, ked) !------------------------------------------------------------------- ! subroutine gwfc computes the gravity wave-driven-forcing on the -! zonal wind given vertical profiles of wind, density, and buoyancy -! frequency. -! Based on version implemented in SKYHI -- 27 Oct 1998 by M.J. +! zonal wind given vertical profiles of wind, density, and buoyancy +! frequency. +! Based on version implemented in SKYHI -- 27 Oct 1998 by M.J. ! Alexander and L. Bruhwiler. !------------------------------------------------------------------- !------------------------------------------------------------------- integer, intent(in) :: is, ie, js, je -integer, dimension(:,:), intent(in) :: source_level -real, dimension(:,:), intent(in) :: source_amp +integer, dimension(:,:), intent(in) :: source_level, damp_level +real, dimension(:,:), intent(in) :: source_amp, lat real, dimension(:,:,0:), intent(in) :: rho, u, bf, z real, dimension(:,:,0:), intent(out) :: gwf real, dimension(:,:,0:), intent(out) :: ked @@ -1216,9 +1275,10 @@ subroutine gwfc (is, ie, js, je, source_level, source_amp, rho, u, & ! in the physics_window being integrated ! source_level k index of model level serving as gravity wave ! source +! damp_level k index of the lowest model level at which all drag that reaches the model top is partially dumped ! source_amp amplitude of gravity wave source [Pa] -! -! rho atmospheric density [ kg/m^3 ] +! +! rho atmospheric density [ kg/m^3 ] ! u zonal wind component [ m/s ] ! bf buoyancy frequency [ /s ] ! z height of model levels [ m ] @@ -1229,7 +1289,7 @@ subroutine gwfc (is, ie, js, je, source_level, source_amp, rho, u, & ! ! intent(out), optional variables: ! -! ked eddy diffusion coefficient from gravity wave +! ked eddy diffusion coefficient from gravity wave ! forcing [ m^2/s ] ! !------------------------------------------------------------------ @@ -1244,28 +1304,28 @@ subroutine gwfc (is, ie, js, je, source_level, source_amp, rho, u, & real , dimension (nc) :: c0mu0, B0 real :: fm, fe, Hb, alp2, Foc, c, test, rbh,& eps, Bsum - integer :: iz0 + integer :: iz0, iztop integer :: i, j, k, ink, n - real :: ampl + real :: ampl, cwthis, Bnthis, flagthis !------------------------------------------------------------------ ! local variables: -! +! ! wv_frcng gravity wave forcing tendency [ m/s^2 ] ! diff_coeff eddy diffusion coefficient [ m2/s ] -! c0mu difference between phase speed of wave n and u +! c0mu difference between phase speed of wave n and u ! [ m/s ] ! dz delta z between model levels [ m ] -! fac factor used in determining if wave is breaking +! fac factor used in determining if wave is breaking ! [ s/m ] -! omc critical frequency that marks total internal +! omc critical frequency that marks total internal ! reflection [ /s ] -! msk indicator as to whether wave n is still propagating -! upwards (msk=1), or has been removed from the +! msk indicator as to whether wave n is still propagating +! upwards (msk=1), or has been removed from the ! spectrum because of breaking or reflection (msk=0) ! c0mu0 difference between phase speed of wave n and u at the ! source level [ m/s ] ! B0 wave momentum flux amplitude for wave n [ (m/s)^2 ] -! fm used to sum up momentum flux from all waves n +! fm used to sum up momentum flux from all waves n ! deposited at a level [ (m/s)^2 ] ! fe used to sum up contributions to diffusion coefficient ! from all waves n at a level [ (m/s)^3 ] @@ -1278,18 +1338,18 @@ subroutine gwfc (is, ie, js, je, source_level, source_amp, rho, u, & ! rbh atmospheric density at half-level (geometric mean) ! [ kg/m^3 ] ! eps intermittency factor -! Bsum total mag of gravity wave momentum flux at source +! Bsum total mag of gravity wave momentum flux at source ! level, divided by the density [ m^2/s^2 ] ! iz0 source level vertical index for the given column ! i,j,k spatial do loop indices -! ink wavenumber loop index -! n phase speed loop index -! ampl gravity wave stress [Pa] +! ink wavenumber loop index +! n phase speed loop index +! ampl gravity wave stress [Pa] ! !-------------------------------------------------------------------- !------------------------------------------------------------------- -! initialize the output arrays. these will hold values at each +! initialize the output arrays. these will hold values at each ! (i,j,k) point, summed over the wavelengths and phase speeds ! defining the gravity wave spectrum. !------------------------------------------------------------------- @@ -1297,19 +1357,39 @@ subroutine gwfc (is, ie, js, je, source_level, source_amp, rho, u, & ked = 0.0 do j=1,size(u,2) - do i=1,size(u,1) -! The following index-offsets are needed in case a physics_window is being used. - iz0 = source_level(i+is-1,j+js-1) - ampl= source_amp(i+is-1,j+js-1) + + do i=1,size(u,1) +!added by cig, january 2017, Changed indicies RC 2019 + !if ((lat(i+is-1,j+js-1) <= dphin) .and. (lat(i+is-1,j+js-1) >= dphis)) then + if ((lat(i,j) <= dphin) .and. (lat(i,j) >= dphis)) then + cwthis=cwtropics + Bnthis=0. + flag=0 + else + cwthis=cw + Bnthis=Bn + flagthis=flag + endif + +! RC Changed Indicies The following index-offsets are needed in case a physics_window is being used. +! iz0 = source_level(i+is-1,j+js-1) +! iztop = damp_level(i+is-1,j+js-1) +! ampl= source_amp(i+is-1,j+js-1) + +! The following index-offsets are needed in case a physics_window is being used. + iz0 = source_level(i,j) + ampl= source_amp(i,j) +! DO I ADD iztop=damp_level(i,j)? + iztop=damp_level(i,j) !-------------------------------------------------------------------- -! define wave momentum flux (B0) at source level for each phase -! speed n, and the sum over all phase speeds (Bsum), which is needed -! to calculate the intermittency. +! define wave momentum flux (B0) at source level for each phase +! speed n, and the sum over all phase speeds (Bsum), which is needed +! to calculate the intermittency. !------------------------------------------------------------------- Bsum = 0. do n=1,nc - c0mu0(n) = c0(n) - u(i,j,iz0) + c0mu0(n) = c0(n) - u(i,j,iz0) !--------------------------------------------------------------------- ! when the wave phase speed is same as wind speed, there is no @@ -1317,19 +1397,21 @@ subroutine gwfc (is, ie, js, je, source_level, source_amp, rho, u, & !--------------------------------------------------------------------- if (c0mu0(n) == 0.0) then B0(n) = 0.0 - else + else !--------------------------------------------------------------------- ! define wave momentum flux at source level for phase speed n. Add ! the contribution from this phase speed to the previous sum. !--------------------------------------------------------------------- - c = c0(n)*flag + c0mu0(n)*(1 - flag) + c = c0(n)*flagthis + c0mu0(n)*(1 - flagthis) if (c0mu0(n) < 0.0) then - B0(n) = -1.0*(Bw*exp(-alog(2.0)*(c/cw)**2) + & - Bn*exp(-alog(2.0)*(c/cn)**2)) - else - B0(n) = (Bw*exp(-alog(2.0)*(c/cw)**2) + & - Bn*exp(-alog(2.0)*(c/cn)**2)) + B0(n) = -1.0*(Bw*exp(-alog(2.0)*(c/cwthis)**2) + & + Bnthis*exp(-alog(2.0)*(c/cn)**2)) + B0(n) =B0(n)*kelvin_kludge + else + B0(n) = (Bw*exp(-alog(2.0)*(c/cwthis)**2) + & + Bnthis*exp(-alog(2.0)*(c/cn)**2)) + endif Bsum = Bsum + abs(B0(n)) endif @@ -1364,7 +1446,7 @@ subroutine gwfc (is, ie, js, je, source_level, source_amp, rho, u, & fac(k) = 0.5*(rho(i,j,k)/rho(i,j,iz0))*kwv(ink)/bf(i,j,k) end do - do k=0,iz0 + do k=0,iz0 dz(k) = z(i,j,k) - z(i,j,k+1) Hb = -(dz(k))/alog(rho(i,j,k)/rho(i,j,k+1)) alp2 = 0.25/(Hb*Hb) @@ -1373,14 +1455,14 @@ subroutine gwfc (is, ie, js, je, source_level, source_amp, rho, u, & end do !--------------------------------------------------------------------- -! initialize a flag which will indicate which waves are still +! initialize a flag which will indicate which waves are still ! propagating upwards. !--------------------------------------------------------------------- msk = 1 !---------------------------------------------------------------------- -! integrate upwards from the source level. define variables over -! which to sum the deposited flux and effective eddy diffusivity +! integrate upwards from the source level. define variables over +! which to sum the deposited flux and effective eddy diffusivity ! from all waves breaking at a given level. !---------------------------------------------------------------------- do k=iz0, 0, -1 @@ -1395,7 +1477,7 @@ subroutine gwfc (is, ie, js, je, source_level, source_amp, rho, u, & c0mu(k) = c0(n) - u(i,j,k) !---------------------------------------------------------------------- -! if phase speed matches the wind speed, remove c0(n) from the +! if phase speed matches the wind speed, remove c0(n) from the ! set of propagating waves. ! epg: This seems to be an unphysical decision, as the wave should ! break, having reached a critical level. But it's extremely @@ -1408,7 +1490,7 @@ subroutine gwfc (is, ie, js, je, source_level, source_amp, rho, u, & else !--------------------------------------------------------------------- -! define the criterion which determines if wave is reflected at this +! define the criterion which determines if wave is reflected at this ! level (test). !--------------------------------------------------------------------- test = abs(c0mu(k))*kwv(ink) - omc(k) @@ -1419,22 +1501,22 @@ subroutine gwfc (is, ie, js, je, source_level, source_amp, rho, u, & ! propagating set. !--------------------------------------------------------------------- msk(n) = 0 - else + else !--------------------------------------------------------------------- -! if wave is not reflected at this level, determine if it is -! breaking at this level (Foc >= 0), or if wave speed relative to -! windspeed has changed sign from its value at the source level +! if wave is not reflected at this level, determine if it is +! breaking at this level (Foc >= 0), or if wave speed relative to +! windspeed has changed sign from its value at the source level ! (c0mu0(n)*c0mu <= 0). if it is above the source level and is -! breaking, then add its momentum flux to the accumulated sum at -! this level, and increase the effective diffusivity accordingly. +! breaking, then add its momentum flux to the accumulated sum at +! this level, and increase the effective diffusivity accordingly. ! set flag to remove phase speed c0(n) from the set of active waves ! moving upwards to the next level. !--------------------------------------------------------------------- -! epg: if you are at the model top, deposit all momentum here, to +! epg: if you are at the model top, deposit all momentum here, to ! prevent waves from escaping the top of the model. See -! Shaw et al. 2010? for details on why this in important. +! Shaw et al. 2010? for details on why this in important. if ( k==0 ) then msk(n) = 0 if (k < iz0) then @@ -1442,7 +1524,7 @@ subroutine gwfc (is, ie, js, je, source_level, source_amp, rho, u, & fe = fe + c0mu(k)*B0(n) endif else - + Foc = B0(n)/(c0mu(k) )**3 - fac(k) if ((Foc >= 0.0) .or. & (c0mu0(n)*c0mu(k) <= 0.0)) then @@ -1459,51 +1541,77 @@ subroutine gwfc (is, ie, js, je, source_level, source_amp, rho, u, & end do ! phase speed loop !---------------------------------------------------------------------- -! compute the gravity wave momentum flux forcing and eddy +! compute the gravity wave momentum flux forcing and eddy ! diffusion coefficient obtained across the entire wave spectrum ! at this level. !---------------------------------------------------------------------- if ( k < iz0) then rbh = sqrt(rho(i,j,k)*rho(i,j,k+1)) wv_frcng(k) = ( rho(i,j,iz0)/rbh)*fm*eps/dz(k) - - !epg: enforce momentum conservation at model top; in this case, all the momentum + + !epg: enforce momentum conservation at model top; in this case, all the momentum ! deposited in the uppermost layer (which exist above the top model level, ! as explained in cg_drag_calc, must be added to the level below, which is ! the actual top level of the model. + !cig: place the extra momentum flux in the top 3 layers instead of all in the top layer if ( k==0 ) then - wv_frcng(k+1) = 0.5*wv_frcng(k+1) + wv_frcng(k) + wv_frcng(k+1) = 0.5*wv_frcng(k+1) !+ weighttop*wv_frcng(k) cig commented out else - wv_frcng(k+1) = 0.5*(wv_frcng(k+1) + wv_frcng(k)) + wv_frcng(k+1) = 0.5*(wv_frcng(k+1) + wv_frcng(k)) endif - + + + + diff_coeff(k) = (rho(i,j,iz0)/rbh)*fe*eps/(dz(k)* & bf(i,j,k)*bf(i,j,k)) - + !epg: following what we did above... + !cig: place the extra momentum flux in the top 3 layers instead of all in the top layer if ( k==0 ) then - diff_coeff(k+1) = 0.5*diff_coeff(k+1) + diff_coeff(k) + diff_coeff(k+1) = 0.5*diff_coeff(k+1) !+ weighttop*diff_coeff(k) cig commented out else diff_coeff(k+1) = 0.5*(diff_coeff(k+1) + diff_coeff(k)) endif - - else + + !cig: following what we did above... + + + else wv_frcng(iz0) = 0.0 diff_coeff(iz0) = 0.0 endif - end do ! (k loop) + end do ! (k loop) + +!cig: place the extra momentum flux in the layers above a specific threshold instead of all in the top layer (k=0 isn't a real model level) +! write (*,*) "iztop",iztop, ' ', damp_level(i,j), ' ', damp_level_pressure, ' ', iz0, ' ', source_level_pressure + do k=1,iztop + wv_frcng(k) = wv_frcng(k) + wv_frcng(0)/REAL(iztop) + diff_coeff(k) = diff_coeff(k) + diff_coeff(0)/REAL(iztop) + end do + +!cig: place the extra momentum flux in the top 3 layers instead of all in the top layer +! wv_frcng(1) = wv_frcng(1) + weighttop*wv_frcng(0) +! wv_frcng(2) = wv_frcng(2) + weightminus1*wv_frcng(0) +! wv_frcng(3) = wv_frcng(3) + weightminus2*wv_frcng(0) + +! diff_coeff(1) = diff_coeff(1) + weighttop*diff_coeff(0) +! diff_coeff(2) = diff_coeff(2) + weightminus1*diff_coeff(0) +! diff_coeff(3) = diff_coeff(3) + weightminus2*diff_coeff(0) !--------------------------------------------------------------------- ! increment the total forcing at each point with that obtained from ! the set of waves with the current wavenumber. !--------------------------------------------------------------------- - do k=0,iz0 + + + do k=0,iz0 gwf(i,j,k) = gwf(i,j,k) + wv_frcng(k) ked(i,j,k) = ked(i,j,k) + diff_coeff(k) - end do + end do end do ! wavelength loop - end do ! i loop - end do ! j loop + end do ! i loop + end do ! j loop !-------------------------------------------------------------------- @@ -1517,5 +1625,3 @@ end subroutine gwfc end module cg_drag_mod - - diff --git a/src/atmos_param/damping_driver/damping_driver.f90 b/src/atmos_param/damping_driver/damping_driver.f90 index f1829d5d3..6c3881595 100644 --- a/src/atmos_param/damping_driver/damping_driver.f90 +++ b/src/atmos_param/damping_driver/damping_driver.f90 @@ -44,7 +44,7 @@ module damping_driver_mod integer :: nlev_rayfric real :: sponge_pbottom = 50. ! [Pa] logical :: do_mg_drag = .false. -!epg: Use cg_drag.f90, GFDL's version of the Alexander and Dunkerton 1999 +!epg: Use cg_drag.f90, GFDL's version of the Alexander and Dunkerton 1999 ! Non-orographic gravity wave parameterization, updated as for Cohen et al. 2013 ! mj actively choose rayleigh friction logical :: do_rayleigh = .false. @@ -101,7 +101,7 @@ module damping_driver_mod real :: rfactr -! note: +! note: ! rfactr = coeff. for damping momentum at the top level character(len=128) :: version = '$Id: damping_driver.f90,v 10.0 2003/10/24 22:00:25 fms Exp $' @@ -119,7 +119,7 @@ subroutine damping_driver (is, js, lat, Time, delt, pfull, phalf, zfull, zhalf, u, v, t, q, r, udt, vdt, tdt, qdt, rdt, & ! mask, kbot) z_pbl, mask, kbot) - + !----------------------------------------------------------------------- integer, intent(in) :: is, js real, dimension(:,:), intent(in) :: lat @@ -265,6 +265,7 @@ subroutine damping_driver (is, js, lat, Time, delt, pfull, phalf, zfull, zhalf, !mj updating call to riga version of cg_drag !call cg_drag_calc (is, js, lat, pfull, zfull, t, u, Time, & ! delt, utnd) + call cg_drag_calc (is, js, lat, pfull, zfull, t, u, v, Time, delt, utnd, vtnd) udt = udt + utnd vdt = vdt + vtnd !mj @@ -274,9 +275,9 @@ subroutine damping_driver (is, js, lat, Time, delt, pfull, phalf, zfull, zhalf, if ( id_udt_cgwd > 0 ) then ! used = send_data ( id_udt_cgwd, utnd, Time, is, js, 1, & ! rmask=mask ) - used = send_data ( id_udt_cgwd, utnd, Time, rmask=mask ) + used = send_data ( id_udt_cgwd, utnd, Time) endif - + endif ! constant drag, modeled on Alexander-Dunkerton winter average @@ -298,7 +299,7 @@ subroutine damping_driver (is, js, lat, Time, delt, pfull, phalf, zfull, zhalf, *( -1.65*abs(lat)**3 +2.5*lat**2 +0.17*abs(lat) +const_drag_off ) end where enddo - udt = udt + utnd + udt = udt + utnd !----- diagnostics ----- @@ -366,12 +367,13 @@ end subroutine damping_driver !####################################################################### - subroutine damping_driver_init ( lonb, latb, pref, axes, Time, sgsmtn) + subroutine damping_driver_init (is, ie, js, je, lonb, latb, pref, axes, Time, sgsmtn) real, intent(in) :: lonb(:), latb(:), pref(:) - integer, intent(in) :: axes(4) + integer, dimension(4), intent(in) :: axes type(time_type), intent(in) :: Time real, dimension(:,:), intent(out) :: sgsmtn + integer, intent(in) :: is, ie, js, je !----------------------------------------------------------------------- ! lonb = longitude in radians of the grid box edges ! latb = latitude in radians of the grid box edges @@ -426,13 +428,13 @@ subroutine damping_driver_init ( lonb, latb, pref, axes, Time, sgsmtn) !----------------------------------------------------------------------- !----- mountain gravity wave drag ----- - if (do_mg_drag) call mg_drag_init (lonb, latb, sgsmtn) + if (do_mg_drag) call mg_drag_init (is, ie, js, je, lonb, latb, sgsmtn) !-------------------------------------------------------------------- !----- Alexander-Dunkerton gravity wave drag ----- - + if (do_cg_drag) then - call cg_drag_init (lonb, latb, pref, Time=Time, axes=axes) + call cg_drag_init (is, ie, js, je, lonb, latb, pref, Time=Time, axes=axes) endif !----------------------------------------------------------------------- @@ -454,7 +456,7 @@ subroutine damping_driver_init ( lonb, latb, pref, axes, Time, sgsmtn) register_diag_field ( mod_name, 'tdt_diss_rdamp', axes(1:3), Time, & 'Dissipative heating from Rayleigh damping',& 'deg_k/s', missing_value=missing_value ) - + id_diss_heat_rdamp = & register_diag_field ( mod_name, 'diss_heat_rdamp', axes(1:2), Time, & 'Integrated dissipative heating from Rayleigh damping',& @@ -499,7 +501,7 @@ subroutine damping_driver_init ( lonb, latb, pref, axes, Time, sgsmtn) register_diag_field ( mod_name, 'tdt_diss_gwd', axes(1:3), Time, & 'Dissipative heating from gravity wave drag',& 'deg_k/s', missing_value=missing_value ) - + id_diss_heat_gwd = & register_diag_field ( mod_name, 'diss_heat_gwd', axes(1:2), Time, & 'Integrated dissipative heating from gravity wave drag',& @@ -521,7 +523,7 @@ subroutine damping_driver_init ( lonb, latb, pref, axes, Time, sgsmtn) 'u wind tendency for constant drag', 'm/s2', & missing_value=missing_value ) endif - + !----------------------------------------------------------------------- !----- topo wave drag ----- @@ -639,4 +641,3 @@ end subroutine rayleigh !####################################################################### end module damping_driver_mod - diff --git a/src/atmos_param/mg_drag/mg_drag.f90 b/src/atmos_param/mg_drag/mg_drag.F90 similarity index 84% rename from src/atmos_param/mg_drag/mg_drag.f90 rename to src/atmos_param/mg_drag/mg_drag.F90 index 1279a852c..a6bd852a4 100644 --- a/src/atmos_param/mg_drag/mg_drag.f90 +++ b/src/atmos_param/mg_drag/mg_drag.F90 @@ -6,16 +6,22 @@ module mg_drag_mod !------------------------------------------------------------------- ! Calculates partial tendencies for the zonal and meridional winds -! due to the effect of mountain gravity wave drag +! due to the effect of mountain gravity wave drag !------------------------------------------------------------------- use topography_mod, only: get_topog_stdev use fms_mod, only: mpp_npes, field_size, file_exist, write_version_number, stdlog, & mpp_pe, mpp_root_pe, error_mesg, FATAL, NOTE, read_data, write_data, & - open_namelist_file, close_file, check_nml_error, open_restart_file, mpp_error + open_namelist_file, close_file, check_nml_error, open_restart_file, mpp_error, & + set_domain use fms_io_mod, only: get_restart_io_mode use constants_mod, only: Grav, Kappa, RDgas, cp_air +#ifdef COLUMN_MODEL + use spec_mpp_mod, only: grid_domain +#else + use transforms_mod, only: grid_domain +#endif !----------------------------------------------------------------------- implicit none @@ -47,13 +53,13 @@ module mg_drag_mod !--------------------------------------------------------------------- ! xl_mtn effective mountain length ( set currently to 100km) ! acoef order unity "tunable" parameter -! gmax order unity "tunable" parameter +! gmax order unity "tunable" parameter ! (may be enhanced to increase drag) ! rho stand value for density of the air at sea-level (1.13 KG/M**3) ! low_lev_frac - fraction of atmosphere (from bottom up) considered ! to be "low-level-layer for base flux calc. and where no ! wave breaking is allowed. -! flux_cut_level pressure level (Pa) above which flux divergence is set to zero +! flux_cut_level pressure level (Pa) above which flux divergence is set to zero !----------------------------------------------------------------------- real :: & @@ -65,22 +71,34 @@ module mg_drag_mod ,low_lev_frac = .23 real :: flux_cut_level= 0.0 +real :: tp_level = .9 + +integer :: ideal_lat_cent_NH = 48 +integer :: ideal_lat_cent_SH = 17 +integer :: ideal_pfull_cent = 19 +real :: ideal_strength = -0.00002 logical :: do_netcdf_restart = .true. logical :: do_conserve_energy = .false. logical :: do_mcm_mg_drag = .false. +logical :: no_strat_drag = .false. +logical :: no_trop_drag = .false. +logical :: do_idealised_NH = .false. +logical :: do_idealised_SH = .false. character(len=128) :: source_of_sgsmtn = 'input' namelist / mg_drag_nml / do_netcdf_restart, & xl_mtn, gmax, acoef, rho, low_lev_frac, & do_conserve_energy, do_mcm_mg_drag, & - source_of_sgsmtn, flux_cut_level + source_of_sgsmtn, flux_cut_level, & + no_strat_drag, no_trop_drag, tp_level, & + do_idealised_NH, do_idealised_SH, ideal_lat_cent_NH, ideal_lat_cent_SH, ideal_pfull_cent, ideal_strength public mg_drag, mg_drag_init, mg_drag_end contains -!############################################################################# +!############################################################################# subroutine mg_drag (is, js, delt, uwnd, vwnd, temp, pfull, phalf, & zfull,zhalf,dtaux,dtauy,dtemp,taubx, tauby, tausf,& @@ -129,14 +147,14 @@ subroutine mg_drag (is, js, delt, uwnd, vwnd, temp, pfull, phalf, & ! (dimensioned IDIM x JDIM)-kg/m/s**2 ! = -(RHO*U**3/(N*XL))*G(FR) FOR N**2 > 0 ! = 0 FOR N**2 <=0 -! DTAUX Tendency of the zonal wind component deceleration +! DTAUX Tendency of the zonal wind component deceleration ! (dimensioned IDIM x JDIM x KDIM) -! DTAUY Tendency of the meridional wind component deceleration +! DTAUY Tendency of the meridional wind component deceleration ! (dimensioned IDIM x JDIM x KDIM) ! dtemp Tendency of temperature due to dissipation of ke ! ! TAUSF = "CLIPPED" SAT MOMENTUM FLUX ( AT HALF LEVELS below top) -! +! !=================================================================== !----------------------------------------------------------------------- @@ -178,7 +196,7 @@ subroutine mg_drag (is, js, delt, uwnd, vwnd, temp, pfull, phalf, & !======================================================================= ! (Intent local) real , dimension(size(uwnd,1),size(uwnd,2)) :: xn, yn, psurf,ptop,taub - real , dimension(size(uwnd,1),size(uwnd,2),size(uwnd,3)) :: theta + real , dimension(size(uwnd,1),size(uwnd,2),size(uwnd,3)) :: theta real , dimension(size(uwnd,1),size(uwnd,2),size(uwnd,3)+1) :: taus real vsamp integer, dimension (size(uwnd,1),size(uwnd,2)) :: ktop, kbtm @@ -195,18 +213,18 @@ subroutine mg_drag (is, js, delt, uwnd, vwnd, temp, pfull, phalf, & ! -> -AETA(L)*PS*UMAG(L)*D(L)*GMAX/XL ! THETA POTENTIAL temperature at full model levels ! (dimensioned IDIM x JDIM x KDIM) -! PSURF Surface pressure +! PSURF Surface pressure ! (dimensioned IDIM x JDIM) ! PTOP Pressure at top of low-level layer ! (dimensioned IDIM x JDIM) ! KTOP Top model level index included in low-level layer ! (dimensioned IDIM x JDIM) ! KBTM Bottom model level index included in low-level layer -! usually the lowest level +! usually the lowest level ! (dimensioned IDIM x JDIM) !----------------------------------------------------------------------- ! type loop indicies - integer i, j, k, kd, kb, kt, kbp1, ktm1 + integer i, j, k, kd, kb, kt, kbp1, ktm1 !----------------------------------------------------------------------- ! Local variables needed only for code that ! implements supersource-like gravity wave drag. @@ -217,7 +235,7 @@ subroutine mg_drag (is, js, delt, uwnd, vwnd, temp, pfull, phalf, & real, dimension(size(uwnd,1),size(uwnd,2)) :: ulow, vlow, tlow, thlow real, dimension(size(uwnd,1),size(uwnd,2)) :: rlow, zsvar, bvfreq, x real, dimension(size(uwnd,1),size(uwnd,2)) :: depth, ave_p -integer, dimension(size(uwnd,1),size(uwnd,2)) :: ntop +integer, dimension(size(uwnd,1),size(uwnd,2)) :: ntop real, dimension(size(uwnd,1),size(uwnd,2),size(uwnd,3)) :: th, sh_ang, test real, dimension(size(uwnd,1),size(uwnd,2),size(uwnd,3)) :: sigma, del_sigma !real, dimension(size(uwnd,1),size(uwnd,2),size(uwnd,3)+1) :: sigma_half @@ -229,7 +247,6 @@ subroutine mg_drag (is, js, delt, uwnd, vwnd, temp, pfull, phalf, & kdim = size( uwnd, 3 ) kdimm1 = kdim - 1 kdimp1 = kdim + 1 - !----------------------------------------------------------------------- ! CODE VARIABLES DESCRIPTION @@ -259,7 +276,7 @@ subroutine mg_drag (is, js, delt, uwnd, vwnd, temp, pfull, phalf, & ! BNV,BNVK = "LOW-LEVEL",V. PROFILE - BRUNT VAISALA FREQ(1 ! (= N,N ! BNV2,BNVK2 = N**2, N(L)**2 -! HPRIME = Sub-grid scale mountain height +! HPRIME = Sub-grid scale mountain height ! over local domain (IDIM x JDIM) ! XL = EFFECTIVE MOUNTAIN LENGTH = (100KM EVERYWHERE) ! SIGTOP = HIGHEST LEVEL TO WHICH GRAVITY WAVE @@ -305,8 +322,8 @@ subroutine mg_drag (is, js, delt, uwnd, vwnd, temp, pfull, phalf, & ! (input via namelist), find highest model level. ptop(:,:) = (1.-low_lev_frac)*psurf(:,:) - do kd=kdim,1,-1 - where (pfull(:,:,kd) .ge. ptop(:,:)) + do kd=kdim,1,-1 + where (pfull(:,:,kd) .ge. ptop(:,:)) ktop(:,:) = kd end where end do @@ -328,7 +345,7 @@ subroutine mg_drag (is, js, delt, uwnd, vwnd, temp, pfull, phalf, & & xn,yn,taub,pfull, phalf,zfull,zhalf,vsamp,taus) ! calculate mountain gravity wave drag tendency contributions - call mgwd_tend (is,js,xn,yn,taub,phalf,taus,dtaux,dtauy, tausf) + call mgwd_tend (is,js,xn,yn,taub,pfull,phalf,taus,dtaux,dtauy,tausf) else if ( do_mcm_mg_drag ) then @@ -378,9 +395,9 @@ subroutine mg_drag (is, js, delt, uwnd, vwnd, temp, pfull, phalf, & do j = 1,jdim do i = 1,idim do k = ntop(i,j), kdim - ulow(i,j) = ulow(i,j) + del_sigma(i,j,k)* uwnd(i,j,k) - vlow(i,j) = vlow(i,j) + del_sigma(i,j,k)* vwnd(i,j,k) - tlow(i,j) = tlow(i,j) + del_sigma(i,j,k)* temp(i,j,k) + ulow(i,j) = ulow(i,j) + del_sigma(i,j,k)* uwnd(i,j,k) + vlow(i,j) = vlow(i,j) + del_sigma(i,j,k)* vwnd(i,j,k) + tlow(i,j) = tlow(i,j) + del_sigma(i,j,k)* temp(i,j,k) thlow(i,j) = thlow(i,j) + del_sigma(i,j,k)*th(i,j,k) depth(i,j) = depth(i,j) + del_sigma(i,j,k) end do @@ -440,7 +457,7 @@ subroutine mg_drag (is, js, delt, uwnd, vwnd, temp, pfull, phalf, & 20 continue if ( klast /= 0 ) test(i,j,1:klast) = 0.0 kcrit = klast + 1 - x(i,j) = x(i,j)/(1.-sigma(i,j,kcrit)) + x(i,j) = x(i,j)/(1.-sigma(i,j,kcrit)) ! should be x(i,j) = x(i,j)/(1 - sigma_half(klast)/sigma_half(kdim)) end do end do @@ -464,15 +481,20 @@ subroutine mg_drag (is, js, delt, uwnd, vwnd, temp, pfull, phalf, & dtemp = 0.0 endif +! Temperature tendencies crashing model due to size? Add a minimum size at x10-11. +where (dtemp<=1E-11) + dtemp = 0.0 +end where + return end subroutine mg_drag !======================================================================= -!############################################################################# - +!############################################################################# + subroutine mgwd_base_flux (is,js,uwnd,vwnd,temp,pfull,phalf,ktop,kbtm, & theta,xn,yn,taub) - + !------------------------------------------------------------------- @@ -493,12 +515,12 @@ subroutine mgwd_base_flux (is,js,uwnd,vwnd,temp,pfull,phalf,ktop,kbtm, & !======================================================================= ! (Intent local) real , dimension(size(uwnd,1),size(uwnd,2)) :: sumw, delp, ulow, bnv, & - & hprime, fr, g, ubar, vbar, bnv2 + & hprime, fr, g, ubar, vbar, bnv2 real grav2, xli, a, small integer idim, jdim,kdim,ie, je !----------------------------------------------------------------------- ! type loop indicies - integer i, j, k, kb, kt, kbp1, ktm1 + integer i, j, k, kb, kt, kbp1, ktm1 !----------------------------------------------------------------------- !=================================================================== @@ -516,7 +538,7 @@ subroutine mgwd_base_flux (is,js,uwnd,vwnd,temp,pfull,phalf,ktop,kbtm, & ! define local scalar variables xli=1.0/xl_mtn grav2=grav*grav - a = acoef + a = acoef !----------------------------------------------------------------------- @@ -565,7 +587,7 @@ subroutine mgwd_base_flux (is,js,uwnd,vwnd,temp,pfull,phalf,ktop,kbtm, & ! v197 uses p* as reference vlues for theta, in above 1000 hPa is used ! theta(:,:,:)=temp(:,:,:)*(pfull(:,:,:)/ & ! & phalf(:,:,kdim+1))**(-kappa) - + do j=1,jdim do i=1,idim kt=ktop(i,j) @@ -582,7 +604,7 @@ subroutine mgwd_base_flux (is,js,uwnd,vwnd,temp,pfull,phalf,ktop,kbtm, & ! ----------------------------------------- small = epsilon(ulow) - where (bnv2(:,:) .gt. 0.0) + where (bnv2(:,:) .gt. 0.0) bnv(:,:) = sqrt(bnv2(:,:)) fr (:,:) = bnv(:,:)*hprime(:,:)/(ulow(:,:) + small) g (:,:) = gmax*fr(:,:)*fr(:,:)/(fr(:,:)*fr(:,:)+a*a) @@ -596,7 +618,7 @@ subroutine mgwd_base_flux (is,js,uwnd,vwnd,temp,pfull,phalf,ktop,kbtm, & end subroutine mgwd_base_flux -!############################################################################# +!############################################################################# subroutine mgwd_satur_flux (uwnd,vwnd,temp,theta,ktop,kbtm, & xn,yn,taub,pfull,phalf,zfull,zhalf,vsamp,taus) @@ -606,7 +628,7 @@ subroutine mgwd_satur_flux (uwnd,vwnd,temp,theta,ktop,kbtm, & real, intent(in), dimension (:,:,:) :: & & uwnd, vwnd, temp, theta, pfull, phalf,zfull, zhalf real, intent(in), dimension (:,:) :: xn, yn, taub - real, intent(in) :: vsamp + real, intent(in) :: vsamp integer, intent(in), dimension (:,:) :: ktop, kbtm !=================================================================== ! Arguments (intent out) @@ -614,16 +636,16 @@ subroutine mgwd_satur_flux (uwnd,vwnd,temp,theta,ktop,kbtm, & !======================================================================= ! (Intent local) real , dimension(size(uwnd,1),size(uwnd,2),size(uwnd,3)) :: & - & dterm, dudz + & dterm, dudz real , dimension(size(uwnd,1),size(uwnd,2),size(uwnd,3)+1) :: & & umag, bnvk2, d,d2, d2i, d2udz2, extend real grav2, xli, small integer :: idim, jdim, kdim, kdimm1, kdimp1 !----------------------------------------------------------------------- ! type loop indicies - integer i, j, k, kb, kt, kbp1, ktm1 + integer i, j, k, kb, kt, kbp1, ktm1 !----------------------------------------------------------------------- -! type flux cutoff +! type flux cutoff integer kcut !======================================================================= @@ -740,7 +762,7 @@ subroutine mgwd_satur_flux (uwnd,vwnd,temp,theta,ktop,kbtm, & ! compute wkb extension term for umag > 0 ! --------------------------------------- - where (umag(:,:,:).gt.0.0) + where (umag(:,:,:).gt.0.0) extend(:,:,:) = vsamp*d2udz2(:,:,:)/umag(:,:,:) elsewhere extend(:,:,:) = 0.0 @@ -784,7 +806,7 @@ subroutine mgwd_satur_flux (uwnd,vwnd,temp,theta,ktop,kbtm, & ! initialize d2i to a large number, which will result in a very ! small vertical wavelength (d) where umag = 0. - where (umag(:,:,:).gt.0.0) + where (umag(:,:,:).gt.0.0) d2i(:,:,:) = (bnvk2(:,:,:)/(umag(:,:,:)* & & umag(:,:,:)) - extend(:,:,:) ) elsewhere @@ -797,7 +819,7 @@ subroutine mgwd_satur_flux (uwnd,vwnd,temp,theta,ktop,kbtm, & ! for 1/d**2 approaching 0 calculate d by dividing by a ! very small but finite number - where (d2i(:,:,:) .lt. 1.e-30) + where (d2i(:,:,:) .lt. 1.e-30) d(:,:,:) = 1.e+30 elsewhere d2(:,:,:) = 1./d2i(:,:,:) @@ -806,7 +828,7 @@ subroutine mgwd_satur_flux (uwnd,vwnd,temp,theta,ktop,kbtm, & ! set d=0 for umag=0. - where (umag(:,:,:).eq.0.0) + where (umag(:,:,:).eq.0.0) d(:,:,:) = 0.0 endwhere @@ -840,16 +862,16 @@ subroutine mgwd_satur_flux (uwnd,vwnd,temp,theta,ktop,kbtm, & ! keep taus profile constant across top model layer (original code) ! calculate taus profile in top model layer (new code) - + taus(:,:,1) = taus(:,:,2) !del taus(:,:,1) = -phalf(:,:,1)*umag(:,:,1)*umag(:,:,1) & !del & *d(:,:,1)*xli*gmax / (temp(:,:,1)*rdgas) ! do not allow wave breaking for unstable layers - + do k = 1,kdimp1 - where ( bnvk2(:,:,k) .lt. 0.0) + where ( bnvk2(:,:,k) .lt. 0.0) taus(:,:,k) = taub(:,:) endwhere end do @@ -858,8 +880,8 @@ subroutine mgwd_satur_flux (uwnd,vwnd,temp,theta,ktop,kbtm, & ! ------------------------------------------- ! tausat(:,:,1) = 0. ! use all forcing ! Instead, let remaining flux escape above flux_cut_level - if( flux_cut_level > 0.0 ) then - kcut= 1 + if( flux_cut_level > 0.0 ) then + kcut= 1 do while( phalf(1,1,kcut) < flux_cut_level ) kcut= kcut+1 enddo @@ -871,13 +893,13 @@ subroutine mgwd_satur_flux (uwnd,vwnd,temp,theta,ktop,kbtm, & end subroutine mgwd_satur_flux -!############################################################################# +!############################################################################# -subroutine mgwd_tend (is,js,xn,yn,taub,phalf,taus,dtaux,dtauy,tausf) +subroutine mgwd_tend (is,js,xn,yn,taub,pfull,phalf,taus,dtaux,dtauy,tausf) !=================================================================== ! Arguments (intent in) - real, intent(in), dimension (:,:,:) :: phalf, taus + real, intent(in), dimension (:,:,:) :: pfull,phalf, taus real, intent(in), dimension (:,:) :: xn, yn, taub integer, intent(in) :: is, js !=================================================================== @@ -885,9 +907,15 @@ subroutine mgwd_tend (is,js,xn,yn,taub,phalf,taus,dtaux,dtauy,tausf) real, intent(out), dimension (:,:,:) :: dtaux, dtauy, tausf !======================================================================= ! (Intent local) - real , dimension(size(phalf,1),size(phalf,2),size(phalf,3)) :: dterm + real , dimension(size(phalf,1),size(phalf,2),size(phalf,3)) :: dterm real , dimension(size(phalf,1),size(phalf,2),size(phalf,3)+1) :: taup integer kdim, kdimp1 + + ! (Intent local) (trop/strat calc) + real , dimension(size(dtaux,1),size(dtaux,2)) :: psurf,ptop_tp + integer, dimension (size(dtaux,1),size(dtaux,2)) :: kbtm + real pfull_av, ptop_av + integer ktop_tp !----------------------------------------------------------------------- ! type loop indicies integer k, kd @@ -919,32 +947,128 @@ subroutine mgwd_tend (is,js,xn,yn,taub,phalf,taus,dtaux,dtauy,tausf) ! <><><><><><><><><><><><><><><><><><><><><><><><><><><><> ! <><><><><><><><> DE-CELERATION CODE <><><><><><><><> - + ! CALCULATE DECELERATION TERMS - DTAUX,DTAUY ! ------------------------------------------ do k=1,kdim dterm(:,:,k) = grav*(taup (:,:,k+1)-taup (:,:,k)) & & /(phalf(:,:,k+1)-phalf(:,:,k)) - + dtaux(:,:,k) = xn(:,:)*dterm(:,:,k) dtauy(:,:,k) = yn(:,:)*dterm(:,:,k) end do - + ! ! RC Added to make stable at higher values of gmax (above 1) REMOVE AS UNREALISTIC!! + ! where (dtaux>0.001) + ! dtaux = 0.001 + ! end where + + ! where (dtaux<-0.001) + ! dtaux = -0.001 + ! end where ! print sample output ! print*, ' mgdrag output for i,j=', is,js -! print *,'taub = ', taub(is,js) -! print *,'taus = ', taus(is,js,:) -! print *,'taup = ', taup(is,js,:) +! print *,'taub = ', taub(is,js) +! print *,'taus = ', taus(is,js,:) +! print *,'taup = ', taup(is,js,:) -! *********************************************************** +!----------------------------------------------------------------------- +!! Turning off drag for stratosphere/troposphere +!----------------------------------------------------------------------- +! calculate bottom of low-leveles layer = lowest level unless kbot is present +! copied this code from the low level fraction calculation in mg_drag + ! if (present(kbot)) then + ! kbtm(:,:) = kbot(:,:) + ! else + ! kbtm(:,:) = kdim + ! endif + ! ! calculate top of low-level layer, first get surface p from phalf + ! if (present(kbot)) then + ! do j=1,jdim + ! do i=1,idim + ! psurf(i,j) = phalf(i,j,kbtm(i,j)+1) + ! end do + ! end do + ! else + ! psurf(:,:) = phalf(:,:,kdimp1) + ! endif + psurf(:,:)=pfull(:,:,kdim) +! Based on fraction of model atmosphere to be considered below the "tropopause" +! (input via namelist), find highest model level in the "troposphere". +! This definition is to do with the distribution of drag, not the actual tropopause + + ptop_tp(:,:) = (1.-tp_level)*psurf(:,:) + ptop_av=sum(ptop_tp(:,:))/size(ptop_tp(:,:)) + do kd=kdim,1,-1 + pfull_av=sum(pfull(:,:,kd))/size(pfull(:,:,kd)) + if (pfull_av .ge. ptop_av) then + ktop_tp = kd + endif + end do +! Turn off drag below the "tropopause" + if (no_trop_drag) then + do k=kdim,ktop_tp,-1 + dtaux(:,:,k) = 0 + dtauy(:,:,k) = 0 + end do + endif + +! Turn off drag above the "tropopause" + if (no_strat_drag) then + do k=ktop_tp,1,-1 + dtaux(:,:,k) = 0 + dtauy(:,:,k) = 0 + end do + endif + +!! Idealised Drag + if (do_idealised_NH) then + dtaux(:,:,:) = 0 + dtauy(:,:,:) = 0 + + if (ideal_lat_cent_NH .GE. js .AND. ideal_lat_cent_NH .LE. js+size(dtaux, 2 )-1) then + dtaux(:, MOD(ideal_lat_cent_NH-1,size(dtaux,2))+1, ideal_pfull_cent) = ideal_strength + dtaux(:, MOD(ideal_lat_cent_NH-1,size(dtaux,2))+1, ideal_pfull_cent-1) = ideal_strength + + if (MOD(ideal_lat_cent_NH,size(dtaux,2)) .NE. 0) then + dtaux(:, MOD(ideal_lat_cent_NH-1,size(dtaux,2))+2, ideal_pfull_cent) = ideal_strength + dtaux(:, MOD(ideal_lat_cent_NH-1,size(dtaux,2))+2, ideal_pfull_cent-1) = ideal_strength + endif + endif + + if (ideal_lat_cent_NH .EQ. js-1 ) then + dtaux(:,1,ideal_pfull_cent) = ideal_strength + dtaux(:,1,ideal_pfull_cent-1) = ideal_strength + endif + endif + + if (do_idealised_SH) then + dtaux(:,:,:) = 0 + dtauy(:,:,:) = 0 + + if (ideal_lat_cent_SH .GE. js .AND. ideal_lat_cent_SH .LE. js+size(dtaux, 2 )-1) then + dtaux(:, MOD(ideal_lat_cent_SH-1,size(dtaux,2))+1, ideal_pfull_cent) = ideal_strength + dtaux(:, MOD(ideal_lat_cent_SH-1,size(dtaux,2))+1, ideal_pfull_cent-1) = ideal_strength + + if (MOD(ideal_lat_cent_SH,size(dtaux,2)) .NE. 1) then + dtaux(:, MOD(ideal_lat_cent_SH-1,size(dtaux,2)), ideal_pfull_cent) = ideal_strength + dtaux(:, MOD(ideal_lat_cent_SH-1,size(dtaux,2)), ideal_pfull_cent-1) = ideal_strength + endif + endif + + if (ideal_lat_cent_SH .EQ. js+size(dtaux,2)) then + dtaux(:,size(dtaux,2),ideal_pfull_cent) = ideal_strength + dtaux(:,size(dtaux,2),ideal_pfull_cent-1) = ideal_strength + endif + endif +!----------------------------------------------------------------------- end subroutine mgwd_tend !####################################################################### - subroutine mg_drag_init( lonb, latb, hprime ) + subroutine mg_drag_init(is, ie, js, je, lonb, latb, hprime ) !======================================================================= ! ***** INITIALIZE Mountain Gravity Wave Drag @@ -955,14 +1079,15 @@ subroutine mg_drag_init( lonb, latb, hprime ) ! lonb = longitude in radians of the grid box edges ! latb = latitude in radians of the grid box edges !--------------------------------------------------------------------- + integer, intent(in) :: is, ie, js, je real, intent(in), dimension(:) :: lonb, latb - + !--------------------------------------------------------------------- ! Arguments (Intent out - optional) ! hprime = array of sub-grid scale mountain heights !--------------------------------------------------------------------- real, intent(out), dimension(:,:), optional :: hprime - + !--------------------------------------------------------------------- ! (Intent local) !--------------------------------------------------------------------- @@ -983,7 +1108,7 @@ subroutine mg_drag_init( lonb, latb, hprime ) unit = open_namelist_file() ierr = 1 do while( ierr .ne. 0 ) - read ( unit, nml = mg_drag_nml, iostat = io, end = 10 ) + read ( unit, nml = mg_drag_nml, iostat = io, end = 10 ) ierr = check_nml_error(io,'mg_drag_nml') end do 10 continue @@ -1007,8 +1132,8 @@ subroutine mg_drag_init( lonb, latb, hprime ) ix = size(lonb(:)) - 1 iy = size(latb(:)) - 1 - allocate( Ghprime(ix,iy) ) ; Ghprime = 0.0 - + allocate( Ghprime(is:ie,js:je) ) ; Ghprime = 0.0 + !------------------------------------------------------------------- module_is_initialized = .true. !--------------------------------------------------------------------- @@ -1025,7 +1150,8 @@ subroutine mg_drag_init( lonb, latb, hprime ) if ( file_exist( 'INPUT/mg_drag.res.nc' ) ) then if (mpp_pe() == mpp_root_pe()) call mpp_error ('mg_drag_mod', & 'Reading NetCDF formatted restart file: INPUT/mg_drag.res.nc', NOTE) - call read_data ('INPUT/mg_drag.res.nc', 'ghprime', Ghprime) + call set_domain(grid_domain) + call read_data ('INPUT/mg_drag.res.nc', 'ghprime', Ghprime, grid_domain) else if ( file_exist( 'INPUT/mg_drag.res' ) ) then if (mpp_pe() == mpp_root_pe()) call mpp_error ('mg_drag_mod', & 'Reading native formatted restart file.', NOTE) @@ -1041,9 +1167,13 @@ subroutine mg_drag_init( lonb, latb, hprime ) ' is not a valid value for source_of_sgsmtn', FATAL) endif + ! ! RC Added to make stable at higher values of gmax (above 1) 1500 STILL HAD ISSUES, TRY 1200 + ! where (Ghprime>700) + ! Ghprime = 700 + ! end where ! return sub-grid scale topography? if (present(hprime)) hprime = Ghprime - + !===================================================================== end subroutine mg_drag_init @@ -1056,7 +1186,8 @@ subroutine mg_drag_end if(do_netcdf_restart) then if (mpp_pe() == mpp_root_pe()) call mpp_error ('mg_drag_mod', & 'Writing NetCDF formatted restart file: RESTART/mg_drag.res.nc', NOTE) - call write_data('RESTART/mg_drag.res.nc', 'ghprime', ghprime) + call set_domain(grid_domain) + call write_data('RESTART/mg_drag.res.nc', 'ghprime', ghprime, grid_domain) else if (mpp_pe() == mpp_root_pe()) call mpp_error ('mg_drag_mod', & 'Writing native formatted restart file.', NOTE) diff --git a/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 b/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 index cfcbb4606..dc1cd44a8 100644 --- a/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 +++ b/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 @@ -653,7 +653,7 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l if(do_damping) then call pressure_variables(p_half_1d,ln_p_half_1d,pref(1:num_levels),ln_p_full_1d,PSTD_MKS) pref(num_levels+1) = PSTD_MKS - call damping_driver_init(rad_lonb_2d(:,1),rad_latb_2d(1,:), pref(:), get_axis_id(), Time, & + call damping_driver_init(is, ie, js, je, rad_lonb_2d(:,1),rad_latb_2d(1,:), pref(:), get_axis_id(), Time, & sgsmtn) !note that in the original this is pref(:,1), which is the full model pressure levels and !the surface pressure at the bottom. There is pref(:2) in this version with 81060 as diff --git a/src/extra/model/column/path_names b/src/extra/model/column/path_names index 6912d90f4..62758fe8b 100644 --- a/src/extra/model/column/path_names +++ b/src/extra/model/column/path_names @@ -105,7 +105,7 @@ atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_spcvmc.f90 atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_taumol.f90 atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_vrtqdr.f90 atmos_param/damping_driver/damping_driver.f90 -atmos_param/mg_drag/mg_drag.f90 +atmos_param/mg_drag/mg_drag.F90 atmos_param/cg_drag/cg_drag.f90 atmos_param/topo_drag/topo_drag.f90 atmos_param/vert_diff/vert_diff.F90 diff --git a/src/extra/model/dry/path_names b/src/extra/model/dry/path_names index f7caddc6c..890ac6372 100644 --- a/src/extra/model/dry/path_names +++ b/src/extra/model/dry/path_names @@ -25,7 +25,7 @@ atmos_param/frierson_monin_obukhov/frierson_monin_obukhov.F90 atmos_param/dry_convection/dry_convection.f90 atmos_param/rayleigh_bottom_drag/rayleigh_bottom_drag.F90 atmos_param/damping_driver/damping_driver.f90 -atmos_param/mg_drag/mg_drag.f90 +atmos_param/mg_drag/mg_drag.F90 atmos_param/cg_drag/cg_drag.f90 atmos_param/topo_drag/topo_drag.f90 atmos_param/vert_diff/vert_diff.F90 diff --git a/src/extra/model/grey/path_names b/src/extra/model/grey/path_names index 8edec5d9f..d5a045000 100644 --- a/src/extra/model/grey/path_names +++ b/src/extra/model/grey/path_names @@ -26,7 +26,7 @@ atmos_param/dry_convection/dry_convection.f90 atmos_param/ras/ras.f90 atmos_param/rayleigh_bottom_drag/rayleigh_bottom_drag.F90 atmos_param/damping_driver/damping_driver.f90 -atmos_param/mg_drag/mg_drag.f90 +atmos_param/mg_drag/mg_drag.F90 atmos_param/cg_drag/cg_drag.f90 atmos_param/topo_drag/topo_drag.f90 atmos_param/vert_diff/vert_diff.F90 diff --git a/src/extra/model/isca/path_names b/src/extra/model/isca/path_names index abcba6b80..74046e3f8 100644 --- a/src/extra/model/isca/path_names +++ b/src/extra/model/isca/path_names @@ -101,7 +101,7 @@ atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_spcvmc.f90 atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_taumol.f90 atmos_param/rrtm_radiation/rrtmg_sw/gcm_model/src/rrtmg_sw_vrtqdr.f90 atmos_param/damping_driver/damping_driver.f90 -atmos_param/mg_drag/mg_drag.f90 +atmos_param/mg_drag/mg_drag.F90 atmos_param/cg_drag/cg_drag.f90 atmos_param/topo_drag/topo_drag.f90 atmos_param/vert_diff/vert_diff.F90 diff --git a/src/extra/model/socrates/path_names b/src/extra/model/socrates/path_names index cab81e74d..f7762fb85 100644 --- a/src/extra/model/socrates/path_names +++ b/src/extra/model/socrates/path_names @@ -26,7 +26,7 @@ atmos_param/frierson_monin_obukhov/frierson_monin_obukhov.F90 atmos_param/dry_convection/dry_convection.f90 atmos_param/rayleigh_bottom_drag/rayleigh_bottom_drag.F90 atmos_param/damping_driver/damping_driver.f90 -atmos_param/mg_drag/mg_drag.f90 +atmos_param/mg_drag/mg_drag.F90 atmos_param/cg_drag/cg_drag.f90 atmos_param/topo_drag/topo_drag.f90 atmos_param/vert_diff/vert_diff.F90 diff --git a/src/extra/model/socrates_column/path_names b/src/extra/model/socrates_column/path_names index 20fae2c23..f3d3f4ef7 100644 --- a/src/extra/model/socrates_column/path_names +++ b/src/extra/model/socrates_column/path_names @@ -30,7 +30,7 @@ atmos_param/frierson_monin_obukhov/frierson_monin_obukhov.F90 atmos_param/dry_convection/dry_convection.f90 atmos_param/rayleigh_bottom_drag/rayleigh_bottom_drag.F90 atmos_param/damping_driver/damping_driver.f90 -atmos_param/mg_drag/mg_drag.f90 +atmos_param/mg_drag/mg_drag.F90 atmos_param/cg_drag/cg_drag.f90 atmos_param/topo_drag/topo_drag.f90 atmos_param/vert_diff/vert_diff.F90 diff --git a/src/extra/python/scripts/regrid_era5_sdor_to_t42.py b/src/extra/python/scripts/regrid_era5_sdor_to_t42.py new file mode 100644 index 000000000..a4e0890f9 --- /dev/null +++ b/src/extra/python/scripts/regrid_era5_sdor_to_t42.py @@ -0,0 +1,111 @@ +""" +Regrids ERA5 "standard deviation of sub-gridscale orography" (sdor) to Isca's +T42 Gaussian grid (64 lat x 128 lon, same grid as era-spectral7_T42_64x128.out.nc +elsewhere in this repo), producing a direct mg_drag "ghprime" input field. +Writes exp/test_cases/mg_drag/input/mg_drag.res.nc, the input file used by +exp/test_cases/mg_drag/mg_drag_socrates_topo_test_case.py. + +Source data: ERA5 single-level invariant parameter "Standard deviation of +orography" (sdor, ECMWF param 160), downloaded from the Copernicus Climate +Data Store (CDS) -- see https://cds.climate.copernicus.eu, reanalysis-era5-single-levels, +variable "standard_deviation_of_orography". Not committed to git (large binary, +regeneratable from CDS) -- edit ERA5_SDOR_PATH below to point at your own copy. + +Why 'input' mode instead of 'computed': mg_drag_nml's source_of_sgsmtn='computed' +path (src/atmos_param/mg_drag/mg_drag.F90 -> src/shared/topography/topography.F90) +expects *raw high-resolution elevation* and computes subgrid standard deviation +itself by binning many fine-grid points into each destination grid cell. ERA5's +sdor is already a precomputed standard-deviation field (ECMWF derives it the same +way, from a much finer DEM), so feeding it through 'computed' would compute the +"stdev of a stdev field" -- meaningless. Instead this script regrids sdor +directly onto the model grid and mg_drag reads it verbatim via +source_of_sgsmtn='input' (mg_drag_init's INPUT/mg_drag.res.nc restart-style +read, variable name 'ghprime'). + +Regridding: bilinear interpolation onto Isca's Gaussian T42 latitudes (via +numpy.polynomial.legendre.leggauss). Interpolation is done by hand with two +separable 1D numpy.interp passes (longitude then latitude) rather than +xarray's/scipy's interp -- this works cleanly here because both the ERA5 +source grid and the T42 target grid are regular in longitude and monotonic in +latitude, so a separable linear interpolation is exact for bilinear +interpolation on a rectilinear grid. Spherical-harmonic filtering is +deliberately not used: sdor is a standard-deviation field and must stay >= 0 +everywhere, and spectral truncation can ring negative near sharp gradients +(e.g. the Himalayas); linear interpolation of a non-negative field stays +non-negative. + +Latitude convention: mg_drag_init reads this file via FMS's read_data(), which +loads the raw array by index position into the model's internal field with no +coordinate-based matching -- it just trusts the array is already ordered to +match the model's own domain-decomposed grid, exactly like every other Isca +input file (era-spectral7_T42_64x128.out.nc is ascending, -90 -> 90). Writing +descending here would silently mirror the whole field north-south when +mg_drag reads it at runtime, so gaussian_grid() below returns ascending +latitudes to match. +""" +import os + +import numpy as np +import xarray as xr +from numpy.polynomial.legendre import leggauss + +from isca import GFDL_BASE + +ERA5_SDOR_PATH = '/home/links/sit204/gwd_tests/69eb30085e940a6d80b1ce305b746b0b.nc' +OUT_PATH = os.path.join(GFDL_BASE, 'exp/test_cases/mg_drag/input/mg_drag.res.nc') + +NLAT_T42 = 64 +NLON_T42 = 128 + + +def gaussian_grid(n_lat: int, n_lon: int) -> tuple[np.ndarray, np.ndarray]: + """Ascending (-90 -> 90) Gaussian latitudes, matching Isca's own grid + convention -- see the 'Latitude convention' note above for why this matters.""" + gauss_points, _ = leggauss(n_lat) + lats = np.arcsin(gauss_points) * (180.0 / np.pi) + lons = np.linspace(0, 360, n_lon, endpoint=False) + return lats, lons + + +def bilinear_regrid(data, lat_src, lon_src, lat_tgt, lon_tgt): + """Separable linear regrid (no scipy). lat_src must be ascending; + lon_src must be ascending and span >= [lon_tgt.min(), lon_tgt.max()].""" + # interpolate along longitude for every source latitude row + by_lon = np.empty((data.shape[0], len(lon_tgt))) + for i in range(data.shape[0]): + by_lon[i, :] = np.interp(lon_tgt, lon_src, data[i, :]) + # then interpolate along latitude for every target longitude column + out = np.empty((len(lat_tgt), len(lon_tgt))) + for j in range(len(lon_tgt)): + out[:, j] = np.interp(lat_tgt, lat_src, by_lon[:, j]) + return out + + +def main(): + era5 = xr.open_dataset(ERA5_SDOR_PATH).squeeze() + sdor = era5['sdor'].sortby('latitude') # ERA5 ships lat descending; need ascending for np.interp + + lat_src = sdor['latitude'].values + lon_src = sdor['longitude'].values + data_src = sdor.values + + target_lat, target_lon = gaussian_grid(NLAT_T42, NLON_T42) + regridded = bilinear_regrid(data_src, lat_src, lon_src, target_lat, target_lon) + + assert regridded.shape == (NLAT_T42, NLON_T42) + assert regridded.min() >= 0.0, 'regridded sdor went negative -- unexpected for linear interp of a >=0 field' + + ghprime = xr.DataArray( + regridded.astype('float64'), + dims=('lat', 'lon'), + coords={'lat': target_lat, 'lon': target_lon}, + name='ghprime', + attrs={'long_name': 'subgrid-scale orography standard deviation (from ERA5 sdor)', 'units': 'm'}, + ) + ghprime.to_netcdf(OUT_PATH) + print(f'Wrote {OUT_PATH}: ghprime min={float(ghprime.min()):.1f} max={float(ghprime.max()):.1f} ' + f'mean={float(ghprime.mean()):.2f} m, shape={ghprime.shape}') + + +if __name__ == '__main__': + main()