From 0b8f36b063372ab6b164ffa05653605e2f3b83e3 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 31 Jul 2022 13:35:27 -0400 Subject: [PATCH 1/2] (*)Set defaults in USE_REGRIDDING get_param calls Set default values in all get_param calls for USE_REGRIDDING. Previously, there had been 4 calls where this was missing, which led to the problems noted at https://github.com/mom-ocean/MOM6/issues/1576. This PR will allow that issue to be closed. Also used the default argument in a get_param call for INPUTDIR, although that case would not change any behavior because the value was set before the get_param call. A fail_in_missing argument was added to the FMS_cap call to get_param for GUST_2D_FILE, mirroring what is done for the solo_driver code, but cases where this was actually missing were very likely to have failed later anyway, but without an explicit error message. This PR could change unpredictable behavior in cases where USE_REGRIDDING is not explicitly set, but all answers are bitwise identical in the MOM6-examples test suite. --- .../FMS_cap/MOM_surface_forcing_gfdl.F90 | 10 +++++----- .../MOM_state_initialization.F90 | 18 ++++++++--------- src/user/dumbbell_initialization.F90 | 20 +++++++++---------- 3 files changed, 24 insertions(+), 24 deletions(-) diff --git a/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 b/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 index faa74a7fe0..f38c3fa7d2 100644 --- a/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 +++ b/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 @@ -1387,7 +1387,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, wind_stagger) call get_param(param_file, mdl, "FLUXCONST_SALT", CS%Flux_const_salt, & "The constant that relates the restoring surface salt fluxes to the relative "//& "surface anomalies (akin to a piston velocity). Note the non-MKS units.", & - fail_if_missing=.false.,default=unscaled_fluxconst, units="m day-1", scale=US%m_to_Z*US%T_to_s) + fail_if_missing=.false., default=unscaled_fluxconst, units="m day-1", scale=US%m_to_Z*US%T_to_s) ! Finish converting CS%Flux_const from m day-1 to [Z T-1 ~> m s-1]. CS%Flux_const = CS%Flux_const / 86400.0 CS%Flux_const_salt = CS%Flux_const_salt / 86400.0 @@ -1435,11 +1435,11 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, wind_stagger) call get_param(param_file, mdl, "FLUXCONST", CS%Flux_const, & "The constant that relates the restoring surface fluxes to the relative "//& "surface anomalies (akin to a piston velocity). Note the non-MKS units.", & - default=0.0, units="m day-1", scale=US%m_to_Z*US%T_to_s,unscaled=unscaled_fluxconst) + default=0.0, units="m day-1", scale=US%m_to_Z*US%T_to_s, unscaled=unscaled_fluxconst) call get_param(param_file, mdl, "FLUXCONST_TEMP", CS%Flux_const_temp, & "The constant that relates the restoring surface temperature fluxes to the relative "//& "surface anomalies (akin to a piston velocity). Note the non-MKS units.", & - fail_if_missing=.false.,default=unscaled_fluxconst, units="m day-1", scale=US%m_to_Z*US%T_to_s) + fail_if_missing=.false., default=unscaled_fluxconst, units="m day-1", scale=US%m_to_Z*US%T_to_s) ! Convert CS%Flux_const from m day-1 to m s-1. CS%Flux_const = CS%Flux_const / 86400.0 CS%Flux_const_temp = CS%Flux_const_temp / 86400.0 @@ -1523,7 +1523,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, wind_stagger) if (CS%read_gust_2d) then call get_param(param_file, mdl, "GUST_2D_FILE", gust_file, & "The file in which the wind gustiness is found in "//& - "variable gustiness.") + "variable gustiness.", fail_if_missing=.true.) call safe_alloc_ptr(CS%gust,isd,ied,jsd,jed) gust_file = trim(CS%inputdir) // trim(gust_file) @@ -1549,7 +1549,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, wind_stagger) if (CS%rigid_sea_ice) then call get_param(param_file, mdl, "G_EARTH", CS%g_Earth, & "The gravitational acceleration of the Earth.", & - units="m s-2", default = 9.80, scale=US%Z_to_m*US%m_s_to_L_T**2) + units="m s-2", default=9.80, scale=US%Z_to_m*US%m_s_to_L_T**2) call get_param(param_file, mdl, "SEA_ICE_MEAN_DENSITY", CS%density_sea_ice, & "A typical density of sea ice, used with the kinematic "//& "viscosity, when USE_RIGID_SEA_ICE is true.", & diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 257d25dad0..d21c13a3e5 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -1087,11 +1087,11 @@ subroutine depress_surface(h, G, GV, US, param_file, tv, just_read, z_top_shelf) call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") inputdir = slasher(inputdir) - call get_param(param_file, mdl, "SURFACE_HEIGHT_IC_FILE", eta_srf_file,& + call get_param(param_file, mdl, "SURFACE_HEIGHT_IC_FILE", eta_srf_file, & "The initial condition file for the surface height.", & fail_if_missing=.not.just_read, do_not_log=just_read) call get_param(param_file, mdl, "SURFACE_HEIGHT_IC_VAR", eta_srf_var, & - "The initial condition variable for the surface height.",& + "The initial condition variable for the surface height.", & default="SSH", do_not_log=just_read) filename = trim(inputdir)//trim(eta_srf_file) if (.not.just_read) & @@ -1263,7 +1263,7 @@ subroutine calc_sfc_displacement(PF, G, GV, US, mass_shelf, tv, h) call get_param(PF, mdl, "ICE_SHELF_INITIALIZATION_Z_TOLERANCE", tol, & "A initialization tolerance for the calculation of the static "// & - "ice shelf displacement (m) using initial temperature and salinity profile.",& + "ice shelf displacement (m) using initial temperature and salinity profile.", & default=0.001, units="m", scale=US%m_to_Z) max_iter = 1e3 call MOM_mesg("Started calculating initial interface position under ice shelf ") @@ -1949,13 +1949,13 @@ subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, depth_t "The name of the inverse damping rate variable in "//& "SPONGE_UV_DAMPING_FILE for the velocities.", default=Idamp_var) endif - call get_param(param_file, mdl, "USE_REGRIDDING", use_ALE, do_not_log=.true.) + call get_param(param_file, mdl, "USE_REGRIDDING", use_ALE, default=.false., do_not_log=.true.) !### NEW_SPONGES should be obsoleted properly, rather than merely deprecated, at which ! point only the else branch of the new_sponge_param block would be retained. call get_param(param_file, mdl, "NEW_SPONGES", new_sponge_param, & "Set True if using the newer sponging code which "//& - "performs on-the-fly regridding in lat-lon-time.",& + "performs on-the-fly regridding in lat-lon-time"//& "of sponge restoring data.", default=.false., do_not_log=.true.) if (new_sponge_param) then call get_param(param_file, mdl, "INTERPOLATE_SPONGE_TIME_SPACE", time_space_interp_sponge, & @@ -2230,7 +2230,7 @@ subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, h, u, v, p default=.false.) endif call get_param(param_file, mdl, "ODA_INCUPD_RESET_NCOUNT", reset_ncount, & - "If True, reinitialize number of updates already done, ncount.",& + "If True, reinitialize number of updates already done, ncount.", & default=.true.) if (.not.oda_inc .and. .not.reset_ncount) & call MOM_error(FATAL, " initialize_oda_incupd: restarting during update "// & @@ -2258,7 +2258,7 @@ subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, h, u, v, p "The name of the meridional vel. inc. variable in "//& "ODA_INCUPD_FILE.", default="v_inc") -! call get_param(param_file, mdl, "USE_REGRIDDING", use_ALE, do_not_log=.true.) +! call get_param(param_file, mdl, "USE_REGRIDDING", use_ALE, default=.false., do_not_log=.true.) ! Read in incremental update for tracers filename = trim(inputdir)//trim(inc_file) @@ -2486,7 +2486,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just if (.not.just_read) call callTree_enter(trim(mdl)//"(), MOM_state_initialization.F90") if (.not.just_read) call log_version(PF, mdl, version, "") - inputdir = "." ; call get_param(PF, mdl, "INPUTDIR", inputdir) + call get_param(PF, mdl, "INPUTDIR", inputdir, default=".") inputdir = slasher(inputdir) eos => tv%eqn_of_state @@ -2525,7 +2525,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just "is True.", default="PPM_IH4", do_not_log=just_read) call get_param(PF, mdl, "Z_INIT_REMAP_GENERAL", remap_general, & "If false, only initializes to z* coordinates. "//& - "If true, allows initialization directly to general coordinates.",& + "If true, allows initialization directly to general coordinates.", & default=.false., do_not_log=just_read) call get_param(PF, mdl, "Z_INIT_REMAP_FULL_COLUMN", remap_full_column, & "If false, only reconstructs profiles for valid data points. "//& diff --git a/src/user/dumbbell_initialization.F90 b/src/user/dumbbell_initialization.F90 index 570e638465..e4ce7e77f5 100644 --- a/src/user/dumbbell_initialization.F90 +++ b/src/user/dumbbell_initialization.F90 @@ -51,13 +51,13 @@ subroutine dumbbell_initialize_topography( D, G, param_file, max_depth ) logical :: dbrotate call get_param(param_file, mdl, "DUMBBELL_LEN",dblen, & - 'Lateral Length scale for dumbbell.',& + 'Lateral Length scale for dumbbell.', & units='km', default=600., do_not_log=.false.) call get_param(param_file, mdl, "DUMBBELL_FRACTION",dbfrac, & - 'Meridional fraction for narrow part of dumbbell.',& + 'Meridional fraction for narrow part of dumbbell.', & units='nondim', default=0.5, do_not_log=.false.) call get_param(param_file, mdl, "DUMBBELL_ROTATION", dbrotate, & - 'Logical for rotation of dumbbell domain.',& + 'Logical for rotation of dumbbell domain.', & units='nondim', default=.false., do_not_log=.false.) if (G%x_axis_units == 'm') then @@ -128,11 +128,11 @@ subroutine dumbbell_initialize_thickness ( h, depth_tot, G, GV, US, param_file, if (.not.just_read) call log_version(param_file, mdl, version, "") call get_param(param_file, mdl,"MIN_THICKNESS", min_thickness, & - 'Minimum thickness for layer',& + 'Minimum thickness for layer', & units='m', default=1.0e-3, scale=US%m_to_Z, do_not_log=just_read) call get_param(param_file, mdl,"REGRIDDING_COORDINATE_MODE", verticalCoordinate, & default=DEFAULT_COORDINATE_MODE, do_not_log=just_read) - call get_param(param_file, mdl, "USE_REGRIDDING", use_ALE, do_not_log = .true.) + call get_param(param_file, mdl, "USE_REGRIDDING", use_ALE, default=.false., do_not_log=.true.) if (.not. use_ALE) verticalCoordinate = "LAYER" ! WARNING: this routine specifies the interface heights so that the last layer @@ -149,7 +149,7 @@ subroutine dumbbell_initialize_thickness ( h, depth_tot, G, GV, US, param_file, select case ( coordinateMode(verticalCoordinate) ) case ( REGRIDDING_LAYER) ! Initial thicknesses for isopycnal coordinates call get_param(param_file, mdl, "DUMBBELL_ROTATION", dbrotate, & - 'Logical for rotation of dumbbell domain.',& + 'Logical for rotation of dumbbell domain.', & units='nondim', default=.false., do_not_log=just_read) do j=js,je do i=is,ie @@ -273,7 +273,7 @@ subroutine dumbbell_initialize_temperature_salinity ( T, S, h, G, GV, US, param_ T_surf = 20.0*US%degC_to_C ! layer mode - call get_param(param_file, mdl, "USE_REGRIDDING", use_ALE, do_not_log = .true.) + call get_param(param_file, mdl, "USE_REGRIDDING", use_ALE, default=.false., do_not_log=.true.) if (.not. use_ALE) call MOM_error(FATAL, "dumbbell_initialize_temperature_salinity: "//& "Please use 'fit' for 'TS_CONFIG' in the LAYER mode.") @@ -357,10 +357,10 @@ subroutine dumbbell_initialize_sponges(G, GV, US, tv, h_in, depth_tot, param_fil logical :: dbrotate ! If true, rotate the domain. call get_param(param_file, mdl,"DUMBBELL_LEN",dblen, & - 'Lateral Length scale for dumbbell ',& + 'Lateral Length scale for dumbbell ', & units='km', default=600., do_not_log=.true.) call get_param(param_file, mdl, "DUMBBELL_ROTATION", dbrotate, & - 'Logical for rotation of dumbbell domain.',& + 'Logical for rotation of dumbbell domain.', & units='nondim', default=.false., do_not_log=.true.) if (G%x_axis_units == 'm') then @@ -379,7 +379,7 @@ subroutine dumbbell_initialize_sponges(G, GV, US, tv, h_in, depth_tot, param_fil 'DUMBBELL salinity range (right-left)', & units='1e-3', default=2., scale=US%ppt_to_S, do_not_log=.true.) call get_param(param_file, mdl,"MIN_THICKNESS", min_thickness, & - 'Minimum thickness for layer',& + 'Minimum thickness for layer', & units='m', default=1.0e-3, scale=US%m_to_Z, do_not_log=.true.) ! no active sponges From acd8c7b4ee893d80abbf5fd578d51b6d450fda3e Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 2 Aug 2022 07:04:03 -0400 Subject: [PATCH 2/2] Pass fail_if_missing to 4 user get_param calls Added fail_if_missing or default arguments to 5 get_param calls in user code, with values set consistently with other calls for the same parameters. Also replaced 4 get_param calls with copies of equivalent fields from the grid type in one use initialization routine. All answers are bitwise identical. --- src/user/BFB_initialization.F90 | 18 +++++++++--------- src/user/DOME2d_initialization.F90 | 2 +- src/user/Kelvin_initialization.F90 | 2 +- src/user/dense_water_initialization.F90 | 6 ++++-- src/user/shelfwave_initialization.F90 | 2 +- 5 files changed, 16 insertions(+), 14 deletions(-) diff --git a/src/user/BFB_initialization.F90 b/src/user/BFB_initialization.F90 index 22d3156723..68a6b6530b 100644 --- a/src/user/BFB_initialization.F90 +++ b/src/user/BFB_initialization.F90 @@ -93,7 +93,11 @@ subroutine BFB_initialize_sponges_southonly(G, GV, US, use_temperature, tv, dept real :: Idamp(SZI_(G),SZJ_(G)) ! The sponge damping rate [T-1 ~> s-1] real :: H0(SZK_(GV)) ! Resting layer thicknesses in depth units [Z ~> m]. real :: min_depth ! The minimum ocean depth in depth units [Z ~> m]. - real :: slat, wlon, lenlat, lenlon, nlat + real :: slat ! The southern latitude of the domain [degrees_N] + real :: wlon ! The western longitude of the domain [degrees_E] + real :: lenlat ! The latitudinal length of the domain [degrees_N] + real :: lenlon ! The longitudinal length of the domain [degrees_E] + real :: nlat ! The northern latitude of the domain [degrees_N] real :: max_damping ! The maximum damping rate [T-1 ~> s-1] character(len=40) :: mdl = "BFB_initialize_sponges_southonly" ! This subroutine's name. integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz @@ -112,14 +116,10 @@ subroutine BFB_initialize_sponges_southonly(G, GV, US, use_temperature, tv, dept call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, & "The minimum depth of the ocean.", units="m", default=0.0, scale=US%m_to_Z) - call get_param(param_file, mdl, "SOUTHLAT", slat, & - "The southern latitude of the domain.", units="degrees") - call get_param(param_file, mdl, "LENLAT", lenlat, & - "The latitudinal length of the domain.", units="degrees") - call get_param(param_file, mdl, "WESTLON", wlon, & - "The western longitude of the domain.", units="degrees", default=0.0) - call get_param(param_file, mdl, "LENLON", lenlon, & - "The longitudinal length of the domain.", units="degrees") + slat = G%south_lat + lenlat = G%len_lat + wlon = G%west_lon + lenlon = G%len_lon nlat = slat + lenlat do k=1,nz ; H0(k) = -G%max_depth * real(k-1) / real(nz) ; enddo diff --git a/src/user/DOME2d_initialization.F90 b/src/user/DOME2d_initialization.F90 index 393347d1f2..d0ed88c128 100644 --- a/src/user/DOME2d_initialization.F90 +++ b/src/user/DOME2d_initialization.F90 @@ -411,7 +411,7 @@ subroutine DOME2d_initialize_sponges(G, GV, US, tv, depth_tot, param_file, use_A call get_param(param_file, mdl, "DOME2D_SHELF_DEPTH", dome2d_depth_bay, & default=0.2, do_not_log=.true.) call get_param(param_file, mdl, "S_REF", S_ref, default=35.0, scale=US%ppt_to_S) - call get_param(param_file, mdl, "T_REF", T_ref, scale=US%degC_to_C) + call get_param(param_file, mdl, "T_REF", T_ref, scale=US%degC_to_C, fail_if_missing=.false.) call get_param(param_file, mdl, "S_RANGE", S_range, default=2.0, scale=US%ppt_to_S) call get_param(param_file, mdl, "T_RANGE", T_range, default=0.0, scale=US%degC_to_C) diff --git a/src/user/Kelvin_initialization.F90 b/src/user/Kelvin_initialization.F90 index a65dc45e73..595736540e 100644 --- a/src/user/Kelvin_initialization.F90 +++ b/src/user/Kelvin_initialization.F90 @@ -75,7 +75,7 @@ function register_Kelvin_OBC(param_file, CS, US, OBC_Reg) default=0) call get_param(param_file, mdl, "F_0", CS%F_0, & default=0.0, units="s-1", scale=US%T_to_s, do_not_log=.true.) - call get_param(param_file, mdl, "TOPO_CONFIG", config, do_not_log=.true.) + call get_param(param_file, mdl, "TOPO_CONFIG", config, fail_if_missing=.true., do_not_log=.true.) if (trim(config) == "Kelvin") then call get_param(param_file, mdl, "ROTATED_COAST_OFFSET_1", CS%coast_offset1, & "The distance along the southern and northern boundaries "//& diff --git a/src/user/dense_water_initialization.F90 b/src/user/dense_water_initialization.F90 index 1c372bf1b7..fa44a78604 100644 --- a/src/user/dense_water_initialization.F90 +++ b/src/user/dense_water_initialization.F90 @@ -197,8 +197,10 @@ subroutine dense_water_initialize_sponges(G, GV, US, tv, depth_tot, param_file, call get_param(param_file, mdl, "DENSE_WATER_SILL_HEIGHT", sill_height, default=default_sill, do_not_log=.true.) call get_param(param_file, mdl, "S_REF", S_ref, default=35.0, scale=US%ppt_to_S, do_not_log=.true.) - call get_param(param_file, mdl, "S_RANGE", S_range, scale=US%ppt_to_S, do_not_log=.true.) - call get_param(param_file, mdl, "T_REF", T_ref, scale=US%degC_to_C, do_not_log=.true.) + call get_param(param_file, mdl, "S_RANGE", S_range, & + units='1e-3', default=2.0, scale=US%ppt_to_S, do_not_log=.true.) + call get_param(param_file, mdl, "T_REF", T_ref, & + units='degC', scale=US%degC_to_C, fail_if_missing=.true., do_not_log=.true.) ! no active sponges if (west_sponge_time_scale <= 0. .and. east_sponge_time_scale <= 0.) return diff --git a/src/user/shelfwave_initialization.F90 b/src/user/shelfwave_initialization.F90 index 3bb031bbb6..a9c1914356 100644 --- a/src/user/shelfwave_initialization.F90 +++ b/src/user/shelfwave_initialization.F90 @@ -66,7 +66,7 @@ function register_shelfwave_OBC(param_file, CS, US, OBC_Reg) call get_param(param_file, mdl, "F_0", CS%f0, & default=0.0, units="s-1", scale=US%T_to_s, do_not_log=.true.) call get_param(param_file, mdl, "LENLAT", len_lat, & - do_not_log=.true.) + do_not_log=.true., fail_if_missing=.true.) call get_param(param_file, mdl,"SHELFWAVE_X_WAVELENGTH",CS%Lx, & "Length scale of shelfwave in x-direction.",& units="Same as x,y", default=100.)