From e48f4a7f42f4e28cc5ac2b2dee5327617eeedec2 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 6 Dec 2021 10:58:13 -0500 Subject: [PATCH 1/3] +Add the new routine unit_no_scaling_init Added the public interface unit_no_scaling_init() to the MOM_unit_scaling module to initialize a non-scaling unit_scale_type without requiring a param_file_type argument. This should be useful for certain trivial unit tests, and it makes it more plausible to make the unit_scale_type arguments mandatory. As a part of this change, the new internal subroutine set_unit_scaling_combos() was carved out of unit_scaling_init to avoid code duplication. Also added comments describing the effective units of the various scaling factors with the standard MOM6 notation. All answers and output are bitwise identical. --- src/framework/MOM_unit_scaling.F90 | 93 ++++++++++++++++++++---------- 1 file changed, 63 insertions(+), 30 deletions(-) diff --git a/src/framework/MOM_unit_scaling.F90 b/src/framework/MOM_unit_scaling.F90 index dbcd2405ec..cd339f410c 100644 --- a/src/framework/MOM_unit_scaling.F90 +++ b/src/framework/MOM_unit_scaling.F90 @@ -8,39 +8,47 @@ module MOM_unit_scaling implicit none ; private -public unit_scaling_init, unit_scaling_end, fix_restart_unit_scaling +! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional +! consistency testing. These are noted in comments with units like Z, H, L, T, R and Q, along with +! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the rescaled +! combination is a nondimensional variable, the notation would be "a slope [Z L-1 ~> nondim]", +! but if (as the case for the variables here), the rescaled combination is exactly 1, the right +! notation would be something like "a dimensional scaling factor [Z m-1 ~> 1]". If the units +! vary with the Boussinesq approximation, the Boussinesq variant is given first. + +public unit_scaling_init, unit_no_scaling_init, unit_scaling_end, fix_restart_unit_scaling !> Describes various unit conversion factors type, public :: unit_scale_type - real :: m_to_Z !< A constant that translates distances in meters to the units of depth. - real :: Z_to_m !< A constant that translates distances in the units of depth to meters. - real :: m_to_L !< A constant that translates lengths in meters to the units of horizontal lengths. - real :: L_to_m !< A constant that translates lengths in the units of horizontal lengths to meters. - real :: s_to_T !< A constant that translates time intervals in seconds to the units of time. - real :: T_to_s !< A constant that translates the units of time to seconds. - real :: R_to_kg_m3 !< A constant that translates the units of density to kilograms per meter cubed. - real :: kg_m3_to_R !< A constant that translates kilograms per meter cubed to the units of density. - real :: Q_to_J_kg !< A constant that translates the units of enthalpy to Joules per kilogram. - real :: J_kg_to_Q !< A constant that translates Joules per kilogram to the units of enthalpy. + real :: m_to_Z !< A constant that translates distances in meters to the units of depth [Z m-1 ~> 1] + real :: Z_to_m !< A constant that translates distances in the units of depth to meters [m Z-1 ~> 1] + real :: m_to_L !< A constant that translates lengths in meters to the units of horizontal lengths [L m-1 ~> 1] + real :: L_to_m !< A constant that translates lengths in the units of horizontal lengths to meters [m L-1 ~> 1] + real :: s_to_T !< A constant that translates time intervals in seconds to the units of time [T s-1 ~> 1] + real :: T_to_s !< A constant that translates the units of time to seconds [s T-1 ~> 1] + real :: R_to_kg_m3 !< A constant that translates the units of density to kilograms per meter cubed [kg m-3 R-1 ~> 1] + real :: kg_m3_to_R !< A constant that translates kilograms per meter cubed to the units of density [R m3 kg-1 ~> 1] + real :: Q_to_J_kg !< A constant that translates the units of enthalpy to Joules per kilogram [J kg-1 Q-1 ~> 1] + real :: J_kg_to_Q !< A constant that translates Joules per kilogram to the units of enthalpy [Q kg J-1 ~> 1] ! These are useful combinations of the fundamental scale conversion factors above. - real :: Z_to_L !< Convert vertical distances to lateral lengths - real :: L_to_Z !< Convert lateral lengths to vertical distances - real :: L_T_to_m_s !< Convert lateral velocities from L T-1 to m s-1. - real :: m_s_to_L_T !< Convert lateral velocities from m s-1 to L T-1. - real :: L_T2_to_m_s2 !< Convert lateral accelerations from L T-2 to m s-2. - real :: Z2_T_to_m2_s !< Convert vertical diffusivities from Z2 T-1 to m2 s-1. - real :: m2_s_to_Z2_T !< Convert vertical diffusivities from m2 s-1 to Z2 T-1. - real :: W_m2_to_QRZ_T !< Convert heat fluxes from W m-2 to Q R Z T-1. - real :: QRZ_T_to_W_m2 !< Convert heat fluxes from Q R Z T-1 to W m-2. - ! Not used enough: real :: kg_m2_to_RZ !< Convert mass loads from kg m-2 to R Z. - real :: RZ_to_kg_m2 !< Convert mass loads from R Z to kg m-2. - real :: kg_m2s_to_RZ_T !< Convert mass fluxes from kg m-2 s-1 to R Z T-1. - real :: RZ_T_to_kg_m2s !< Convert mass fluxes from R Z T-1 to kg m-2 s-1. - real :: RZ3_T3_to_W_m2 !< Convert turbulent kinetic energy fluxes from R Z3 T-3 to W m-2. - real :: W_m2_to_RZ3_T3 !< Convert turbulent kinetic energy fluxes from W m-2 to R Z3 T-3. - real :: RL2_T2_to_Pa !< Convert pressures from R L2 T-2 to Pa. - ! Not used enough: real :: Pa_to_RL2_T2 !< Convert pressures from Pa to R L2 T-2. + real :: Z_to_L !< Convert vertical distances to lateral lengths [L Z-1 ~> 1] + real :: L_to_Z !< Convert lateral lengths to vertical distances [Z L-1 ~> 1] + real :: L_T_to_m_s !< Convert lateral velocities from L T-1 to m s-1 [T m L-1 s-1 ~> 1] + real :: m_s_to_L_T !< Convert lateral velocities from m s-1 to L T-1 [L s T-1 m-1 ~> 1] + real :: L_T2_to_m_s2 !< Convert lateral accelerations from L T-2 to m s-2 [L s2 T-2 m-1 ~> 1] + real :: Z2_T_to_m2_s !< Convert vertical diffusivities from Z2 T-1 to m2 s-1 [T1 m2 Z-2 s-1 ~> 1] + real :: m2_s_to_Z2_T !< Convert vertical diffusivities from m2 s-1 to Z2 T-1 [Z2 s T-1 m-2 ~> 1] + real :: W_m2_to_QRZ_T !< Convert heat fluxes from W m-2 to Q R Z T-1 [Q R Z m2 T-1 W-1 ~> 1] + real :: QRZ_T_to_W_m2 !< Convert heat fluxes from Q R Z T-1 to W m-2 [W T Q-1 R-1 Z-1 m-2 ~> 1] + ! Not used enough: real :: kg_m2_to_RZ !< Convert mass loads from kg m-2 to R Z [R Z m2 kg-1 ~> 1] + real :: RZ_to_kg_m2 !< Convert mass loads from R Z to kg m-2 [kg R-1 Z-1 m-2 ~> 1] + real :: kg_m2s_to_RZ_T !< Convert mass fluxes from kg m-2 s-1 to R Z T-1 [R Z m2 s T-1 kg-1 ~> 1] + real :: RZ_T_to_kg_m2s !< Convert mass fluxes from R Z T-1 to kg m-2 s-1 [T kg R-1 Z-1 m-2 s-1 ~> 1] + real :: RZ3_T3_to_W_m2 !< Convert turbulent kinetic energy fluxes from R Z3 T-3 to W m-2 [W T3 R-1 Z-3 m-2 ~> 1] + real :: W_m2_to_RZ3_T3 !< Convert turbulent kinetic energy fluxes from W m-2 to R Z3 T-3 [R Z3 m2 T-3 W-1 ~> 1] + real :: RL2_T2_to_Pa !< Convert pressures from R L2 T-2 to Pa [Pa T2 R-1 L-2 ~> 1] + ! Not used enough: real :: Pa_to_RL2_T2 !< Convert pressures from Pa to R L2 T-2 [R L2 T-2 Pa-1 ~> 1] ! These are used for changing scaling across restarts. real :: m_to_Z_restart = 0.0 !< A copy of the m_to_Z that is used in restart files. @@ -130,7 +138,32 @@ subroutine unit_scaling_init( param_file, US ) US%Q_to_J_kg = 1.0 * Q_Rescale_factor US%J_kg_to_Q = 1.0 / Q_Rescale_factor - ! These are useful combinations of the fundamental scale conversion factors set above. + call set_unit_scaling_combos(US) +end subroutine unit_scaling_init + +!> Allocates and initializes the ocean model unit scaling type to unscaled values. +subroutine unit_no_scaling_init(US) + type(unit_scale_type), pointer :: US !< A dimensional unit scaling type + + if (associated(US)) call MOM_error(FATAL, & + 'unit_scaling_init: called with an associated US pointer.') + allocate(US) + + US%Z_to_m = 1.0 ; US%m_to_Z = 1.0 + US%L_to_m = 1.0 ; US%m_to_L = 1.0 + US%T_to_s = 1.0 ; US%s_to_T = 1.0 + US%R_to_kg_m3 = 1.0 ; US%kg_m3_to_R = 1.0 + US%Q_to_J_kg = 1.0 ; US%J_kg_to_Q = 1.0 + + call set_unit_scaling_combos(US) +end subroutine unit_no_scaling_init + +!> This subroutine sets useful combinations of the fundamental scale conversion factors +!! in the unit scaling type. +subroutine set_unit_scaling_combos(US) + type(unit_scale_type), intent(inout) :: US !< A dimensional unit scaling type + + ! Convert vertical to horizontal length scales and the reverse: US%Z_to_L = US%Z_to_m * US%m_to_L US%L_to_Z = US%L_to_m * US%m_to_Z ! Horizontal velocities: @@ -159,7 +192,7 @@ subroutine unit_scaling_init( param_file, US ) ! It does not seem like US%Pa_to_RL2_T2 would be used enough in MOM6 to justify its existence. ! US%Pa_to_RL2_T2 = US%kg_m3_to_R * US%m_s_to_L_T**2 -end subroutine unit_scaling_init +end subroutine set_unit_scaling_combos !> Set the unit scaling factors for output to restart files to the unit scaling !! factors for this run. From 59c592649bc404cf933b88b988077f7ecdf9bd65 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 6 Dec 2021 11:50:39 -0500 Subject: [PATCH 2/3] (*)Provide US arguments to 4 existing calls Provide optional US arguments to 4 existing calls to set_grid_metrics() and MOM_initialize_topography(), enabling these arguments to become mandatory in the next commit. In some cases this requires a call to unit_no_scaling_init() or the removal of a call to rescale_dyn_horgrid_bathymetry(). In addition, a mis-scaled value was corrected in the initialization of the ODA control structures private thickness array; this latter issue is not a problem for Boussinesq models without dimensional rescaling (i.e. the tested cases), but could change answers in other cases with data assimilation. All answers in the MOM6-examples or TC test cases are bitwise identical. --- config_src/drivers/unit_drivers/MOM_sum_driver.F90 | 7 ++++++- src/ice_shelf/MOM_ice_shelf.F90 | 6 ++---- src/ocean_data_assim/MOM_oda_driver.F90 | 8 ++++---- 3 files changed, 12 insertions(+), 9 deletions(-) diff --git a/config_src/drivers/unit_drivers/MOM_sum_driver.F90 b/config_src/drivers/unit_drivers/MOM_sum_driver.F90 index 7291eb913a..9f3950ac7f 100644 --- a/config_src/drivers/unit_drivers/MOM_sum_driver.F90 +++ b/config_src/drivers/unit_drivers/MOM_sum_driver.F90 @@ -29,6 +29,7 @@ program MOM_main use MOM_io, only : MOM_io_init, file_exists, open_file, close_file use MOM_io, only : check_nml_error, io_infra_init, io_infra_end use MOM_io, only : APPEND_FILE, ASCII_FILE, READONLY_FILE, SINGLE_FILE + use MOM_unit_scaling, only : unit_scale_type, unit_no_scaling_init, unit_scaling_end implicit none @@ -39,6 +40,8 @@ program MOM_main type(hor_index_type) :: HI ! A hor_index_type for array extents type(param_file_type) :: param_file ! The structure indicating the file(s) ! containing all run-time parameters. + type(unit_scale_type), pointer :: US => NULL() !< A structure containing various unit + ! conversion factors, but in this case all are 1. real :: max_depth ! The maximum ocean depth [m] integer :: verbosity integer :: num_sums @@ -104,7 +107,8 @@ program MOM_main allocate(depth_tot_fastR(num_sums)) ; depth_tot_fastR(:) = 0.0 ! Set up the parameters of the physical grid - call set_grid_metrics(grid, param_file) + call unit_no_scaling_init(US) + call set_grid_metrics(grid, param_file, US) ! Set up the bottom depth, grid%bathyT either analytically or from file call get_param(param_file, "MOM", "MAXIMUM_DEPTH", max_depth, & @@ -162,6 +166,7 @@ program MOM_main enddo call destroy_dyn_horgrid(grid) + call unit_scaling_end(US) call io_infra_end ; call MOM_infra_end contains diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 77166cece0..674b84807d 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -21,7 +21,6 @@ module MOM_ice_shelf use MOM_domains, only : MOM_domains_init, pass_var, pass_vector, clone_MOM_domain use MOM_domains, only : TO_ALL, CGRID_NE, BGRID_NE, CORNER use MOM_dyn_horgrid, only : dyn_horgrid_type, create_dyn_horgrid, destroy_dyn_horgrid -use MOM_dyn_horgrid, only : rescale_dyn_horgrid_bathymetry use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe use MOM_file_parser, only : read_param, get_param, log_param, log_version, param_file_type use MOM_grid, only : MOM_grid_init, ocean_grid_type @@ -1306,9 +1305,8 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, call create_dyn_horgrid(dG, CS%Grid%HI) call clone_MOM_domain(CS%Grid%Domain,dG%Domain) call set_grid_metrics(dG,param_file,CS%US) - ! Set up the bottom depth, G%D either analytically or from file - call MOM_initialize_topography(dG%bathyT, CS%Grid%max_depth, dG, param_file) - call rescale_dyn_horgrid_bathymetry(dG, CS%US%Z_to_m) + ! Set up the bottom depth, dG%bathyT, either analytically or from file + call MOM_initialize_topography(dG%bathyT, CS%Grid%max_depth, dG, param_file, CS%US) call copy_dyngrid_to_MOM_grid(dG, CS%Grid, CS%US) call destroy_dyn_horgrid(dG) ! endif diff --git a/src/ocean_data_assim/MOM_oda_driver.F90 b/src/ocean_data_assim/MOM_oda_driver.F90 index 161cf16115..d5259d760a 100644 --- a/src/ocean_data_assim/MOM_oda_driver.F90 +++ b/src/ocean_data_assim/MOM_oda_driver.F90 @@ -278,8 +278,8 @@ subroutine init_oda(Time, G, GV, diag_CS, CS) allocate(dG) call create_dyn_horgrid(dG, HI) call clone_MOM_domain(CS%Grid%Domain, dG%Domain,symmetric=.false.) - call set_grid_metrics(dG,PF) - call MOM_initialize_topography(dg%bathyT,dG%max_depth,dG,PF) + call set_grid_metrics(dG, PF, CS%US) + call MOM_initialize_topography(dG%bathyT, dG%max_depth, dG, PF, CS%US) call MOM_initialize_coord(CS%GV, CS%US, PF, .false., & dirs%output_directory, tv_dummy, dG%max_depth) call ALE_init(PF, CS%GV, CS%US, dG%max_depth, CS%ALE_CS) @@ -313,9 +313,9 @@ subroutine init_oda(Time, G, GV, diag_CS, CS) !jsd=jsd+jdg_offset; jed=jed+jdg_offset ! TODO: switch to local indexing? (mjh) if (.not. associated(CS%h)) then - allocate(CS%h(isd:ied,jsd:jed,CS%GV%ke), source=CS%GV%Angstrom_m*CS%GV%H_to_m) + allocate(CS%h(isd:ied,jsd:jed,CS%GV%ke), source=CS%GV%Angstrom_H) ! assign thicknesses - call ALE_initThicknessToCoord(CS%ALE_CS,G,CS%GV,CS%h) + call ALE_initThicknessToCoord(CS%ALE_CS, G, CS%GV, CS%h) endif allocate(CS%tv) allocate(CS%tv%T(isd:ied,jsd:jed,CS%GV%ke), source=0.0) From 3162bd08690b518633ab66bd8c37ebe3125616d6 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 6 Dec 2021 11:57:20 -0500 Subject: [PATCH 3/3] +Make US arguments non-optional for 28 routines Made the unit_scale_type arguments non-optional for 28 routines. These arguments had been optional in the first place to manage the coordination between the MOM6 and SIS2 repositories, but SIS2 has been using these optional arguments for several years now, and they can be made mandatory without imposing any disruptions. This change simplifies and clarifies the code. All answers and output are bitwise identical. --- src/core/MOM_checksum_packages.F90 | 8 +- .../MOM_fixed_initialization.F90 | 15 +- src/initialization/MOM_grid_initialize.F90 | 105 +++----- .../MOM_shared_initialization.F90 | 252 +++++++----------- .../MOM_state_initialization.F90 | 2 +- src/user/DOME_initialization.F90 | 15 +- src/user/ISOMIP_initialization.F90 | 15 +- src/user/Kelvin_initialization.F90 | 11 +- src/user/Phillips_initialization.F90 | 10 +- src/user/benchmark_initialization.F90 | 13 +- src/user/shelfwave_initialization.F90 | 15 +- src/user/user_initialization.F90 | 6 +- 12 files changed, 189 insertions(+), 278 deletions(-) diff --git a/src/core/MOM_checksum_packages.F90 b/src/core/MOM_checksum_packages.F90 index be00de8779..917a4afdc3 100644 --- a/src/core/MOM_checksum_packages.F90 +++ b/src/core/MOM_checksum_packages.F90 @@ -92,23 +92,21 @@ subroutine MOM_state_chksum_3arg(mesg, u, v, h, G, GV, US, haloshift, symmetric) intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1] or [m s-1].. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]. - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type, which is + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type, which is !! used to rescale u and v if present. integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0). logical, optional, intent(in) :: symmetric !< If true, do checksums on the fully !! symmetric computational domain. - real :: L_T_to_m_s ! A rescaling factor for velocities [m T s-1 L-1 ~> 1] or [1] + integer :: hs logical :: sym - L_T_to_m_s = 1.0 ; if (present(US)) L_T_to_m_s = US%L_T_to_m_s - ! Note that for the chksum calls to be useful for reproducing across PE ! counts, there must be no redundant points, so all variables use is..ie ! and js...je as their extent. hs = 1 ; if (present(haloshift)) hs = haloshift sym = .false. ; if (present(symmetric)) sym = symmetric - call uvchksum(mesg//" u", u, v, G%HI, haloshift=hs, symmetric=sym, scale=L_T_to_m_s) + call uvchksum(mesg//" u", u, v, G%HI, haloshift=hs, symmetric=sym, scale=US%L_T_to_m_s) call hchksum(h, mesg//" h",G%HI, haloshift=hs, scale=GV%H_to_m) end subroutine MOM_state_chksum_3arg diff --git a/src/initialization/MOM_fixed_initialization.F90 b/src/initialization/MOM_fixed_initialization.F90 index b67d21ebcb..f0fb1d23f9 100644 --- a/src/initialization/MOM_fixed_initialization.F90 +++ b/src/initialization/MOM_fixed_initialization.F90 @@ -6,7 +6,7 @@ module MOM_fixed_initialization use MOM_debugging, only : hchksum, qchksum, uvchksum use MOM_domains, only : pass_var -use MOM_dyn_horgrid, only : dyn_horgrid_type, rescale_dyn_horgrid_bathymetry +use MOM_dyn_horgrid, only : dyn_horgrid_type use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint use MOM_file_parser, only : get_param, read_param, log_param, param_file_type @@ -82,7 +82,6 @@ subroutine MOM_initialize_fixed(G, US, OBC, PF, write_geom, output_dir) ! This also sets G%max_depth based on the input parameter MAXIMUM_DEPTH, ! or, if absent, is diagnosed as G%max_depth = max( G%D(:,:) ) call MOM_initialize_topography(G%bathyT, G%max_depth, G, PF, US) -! call rescale_dyn_horgrid_bathymetry(G, US%Z_to_m) ! To initialize masks, the bathymetry in halo regions must be filled in call pass_var(G%bathyT, G%Domain) @@ -174,20 +173,16 @@ subroutine MOM_initialize_topography(D, max_depth, G, PF, US) intent(out) :: D !< Ocean bottom depth [Z ~> m] or [m] type(param_file_type), intent(in) :: PF !< Parameter file structure real, intent(out) :: max_depth !< Maximum depth of model [Z ~> m] or [m] - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! This subroutine makes the appropriate call to set up the bottom depth. ! This is a separate subroutine so that it can be made public and shared with ! the ice-sheet code or other components. ! Local variables - real :: m_to_Z, Z_to_m ! Dimensional rescaling factors character(len=40) :: mdl = "MOM_initialize_topography" ! This subroutine's name. character(len=200) :: config - m_to_Z = 1.0 ; if (present(US)) m_to_Z = US%m_to_Z - Z_to_m = 1.0 ; if (present(US)) Z_to_m = US%Z_to_m - call get_param(PF, mdl, "TOPO_CONFIG", config, & "This specifies how bathymetry is specified: \n"//& " \t file - read bathymetric information from the file \n"//& @@ -216,7 +211,7 @@ subroutine MOM_initialize_topography(D, max_depth, G, PF, US) " \t dense - Denmark Strait-like dense water formation and overflow.\n"//& " \t USER - call a user modified routine.", & fail_if_missing=.true.) - max_depth = -1.e9*m_to_Z ; call read_param(PF, "MAXIMUM_DEPTH", max_depth, scale=m_to_Z) + max_depth = -1.e9*US%m_to_Z ; call read_param(PF, "MAXIMUM_DEPTH", max_depth, scale=US%m_to_Z) select case ( trim(config) ) case ("file"); call initialize_topography_from_file(D, G, PF, US) case ("flat"); call initialize_topography_named(D, G, PF, config, max_depth, US) @@ -241,11 +236,11 @@ subroutine MOM_initialize_topography(D, max_depth, G, PF, US) "Unrecognized topography setup '"//trim(config)//"'") end select if (max_depth>0.) then - call log_param(PF, mdl, "MAXIMUM_DEPTH", max_depth*Z_to_m, & + call log_param(PF, mdl, "MAXIMUM_DEPTH", max_depth*US%Z_to_m, & "The maximum depth of the ocean.", units="m") else max_depth = diagnoseMaximumDepth(D,G) - call log_param(PF, mdl, "!MAXIMUM_DEPTH", max_depth*Z_to_m, & + call log_param(PF, mdl, "!MAXIMUM_DEPTH", max_depth*US%Z_to_m, & "The (diagnosed) maximum depth of the ocean.", units="m", like_default=.true.) endif if (trim(config) /= "DOME") then diff --git a/src/initialization/MOM_grid_initialize.F90 b/src/initialization/MOM_grid_initialize.F90 index f67a977d27..498e1915ba 100644 --- a/src/initialization/MOM_grid_initialize.F90 +++ b/src/initialization/MOM_grid_initialize.F90 @@ -58,19 +58,14 @@ module MOM_grid_initialize subroutine set_grid_metrics(G, param_file, US) type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type type(param_file_type), intent(in) :: param_file !< Parameter file structure - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! Local variables - real :: m_to_L ! A length unit conversion factor [L m-1 ~> 1] - real :: L_to_m ! A length unit conversion factor [m L-1 ~> 1] ! This include declares and sets the variable "version". # include "version_variable.h" logical :: debug character(len=256) :: config - m_to_L = 1.0 ; if (present(US)) m_to_L = US%m_to_L - L_to_m = 1.0 ; if (present(US)) L_to_m = US%L_to_m - call callTree_enter("set_grid_metrics(), MOM_grid_initialize.F90") call log_version(param_file, "MOM_grid_init", version, "") call get_param(param_file, "MOM_grid_init", "GRID_CONFIG", config, & @@ -88,7 +83,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%Rad_Earth_L = -1.0*m_to_L ; G%len_lat = 0.0 ; G%len_lon = 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) case ("cartesian"); call set_grid_metrics_cartesian(G, param_file, US) @@ -104,11 +99,11 @@ subroutine set_grid_metrics(G, param_file, US) ! The grid metrics were set with an option that does not explicitly initialize Rad_Earth. ! ### Rad_Earth should be read as in: ! call get_param(param_file, mdl, "RAD_EARTH", G%Rad_Earth_L, & - ! "The radius of the Earth.", units="m", default=6.378e6, scale=m_to_L) + ! "The radius of the Earth.", units="m", default=6.378e6, scale=US%m_to_L) ! but for now it is being set via a hard-coded value to reproduce current behavior. - G%Rad_Earth_L = 6.378e6*m_to_L + G%Rad_Earth_L = 6.378e6*US%m_to_L endif - G%Rad_Earth = L_to_m*G%Rad_Earth_L + G%Rad_Earth = US%L_to_m*G%Rad_Earth_L ! Calculate derived metrics (i.e. reciprocals and products) call callTree_enter("set_derived_metrics(), MOM_grid_initialize.F90") @@ -127,39 +122,35 @@ end subroutine set_grid_metrics subroutine grid_metrics_chksum(parent, G, US) character(len=*), intent(in) :: parent !< A string identifying the caller type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - real :: m_to_L ! A unit conversion factor [L m-1 ~> 1] - real :: L_to_m ! A unit conversion factor [m L-1 ~> 1] integer :: halo - m_to_L = 1.0 ; if (present(US)) m_to_L = US%m_to_L - L_to_m = 1.0 ; if (present(US)) L_to_m = US%L_to_m halo = min(G%ied-G%iec, G%jed-G%jec, 1) call hchksum_pair(trim(parent)//': d[xy]T', G%dxT, G%dyT, G%HI, & - haloshift=halo, scale=L_to_m, scalar_pair=.true.) + haloshift=halo, scale=US%L_to_m, scalar_pair=.true.) - call uvchksum(trim(parent)//': dxC[uv]', G%dxCu, G%dyCv, G%HI, haloshift=halo, scale=L_to_m) + call uvchksum(trim(parent)//': dxC[uv]', G%dxCu, G%dyCv, G%HI, haloshift=halo, scale=US%L_to_m) - call uvchksum(trim(parent)//': dxC[uv]', G%dyCu, G%dxCv, G%HI, haloshift=halo, scale=L_to_m) + call uvchksum(trim(parent)//': dxC[uv]', G%dyCu, G%dxCv, G%HI, haloshift=halo, scale=US%L_to_m) - call Bchksum_pair(trim(parent)//': dxB[uv]', G%dxBu, G%dyBu, G%HI, haloshift=halo, scale=L_to_m) + call Bchksum_pair(trim(parent)//': dxB[uv]', G%dxBu, G%dyBu, G%HI, haloshift=halo, scale=US%L_to_m) call hchksum_pair(trim(parent)//': Id[xy]T', G%IdxT, G%IdyT, G%HI, & - haloshift=halo, scale=m_to_L, scalar_pair=.true.) + haloshift=halo, scale=US%m_to_L, scalar_pair=.true.) - call uvchksum(trim(parent)//': Id[xy]C[uv]', G%IdxCu, G%IdyCv, G%HI, haloshift=halo, scale=m_to_L) + call uvchksum(trim(parent)//': Id[xy]C[uv]', G%IdxCu, G%IdyCv, G%HI, haloshift=halo, scale=US%m_to_L) - call uvchksum(trim(parent)//': Id[xy]C[uv]', G%IdyCu, G%IdxCv, G%HI, haloshift=halo, scale=m_to_L) + call uvchksum(trim(parent)//': Id[xy]C[uv]', G%IdyCu, G%IdxCv, G%HI, haloshift=halo, scale=US%m_to_L) - call Bchksum_pair(trim(parent)//': Id[xy]B[uv]', G%IdxBu, G%IdyBu, G%HI, haloshift=halo, scale=m_to_L) + call Bchksum_pair(trim(parent)//': Id[xy]B[uv]', G%IdxBu, G%IdyBu, G%HI, haloshift=halo, scale=US%m_to_L) - call hchksum(G%areaT, trim(parent)//': areaT',G%HI, haloshift=halo, scale=L_to_m**2) - call Bchksum(G%areaBu, trim(parent)//': areaBu',G%HI, haloshift=halo, scale=L_to_m**2) + call hchksum(G%areaT, trim(parent)//': areaT',G%HI, haloshift=halo, scale=US%L_to_m**2) + call Bchksum(G%areaBu, trim(parent)//': areaBu',G%HI, haloshift=halo, scale=US%L_to_m**2) - call hchksum(G%IareaT, trim(parent)//': IareaT',G%HI, haloshift=halo, scale=m_to_L**2) - call Bchksum(G%IareaBu, trim(parent)//': IareaBu',G%HI, haloshift=halo, scale=m_to_L**2) + call hchksum(G%IareaT, trim(parent)//': IareaT',G%HI, haloshift=halo, scale=US%m_to_L**2) + call Bchksum(G%IareaBu, trim(parent)//': IareaBu',G%HI, haloshift=halo, scale=US%m_to_L**2) call hchksum(G%geoLonT,trim(parent)//': geoLonT',G%HI, haloshift=halo) call hchksum(G%geoLatT,trim(parent)//': geoLatT',G%HI, haloshift=halo) @@ -178,8 +169,8 @@ end subroutine grid_metrics_chksum !> Sets the grid metrics from a mosaic file. subroutine set_grid_metrics_from_mosaic(G, param_file, US) type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type - type(param_file_type), intent(in) :: param_file !< Parameter file structure - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + type(param_file_type), intent(in) :: param_file !< Parameter file structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! Local variables real, dimension(G%isd :G%ied ,G%jsd :G%jed ) :: tempH1, tempH2, tempH3, tempH4 real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB) :: tempQ1, tempQ2, tempQ3, tempQ4 @@ -197,7 +188,6 @@ subroutine set_grid_metrics_from_mosaic(G, param_file, US) real, dimension(2*G%isd-2:2*G%ied+1,2*G%jsd-3:2*G%jed+1) :: tmpV real, dimension(2*G%isd-3:2*G%ied+1,2*G%jsd-3:2*G%jed+1) :: tmpZ real, dimension(:,:), allocatable :: tmpGlbl - real :: m_to_L ! A unit conversion factor [L m-1 ~> 1] character(len=200) :: filename, grid_file, inputdir character(len=64) :: mdl = "MOM_grid_init set_grid_metrics_from_mosaic" type(MOM_domain_type), pointer :: SGdom => NULL() ! Supergrid domain @@ -207,7 +197,6 @@ subroutine set_grid_metrics_from_mosaic(G, param_file, US) call callTree_enter("set_grid_metrics_from_mosaic(), MOM_grid_initialize.F90") - m_to_L = 1.0 ; if (present(US)) m_to_L = US%m_to_L call get_param(param_file, mdl, "GRID_FILE", grid_file, & "Name of the file from which to read horizontal grid data.", & fail_if_missing=.true.) @@ -331,16 +320,16 @@ subroutine set_grid_metrics_from_mosaic(G, param_file, US) call pass_var(areaBu, G%Domain, position=CORNER) do i=G%isd,G%ied ; do j=G%jsd,G%jed - G%dxT(i,j) = m_to_L*dxT(i,j) ; G%dyT(i,j) = m_to_L*dyT(i,j) ; G%areaT(i,j) = m_to_L**2*areaT(i,j) + G%dxT(i,j) = US%m_to_L*dxT(i,j) ; G%dyT(i,j) = US%m_to_L*dyT(i,j) ; G%areaT(i,j) = US%m_to_L**2*areaT(i,j) enddo ; enddo do I=G%IsdB,G%IedB ; do j=G%jsd,G%jed - G%dxCu(I,j) = m_to_L*dxCu(I,j) ; G%dyCu(I,j) = m_to_L*dyCu(I,j) + G%dxCu(I,j) = US%m_to_L*dxCu(I,j) ; G%dyCu(I,j) = US%m_to_L*dyCu(I,j) enddo ; enddo do i=G%isd,G%ied ; do J=G%JsdB,G%JedB - G%dxCv(i,J) = m_to_L*dxCv(i,J) ; G%dyCv(i,J) = m_to_L*dyCv(i,J) + G%dxCv(i,J) = US%m_to_L*dxCv(i,J) ; G%dyCv(i,J) = US%m_to_L*dyCv(i,J) enddo ; enddo do I=G%IsdB,G%IedB ; do J=G%JsdB,G%JedB - G%dxBu(I,J) = m_to_L*dxBu(I,J) ; G%dyBu(I,J) = m_to_L*dyBu(I,J) ; G%areaBu(I,J) = m_to_L**2*areaBu(I,J) + G%dxBu(I,J) = US%m_to_L*dxBu(I,J) ; G%dyBu(I,J) = US%m_to_L*dyBu(I,J) ; G%areaBu(I,J) = US%m_to_L**2*areaBu(I,J) enddo ; enddo ! Construct axes for diagnostic output (only necessary because "ferret" uses @@ -395,8 +384,8 @@ end subroutine set_grid_metrics_from_mosaic !! sets of points. subroutine set_grid_metrics_cartesian(G, param_file, US) type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type - type(param_file_type), intent(in) :: param_file !< Parameter file structure - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + type(param_file_type), intent(in) :: param_file !< Parameter file structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! Local variables integer :: i, j, isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, I1off, J1off integer :: niglobal, njglobal @@ -405,7 +394,6 @@ subroutine set_grid_metrics_cartesian(G, param_file, US) real :: dx_everywhere, dy_everywhere ! Grid spacings [L ~> m]. real :: I_dx, I_dy ! Inverse grid spacings [L-1 ~> m-1]. real :: PI - real :: m_to_L ! A unit conversion factor [L m-1 ~> 1] character(len=80) :: units_temp character(len=48) :: mdl = "MOM_grid_init set_grid_metrics_cartesian" @@ -416,7 +404,6 @@ subroutine set_grid_metrics_cartesian(G, param_file, US) call callTree_enter("set_grid_metrics_cartesian(), MOM_grid_initialize.F90") - m_to_L = 1.0 ; if (present(US)) m_to_L = US%m_to_L PI = 4.0*atan(1.0) call get_param(param_file, mdl, "AXIS_UNITS", units_temp, & @@ -438,7 +425,7 @@ subroutine set_grid_metrics_cartesian(G, param_file, US) "The longitudinal or x-direction length of the domain.", & units=units_temp, fail_if_missing=.true.) call get_param(param_file, mdl, "RAD_EARTH", G%Rad_Earth_L, & - "The radius of the Earth.", units="m", default=6.378e6, scale=m_to_L) + "The radius of the Earth.", units="m", default=6.378e6, scale=US%m_to_L) if (units_temp(1:1) == 'k') then G%x_axis_units = "kilometers" ; G%y_axis_units = "kilometers" @@ -476,11 +463,11 @@ subroutine set_grid_metrics_cartesian(G, param_file, US) enddo if (units_temp(1:1) == 'k') then ! Axes are measured in km. - dx_everywhere = 1000.0*m_to_L * G%len_lon / (REAL(niglobal)) - dy_everywhere = 1000.0*m_to_L * G%len_lat / (REAL(njglobal)) + 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. - dx_everywhere = m_to_L*G%len_lon / (REAL(niglobal)) - dy_everywhere = m_to_L*G%len_lat / (REAL(njglobal)) + 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. dx_everywhere = G%Rad_Earth_L * G%len_lon * PI / (180.0 * niglobal) dy_everywhere = G%Rad_Earth_L * G%len_lat * PI / (180.0 * njglobal) @@ -531,8 +518,8 @@ end subroutine set_grid_metrics_cartesian !! sets of points. subroutine set_grid_metrics_spherical(G, param_file, US) type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type - type(param_file_type), intent(in) :: param_file !< Parameter file structure - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + type(param_file_type), intent(in) :: param_file !< Parameter file structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! Local variables real :: PI, PI_180! PI = 3.1415926... as 4*atan(1) integer :: i, j, isd, ied, jsd, jed @@ -541,7 +528,6 @@ subroutine set_grid_metrics_spherical(G, param_file, US) real :: grid_latT(G%jsd:G%jed), grid_latB(G%JsdB:G%JedB) real :: grid_lonT(G%isd:G%ied), grid_lonB(G%IsdB:G%IedB) real :: dLon,dLat,latitude,longitude,dL_di - real :: m_to_L ! A unit conversion factor [L m-1 ~> 1] character(len=48) :: mdl = "MOM_grid_init set_grid_metrics_spherical" is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec @@ -551,7 +537,6 @@ subroutine set_grid_metrics_spherical(G, param_file, US) i_offset = G%idg_offset ; j_offset = G%jdg_offset call callTree_enter("set_grid_metrics_spherical(), MOM_grid_initialize.F90") - m_to_L = 1.0 ; if (present(US)) m_to_L = US%m_to_L ! Calculate the values of the metric terms that might be used ! and save them in arrays. @@ -570,7 +555,7 @@ subroutine set_grid_metrics_spherical(G, param_file, US) "The longitudinal length of the domain.", units="degrees", & fail_if_missing=.true.) call get_param(param_file, mdl, "RAD_EARTH", G%Rad_Earth_L, & - "The radius of the Earth.", units="m", default=6.378e6, scale=m_to_L) + "The radius of the Earth.", units="m", default=6.378e6, scale=US%m_to_L) dLon = G%len_lon/G%Domain%niglobal dLat = G%len_lat/G%Domain%njglobal @@ -670,8 +655,8 @@ end subroutine set_grid_metrics_spherical !! sets of points. subroutine set_grid_metrics_mercator(G, param_file, US) type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type - type(param_file_type), intent(in) :: param_file !< Parameter file structure - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + type(param_file_type), intent(in) :: param_file !< Parameter file structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! Local variables integer :: i, j, isd, ied, jsd, jed integer :: I_off, J_off @@ -691,7 +676,6 @@ subroutine set_grid_metrics_mercator(G, param_file, US) real :: fnRef ! fnRef is the value of Int_dj_dy or ! Int_dj_dy at a latitude or longitude that is real :: jRef, iRef ! being set to be at grid index jRef or iRef. - real :: m_to_L ! A unit conversion factor [L m-1 ~> 1] integer :: itt1, itt2 logical :: debug = .FALSE., simple_area = .true. integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, IsdB, IedB, JsdB, JedB @@ -710,7 +694,6 @@ subroutine set_grid_metrics_mercator(G, param_file, US) call callTree_enter("set_grid_metrics_mercator(), MOM_grid_initialize.F90") - m_to_L = 1.0 ; if (present(US)) m_to_L = US%m_to_L ! Calculate the values of the metric terms that might be used ! and save them in arrays. PI = 4.0*atan(1.0) ; PI_2 = 0.5*PI @@ -728,7 +711,7 @@ subroutine set_grid_metrics_mercator(G, param_file, US) "The longitudinal length of the domain.", units="degrees", & fail_if_missing=.true.) call get_param(param_file, mdl, "RAD_EARTH", GP%Rad_Earth_L, & - "The radius of the Earth.", units="m", default=6.378e6, scale=m_to_L) + "The radius of the Earth.", units="m", default=6.378e6, scale=US%m_to_L) G%south_lat = GP%south_lat ; G%len_lat = GP%len_lat G%west_lon = GP%west_lon ; G%len_lon = GP%len_lon G%Rad_Earth_L = GP%Rad_Earth_L @@ -1210,11 +1193,10 @@ end function Adcroft_reciprocal !! any land or boundary point. For points in the interior, mask2dCu, !! mask2dCv, and mask2dBu are all 1.0. subroutine initialize_masks(G, PF, US) - type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type - type(param_file_type), intent(in) :: PF !< Parameter file structure - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type + type(param_file_type), intent(in) :: PF !< Parameter file structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! Local variables - real :: m_to_Z_scale ! A unit conversion factor from m to Z [Z m-1 ~> 1] real :: Dmask ! The depth for masking in the same units as G%bathyT [Z ~> m]. real :: min_depth ! The minimum ocean depth in the same units as G%bathyT [Z ~> m]. real :: mask_depth ! The depth shallower than which to mask a point as land [Z ~> m]. @@ -1222,22 +1204,21 @@ subroutine initialize_masks(G, PF, US) integer :: i, j call callTree_enter("initialize_masks(), MOM_grid_initialize.F90") - m_to_Z_scale = 1.0 ; if (present(US)) m_to_Z_scale = US%m_to_Z call get_param(PF, mdl, "MINIMUM_DEPTH", min_depth, & "If MASKING_DEPTH is unspecified, then anything shallower than "//& "MINIMUM_DEPTH is assumed to be land and all fluxes are masked out. "//& "If MASKING_DEPTH is specified, then all depths shallower than "//& "MINIMUM_DEPTH but deeper than MASKING_DEPTH are rounded to MINIMUM_DEPTH.", & - units="m", default=0.0, scale=m_to_Z_scale) + units="m", default=0.0, scale=US%m_to_Z) call get_param(PF, mdl, "MASKING_DEPTH", mask_depth, & "The depth below which to mask points as land points, for which all "//& "fluxes are zeroed out. MASKING_DEPTH is ignored if it has the special "//& "default value.", & - units="m", default=-9999.0, scale=m_to_Z_scale) + units="m", default=-9999.0, scale=US%m_to_Z) Dmask = mask_depth - if (mask_depth == -9999.*m_to_Z_scale) Dmask = min_depth + if (mask_depth == -9999.0*US%m_to_Z) Dmask = min_depth G%mask2dCu(:,:) = 0.0 ; G%mask2dCv(:,:) = 0.0 ; G%mask2dBu(:,:) = 0.0 diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90 index 9973b64a21..bb5a84033b 100644 --- a/src/initialization/MOM_shared_initialization.F90 +++ b/src/initialization/MOM_shared_initialization.F90 @@ -57,7 +57,7 @@ subroutine MOM_initialize_rotation(f, G, PF, US) type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), intent(out) :: f !< The Coriolis parameter [T-1 ~> s-1] type(param_file_type), intent(in) :: PF !< Parameter file structure - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! This subroutine makes the appropriate call to set up the Coriolis parameter. ! This is a separate subroutine so that it can be made public and shared with @@ -95,11 +95,8 @@ subroutine MOM_calculate_grad_Coriolis(dF_dx, dF_dy, G, US) type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type ! Local variables integer :: i,j - real :: m_to_L ! A unit conversion factor [L m-1 ~> 1] real :: f1, f2 - m_to_L = 1.0 ; if (present(US)) m_to_L = US%m_to_L - if ((LBOUND(G%CoriolisBu,1) > G%isc-1) .or. & (LBOUND(G%CoriolisBu,2) > G%isc-1)) then ! The gradient of the Coriolis parameter can not be calculated with this grid. @@ -139,19 +136,16 @@ end function diagnoseMaximumDepth subroutine initialize_topography_from_file(D, G, param_file, US) type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type real, dimension(G%isd:G%ied,G%jsd:G%jed), & - intent(out) :: D !< Ocean bottom depth in m or Z if US is present + intent(out) :: D !< Ocean bottom depth [Z ~> m] type(param_file_type), intent(in) :: param_file !< Parameter file structure - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! Local variables - real :: m_to_Z ! A dimensional rescaling factor. character(len=200) :: filename, topo_file, inputdir ! Strings for file/path character(len=200) :: topo_varname ! Variable name in file character(len=40) :: mdl = "initialize_topography_from_file" ! This subroutine's name. call callTree_enter(trim(mdl)//"(), MOM_shared_initialization.F90") - m_to_Z = 1.0 ; if (present(US)) m_to_Z = US%m_to_Z - call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") inputdir = slasher(inputdir) call get_param(param_file, mdl, "TOPO_FILE", topo_file, & @@ -167,13 +161,13 @@ subroutine initialize_topography_from_file(D, G, param_file, US) if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, & " initialize_topography_from_file: Unable to open "//trim(filename)) - D(:,:) = -9.e30*m_to_Z ! Initializing to a very large negative depth (tall mountains) everywhere + D(:,:) = -9.0e30*US%m_to_Z ! Initializing to a very large negative depth (tall mountains) everywhere ! before reading from a file should do nothing. However, in the instance of ! masked-out PEs, halo regions are not updated when a processor does not ! exist. We need to ensure the depth in masked-out PEs appears to be that ! of land so this line does that in the halo regions. For non-masked PEs ! the halo region is filled properly with a later pass_var(). - call MOM_read_data(filename, trim(topo_varname), D, G%Domain, scale=m_to_Z) + call MOM_read_data(filename, trim(topo_varname), D, G%Domain, scale=US%m_to_Z) call apply_topography_edits_from_file(D, G, param_file, US) @@ -187,10 +181,9 @@ subroutine apply_topography_edits_from_file(D, G, param_file, US) intent(inout) :: D !< Ocean bottom depth [m] or [Z ~> m] if !! US is present type(param_file_type), intent(in) :: param_file !< Parameter file structure - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! Local variables - real :: m_to_Z ! A dimensional rescaling factor. real, dimension(:), allocatable :: new_depth ! The new values of the depths [m] integer, dimension(:), allocatable :: ig, jg ! The global indicies of the points to modify character(len=200) :: topo_edits_file, inputdir ! Strings for file/path @@ -202,8 +195,6 @@ subroutine apply_topography_edits_from_file(D, G, param_file, US) call callTree_enter(trim(mdl)//"(), MOM_shared_initialization.F90") - m_to_Z = 1.0 ; if (present(US)) m_to_Z = US%m_to_Z - call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") inputdir = slasher(inputdir) call get_param(param_file, mdl, "TOPO_EDITS_FILE", topo_edits_file, & @@ -217,13 +208,13 @@ subroutine apply_topography_edits_from_file(D, G, param_file, US) "MINIMUM_DEPTH is assumed to be land and all fluxes are masked out. "//& "If MASKING_DEPTH is specified, then all depths shallower than "//& "MINIMUM_DEPTH but deeper than MASKING_DEPTH are rounded to MINIMUM_DEPTH.", & - units="m", default=0.0, scale=m_to_Z) + units="m", default=0.0, scale=US%m_to_Z) call get_param(param_file, mdl, "MASKING_DEPTH", mask_depth, & "The depth below which to mask points as land points, for which all "//& "fluxes are zeroed out. MASKING_DEPTH is ignored if it has the special "//& "default value.", & - units="m", default=-9999.0, scale=m_to_Z) - if (mask_depth == -9999.*m_to_Z) mask_depth = min_depth + units="m", default=-9999.0, scale=US%m_to_Z) + if (mask_depth == -9999.*US%m_to_Z) mask_depth = min_depth if (len_trim(topo_edits_file)==0) return @@ -263,15 +254,15 @@ subroutine apply_topography_edits_from_file(D, G, param_file, US) i = ig(n) - G%isd_global + 2 ! +1 for python indexing and +1 for ig-isd_global+1 j = jg(n) - G%jsd_global + 2 if (i>=G%isc .and. i<=G%iec .and. j>=G%jsc .and. j<=G%jec) then - if (new_depth(n)*m_to_Z /= mask_depth) then + if (new_depth(n)*US%m_to_Z /= mask_depth) then write(stdout,'(a,3i5,f8.2,a,f8.2,2i4)') & - 'Ocean topography edit: ', n, ig(n), jg(n), D(i,j)/m_to_Z, '->', abs(new_depth(n)), i, j - D(i,j) = abs(m_to_Z*new_depth(n)) ! Allows for height-file edits (i.e. converts negatives) + 'Ocean topography edit: ', n, ig(n), jg(n), D(i,j)*US%Z_to_m, '->', abs(new_depth(n)), i, j + D(i,j) = abs(US%m_to_Z*new_depth(n)) ! Allows for height-file edits (i.e. converts negatives) else if (topo_edits_change_mask) then write(stdout,'(a,3i5,f8.2,a,f8.2,2i4)') & - 'Ocean topography edit: ',n,ig(n),jg(n),D(i,j)/m_to_Z,'->',abs(new_depth(n)),i,j - D(i,j) = abs(m_to_Z*new_depth(n)) ! Allows for height-file edits (i.e. converts negatives) + 'Ocean topography edit: ',n,ig(n),jg(n),D(i,j)*US%Z_to_m,'->',abs(new_depth(n)),i,j + D(i,j) = abs(US%m_to_Z*new_depth(n)) ! Allows for height-file edits (i.e. converts negatives) else call MOM_error(FATAL, ' apply_topography_edits_from_file: '//& "A zero depth edit would change the land mask and is not allowed in"//trim(topo_edits_file)) @@ -289,18 +280,16 @@ end subroutine apply_topography_edits_from_file subroutine initialize_topography_named(D, G, param_file, topog_config, max_depth, US) type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type real, dimension(G%isd:G%ied,G%jsd:G%jed), & - intent(out) :: D !< Ocean bottom depth in m or Z if US is present + intent(out) :: D !< Ocean bottom depth [Z ~> m] type(param_file_type), intent(in) :: param_file !< Parameter file structure character(len=*), intent(in) :: topog_config !< The name of an idealized !! topographic configuration - real, intent(in) :: max_depth !< Maximum depth of model in the units of D - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + real, intent(in) :: max_depth !< Maximum depth [Z ~> m] + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! This subroutine places the bottom depth in m into D(:,:), shaped according to the named config. ! Local variables - real :: m_to_Z ! A dimensional rescaling factor [Z m-1 ~> 1] - real :: m_to_L ! A dimensional rescaling factor [L m-1 ~> 1] real :: min_depth ! The minimum depth [Z ~> m]. real :: PI ! 3.1415926... calculated as 4*atan(1) real :: D0 ! A constant to make the maximum basin depth MAXIMUM_DEPTH [Z ~> m] @@ -315,21 +304,18 @@ subroutine initialize_topography_named(D, G, param_file, topog_config, max_depth call MOM_mesg(" MOM_shared_initialization.F90, initialize_topography_named: "//& "TOPO_CONFIG = "//trim(topog_config), 5) - m_to_Z = 1.0 ; if (present(US)) m_to_Z = US%m_to_Z - m_to_L = 1.0 ; if (present(US)) m_to_L = US%m_to_L - call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, & - "The minimum depth of the ocean.", units="m", default=0.0, scale=m_to_Z) + "The minimum depth of the ocean.", units="m", default=0.0, scale=US%m_to_Z) if (max_depth<=0.) call MOM_error(FATAL,"initialize_topography_named: "// & "MAXIMUM_DEPTH has a non-sensical value! Was it set?") if (trim(topog_config) /= "flat") then call get_param(param_file, mdl, "EDGE_DEPTH", Dedge, & "The depth at the edge of one of the named topographies.", & - units="m", default=100.0, scale=m_to_Z) + units="m", default=100.0, scale=US%m_to_Z) call get_param(param_file, mdl, "TOPOG_SLOPE_SCALE", expdecay, & "The exponential decay scale used in defining some of "//& - "the named topographies.", units="m", default=400000.0, scale=m_to_L) + "the named topographies.", units="m", default=400000.0, scale=US%m_to_L) endif @@ -389,13 +375,12 @@ end subroutine initialize_topography_named subroutine limit_topography(D, G, param_file, max_depth, US) type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type real, dimension(G%isd:G%ied,G%jsd:G%jed), & - intent(inout) :: D !< Ocean bottom depth in m or Z if US is present + intent(inout) :: D !< Ocean bottom depth [Z ~> m] type(param_file_type), intent(in) :: param_file !< Parameter file structure - real, intent(in) :: max_depth !< Maximum depth of model in the units of D - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + real, intent(in) :: max_depth !< Maximum depth of model [Z ~> m] + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! Local variables - real :: m_to_Z ! A dimensional rescaling factor. integer :: i, j character(len=40) :: mdl = "limit_topography" ! This subroutine's name. real :: min_depth ! The shallowest value of wet points [Z ~> m] @@ -403,24 +388,22 @@ subroutine limit_topography(D, G, param_file, max_depth, US) call callTree_enter(trim(mdl)//"(), MOM_shared_initialization.F90") - m_to_Z = 1.0 ; if (present(US)) m_to_Z = US%m_to_Z - call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, & "If MASKING_DEPTH is unspecified, then anything shallower than "//& "MINIMUM_DEPTH is assumed to be land and all fluxes are masked out. "//& "If MASKING_DEPTH is specified, then all depths shallower than "//& "MINIMUM_DEPTH but deeper than MASKING_DEPTH are rounded to MINIMUM_DEPTH.", & - units="m", default=0.0, scale=m_to_Z) + units="m", default=0.0, scale=US%m_to_Z) call get_param(param_file, mdl, "MASKING_DEPTH", mask_depth, & "The depth below which to mask points as land points, for which all "//& "fluxes are zeroed out. MASKING_DEPTH is ignored if it has the special "//& "default value.", & - units="m", default=-9999.0, scale=m_to_Z, do_not_log=.true.) + units="m", default=-9999.0, scale=US%m_to_Z, do_not_log=.true.) ! Make sure that min_depth < D(x,y) < max_depth for ocean points ! TBD: The following f.p. equivalence uses a special value. Originally, any negative value ! indicated the branch. We should create a logical flag to indicate this branch. - if (mask_depth == -9999.*m_to_Z) then + if (mask_depth == -9999.*US%m_to_Z) then if (min_depth<0.) then call MOM_error(FATAL, trim(mdl)//": MINIMUM_DEPTH<0 does not work as expected "//& "unless MASKING_DEPTH has been set appropriately. Set a meaningful "//& @@ -460,22 +443,19 @@ subroutine set_rotation_planetary(f, G, param_file, US) real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), & intent(out) :: f !< Coriolis parameter (vertical component) [T-1 ~> s-1] type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! This subroutine sets up the Coriolis parameter for a sphere character(len=30) :: mdl = "set_rotation_planetary" ! This subroutine's name. integer :: I, J real :: PI real :: omega ! The planetary rotation rate [T-1 ~> s-1] - real :: T_to_s ! A time unit conversion factor call callTree_enter(trim(mdl)//"(), MOM_shared_initialization.F90") - T_to_s = 1.0 ; if (present(US)) T_to_s = US%T_to_s - call get_param(param_file, "set_rotation_planetary", "OMEGA", omega, & "The rotation rate of the earth.", units="s-1", & - default=7.2921e-5, scale=T_to_s) + default=7.2921e-5, scale=US%T_to_s) PI = 4.0*atan(1.0) do I=G%IsdB,G%IedB ; do J=G%JsdB,G%JedB @@ -493,7 +473,7 @@ subroutine set_rotation_beta_plane(f, G, param_file, US) real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), & intent(out) :: f !< Coriolis parameter (vertical component) [T-1 ~> s-1] type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! This subroutine sets up the Coriolis parameter for a beta-plane integer :: I, J @@ -502,9 +482,6 @@ subroutine set_rotation_beta_plane(f, G, param_file, US) real :: beta_lat_ref ! The reference latitude for the beta plane [degrees/km/m/cm] 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 :: T_to_s ! A time unit conversion factor [s T-1 ~> 1] - real :: m_to_L ! A length unit conversion factor [L m-1 ~> 1] - real :: L_to_m ! A length unit conversion factor [m L-1 ~> 1] real :: PI character(len=40) :: mdl = "set_rotation_beta_plane" ! This subroutine's name. character(len=200) :: axis_units @@ -512,31 +489,27 @@ subroutine set_rotation_beta_plane(f, G, param_file, US) call callTree_enter(trim(mdl)//"(), MOM_shared_initialization.F90") - T_to_s = 1.0 ; if (present(US)) T_to_s = US%T_to_s - m_to_L = 1.0 ; if (present(US)) m_to_L = US%m_to_L - L_to_m = 1.0 ; if (present(US)) L_to_m = US%L_to_m - call get_param(param_file, mdl, "F_0", f_0, & "The reference value of the Coriolis parameter with the "//& - "betaplane option.", units="s-1", default=0.0, scale=T_to_s) + "betaplane option.", units="s-1", default=0.0, scale=US%T_to_s) call get_param(param_file, mdl, "BETA", beta, & "The northward gradient of the Coriolis parameter with "//& - "the betaplane option.", units="m-1 s-1", default=0.0, scale=T_to_s*L_to_m) + "the betaplane option.", units="m-1 s-1", default=0.0, scale=US%T_to_s*US%L_to_m) call get_param(param_file, mdl, "AXIS_UNITS", axis_units, default="degrees") PI = 4.0*atan(1.0) 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=m_to_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.E3 * m_to_L + y_scl = 1.0e3 * US%m_to_L case ("m") beta_lat_ref_units = "meters" - y_scl = 1. * m_to_L + 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 @@ -633,90 +606,89 @@ subroutine reset_face_lengths_named(G, param_file, name, US) type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters character(len=*), intent(in) :: name !< The name for the set of face lengths. Only "global_1deg" !! is currently implemented. - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! Local variables character(len=256) :: mesg ! Message for error messages. - real :: m_to_L ! A unit conversion factor [L m-1 ~> 1] - real :: L_to_m ! A unit conversion factor [m L-1 ~> 1] - real :: dx_2 = -1.0, dy_2 = -1.0 + real :: dx_2 ! Half the local zonal grid spacing [degreesE] + real :: dy_2 ! Half the local meridional grid spacing [degreesN] real :: pi_180 - integer :: option = -1 + integer :: option integer :: i, j, isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB pi_180 = (4.0*atan(1.0))/180.0 + dx_2 = -1.0 ; dy_2 = -1.0 + option = -1 + select case ( trim(name) ) case ("global_1deg") ; option = 1 ; dx_2 = 0.5*1.0 case default ; call MOM_error(FATAL, "reset_face_lengths_named: "//& "Unrecognized channel configuration name "//trim(name)) end select - m_to_L = 1.0 ; if (present(US)) m_to_L = US%m_to_L - L_to_m = 1.0 ; if (present(US)) L_to_m = US%L_to_m - if (option==1) then ! 1-degree settings. do j=jsd,jed ; do I=IsdB,IedB ! Change any u-face lengths within this loop. dy_2 = dx_2 * G%dyCu(I,j)*G%IdxCu(I,j) * cos(pi_180 * G%geoLatCu(I,j)) if ((abs(G%geoLatCu(I,j)-35.5) < dy_2) .and. (G%geoLonCu(I,j) < -4.5) .and. & (G%geoLonCu(I,j) > -6.5)) & - G%dy_Cu(I,j) = G%mask2dCu(I,j)*12000.0*m_to_L ! Gibraltar + G%dy_Cu(I,j) = G%mask2dCu(I,j)*12000.0*US%m_to_L ! Gibraltar if ((abs(G%geoLatCu(I,j)-12.5) < dy_2) .and. (abs(G%geoLonCu(I,j)-43.0) < dx_2)) & - G%dy_Cu(I,j) = G%mask2dCu(I,j)*10000.0*m_to_L ! Red Sea + G%dy_Cu(I,j) = G%mask2dCu(I,j)*10000.0*US%m_to_L ! Red Sea if ((abs(G%geoLatCu(I,j)-40.5) < dy_2) .and. (abs(G%geoLonCu(I,j)-26.0) < dx_2)) & - G%dy_Cu(I,j) = G%mask2dCu(I,j)*5000.0*m_to_L ! Dardanelles + G%dy_Cu(I,j) = G%mask2dCu(I,j)*5000.0*US%m_to_L ! Dardanelles if ((abs(G%geoLatCu(I,j)-41.5) < dy_2) .and. (abs(G%geoLonCu(I,j)+220.0) < dx_2)) & - G%dy_Cu(I,j) = G%mask2dCu(I,j)*35000.0*m_to_L ! Tsugaru strait at 140.0e + G%dy_Cu(I,j) = G%mask2dCu(I,j)*35000.0*US%m_to_L ! Tsugaru strait at 140.0e if ((abs(G%geoLatCu(I,j)-45.5) < dy_2) .and. (abs(G%geoLonCu(I,j)+217.5) < 0.9)) & - G%dy_Cu(I,j) = G%mask2dCu(I,j)*15000.0*m_to_L ! Betw Hokkaido and Sakhalin at 217&218 = 142e + G%dy_Cu(I,j) = G%mask2dCu(I,j)*15000.0*US%m_to_L ! Betw Hokkaido and Sakhalin at 217&218 = 142e ! Greater care needs to be taken in the tripolar region. if ((abs(G%geoLatCu(I,j)-80.84) < 0.2) .and. (abs(G%geoLonCu(I,j)+64.9) < 0.8)) & - G%dy_Cu(I,j) = G%mask2dCu(I,j)*38000.0*m_to_L ! Smith Sound in Canadian Arch - tripolar region + G%dy_Cu(I,j) = G%mask2dCu(I,j)*38000.0*US%m_to_L ! Smith Sound in Canadian Arch - tripolar region enddo ; enddo do J=JsdB,JedB ; do i=isd,ied ! Change any v-face lengths within this loop. dy_2 = dx_2 * G%dyCv(i,J)*G%IdxCv(i,J) * cos(pi_180 * G%geoLatCv(i,J)) if ((abs(G%geoLatCv(i,J)-41.0) < dy_2) .and. (abs(G%geoLonCv(i,J)-28.5) < dx_2)) & - G%dx_Cv(i,J) = G%mask2dCv(i,J)*2500.0*m_to_L ! Bosporus - should be 1000.0 m wide. + G%dx_Cv(i,J) = G%mask2dCv(i,J)*2500.0*US%m_to_L ! Bosporus - should be 1000.0 m wide. if ((abs(G%geoLatCv(i,J)-13.0) < dy_2) .and. (abs(G%geoLonCv(i,J)-42.5) < dx_2)) & - G%dx_Cv(i,J) = G%mask2dCv(i,J)*10000.0*m_to_L ! Red Sea + G%dx_Cv(i,J) = G%mask2dCv(i,J)*10000.0*US%m_to_L ! Red Sea if ((abs(G%geoLatCv(i,J)+2.8) < 0.8) .and. (abs(G%geoLonCv(i,J)+241.5) < dx_2)) & - G%dx_Cv(i,J) = G%mask2dCv(i,J)*40000.0*m_to_L ! Makassar Straits at 241.5 W = 118.5 E + G%dx_Cv(i,J) = G%mask2dCv(i,J)*40000.0*US%m_to_L ! Makassar Straits at 241.5 W = 118.5 E if ((abs(G%geoLatCv(i,J)-0.56) < 0.5) .and. (abs(G%geoLonCv(i,J)+240.5) < dx_2)) & - G%dx_Cv(i,J) = G%mask2dCv(i,J)*80000.0*m_to_L ! entry to Makassar Straits at 240.5 W = 119.5 E + G%dx_Cv(i,J) = G%mask2dCv(i,J)*80000.0*US%m_to_L ! entry to Makassar Straits at 240.5 W = 119.5 E if ((abs(G%geoLatCv(i,J)-0.19) < 0.5) .and. (abs(G%geoLonCv(i,J)+230.5) < dx_2)) & - G%dx_Cv(i,J) = G%mask2dCv(i,J)*25000.0*m_to_L ! Channel betw N Guinea and Halmahara 230.5 W = 129.5 E + G%dx_Cv(i,J) = G%mask2dCv(i,J)*25000.0*US%m_to_L ! Channel betw N Guinea and Halmahara 230.5 W = 129.5 E if ((abs(G%geoLatCv(i,J)-0.19) < 0.5) .and. (abs(G%geoLonCv(i,J)+229.5) < dx_2)) & - G%dx_Cv(i,J) = G%mask2dCv(i,J)*25000.0*m_to_L ! Channel betw N Guinea and Halmahara 229.5 W = 130.5 E + G%dx_Cv(i,J) = G%mask2dCv(i,J)*25000.0*US%m_to_L ! Channel betw N Guinea and Halmahara 229.5 W = 130.5 E if ((abs(G%geoLatCv(i,J)-0.0) < 0.25) .and. (abs(G%geoLonCv(i,J)+228.5) < dx_2)) & - G%dx_Cv(i,J) = G%mask2dCv(i,J)*25000.0*m_to_L ! Channel betw N Guinea and Halmahara 228.5 W = 131.5 E + G%dx_Cv(i,J) = G%mask2dCv(i,J)*25000.0*US%m_to_L ! Channel betw N Guinea and Halmahara 228.5 W = 131.5 E if ((abs(G%geoLatCv(i,J)+8.5) < 0.5) .and. (abs(G%geoLonCv(i,J)+244.5) < dx_2)) & - G%dx_Cv(i,J) = G%mask2dCv(i,J)*20000.0*m_to_L ! Lombok Straits at 244.5 W = 115.5 E + G%dx_Cv(i,J) = G%mask2dCv(i,J)*20000.0*US%m_to_L ! Lombok Straits at 244.5 W = 115.5 E if ((abs(G%geoLatCv(i,J)+8.5) < 0.5) .and. (abs(G%geoLonCv(i,J)+235.5) < dx_2)) & - G%dx_Cv(i,J) = G%mask2dCv(i,J)*20000.0*m_to_L ! Timor Straits at 235.5 W = 124.5 E + G%dx_Cv(i,J) = G%mask2dCv(i,J)*20000.0*US%m_to_L ! Timor Straits at 235.5 W = 124.5 E if ((abs(G%geoLatCv(i,J)-52.5) < dy_2) .and. (abs(G%geoLonCv(i,J)+218.5) < dx_2)) & - G%dx_Cv(i,J) = G%mask2dCv(i,J)*2500.0*m_to_L ! Russia and Sakhalin Straits at 218.5 W = 141.5 E + G%dx_Cv(i,J) = G%mask2dCv(i,J)*2500.0*US%m_to_L ! Russia and Sakhalin Straits at 218.5 W = 141.5 E ! Greater care needs to be taken in the tripolar region. if ((abs(G%geoLatCv(i,J)-76.8) < 0.06) .and. (abs(G%geoLonCv(i,J)+88.7) < dx_2)) & - G%dx_Cv(i,J) = G%mask2dCv(i,J)*8400.0*m_to_L ! Jones Sound in Canadian Arch - tripolar region + G%dx_Cv(i,J) = G%mask2dCv(i,J)*8400.0*US%m_to_L ! Jones Sound in Canadian Arch - tripolar region enddo ; enddo endif @@ -724,10 +696,10 @@ subroutine reset_face_lengths_named(G, param_file, name, US) ! These checks apply regardless of the chosen option. do j=jsd,jed ; do I=IsdB,IedB - if (L_to_m*G%dy_Cu(I,j) > L_to_m*G%dyCu(I,j)) then + if (G%dy_Cu(I,j) > G%dyCu(I,j)) then write(mesg,'("dy_Cu of ",ES11.4," exceeds unrestricted width of ",ES11.4,& &" by ",ES11.4," at lon/lat of ", ES11.4, ES11.4)') & - L_to_m*G%dy_Cu(I,j), L_to_m*G%dyCu(I,j), L_to_m*G%dy_Cu(I,j)-L_to_m*G%dyCu(I,j), & + US%L_to_m*G%dy_Cu(I,j), US%L_to_m*G%dyCu(I,j), US%L_to_m*(G%dy_Cu(I,j)-G%dyCu(I,j)), & G%geoLonCu(I,j), G%geoLatCu(I,j) call MOM_error(FATAL,"reset_face_lengths_named "//mesg) endif @@ -737,10 +709,10 @@ subroutine reset_face_lengths_named(G, param_file, name, US) enddo ; enddo do J=JsdB,JedB ; do i=isd,ied - if (L_to_m*G%dx_Cv(i,J) > L_to_m*G%dxCv(i,J)) then + if (G%dx_Cv(i,J) > G%dxCv(i,J)) then write(mesg,'("dx_Cv of ",ES11.4," exceeds unrestricted width of ",ES11.4,& &" by ",ES11.4, " at lon/lat of ", ES11.4, ES11.4)') & - L_to_m*G%dx_Cv(i,J), L_to_m*G%dxCv(i,J), L_to_m*G%dx_Cv(i,J)-L_to_m*G%dxCv(i,J), & + US%L_to_m*G%dx_Cv(i,J), US%L_to_m*G%dxCv(i,J), US%L_to_m*(G%dx_Cv(i,J)-G%dxCv(i,J)), & G%geoLonCv(i,J), G%geoLatCv(i,J) call MOM_error(FATAL,"reset_face_lengths_named "//mesg) @@ -759,22 +731,18 @@ end subroutine reset_face_lengths_named subroutine reset_face_lengths_file(G, param_file, US) type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! Local variables character(len=40) :: mdl = "reset_face_lengths_file" ! This subroutine's name. character(len=256) :: mesg ! Message for error messages. character(len=200) :: filename, chan_file, inputdir ! Strings for file/path - real :: m_to_L ! A unit conversion factor [L m-1 ~> 1] - real :: L_to_m ! A unit conversion factor [m L-1 ~> 1] integer :: i, j, isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB ! These checks apply regardless of the chosen option. call callTree_enter(trim(mdl)//"(), MOM_shared_initialization.F90") - m_to_L = 1.0 ; if (present(US)) m_to_L = US%m_to_L - L_to_m = 1.0 ; if (present(US)) L_to_m = US%L_to_m call get_param(param_file, mdl, "CHANNEL_WIDTH_FILE", chan_file, & "The file from which the list of narrowed channels is read.", & @@ -789,14 +757,14 @@ subroutine reset_face_lengths_file(G, param_file, US) trim(filename)) endif - call MOM_read_vector(filename, "dyCuo", "dxCvo", G%dy_Cu, G%dx_Cv, G%Domain, scale=m_to_L) + call MOM_read_vector(filename, "dyCuo", "dxCvo", G%dy_Cu, G%dx_Cv, G%Domain, scale=US%m_to_L) call pass_vector(G%dy_Cu, G%dx_Cv, G%Domain, To_All+SCALAR_PAIR, CGRID_NE) do j=jsd,jed ; do I=IsdB,IedB - if (L_to_m*G%dy_Cu(I,j) > L_to_m*G%dyCu(I,j)) then + if (G%dy_Cu(I,j) > G%dyCu(I,j)) then write(mesg,'("dy_Cu of ",ES11.4," exceeds unrestricted width of ",ES11.4,& &" by ",ES11.4," at lon/lat of ", ES11.4, ES11.4)') & - L_to_m*G%dy_Cu(I,j), L_to_m*G%dyCu(I,j), L_to_m*G%dy_Cu(I,j)-L_to_m*G%dyCu(I,j), & + US%L_to_m*G%dy_Cu(I,j), US%L_to_m*G%dyCu(I,j), US%L_to_m*(G%dy_Cu(I,j)-G%dyCu(I,j)), & G%geoLonCu(I,j), G%geoLatCu(I,j) call MOM_error(FATAL,"reset_face_lengths_file "//mesg) endif @@ -806,10 +774,10 @@ subroutine reset_face_lengths_file(G, param_file, US) enddo ; enddo do J=JsdB,JedB ; do i=isd,ied - if (L_to_m*G%dx_Cv(i,J) > L_to_m*G%dxCv(i,J)) then + if (G%dx_Cv(i,J) > G%dxCv(i,J)) then write(mesg,'("dx_Cv of ",ES11.4," exceeds unrestricted width of ",ES11.4,& &" by ",ES11.4, " at lon/lat of ", ES11.4, ES11.4)') & - L_to_m*G%dx_Cv(i,J), L_to_m*G%dxCv(i,J), L_to_m*G%dx_Cv(i,J)-L_to_m*G%dxCv(i,J), & + US%L_to_m*G%dx_Cv(i,J), US%L_to_m*G%dxCv(i,J), US%L_to_m*(G%dx_Cv(i,J)-G%dxCv(i,J)), & G%geoLonCv(i,J), G%geoLatCv(i,J) call MOM_error(FATAL,"reset_face_lengths_file "//mesg) @@ -829,7 +797,7 @@ end subroutine reset_face_lengths_file subroutine reset_face_lengths_list(G, param_file, US) type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! Local variables character(len=120), pointer, dimension(:) :: lines => NULL() @@ -847,9 +815,6 @@ subroutine reset_face_lengths_list(G, param_file, US) Dmin_u, Dmax_u, Davg_u ! Porous barrier monomial fit params [m] real, allocatable, dimension(:) :: & Dmin_v, Dmax_v, Davg_v - real :: m_to_L ! A unit conversion factor [L m-1 ~> 1] - real :: L_to_m ! A unit conversion factor [m L-1 ~> 1] - real :: m_to_Z ! A unit conversion factor [Z m-1 ~> 1] real :: lat, lon ! The latitude and longitude of a point. real :: len_lon ! The periodic range of longitudes, usually 360 degrees. real :: len_lat ! The range of latitudes, usually 180 degrees. @@ -870,9 +835,6 @@ subroutine reset_face_lengths_list(G, param_file, US) IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB call callTree_enter(trim(mdl)//"(), MOM_shared_initialization.F90") - m_to_L = 1.0 ; if (present(US)) m_to_L = US%m_to_L - L_to_m = 1.0 ; if (present(US)) L_to_m = US%L_to_m - m_to_Z = 1.0 ; if (present(US)) m_to_Z = US%m_to_Z call get_param(param_file, mdl, "CHANNEL_LIST_FILE", chan_file, & "The file from which the list of narrowed channels is read.", & @@ -1053,10 +1015,10 @@ subroutine reset_face_lengths_list(G, param_file, US) ((lon_p >= u_lon(1,npt)) .and. (lon_p <= u_lon(2,npt))) .or. & ((lon_m >= u_lon(1,npt)) .and. (lon_m <= u_lon(2,npt)))) ) then - G%dy_Cu(I,j) = G%mask2dCu(I,j) * min(G%dyCu(I,j), max(m_to_L*u_width(npt), 0.0)) - G%porous_DminU(I,j) = m_to_Z*Dmin_u(npt) - G%porous_DmaxU(I,j) = m_to_Z*Dmax_u(npt) - G%porous_DavgU(I,j) = m_to_Z*Davg_u(npt) + G%dy_Cu(I,j) = G%mask2dCu(I,j) * min(G%dyCu(I,j), max(US%m_to_L*u_width(npt), 0.0)) + G%porous_DminU(I,j) = US%m_to_Z*Dmin_u(npt) + G%porous_DmaxU(I,j) = US%m_to_Z*Dmax_u(npt) + G%porous_DavgU(I,j) = US%m_to_Z*Davg_u(npt) if (j>=G%jsc .and. j<=G%jec .and. I>=G%isc .and. I<=G%iec) then ! Limit messages/checking to compute domain if ( G%mask2dCu(I,j) == 0.0 ) then @@ -1066,7 +1028,7 @@ subroutine reset_face_lengths_list(G, param_file, US) u_line_used(npt) = u_line_used(npt) + 1 write(stdout,'(A,2F8.2,A,4F8.2,A5,F9.2,A1)') & "read_face_lengths_list : Modifying dy_Cu gridpoint at ",lat,lon," (",& - u_lat(1,npt), u_lat(2,npt), u_lon(1,npt), u_lon(2,npt),") to ",L_to_m*G%dy_Cu(I,j),"m" + u_lat(1,npt), u_lat(2,npt), u_lon(1,npt), u_lon(2,npt),") to ",US%L_to_m*G%dy_Cu(I,j),"m" write(stdout,'(A,3F8.2,A)') & "read_face_lengths_list : Porous Topography parameters: Dmin, Dmax, Davg (",G%porous_DminU(I,j),& G%porous_DmaxU(I,j), G%porous_DavgU(I,j),")m" @@ -1090,10 +1052,10 @@ subroutine reset_face_lengths_list(G, param_file, US) (((lon >= v_lon(1,npt)) .and. (lon <= v_lon(2,npt))) .or. & ((lon_p >= v_lon(1,npt)) .and. (lon_p <= v_lon(2,npt))) .or. & ((lon_m >= v_lon(1,npt)) .and. (lon_m <= v_lon(2,npt)))) ) then - G%dx_Cv(i,J) = G%mask2dCv(i,J) * min(G%dxCv(i,J), max(m_to_L*v_width(npt), 0.0)) - G%porous_DminV(i,J) = m_to_Z*Dmin_v(npt) - G%porous_DmaxV(i,J) = m_to_Z*Dmax_v(npt) - G%porous_DavgV(i,J) = m_to_Z*Davg_v(npt) + G%dx_Cv(i,J) = G%mask2dCv(i,J) * min(G%dxCv(i,J), max(US%m_to_L*v_width(npt), 0.0)) + G%porous_DminV(i,J) = US%m_to_Z*Dmin_v(npt) + G%porous_DmaxV(i,J) = US%m_to_Z*Dmax_v(npt) + G%porous_DavgV(i,J) = US%m_to_Z*Davg_v(npt) if (i>=G%isc .and. i<=G%iec .and. J>=G%jsc .and. J<=G%jec) then ! Limit messages/checking to compute domain if ( G%mask2dCv(i,J) == 0.0 ) then @@ -1103,7 +1065,7 @@ subroutine reset_face_lengths_list(G, param_file, US) v_line_used(npt) = v_line_used(npt) + 1 write(stdout,'(A,2F8.2,A,4F8.2,A5,F9.2,A1)') & "read_face_lengths_list : Modifying dx_Cv gridpoint at ",lat,lon," (",& - v_lat(1,npt), v_lat(2,npt), v_lon(1,npt), v_lon(2,npt),") to ",L_to_m*G%dx_Cv(I,j),"m" + v_lat(1,npt), v_lat(2,npt), v_lon(1,npt), v_lon(2,npt),") to ",US%L_to_m*G%dx_Cv(I,j),"m" write(stdout,'(A,3F8.2,A)') & "read_face_lengths_list : Porous Topography parameters: Dmin, Dmax, Davg (",G%porous_DminV(i,J),& G%porous_DmaxV(i,J), G%porous_DavgV(i,J),")m" @@ -1247,15 +1209,15 @@ end subroutine set_velocity_depth_min !> Pre-compute global integrals of grid quantities (like masked ocean area) for !! later use in reporting diagnostics subroutine compute_global_grid_integrals(G, US) - type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! Local variables real, dimension(G%isc:G%iec, G%jsc:G%jec) :: tmpForSumming - real :: area_scale ! A scaling factor for area into MKS units + real :: area_scale ! A scaling factor for area into MKS units [m2 L-2 ~> 1] integer :: i,j - area_scale = 1.0 ; if (present(US)) area_scale = US%L_to_m**2 + area_scale = US%L_to_m**2 tmpForSumming(:,:) = 0. G%areaT_global = 0.0 ; G%IareaT_global = 0.0 @@ -1275,13 +1237,13 @@ end subroutine compute_global_grid_integrals ! ----------------------------------------------------------------------------- !> Write out a file describing the topography, Coriolis parameter, grid locations !! and various other fixed fields from the grid. -subroutine write_ocean_geometry_file(G, param_file, directory, geom_file, US) +subroutine write_ocean_geometry_file(G, param_file, directory, US, geom_file) type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type(param_file_type), intent(in) :: param_file !< Parameter file structure character(len=*), intent(in) :: directory !< The directory into which to place the geometry file. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type character(len=*), optional, intent(in) :: geom_file !< If present, the name of the geometry file !! (otherwise the file is "ocean_geometry") - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type ! Local variables. character(len=240) :: filepath ! The full path to the file to write @@ -1290,9 +1252,6 @@ subroutine write_ocean_geometry_file(G, param_file, directory, geom_file, US) vars ! Types with metadata about the variables and their staggering type(fieldtype), dimension(:), allocatable :: & fields ! Opaque types used by MOM_io to store variable metadata information - real :: Z_to_m_scale ! A unit conversion factor from Z to m - real :: s_to_T_scale ! A unit conversion factor from T-1 to s-1 - real :: L_to_m_scale ! A unit conversion factor from L to m type(file_type) :: IO_handle ! The I/O handle of the fileset integer :: nFlds ! The number of variables in this file integer :: file_threading @@ -1300,11 +1259,6 @@ subroutine write_ocean_geometry_file(G, param_file, directory, geom_file, US) call callTree_enter('write_ocean_geometry_file()') - Z_to_m_scale = 1.0 ; if (present(US)) Z_to_m_scale = US%Z_to_m - s_to_T_scale = 1.0 ; if (present(US)) s_to_T_scale = US%s_to_T - L_to_m_scale = 1.0 ; if (present(US)) L_to_m_scale = US%L_to_m - - nFlds = 19 ; if (G%bathymetry_at_vel) nFlds = 23 allocate(vars(nFlds)) @@ -1369,30 +1323,30 @@ subroutine write_ocean_geometry_file(G, param_file, directory, geom_file, US) call MOM_write_field(IO_handle, fields(3), G%Domain, G%geoLatT) call MOM_write_field(IO_handle, fields(4), G%Domain, G%geoLonT) - call MOM_write_field(IO_handle, fields(5), G%Domain, G%bathyT, scale=Z_to_m_scale) - call MOM_write_field(IO_handle, fields(6), G%Domain, G%CoriolisBu, scale=s_to_T_scale) + call MOM_write_field(IO_handle, fields(5), G%Domain, G%bathyT, scale=US%Z_to_m) + call MOM_write_field(IO_handle, fields(6), G%Domain, G%CoriolisBu, scale=US%s_to_T) - call MOM_write_field(IO_handle, fields(7), G%Domain, G%dxCv, scale=L_to_m_scale) - call MOM_write_field(IO_handle, fields(8), G%Domain, G%dyCu, scale=L_to_m_scale) - call MOM_write_field(IO_handle, fields(9), G%Domain, G%dxCu, scale=L_to_m_scale) - call MOM_write_field(IO_handle, fields(10), G%Domain, G%dyCv, scale=L_to_m_scale) - call MOM_write_field(IO_handle, fields(11), G%Domain, G%dxT, scale=L_to_m_scale) - call MOM_write_field(IO_handle, fields(12), G%Domain, G%dyT, scale=L_to_m_scale) - call MOM_write_field(IO_handle, fields(13), G%Domain, G%dxBu, scale=L_to_m_scale) - call MOM_write_field(IO_handle, fields(14), G%Domain, G%dyBu, scale=L_to_m_scale) + call MOM_write_field(IO_handle, fields(7), G%Domain, G%dxCv, scale=US%L_to_m) + call MOM_write_field(IO_handle, fields(8), G%Domain, G%dyCu, scale=US%L_to_m) + call MOM_write_field(IO_handle, fields(9), G%Domain, G%dxCu, scale=US%L_to_m) + call MOM_write_field(IO_handle, fields(10), G%Domain, G%dyCv, scale=US%L_to_m) + call MOM_write_field(IO_handle, fields(11), G%Domain, G%dxT, scale=US%L_to_m) + call MOM_write_field(IO_handle, fields(12), G%Domain, G%dyT, scale=US%L_to_m) + call MOM_write_field(IO_handle, fields(13), G%Domain, G%dxBu, scale=US%L_to_m) + call MOM_write_field(IO_handle, fields(14), G%Domain, G%dyBu, scale=US%L_to_m) - call MOM_write_field(IO_handle, fields(15), G%Domain, G%areaT, scale=L_to_m_scale**2) - call MOM_write_field(IO_handle, fields(16), G%Domain, G%areaBu, scale=L_to_m_scale**2) + call MOM_write_field(IO_handle, fields(15), G%Domain, G%areaT, scale=US%L_to_m**2) + call MOM_write_field(IO_handle, fields(16), G%Domain, G%areaBu, scale=US%L_to_m**2) - call MOM_write_field(IO_handle, fields(17), G%Domain, G%dx_Cv, scale=L_to_m_scale) - call MOM_write_field(IO_handle, fields(18), G%Domain, G%dy_Cu, scale=L_to_m_scale) + call MOM_write_field(IO_handle, fields(17), G%Domain, G%dx_Cv, scale=US%L_to_m) + call MOM_write_field(IO_handle, fields(18), G%Domain, G%dy_Cu, scale=US%L_to_m) call MOM_write_field(IO_handle, fields(19), G%Domain, G%mask2dT) if (G%bathymetry_at_vel) then - call MOM_write_field(IO_handle, fields(20), G%Domain, G%Dblock_u, scale=Z_to_m_scale) - call MOM_write_field(IO_handle, fields(21), G%Domain, G%Dopen_u, scale=Z_to_m_scale) - call MOM_write_field(IO_handle, fields(22), G%Domain, G%Dblock_v, scale=Z_to_m_scale) - call MOM_write_field(IO_handle, fields(23), G%Domain, G%Dopen_v, scale=Z_to_m_scale) + call MOM_write_field(IO_handle, fields(20), G%Domain, G%Dblock_u, scale=US%Z_to_m) + call MOM_write_field(IO_handle, fields(21), G%Domain, G%Dopen_u, scale=US%Z_to_m) + call MOM_write_field(IO_handle, fields(22), G%Domain, G%Dblock_v, scale=US%Z_to_m) + call MOM_write_field(IO_handle, fields(23), G%Domain, G%Dopen_v, scale=US%Z_to_m) endif call close_file(IO_handle) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index ce1f9ad92f..2aab378b4a 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -742,7 +742,7 @@ end subroutine initialize_thickness_from_file !> Adjust interface heights to fit the bathymetry and diagnose layer thickness. !! !! If the bottom most interface is below the topography then the bottom-most -!! layers are contracted to GV%Angstrom_m. +!! layers are contracted to ANGSTROM thickness (which may be 0). !! If the bottom most interface is above the topography then the entire column !! is dilated (expanded) to fill the void. !! @remark{There is a (hard-wired) "tolerance" parameter such that the diff --git a/src/user/DOME_initialization.F90 b/src/user/DOME_initialization.F90 index 23ef41be94..248bf6c0f0 100644 --- a/src/user/DOME_initialization.F90 +++ b/src/user/DOME_initialization.F90 @@ -40,13 +40,12 @@ module DOME_initialization subroutine DOME_initialize_topography(D, G, param_file, max_depth, US) type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type real, dimension(G%isd:G%ied,G%jsd:G%jed), & - intent(out) :: D !< Ocean bottom depth in [m] or [Z ~> m] if US is present + intent(out) :: D !< Ocean bottom depth [Z ~> m] type(param_file_type), intent(in) :: param_file !< Parameter file structure - real, intent(in) :: max_depth !< Maximum model depth [m] or [Z ~> m] - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + real, intent(in) :: max_depth !< Maximum model depth [Z ~> m] + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! Local variables - real :: m_to_Z ! A dimensional rescaling factor. real :: min_depth ! The minimum and maximum depths [Z ~> m]. ! This include declares and sets the variable "version". # include "version_variable.h" @@ -57,22 +56,20 @@ subroutine DOME_initialize_topography(D, G, param_file, max_depth, US) call MOM_mesg(" DOME_initialization.F90, DOME_initialize_topography: setting topography", 5) - m_to_Z = 1.0 ; if (present(US)) m_to_Z = US%m_to_Z - call log_version(param_file, mdl, version, "") call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, & - "The minimum depth of the ocean.", units="m", default=0.0, scale=m_to_Z) + "The minimum depth of the ocean.", units="m", default=0.0, scale=US%m_to_Z) do j=js,je ; do i=is,ie if (G%geoLatT(i,j) < 600.0) then if (G%geoLatT(i,j) < 300.0) then D(i,j) = max_depth else - D(i,j) = max_depth - 10.0*m_to_Z * (G%geoLatT(i,j)-300.0) + D(i,j) = max_depth - 10.0*US%m_to_Z * (G%geoLatT(i,j)-300.0) endif else if ((G%geoLonT(i,j) > 1000.0) .AND. (G%geoLonT(i,j) < 1100.0)) then - D(i,j) = 600.0*m_to_Z + D(i,j) = 600.0*US%m_to_Z else D(i,j) = 0.5*min_depth endif diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index 617ac0da3d..2ebac05a68 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -42,14 +42,13 @@ module ISOMIP_initialization subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth, US) type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type real, dimension(G%isd:G%ied,G%jsd:G%jed), & - intent(out) :: D !< Ocean bottom depth in m or Z if US is present + intent(out) :: D !< Ocean bottom depth [m ~> Z] type(param_file_type), intent(in) :: param_file !< Parameter file structure - real, intent(in) :: max_depth !< Maximum model depth in the units of D - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + real, intent(in) :: max_depth !< Maximum model depth [m ~> Z] + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! Local variables real :: min_depth ! The minimum and maximum depths [Z ~> m]. - real :: m_to_Z ! A dimensional rescaling factor. ! The following variables are used to set up the bathymetry in the ISOMIP example. real :: bmax ! max depth of bedrock topography [Z ~> m] real :: b0,b2,b4,b6 ! first, second, third and fourth bedrock topography coeffs [Z ~> m] @@ -70,16 +69,14 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth, US) call MOM_mesg(" ISOMIP_initialization.F90, ISOMIP_initialize_topography: setting topography", 5) - m_to_Z = 1.0 ; if (present(US)) m_to_Z = US%m_to_Z - call log_version(param_file, mdl, version, "") call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, & - "The minimum depth of the ocean.", units="m", default=0.0, scale=m_to_Z) + "The minimum depth of the ocean.", units="m", default=0.0, scale=US%m_to_Z) call get_param(param_file, mdl, "ISOMIP_2D",is_2D,'If true, use a 2D setup.', default=.false.) ! The following variables should be transformed into runtime parameters? - bmax = 720.0*m_to_Z ; dc = 500.0*m_to_Z - b0 = -150.0*m_to_Z ; b2 = -728.8*m_to_Z ; b4 = 343.91*m_to_Z ; b6 = -50.57*m_to_Z + bmax = 720.0*US%m_to_Z ; dc = 500.0*US%m_to_Z + 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 xbar = 300.0e3 ; fc = 4.0e3 ; wc = 24.0e3 ; ly = 80.0e3 bx = 0.0 ; by = 0.0 ; xtil = 0.0 diff --git a/src/user/Kelvin_initialization.F90 b/src/user/Kelvin_initialization.F90 index 9bdf9b45c3..4c0c55f746 100644 --- a/src/user/Kelvin_initialization.F90 +++ b/src/user/Kelvin_initialization.F90 @@ -119,14 +119,13 @@ end subroutine Kelvin_OBC_end subroutine Kelvin_initialize_topography(D, G, param_file, max_depth, US) type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type real, dimension(G%isd:G%ied,G%jsd:G%jed), & - intent(out) :: D !< Ocean bottom depth in m or Z if US is present + intent(out) :: D !< Ocean bottom depth [m ~> Z] type(param_file_type), intent(in) :: param_file !< Parameter file structure - real, intent(in) :: max_depth !< Maximum model depth in the units of D [Z ~> m or m] - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + real, intent(in) :: max_depth !< Maximum model depth [Z ~> m] + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! Local variables character(len=40) :: mdl = "Kelvin_initialize_topography" ! This subroutine's name. - real :: m_to_Z ! A dimensional rescaling factor. real :: min_depth ! The minimum and maximum depths [Z ~> m]. real :: PI ! 3.1415... real :: coast_offset1, coast_offset2, coast_angle, right_angle @@ -134,10 +133,8 @@ subroutine Kelvin_initialize_topography(D, G, param_file, max_depth, US) call MOM_mesg(" Kelvin_initialization.F90, Kelvin_initialize_topography: setting topography", 5) - m_to_Z = 1.0 ; if (present(US)) m_to_Z = US%m_to_Z - call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, & - "The minimum depth of the ocean.", units="m", default=0.0, scale=m_to_Z) + "The minimum depth of the ocean.", units="m", default=0.0, scale=US%m_to_Z) call get_param(param_file, mdl, "ROTATED_COAST_OFFSET_1", coast_offset1, & default=100.0, do_not_log=.true.) call get_param(param_file, mdl, "ROTATED_COAST_OFFSET_2", coast_offset2, & diff --git a/src/user/Phillips_initialization.F90 b/src/user/Phillips_initialization.F90 index ed7bc07ba3..110a12c5f5 100644 --- a/src/user/Phillips_initialization.F90 +++ b/src/user/Phillips_initialization.F90 @@ -325,13 +325,12 @@ end function sech subroutine Phillips_initialize_topography(D, G, param_file, max_depth, US) type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type real, dimension(G%isd:G%ied,G%jsd:G%jed), & - intent(out) :: D !< Ocean bottom depth in m or Z if US is present + intent(out) :: D !< Ocean bottom depth [m ~> Z] type(param_file_type), intent(in) :: param_file !< Parameter file structure - real, intent(in) :: max_depth !< Maximum model depth in the units of D - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + real, intent(in) :: max_depth !< Maximum model depth [Z ~> m] + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! Local variables - real :: m_to_Z ! A dimensional rescaling factor. real :: PI, Htop, Wtop, Ltop, offset, dist real :: x1, x2, x3, x4, y1, y2 integer :: i,j,is,ie,js,je @@ -340,10 +339,9 @@ subroutine Phillips_initialize_topography(D, G, param_file, max_depth, US) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec PI = 4.0*atan(1.0) - m_to_Z = 1.0 ; if (present(US)) m_to_Z = US%m_to_Z call get_param(param_file, mdl, "PHILLIPS_HTOP", Htop, & - "The maximum height of the topography.", units="m", scale=m_to_Z, & + "The maximum height of the topography.", units="m", scale=US%m_to_Z, & fail_if_missing=.true.) ! Htop=0.375*max_depth ! max height of topog. above max_depth Wtop = 0.5*G%len_lat ! meridional width of drake and mount diff --git a/src/user/benchmark_initialization.F90 b/src/user/benchmark_initialization.F90 index b955f75a32..7b46295c20 100644 --- a/src/user/benchmark_initialization.F90 +++ b/src/user/benchmark_initialization.F90 @@ -34,17 +34,16 @@ module benchmark_initialization subroutine benchmark_initialize_topography(D, G, param_file, max_depth, US) type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type real, dimension(G%isd:G%ied,G%jsd:G%jed), & - intent(out) :: D !< Ocean bottom depth in m or [Z ~> m] if US is present + intent(out) :: D !< Ocean bottom depth [Z ~> m] type(param_file_type), intent(in) :: param_file !< Parameter file structure - real, intent(in) :: max_depth !< Maximum model depth in the units of D - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + real, intent(in) :: max_depth !< Maximum model depth [Z ~> m] + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! Local variables - real :: min_depth ! The minimum and maximum depths [Z ~> m]. + real :: min_depth ! The minimum depth [Z ~> m] real :: PI ! 3.1415926... calculated as 4*atan(1) real :: D0 ! A constant to make the maximum ! ! basin depth MAXIMUM_DEPTH. ! - real :: m_to_Z ! A dimensional rescaling factor. real :: x, y ! This include declares and sets the variable "version". # include "version_variable.h" @@ -55,11 +54,9 @@ subroutine benchmark_initialize_topography(D, G, param_file, max_depth, US) call MOM_mesg(" benchmark_initialization.F90, benchmark_initialize_topography: setting topography", 5) - m_to_Z = 1.0 ; if (present(US)) m_to_Z = US%m_to_Z - call log_version(param_file, mdl, version, "") call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, & - "The minimum depth of the ocean.", units="m", default=0.0, scale=m_to_Z) + "The minimum depth of the ocean.", units="m", default=0.0, scale=US%m_to_Z) PI = 4.0*atan(1.0) D0 = max_depth / 0.5 diff --git a/src/user/shelfwave_initialization.F90 b/src/user/shelfwave_initialization.F90 index 041d77d9f9..2c84a6040c 100644 --- a/src/user/shelfwave_initialization.F90 +++ b/src/user/shelfwave_initialization.F90 @@ -28,8 +28,8 @@ module shelfwave_initialization !> Control structure for shelfwave open boundaries. type, public :: shelfwave_OBC_CS ; private - real :: Lx = 100.0 !< Long-shore length scale of bathymetry. - real :: Ly = 50.0 !< Cross-shore length scale. + real :: Lx = 100.0 !< Long-shore length scale of bathymetry [km] + real :: Ly = 50.0 !< Cross-shore length scale [km] real :: f0 = 1.e-4 !< Coriolis parameter [T-1 ~> s-1] real :: jj = 1 !< Cross-shore wave mode. real :: kk !< Parameter. @@ -101,22 +101,19 @@ end subroutine shelfwave_OBC_end subroutine shelfwave_initialize_topography( D, G, param_file, max_depth, US ) type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type real, dimension(G%isd:G%ied,G%jsd:G%jed), & - intent(out) :: D !< Ocean bottom depth in m or Z if US is present + intent(out) :: D !< Ocean bottom depth [Z ~> m] type(param_file_type), intent(in) :: param_file !< Parameter file structure - real, intent(in) :: max_depth !< Maximum model depth in the units of D - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + real, intent(in) :: max_depth !< Maximum model depth [Z ~> m] + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! Local variables - real :: m_to_Z ! A dimensional rescaling factor. integer :: i, j real :: y, rLy, Ly, H0 - m_to_Z = 1.0 ; if (present(US)) m_to_Z = US%m_to_Z - call get_param(param_file, mdl,"SHELFWAVE_Y_LENGTH_SCALE",Ly, & default=50., do_not_log=.true.) call get_param(param_file, mdl,"MINIMUM_DEPTH", H0, & - default=10., units="m", scale=m_to_Z, do_not_log=.true.) + default=10., units="m", scale=US%m_to_Z, do_not_log=.true.) rLy = 0. ; if (Ly>0.) rLy = 1. / Ly diff --git a/src/user/user_initialization.F90 b/src/user/user_initialization.F90 index 18b1fa5225..d59d271471 100644 --- a/src/user/user_initialization.F90 +++ b/src/user/user_initialization.F90 @@ -61,10 +61,10 @@ end subroutine USER_set_coord subroutine USER_initialize_topography(D, G, param_file, max_depth, US) type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type real, dimension(G%isd:G%ied,G%jsd:G%jed), & - intent(out) :: D !< Ocean bottom depth in m or Z if US is present + intent(out) :: D !< Ocean bottom depth [Z ~> m] type(param_file_type), intent(in) :: param_file !< Parameter file structure - real, intent(in) :: max_depth !< Maximum model depth in the units of D - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + real, intent(in) :: max_depth !< Maximum model depth [Z ~> m] + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type call MOM_error(FATAL, & "USER_initialization.F90, USER_initialize_topography: " // &