diff --git a/cime_config/namelist_definition_cam.xml b/cime_config/namelist_definition_cam.xml index 52d8d3318..7379936ba 100644 --- a/cime_config/namelist_definition_cam.xml +++ b/cime_config/namelist_definition_cam.xml @@ -377,6 +377,19 @@ + + + char*256 + gravity_wave_drag + gw_drag_input_nl + + Full pathname of boundary dataset for meso-gamma ridges. + + + UNSET_PATH + + + diff --git a/cime_config/testdefs/testlist_cam.xml b/cime_config/testdefs/testlist_cam.xml index da0b7dad9..b3e021999 100644 --- a/cime_config/testdefs/testlist_cam.xml +++ b/cime_config/testdefs/testlist_cam.xml @@ -88,6 +88,24 @@ + + + + + + + + + + + + + + + + + + diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq_gw_cam4_derecho/shell_commands b/cime_config/testdefs/testmods_dirs/cam/outfrq_gw_cam4_derecho/shell_commands new file mode 100644 index 000000000..0d4b3aa45 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq_gw_cam4_derecho/shell_commands @@ -0,0 +1 @@ + ./xmlchange CAM_CONFIG_OPTS="--dyn none --physics-suites gw_cam4" diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq_gw_cam4_derecho/user_nl_cam b/cime_config/testdefs/testmods_dirs/cam/outfrq_gw_cam4_derecho/user_nl_cam new file mode 100644 index 000000000..4f071621c --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq_gw_cam4_derecho/user_nl_cam @@ -0,0 +1,21 @@ +! cam4 snapshot +pver = 26 +ncdata = '/glade/campaign/cesm/community/amwg/sima_baselines/cam_sima_test_snapshots/cam_ne3pg3_fhistc4_gw_drag_cam4_oro_snapshot_derecho_gnu_before_c20250826.nc' +ncdata_check = '/glade/campaign/cesm/community/amwg/sima_baselines/cam_sima_test_snapshots/cam_ne3pg3_fhistc4_gw_drag_cam4_oro_snapshot_derecho_gnu_after_c20250826.nc' + +debug_output = 0 + +! tolerances for testing +ncdata_check_err = .true. +min_difference = 2e-15 + +effgw_oro = 0.125D0 +gw_apply_tndmax = .true. +gw_dc = 0.D0 +gw_lndscl_sgh = .true. +gw_oro_south_fac = 1.d0 +gw_prndl = 0.25d0 +gw_qbo_hdepth_scaling = 1.0D0 +gw_top_taper = .false. +pgwv = 0 +tau_0_ubc = .false. diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq_gw_cam7_derecho/shell_commands b/cime_config/testdefs/testmods_dirs/cam/outfrq_gw_cam7_derecho/shell_commands new file mode 100644 index 000000000..5c619d67a --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq_gw_cam7_derecho/shell_commands @@ -0,0 +1 @@ + ./xmlchange CAM_CONFIG_OPTS="--dyn none --physics-suites gw_cam7_se" diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq_gw_cam7_derecho/user_nl_cam b/cime_config/testdefs/testmods_dirs/cam/outfrq_gw_cam7_derecho/user_nl_cam new file mode 100644 index 000000000..6243140e6 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq_gw_cam7_derecho/user_nl_cam @@ -0,0 +1,43 @@ +! fhistc_ltso snapshot +pver = 58 +ncdata = '/glade/campaign/cesm/community/amwg/sima_baselines/cam_sima_test_snapshots/cam_ne3pg3_fhistc_ltso_gw_drag_cam7se_snapshot_derecho_gnu_before_c20251013.nc' +ncdata_check = '/glade/campaign/cesm/community/amwg/sima_baselines/cam_sima_test_snapshots/cam_ne3pg3_fhistc_ltso_gw_drag_cam7se_snapshot_derecho_gnu_after_c20251013.nc' + +debug_output = 0 + +! tolerances for testing +ncdata_check_err = .true. +min_difference = 2e-15 + +! topo file for beta ridge scheme input +bnd_topo = '/glade/campaign/cesm/cesmdata/inputdata/atm/cam/topo/se/ne3pg3_gmted2010_modis_bedmachine_nc0540_Laplace1000_noleak_20230209.nc' + +gw_movmtn_use_pbl_source = .true. ! se dycore +gw_movmtn_alpha = 0.01d0 +effgw_beres_dp = 0.15D0 +effgw_cm = 1.D0 +effgw_movmtn_pbl = 1.0d0 +effgw_rdg_beta = 0.5D0 +effgw_rdg_beta_max = 0.5D0 +front_gaussian_width = 30.D0 +gw_frontgfc = 1.25D-15 +gw_apply_tndmax = .false. +gw_dc = 2.5D0 +gw_dc_long = 0.D0 +gw_drag_file_dp = '/glade/campaign/cesm/cesmdata/inputdata/atm/waccm/gw/newmfspectra40_dc25.nc' +gw_drag_file_movmtn = '/glade/campaign/cesm/cesmdata/inputdata/atm/waccm/gw/mfc0lookup_mm.nc' +gw_lndscl_sgh = .true. +gw_oro_south_fac = 1.d0 +gw_polar_taper = .false. +gw_prndl = 0.5D0 +gw_qbo_hdepth_scaling = 0.32D0 +gw_movmtn_plaunch = 32500.0d0 +gw_movmtn_psteer = 65000.0d0 +gw_movmtn_source = 1 +gw_rdg_beta_nrdg = 10 +pgwv = 32 +pgwv_long = 0 +gw_rdg_beta_cd_llb = 1.0D0 +tau_0_ubc = .true. +gw_front_taubgnd = 1.5D-3 +gw_rdg_beta_trpd_leewv = .false. diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq_se_cslam_analy_ic/user_nl_cam b/cime_config/testdefs/testmods_dirs/cam/outfrq_se_cslam_analy_ic/user_nl_cam index ddd9f5800..69eb4b4f1 100644 --- a/cime_config/testdefs/testmods_dirs/cam/outfrq_se_cslam_analy_ic/user_nl_cam +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq_se_cslam_analy_ic/user_nl_cam @@ -2,3 +2,11 @@ hist_add_inst_fields;h1: T, Q, U, V, PS hist_output_frequency;h1: 1*nsteps hist_write_nstep0;h1: .true. analytic_ic_type=held_suarez_1994 + +! generally, analytic_ic runs do not need a topo file (see namelist_definition_cam.xml) +! however, a topo file needs to be supplied to run the Beta Ridge gravity wave scheme. +! +! in the future, the FCAM7 test can be migrated to using 'real' initial conditions, +! and this line can be removed, and +! analytic ic runs can go back to having no topo file supplied as it is by default. +bnd_topo='/glade/campaign/cesm/cesmdata/inputdata/atm/cam/topo/se/ne3pg3_gmted2010_modis_bedmachine_nc0540_Laplace1000_noleak_20230209.nc' diff --git a/src/control/cam_comp.F90 b/src/control/cam_comp.F90 index 3af070919..4c23eb524 100644 --- a/src/control/cam_comp.F90 +++ b/src/control/cam_comp.F90 @@ -103,6 +103,7 @@ subroutine cam_init(caseid, ctitle, model_doi_url, & use vert_coord, only: pver use phys_vars_init_check, only: mark_as_initialized use tropopause_climo_read, only: tropopause_climo_read_file + use gravity_wave_drag_ridge_read, only: gravity_wave_drag_ridge_read_file use orbital_data, only: orbital_data_init use ccpp_kinds, only: kind_phys use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t @@ -248,6 +249,9 @@ subroutine cam_init(caseid, ctitle, model_doi_url, & ! Read tropopause climatology call tropopause_climo_read_file() + ! Read gravity wave drag data for ridge parameterization + call gravity_wave_drag_ridge_read_file() + ! TEMPORARY: Prescribe realistic but inaccurate physical quantities ! necessary for MUSICA that are currently unavailable in CAM-SIMA. ! diff --git a/src/control/runtime_obj.F90 b/src/control/runtime_obj.F90 index 157545a8b..ff27e163e 100644 --- a/src/control/runtime_obj.F90 +++ b/src/control/runtime_obj.F90 @@ -20,10 +20,6 @@ module runtime_obj type, public :: runtime_options character(len=CS), private :: phys_suite = unset_str character(len=16), private :: waccmx_opt = unset_str - ! use_gw_front: Frontogenesis - logical, private :: use_gw_front = .false. - ! use_gw_front_igw: Frontogenesis to inertial spectrum. - logical, private :: use_gw_front_igw = .false. ! update_thermo_variables: update thermo "constants" to composition-dependent thermo variables logical, private :: update_thermo_variables = .false. contains @@ -33,8 +29,6 @@ module runtime_obj ! Runtime parameters of interest to dycore procedure, public :: waccmx_on procedure, public :: waccmx_option - procedure, public :: gw_front - procedure, public :: gw_front_igw procedure, public :: update_thermodynamic_variables end type runtime_options @@ -74,20 +68,6 @@ pure character(len=16) function waccmx_option(self) end function waccmx_option - pure logical function gw_front(self) - class(runtime_options), intent(in) :: self - - gw_front = self%use_gw_front - - end function gw_front - - pure logical function gw_front_igw(self) - class(runtime_options), intent(in) :: self - - gw_front_igw = self%use_gw_front_igw - - end function gw_front_igw - pure logical function update_thermodynamic_variables(self) class(runtime_options), intent(in) :: self @@ -95,14 +75,12 @@ pure logical function update_thermodynamic_variables(self) end function update_thermodynamic_variables - subroutine cam_set_runtime_opts(phys_suite, waccmx_opt, & - gw_front, gw_front_igw) + subroutine cam_set_runtime_opts(phys_suite, waccmx_opt) use cam_abortutils, only: endrun + ! Initialize the CAM runtime object character(len=CS), intent(in) :: phys_suite character(len=16), intent(in) :: waccmx_opt - logical, intent(in) :: gw_front - logical, intent(in) :: gw_front_igw if (runtime_configured) then ! We might need more action to reset this so do not allow it now @@ -111,8 +89,6 @@ subroutine cam_set_runtime_opts(phys_suite, waccmx_opt, & cam_runtime_opts%phys_suite = trim(phys_suite) cam_runtime_opts%waccmx_opt = trim(waccmx_opt) - cam_runtime_opts%use_gw_front = gw_front - cam_runtime_opts%use_gw_front_igw = gw_front_igw cam_runtime_opts%update_thermo_variables = (trim(waccmx_opt) == 'ionosphere' .or. & trim(waccmx_opt) == 'neutral') diff --git a/src/control/runtime_opts.F90 b/src/control/runtime_opts.F90 index 7eb79f322..2a88f1754 100644 --- a/src/control/runtime_opts.F90 +++ b/src/control/runtime_opts.F90 @@ -45,6 +45,7 @@ subroutine read_namelist(nlfilename, single_column, scmlat, scmlon) use inic_analytic_utils, only: analytic_ic_readnl use tropopause_climo_read, only: tropopause_climo_readnl + use gravity_wave_drag_ridge_read, only: gravity_wave_drag_ridge_read_readnl ! use tracers, only: tracers_readnl ! use nudging, only: nudging_readnl @@ -64,15 +65,11 @@ subroutine read_namelist(nlfilename, single_column, scmlat, scmlon) character(len=512) :: errmsg character(len=64), allocatable :: schemes(:) integer :: errflg - logical :: use_gw_front - logical :: use_gw_front_igw character(len=*), parameter :: subname = "read_namelist" ! Initialize system-wide runtime configuration variables waccmx_opt = unset_str - use_gw_front = .false. - use_gw_front_igw = .false. !----------------------------------------------------------------------- ! Call subroutines for modules to read their own namelist. @@ -104,7 +101,7 @@ subroutine read_namelist(nlfilename, single_column, scmlat, scmlon) call tropopause_climo_readnl(nlfilename) ! call scam_readnl(nlfilename, single_column, scmlat, scmlon) ! call nudging_readnl(nlfilename) - + call gravity_wave_drag_ridge_read_readnl(nlfilename) call dyn_readnl(nlfilename) ! Read the namelists for active physics schemes @@ -117,12 +114,8 @@ subroutine read_namelist(nlfilename, single_column, scmlat, scmlon) call cam_read_ccpp_scheme_namelists(nlfilename, schemes, & mpicom, masterprocid, masterproc, iulog) - ! XXgoldyXX: Need to figure out how to grab gravity-wave (and other) - ! scheme variables if and only if they are compiled in. - ! Finally, set the system-wide runtime configuration object - call cam_set_runtime_opts(phys_suite_name, waccmx_opt, & - use_gw_front, use_gw_front_igw) + call cam_set_runtime_opts(phys_suite_name, waccmx_opt) end subroutine read_namelist diff --git a/src/data/ref_pres.F90 b/src/data/ref_pres.F90 index bf4bda458..58c1c3e93 100644 --- a/src/data/ref_pres.F90 +++ b/src/data/ref_pres.F90 @@ -56,6 +56,12 @@ module ref_pres logical, protected :: do_molec_diff = .false. integer, protected :: nbot_molec = 0 + ! Pressure limit used to set gravity wave tapering at top of model (Pa) + real(kind_phys), parameter :: gravity_wave_taper_bot_press = 0.6E-02_kind_phys + + ! Bottom level for tapering gravity waves at top of model + integer, protected :: nbot_gravity_wave_top_taper = 0 + !==================================================================================== contains !==================================================================================== @@ -165,6 +171,11 @@ subroutine ref_pres_init(pref_edge_in, pref_mid_in, num_pr_lev_in) top=.false.) end if + ! Find level corresponding to the bottom for tapering gravity waves + ! at the top of model. + nbot_gravity_wave_top_taper = press_lim_idx(gravity_wave_taper_bot_press, & + top=.true.) + ! Tell rest of model that variables have been initialized: ! pref_edge_in call mark_as_initialized("reference_pressure_at_interface") @@ -185,9 +196,13 @@ subroutine ref_pres_init(pref_edge_in, pref_mid_in, num_pr_lev_in) ! molec_diff_bot_press call mark_as_initialized("pressure_at_bottom_of_molecular_diffusion") ! do_molec_diff - call mark_as_initialized("flag_for_molecular_diffusion") + call mark_as_initialized("do_molecular_diffusion") ! nbot_molec call mark_as_initialized("index_of_pressure_at_bottom_of_molecular_diffusion") + ! gravity_wave_taper_bot_press + call mark_as_initialized("largest_model_top_pressure_that_allows_tapering_gravity_wave_drag_at_model_top") + ! nbot_gravity_wave_top_taper + call mark_as_initialized("vertical_index_of_bottom_limit_for_tapering_gravity_wave_drag_at_model_top") end subroutine ref_pres_init diff --git a/src/data/ref_pres.meta b/src/data/ref_pres.meta index 79f9f3198..953c251f4 100644 --- a/src/data/ref_pres.meta +++ b/src/data/ref_pres.meta @@ -59,7 +59,7 @@ dimensions = () protected = True [ do_molec_diff ] - standard_name = flag_for_molecular_diffusion + standard_name = do_molecular_diffusion units = flag type = logical dimensions = () @@ -70,3 +70,14 @@ type = integer dimensions = () protected = True +[ gravity_wave_taper_bot_press ] + standard_name = largest_model_top_pressure_that_allows_tapering_gravity_wave_drag_at_model_top + units = Pa + type = real | kind = kind_phys + dimensions = () +[ nbot_gravity_wave_top_taper ] + standard_name = vertical_index_of_bottom_limit_for_tapering_gravity_wave_drag_at_model_top + units = index + type = integer + dimensions = () + protected = True diff --git a/src/data/registry.xml b/src/data/registry.xml index 9985d559a..ed41e7345 100644 --- a/src/data/registry.xml +++ b/src/data/registry.xml @@ -14,6 +14,7 @@ $SRCROOT/src/physics/utils/physics_grid.meta $SRCROOT/src/physics/utils/cam_constituents.meta $SRCROOT/src/physics/utils/tropopause_climo_read.meta + $SRCROOT/src/physics/utils/gravity_wave_drag_ridge_read.meta $SRCROOT/src/data/air_composition.meta $SRCROOT/src/data/cam_thermo.meta $SRCROOT/src/data/cam_thermo_formula.meta @@ -126,17 +127,24 @@ horizontal_dimension vertical_layer_dimension - frontgf pbuf_frontgf + frontgf pbuf_FRONTGF horizontal_dimension vertical_layer_dimension - frontga pbuf_frontga + frontga pbuf_FRONTGA + + + horizontal_dimension vertical_layer_dimension + pbuf_VORT4GW reciprocal_of_dimensionless_exner_function_wrt_surface_air_pressure frontogenesis_function frontogenesis_angle + relative_vorticity vertically_integrated_total_energy_using_dycore_energy_formula_at_start_of_physics_timestep vertically_integrated_total_energy_using_dycore_energy_formula vertically_integrated_total_water_at_start_of_physics_timestep @@ -1324,6 +1333,7 @@ units="None" type="cam_out_t"> + tpert pbuf_tpert + + + horizontal_dimension vertical_layer_dimension + 0.0_kind_phys + pbuf_TTEND_CLUBB + + + horizontal_dimension vertical_interface_dimension + 0.0_kind_phys + pbuf_UPWP_CLUBB_GW + + + horizontal_dimension vertical_interface_dimension + 0.0_kind_phys + pbuf_VPWP_CLUBB_GW + + + + + horizontal_dimension vertical_interface_dimension + 0.0_kind_phys + pbuf_kvt + + + horizontal_dimension + 0.0_kind_phys + pbuf_SGH + + should cloud fraction use shallow convection calculated convective cloud area fraction - + .false. + allocatable="parameter"> .true. + allocatable="parameter"> 1 horizontal_dimension vertical_layer_dimension number_of_ccpp_constituents 1 + + horizontal_dimension vertical_layer_dimension + 0.0_kind_phys + pbuf_TTEND_DP + + units="flag" type="logical" + allocatable="parameter"> .true. .false. diff --git a/src/dynamics/se/dp_coupling.F90 b/src/dynamics/se/dp_coupling.F90 index 572f663f0..de4f7fe06 100644 --- a/src/dynamics/se/dp_coupling.F90 +++ b/src/dynamics/se/dp_coupling.F90 @@ -53,6 +53,7 @@ subroutine d_p_coupling(cam_runtime_opts, phys_state, phys_tend, dyn_out) ! dry air mass. use gravity_waves_sources, only: gws_src_fnct + use gravity_waves_sources, only: gws_src_vort use hycoef, only: hyai, ps0 use test_fvm_mapping, only: test_mapping_overwrite_dyn_state, test_mapping_output_phys_state use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t @@ -64,6 +65,8 @@ subroutine d_p_coupling(cam_runtime_opts, phys_state, phys_tend, dyn_out) use time_mod, only: timelevel_qdp use control_mod, only: qsplit + use shr_kind_mod, only: shr_kind_cl + ! arguments type(runtime_options), intent(in) :: cam_runtime_opts ! Runtime settings object type(dyn_export_t), intent(inout) :: dyn_out ! dynamics export @@ -92,6 +95,9 @@ subroutine d_p_coupling(cam_runtime_opts, phys_state, phys_tend, dyn_out) real (kind=r8), allocatable :: frontgf(:,:,:) ! temp arrays to hold frontogenesis real (kind=r8), allocatable :: frontga(:,:,:) ! function (frontgf) and angle (frontga) + ! Vorticity + real (kind=r8), allocatable :: vort4gw(:,:,:) ! temp arrays to hold vorticity + integer :: ncols,ierr integer :: blk_ind(1), m, m_cnst integer :: nphys @@ -102,6 +108,7 @@ subroutine d_p_coupling(cam_runtime_opts, phys_state, phys_tend, dyn_out) character(len=*), parameter :: subname = 'd_p_coupling' character(len=200) :: stdname_test + character(len=shr_kind_cl) :: errmsg !---------------------------------------------------------------------------- @@ -116,9 +123,9 @@ subroutine d_p_coupling(cam_runtime_opts, phys_state, phys_tend, dyn_out) if (fv_nphys > 0) then nphys = fv_nphys else - allocate(qgll(np,np,nlev,num_advected), stat=ierr) + allocate(qgll(np,np,nlev,num_advected), stat=ierr, errmsg=errmsg) call check_allocate(ierr, subname, 'qgll(np,np,nlev,num_advected)', & - file=__FILE__, line=__LINE__) + file=__FILE__, line=__LINE__, errmsg=errmsg) nphys = np end if @@ -126,60 +133,59 @@ subroutine d_p_coupling(cam_runtime_opts, phys_state, phys_tend, dyn_out) const_data_ptr => cam_constituents_array() ! Allocate temporary arrays to hold data for physics decomposition - allocate(ps_tmp(nphys_pts,nelemd), stat=ierr) + allocate(ps_tmp(nphys_pts,nelemd), stat=ierr, errmsg=errmsg) call check_allocate(ierr, subname, 'ps_tmp(nphys_pts,nelemd)', & - file=__FILE__, line=__LINE__) + file=__FILE__, line=__LINE__, errmsg=errmsg) - allocate(dp3d_tmp(nphys_pts,pver,nelemd), stat=ierr) + allocate(dp3d_tmp(nphys_pts,pver,nelemd), stat=ierr, errmsg=errmsg) call check_allocate(ierr, subname, 'dp3d_tmp(nphys_pts,pver,nelemd)', & - file=__FILE__, line=__LINE__) + file=__FILE__, line=__LINE__, errmsg=errmsg) - allocate(dp3d_tmp_tmp(nphys_pts,pver), stat=ierr) + allocate(dp3d_tmp_tmp(nphys_pts,pver), stat=ierr, errmsg=errmsg) call check_allocate(ierr, subname, 'dp3d_tmp_tmp(nphys_pts,pver)', & - file=__FILE__, line=__LINE__) + file=__FILE__, line=__LINE__, errmsg=errmsg) - allocate(phis_tmp(nphys_pts,nelemd), stat=ierr) + allocate(phis_tmp(nphys_pts,nelemd), stat=ierr, errmsg=errmsg) call check_allocate(ierr, subname, 'phis_tmp(nphys_pts,nelemd)', & - file=__FILE__, line=__LINE__) + file=__FILE__, line=__LINE__, errmsg=errmsg) - allocate(T_tmp(nphys_pts,pver,nelemd), stat=ierr) + allocate(T_tmp(nphys_pts,pver,nelemd), stat=ierr, errmsg=errmsg) call check_allocate(ierr, subname, 'T_tmp(nphys_pts,pver,nelemd)', & - file=__FILE__, line=__LINE__) + file=__FILE__, line=__LINE__, errmsg=errmsg) - allocate(uv_tmp(nphys_pts,2,pver,nelemd), stat=ierr) + allocate(uv_tmp(nphys_pts,2,pver,nelemd), stat=ierr, errmsg=errmsg) call check_allocate(ierr, subname, 'uv_tmp(nphys_pts,2,pver,nelemd)', & - file=__FILE__, line=__LINE__) + file=__FILE__, line=__LINE__, errmsg=errmsg) - allocate(q_tmp(nphys_pts,pver,num_advected,nelemd), stat=ierr) + allocate(q_tmp(nphys_pts,pver,num_advected,nelemd), stat=ierr, errmsg=errmsg) call check_allocate(ierr, subname, 'q_tmp(nphys_pts,pver,num_advected,nelemd)', & - file=__FILE__, line=__LINE__) + file=__FILE__, line=__LINE__, errmsg=errmsg) - allocate(omega_tmp(nphys_pts,pver,nelemd), stat=ierr) + allocate(omega_tmp(nphys_pts,pver,nelemd), stat=ierr, errmsg=errmsg) call check_allocate(ierr, subname, 'omega_tmp(nphys_pts,pver,nelemd)', & - file=__FILE__, line=__LINE__) + file=__FILE__, line=__LINE__, errmsg=errmsg) - if (cam_runtime_opts%gw_front() .or. & - cam_runtime_opts%gw_front_igw()) then + allocate(frontgf(nphys_pts,pver,nelemd), stat=ierr, errmsg=errmsg) + call check_allocate(ierr, subname, 'frontgf(nphys_pts,pver,nelemd)', & + file=__FILE__, line=__LINE__, errmsg=errmsg) - allocate(frontgf(nphys_pts,pver,nelemd), stat=ierr) - call check_allocate(ierr, subname, 'frontgf(nphys_pts,pver,nelemd)', & - file=__FILE__, line=__LINE__) + allocate(frontga(nphys_pts,pver,nelemd), stat=ierr, errmsg=errmsg) + call check_allocate(ierr, subname, 'frontga(nphys_pts,pver,nelemd)', & + file=__FILE__, line=__LINE__, errmsg=errmsg) - allocate(frontga(nphys_pts,pver,nelemd), stat=ierr) - call check_allocate(ierr, subname, 'frontga(nphys_pts,pver,nelemd)', & - file=__FILE__, line=__LINE__) - end if + allocate(vort4gw(nphys_pts,pver,nelemd), stat=ierr, errmsg=errmsg) + call check_allocate(ierr, subname, 'vort4gw(nphys_pts,pver,nelemd)', & + file=__FILE__, line=__LINE__, errmsg=errmsg) - if (iam < par%nprocs) then - ! Gravity Waves - if (cam_runtime_opts%gw_front() .or. & - cam_runtime_opts%gw_front_igw()) then + if (iam < par%nprocs) then - ! Calculate frontogenesis function and angle - call gws_src_fnct(elem, tl_f, tl_qdp_np0, frontgf, frontga, nphys) + ! Calculate frontogenesis function and angle + ! for gravity wave parameterization. + call gws_src_fnct(elem, tl_f, tl_qdp_np0, frontgf, frontga, nphys) - end if + ! Calculate vorticity for moving mountain gravity wave parameterization. + call gws_src_vort(elem, tl_f, tl_qdp_np0, vort4gw, nphys) if (fv_nphys > 0) then call test_mapping_overwrite_dyn_state(elem,dyn_out%fvm) @@ -240,14 +246,9 @@ subroutine d_p_coupling(cam_runtime_opts, phys_state, phys_tend, dyn_out) omega_tmp(:,:,:) = 0._r8 phis_tmp(:,:) = 0._r8 q_tmp(:,:,:,:) = 0._r8 - - if (cam_runtime_opts%gw_front() .or. & - cam_runtime_opts%gw_front_igw()) then - - frontgf(:,:,:) = 0._r8 - frontga(:,:,:) = 0._r8 - - end if + frontgf(:,:,:) = 0._r8 + frontga(:,:,:) = 0._r8 + vort4gw(:,:,:) = 0._r8 endif ! iam < par%nprocs @@ -275,11 +276,9 @@ subroutine d_p_coupling(cam_runtime_opts, phys_state, phys_tend, dyn_out) phys_state%u(icol, ilyr) = real(uv_tmp(blk_ind(1), 1, ilyr, ie), kind_phys) phys_state%v(icol, ilyr) = real(uv_tmp(blk_ind(1), 2, ilyr, ie), kind_phys) phys_state%omega(icol, ilyr) = real(omega_tmp(blk_ind(1), ilyr, ie), kind_phys) - - if (cam_runtime_opts%gw_front() .or. cam_runtime_opts%gw_front_igw()) then - phys_state%frontgf(icol, ilyr) = real(frontgf(blk_ind(1), ilyr, ie), kind_phys) - phys_state%frontga(icol, ilyr) = real(frontga(blk_ind(1), ilyr, ie), kind_phys) - end if + phys_state%frontgf(icol, ilyr) = real(frontgf(blk_ind(1), ilyr, ie), kind_phys) + phys_state%frontga(icol, ilyr) = real(frontga(blk_ind(1), ilyr, ie), kind_phys) + phys_state%vorticity(icol, ilyr) = real(vort4gw(blk_ind(1), ilyr, ie), kind_phys) end do do m = 1, num_advected diff --git a/src/dynamics/se/dyn_comp.F90 b/src/dynamics/se/dyn_comp.F90 index bf5bfcab6..53d54af7e 100644 --- a/src/dynamics/se/dyn_comp.F90 +++ b/src/dynamics/se/dyn_comp.F90 @@ -838,7 +838,8 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) call prim_init2(elem, fvm, hybrid, nets, nete, TimeLevel, hvcoord) !$OMP END PARALLEL - if (cam_runtime_opts%gw_front() .or. cam_runtime_opts%gw_front_igw()) call gws_init(elem) + ! initialize gravity wave sources + call gws_init(elem) end if ! iam < par%nprocs !Remove/replace after CAMDEN history output is enabled -JN: @@ -1864,6 +1865,9 @@ subroutine read_inidat(dyn_in) call mark_as_initialized("tendency_of_eastward_wind_due_to_model_physics") call mark_as_initialized("tendency_of_northward_wind_due_to_model_physics") call mark_as_initialized("specific_heat_of_air_used_in_dycore") + call mark_as_initialized("frontogenesis_function") + call mark_as_initialized("frontogenesis_angle") + call mark_as_initialized("relative_vorticity") ! These energy variables are calculated by check_energy_timestep_init ! but need to be marked here @@ -2320,7 +2324,7 @@ subroutine read_dyn_field_2d(fieldname, fh, dimname, buffer) end if ! This code allows use of compiler option to set uninitialized values - ! to NaN. In that case cam_read_feild can return NaNs where the element + ! to NaN. In that case cam_read_field can return NaNs where the element ! GLL points are not "unique columns". ! Set NaNs or fillvalue points to zero: where (shr_infnan_isnan(buffer) .or. (buffer==fillvalue)) buffer = 0.0_r8 @@ -2350,7 +2354,7 @@ subroutine read_dyn_field_3d(fieldname, fh, dimname, buffer) end if ! This code allows use of compiler option to set uninitialized values - ! to NaN. In that case infld can return NaNs where the element GLL + ! to NaN. In that case cam_read_field can return NaNs where the element GLL ! points are not "unique columns". ! Set NaNs or fillvalue points to zero: where (shr_infnan_isnan(buffer) .or. (buffer==fillvalue)) buffer = 0.0_r8 diff --git a/src/dynamics/se/gravity_waves_sources.F90 b/src/dynamics/se/gravity_waves_sources.F90 index b3484c7cc..0591c3ef5 100644 --- a/src/dynamics/se/gravity_waves_sources.F90 +++ b/src/dynamics/se/gravity_waves_sources.F90 @@ -18,10 +18,12 @@ module gravity_waves_sources !! for use by WACCM (via dp_coupling) public :: gws_src_fnct + public :: gws_src_vort public :: gws_init private :: compute_frontogenesis + private :: compute_vorticity_4gw - type (EdgeBuffer_t) :: edge3 + type (EdgeBuffer_t) :: edge3,edge1 type (derivative_t) :: deriv real(r8) :: psurf_ref @@ -43,6 +45,7 @@ subroutine gws_init(elem) ! Set up variables similar to dyn_comp and prim_driver_mod initializations call initEdgeBuffer(par, edge3, elem, 3*nlev,nthreads=1) + call initEdgeBuffer(par, edge1, elem, nlev,nthreads=1) psurf_ref = hypi(nlevp) @@ -52,6 +55,7 @@ subroutine gws_src_fnct(elem, tl, tlq, frontgf, frontga,nphys) use vert_coord, only: pver use cam_abortutils, only: check_allocate + use shr_kind_mod, only: shr_kind_cl !SE dycore: use derivative_mod, only: derivinit @@ -70,8 +74,9 @@ subroutine gws_src_fnct(elem, tl, tlq, frontgf, frontga,nphys) ! Local variables type (hybrid_t) :: hybrid integer :: nets, nete, ithr, ncols, ie, iret - real(kind=r8), allocatable :: frontgf_thr(:,:,:,:) - real(kind=r8), allocatable :: frontga_thr(:,:,:,:) + real(kind=r8), allocatable :: frontgf_thr(:,:,:,:) + real(kind=r8), allocatable :: frontga_thr(:,:,:,:) + character(len=shr_kind_cl) :: errmsg character(len=*), parameter :: subname = 'gws_src_fnct' @@ -82,15 +87,15 @@ subroutine gws_src_fnct(elem, tl, tlq, frontgf, frontga,nphys) hybrid = config_thread_region(par,'serial') call get_loop_ranges(hybrid,ibeg=nets,iend=nete) - allocate(frontgf_thr(nphys,nphys,nlev,nets:nete), stat=iret) + allocate(frontgf_thr(nphys,nphys,nlev,nets:nete), stat=iret, errmsg=errmsg) call check_allocate(iret, subname, & 'frontgf_thr(nphys,nphys,nlev,nets:nete)', & - file=__FILE__, line=__LINE__) + file=__FILE__, line=__LINE__, errmsg=errmsg) - allocate(frontga_thr(nphys,nphys,nlev,nets:nete), stat=iret) + allocate(frontga_thr(nphys,nphys,nlev,nets:nete), stat=iret, errmsg=errmsg) call check_allocate(iret, subname, & 'frontga_thr(nphys,nphys,nlev,nets:nete)', & - file=__FILE__, line=__LINE__) + file=__FILE__, line=__LINE__, errmsg=errmsg) call compute_frontogenesis(frontgf_thr,frontga_thr,tl,tlq,elem,deriv,hybrid,nets,nete,nphys) if (fv_nphys>0) then @@ -111,6 +116,137 @@ subroutine gws_src_fnct(elem, tl, tlq, frontgf, frontga,nphys) end subroutine gws_src_fnct +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine gws_src_vort(elem, tl, tlq, vort4gw, nphys) + use derivative_mod, only : derivinit + use dimensions_mod, only : nelemd + use dof_mod, only : UniquePoints + use hybrid_mod, only : config_thread_region, get_loop_ranges + use parallel_mod, only : par + use vert_coord, only : pver + use thread_mod, only : horz_num_threads + use dimensions_mod, only : fv_nphys + use cam_abortutils, only : check_allocate + use shr_kind_mod, only : shr_kind_cl + + implicit none + type (element_t), intent(in), dimension(:) :: elem + integer, intent(in) :: tl, nphys, tlq + + ! + real (kind=r8), intent(out) :: vort4gw(nphys*nphys,pver,nelemd) + + ! Local variables + type (hybrid_t) :: hybrid + integer :: nets, nete, ithr, ncols, ie, ierr + + ! + real(kind=r8), allocatable :: vort4gw_thr(:,:,:,:) + character(len=shr_kind_cl) :: errmsg + + character(len=*), parameter :: subname = 'gws_src_vort' + + ! This does not need to be a thread private data-structure + call derivinit(deriv) + !!$OMP PARALLEL NUM_THREADS(horz_num_threads), DEFAULT(SHARED), PRIVATE(nets,nete,hybrid,ie,ncols,vort4gw_thr) + hybrid = config_thread_region(par,'serial') + call get_loop_ranges(hybrid,ibeg=nets,iend=nete) + + allocate(vort4gw_thr(nphys,nphys,nlev,nets:nete), stat=ierr, errmsg=errmsg) + call check_allocate(ierr, subname, & + 'vort4gw_thr(nphys,nphys,nlev,nets:nete)', & + file=__FILE__, line=__LINE__, errmsg=errmsg) + + call compute_vorticity_4gw(vort4gw_thr,tl,tlq,elem,deriv,hybrid,nets,nete,nphys) + + if (fv_nphys>0) then + do ie=nets,nete + vort4gw(:,:,ie) = RESHAPE(vort4gw_thr(:,:,:,ie),(/nphys*nphys,nlev/)) + end do + else + do ie=nets,nete + ncols = elem(ie)%idxP%NumUniquePts + call UniquePoints(elem(ie)%idxP, nlev, vort4gw_thr(:,:,:,ie), vort4gw(1:ncols,:,ie)) + end do + end if + deallocate(vort4gw_thr) + + !!$OMP END PARALLEL + + end subroutine gws_src_vort + + subroutine compute_vorticity_4gw(vort4gw,tl,tlq,elem,ederiv,hybrid,nets,nete,nphys) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! compute vorticity for use in gw params + ! F = ( curl ) [U,V] + ! + ! Original by Peter Lauritzen, Julio Bacmeister*, Dec 2024 + ! Patterned on 'compute_frontogenesis' + ! + ! * corresponding/blame-able + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + use derivative_mod, only: vorticity_sphere + use edge_mod, only: edgevpack, edgevunpack + use bndry_mod, only: bndry_exchange + use dimensions_mod, only: fv_nphys + use fvm_mapping, only: dyn2phys + + type(hybrid_t), intent(in) :: hybrid + type(element_t), intent(in) :: elem(:) + type(derivative_t), intent(in) :: ederiv + integer, intent(in) :: nets,nete,nphys + integer, intent(in) :: tl,tlq + real(r8), intent(out) :: vort4gw(nphys,nphys,nlev,nets:nete) + + ! local + real(r8) :: area_inv(fv_nphys,fv_nphys), tmp(np,np) + real(r8) :: vort_gll(np,np,nlev,nets:nete) + integer :: k,kptr,i,j,ie,component,h,nq,m_cnst,n0 + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! First calculate vorticity on GLL grid + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! set timelevel=1 for velocities + n0=tl + do ie=nets,nete + do k=1,nlev + call vorticity_sphere(elem(ie)%state%v(:,:,:,k,n0),ederiv,elem(ie),vort_gll(:,:,k,ie)) + end do + do k=1,nlev + vort_gll(:,:,k,ie) = vort_gll(:,:,k,ie)*elem(ie)%spheremp(:,:) + end do + ! pack + call edgeVpack(edge1, vort_gll(:,:,:,ie),nlev,0,ie) + end do + call bndry_exchange(hybrid,edge1,location='compute_vorticity_4gw') + do ie=nets,nete + call edgeVunpack(edge1, vort_gll(:,:,:,ie),nlev,0,ie) + ! apply inverse mass matrix, + do k=1,nlev + vort_gll(:,:,k,ie) = vort_gll(:,:,k,ie)*elem(ie)%rspheremp(:,:) + end do + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Now regrid from GLL to PhysGrid if necessary + ! otherwise just return vorticity on GLL grid + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + if (fv_nphys>0) then + tmp = 1.0_r8 + area_inv = dyn2phys(tmp,elem(ie)%metdet) + area_inv = 1.0_r8/area_inv + do k=1,nlev + vort4gw(:,:,k,ie) = dyn2phys( vort_gll(:,:,k,ie) , elem(ie)%metdet , area_inv ) + end do + else + do k=1,nlev + vort4gw(:,:,k,ie) = vort_gll(:,:,k,ie) + end do + end if + end do + + + end subroutine compute_vorticity_4gw + subroutine compute_frontogenesis(frontgf,frontga,tl,tlq,elem,ederiv,hybrid,nets,nete,nphys) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! compute frontogenesis function F @@ -125,6 +261,8 @@ subroutine compute_frontogenesis(frontgf,frontga,tl,tlq,elem,ederiv,hybrid,nets, ! Integrated into gravity_waves_sources module, several arguments made global ! to prevent repeated allocation/initialization ! + ! Frontogenesis function correction by Walter Hannah, Mark Taylor, and Jack Chen. October 2025 + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! use physconst, only: cappa use air_composition, only: dry_air_species_num,thermodynamic_active_species_num @@ -153,14 +291,45 @@ subroutine compute_frontogenesis(frontgf,frontga,tl,tlq,elem,ederiv,hybrid,nets, real(r8) :: frontga_gll(np,np,nlev,nets:nete) integer :: k,kptr,i,j,ie,component,h,nq,m_cnst real(r8) :: gradth(np,np,2,nlev,nets:nete) ! grad(theta) - real(r8) :: p(np,np) ! pressure at mid points - real(r8) :: pint(np,np) ! pressure at interface points - real(r8) :: theta(np,np) ! potential temperature at mid points + real(r8) :: p(np,np,nlev) ! pressure at mid points + real(r8) :: pint(np,np,nlev+1) ! pressure at interface points + real(r8) :: gradp(np,np,2) ! grad(pressure) + real(r8) :: theta(np,np,nlev) ! potential temperature at mid points + real(r8) :: dtheta_dp(np,np,nlev) ! d(theta)/dp for eta to pressure surface correction + real(r8) :: grad_wind_cart(np,np,2) ! horizontal gradient of zonal and meridional wind on cartesian coordinate on isobaric surface + real(r8) :: wind_cart(np,np,3,nlev) ! zonal & meridional wind on cartesian coordinate + real(r8) :: ddp_wind_cart(np,np,3,nlev) ! vertical gradient of zonal & meridional wind on cartesian coordinate real(r8) :: C(np,np,2), sum_water(np,np) + ! By Mark Taylor + ! For a vector velocity "v", a tensor "grad(v)", and a vector "grad(theta)", + ! this loop computes the vector "grad(theta)*grad(v)" + ! + ! Representing the tensor "grad(v)" in spherical coordinates is difficult. This routine + ! avoids this by computing a mathematically equivalent form using a mixture of + ! Cartesian and spherical coordinates + ! + ! This routine is a modified version of derivative_mod.F90:ugradv_sphere() in that the + ! grad(v) term is modified to compute grad_p(v) - the gradient on p-surfaces expressed + ! in terms of the gradient on model surfaces and a vertical pressure gradient. + ! + ! First, v is represented in cartesian coordinates v(c) for c=1,2,3 + ! For each v(c), we compute its gradient on p-surfaces via: + ! grad(v(c)) - d(v(c))/dz grad(p) + ! Each of these gradients is represented in *spherical* coordinates (i=1,2) + ! + ! We then dot each of these vectors with grad(theta). This dot product is computed + ! in spherical coordinates. The end result is wind_cart(c), for c=1,2,3 + ! These three scalars are the three Cartesian coefficients of + ! the vector "grad(theta)*grad(v)" + ! + ! This Cartesian vector is then transformed back to spherical coordinates + ! + do ie=nets,nete ! pressure at model top - pint(:,:) = hvcoord%hyai(1) + pint(:,:,1) = hvcoord%hyai(1)*hvcoord%ps0 + do k=1,nlev ! moist pressure at mid points sum_water(:,:) = 1.0_r8 @@ -171,15 +340,49 @@ subroutine compute_frontogenesis(frontgf,frontga,tl,tlq,elem,ederiv,hybrid,nets, ! sum_water(:,:) = sum_water(:,:) + elem(ie)%state%Qdp(:,:,k,m_cnst,tlq)/elem(ie)%state%dp3d(:,:,k,tl) end do - p(:,:) = pint(:,:) + 0.5_r8*sum_water(:,:)*elem(ie)%state%dp3d(:,:,k,tl) + p(:,:,k) = pint(:,:,k) + 0.5_r8*sum_water(:,:)*elem(ie)%state%dp3d(:,:,k,tl) ! moist pressure at interface for next iteration - pint(:,:) = pint(:,:)+elem(ie)%state%dp3d(:,:,k,tl) + pint(:,:,k+1) = pint(:,:,k)+elem(ie)%state%dp3d(:,:,k,tl) ! - theta(:,:) = elem(ie)%state%T(:,:,k,tl)*(psurf_ref / p(:,:))**cappa - ! gradth(:,:,:,k,ie) = gradient_sphere(theta,ederiv,elem(ie)%Dinv) - call gradient_sphere(theta,ederiv,elem(ie)%Dinv,gradth(:,:,:,k,ie)) - ! compute C = (grad(theta) dot grad ) u - C(:,:,:) = ugradv_sphere(gradth(:,:,:,k,ie), elem(ie)%state%v(:,:,:,k,tl),ederiv,elem(ie)) + theta(:,:,k) = elem(ie)%state%T(:,:,k,tl)*(psurf_ref / p(:,:,k))**cappa + end do + + call compute_vertical_derivative(pint,p,theta,dtheta_dp) + + do k=1,nlev + call gradient_sphere(theta(:,:,k),ederiv,elem(ie)%Dinv,gradth(:,:,:,k,ie)) + + call gradient_sphere(p(:,:,k),ederiv,elem(ie)%Dinv,gradp) + + do component=1,2 + gradth(:,:,component,k,ie) = gradth(:,:,component,k,ie) - dtheta_dp(:,:,k) * gradp(:,:,component) + end do + end do + + do k=1,nlev + do component=1,3 + wind_cart(:,:,component,k) = sum( elem(ie)%vec_sphere2cart(:,:,component,:) * elem(ie)%state%v(:,:,:,k,tl),3 ) + end do + end do + + do component=1,3 + call compute_vertical_derivative(pint,p,wind_cart(:,:,component,:),ddp_wind_cart(:,:,component,:)) + end do + do k=1,nlev + call gradient_sphere(p(:,:,k),ederiv,elem(ie)%Dinv,gradp) + + do component=1,3 + call gradient_sphere(wind_cart(:,:,component,k),ederiv,elem(ie)%Dinv,grad_wind_cart) + do i=1,2 + grad_wind_cart(:,:,i) = grad_wind_cart(:,:,i) - ddp_wind_cart(:,:,component,k) * gradp(:,:,i) + end do + wind_cart(:,:,component,k) = sum( gradth(:,:,:,k,ie) * grad_wind_cart , 3 ) + end do + + do i=1,2 + C(:,:,i) = sum(wind_cart(:,:,:,k)*elem(ie)%vec_sphere2cart(:,:,:,i), 3) + end do + ! gradth dot C frontgf_gll(:,:,k,ie) = -( C(:,:,1)*gradth(:,:,1,k,ie) + C(:,:,2)*gradth(:,:,2,k,ie) ) ! apply mass matrix @@ -231,5 +434,38 @@ subroutine compute_frontogenesis(frontgf,frontga,tl,tlq,elem,ederiv,hybrid,nets, enddo end subroutine compute_frontogenesis + subroutine compute_vertical_derivative(pint,pmid,data,ddata_dp) + !--------------------------------------------------------------------------- + real(r8), intent(in ) :: pint(np,np,nlev+1) + real(r8), intent(in ) :: pmid(np,np,nlev) + real(r8), intent(in ) :: data(np,np,nlev) + real(r8), intent(out) :: ddata_dp(np,np,nlev) + !--------------------------------------------------------------------------- + integer :: k + real(r8) :: pint_above(np,np) ! pressure interpolated to interface above the current k mid-point + real(r8) :: pint_below(np,np) ! pressure interpolated to interface below the current k mid-point + real(r8) :: dint_above(np,np) ! data interpolated to interface above the current k mid-point + real(r8) :: dint_below(np,np) ! data interpolated to interface below the current k mid-point + !--------------------------------------------------------------------------- + do k = 1,nlev + if (k==1) then + pint_above = pmid(:,:,k) + pint_below = pint(:,:,k+1) + dint_above = data(:,:,k) + dint_below = ( data(:,:,k+1) + data(:,:,k) ) / 2.0_r8 + elseif (k==nlev) then + pint_above = pint(:,:,k) + pint_below = pmid(:,:,k) + dint_above = ( data(:,:,k-1) + data(:,:,k) ) / 2.0_r8 + dint_below = data(:,:,k) + else + pint_above = pint(:,:,k) + pint_below = pint(:,:,k+1) + dint_above = ( data(:,:,k-1) + data(:,:,k) ) / 2.0_r8 + dint_below = ( data(:,:,k+1) + data(:,:,k) ) / 2.0_r8 + end if + ddata_dp(:,:,k) = ( dint_above - dint_below ) / ( pint_above - pint_below ) + end do + end subroutine compute_vertical_derivative end module gravity_waves_sources diff --git a/src/physics/utils/gravity_wave_drag_ridge_read.F90 b/src/physics/utils/gravity_wave_drag_ridge_read.F90 new file mode 100644 index 000000000..61ec571a1 --- /dev/null +++ b/src/physics/utils/gravity_wave_drag_ridge_read.F90 @@ -0,0 +1,360 @@ +! Support module for CAM-SIMA gravity_wave_drag_ridge parameterization +! to read in ridge data from the topo file. +! +! Remarks: this module is not CCPP-ized but is written specifically for +! grid decomposition-aware I/O in CAM-SIMA; it also has variables that +! are not a grid dimension (prdg = 16) that are read by specifying the +! custom dimension. It also makes that dimension available to the CCPP +! framework (prdg). +! This module can be a useful reference for how to provide +! gridded data to underlying CCPP schemes via the CAM PIO decomposition. +! (hplin, 8/28/25) +module gravity_wave_drag_ridge_read + use ccpp_kinds, only: kind_phys + use shr_kind_mod, only: shr_kind_cl + + implicit none + private + + public :: gravity_wave_drag_ridge_read_readnl + public :: gravity_wave_drag_ridge_read_file + + ! Topography file (rel pathname) for meso-Beta ridge data. + character(len=shr_kind_cl) :: bnd_topo = 'UNSET_PATH' + ! Resolved pathname. + character(len=shr_kind_cl) :: bnd_topo_loc = 'UNSET_PATH' + + ! Topography file (rel pathname) for meso-Gamma ridge data. + character(len=shr_kind_cl) :: bnd_rdggm = 'UNSET_PATH' + ! Resolved pathname. + character(len=shr_kind_cl) :: bnd_rdggm_loc = 'UNSET_PATH' + + ! Below data is provided externally to CCPP schemes. +!> \section arg_table_gravity_wave_drag_ridge_read Argument Table +!! \htmlinclude gravity_wave_drag_ridge_read.html + integer, parameter, public :: prdg = 16 ! # of ridges + + ! Meso-Beta ridges: + real(kind_phys), allocatable, public :: rdg_gbxar (:) ! Grid box area (ncol) [km2] + real(kind_phys), allocatable, public :: rdg_isovar(:) ! (B) Isotropic variance (ncol) [m] + real(kind_phys), allocatable, public :: rdg_isowgt(:) ! (B) Isotropic weight (ncol) [1] + real(kind_phys), allocatable, public :: rdg_hwdth (:, :) ! (B) Ridge Half-widths (ncol,prdg) [km] + real(kind_phys), allocatable, public :: rdg_clngt (:, :) ! (B) Ridge length (ncol,prdg) [km] + real(kind_phys), allocatable, public :: rdg_mxdis (:, :) ! (B) Ridge/obstacle height (ncol,prdg) [m] + real(kind_phys), allocatable, public :: rdg_anixy (:, :) ! (B) Ridge anisotropy (ncol,prdg) [1] + real(kind_phys), allocatable, public :: rdg_angll (:, :) ! (B) Ridge clockwise angle w.r.t. N-S direction (ncol,prdg) [degrees] + + ! Meso-Gamma ridges: + real(kind_phys), allocatable, public :: rdg_gbxarg(:) ! Grid box area (ncol) [km2] + real(kind_phys), allocatable, public :: rdg_hwdthg(:, :) ! (G) Ridge Half-widths (ncol,prdg) [km] + real(kind_phys), allocatable, public :: rdg_clngtg(:, :) ! (G) Ridge length (ncol,prdg) [km] + real(kind_phys), allocatable, public :: rdg_mxdisg(:, :) ! (G) Ridge/obstacle height (ncol,prdg) [m] + real(kind_phys), allocatable, public :: rdg_anixyg(:, :) ! (G) Ridge anisotropy (ncol,prdg) [1] + real(kind_phys), allocatable, public :: rdg_angllg(:, :) ! (G) Ridge clockwise angle w.r.t. N-S direction (ncol,prdg) [degrees] + +contains + + subroutine gravity_wave_drag_ridge_read_readnl(nlfile) + use shr_nl_mod, only: find_group_name => shr_nl_find_group_name + use shr_kind_mod, only: shr_kind_cm + use mpi, only: mpi_character + use spmd_utils, only: mpicom + use cam_logfile, only: iulog + use cam_abortutils, only: endrun + use spmd_utils, only: masterproc + use cam_initfiles, only: unset_path_str + + ! topography dataset is available in this module: + use cam_initfiles, only: cam_bnd_topo => bnd_topo + ! and does not need to be re-read from namelist. + + ! filepath for file containing namelist input + character(len=*), intent(in) :: nlfile + + ! Local variables + integer :: unitn, errflg + character(len=*), parameter :: subname = 'gravity_wave_drag_ridge_read_readnl' + character(len=shr_kind_cm) :: errmsg + + namelist /gw_drag_input_nl/ bnd_rdggm + + errmsg = '' + errflg = 0 + + if (masterproc) then + open(newunit=unitn, file=trim(nlfile), status='old') + call find_group_name(unitn, 'gw_drag_input_nl', status=errflg) + if (errflg == 0) then + read(unitn, gw_drag_input_nl, iostat=errflg, iomsg=errmsg) + if (errflg /= 0) then + call endrun(subname // ':: ERROR reading namelist:' // errmsg) + end if + end if + close(unitn) + end if + + ! Broadcast namelist variables + call mpi_bcast(bnd_rdggm, len(bnd_rdggm), mpi_character, 0, mpicom, errflg) + + ! Retrieve topo file location: + bnd_topo = cam_bnd_topo + + ! Print out namelist variables + if (masterproc) then + write(iulog,*) subname, ' options:' + if(bnd_topo /= unset_path_str) then + write(iulog,*) ' Gravity wave meso-Beta ridge topo file: ', trim(bnd_topo) + else + write(iulog,*) ' Gravity wave meso-Beta ridge input data unavailable.' + endif + + if(bnd_rdggm /= unset_path_str) then + write(iulog,*) ' Gravity wave meso-Gamma ridge input file: ', trim(bnd_rdggm) + else + write(iulog,*) ' Gravity wave meso-Gamma ridge input data unavailable.' + endif + endif + end subroutine gravity_wave_drag_ridge_read_readnl + + subroutine gravity_wave_drag_ridge_read_file() + use spmd_utils, only: masterproc + use cam_logfile, only: iulog + use cam_abortutils, only: endrun, check_allocate + use pio, only: file_desc_t, pio_nowrite + use cam_pio_utils, only: cam_pio_openfile, cam_pio_closefile + use ioFileMod, only: cam_get_file + use cam_initfiles, only: topo_file_get_id + use cam_field_read, only: cam_read_field + use physconst, only: rearth + use cam_initfiles, only: unset_path_str + use physics_grid, only: ncol => columns_on_task + use phys_vars_init_check, only: mark_as_initialized + + ! Local variables + type(file_desc_t), pointer :: fh_topo + type(file_desc_t), pointer :: fh_rdggm + integer :: errflg + character(len=512) :: errmsg + character(len=*), parameter :: subname = 'gravity_wave_drag_ridge_read_file' + + logical :: found + logical :: has_gbxar_from_topo + + errmsg = '' + errflg = 0 + + has_gbxar_from_topo = .false. + call mark_as_initialized('number_of_ridges_in_ridge_gravity_wave_drag') + + ! Do we have meso-Beta file? + if(bnd_topo /= unset_path_str) then + call cam_get_file(bnd_topo, bnd_topo_loc) + + ! Try getting from initfiles first. + fh_topo => topo_file_get_id() + if(.not. associated(fh_topo)) then + ! I think this case will never be hit in SIMA. + ! There was a fallback in CAM. + call endrun(trim(subname) // ": fh_topo from cam_initfiles is not available and this is not implemented.") + endif + + if(masterproc) then + write (iulog,*) trim(subname)//': Reading meso-Beta ridge data from ', trim(bnd_topo_loc) + endif + + ! Allocate and initialize data to zeros. + allocate(rdg_gbxar(ncol), stat=errflg, errmsg=errmsg) + call check_allocate(errflg, subname, 'rdg_gbxar', errmsg=errmsg) + allocate(rdg_isovar(ncol), stat=errflg, errmsg=errmsg) + call check_allocate(errflg, subname, 'rdg_isovar', errmsg=errmsg) + allocate(rdg_isowgt(ncol), stat=errflg, errmsg=errmsg) + call check_allocate(errflg, subname, 'rdg_isowgt', errmsg=errmsg) + allocate(rdg_hwdth(ncol, prdg), stat=errflg, errmsg=errmsg) + call check_allocate(errflg, subname, 'rdg_hwdth', errmsg=errmsg) + allocate(rdg_clngt(ncol, prdg), stat=errflg, errmsg=errmsg) + call check_allocate(errflg, subname, 'rdg_clngt', errmsg=errmsg) + allocate(rdg_mxdis(ncol, prdg), stat=errflg, errmsg=errmsg) + call check_allocate(errflg, subname, 'rdg_mxdis', errmsg=errmsg) + allocate(rdg_anixy(ncol, prdg), stat=errflg, errmsg=errmsg) + call check_allocate(errflg, subname, 'rdg_anixy', errmsg=errmsg) + allocate(rdg_angll(ncol, prdg), stat=errflg, errmsg=errmsg) + call check_allocate(errflg, subname, 'rdg_angll', errmsg=errmsg) + + rdg_gbxar(:) = 0._kind_phys + rdg_isovar(:) = 0._kind_phys + rdg_isowgt(:) = 0._kind_phys + rdg_hwdth(:,:) = 0._kind_phys + rdg_clngt(:,:) = 0._kind_phys + rdg_mxdis(:,:) = 0._kind_phys + rdg_anixy(:,:) = 0._kind_phys + rdg_angll(:,:) = 0._kind_phys + + ! Use cam_field_read to read the file data: + + ! Read required 1D field: GBXAR (grid box area) + call cam_read_field('GBXAR', fh_topo, rdg_gbxar, found) + if(.not. found) then + call endrun(trim(subname) // ': GBXAR not found in input file') + endif + ! Convert from m2 to km2 + rdg_gbxar = rdg_gbxar * (rearth/1000._kind_phys) * (rearth/1000._kind_phys) + has_gbxar_from_topo = .true. + + ! Read optional 1D field: ISOVAR (isotropic variance) + call cam_read_field('ISOVAR', fh_topo, rdg_isovar, found) + if(.not. found) then + if(masterproc) then + write(iulog,*) trim(subname) // ': ISOVAR not found in topo file, using zero values' + endif + ! rdg_isovar already initialized to zero above + endif + + ! Read optional 1D field: ISOWGT (isotropic weight) + call cam_read_field('ISOWGT', fh_topo, rdg_isowgt, found) + if(.not. found) then + if(masterproc) then + write(iulog,*) trim(subname) // ': ISOWGT not found in topo file, using zero values' + endif + ! rdg_isowgt already initialized to zero above + endif + + ! Read required 2D field: HWDTH (ridge half-widths) + call cam_read_field('HWDTH', fh_topo, rdg_hwdth, found, dim3name='nrdg', dim3_bnds=(/1, prdg/)) + if(.not. found) then + call endrun(trim(subname) // ': HWDTH not found in input file') + endif + + ! Read required 2D field: CLNGT (ridge length) + call cam_read_field('CLNGT', fh_topo, rdg_clngt, found, dim3name='nrdg', dim3_bnds=(/1, prdg/)) + if(.not. found) then + call endrun(trim(subname) // ': CLNGT not found in input file') + endif + + ! Read required 2D field: MXDIS (ridge/obstacle height) + call cam_read_field('MXDIS', fh_topo, rdg_mxdis, found, dim3name='nrdg', dim3_bnds=(/1, prdg/)) + if(.not. found) then + call endrun(trim(subname) // ': MXDIS not found in input file') + endif + + ! Read required 2D field: ANIXY (ridge anisotropy) + call cam_read_field('ANIXY', fh_topo, rdg_anixy, found, dim3name='nrdg', dim3_bnds=(/1, prdg/)) + if(.not. found) then + call endrun(trim(subname) // ': ANIXY not found in input file') + endif + + ! Read required 2D field: ANGLL (ridge clockwise angle w.r.t. N-S direction) + call cam_read_field('ANGLL', fh_topo, rdg_angll, found, dim3name='nrdg', dim3_bnds=(/1, prdg/)) + if(.not. found) then + call endrun(trim(subname) // ': ANGLL not found in input file') + endif + + ! Mark variables as initialized so they are not read from ic file. + call mark_as_initialized('grid_box_area_for_beta_ridge_gravity_wave_drag') + call mark_as_initialized('isotropic_variance_for_beta_ridge_gravity_wave_drag') + call mark_as_initialized('isotropic_weight_for_beta_ridge_gravity_wave_drag') + call mark_as_initialized('ridge_half_width_for_beta_ridge_gravity_wave_drag') + call mark_as_initialized('ridge_length_for_beta_ridge_gravity_wave_drag') + call mark_as_initialized('ridge_obstacle_height_for_beta_ridge_gravity_wave_drag') + call mark_as_initialized('ridge_anisotropy_for_beta_ridge_gravity_wave_drag') + call mark_as_initialized('ridge_clockwise_angle_from_north_for_beta_ridge_gravity_wave_drag') + endif + + ! Do we have meso-Gamma file? + if(bnd_rdggm /= unset_path_str) then + call cam_get_file(bnd_rdggm, bnd_rdggm_loc) + call cam_pio_openfile(fh_rdggm, bnd_rdggm_loc, pio_nowrite) + + if(masterproc) then + write (iulog,*) trim(subname)//': Reading meso-Gamma ridge data from ', trim(bnd_topo_loc) + endif + + ! Allocate meso-Gamma ridge data arrays + allocate(rdg_gbxarg(ncol), stat=errflg, errmsg=errmsg) + call check_allocate(errflg, subname, 'rdg_gbxarg', errmsg=errmsg) + allocate(rdg_hwdthg(ncol, prdg), stat=errflg, errmsg=errmsg) + call check_allocate(errflg, subname, 'rdg_hwdthg', errmsg=errmsg) + allocate(rdg_clngtg(ncol, prdg), stat=errflg, errmsg=errmsg) + call check_allocate(errflg, subname, 'rdg_clngtg', errmsg=errmsg) + allocate(rdg_mxdisg(ncol, prdg), stat=errflg, errmsg=errmsg) + call check_allocate(errflg, subname, 'rdg_mxdisg', errmsg=errmsg) + allocate(rdg_anixyg(ncol, prdg), stat=errflg, errmsg=errmsg) + call check_allocate(errflg, subname, 'rdg_anixyg', errmsg=errmsg) + allocate(rdg_angllg(ncol, prdg), stat=errflg, errmsg=errmsg) + call check_allocate(errflg, subname, 'rdg_angllg', errmsg=errmsg) + + ! Initialize all gamma fields to zero first + rdg_gbxarg(:) = 0._kind_phys + rdg_hwdthg(:,:) = 0._kind_phys + rdg_clngtg(:,:) = 0._kind_phys + rdg_mxdisg(:,:) = 0._kind_phys + rdg_anixyg(:,:) = 0._kind_phys + rdg_angllg(:,:) = 0._kind_phys + + if(.not. has_gbxar_from_topo) then + call cam_read_field('GBXAR', fh_rdggm, rdg_gbxarg, found) + if(.not. found) then + call endrun(trim(subname) // ': GBXAR not found in neither topo or gamma ridge file') + endif + ! Convert from m2 to km2 + rdg_gbxarg = rdg_gbxarg * (rearth/1000._kind_phys) * (rearth/1000._kind_phys) + else + if(masterproc) then + write(iulog,*) trim(subname) // ': Using GBXAR from topo file, skipping gamma file GBXAR' + endif + endif + + ! ISOVAR and ISOWGT are intentionally not read from gamma file, + ! because (1) they are unavailable, and (2) the original code did not read them in + ! and left them as dangling pointers. See gw_drag_cam CAM interface, + ! gravity_wave_drag_ridge_gamma_run implementation. + + ! Read required 2D field: HWDTH (ridge half-widths gamma) + call cam_read_field('HWDTH', fh_rdggm, rdg_hwdthg, found, dim3name='nrdg', dim3_bnds=(/1, prdg/)) + if(.not. found) then + call endrun(trim(subname) // ': HWDTH not found in gamma ridge file') + endif + + ! Read required 2D field: CLNGT (ridge length gamma) + call cam_read_field('CLNGT', fh_rdggm, rdg_clngtg, found, dim3name='nrdg', dim3_bnds=(/1, prdg/)) + if(.not. found) then + call endrun(trim(subname) // ': CLNGT not found in gamma ridge file') + endif + + ! Read required 2D field: MXDIS (ridge/obstacle height gamma) + call cam_read_field('MXDIS', fh_rdggm, rdg_mxdisg, found, dim3name='nrdg', dim3_bnds=(/1, prdg/)) + if(.not. found) then + call endrun(trim(subname) // ': MXDIS not found in gamma ridge file') + endif + + ! Apply negative value correction for gamma ridge maximum displacement + where (rdg_mxdisg < 0._kind_phys) + rdg_mxdisg = 0._kind_phys + end where + + ! Read required 2D field: ANIXY (ridge anisotropy gamma) + call cam_read_field('ANIXY', fh_rdggm, rdg_anixyg, found, dim3name='nrdg', dim3_bnds=(/1, prdg/)) + if(.not. found) then + call endrun(trim(subname) // ': ANIXY not found in gamma ridge file') + endif + + ! Read required 2D field: ANGLL (ridge clockwise angle w.r.t. N-S direction gamma) + call cam_read_field('ANGLL', fh_rdggm, rdg_angllg, found, dim3name='nrdg', dim3_bnds=(/1, prdg/)) + if(.not. found) then + call endrun(trim(subname) // ': ANGLL not found in gamma ridge file') + endif + + call cam_pio_closefile(fh_rdggm) + deallocate(fh_rdggm) + nullify(fh_rdggm) + + ! Mark variables as initialized so they are not read from ic file. + call mark_as_initialized('grid_box_area_for_gamma_ridge_gravity_wave_drag') + call mark_as_initialized('ridge_half_width_for_gamma_ridge_gravity_wave_drag') + call mark_as_initialized('ridge_length_for_gamma_ridge_gravity_wave_drag') + call mark_as_initialized('ridge_obstacle_height_for_gamma_ridge_gravity_wave_drag') + call mark_as_initialized('ridge_anisotropy_for_gamma_ridge_gravity_wave_drag') + call mark_as_initialized('ridge_clockwise_angle_from_north_for_gamma_ridge_gravity_wave_drag') + endif + end subroutine gravity_wave_drag_ridge_read_file + +end module gravity_wave_drag_ridge_read diff --git a/src/physics/utils/gravity_wave_drag_ridge_read.meta b/src/physics/utils/gravity_wave_drag_ridge_read.meta new file mode 100644 index 000000000..e3dd480b6 --- /dev/null +++ b/src/physics/utils/gravity_wave_drag_ridge_read.meta @@ -0,0 +1,82 @@ +[ccpp-table-properties] + name = gravity_wave_drag_ridge_read + type = module + +[ccpp-arg-table] + name = gravity_wave_drag_ridge_read + type = module +[ prdg ] + standard_name = number_of_ridges_in_ridge_gravity_wave_drag + units = count + type = integer + dimensions = () +[ rdg_gbxar ] + standard_name = grid_box_area_for_beta_ridge_gravity_wave_drag + units = km2 + type = real | kind = kind_phys + dimensions = (horizontal_dimension) +[ rdg_isovar ] + standard_name = isotropic_variance_for_beta_ridge_gravity_wave_drag + units = m + type = real | kind = kind_phys + dimensions = (horizontal_dimension) +[ rdg_isowgt ] + standard_name = isotropic_weight_for_beta_ridge_gravity_wave_drag + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_dimension) +[ rdg_hwdth ] + standard_name = ridge_half_width_for_beta_ridge_gravity_wave_drag + units = km + type = real | kind = kind_phys + dimensions = (horizontal_dimension, number_of_ridges_in_ridge_gravity_wave_drag) +[ rdg_clngt ] + standard_name = ridge_length_for_beta_ridge_gravity_wave_drag + units = km + type = real | kind = kind_phys + dimensions = (horizontal_dimension, number_of_ridges_in_ridge_gravity_wave_drag) +[ rdg_mxdis ] + standard_name = ridge_obstacle_height_for_beta_ridge_gravity_wave_drag + units = m + type = real | kind = kind_phys + dimensions = (horizontal_dimension, number_of_ridges_in_ridge_gravity_wave_drag) +[ rdg_anixy ] + standard_name = ridge_anisotropy_for_beta_ridge_gravity_wave_drag + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_dimension, number_of_ridges_in_ridge_gravity_wave_drag) +[ rdg_angll ] + standard_name = ridge_clockwise_angle_from_north_for_beta_ridge_gravity_wave_drag + units = degree + type = real | kind = kind_phys + dimensions = (horizontal_dimension, number_of_ridges_in_ridge_gravity_wave_drag) +[ rdg_gbxarg ] + standard_name = grid_box_area_for_gamma_ridge_gravity_wave_drag + units = km2 + type = real | kind = kind_phys + dimensions = (horizontal_dimension) +[ rdg_hwdthg ] + standard_name = ridge_half_width_for_gamma_ridge_gravity_wave_drag + units = km + type = real | kind = kind_phys + dimensions = (horizontal_dimension, number_of_ridges_in_ridge_gravity_wave_drag) +[ rdg_clngtg ] + standard_name = ridge_length_for_gamma_ridge_gravity_wave_drag + units = km + type = real | kind = kind_phys + dimensions = (horizontal_dimension, number_of_ridges_in_ridge_gravity_wave_drag) +[ rdg_mxdisg ] + standard_name = ridge_obstacle_height_for_gamma_ridge_gravity_wave_drag + units = m + type = real | kind = kind_phys + dimensions = (horizontal_dimension, number_of_ridges_in_ridge_gravity_wave_drag) +[ rdg_anixyg ] + standard_name = ridge_anisotropy_for_gamma_ridge_gravity_wave_drag + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_dimension, number_of_ridges_in_ridge_gravity_wave_drag) +[ rdg_angllg ] + standard_name = ridge_clockwise_angle_from_north_for_gamma_ridge_gravity_wave_drag + units = degree + type = real | kind = kind_phys + dimensions = (horizontal_dimension, number_of_ridges_in_ridge_gravity_wave_drag) diff --git a/test/unit/python/sample_files/ref_pres.meta b/test/unit/python/sample_files/ref_pres.meta index 4afd23c37..ba4140bd8 100644 --- a/test/unit/python/sample_files/ref_pres.meta +++ b/test/unit/python/sample_files/ref_pres.meta @@ -65,7 +65,7 @@ type = real | kind = kind_phys dimensions = () [ do_molec_diff ] - standard_name = flag_for_molecular_diffusion + standard_name = do_molecular_diffusion units = flag type = logical dimensions = () diff --git a/test/unit/python/sample_files/ref_pres_SourceMods.meta b/test/unit/python/sample_files/ref_pres_SourceMods.meta index 7a42dc996..141eb2fde 100644 --- a/test/unit/python/sample_files/ref_pres_SourceMods.meta +++ b/test/unit/python/sample_files/ref_pres_SourceMods.meta @@ -65,7 +65,7 @@ type = real | kind = kind_phys dimensions = () [ do_molec_diff ] - standard_name = flag_for_molecular_diffusion + standard_name = do_molecular_diffusion units = flag type = logical dimensions = ()