Skip to content

+Set input file variable names via runtime params #231

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
marshallward merged 2 commits into from
Nov 16, 2022
Merged

+Set input file variable names via runtime params #231

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
47a6ad1
+Set input file variable names via runtime params
Hallberg-NOAA Nov 7, 2022
bf9e191
+Set more input file variable names at runtime
Hallberg-NOAA Nov 14, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
41 changes: 22 additions & 19 deletions src/ice_shelf/MOM_ice_shelf.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -139,7 +139,7 @@ module MOM_ice_shelf
real :: time_step !< this is the shortest timestep that the ice shelf sees [T ~> s], and
!! is equal to the forcing timestep (it is passed in when the shelf
!! is initialized - so need to reorganize MOM driver.
!! it will be the prognistic timestep ... maybe.
!! it will be the prognostic timestep ... maybe.

logical :: solo_ice_sheet !< whether the ice model is running without being
!! coupled to the ocean
Expand Down Expand Up @@ -288,13 +288,13 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS)
real :: n_star_term ! A term in the expression for nstar [T3 Z-2 ~> s3 m-2]
real :: absf ! The absolute value of the Coriolis parameter [T-1 ~> s-1]
real :: dIns_dwB !< The partial derivative of I_n_star with wB_flux, in [T3 Z-2 ~> s3 m-2]
real :: dT_ustar ! The difference between the the freezing point and the ocean boundary layer
real :: dT_ustar ! The difference between the freezing point and the ocean boundary layer
! temperature times the friction velocity [C Z T-1 ~> degC m s-1]
real :: dS_ustar ! The difference between the salinity at the ice-ocean interface and the ocean
! boundary layer salinity times the friction velocity [S Z T-1 ~> ppt m s-1]
real :: ustar_h ! The friction velocity in the water below the ice shelf [Z T-1 ~> m s-1]
real :: Gam_turb ! [nondim]
real :: Gam_mol_t, Gam_mol_s ! Relative coefficients of molecular diffusivites [nondim]
real :: Gam_mol_t, Gam_mol_s ! Relative coefficients of molecular diffusivities [nondim]
real :: RhoCp ! A typical ocean density times the heat capacity of water [Q R C-1 ~> J m-3 degC-1]
real :: ln_neut
real :: mass_exch ! A mass exchange rate [R Z T-1 ~> kg m-2 s-1]
Expand All @@ -312,7 +312,7 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS)
real :: Irho0 ! The inverse of the mean density times a unit conversion factor [R-1 L Z-1 ~> m3 kg-1]
logical :: Sb_min_set, Sb_max_set
logical :: update_ice_vel ! If true, it is time to update the ice shelf velocities.
logical :: coupled_GL ! If true, the grouding line position is determined based on
logical :: coupled_GL ! If true, the grounding line position is determined based on
! coupled ice-ocean dynamics.

real, parameter :: c2_3 = 2.0/3.0
Expand Down Expand Up @@ -524,7 +524,7 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS)
wB_flux = dB_dS * (dS_ustar * I_Gam_S) + dB_dT * wT_flux

if (wB_flux < 0.0) then
! The buoyancy flux is stabilizing and will reduce the tubulent
! The buoyancy flux is stabilizing and will reduce the turbulent
! fluxes, and iteration is required.
n_star_term = (ZETA_N/RC) * (hBL_neut * VK) / (ustar_h)**3
do it3 = 1,30
Expand Down Expand Up @@ -572,9 +572,9 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS)
exch_vel_s(i,j) = ustar_h * I_Gam_S

! Calculate the heat flux inside the ice shelf.
! Vertical adv/diff as in H+J 1999, eqns (26) & approx from (31).
! Vertical adv/diff as in H+J 1999, equations (26) & approx from (31).
! Q_ice = density_ice * CS%Cp_ice * K_ice * dT/dz (at interface)
! vertical adv/diff as in H+J 1999, eqs (31) & (26)...
! vertical adv/diff as in H+J 1999, equations (31) & (26)...
! dT/dz ~= min( (lprec/(density_ice*K_ice))*(CS%Temp_Ice-T_freeze) , 0.0 )
! If this approximation is not made, iterations are required... See H+J Fig 3.

Expand Down Expand Up @@ -1012,7 +1012,7 @@ subroutine add_shelf_flux(G, US, CS, sfc_state, fluxes)
real :: balancing_area !< total area where the balancing flux is applied [m2]
type(time_type) :: dTime !< The time step as a time_type
type(time_type) :: Time0 !< The previous time (Time-dt)
real, dimension(SZDI_(G),SZDJ_(G)) :: bal_frac !< Fraction of the cel1 where the mass flux
real, dimension(SZDI_(G),SZDJ_(G)) :: bal_frac !< Fraction of the cell where the mass flux
!! balancing the net melt flux occurs, 0 to 1 [nondim]
real, dimension(SZDI_(G),SZDJ_(G)) :: last_mass_shelf !< Ice shelf mass
!! at at previous time (Time-dt) [R Z ~> kg m-2]
Expand Down Expand Up @@ -1235,13 +1235,14 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in,
logical :: new_sim, save_IC
!This include declares and sets the variable "version".
# include "version_variable.h"
character(len=200) :: IC_file, inputdir
character(len=200) :: IC_file, inputdir ! Input file names or paths
character(len=40) :: mdl = "MOM_ice_shelf" ! This module's name.
integer :: i, j, is, ie, js, je, isd, ied, jsd, jed, Isdq, Iedq, Jsdq, Jedq
integer :: wd_halos(2)
logical :: read_TideAmp, shelf_mass_is_dynamic, debug
logical :: global_indexing
character(len=240) :: Tideamp_file
character(len=240) :: Tideamp_file ! Input file names
character(len=80) :: tideamp_var ! Input file variable names
real :: utide ! A tidal velocity [L T-1 ~> m s-1]
real :: col_thick_melt_thresh ! An ocean column thickness below which iceshelf melting
! does not occur [Z ~> m]
Expand Down Expand Up @@ -1397,7 +1398,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in,
"consistency to calculate the fluxes at the ice-ocean "//&
"interface.", default=.true.)
call get_param(param_file, mdl, "SHELF_INSULATOR", CS%insulator, &
"If true, the ice shelf is a perfect insulatior "//&
"If true, the ice shelf is a perfect insulator "//&
"(no conduction).", default=.false.)
call get_param(param_file, mdl, "MELTING_CUTOFF_DEPTH", CS%cutoff_depth, &
"Depth above which the melt is set to zero (it must be >= 0) "//&
Expand Down Expand Up @@ -1491,7 +1492,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in,
"The viscosity of the ice.", &
units="m2 s-1", default=1.0e10, scale=US%Z_to_L**2*US%m_to_L**2*US%T_to_s)
call get_param(param_file, mdl, "KV_MOLECULAR", CS%kv_molec, &
"The molecular kinimatic viscosity of sea water at the freezing temperature.", &
"The molecular kinematic viscosity of sea water at the freezing temperature.", &
units="m2 s-1", default=1.95e-6, scale=US%m2_s_to_Z2_T)
call get_param(param_file, mdl, "ICE_SHELF_SALINITY", CS%Salin_ice, &
"The salinity of the ice inside the ice shelf.", &
Expand Down Expand Up @@ -1537,19 +1538,21 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in,

if (read_TIDEAMP) then
call get_param(param_file, mdl, "TIDEAMP_FILE", TideAmp_file, &
"The path to the file containing the spatially varying "//&
"tidal amplitudes.", &
"The path to the file containing the spatially varying tidal amplitudes.", &
default="tideamp.nc")
call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
call get_param(param_file, mdl, "TIDEAMP_VARNAME", tideamp_var, &
"The name of the tidal amplitude variable in the input file.", &
default="tideamp")
call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
inputdir = slasher(inputdir)
TideAmp_file = trim(inputdir) // trim(TideAmp_file)
if (CS%rotate_index) then
allocate(tmp2d(CS%Grid_in%isd:CS%Grid_in%ied,CS%Grid_in%jsd:CS%Grid_in%jed), source=0.0)
call MOM_read_data(TideAmp_file, 'tideamp', tmp2d, CS%Grid_in%domain, timelevel=1, scale=US%m_s_to_L_T)
call MOM_read_data(TideAmp_file, tideamp_var, tmp2d, CS%Grid_in%domain, timelevel=1, scale=US%m_s_to_L_T)
call rotate_array(tmp2d, CS%turns, CS%utide)
deallocate(tmp2d)
else
call MOM_read_data(TideAmp_file, 'tideamp', CS%utide, CS%Grid%domain, timelevel=1, scale=US%m_s_to_L_T)
call MOM_read_data(TideAmp_file, tideamp_var, CS%utide, CS%Grid%domain, timelevel=1, scale=US%m_s_to_L_T)
endif
else
call get_param(param_file, mdl, "UTIDE", utide, &
Expand Down Expand Up @@ -2192,7 +2195,7 @@ subroutine solo_step_ice_shelf(CS, time_interval, nsteps, Time, min_time_step_in
real :: min_time_step ! The minimal required timestep that would indicate a fatal problem [T ~> s]
character(len=240) :: mesg
logical :: update_ice_vel ! If true, it is time to update the ice shelf velocities.
logical :: coupled_GL ! If true the grouding line position is determined based on
logical :: coupled_GL ! If true the grounding line position is determined based on
! coupled ice-ocean dynamics.
integer :: is, iec, js, jec

Expand Down Expand Up @@ -2261,7 +2264,7 @@ end subroutine solo_step_ice_shelf
!! update_shelf_mass - updates ice shelf mass via netCDF file
!! USER_update_shelf_mass (TODO).
!! solo_step_ice_shelf - called only in ice-only mode.
!! shelf_calc_flux - after melt rate & fluxes are calculated, ice dynamics are done. currently mass_shelf is
!! shelf_calc_flux - after melt rate & fluxes are calculated, ice dynamics are done. Currently mass_shelf is
!! updated immediately after ice_shelf_advect in fully dynamic mode.
!!
!! NOTES: be aware that hmask(:,:) has a number of functions; it is used for front advancement,
Expand Down
20 changes: 14 additions & 6 deletions src/initialization/MOM_coord_initialization.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -261,6 +261,8 @@ subroutine set_coord_from_TS_profile(Rlay, g_prime, GV, US, param_file, eqn_of_s
integer :: k, nz
character(len=40) :: mdl = "set_coord_from_TS_profile" ! This subroutine's name.
character(len=200) :: filename, coord_file, inputdir ! Strings for file/path
character(len=64) :: temp_var, salt_var ! Temperature and salinity names in files

nz = GV%ke

call callTree_enter(trim(mdl)//"(), MOM_coord_initialization.F90")
Expand All @@ -269,15 +271,21 @@ subroutine set_coord_from_TS_profile(Rlay, g_prime, GV, US, param_file, eqn_of_s
"The reduced gravity at the free surface.", units="m s-2", &
default=GV%g_Earth*US%L_T_to_m_s**2*US%m_to_Z, scale=US%m_s_to_L_T**2*US%Z_to_m)
call get_param(param_file, mdl, "COORD_FILE", coord_file, &
"The file from which the coordinate temperatures and "//&
"salinities are read.", fail_if_missing=.true.)
"The file from which the coordinate temperatures and salinities are read.", &
fail_if_missing=.true.)
call get_param(param_file, mdl, "TEMP_COORD_VAR", temp_var, &
"The coordinate reference profile variable name for potential temperature.", &
default="PTEMP")
call get_param(param_file, mdl, "SALT_COORD_VAR", salt_var, &
"The coordinate reference profile variable name for salinity.", &
default="SALT")

call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
filename = trim(slasher(inputdir))//trim(coord_file)
call log_param(param_file, mdl, "INPUTDIR/COORD_FILE", filename)

call MOM_read_data(filename, "PTEMP", T0(:), scale=US%degC_to_C)
call MOM_read_data(filename, "SALT", S0(:), scale=US%ppt_to_S)
call MOM_read_data(filename, temp_var, T0(:), scale=US%degC_to_C)
call MOM_read_data(filename, salt_var, S0(:), scale=US%ppt_to_S)

if (.not.file_exists(filename)) call MOM_error(FATAL, &
" set_coord_from_TS_profile: Unable to open " //trim(filename))
Expand Down Expand Up @@ -357,7 +365,7 @@ subroutine set_coord_from_TS_range(Rlay, g_prime, GV, US, param_file, eqn_of_sta

k_light = GV%nk_rho_varies + 1

! Set T0(k) to range from T_LIGHT to T_DENSE, and simliarly for S0(k).
! Set T0(k) to range from T_LIGHT to T_DENSE, and similarly for S0(k).
T0(k_light) = T_Light ; S0(k_light) = S_Light
a1 = 2.0 * res_rat / (1.0 + res_rat)
do k=k_light+1,nz
Expand Down Expand Up @@ -458,7 +466,7 @@ subroutine set_coord_linear(Rlay, g_prime, GV, US, param_file)
"The reduced gravity at the free surface.", units="m s-2", &
default=GV%g_Earth*US%L_T_to_m_s**2*US%m_to_Z, scale=US%m_s_to_L_T**2*US%Z_to_m)

! This following sets the target layer densities such that a the
! This following sets the target layer densities such that the
! surface interface has density Rlay_ref and the bottom
! is Rlay_range larger
do k=1,nz
Expand Down
Loading