From dd1f3f52f5fef0e78ae7216e64d3489f28123cc9 Mon Sep 17 00:00:00 2001 From: claireyung <61528379+claireyung@users.noreply.github.com> Date: Wed, 15 Jan 2025 09:06:36 +1100 Subject: [PATCH 1/4] Add ice shelf pressure initialisation bug fix (#800) * Add ice shelf pressure initialisation bug fix This commit adds a new parameter FRAC_DP_AT_POS_NEGATIVE_P_BUGFIX and associated bug fix. This bug occurs when TRIM_IC_FOR_P_SURF pressure initialisation is used for ice shelf, whilst in Boussinesq mode and a nonlinear equation of state. The subroutine trim_for_ice calls cut_off_column_top. This routine in Boussinesq mode calls find_depth_of_pressure_in_cell in MOM_density_integrals.F90. This subsequently calls the function frac_dp_at_pos which calls the density equation of state based on given salinity, temperature and depth, which is incorrectly converted into pressure by z position (which is negative) multiplied by g and rho0. The bug results in incorrect densities from the equation of state and therefore an imperfect initialisation of layer thicknesses and sea surface height due to the displacement of water by the ice. The bug is fixed by multiplying the pressure by -1. This commit introduces parameter FRAC_DP_AT_POS_NEGATIVE_P_BUGFIX with default value False to preserve answers. If true, the ice shelf initialisation is accurate to a high degree with a nonlinear equation of state. Answers should not change, except for the added parameter in MOM_parameter_doc. The same functions are called by a unit test in MOM_state_initialization; in this commit I set the MOM_state_initialization unit test to use the existing bug (with a false flag). * Resolve indenting and white space inconsitencies with MOM6 style for ice shelf pressure initialisation bug fix FRAC_DP_AT_POS_NEGATIVE_P_BUGFIX --- src/core/MOM_density_integrals.F90 | 16 +++++++++----- .../MOM_state_initialization.F90 | 21 +++++++++++++------ 2 files changed, 26 insertions(+), 11 deletions(-) diff --git a/src/core/MOM_density_integrals.F90 b/src/core/MOM_density_integrals.F90 index 90994dd073..151bc4ef3c 100644 --- a/src/core/MOM_density_integrals.F90 +++ b/src/core/MOM_density_integrals.F90 @@ -2001,7 +2001,7 @@ end subroutine diagnose_mass_weight_p !> Find the depth at which the reconstructed pressure matches P_tgt subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_tgt, & - rho_ref, G_e, EOS, US, P_b, z_out, z_tol) + rho_ref, G_e, EOS, US, P_b, z_out, z_tol, frac_dp_bugfix) real, intent(in) :: T_t !< Potential temperature at the cell top [C ~> degC] real, intent(in) :: T_b !< Potential temperature at the cell bottom [C ~> degC] real, intent(in) :: S_t !< Salinity at the cell top [S ~> ppt] @@ -2020,6 +2020,7 @@ subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_t real, intent(out) :: P_b !< Pressure at the bottom of the cell [R L2 T-2 ~> Pa] real, intent(out) :: z_out !< Absolute depth at which anomalous pressure = p_tgt [Z ~> m] real, intent(in) :: z_tol !< The tolerance in finding z_out [Z ~> m] + logical, intent(in) :: frac_dp_bugfix !< If true, use bugfix in frac_dp_at_pos ! Local variables real :: dp ! Pressure thickness of the layer [R L2 T-2 ~> Pa] @@ -2032,7 +2033,7 @@ subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_t GxRho = G_e * rho_ref ! Anomalous pressure difference across whole cell - dp = frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, 1.0, EOS) + dp = frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, 1.0, EOS, frac_dp_bugfix) P_b = P_t + dp ! Anomalous pressure at bottom of cell @@ -2063,7 +2064,7 @@ subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_t call MOM_error(FATAL, 'find_depth_of_pressure_in_cell completes too many iterations: '//msg) endif z_out = z_t + ( z_b - z_t ) * F_guess - Pa = frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, F_guess, EOS) - ( P_tgt - P_t ) + Pa = frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, F_guess, EOS, frac_dp_bugfix) - ( P_tgt - P_t ) if (Pa Returns change in anomalous pressure change from top to non-dimensional !! position pos between z_t and z_b [R L2 T-2 ~> Pa] -real function frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, pos, EOS) +real function frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, pos, EOS, frac_dp_bugfix) real, intent(in) :: T_t !< Potential temperature at the cell top [C ~> degC] real, intent(in) :: T_b !< Potential temperature at the cell bottom [C ~> degC] real, intent(in) :: S_t !< Salinity at the cell top [S ~> ppt] @@ -2131,6 +2132,7 @@ real function frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, pos, EO real, intent(in) :: G_e !< The Earth's gravitational acceleration [L2 Z-1 T-2 ~> m s-2] real, intent(in) :: pos !< The fractional vertical position, 0 to 1 [nondim] type(EOS_type), intent(in) :: EOS !< Equation of state structure + logical, intent(in) :: frac_dp_bugfix !< If true, use bugfix in frac_dp_at_pos ! Local variables real, parameter :: C1_90 = 1.0/90.0 ! A rational constant [nondim] @@ -2150,7 +2152,11 @@ real function frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, pos, EO ! Salinity and temperature points are linearly interpolated S5(n) = top_weight * S_t + bottom_weight * S_b T5(n) = top_weight * T_t + bottom_weight * T_b - p5(n) = ( top_weight * z_t + bottom_weight * z_b ) * ( G_e * rho_ref ) + if (frac_dp_bugfix) then + p5(n) = (-1) * ( top_weight * z_t + bottom_weight * z_b ) * ( G_e * rho_ref ) + else + p5(n) = ( top_weight * z_t + bottom_weight * z_b ) * ( G_e * rho_ref ) + endif !bugfix enddo call calculate_density(T5, S5, p5, rho5, EOS) rho5(:) = rho5(:) !- rho_ref ! Work with anomalies relative to rho_ref diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index d908ec23a0..7909bddb80 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -1205,6 +1205,7 @@ subroutine trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read) ! answers from 2018, while higher values use more robust ! forms of the same remapping expressions. logical :: use_remapping ! If true, remap the initial conditions. + logical :: use_frac_dp_bugfix ! If true, use bugfix. Otherwise, pressure input to EOS is negative. type(remapping_CS), pointer :: remap_CS => NULL() call get_param(PF, mdl, "SURFACE_PRESSURE_FILE", p_surf_file, & @@ -1227,7 +1228,10 @@ subroutine trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read) "The tolerance with which to find the depth matching the specified "//& "surface pressure with TRIM_IC_FOR_P_SURF.", & units="m", default=1.0e-5, scale=US%m_to_Z, do_not_log=just_read) - + call get_param(PF, mdl, "FRAC_DP_AT_POS_NEGATIVE_P_BUGFIX", use_frac_dp_bugfix, & + "If true, use bugfix in ice shelf TRIM_IC initialization. "//& + "Otherwise, pressure input to density EOS is negative.", & + default=.false., do_not_log=just_read) call get_param(PF, mdl, "TRIMMING_USES_REMAPPING", use_remapping, & 'When trimming the column, also remap T and S.', & default=.false., do_not_log=just_read) @@ -1277,7 +1281,8 @@ subroutine trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read) do j=G%jsc,G%jec ; do i=G%isc,G%iec call cut_off_column_top(GV%ke, tv, GV, US, GV%g_Earth, G%bathyT(i,j)+G%Z_ref, min_thickness, & tv%T(i,j,:), T_t(i,j,:), T_b(i,j,:), tv%S(i,j,:), S_t(i,j,:), S_b(i,j,:), & - p_surf(i,j), h(i,j,:), remap_CS, z_tol=z_tolerance) + p_surf(i,j), h(i,j,:), remap_CS, z_tol=z_tolerance, & + frac_dp_bugfix=use_frac_dp_bugfix) enddo ; enddo end subroutine trim_for_ice @@ -1368,7 +1373,7 @@ end subroutine calc_sfc_displacement !> Adjust the layer thicknesses by removing the top of the water column above the !! depth where the hydrostatic pressure matches p_surf subroutine cut_off_column_top(nk, tv, GV, US, G_earth, depth, min_thickness, T, T_t, T_b, & - S, S_t, S_b, p_surf, h, remap_CS, z_tol) + S, S_t, S_b, p_surf, h, remap_CS, z_tol, frac_dp_bugfix) integer, intent(in) :: nk !< Number of layers type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. @@ -1388,6 +1393,7 @@ subroutine cut_off_column_top(nk, tv, GV, US, G_earth, depth, min_thickness, T, !! if associated real, intent(in) :: z_tol !< The tolerance with which to find the depth !! matching the specified pressure [Z ~> m]. + logical, intent(in) :: frac_dp_bugfix !< If true, use bugfix in frac_dp_at_pos ! Local variables real, dimension(nk+1) :: e ! Top and bottom edge positions for reconstructions [Z ~> m] @@ -1416,7 +1422,8 @@ subroutine cut_off_column_top(nk, tv, GV, US, G_earth, depth, min_thickness, T, do k=1,nk call find_depth_of_pressure_in_cell(T_t(k), T_b(k), S_t(k), S_b(k), e(K), e(K+1), & P_t, p_surf, GV%Rho0, G_earth, tv%eqn_of_state, & - US, P_b, z_out, z_tol=z_tol) + US, P_b, z_out, z_tol=z_tol, & + frac_dp_bugfix=frac_dp_bugfix) if (z_out>=e(K)) then ! Imposed pressure was less that pressure at top of cell exit @@ -3139,7 +3146,8 @@ subroutine MOM_state_init_tests(G, GV, US, tv) P_t = 0. do k = 1, nk call find_depth_of_pressure_in_cell(T_t(k), T_b(k), S_t(k), S_b(k), e(K), e(K+1), P_t, 0.5*P_tot, & - GV%Rho0, GV%g_Earth, tv%eqn_of_state, US, P_b, z_out, z_tol=z_tol) + GV%Rho0, GV%g_Earth, tv%eqn_of_state, US, P_b, z_out, z_tol=z_tol, & + frac_dp_bugfix=.false.) write(0,*) k, US%RL2_T2_to_Pa*P_t, US%RL2_T2_to_Pa*P_b, 0.5*US%RL2_T2_to_Pa*P_tot, & US%Z_to_m*e(K), US%Z_to_m*e(K+1), US%Z_to_m*z_out P_t = P_b @@ -3158,7 +3166,8 @@ subroutine MOM_state_init_tests(G, GV, US, tv) ! h_neglect=GV%H_subroundoff, h_neglect_edge=GV%H_subroundoff) ! endif call cut_off_column_top(nk, tv, GV, US, GV%g_Earth, -e(nk+1), GV%Angstrom_H, & - T, T_t, T_b, S, S_t, S_b, 0.5*P_tot, h, remap_CS, z_tol=z_tol) + T, T_t, T_b, S, S_t, S_b, 0.5*P_tot, h, remap_CS, z_tol=z_tol, & + frac_dp_bugfix=.false.) write(0,*) GV%H_to_m*h(:) if (associated(remap_CS)) deallocate(remap_CS) From d8da512c48a4445114d4b371144f625feace89cd Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 5 Jan 2025 06:22:40 -0500 Subject: [PATCH 2/4] +Add the optional argument old_name to get_param Added the new optional argument old_name to the 8 get_param() routines. This new capability allows for an archaic parameter name to be specified and for appropriate warnings encouraging the user to migrate to using the new name while still setting the parameter as intended, or error messages in the case of inconsistent setting via the archaic name and the correct name. The logging inside of the MOM_parameter_doc files only uses the correct parameter name. Also added the new optional argument set to the 8 read_param routines, to indicate whether a parameter has been found and successfully set. The new set argument is now being used in read_param() calls in obsolete_int(), obsolete_real(), obsolete_char() and obsolete_logical(). Obsolete_logical() in particular was substantially simplified by the use of this new argument, and is now only about half as long as it was. The read_param() set argument is also used in all of the get_param() routines when they are given an old_name argument. The new old_name argument to get_param() is not yet being used in the version of the MOM6 code that is being checked in, but it has been tested extensively by adding or modifying get_param calls in a variant of the initialization code, and it will be used in an updated version of github.com/NOAA-GFDL/MOM6/pull/725 to gracefully handle the deprecation of 4 parameter names. All answers are bitwise identical, but there are new optional arguments to two widely used interfaces. --- src/diagnostics/MOM_obsolete_params.F90 | 37 ++- src/framework/MOM_file_parser.F90 | 288 +++++++++++++++++++++--- 2 files changed, 274 insertions(+), 51 deletions(-) diff --git a/src/diagnostics/MOM_obsolete_params.F90 b/src/diagnostics/MOM_obsolete_params.F90 index a590ae3893..a0401d7199 100644 --- a/src/diagnostics/MOM_obsolete_params.F90 +++ b/src/diagnostics/MOM_obsolete_params.F90 @@ -168,29 +168,15 @@ subroutine obsolete_logical(param_file, varname, warning_val, hint) character(len=*), optional, intent(in) :: hint !< A hint to the user about what to do. ! Local variables logical :: test_logic, fatal_err + logical :: var_is_set ! True if this value was read by read_param. character(len=128) :: hint_msg - test_logic = .false. ; call read_param(param_file, varname, test_logic) + test_logic = .false. ; call read_param(param_file, varname, test_logic, set=var_is_set) fatal_err = .true. - if (present(warning_val)) fatal_err = (warning_val .neqv. .true.) + if (var_is_set .and. present(warning_val)) fatal_err = (warning_val .neqv. test_logic) hint_msg = " " ; if (present(hint)) hint_msg = hint - if (test_logic) then - if (fatal_err) then - call MOM_ERROR(FATAL, "MOM_obsolete_params: "//trim(varname)// & - " is an obsolete run-time flag, and should not be used. "// & - trim(hint_msg)) - else - call MOM_ERROR(WARNING, "MOM_obsolete_params: "//trim(varname)// & - " is an obsolete run-time flag. "//trim(hint_msg)) - endif - endif - - test_logic = .true. ; call read_param(param_file, varname, test_logic) - fatal_err = .true. - if (present(warning_val)) fatal_err = (warning_val .neqv. .false.) - - if (.not.test_logic) then + if (var_is_set) then if (fatal_err) then call MOM_ERROR(FATAL, "MOM_obsolete_params: "//trim(varname)// & " is an obsolete run-time flag, and should not be used. "// & @@ -211,12 +197,13 @@ subroutine obsolete_char(param_file, varname, warning_val, hint) character(len=*), optional, intent(in) :: hint !< A hint to the user about what to do. ! Local variables character(len=200) :: test_string, hint_msg + logical :: var_is_set ! True if this value was read by read_param. logical :: only_warn - test_string = ''; call read_param(param_file, varname, test_string) + test_string = ''; call read_param(param_file, varname, test_string, set=var_is_set) hint_msg = " " ; if (present(hint)) hint_msg = hint - if (len_trim(test_string) > 0) then + if (var_is_set) then only_warn = .false. if (present(warning_val)) then ! Check if test_string and warning_val are the same. if (len_trim(warning_val) == len_trim(test_string)) then @@ -246,15 +233,16 @@ subroutine obsolete_real(param_file, varname, warning_val, hint, only_warn) ! Local variables real :: test_val, warn_val + logical :: var_is_set ! True if this value was read by read_param. logical :: issue_warning character(len=128) :: hint_msg - test_val = -9e35; call read_param(param_file, varname, test_val) + test_val = -9e35; call read_param(param_file, varname, test_val, set=var_is_set) warn_val = -9e35; if (present(warning_val)) warn_val = warning_val hint_msg = " " ; if (present(hint)) hint_msg = hint issue_warning = .false. ; if (present(only_warn)) issue_warning = only_warn - if (test_val /= -9e35) then + if (var_is_set) then if ((test_val == warn_val) .or. issue_warning) then call MOM_ERROR(WARNING, "MOM_obsolete_params: "//trim(varname)// & " is an obsolete run-time flag. "//trim(hint_msg)) @@ -273,14 +261,15 @@ subroutine obsolete_int(param_file, varname, warning_val, hint) integer, optional, intent(in) :: warning_val !< An allowed value that causes a warning instead of an error. character(len=*), optional, intent(in) :: hint !< A hint to the user about what to do. ! Local variables + logical :: var_is_set ! True if this value was read by read_param. integer :: test_val, warn_val character(len=128) :: hint_msg - test_val = -123456788; call read_param(param_file, varname, test_val) + test_val = -123456788; call read_param(param_file, varname, test_val, set=var_is_set) warn_val = -123456788; if (present(warning_val)) warn_val = warning_val hint_msg = " " ; if (present(hint)) hint_msg = hint - if (test_val /= -123456788) then + if (var_is_set) then if (test_val == warn_val) then call MOM_ERROR(WARNING, "MOM_obsolete_params: "//trim(varname)// & " is an obsolete run-time flag. "//trim(hint_msg)) diff --git a/src/framework/MOM_file_parser.F90 b/src/framework/MOM_file_parser.F90 index 7d3337ea24..291d44492d 100644 --- a/src/framework/MOM_file_parser.F90 +++ b/src/framework/MOM_file_parser.F90 @@ -8,7 +8,7 @@ module MOM_file_parser use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, assert use MOM_error_handler, only : is_root_pe, stdlog, stdout use MOM_time_manager, only : get_time, time_type, get_ticks_per_second -use MOM_time_manager, only : set_date, get_date, real_to_time, operator(-), set_time +use MOM_time_manager, only : set_date, get_date, real_to_time, operator(-), operator(==), set_time use MOM_document, only : doc_param, doc_module, doc_init, doc_end, doc_type use MOM_document, only : doc_openBlock, doc_closeBlock use MOM_string_functions, only : left_int, left_ints, slasher @@ -618,7 +618,7 @@ function simplifyWhiteSpace(string) end function simplifyWhiteSpace !> This subroutine reads the value of an integer model parameter from a parameter file. -subroutine read_param_int(CS, varname, value, fail_if_missing) +subroutine read_param_int(CS, varname, value, fail_if_missing, set) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: varname !< The case-sensitive name of the parameter to read @@ -626,6 +626,8 @@ subroutine read_param_int(CS, varname, value, fail_if_missing) !! read from the parameter file logical, optional, intent(in) :: fail_if_missing !< If present and true, a fatal error occurs !! if this variable is not found in the parameter file + logical, optional, intent(out) :: set !< If present, this indicates whether this parameter + !! has been found and successfully set in the input files. ! Local variables character(len=CS%max_line_len) :: value_string(1) logical :: found, defined @@ -633,6 +635,7 @@ subroutine read_param_int(CS, varname, value, fail_if_missing) call get_variable_line(CS, varname, found, defined, value_string) if (found .and. defined .and. (LEN_TRIM(value_string(1)) > 0)) then read(value_string(1),*,err = 1001) value + if (present(set)) set = .true. else if (present(fail_if_missing)) then ; if (fail_if_missing) then if (.not.found) then @@ -643,6 +646,7 @@ subroutine read_param_int(CS, varname, value, fail_if_missing) ' found but not set in input files.') endif endif ; endif + if (present(set)) set = .false. endif return 1001 call MOM_error(FATAL,'read_param_int: read error for integer variable '//trim(varname)// & @@ -650,7 +654,7 @@ subroutine read_param_int(CS, varname, value, fail_if_missing) end subroutine read_param_int !> This subroutine reads the values of an array of integer model parameters from a parameter file. -subroutine read_param_int_array(CS, varname, value, fail_if_missing) +subroutine read_param_int_array(CS, varname, value, fail_if_missing, set) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: varname !< The case-sensitive name of the parameter to read @@ -658,12 +662,15 @@ subroutine read_param_int_array(CS, varname, value, fail_if_missing) !! read from the parameter file logical, optional, intent(in) :: fail_if_missing !< If present and true, a fatal error occurs !! if this variable is not found in the parameter file + logical, optional, intent(out) :: set !< If present, this indicates whether this parameter + !! has been found and successfully set in the input files. ! Local variables character(len=CS%max_line_len) :: value_string(1) logical :: found, defined call get_variable_line(CS, varname, found, defined, value_string) if (found .and. defined .and. (LEN_TRIM(value_string(1)) > 0)) then + if (present(set)) set = .true. read(value_string(1),*,end=991,err=1002) value 991 return else @@ -676,6 +683,7 @@ subroutine read_param_int_array(CS, varname, value, fail_if_missing) ' found but not set in input files.') endif endif ; endif + if (present(set)) set = .false. endif return 1002 call MOM_error(FATAL,'read_param_int_array: read error for integer array '//trim(varname)// & @@ -683,7 +691,7 @@ subroutine read_param_int_array(CS, varname, value, fail_if_missing) end subroutine read_param_int_array !> This subroutine reads the value of a real model parameter from a parameter file. -subroutine read_param_real(CS, varname, value, fail_if_missing, scale) +subroutine read_param_real(CS, varname, value, fail_if_missing, scale, set) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: varname !< The case-sensitive name of the parameter to read @@ -693,6 +701,8 @@ subroutine read_param_real(CS, varname, value, fail_if_missing, scale) !! if this variable is not found in the parameter file real, optional, intent(in) :: scale !< A scaling factor that the parameter is multiplied !! by before it is returned. + logical, optional, intent(out) :: set !< If present, this indicates whether this parameter + !! has been found and successfully set in the input files. ! Local variables character(len=CS%max_line_len) :: value_string(1) @@ -702,6 +712,7 @@ subroutine read_param_real(CS, varname, value, fail_if_missing, scale) if (found .and. defined .and. (LEN_TRIM(value_string(1)) > 0)) then read(value_string(1),*,err=1003) value if (present(scale)) value = scale*value + if (present(set)) set = .true. else if (present(fail_if_missing)) then ; if (fail_if_missing) then if (.not.found) then @@ -712,6 +723,7 @@ subroutine read_param_real(CS, varname, value, fail_if_missing, scale) ' found but not set in input files.') endif endif ; endif + if (present(set)) set = .false. endif return 1003 call MOM_error(FATAL,'read_param_real: read error for real variable '//trim(varname)// & @@ -719,7 +731,7 @@ subroutine read_param_real(CS, varname, value, fail_if_missing, scale) end subroutine read_param_real !> This subroutine reads the values of an array of real model parameters from a parameter file. -subroutine read_param_real_array(CS, varname, value, fail_if_missing, scale) +subroutine read_param_real_array(CS, varname, value, fail_if_missing, scale, set) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: varname !< The case-sensitive name of the parameter to read @@ -729,6 +741,8 @@ subroutine read_param_real_array(CS, varname, value, fail_if_missing, scale) !! if this variable is not found in the parameter file real, optional, intent(in) :: scale !< A scaling factor that the parameter is multiplied !! by before it is returned. + logical, optional, intent(out) :: set !< If present, this indicates whether this parameter + !! has been found and successfully set in the input files. ! Local variables character(len=CS%max_line_len) :: value_string(1) @@ -739,7 +753,7 @@ subroutine read_param_real_array(CS, varname, value, fail_if_missing, scale) read(value_string(1),*,end=991,err=1004) value 991 continue if (present(scale)) value(:) = scale*value(:) - return + if (present(set)) set = .true. else if (present(fail_if_missing)) then ; if (fail_if_missing) then if (.not.found) then @@ -750,6 +764,7 @@ subroutine read_param_real_array(CS, varname, value, fail_if_missing, scale) ' found but not set in input files.') endif endif ; endif + if (present(set)) set = .false. endif return 1004 call MOM_error(FATAL,'read_param_real_array: read error for real array '//trim(varname)// & @@ -757,7 +772,7 @@ subroutine read_param_real_array(CS, varname, value, fail_if_missing, scale) end subroutine read_param_real_array !> This subroutine reads the value of a character string model parameter from a parameter file. -subroutine read_param_char(CS, varname, value, fail_if_missing) +subroutine read_param_char(CS, varname, value, fail_if_missing, set) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: varname !< The case-sensitive name of the parameter to read @@ -765,6 +780,8 @@ subroutine read_param_char(CS, varname, value, fail_if_missing) !! read from the parameter file logical, optional, intent(in) :: fail_if_missing !< If present and true, a fatal error occurs !! if this variable is not found in the parameter file + logical, optional, intent(out) :: set !< If present, this indicates whether this parameter + !! has been found and successfully set in the input files. ! Local variables character(len=CS%max_line_len) :: value_string(1) logical :: found, defined @@ -776,10 +793,12 @@ subroutine read_param_char(CS, varname, value, fail_if_missing) call MOM_error(FATAL, 'Unable to find variable '//trim(varname)//' in any input files.') endif ; endif + if (present(set)) set = found + end subroutine read_param_char !> This subroutine reads the values of an array of character string model parameters from a parameter file. -subroutine read_param_char_array(CS, varname, value, fail_if_missing) +subroutine read_param_char_array(CS, varname, value, fail_if_missing, set) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: varname !< The case-sensitive name of the parameter to read @@ -787,6 +806,8 @@ subroutine read_param_char_array(CS, varname, value, fail_if_missing) !! read from the parameter file logical, optional, intent(in) :: fail_if_missing !< If present and true, a fatal error occurs !! if this variable is not found in the parameter file + logical, optional, intent(out) :: set !< If present, this indicates whether this parameter + !! has been found and successfully set in the input files. ! Local variables character(len=CS%max_line_len) :: value_string(1), loc_string @@ -813,10 +834,12 @@ subroutine read_param_char_array(CS, varname, value, fail_if_missing) call MOM_error(FATAL, 'Unable to find variable '//trim(varname)//' in any input files.') endif ; endif + if (present(set)) set = found + end subroutine read_param_char_array !> This subroutine reads the value of a logical model parameter from a parameter file. -subroutine read_param_logical(CS, varname, value, fail_if_missing) +subroutine read_param_logical(CS, varname, value, fail_if_missing, set) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: varname !< The case-sensitive name of the parameter to read @@ -824,6 +847,8 @@ subroutine read_param_logical(CS, varname, value, fail_if_missing) !! read from the parameter file logical, optional, intent(in) :: fail_if_missing !< If present and true, a fatal error occurs !! if this variable is not found in the parameter file + logical, optional, intent(out) :: set !< If present, this indicates whether this parameter + !! has been found and successfully set in the input files. ! Local variables character(len=CS%max_line_len) :: value_string(1) @@ -835,10 +860,13 @@ subroutine read_param_logical(CS, varname, value, fail_if_missing) elseif (present(fail_if_missing)) then ; if (fail_if_missing) then call MOM_error(FATAL, 'Unable to find variable '//trim(varname)//' in any input files.') endif ; endif + + if (present(set)) set = found + end subroutine read_param_logical !> This subroutine reads the value of a time_type model parameter from a parameter file. -subroutine read_param_time(CS, varname, value, timeunit, fail_if_missing, date_format) +subroutine read_param_time(CS, varname, value, timeunit, fail_if_missing, date_format, set) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: varname !< The case-sensitive name of the parameter to read @@ -850,6 +878,8 @@ subroutine read_param_time(CS, varname, value, timeunit, fail_if_missing, date_f logical, optional, intent(out) :: date_format !< If present, this indicates whether this !! parameter was read in a date format, so that it can !! later be logged in the same format. + logical, optional, intent(out) :: set !< If present, this indicates whether this parameter + !! has been found and successfully set in the input files. ! Local variables character(len=CS%max_line_len) :: value_string(1) @@ -891,6 +921,7 @@ subroutine read_param_time(CS, varname, value, timeunit, fail_if_missing, date_f read( value_string(1), *) real_time value = real_to_time(real_time*time_unit) endif + if (present(set)) set = .true. else if (present(fail_if_missing)) then ; if (fail_if_missing) then if (.not.found) then @@ -899,6 +930,7 @@ subroutine read_param_time(CS, varname, value, timeunit, fail_if_missing, date_f call MOM_error(FATAL, 'Variable '//trim(varname)//' found but not set in input files.') endif endif ; endif + if (present(set)) set = .false. endif return @@ -1704,7 +1736,7 @@ end function convert_date_to_string !! and logs it in documentation files. subroutine get_param_int(CS, modulename, varname, value, desc, units, & default, fail_if_missing, do_not_read, do_not_log, & - layoutParam, debuggingParam) + layoutParam, debuggingParam, old_name) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: modulename !< The name of the calling module @@ -1725,15 +1757,37 @@ subroutine get_param_int(CS, modulename, varname, value, desc, units, & !! logged in the layout parameter file logical, optional, intent(in) :: debuggingParam !< If present and true, this parameter is !! logged in the debugging parameter file + character(len=*), optional, intent(in) :: old_name !< A case-sensitive archaic name of the parameter + !! to read. Errors or warnings are issued if the old name + !! is being used. + ! Local variables logical :: do_read, do_log + logical :: new_name_used, old_name_used, same_value + integer :: new_name_value ! The value that is set when the standard name is used. do_read = .true. ; if (present(do_not_read)) do_read = .not.do_not_read do_log = .true. ; if (present(do_not_log)) do_log = .not.do_not_log if (do_read) then if (present(default)) value = default - call read_param_int(CS, varname, value, fail_if_missing) + + old_name_used = .false. + if (present(old_name)) then + new_name_value = value + call read_param_int(CS, old_name, value, set=old_name_used) + if (old_name_used) then + call read_param_int(CS, varname, new_name_value, set=new_name_used) + + ! Issue appropriate warnings or error messages. + same_value = (value == new_name_value) + call archaic_param_name_message(varname, old_name, new_name_used, same_value) + endif + endif + + if (.not.old_name_used) then ! Old name is either not present or not set. + call read_param_int(CS, varname, value, fail_if_missing) + endif endif if (do_log) then @@ -1747,7 +1801,7 @@ end subroutine get_param_int !! and logs them in documentation files. subroutine get_param_int_array(CS, modulename, varname, value, desc, units, & default, defaults, fail_if_missing, do_not_read, do_not_log, & - layoutParam, debuggingParam) + layoutParam, debuggingParam, old_name) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: modulename !< The name of the calling module @@ -1769,8 +1823,15 @@ subroutine get_param_int_array(CS, modulename, varname, value, desc, units, & !! logged in the layout parameter file logical, optional, intent(in) :: debuggingParam !< If present and true, this parameter is !! logged in the debugging parameter file + character(len=*), optional, intent(in) :: old_name !< A case-sensitive archaic name of the parameter + !! to read. Errors or warnings are issued if the old name + !! is being used. + ! Local variables logical :: do_read, do_log + logical :: new_name_used, old_name_used, same_value + integer :: new_name_value(size(value)) ! The values that are set when the old name is used. + integer :: m do_read = .true. ; if (present(do_not_read)) do_read = .not.do_not_read do_log = .true. ; if (present(do_not_log)) do_log = .not.do_not_log @@ -1785,7 +1846,24 @@ subroutine get_param_int_array(CS, modulename, varname, value, desc, units, & if (do_read) then if (present(default)) value(:) = default if (present(defaults)) value(:) = defaults(:) - call read_param_int_array(CS, varname, value, fail_if_missing) + + old_name_used = .false. + if (present(old_name)) then + new_name_value(:) = value(:) + call read_param_int_array(CS, old_name, value, set=old_name_used) + if (old_name_used) then + call read_param_int_array(CS, varname, new_name_value, set=new_name_used) + + ! Issue appropriate warnings or error messages. + same_value = .true. + do m=1,size(value) ; if (value(m) /= new_name_value(m)) same_value = .false. ; enddo + call archaic_param_name_message(varname, old_name, new_name_used, same_value) + endif + endif + + if (.not.old_name_used) then ! Old name is either not present or not set. + call read_param_int_array(CS, varname, value, fail_if_missing) + endif endif if (do_log) then @@ -1799,7 +1877,7 @@ end subroutine get_param_int_array !! and logs it in documentation files. subroutine get_param_real(CS, modulename, varname, value, desc, units, & default, fail_if_missing, do_not_read, do_not_log, & - debuggingParam, scale, unscaled) + debuggingParam, scale, unscaled, old_name) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: modulename !< The name of the calling module @@ -1822,15 +1900,37 @@ subroutine get_param_real(CS, modulename, varname, value, desc, units, & !! multiplied by before it is returned. real, optional, intent(out) :: unscaled !< The value of the parameter that would be !! returned without any multiplication by a scaling factor. + character(len=*), optional, intent(in) :: old_name !< A case-sensitive archaic name of the parameter + !! to read. Errors or warnings are issued if the old name + !! is being used. + ! Local variables logical :: do_read, do_log + logical :: new_name_used, old_name_used, same_value + real :: new_name_value ! The value that is set when the old name is used. do_read = .true. ; if (present(do_not_read)) do_read = .not.do_not_read do_log = .true. ; if (present(do_not_log)) do_log = .not.do_not_log if (do_read) then if (present(default)) value = default - call read_param_real(CS, varname, value, fail_if_missing) + + old_name_used = .false. + if (present(old_name)) then + new_name_value = value + call read_param_real(CS, old_name, value, set=old_name_used) + if (old_name_used) then + call read_param_real(CS, varname, new_name_value, set=new_name_used) + + ! Issue appropriate warnings or error messages. + same_value = (new_name_used .and. old_name_used .and. (value == new_name_value)) + call archaic_param_name_message(varname, old_name, new_name_used, same_value) + endif + endif + + if (.not.old_name_used) then ! Old name is either not present or not set. + call read_param_real(CS, varname, value, fail_if_missing) + endif endif if (do_log) then @@ -1847,7 +1947,7 @@ end subroutine get_param_real !! and logs them in documentation files. subroutine get_param_real_array(CS, modulename, varname, value, desc, units, & default, defaults, fail_if_missing, do_not_read, do_not_log, debuggingParam, & - scale, unscaled) + scale, unscaled, old_name) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: modulename !< The name of the calling module @@ -1871,8 +1971,15 @@ subroutine get_param_real_array(CS, modulename, varname, value, desc, units, & !! multiplied by before it is returned. real, dimension(:), optional, intent(out) :: unscaled !< The value of the parameter that would be !! returned without any multiplication by a scaling factor. + character(len=*), optional, intent(in) :: old_name !< A case-sensitive archaic name of the parameter + !! to read. Errors or warnings are issued if the old name + !! is being used. + ! Local variables logical :: do_read, do_log + logical :: new_name_used, old_name_used, same_value + real :: new_name_value(size(value)) ! The values that are set when the standard name is used. + integer :: m do_read = .true. ; if (present(do_not_read)) do_read = .not.do_not_read do_log = .true. ; if (present(do_not_log)) do_log = .not.do_not_log @@ -1887,7 +1994,24 @@ subroutine get_param_real_array(CS, modulename, varname, value, desc, units, & if (do_read) then if (present(default)) value(:) = default if (present(defaults)) value(:) = defaults(:) - call read_param_real_array(CS, varname, value, fail_if_missing) + + old_name_used = .false. + if (present(old_name)) then + new_name_value(:) = value(:) + call read_param_real_array(CS, old_name, value, set=old_name_used) + if (old_name_used) then + call read_param_real_array(CS, varname, new_name_value, set=new_name_used) + + ! Issue appropriate warnings or error messages. + same_value = .true. + do m=1,size(value) ; if (value(m) /= new_name_value(m)) same_value = .false. ; enddo + call archaic_param_name_message(varname, old_name, new_name_used, same_value) + endif + endif + + if (.not.old_name_used) then ! Old name is either not present or not set. + call read_param_real_array(CS, varname, value, fail_if_missing) + endif endif if (do_log) then @@ -1904,7 +2028,7 @@ end subroutine get_param_real_array !! and logs it in documentation files. subroutine get_param_char(CS, modulename, varname, value, desc, units, & default, fail_if_missing, do_not_read, do_not_log, & - layoutParam, debuggingParam) + layoutParam, debuggingParam, old_name) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: modulename !< The name of the calling module @@ -1925,15 +2049,37 @@ subroutine get_param_char(CS, modulename, varname, value, desc, units, & !! logged in the layout parameter file logical, optional, intent(in) :: debuggingParam !< If present and true, this parameter is !! logged in the debugging parameter file + character(len=*), optional, intent(in) :: old_name !< A case-sensitive archaic name of the parameter + !! to read. Errors or warnings are issued if the old name + !! is being used. + ! Local variables logical :: do_read, do_log + logical :: new_name_used, old_name_used, same_value + character(len=:), allocatable :: new_name_value ! The value that is set when the standard name is used. do_read = .true. ; if (present(do_not_read)) do_read = .not.do_not_read do_log = .true. ; if (present(do_not_log)) do_log = .not.do_not_log if (do_read) then if (present(default)) value = default - call read_param_char(CS, varname, value, fail_if_missing) + + old_name_used = .false. + if (present(old_name)) then + new_name_value = value + call read_param_char(CS, old_name, value, set=old_name_used) + if (old_name_used) then + call read_param_char(CS, varname, new_name_value, set=new_name_used) + + ! Issue appropriate warnings or error messages. + same_value = (trim(value) == trim(new_name_value)) + call archaic_param_name_message(varname, old_name, new_name_used, same_value) + endif + endif + + if (.not.old_name_used) then ! Old name is either not present or not set. + call read_param_char(CS, varname, value, fail_if_missing) + endif endif if (do_log) then @@ -1946,7 +2092,7 @@ end subroutine get_param_char !> This subroutine reads the values of an array of character string model parameters !! from a parameter file and logs them in documentation files. subroutine get_param_char_array(CS, modulename, varname, value, desc, units, & - default, fail_if_missing, do_not_read, do_not_log) + default, fail_if_missing, do_not_read, do_not_log, old_name) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: modulename !< The name of the calling module @@ -1963,18 +2109,40 @@ subroutine get_param_char_array(CS, modulename, varname, value, desc, units, & !! value for this parameter, although it might be logged. logical, optional, intent(in) :: do_not_log !< If present and true, do not log this !! parameter to the documentation files + character(len=*), optional, intent(in) :: old_name !< A case-sensitive archaic name of the parameter + !! to read. Errors or warnings are issued if the old name + !! is being used. ! Local variables logical :: do_read, do_log - integer :: i, len_tot, len_val + logical :: new_name_used, old_name_used, same_value + integer :: i, m, len_tot, len_val character(len=:), allocatable :: cat_val + character(len=:), allocatable :: new_name_value(:) ! The value that is set when the standard name is used. do_read = .true. ; if (present(do_not_read)) do_read = .not.do_not_read do_log = .true. ; if (present(do_not_log)) do_log = .not.do_not_log if (do_read) then if (present(default)) value(:) = default - call read_param_char_array(CS, varname, value, fail_if_missing) + + old_name_used = .false. + if (present(old_name)) then + new_name_value(:) = value(:) + call read_param_char_array(CS, old_name, value, set=old_name_used) + if (old_name_used) then + call read_param_char_array(CS, varname, new_name_value, set=new_name_used) + + ! Issue appropriate warnings or error messages. + same_value = .true. + do m=1,size(value) ; if (trim(value(m)) /= trim(new_name_value(m))) same_value = .false. ; enddo + call archaic_param_name_message(varname, old_name, new_name_used, same_value) + endif + endif + + if (.not.old_name_used) then ! Old name is either not present or not set. + call read_param_char_array(CS, varname, value, fail_if_missing) + endif endif if (do_log) then @@ -1996,7 +2164,7 @@ end subroutine get_param_char_array !! and logs it in documentation files. subroutine get_param_logical(CS, modulename, varname, value, desc, units, & default, fail_if_missing, do_not_read, do_not_log, & - layoutParam, debuggingParam) + layoutParam, debuggingParam, old_name) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: modulename !< The name of the calling module @@ -2017,15 +2185,37 @@ subroutine get_param_logical(CS, modulename, varname, value, desc, units, & !! logged in the layout parameter file logical, optional, intent(in) :: debuggingParam !< If present and true, this parameter is !! logged in the debugging parameter file + character(len=*), optional, intent(in) :: old_name !< A case-sensitive archaic name of the parameter + !! to read. Errors or warnings are issued if the old name + !! is being used. + ! Local variables logical :: do_read, do_log + logical :: new_name_used, old_name_used, same_value + logical :: new_name_value ! The value that is set when the standard name is used. do_read = .true. ; if (present(do_not_read)) do_read = .not.do_not_read do_log = .true. ; if (present(do_not_log)) do_log = .not.do_not_log if (do_read) then if (present(default)) value = default - call read_param_logical(CS, varname, value, fail_if_missing) + + old_name_used = .false. + if (present(old_name)) then + new_name_value = value + call read_param_logical(CS, old_name, value, set=old_name_used) + if (old_name_used) then + call read_param_logical(CS, varname, new_name_value, set=new_name_used) + + ! Issue appropriate warnings or error messages. + same_value = (value .eqv. new_name_value) + call archaic_param_name_message(varname, old_name, new_name_used, same_value) + endif + endif + + if (.not.old_name_used) then ! Old name is either not present or not set. + call read_param_logical(CS, varname, value, fail_if_missing) + endif endif if (do_log) then @@ -2040,7 +2230,7 @@ end subroutine get_param_logical subroutine get_param_time(CS, modulename, varname, value, desc, units, & default, fail_if_missing, do_not_read, do_not_log, & timeunit, layoutParam, debuggingParam, & - log_as_date) + log_as_date, old_name) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: modulename !< The name of the calling module @@ -2065,8 +2255,14 @@ subroutine get_param_time(CS, modulename, varname, value, desc, units, & !! logged in the debugging parameter file logical, optional, intent(in) :: log_as_date !< If true, log the time_type in date !! format. The default is false. + character(len=*), optional, intent(in) :: old_name !< A case-sensitive archaic name of the parameter + !! to read. Errors or warnings are issued if the old name + !! is being used. + ! Local variables logical :: do_read, do_log, log_date + logical :: new_name_used, old_name_used, same_value + type(time_type) :: new_name_value ! The value that is set when the standard name is used. do_read = .true. ; if (present(do_not_read)) do_read = .not.do_not_read do_log = .true. ; if (present(do_not_log)) do_log = .not.do_not_log @@ -2074,7 +2270,23 @@ subroutine get_param_time(CS, modulename, varname, value, desc, units, & if (do_read) then if (present(default)) value = default - call read_param_time(CS, varname, value, timeunit, fail_if_missing, date_format=log_date) + + old_name_used = .false. + if (present(old_name)) then + new_name_value = value + call read_param_time(CS, old_name, value, timeunit, date_format=log_date, set=old_name_used) + if (old_name_used) then + call read_param_time(CS, varname, new_name_value, timeunit, date_format=log_date, set=new_name_used) + + ! Issue appropriate warnings or error messages. + same_value = (value == new_name_value) + call archaic_param_name_message(varname, old_name, new_name_used, same_value) + endif + endif + + if (.not.old_name_used) then ! Old name is either not present or not set. + call read_param_time(CS, varname, value, timeunit, fail_if_missing, date_format=log_date) + endif endif if (do_log) then @@ -2086,6 +2298,28 @@ subroutine get_param_time(CS, modulename, varname, value, desc, units, & end subroutine get_param_time +!> Issue error messages or warnings about the use of an archaic parameter name. +subroutine archaic_param_name_message(varname, old_name, new_name_used, same_value) + character(len=*), intent(in) :: varname !< The case-sensitive name of the parameter to read + character(len=*), intent(in) :: old_name !< The case-sensitive archaic name of the parameter + logical, intent(in) :: new_name_used !< True if varname is used in the parameter file. + logical, intent(in) :: same_value !< True if varname and old_name give the same values. + + if (new_name_used .and. same_value) then + call MOM_error(WARNING, "The runtime parameter "//trim(varname)//& + " is also being set consistently via its older name of "//trim(old_name)//& + ". Please migrate to only using "//trim(varname)//".") + elseif (new_name_used .and. .not.same_value) then + call MOM_error(FATAL, "The runtime parameter "//trim(varname)//& + " is also being set inconsistently via its older name of "//trim(old_name)//& + ". Only use "//trim(varname)//".") + else + call MOM_error(WARNING, "The runtime parameter "//trim(varname)//& + " is being set via its soon to be obsolete name of "//trim(old_name)//& + ". Please migrate to using "//trim(varname)//".") + endif +end subroutine archaic_param_name_message + ! ----------------------------------------------------------------------------- !> Resets the parameter block name to blank From 400c982b003d49e156951ffc2d97f364352322cd Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 4 Dec 2024 05:39:07 -0500 Subject: [PATCH 3/4] +Add grid_unit_to_L to the ocean_grid_type Add the new element grid_unit_to_L to the ocean_grid_type and the dyn_horgrid_type, which can be used to convert the units of the geoLat and geoLon variables to rescaled horizontal distance units ([L ~> m]) when they are Cartesian coordinates. When Cartesian coordinates are not in use, G%grid_unit_to_L is set to 0. This new element of the grid type is used to test for inconsistent grids or to eliminate rescaling variables in set_rotation_beta_plane(), initialize_velocity_circular(), DOME_initialize_topography(), DOME_initialize_sponges(), DOME_set_OBC_data(), ISOMIP_initialize_topography(), idealized_hurricane_wind_forcing(), Kelvin_set_OBC_data(), Rossby_front_initialize_velocity(), soliton_initialize_thickness(), and soliton_initialize_velocity(). These are the instances where this new variable could be used and bitwise identical answers are recovered. There are a few other places where they should be used, but where answers would change, and these will be deferred to a subsequent commit. All answers are bitwise identical, but there are new elements in two transparent and widely used types. --- src/core/MOM_grid.F90 | 4 ++ src/core/MOM_transcribe_grid.F90 | 2 + src/framework/MOM_dyn_horgrid.F90 | 5 ++ src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 | 3 + src/initialization/MOM_grid_initialize.F90 | 9 ++- .../MOM_shared_initialization.F90 | 9 +-- .../MOM_state_initialization.F90 | 7 ++- src/user/DOME_initialization.F90 | 60 +++++++++++++------ src/user/ISOMIP_initialization.F90 | 18 +++--- src/user/Idealized_Hurricane.F90 | 13 ++-- src/user/Kelvin_initialization.F90 | 21 ++++--- src/user/Rossby_front_2d_initialization.F90 | 8 +-- src/user/soliton_initialization.F90 | 8 ++- 13 files changed, 110 insertions(+), 57 deletions(-) diff --git a/src/core/MOM_grid.F90 b/src/core/MOM_grid.F90 index 6fb8426395..5e5b307103 100644 --- a/src/core/MOM_grid.F90 +++ b/src/core/MOM_grid.F90 @@ -190,6 +190,10 @@ module MOM_grid ! These parameters are run-time parameters that are used during some ! initialization routines (but not all) + real :: grid_unit_to_L !< A factor that converts a the geoLat and geoLon variables and related + !! variables like len_lat and len_lon into rescaled horizontal distance + !! units on a Cartesian grid, in [L km ~> 1000] or [L m-1 ~> 1] or + !! is 0 for a non-Cartesian grid. real :: south_lat !< The latitude (or y-coordinate) of the first v-line [degrees_N] or [km] or [m] real :: west_lon !< The longitude (or x-coordinate) of the first u-line [degrees_E] or [km] or [m] real :: len_lat !< The latitudinal (or y-coord) extent of physical domain [degrees_N] or [km] or [m] diff --git a/src/core/MOM_transcribe_grid.F90 b/src/core/MOM_transcribe_grid.F90 index f8ae58d9e1..4c1b87f37e 100644 --- a/src/core/MOM_transcribe_grid.F90 +++ b/src/core/MOM_transcribe_grid.F90 @@ -135,6 +135,7 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG, US) ! Copy various scalar variables and strings. oG%x_axis_units = dG%x_axis_units ; oG%y_axis_units = dG%y_axis_units oG%x_ax_unit_short = dG%x_ax_unit_short ; oG%y_ax_unit_short = dG%y_ax_unit_short + oG%grid_unit_to_L = dG%grid_unit_to_L oG%areaT_global = dG%areaT_global ; oG%IareaT_global = dG%IareaT_global oG%south_lat = dG%south_lat ; oG%west_lon = dG%west_lon oG%len_lat = dG%len_lat ; oG%len_lon = dG%len_lon @@ -296,6 +297,7 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG, US) ! Copy various scalar variables and strings. dG%x_axis_units = oG%x_axis_units ; dG%y_axis_units = oG%y_axis_units dG%x_ax_unit_short = oG%x_ax_unit_short ; dG%y_ax_unit_short = oG%y_ax_unit_short + dG%grid_unit_to_L = oG%grid_unit_to_L dG%areaT_global = oG%areaT_global ; dG%IareaT_global = oG%IareaT_global dG%south_lat = oG%south_lat ; dG%west_lon = oG%west_lon dG%len_lat = oG%len_lat ; dG%len_lon = oG%len_lon diff --git a/src/framework/MOM_dyn_horgrid.F90 b/src/framework/MOM_dyn_horgrid.F90 index 987d5bf502..440ea351b9 100644 --- a/src/framework/MOM_dyn_horgrid.F90 +++ b/src/framework/MOM_dyn_horgrid.F90 @@ -181,6 +181,10 @@ module MOM_dyn_horgrid ! These parameters are run-time parameters that are used during some ! initialization routines (but not all) + real :: grid_unit_to_L !< A factor that converts a the geoLat and geoLon variables and related + !! variables like len_lat and len_lon into rescaled horizontal distance + !! units on a Cartesian grid, in [L km ~> 1000] or [L m-1 ~> 1] or + !! is 0 for a non-Cartesian grid. real :: south_lat !< The latitude (or y-coordinate) of the first v-line [degrees_N] or [km] or [m] real :: west_lon !< The longitude (or x-coordinate) of the first u-line [degrees_E] or [km] or [m] real :: len_lat !< The latitudinal (or y-coord) extent of physical domain [degrees_N] or [km] or [m] @@ -400,6 +404,7 @@ subroutine rotate_dyn_horgrid(G_in, G, US, turns) G%len_lon = G_in%len_lon ! Rotation-invariant fields + G%grid_unit_to_L = G_in%grid_unit_to_L G%areaT_global = G_in%areaT_global G%IareaT_global = G_in%IareaT_global G%Rad_Earth = G_in%Rad_Earth diff --git a/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 b/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 index 5ecfb9e788..fe54dd6533 100644 --- a/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 +++ b/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 @@ -132,6 +132,7 @@ subroutine set_IS_axes_info(G, param_file, diag_cs, axes_set_name) G%x_axis_units = "degrees_E" ; G%y_axis_units = "degrees_N" G%x_ax_unit_short = "degrees_E" ; G%y_ax_unit_short = "degrees_N" + G%grid_unit_to_L = 0.0 if (index(lowercase(trim(grid_config)),"cartesian") > 0) then ! This is a cartesian grid, and may have different axis units. @@ -145,9 +146,11 @@ subroutine set_IS_axes_info(G, param_file, diag_cs, axes_set_name) if (units_temp(1:1) == 'k') then G%x_axis_units = "kilometers" ; G%y_axis_units = "kilometers" G%x_ax_unit_short = "km" ; G%y_ax_unit_short = "km" + G%grid_unit_to_L = 1000.0*diag_cs%US%m_to_L elseif (units_temp(1:1) == 'm') then G%x_axis_units = "meters" ; G%y_axis_units = "meters" G%x_ax_unit_short = "m" ; G%y_ax_unit_short = "m" + G%grid_unit_to_L = diag_cs%US%m_to_L endif call log_param(param_file, mdl, "explicit AXIS_UNITS", G%x_axis_units) else diff --git a/src/initialization/MOM_grid_initialize.F90 b/src/initialization/MOM_grid_initialize.F90 index ef78a896c3..429ce24d79 100644 --- a/src/initialization/MOM_grid_initialize.F90 +++ b/src/initialization/MOM_grid_initialize.F90 @@ -84,7 +84,7 @@ subroutine set_grid_metrics(G, param_file, US) ! These are defaults that may be changed in the next select block. G%x_axis_units = "degrees_east" ; G%y_axis_units = "degrees_north" G%x_ax_unit_short = "degrees_E" ; G%y_ax_unit_short = "degrees_N" - + G%grid_unit_to_L = 0.0 G%Rad_Earth_L = -1.0*US%m_to_L ; G%len_lat = 0.0 ; G%len_lon = 0.0 select case (trim(config)) case ("mosaic"); call set_grid_metrics_from_mosaic(G, param_file, US) @@ -251,6 +251,11 @@ subroutine set_grid_metrics_from_mosaic(G, param_file, US) G%geoLatCv(i,J) = tmpZ(i2-1,j2) enddo ; enddo + ! This routine could be modified to support the use of a mosaic using Cartesian grid coordinates, + ! in which case the values of G%x_axis_units, G%y_axis_units and G%grid_unit_to_L would need to be + ! reset appropriately here, but this option has not yet been implemented, and the grid coordinates + ! are assumed to be degrees of longitude and latitude. + ! Read DX,DY from the supergrid tmpU(:,:) = 0. ; tmpV(:,:) = 0. call MOM_read_data(filename, 'dx', tmpV, SGdom, position=NORTH_FACE, scale=US%m_to_L) @@ -440,9 +445,11 @@ subroutine set_grid_metrics_cartesian(G, param_file, US) enddo if (units_temp(1:1) == 'k') then ! Axes are measured in km. + G%grid_unit_to_L = 1000.0*US%m_to_L dx_everywhere = 1000.0*US%m_to_L * G%len_lon / (REAL(niglobal)) dy_everywhere = 1000.0*US%m_to_L * G%len_lat / (REAL(njglobal)) elseif (units_temp(1:1) == 'm') then ! Axes are measured in m. + G%grid_unit_to_L = US%m_to_L dx_everywhere = US%m_to_L*G%len_lon / (REAL(niglobal)) dy_everywhere = US%m_to_L*G%len_lat / (REAL(njglobal)) else ! Axes are measured in degrees of latitude and longitude. diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90 index 0a257b6ca1..1f7a519d9e 100644 --- a/src/initialization/MOM_shared_initialization.F90 +++ b/src/initialization/MOM_shared_initialization.F90 @@ -485,7 +485,6 @@ subroutine set_rotation_beta_plane(f, G, param_file, US) real :: f_0 ! The reference value of the Coriolis parameter [T-1 ~> s-1] real :: beta ! The meridional gradient of the Coriolis parameter [T-1 L-1 ~> s-1 m-1] real :: beta_lat_ref ! The reference latitude for the beta plane [degrees_N] or [km] or [m] - real :: Rad_Earth_L ! The radius of the planet in rescaled units [L ~> m] real :: y_scl ! A scaling factor from the units of latitude [L lat-1 ~> m lat-1] real :: PI ! The ratio of the circumference of a circle to its diameter [nondim] character(len=40) :: mdl = "set_rotation_beta_plane" ! This subroutine's name. @@ -503,18 +502,16 @@ subroutine set_rotation_beta_plane(f, G, param_file, US) call get_param(param_file, mdl, "AXIS_UNITS", axis_units, default="degrees") PI = 4.0*atan(1.0) + y_scl = G%grid_unit_to_L + if (G%grid_unit_to_L <= 0.0) y_scl = PI * G%Rad_Earth_L / 180. + select case (axis_units(1:1)) case ("d") - call get_param(param_file, mdl, "RAD_EARTH", Rad_Earth_L, & - "The radius of the Earth.", units="m", default=6.378e6, scale=US%m_to_L) beta_lat_ref_units = "degrees" - y_scl = PI * Rad_Earth_L / 180. case ("k") beta_lat_ref_units = "kilometers" - y_scl = 1.0e3 * US%m_to_L case ("m") beta_lat_ref_units = "meters" - y_scl = 1.0 * US%m_to_L case default ; call MOM_error(FATAL, & " set_rotation_beta_plane: unknown AXIS_UNITS = "//trim(axis_units)) end select diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 7909bddb80..0f5d06aa4e 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -1636,7 +1636,10 @@ subroutine initialize_velocity_circular(u, v, G, GV, US, param_file, just_read) if (just_read) return ! All run-time parameters have been read, so return. - dpi=acos(0.0)*2.0 ! pi + if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, "MOM_state_initialization.F90: "//& + "initialize_velocity_circular() is only set to work with Cartesian axis units.") + + dpi = acos(0.0)*2.0 ! pi do k=1,nz ; do j=js,je ; do I=Isq,Ieq psi1 = my_psi(I,j) @@ -1663,7 +1666,7 @@ real function my_psi(ig,jg) r = sqrt( (x**2) + (y**2) ) ! Circular stream function is a function of radius only r = min(1.0, r) ! Flatten stream function in corners of box my_psi = 0.5*(1.0 - cos(dpi*r)) - my_psi = my_psi * (circular_max_u * G%US%m_to_L*G%len_lon*1e3 / dpi) ! len_lon is in km + my_psi = my_psi * (circular_max_u * G%len_lon * G%grid_unit_to_L / dpi) ! len_lon is in km end function my_psi end subroutine initialize_velocity_circular diff --git a/src/user/DOME_initialization.F90 b/src/user/DOME_initialization.F90 index 858ca32f93..e608dbd1c2 100644 --- a/src/user/DOME_initialization.F90 +++ b/src/user/DOME_initialization.F90 @@ -49,10 +49,11 @@ subroutine DOME_initialize_topography(D, G, param_file, max_depth, US) real :: min_depth ! The minimum ocean depth [Z ~> m] real :: shelf_depth ! The ocean depth on the shelf in the DOME configuration [Z ~> m] real :: slope ! The bottom slope in the DOME configuration [Z L-1 ~> nondim] - real :: shelf_edge_lat ! The latitude of the edge of the topographic shelf [km] - real :: inflow_lon ! The edge longitude of the DOME inflow [km] - real :: inflow_width ! The longitudinal width of the DOME inflow channel [km] - real :: km_to_L ! The conversion factor from the units of latitude to L [L km-1 ~> 1e3] + real :: shelf_edge_lat ! The latitude of the edge of the topographic shelf in the same units as geolat, often [km] + real :: inflow_lon ! The edge longitude of the DOME inflow in the same units as geolon, often [km] + real :: inflow_width ! The longitudinal width of the DOME inflow channel in the same units as geolat, often [km] + real :: km_to_grid_unit ! The conversion factor from km to the units of latitude often 1 [nondim], + ! but this could be 1000 [m km-1] ! This include declares and sets the variable "version". # include "version_variable.h" character(len=40) :: mdl = "DOME_initialize_topography" ! This subroutine's name. @@ -60,7 +61,16 @@ subroutine DOME_initialize_topography(D, G, param_file, max_depth, US) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed - km_to_L = 1.0e3*US%m_to_L + if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, "DOME_initialization: "//& + "DOME_initialize_topography is only set to work with Cartesian axis units.") + if (abs(G%grid_unit_to_L*US%L_to_m - 1000.0) < 1.0e-3) then ! The grid latitudes are in km. + km_to_grid_unit = 1.0 + elseif (abs(G%grid_unit_to_L*US%L_to_m - 1.0) < 1.0e-6) then ! The grid latitudes are in m. + km_to_grid_unit = 1000.0 + else + call MOM_error(FATAL, "DOME_initialization: "//& + "DOME_initialize_topography is not recognizing the value of G%grid_unit_to_L.") + endif call MOM_mesg(" DOME_initialization.F90, DOME_initialize_topography: setting topography", 5) @@ -75,15 +85,16 @@ subroutine DOME_initialize_topography(D, G, param_file, max_depth, US) default=600.0, units="m", scale=US%m_to_Z) call get_param(param_file, mdl, "DOME_SHELF_EDGE_LAT", shelf_edge_lat, & "The latitude of the shelf edge in the DOME configuration.", & - default=600.0, units="km") + default=600.0, units="km", scale=km_to_grid_unit) call get_param(param_file, mdl, "DOME_INFLOW_LON", inflow_lon, & - "The edge longitude of the DOME inflow.", units="km", default=1000.0) + "The edge longitude of the DOME inflow.", units="km", default=1000.0, scale=km_to_grid_unit) call get_param(param_file, mdl, "DOME_INFLOW_WIDTH", inflow_width, & - "The longitudinal width of the DOME inflow channel.", units="km", default=100.0) + "The longitudinal width of the DOME inflow channel.", & + units="km", default=100.0, scale=km_to_grid_unit) do j=js,je ; do i=is,ie if (G%geoLatT(i,j) < shelf_edge_lat) then - D(i,j) = min(shelf_depth - slope * (G%geoLatT(i,j)-shelf_edge_lat)*km_to_L, max_depth) + D(i,j) = min(shelf_depth - slope * (G%geoLatT(i,j)-shelf_edge_lat)*G%grid_unit_to_L, max_depth) else if ((G%geoLonT(i,j) > inflow_lon) .AND. (G%geoLonT(i,j) < inflow_lon+inflow_width)) then D(i,j) = shelf_depth @@ -177,7 +188,6 @@ subroutine DOME_initialize_sponges(G, GV, US, tv, depth_tot, PF, CSp) real :: min_depth ! The minimum depth at which to apply damping [Z ~> m] real :: damp_W, damp_E ! Damping rates in the western and eastern sponges [T-1 ~> s-1] real :: peak_damping ! The maximum sponge damping rates as the edges [T-1 ~> s-1] - real :: km_to_L ! The conversion factor from the units of longitude to L [L km-1 ~> 1e3] real :: edge_dist ! The distance to an edge [L ~> m] real :: sponge_width ! The width of the sponges [L ~> m] character(len=40) :: mdl = "DOME_initialize_sponges" ! This subroutine's name. @@ -186,7 +196,8 @@ subroutine DOME_initialize_sponges(G, GV, US, tv, depth_tot, PF, CSp) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed - km_to_L = 1.0e3*US%m_to_L + if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, "DOME_initialization: "//& + "DOME_initialize_sponges is only set to work with Cartesian axis units.") ! Set up sponges for the DOME configuration call get_param(PF, mdl, "MINIMUM_DEPTH", min_depth, & @@ -196,7 +207,7 @@ subroutine DOME_initialize_sponges(G, GV, US, tv, depth_tot, PF, CSp) default=10.0, units="day-1", scale=1.0/(86400.0*US%s_to_T)) call get_param(PF, mdl, "DOME_SPONGE_WIDTH", sponge_width, & "The width of the the DOME sponges.", & - default=200.0, units="km", scale=km_to_L) + default=200.0, units="km", scale=1.0e3*US%m_to_L) ! Here the inverse damping time [T-1 ~> s-1], is set. Set Idamp to 0 wherever ! there is no sponge, and the subroutines that are called will automatically @@ -204,7 +215,7 @@ subroutine DOME_initialize_sponges(G, GV, US, tv, depth_tot, PF, CSp) Idamp(:,:) = 0.0 do j=js,je ; do i=is,ie ; if (depth_tot(i,j) > min_depth) then - edge_dist = (G%geoLonT(i,j) - G%west_lon) * km_to_L + edge_dist = (G%geoLonT(i,j) - G%west_lon) * G%grid_unit_to_L if (edge_dist < 0.5*sponge_width) then damp_W = peak_damping elseif (edge_dist < sponge_width) then @@ -213,7 +224,7 @@ subroutine DOME_initialize_sponges(G, GV, US, tv, depth_tot, PF, CSp) damp_W = 0.0 endif - edge_dist = ((G%len_lon + G%west_lon) - G%geoLonT(i,j)) * km_to_L + edge_dist = ((G%len_lon + G%west_lon) - G%geoLonT(i,j)) * G%grid_unit_to_L if (edge_dist < 0.5*sponge_width) then damp_E = peak_damping elseif (edge_dist < sponge_width) then @@ -328,10 +339,12 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, US, PF, tr_Reg) ! properties [T-1 ~> s-1] real :: g_prime_tot ! The reduced gravity across all layers [L2 Z-1 T-2 ~> m s-2] real :: Def_Rad ! The deformation radius, based on fluid of thickness D_edge [L ~> m] - real :: inflow_lon ! The edge longitude of the DOME inflow [km] + real :: inflow_lon ! The edge longitude of the DOME inflow in the same units as geolon, often [km] real :: I_Def_Rad ! The inverse of the deformation radius in the same units as G%geoLon [km-1] real :: Ri_trans ! The shear Richardson number in the transition ! region of the specified shear profile [nondim] + real :: km_to_grid_unit ! The conversion factor from km to the units of latitude often 1 [nondim], + ! but this could be 1000 [m km-1] character(len=32) :: name ! The name of a tracer field. character(len=40) :: mdl = "DOME_set_OBC_data" ! This subroutine's name. integer :: i, j, k, itt, is, ie, js, je, isd, ied, jsd, jed, m, nz, ntherm, ntr_id @@ -343,6 +356,17 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, US, PF, tr_Reg) isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB + if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, "DOME_initialization: "//& + "DOME_initialize_topography is only set to work with Cartesian axis units.") + if (abs(G%grid_unit_to_L*US%L_to_m - 1000.0) < 1.0e-3) then ! The grid latitudes are in km. + km_to_grid_unit = 1.0 + elseif (abs(G%grid_unit_to_L*US%L_to_m - 1.0) < 1.0e-6) then ! The grid latitudes are in m. + km_to_grid_unit = 1000.0 + else + call MOM_error(FATAL, "DOME_initialization: "//& + "DOME_initialize_topography is not recognizing the value of G%grid_unit_to_L.") + endif + call get_param(PF, mdl, "DOME_INFLOW_THICKNESS", D_edge, & "The thickness of the dense DOME inflow at the inner edge.", & default=300.0, units="m", scale=US%m_to_Z) @@ -362,7 +386,7 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, US, PF, tr_Reg) "The value of the Coriolis parameter that is used to determine the DOME "//& "inflow properties.", units="s-1", default=f_0*US%s_to_T, scale=US%T_to_s) call get_param(PF, mdl, "DOME_INFLOW_LON", inflow_lon, & - "The edge longitude of the DOME inflow.", units="km", default=1000.0) + "The edge longitude of the DOME inflow.", units="km", default=1000.0, scale=km_to_grid_unit) if (associated(tv%S) .or. associated(tv%T)) then call get_param(PF, mdl, "S_REF", S_ref, & units="ppt", default=35.0, scale=US%ppt_to_S, do_not_log=.true.) @@ -383,7 +407,9 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, US, PF, tr_Reg) tr_0 = (-D_edge*sqrt(D_edge*g_prime_tot)*0.5*Def_Rad) * (Rlay_Ref + 0.5*Rlay_range) * GV%RZ_to_H endif - I_Def_Rad = 1.0 / (1.0e-3*US%L_to_m*Def_Rad) + I_Def_Rad = 1.0 / ((1.0e-3*US%L_to_m*km_to_grid_unit) * Def_Rad) + ! This is mathematically equivalent to + ! I_Def_Rad = G%grid_unit_to_L / Def_Rad if (OBC%number_of_segments /= 1) then call MOM_error(WARNING, 'Error in DOME OBC segment setup', .true.) diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index d03a07e313..4bf7931856 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -59,7 +59,6 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth, US) real :: ly ! domain width (across ice flow) [L ~> m] real :: bx, by ! The x- and y- contributions to the bathymetric profiles at a point [Z ~> m] real :: xtil ! x-positon normalized by the characteristic along-flow length scale [nondim] - real :: km_to_L ! The conversion factor from the axis units to L [L km-1 ~> 1e3] logical :: is_2D ! If true, use a 2D setup ! This include declares and sets the variable "version". # include "version_variable.h" @@ -93,7 +92,8 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth, US) "Characteristic width of the side walls of the channel in the ISOMIP configuration.", & units="m", default=4.0e3, scale=US%m_to_L) - km_to_L = 1.0e3*US%m_to_L + if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, "ISOMIP_initialization.F90: " //& + "ISOMIP_initialize_topography is only set to work with Cartesian axis units.") ! The following variables should be transformed into runtime parameters. b0 = -150.0*US%m_to_Z ; b2 = -728.8*US%m_to_Z ; b4 = 343.91*US%m_to_Z ; b6 = -50.57*US%m_to_Z @@ -101,8 +101,8 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth, US) if (is_2D) then do j=js,je ; do i=is,ie ! For the 2D setup take a slice through the middle of the domain - xtil = G%geoLonT(i,j)*km_to_L / xbar - !xtil = 450.*km_to_L / xbar + xtil = G%geoLonT(i,j)*G%grid_unit_to_L / xbar + ! xtil = 450.0e3*US%m_to_L / xbar bx = b0 + b2*xtil**2 + b4*xtil**4 + b6*xtil**6 by = 2.0 * dc / (1.0 + exp(2.0*wc / fc)) @@ -117,17 +117,17 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth, US) ! 3D setup ! ===== TEST ===== !if (G%geoLonT(i,j)<500.) then - ! xtil = 500.*km_to_L / xbar + ! xtil = 500.0e3*US%m_to_L / xbar !else - ! xtil = G%geoLonT(i,j)*km_to_L / xbar + ! xtil = G%geoLonT(i,j)*G%grid_unit_to_L / xbar !endif ! ===== TEST ===== - xtil = G%geoLonT(i,j)*km_to_L / xbar + xtil = G%geoLonT(i,j)*G%grid_unit_to_L / xbar bx = b0 + b2*xtil**2 + b4*xtil**4 + b6*xtil**6 - by = (dc / (1.0 + exp(-2.*(G%geoLatT(i,j)*km_to_L - 0.5*ly - wc) / fc))) + & - (dc / (1.0 + exp(2.*(G%geoLatT(i,j)*km_to_L - 0.5*ly + wc) / fc))) + by = (dc / (1.0 + exp(-2.*(G%geoLatT(i,j)*G%grid_unit_to_L - 0.5*ly - wc) / fc))) + & + (dc / (1.0 + exp(2.*(G%geoLatT(i,j)*G%grid_unit_to_L - 0.5*ly + wc) / fc))) D(i,j) = -max(bx+by, -bmax) if (D(i,j) > max_depth) D(i,j) = max_depth diff --git a/src/user/Idealized_Hurricane.F90 b/src/user/Idealized_Hurricane.F90 index aca19e86d0..99cc7b88c5 100644 --- a/src/user/Idealized_Hurricane.F90 +++ b/src/user/Idealized_Hurricane.F90 @@ -349,7 +349,6 @@ subroutine idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, CS) real :: fbench !< The benchmark 'f' value [T-1 ~> s-1] real :: fbench_fac !< A factor that is set to 0 to use the !! benchmark 'f' value [nondim] - real :: km_to_L !< The conversion factor from the units of latitude to L [L km-1 ~> 1e3] real :: rel_tau_fac !< A factor that is set to 0 to disable !! current relative stress calculation [nondim] @@ -359,7 +358,8 @@ subroutine idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, CS) isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB - km_to_L = 1.0e3*US%m_to_L + if ((G%grid_unit_to_L <= 0.) .and. (.not.CS%SCM_mode)) call MOM_error(FATAL, "Idealized_Hurricane.F90: " //& + "idealized_hurricane_wind_forcing is only set to work with Cartesian axis units.") ! Allocate the forcing arrays, if necessary. call allocate_mech_forcing(G, forces, stress=.true., ustar=.true., tau_mag=.true.) @@ -376,7 +376,6 @@ subroutine idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, CS) YC = CS%Hurr_cen_Y0 + (time_type_to_real(day)*US%s_to_T * CS%hurr_translation_spd * & sin(CS%hurr_translation_dir)) - if (CS%BR_Bench) then ! f reset to value used in generated wind for benchmark test fbench = CS%f_column @@ -403,8 +402,8 @@ subroutine idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, CS) YY = CS%dy_from_center - YC XX = -XC else - YY = G%geoLatCu(I,j)*km_to_L - YC - XX = G%geoLonCu(I,j)*km_to_L - XC + YY = G%geoLatCu(I,j) * G%grid_unit_to_L - YC + XX = G%geoLonCu(I,j) * G%grid_unit_to_L - XC endif call idealized_hurricane_wind_profile(CS, US, f_local, YY, XX, Uocn, Vocn, TX, TY) forces%taux(I,j) = G%mask2dCu(I,j) * TX @@ -427,8 +426,8 @@ subroutine idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, CS) YY = CS%dy_from_center - YC XX = -XC else - YY = G%geoLatCv(i,J)*km_to_L - YC - XX = G%geoLonCv(i,J)*km_to_L - XC + YY = G%geoLatCv(i,J) * G%grid_unit_to_L - YC + XX = G%geoLonCv(i,J) * G%grid_unit_to_L - XC endif call idealized_hurricane_wind_profile(CS, US, f_local, YY, XX, Uocn, Vocn, TX, TY) forces%tauy(i,J) = G%mask2dCv(i,J) * TY diff --git a/src/user/Kelvin_initialization.F90 b/src/user/Kelvin_initialization.F90 index 1fc8a2f564..118d76ac93 100644 --- a/src/user/Kelvin_initialization.F90 +++ b/src/user/Kelvin_initialization.F90 @@ -215,10 +215,11 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time) if (.not.associated(OBC)) call MOM_error(FATAL, 'Kelvin_initialization.F90: '// & 'Kelvin_set_OBC_data() was called but OBC type was not initialized!') + if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, 'Kelvin_initialization.F90: '// & + "Kelvin_set_OBC_data() is only set to work with Cartesian axis units.") time_sec = US%s_to_T*time_type_to_real(Time) PI = 4.0*atan(1.0) - km_to_L_scale = 1000.0*US%m_to_L do j=jsd,jed ; do i=isd,ied depth_tot(i,j) = 0.0 @@ -237,6 +238,8 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time) lambda = PI * CS%mode * CS%F_0 / (CS%H0 * N0) ! Two wavelengths in domain omega = (4.0 * CS%H0 * N0) / (CS%mode * US%m_to_L*G%len_lon) + !### There is a bug here when len_lon is in km. This should be + ! omega = (4.0 * CS%H0 * N0) / (CS%mode * G%grid_unit_to_L*G%len_lon) endif sina = sin(CS%coast_angle) @@ -257,8 +260,8 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time) jsd = segment%HI%jsd ; jed = segment%HI%jed JsdB = segment%HI%JsdB ; JedB = segment%HI%JedB do j=jsd,jed ; do I=IsdB,IedB - x1 = km_to_L_scale * G%geoLonCu(I,j) - y1 = km_to_L_scale * G%geoLatCu(I,j) + x1 = G%grid_unit_to_L * G%geoLonCu(I,j) + y1 = G%grid_unit_to_L * G%geoLatCu(I,j) x = (x1 - CS%coast_offset1) * cosa + y1 * sina y = -(x1 - CS%coast_offset1) * sina + y1 * cosa if (CS%mode == 0) then @@ -299,8 +302,8 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time) enddo ; enddo if (allocated(segment%tangential_vel)) then do J=JsdB+1,JedB-1 ; do I=IsdB,IedB - x1 = km_to_L_scale * G%geoLonBu(I,J) - y1 = km_to_L_scale * G%geoLatBu(I,J) + x1 = G%grid_unit_to_L * G%geoLonBu(I,J) + y1 = G%grid_unit_to_L * G%geoLatBu(I,J) x = (x1 - CS%coast_offset1) * cosa + y1 * sina y = - (x1 - CS%coast_offset1) * sina + y1 * cosa cff = sqrt(GV%g_Earth * depth_tot(i+1,j) ) @@ -316,8 +319,8 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time) isd = segment%HI%isd ; ied = segment%HI%ied JsdB = segment%HI%JsdB ; JedB = segment%HI%JedB do J=JsdB,JedB ; do i=isd,ied - x1 = km_to_L_scale * G%geoLonCv(i,J) - y1 = km_to_L_scale * G%geoLatCv(i,J) + x1 = G%grid_unit_to_L * G%geoLonCv(i,J) + y1 = G%grid_unit_to_L * G%geoLatCv(i,J) x = (x1 - CS%coast_offset1) * cosa + y1 * sina y = - (x1 - CS%coast_offset1) * sina + y1 * cosa if (CS%mode == 0) then @@ -355,8 +358,8 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time) enddo ; enddo if (allocated(segment%tangential_vel)) then do J=JsdB,JedB ; do I=IsdB+1,IedB-1 - x1 = km_to_L_scale * G%geoLonBu(I,J) - y1 = km_to_L_scale * G%geoLatBu(I,J) + x1 = G%grid_unit_to_L * G%geoLonBu(I,J) + y1 = G%grid_unit_to_L * G%geoLatBu(I,J) x = (x1 - CS%coast_offset1) * cosa + y1 * sina y = - (x1 - CS%coast_offset1) * sina + y1 * cosa cff = sqrt(GV%g_Earth * depth_tot(i,j+1) ) diff --git a/src/user/Rossby_front_2d_initialization.F90 b/src/user/Rossby_front_2d_initialization.F90 index 33c7641a00..09219eb845 100644 --- a/src/user/Rossby_front_2d_initialization.F90 +++ b/src/user/Rossby_front_2d_initialization.F90 @@ -250,6 +250,8 @@ subroutine Rossby_front_initialize_velocity(u, v, h, G, GV, US, param_file, just if (max_depth <= 0.0) call MOM_error(FATAL, & "Rossby_front_initialize_thickness, Rossby_front_initialize_velocity: "//& "This module requires a positive value of MAXIMUM_DEPTH.") + if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, 'Rossby_front_2d_initialization.F90: '// & + "dTdy() is only set to work with Cartesian axis units.") v(:,:,:) = 0.0 u(:,:,:) = 0.0 @@ -335,18 +337,16 @@ end function Hml real function dTdy( G, dT, lat, US ) type(ocean_grid_type), intent(in) :: G !< Grid structure real, intent(in) :: dT !< Top to bottom temperature difference [C ~> degC] - real, intent(in) :: lat !< Latitude in [km] + real, intent(in) :: lat !< Latitude in the same units as geoLat, often [km] type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! Local real :: PI ! The ratio of the circumference of a circle to its diameter [nondim] real :: dHML ! The range of the mixed layer depths [Z ~> m] real :: dHdy ! The mixed layer depth gradient [Z L-1 ~> m m-1] - real :: km_to_L ! Horizontal axis unit conversion factor when AXIS_UNITS = 'k' (1000 m) [L km-1 ~> 1000] PI = 4.0 * atan(1.0) - km_to_L = 1.0e3*US%m_to_L dHML = 0.5 * ( HMLmax - HMLmin ) * G%max_depth - dHdy = dHML * ( PI / ( frontFractionalWidth * G%len_lat * km_to_L ) ) * cos( yPseudo(G, lat) ) + dHdy = dHML * ( PI / ( frontFractionalWidth * G%len_lat * G%grid_unit_to_L ) ) * cos( yPseudo(G, lat) ) dTdy = -( dT / G%max_depth ) * dHdy end function dTdy diff --git a/src/user/soliton_initialization.F90 b/src/user/soliton_initialization.F90 index 722a41b7e5..a734574995 100644 --- a/src/user/soliton_initialization.F90 +++ b/src/user/soliton_initialization.F90 @@ -79,10 +79,12 @@ subroutine soliton_initialize_thickness(h, depth_tot, G, GV, US, param_file, jus if (abs(beta) <= 0.0) call MOM_error(FATAL, & "soliton_initialization, soliton_initialize_thickness: "//& "This module requires a non-zero value of BETA.") + if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, "soliton_initialization.F90: "//& + "soliton_initialize_thickness() is only set to work with Cartesian axis units.") cg_max = sqrt(GV%g_Earth * max_depth) L_eq = sqrt(cg_max / abs(beta)) - scale_pos = US%m_to_L / L_eq + scale_pos = G%grid_unit_to_L / L_eq I_nz = 1.0 / real(nz) x0 = 2.0*G%len_lon/3.0 @@ -150,10 +152,12 @@ subroutine soliton_initialize_velocity(u, v, G, GV, US, param_file, just_read) if (abs(beta) <= 0.0) call MOM_error(FATAL, & "soliton_initialization, soliton_initialize_velocity: "//& "This module requires a non-zero value of BETA.") + if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, "soliton_initialization.F90: "//& + "soliton_initialize_velocity() is only set to work with Cartesian axis units.") cg_max = sqrt(GV%g_Earth * max_depth) L_eq = sqrt(cg_max / abs(beta)) - scale_pos = US%m_to_L / L_eq + scale_pos = G%grid_unit_to_L / L_eq x0 = 2.0*G%len_lon/3.0 y0 = 0.0 From 14f2c97aa9ee2f5d508aa281c8b2fc4a0747bcbd Mon Sep 17 00:00:00 2001 From: Cory Spencer Jones Date: Mon, 23 Dec 2024 15:59:40 -0600 Subject: [PATCH 4/4] Update particles_run call to allow accurate particle advection using u, v or uh, vh The drifters package is able to run in 2 different modes: one where the particles are advected using u and v, and one where they are advected using uh and vh. When u and v are used (i.e. the drifters are advected using the resolved velocities only, with no sub-grd component), the particles package needs the layer thickness that is consistent with these velocities, so particles_run needs to be called before any subgrid scale parameterizations are applied. This was not implemented correctly in my last PR, and is corrected here. The particles package also needs the correct timestep for particle advection. This is added here. --- .../external/drifters/MOM_particles.F90 | 4 +-- src/core/MOM.F90 | 26 +++++++++++++------ 2 files changed, 20 insertions(+), 10 deletions(-) diff --git a/config_src/external/drifters/MOM_particles.F90 b/config_src/external/drifters/MOM_particles.F90 index 95470e6510..1c41170582 100644 --- a/config_src/external/drifters/MOM_particles.F90 +++ b/config_src/external/drifters/MOM_particles.F90 @@ -28,7 +28,7 @@ subroutine particles_init(parts, Grid, Time, dt, u, v, h) end subroutine particles_init !> The main driver the steps updates particles -subroutine particles_run(parts, time, uo, vo, ho, tv, use_uh, stagger) +subroutine particles_run(parts, time, uo, vo, ho, tv, dt_adv, use_uh) ! Arguments type(particles), pointer :: parts !< Container for all types and memory type(time_type), intent(in) :: time !< Model time @@ -40,8 +40,8 @@ subroutine particles_run(parts, time, uo, vo, ho, tv, use_uh, stagger) !! that are used to advect tracers [H L2 ~> m3 or kg] real, dimension(:,:,:), intent(in) :: ho !< Ocean layer thickness [H ~> m or kg m-2] type(thermo_var_ptrs), intent(in) :: tv !< structure containing pointers to available thermodynamic fields + real, intent(in) :: dt_adv !< timestep for advecting particles [s] logical :: use_uh !< Flag for whether u and v are weighted by thickness - integer, optional, intent(in) :: stagger !< Flag for whether velocities are staggered end subroutine particles_run diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index e5d1f46086..8cfab91cfc 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1283,6 +1283,18 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & endif ! -------------------------------------------------- end SPLIT + if (CS%use_particles .and. CS%do_dynamics .and. (.not. CS%use_uh_particles)) then + if (CS%thickness_diffuse_first) call MOM_error(WARNING,"particles_run: "//& + "Thickness_diffuse_first is true and use_uh_particles is false. "//& + "This is usually a bad combination.") + !Run particles using unweighted velocity + call particles_run(CS%particles, Time_local, CS%u, CS%v, CS%h, & + CS%tv, dt, CS%use_uh_particles) + call particles_to_z_space(CS%particles, h) + endif + + + ! Update the model's current to reflect wind-wave growth if (Waves%Stokes_DDT .and. (.not.Waves%Passive_Stokes_DDT)) then do J=jsq,jeq ; do i=is,ie @@ -1368,23 +1380,21 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & endif call disable_averaging(CS%diag) + ! Advance the dynamics time by dt. + CS%t_dyn_rel_adv = CS%t_dyn_rel_adv + dt + if (CS%use_particles .and. CS%do_dynamics .and. CS%use_uh_particles) then !Run particles using thickness-weighted velocity call particles_run(CS%particles, Time_local, CS%uhtr, CS%vhtr, CS%h, & - CS%tv, CS%use_uh_particles) - elseif (CS%use_particles .and. CS%do_dynamics) then - !Run particles using unweighted velocity - call particles_run(CS%particles, Time_local, CS%u, CS%v, CS%h, & - CS%tv, CS%use_uh_particles) + CS%tv, CS%t_dyn_rel_adv, CS%use_uh_particles) endif - - ! Advance the dynamics time by dt. - CS%t_dyn_rel_adv = CS%t_dyn_rel_adv + dt CS%n_dyn_steps_in_adv = CS%n_dyn_steps_in_adv + 1 if (CS%alternate_first_direction) then call set_first_direction(G, MODULO(G%first_direction+1,2)) CS%first_dir_restart = real(G%first_direction) + elseif (CS%use_particles .and. CS%do_dynamics .and. (.not.CS%use_uh_particles)) then + call particles_to_k_space(CS%particles, h) endif CS%t_dyn_rel_thermo = CS%t_dyn_rel_thermo + dt if (abs(CS%t_dyn_rel_thermo) < 1e-6*dt) CS%t_dyn_rel_thermo = 0.0