diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 35dcdaa819..94d7852851 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -1275,8 +1275,8 @@ subroutine build_rho_grid( G, GV, US, h, tv, dzInterface, remapCS, CS, frac_shel !! [H ~> m or kg m-2] type(remapping_CS), intent(in) :: remapCS !< The remapping control structure type(regridding_CS), intent(in) :: CS !< Regridding control structure - real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: frac_shelf_h !< Fractional - !! ice shelf coverage [nodim] + real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: frac_shelf_h !< Fractional ice + !! shelf coverage [nondim] ! Local variables integer :: nz integer :: i, j, k @@ -1412,14 +1412,14 @@ subroutine build_grid_HyCOM1( G, GV, US, h, tv, h_new, dzInterface, CS, frac_she type(regridding_CS), intent(in) :: CS !< Regridding control structure real, dimension(SZI_(G),SZJ_(G),CS%nk), intent(inout) :: h_new !< New layer thicknesses [H ~> m or kg m-2] real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< Changes in interface position - real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: frac_shelf_h !< Fractional - !! ice shelf coverage [nodim] + real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: frac_shelf_h !< Fractional ice shelf + !! coverage [nondim] ! Local variables real, dimension(SZK_(GV)+1) :: z_col ! Source interface positions relative to the surface [H ~> m or kg m-2] real, dimension(CS%nk+1) :: z_col_new ! New interface positions relative to the surface [H ~> m or kg m-2] real, dimension(SZK_(GV)+1) :: dz_col ! The realized change in z_col [H ~> m or kg m-2] - real, dimension(SZK_(GV)) :: p_col ! Layer center pressure [Pa] + real, dimension(SZK_(GV)) :: p_col ! Layer center pressure [R L2 T-2 ~> Pa] integer :: i, j, k, nki real :: depth, nominalDepth real :: h_neglect, h_neglect_edge diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 3dc12c57e7..d3c3570ca2 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1711,7 +1711,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz integer :: IsdB, IedB, JsdB, JedB - real :: dtbt ! The barotropic timestep [s] + real :: dtbt ! If negative, this specifies the barotropic timestep as a fraction + ! of the maximum stable value [nondim]. real, allocatable, dimension(:,:) :: eta ! free surface height or column mass [H ~> m or kg m-2] real, allocatable, dimension(:,:) :: area_shelf_in ! area occupied by ice shelf [L2 ~> m2] diff --git a/src/core/MOM_CoriolisAdv.F90 b/src/core/MOM_CoriolisAdv.F90 index 08deeb7d0e..77b6305919 100644 --- a/src/core/MOM_CoriolisAdv.F90 +++ b/src/core/MOM_CoriolisAdv.F90 @@ -25,6 +25,7 @@ module MOM_CoriolisAdv !> Control structure for mom_coriolisadv type, public :: CoriolisAdv_CS ; private + logical :: initialized = .false. !< True if this control structure has been initialized. integer :: Coriolis_Scheme !< Selects the discretization for the Coriolis terms. !! Valid values are: !! - SADOURNY75_ENERGY - Sadourny, 1975 @@ -245,6 +246,9 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS) ! v(is-1:ie+2,js-1:je+1), u(is-1:ie+1,js-1:je+2), h(is-1:ie+2,js-1:je+2), ! uh(is-1,ie,js:je+1) and vh(is:ie+1,js-1:je). + if (.not.CS%initialized) call MOM_error(FATAL, & + "MOM_CoriolisAdv: Module must be initialized before it is used.") + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nz = GV%ke vol_neglect = GV%H_subroundoff * (1e-4 * US%m_to_L)**2 @@ -1123,6 +1127,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS) isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = GV%ke IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB + CS%initialized = .true. CS%diag => diag ; CS%Time => Time ! Read all relevant parameters and write them to the model log. diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90 index 112730fb59..8436b92f80 100644 --- a/src/core/MOM_PressureForce_FV.F90 +++ b/src/core/MOM_PressureForce_FV.F90 @@ -34,6 +34,7 @@ module MOM_PressureForce_FV !> Finite volume pressure gradient control structure type, public :: PressureForce_FV_CS ; private + logical :: initialized = .false. !< True if this control structure has been initialized. logical :: tides !< If true, apply tidal momentum forcing. real :: Rho0 !< The density used in the Boussinesq !! approximation [R ~> kg m-3]. @@ -163,6 +164,9 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB EOSdom(1) = Isq - (G%isd-1) ; EOSdom(2) = G%iec+1 - (G%isd-1) + if (.not.CS%initialized) call MOM_error(FATAL, & + "MOM_PressureForce_FV_nonBouss: Module must be initialized before it is used.") + if (CS%Stanley_T2_det_coeff>=0.) call MOM_error(FATAL, & "MOM_PressureForce_FV_nonBouss: The Stanley parameterization is not yet"//& "implemented in non-Boussinesq mode.") @@ -497,6 +501,9 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB EOSdom(1) = Isq - (G%isd-1) ; EOSdom(2) = G%iec+1 - (G%isd-1) + if (.not.CS%initialized) call MOM_error(FATAL, & + "MOM_PressureForce_FV_Bouss: Module must be initialized before it is used.") + use_p_atm = associated(p_atm) use_EOS = associated(tv%eqn_of_state) do i=Isq,Ieq+1 ; p0(i) = 0.0 ; enddo @@ -809,6 +816,7 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, tides_CS character(len=40) :: mdl ! This module's name. logical :: use_ALE + CS%initialized = .true. CS%diag => diag ; CS%Time => Time if (present(tides_CSp)) & CS%tides_CSp => tides_CSp diff --git a/src/core/MOM_PressureForce_Montgomery.F90 b/src/core/MOM_PressureForce_Montgomery.F90 index d77d31484a..6a0831eca9 100644 --- a/src/core/MOM_PressureForce_Montgomery.F90 +++ b/src/core/MOM_PressureForce_Montgomery.F90 @@ -30,6 +30,7 @@ module MOM_PressureForce_Mont !> Control structure for the Montgomery potential form of pressure gradient type, public :: PressureForce_Mont_CS ; private + logical :: initialized = .false. !< True if this control structure has been initialized. logical :: tides !< If true, apply tidal momentum forcing. real :: Rho0 !< The density used in the Boussinesq !! approximation [R ~> kg m-3]. @@ -137,6 +138,9 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb is_split = present(pbce) use_EOS = associated(tv%eqn_of_state) + if (.not.CS%initialized) call MOM_error(FATAL, & + "MOM_PressureForce_Mont: Module must be initialized before it is used.") + if (use_EOS) then if (query_compressible(tv%eqn_of_state)) call MOM_error(FATAL, & "PressureForce_Mont_nonBouss: The Montgomery form of the pressure force "//& @@ -422,6 +426,9 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce, is_split = present(pbce) use_EOS = associated(tv%eqn_of_state) + if (.not.CS%initialized) call MOM_error(FATAL, & + "MOM_PressureForce_Mont: Module must be initialized before it is used.") + if (use_EOS) then if (query_compressible(tv%eqn_of_state)) call MOM_error(FATAL, & "PressureForce_Mont_Bouss: The Montgomery form of the pressure force "//& @@ -829,6 +836,7 @@ subroutine PressureForce_Mont_init(Time, G, GV, US, param_file, diag, CS, tides_ # include "version_variable.h" character(len=40) :: mdl ! This module's name. + CS%initialized = .true. CS%diag => diag ; CS%Time => Time if (present(tides_CSp)) & CS%tides_CSp => tides_CSp diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 32eb036a94..ca1a7d20e5 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -190,7 +190,7 @@ module MOM_barotropic logical :: integral_bt_cont !< If true, use the time-integrated velocity over the barotropic steps !! to determine the integrated transports used to update the continuity !! equation. Otherwise the transports are the sum of the transports - !! based on ]a series of instantaneous velocities and the BT_CONT_TYPE + !! based on a series of instantaneous velocities and the BT_CONT_TYPE !! for transports. This is only valid if a BT_CONT_TYPE is used. logical :: Nonlinear_continuity !< If true, the barotropic continuity equation !! uses the full ocean thickness for transport. @@ -338,9 +338,9 @@ module MOM_barotropic !! the time-integrated barotropic velocity [L ~> m], beyond which the marginal !! open face area is FA_u_EE. uBT_EE must be non-positive. real :: uh_crvW !< The curvature of face area with velocity for flow from the west [H T2 L-1 ~> s2 or kg s2 m-3] - !! or [H L-1 ~> 1 or kg m-3] with INTEGRAL_BT_CONTINUITY. + !! or [H L-1 ~> nondim or kg m-3] with INTEGRAL_BT_CONTINUITY. real :: uh_crvE !< The curvature of face area with velocity for flow from the east [H T2 L-1 ~> s2 or kg s2 m-3] - !! or [H L-1 ~> 1 or kg m-3] with INTEGRAL_BT_CONTINUITY. + !! or [H L-1 ~> nondim or kg m-3] with INTEGRAL_BT_CONTINUITY. real :: uh_WW !< The zonal transport when ubt=ubt_WW [H L2 T-1 ~> m3 s-1 or kg s-1], or the equivalent !! time-integrated transport with INTEGRAL_BT_CONTINUITY [H L2 ~> m3 or kg]. real :: uh_EE !< The zonal transport when ubt=ubt_EE [H L2 T-1 ~> m3 s-1 or kg s-1], or the equivalent @@ -364,9 +364,9 @@ module MOM_barotropic !! the time-integrated barotropic velocity [L ~> m], beyond which the marginal !! open face area is FA_v_NN. vBT_NN must be non-positive. real :: vh_crvS !< The curvature of face area with velocity for flow from the south [H T2 L-1 ~> s2 or kg s2 m-3] - !! or [H L-1 ~> 1 or kg m-3] with INTEGRAL_BT_CONTINUITY. + !! or [H L-1 ~> nondim or kg m-3] with INTEGRAL_BT_CONTINUITY. real :: vh_crvN !< The curvature of face area with velocity for flow from the north [H T2 L-1 ~> s2 or kg s2 m-3] - !! or [H L-1 ~> 1 or kg m-3] with INTEGRAL_BT_CONTINUITY. + !! or [H L-1 ~> nondim or kg m-3] with INTEGRAL_BT_CONTINUITY. real :: vh_SS !< The meridional transport when vbt=vbt_SS [H L2 T-1 ~> m3 s-1 or kg s-1], or the equivalent !! time-integrated transport with INTEGRAL_BT_CONTINUITY [H L2 ~> m3 or kg]. real :: vh_NN !< The meridional transport when vbt=vbt_NN [H L2 T-1 ~> m3 s-1 or kg s-1], or the equivalent @@ -622,24 +622,23 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, vhbt_prev, vhbt_sum_prev, & ! Previous transports stored for OBCs [L2 H T-1 ~> m3 s-1] vbt_int_prev, & ! Previous value of time-integrated velocity stored for OBCs [L ~> m] vhbt_int_prev ! Previous value of time-integrated transport stored for OBCs [L2 H ~> m3] - real :: mass_to_Z ! The depth unit conversion divided by the mean density (Rho0) [Z m-1 R-1 ~> m3 kg-1] !### R-1 - real :: mass_accel_to_Z ! The inverse of the mean density (Rho0) [R-1 ~> m3 kg-1]. - real :: visc_rem ! A work variable that may equal visc_rem_[uv]. Nondim. + real :: mass_to_Z ! The inverse of the the mean density (Rho0) [R-1 ~> m3 kg-1] + real :: visc_rem ! A work variable that may equal visc_rem_[uv] [nondim] real :: vel_prev ! The previous velocity [L T-1 ~> m s-1]. real :: dtbt ! The barotropic time step [T ~> s]. real :: bebt ! A copy of CS%bebt [nondim]. real :: be_proj ! The fractional amount by which velocities are projected - ! when project_velocity is true. For now be_proj is set + ! when project_velocity is true [nondim]. For now be_proj is set ! to equal bebt, as they have similar roles and meanings. real :: Idt ! The inverse of dt [T-1 ~> s-1]. real :: det_de ! The partial derivative due to self-attraction and loading - ! of the reference geopotential with the sea surface height. + ! of the reference geopotential with the sea surface height [nondim]. ! This is typically ~0.09 or less. real :: dgeo_de ! The constant of proportionality between geopotential and - ! sea surface height. It is a nondimensional number of + ! sea surface height [nondim]. It is a nondimensional number of ! order 1. For stability, this may be made larger ! than the physical problem would suggest. - real :: Instep ! The inverse of the number of barotropic time steps to take. + real :: Instep ! The inverse of the number of barotropic time steps to take [nondim]. real :: wt_end ! The weighting of the final value of eta_PF [nondim] integer :: nstep ! The number of barotropic time steps to take. type(time_type) :: & @@ -694,6 +693,9 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, integer :: ioff, joff integer :: l_seg + if (.not.CS%module_is_initialized) call MOM_error(FATAL, & + "btstep: Module MOM_barotropic must be initialized before it is used.") + if (.not.CS%split) return is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB @@ -769,8 +771,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, Idtbt = 1.0 / dtbt bebt = CS%bebt be_proj = CS%bebt - mass_accel_to_Z = 1.0 / GV%Rho0 - mass_to_Z = US%m_to_Z / GV%Rho0 !### THis should be the same as mass_accel_to_Z. + mass_to_Z = 1.0 / GV%Rho0 !--- setup the weight when computing vbt_trans and ubt_trans if (project_velocity) then @@ -1240,7 +1241,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, endif endif - BT_force_u(I,j) = forces%taux(I,j) * mass_accel_to_Z * CS%IDatu(I,j)*visc_rem_u(I,j,1) + BT_force_u(I,j) = forces%taux(I,j) * mass_to_Z * CS%IDatu(I,j)*visc_rem_u(I,j,1) else BT_force_u(I,j) = 0.0 endif ; enddo ; enddo @@ -1266,11 +1267,11 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, endif endif - BT_force_v(i,J) = forces%tauy(i,J) * mass_accel_to_Z * CS%IDatv(i,J)*visc_rem_v(i,J,1) + BT_force_v(i,J) = forces%tauy(i,J) * mass_to_Z * CS%IDatv(i,J)*visc_rem_v(i,J,1) else BT_force_v(i,J) = 0.0 endif ; enddo ; enddo - if (associated(taux_bot) .and. associated(tauy_bot)) then + if (associated(taux_bot) .and. associated(tauy_bot)) then !$OMP parallel do default(shared) do j=js,je ; do I=is-1,ie ; if (G%mask2dCu(I,j) > 0.0) then BT_force_u(I,j) = BT_force_u(I,j) - taux_bot(I,j) * mass_to_Z * CS%IDatu(I,j) @@ -2764,6 +2765,9 @@ subroutine set_dtbt(G, GV, US, CS, eta, pbce, BT_cont, gtot_est, SSH_add) character(len=200) :: mesg integer :: i, j, k, is, ie, js, je, nz + if (.not.CS%module_is_initialized) call MOM_error(FATAL, & + "set_dtbt: Module MOM_barotropic must be initialized before it is used.") + if (.not.CS%split) return is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke MS%isdw = G%isd ; MS%iedw = G%ied ; MS%jsdw = G%jsd ; MS%jedw = G%jed @@ -3298,6 +3302,9 @@ subroutine btcalc(h, G, GV, CS, h_u, h_v, may_use_default, OBC) ! This section interpolates thicknesses onto u & v grid points with the ! second order accurate estimate h = 2*(h+ * h-)/(h+ + h-). + if (.not.CS%module_is_initialized) call MOM_error(FATAL, & + "btcalc: Module MOM_barotropic must be initialized before it is used.") + if (.not.CS%split) return use_default = .false. @@ -4186,6 +4193,9 @@ subroutine bt_mass_source(h, eta, set_cor, G, GV, CS) ! thicknesses [H ~> m or kg m-2]. integer :: is, ie, js, je, nz, i, j, k + if (.not.CS%module_is_initialized) call MOM_error(FATAL, "bt_mass_source: "// & + "Module MOM_barotropic must be initialized before it is used.") + if (.not.CS%split) return is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke diff --git a/src/core/MOM_checksum_packages.F90 b/src/core/MOM_checksum_packages.F90 index a1c300c94f..be00de8779 100644 --- a/src/core/MOM_checksum_packages.F90 +++ b/src/core/MOM_checksum_packages.F90 @@ -97,7 +97,7 @@ subroutine MOM_state_chksum_3arg(mesg, u, v, h, G, GV, US, haloshift, symmetric) 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 ~> nondim] or [nondim] + real :: L_T_to_m_s ! A rescaling factor for velocities [m T s-1 L-1 ~> 1] or [1] integer :: hs logical :: sym diff --git a/src/core/MOM_continuity_PPM.F90 b/src/core/MOM_continuity_PPM.F90 index cc56654d30..b7fb884667 100644 --- a/src/core/MOM_continuity_PPM.F90 +++ b/src/core/MOM_continuity_PPM.F90 @@ -26,6 +26,7 @@ module MOM_continuity_PPM !> Control structure for mom_continuity_ppm type, public :: continuity_PPM_CS ; private + logical :: initialized = .false. !< True if this control structure has been initialized. type(diag_ctrl), pointer :: diag !< Diagnostics control structure. logical :: upwind_1st !< If true, use a first-order upwind scheme. logical :: monotonic !< If true, use the Colella & Woodward monotonic @@ -134,6 +135,9 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, uhbt, vh h_min = GV%Angstrom_H + if (.not.CS%initialized) call MOM_error(FATAL, & + "MOM_continuity_PPM: Module must be initialized before it is used.") + x_first = (MOD(G%first_direction,2) == 0) if (present(visc_rem_u) .neqv. present(visc_rem_v)) call MOM_error(FATAL, & @@ -2202,6 +2206,8 @@ subroutine continuity_PPM_init(Time, G, GV, US, param_file, diag, CS) real :: tol_eta_m ! An unscaled version of tol_eta [m]. character(len=40) :: mdl = "MOM_continuity_PPM" ! This module's name. + CS%initialized = .true. + ! Read all relevant parameters and write them to the model log. call log_version(param_file, mdl, version, "") call get_param(param_file, mdl, "MONOTONIC_CONTINUITY", CS%monotonic, & diff --git a/src/core/MOM_density_integrals.F90 b/src/core/MOM_density_integrals.F90 index 9fb4fdabcc..3bdca94af3 100644 --- a/src/core/MOM_density_integrals.F90 +++ b/src/core/MOM_density_integrals.F90 @@ -802,7 +802,7 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & !! subtracted out to reduce the magnitude of each of the integrals. real, intent(in) :: rho_0 !< A density [R ~> kg m-3] or [kg m-3], that is used to calculate !! the pressure (as p~=-z*rho_0*G_e) used in the equation of state. - real, intent(in) :: G_e !< The Earth's gravitational acceleration [m s-2] + real, intent(in) :: G_e !< The Earth's gravitational acceleration [L2 Z-1 T-2 ~> m s-2] real, intent(in) :: dz_subroundoff !< A minuscule thickness change [Z ~> m] real, dimension(SZI_(HI),SZJ_(HI)), & intent(in) :: bathyT !< The depth of the bathymetry [Z ~> m] @@ -1457,7 +1457,7 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, real :: wt_t(5), wt_b(5) ! Weights of top and bottom values at quadrature points [nondim] real :: T_top, T_bot, S_top, S_bot, P_top, P_bot - real :: alpha_anom ! The depth averaged specific density anomaly [m3 kg-1] + real :: alpha_anom ! The depth averaged specific density anomaly [R-1 ~> m3 kg-1] or [m3 kg-1] real :: dp ! The pressure change through a layer [R L2 T-2 ~> Pa] real :: dp_90(2:4) ! The pressure change through a layer divided by 90 [R L2 T-2 ~> Pa] real :: hWght ! A pressure-thickness below topography [R L2 T-2 ~> Pa] diff --git a/src/core/MOM_grid.F90 b/src/core/MOM_grid.F90 index 7592dc8477..9c9ebf4960 100644 --- a/src/core/MOM_grid.F90 +++ b/src/core/MOM_grid.F90 @@ -171,10 +171,11 @@ module MOM_grid ! initialization routines (but not all) real :: south_lat !< The latitude (or y-coordinate) of the first v-line real :: west_lon !< The longitude (or x-coordinate) of the first u-line - real :: len_lat = 0. !< The latitudinal (or y-coord) extent of physical domain - real :: len_lon = 0. !< The longitudinal (or x-coord) extent of physical domain - real :: Rad_Earth = 6.378e6 !< The radius of the planet [m]. - real :: max_depth !< The maximum depth of the ocean in depth units [Z ~> m]. + real :: len_lat !< The latitudinal (or y-coord) extent of physical domain + real :: len_lon !< The longitudinal (or x-coord) extent of physical domain + real :: Rad_Earth !< The radius of the planet [m] + real :: Rad_Earth_L !< The radius of the planet in rescaled units [L ~> m] + real :: max_depth !< The maximum depth of the ocean in depth units [Z ~> m] end type ocean_grid_type contains diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 581cd5e68e..cb9422a412 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -3730,13 +3730,16 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) integer :: is_obc, ie_obc, js_obc, je_obc ! segment indices within local domain integer :: ishift, jshift ! offsets for staggered locations real, dimension(:,:,:), allocatable, target :: tmp_buffer - real, dimension(:), allocatable :: h_stack + real, dimension(:), allocatable :: h_stack ! Thicknesses at corner points [H ~> m or kg m-2] integer :: is_obc2, js_obc2 - real :: net_H_src, net_H_int, scl_fac - real :: tidal_vel, tidal_elev - real, allocatable :: normal_trans_bt(:,:) ! barotropic transport + real :: net_H_src ! Total thickness of the incoming flow in the source field [H ~> m or kg m-2] + real :: net_H_int ! Total thickness of the incoming flow in the model [H ~> m or kg m-2] + real :: scl_fac ! A nondimensional scaling factor [nondim] + real :: tidal_vel ! Tangential tidal velocity [m s-1] + real :: tidal_elev ! Tidal elevation at an OBC point [m] + real, allocatable :: normal_trans_bt(:,:) ! barotropic transport [H L2 T-1 ~> m3 s-1] integer :: turns ! Number of index quarter turns - real :: time_delta ! Time since tidal reference date + real :: time_delta ! Time since tidal reference date [s] is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed @@ -5110,7 +5113,7 @@ subroutine adjustSegmentEtaToFitBathymetry(G, GV, US, segment,fld) integer :: i, j, k, is, ie, js, je, nz, contractions, dilations integer :: n - real, allocatable, dimension(:,:,:) :: eta ! Segment source data interface heights, [Z -> m] + real, allocatable, dimension(:,:,:) :: eta ! Segment source data interface heights [Z ~> m] real :: hTolerance = 0.1 !< Tolerance to exceed adjustment criteria [Z ~> m] real :: hTmp, eTmp, dilate character(len=100) :: mesg diff --git a/src/core/MOM_transcribe_grid.F90 b/src/core/MOM_transcribe_grid.F90 index d3447f6590..58063f0669 100644 --- a/src/core/MOM_transcribe_grid.F90 +++ b/src/core/MOM_transcribe_grid.F90 @@ -126,7 +126,8 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG, US) 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 - oG%Rad_Earth = dG%Rad_Earth ; oG%max_depth = dG%max_depth + oG%Rad_Earth = dG%Rad_Earth ; oG%Rad_Earth_L = dG%Rad_Earth_L + oG%max_depth = dG%max_depth ! Update the halos in case the dynamic grid has smaller halos than the ocean grid. call pass_var(oG%areaT, oG%Domain) @@ -272,7 +273,8 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG, US) 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 - dG%Rad_Earth = oG%Rad_Earth ; dG%max_depth = oG%max_depth + dG%Rad_Earth = oG%Rad_Earth ; dG%Rad_Earth_L = oG%Rad_Earth_L + dG%max_depth = oG%max_depth ! Update the halos in case the dynamic grid has smaller halos than the ocean grid. call pass_var(dG%areaT, dG%Domain) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 363f3eebfb..c410ce094e 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -67,7 +67,7 @@ module MOM_variables logical :: T_is_conT = .false. !< If true, the temperature variable SST is actually the !! conservative temperature in [degC]. logical :: S_is_absS = .false. !< If true, the salinity variable SSS is actually the - !! absolute salinity in [g/kg]. + !! absolute salinity in [gSalt kg-1]. type(coupler_2d_bc_type) :: tr_fields !< A structure that may contain an !! array of named fields describing tracer-related quantities. !### NOTE: ALL OF THE ARRAYS IN TR_FIELDS USE THE COUPLER'S INDEXING CONVENTION AND HAVE NO @@ -95,7 +95,7 @@ module MOM_variables logical :: T_is_conT = .false. !< If true, the temperature variable tv%T is !! actually the conservative temperature [degC]. logical :: S_is_absS = .false. !< If true, the salinity variable tv%S is - !! actually the absolute salinity in units of [gSalt/kg]. + !! actually the absolute salinity in units of [gSalt kg-1]. real :: min_salinity = 0.01 !< The minimum value of salinity when BOUND_SALINITY=True [ppt]. !! The default is 0.01 for backward compatibility but should be 0. ! These arrays are accumulated fluxes for communication with other components. diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index bcee812c73..439ed242b8 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -49,6 +49,7 @@ module MOM_diagnostics !> The control structure for the MOM_diagnostics module type, public :: diagnostics_CS ; private + logical :: initialized = .false. !< True if this control structure has been initialized. real :: mono_N2_column_fraction = 0. !< The lower fraction of water column over which N2 is limited as !! monotonic for the purposes of calculating the equivalent !! barotropic wave speed. @@ -94,8 +95,8 @@ module MOM_diagnostics ! The following arrays hold diagnostics in the layer-integrated energy budget. real, allocatable :: KE(:,:,:) !< KE per unit mass [L2 T-2 ~> m2 s-2] real, allocatable :: dKE_dt(:,:,:) !< time derivative of the layer KE [H L2 T-3 ~> m3 s-3] - real, allocatable :: PE_to_KE(:,:,:) !< potential energy to KE term [m3 s-3] - real, allocatable :: KE_BT(:,:,:) !< barotropic contribution to KE term [m3 s-3] + real, allocatable :: PE_to_KE(:,:,:) !< potential energy to KE term [H L2 T-3 ~> m3 s-3] + real, allocatable :: KE_BT(:,:,:) !< barotropic contribution to KE term [H L2 T-3 ~> m3 s-3] real, allocatable :: KE_CorAdv(:,:,:) !< KE source from the combined Coriolis and !! advection terms [H L2 T-3 ~> m3 s-3]. !! The Coriolis source should be zero, but is not due to truncation @@ -267,6 +268,9 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & if (loc(CS)==0) call MOM_error(FATAL, & "calculate_diagnostic_fields: Module must be initialized before used.") + if (.not. CS%initialized) call MOM_error(FATAL, & + "calculate_diagnostic_fields: Module must be initialized before used.") + call calculate_derivs(dt, G, CS) if (dt > 0.0) then @@ -418,7 +422,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & call post_data(CS%id_masso, masso, CS%diag) endif - ! diagnose thickness/volumes of grid cells [m] + ! diagnose thickness/volumes of grid cells [Z ~> m] and [m3] if (CS%id_thkcello>0 .or. CS%id_volcello>0) then if (GV%Boussinesq) then ! thkcello = h for Boussinesq if (CS%id_thkcello > 0) then ; if (GV%H_to_Z == 1.0) then @@ -899,7 +903,7 @@ subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS) ! at the ocean surface [R L2 T-2 ~> Pa]. dpress, & ! Change in hydrostatic pressure across a layer [R L2 T-2 ~> Pa]. tr_int ! vertical integral of a tracer times density, - ! (Rho_0 in a Boussinesq model) [TR kg m-2]. + ! (Rho_0 in a Boussinesq model) [Conc R Z ~> Conc kg m-2]. real :: IG_Earth ! Inverse of gravitational acceleration [T2 Z L-2 ~> s2 m-1]. integer :: i, j, k, is, ie, js, je, nz @@ -1245,6 +1249,9 @@ subroutine register_time_deriv(lb, f_ptr, deriv_ptr, CS) integer :: m !< New index of deriv_ptr in CS%deriv integer :: ub(3) !< Upper index bound of f_ptr, based on shape. + if (.not.CS%initialized) call MOM_error(FATAL, & + "register_time_deriv: Module must be initialized before it is used.") + if (CS%num_time_deriv >= MAX_FIELDS_) then call MOM_error(WARNING,"MOM_diagnostics: Attempted to register more than " // & "MAX_FIELDS_ diagnostic time derivatives via register_time_deriv.") @@ -1604,6 +1611,8 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = GV%ke IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB + CS%initialized = .true. + CS%diag => diag use_temperature = associated(tv%T) call get_param(param_file, mdl, "ADIABATIC", adiabatic, default=.false., & diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90 index 602041372b..3b6fb0c510 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/src/diagnostics/MOM_sum_output.F90 @@ -59,6 +59,8 @@ module MOM_sum_output !> The control structure for the MOM_sum_output module type, public :: sum_output_CS ; private + logical :: initialized = .false. !< True if this control structure has been initialized. + type(Depth_List) :: DL !< The sorted depth list. integer, allocatable, dimension(:) :: lH @@ -145,11 +147,10 @@ subroutine MOM_sum_output_init(G, GV, US, param_file, directory, ntrnc, & type(Sum_output_CS), pointer :: CS !< A pointer that is set to point to the !! control structure for this module. ! Local variables - real :: Time_unit ! The time unit in seconds for ENERGYSAVEDAYS. - real :: Rho_0 ! A reference density [kg m-3] + real :: Time_unit ! The time unit in seconds for ENERGYSAVEDAYS [s] real :: maxvel ! The maximum permitted velocity [m s-1] -! This include declares and sets the variable "version". -#include "version_variable.h" + ! This include declares and sets the variable "version". +# include "version_variable.h" character(len=40) :: mdl = "MOM_sum_output" ! This module's name. character(len=200) :: energyfile ! The name of the energy file. character(len=32) :: filename_appendix = '' !fms appendix to filename for ensemble runs @@ -160,6 +161,8 @@ subroutine MOM_sum_output_init(G, GV, US, param_file, directory, ntrnc, & endif allocate(CS) + CS%initialized = .true. + ! Read all relevant parameters and write them to the model log. call log_version(param_file, mdl, version, "") call get_param(param_file, mdl, "CALCULATE_APE", CS%do_APE_calc, & @@ -377,11 +380,10 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci mass_anom_EFP ! The change in fresh water that cannot be accounted for by the surface ! fluxes [kg]. type(EFP_type), dimension(5) :: EFP_list ! An array of EFP types for joint global sums. - real :: CFL_Iarea ! Direction-based inverse area used in CFL test [L-2]. + real :: CFL_Iarea ! Direction-based inverse area used in CFL test [L-2 ~> m-2]. real :: CFL_trans ! A transport-based definition of the CFL number [nondim]. real :: CFL_lin ! A simpler definition of the CFL number [nondim]. real :: max_CFL(2) ! The maxima of the CFL numbers [nondim]. - real :: Irho0 ! The inverse of the reference density [m3 kg-1]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: & tmp1 ! A temporary array real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: & @@ -390,9 +392,9 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci Temp_int, Salt_int ! Layer and cell integrated heat and salt [J] and [g Salt]. real :: HL2_to_kg ! A conversion factor from a thickness-volume to mass [kg H-1 L-2 ~> kg m-3 or 1] real :: KE_scale_factor ! The combination of unit rescaling factors in the kinetic energy - ! calculation [kg T2 H-1 L-2 s-2 ~> kg m-3 or nondim] + ! calculation [kg T2 H-1 L-2 s-2 ~> kg m-3 or 1] real :: PE_scale_factor ! The combination of unit rescaling factors in the potential energy - ! calculation [kg T2 R-1 Z-1 L-2 s-2 ~> nondim] + ! calculation [kg T2 R-1 Z-1 L-2 s-2 ~> 1] integer :: num_nc_fields ! The number of fields that will actually go into ! the NetCDF file. integer :: i, j, k, is, ie, js, je, ns, nz, m, Isq, Ieq, Jsq, Jeq, isr, ier, jsr, jer @@ -490,6 +492,9 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci if (.not.associated(CS)) call MOM_error(FATAL, & "write_energy: Module must be initialized before it is used.") + if (.not.CS%initialized) call MOM_error(FATAL, & + "write_energy: Module must be initialized before it is used.") + do j=js,je ; do i=is,ie areaTm(i,j) = G%mask2dT(i,j)*G%areaT(i,j) enddo ; enddo @@ -734,7 +739,6 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci if (nTr_stocks > 0) call sum_across_PEs(Tr_stocks,nTr_stocks) call max_across_PEs(max_CFL, 2) - Irho0 = 1.0 / (US%R_to_kg_m3*GV%Rho0) if (CS%use_temperature) then if (CS%previous_calls == 0) then @@ -1031,7 +1035,7 @@ subroutine accumulate_net_input(fluxes, sfc_state, tv, dt, G, US, CS) ! enddo ; enddo ; endif if (associated(fluxes%salt_flux)) then ; do j=js,je ; do i=is,ie - ! convert salt_flux from kg (salt)/(m^2 s) to ppt * [m s-1]. + ! integrate salt_flux in [R Z T-1 ~> kgSalt m-2 s-1] to give [ppt kg] salt_in(i,j) = RZL2_to_kg * dt * & G%areaT(i,j)*(1000.0*fluxes%salt_flux(i,j)) enddo ; enddo ; endif diff --git a/src/diagnostics/MOM_wave_speed.F90 b/src/diagnostics/MOM_wave_speed.F90 index acec868561..6765f9aa12 100644 --- a/src/diagnostics/MOM_wave_speed.F90 +++ b/src/diagnostics/MOM_wave_speed.F90 @@ -26,6 +26,7 @@ module MOM_wave_speed !> Control structure for MOM_wave_speed type, public :: wave_speed_CS ; private + logical :: initialized = .false. !< True if this control structure has been initialized. logical :: use_ebt_mode = .false. !< If true, calculate the equivalent barotropic wave speed instead !! of the first baroclinic wave speed. !! This parameter controls the default behavior of wave_speed() which @@ -146,6 +147,9 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_ is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + if (.not. CS%initialized) call MOM_error(FATAL, "MOM_wave_speed: "// & + "Module must be initialized before it is used.") + if (present(full_halos)) then ; if (full_halos) then is = G%isd ; ie = G%ied ; js = G%jsd ; je = G%jed endif ; endif @@ -619,7 +623,7 @@ subroutine tdma6(n, a, c, lam, y) endif yy(k) = y(k) + a(k) * yy(k-1) * I_beta(k-1) enddo - ! The units of y change by a factor of [L2 T-2] in the following lines. + ! The units of y change by a factor of [L2 T-2 ~> m2 s-2] in the following lines. y(n) = yy(n) * I_beta(n) do k = n-1, 1, -1 y(k) = ( yy(k) + c(k) * y(k+1) ) * I_beta(k) @@ -722,6 +726,11 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, full_halos) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + if (present(CS)) then + if (.not. CS%initialized) call MOM_error(FATAL, "MOM_wave_speed: "// & + "Module must be initialized before it is used.") + endif + if (present(full_halos)) then ; if (full_halos) then is = G%isd ; ie = G%ied ; js = G%jsd ; je = G%jed endif ; endif @@ -1185,6 +1194,8 @@ subroutine wave_speed_init(CS, use_ebt_mode, mono_N2_column_fraction, mono_N2_de # include "version_variable.h" character(len=40) :: mdl = "MOM_wave_speed" ! This module's name. + CS%initialized = .true. + ! Write all relevant parameters to the model log. call log_version(mdl, version) diff --git a/src/diagnostics/MOM_wave_structure.F90 b/src/diagnostics/MOM_wave_structure.F90 index 3ae4f218d4..2dd272d409 100644 --- a/src/diagnostics/MOM_wave_structure.F90 +++ b/src/diagnostics/MOM_wave_structure.F90 @@ -36,6 +36,7 @@ module MOM_wave_structure !> The control structure for the MOM_wave_structure module type, public :: wave_structure_CS ; !private + logical :: initialized = .false. !< True if this control structure has been initialized. type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to !! regulate the timing of diagnostic output. real, allocatable, dimension(:,:,:) :: w_strct @@ -191,6 +192,9 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke I_a_int = 1/a_int + if (.not. CS%initialized) call MOM_error(FATAL, "MOM_wave_structure: "// & + "Module must be initialized before it is used.") + if (present(full_halos)) then ; if (full_halos) then is = G%isd ; ie = G%ied ; js = G%jsd ; je = G%jed endif ; endif @@ -728,6 +732,8 @@ subroutine wave_structure_init(Time, G, GV, param_file, diag, CS) isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = GV%ke + CS%initialized = .true. + call get_param(param_file, mdl, "INTERNAL_TIDE_SOURCE_X", CS%int_tide_source_x, & "X Location of generation site for internal tide", default=1.) call get_param(param_file, mdl, "INTERNAL_TIDE_SOURCE_Y", CS%int_tide_source_y, & diff --git a/src/framework/MOM_dyn_horgrid.F90 b/src/framework/MOM_dyn_horgrid.F90 index c7db67ee17..b67056cd80 100644 --- a/src/framework/MOM_dyn_horgrid.F90 +++ b/src/framework/MOM_dyn_horgrid.F90 @@ -163,10 +163,11 @@ module MOM_dyn_horgrid ! initialization routines (but not all) real :: south_lat !< The latitude (or y-coordinate) of the first v-line real :: west_lon !< The longitude (or x-coordinate) of the first u-line - real :: len_lat = 0. !< The latitudinal (or y-coord) extent of physical domain - real :: len_lon = 0. !< The longitudinal (or x-coord) extent of physical domain - real :: Rad_Earth = 6.378e6 !< The radius of the planet [m]. - real :: max_depth !< The maximum depth of the ocean [Z ~> m]. + real :: len_lat !< The latitudinal (or y-coord) extent of physical domain + real :: len_lon !< The longitudinal (or x-coord) extent of physical domain + real :: Rad_Earth !< The radius of the planet [m] + real :: Rad_Earth_L !< The radius of the planet in rescaled units [L ~> m] + real :: max_depth !< The maximum depth of the ocean [Z ~> m] end type dyn_horgrid_type contains @@ -361,6 +362,7 @@ subroutine rotate_dyn_horgrid(G_in, G, US, turns) G%areaT_global = G_in%areaT_global G%IareaT_global = G_in%IareaT_global G%Rad_Earth = G_in%Rad_Earth + G%Rad_Earth_L = G_in%Rad_Earth_L G%max_depth = G_in%max_depth call set_derived_dyn_horgrid(G, US) @@ -406,12 +408,8 @@ subroutine set_derived_dyn_horgrid(G, US) type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type ! Various inverse grid spacings and derived areas are calculated within this ! subroutine. - real :: m_to_L ! A unit conversion factor [L m-1 ~> nondim] - real :: L_to_m ! A unit conversion factor [L m-1 ~> nondim] integer :: i, j, isd, ied, jsd, jed integer :: IsdB, IedB, JsdB, JedB - 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 isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB diff --git a/src/framework/MOM_restart.F90 b/src/framework/MOM_restart.F90 index 019cfe135c..1328fd676c 100644 --- a/src/framework/MOM_restart.F90 +++ b/src/framework/MOM_restart.F90 @@ -69,6 +69,7 @@ module MOM_restart !> A restart registry and the control structure for restarts type, public :: MOM_restart_CS ; private + logical :: initialized = .false. !< True if this control structure has been initialized. logical :: restart !< restart is set to .true. if the run has been started from a full restart !! file. Otherwise some fields must be initialized approximately. integer :: novars = 0 !< The number of restart fields that have been registered. @@ -153,6 +154,9 @@ subroutine register_restart_field_ptr3d(f_ptr, var_desc, mandatory, CS) !! successfully read from the restart file. type(MOM_restart_CS), intent(inout) :: CS !< MOM restart control struct + if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // & + "register_restart_field: Module must be initialized before it is used.") + call lock_check(CS, var_desc) CS%novars = CS%novars+1 @@ -183,6 +187,9 @@ subroutine register_restart_field_ptr4d(f_ptr, var_desc, mandatory, CS) !! successfully read from the restart file. type(MOM_restart_CS), intent(inout) :: CS !< MOM restart control struct + if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // & + "register_restart_field: Module must be initialized before it is used.") + call lock_check(CS, var_desc) CS%novars = CS%novars+1 @@ -213,6 +220,9 @@ subroutine register_restart_field_ptr2d(f_ptr, var_desc, mandatory, CS) !! successfully read from the restart file. type(MOM_restart_CS), intent(inout) :: CS !< MOM restart control struct + if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // & + "register_restart_field: Module must be initialized before it is used.") + call lock_check(CS, var_desc) CS%novars = CS%novars+1 @@ -242,6 +252,9 @@ subroutine register_restart_field_ptr1d(f_ptr, var_desc, mandatory, CS) !! successfully read from the restart file. type(MOM_restart_CS), intent(inout) :: CS !< MOM restart control struct + if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // & + "register_restart_field: Module must be initialized before it is used.") + call lock_check(CS, var_desc) CS%novars = CS%novars+1 @@ -271,6 +284,9 @@ subroutine register_restart_field_ptr0d(f_ptr, var_desc, mandatory, CS) !! successfully read from the restart file. type(MOM_restart_CS), intent(inout) :: CS !< MOM restart control struct + if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // & + "register_restart_field: Module must be initialized before it is used.") + call lock_check(CS, var_desc) CS%novars = CS%novars+1 @@ -378,6 +394,10 @@ subroutine register_restart_field_4d(f_ptr, name, mandatory, CS, longname, units type(vardesc) :: vd + if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart: " // & + "register_restart_field_4d: Module must be initialized before "//& + "it is used to register "//trim(name)) + call lock_check(CS, name=name) vd = var_desc(name, units=units, longname=longname, hor_grid=hor_grid, & @@ -404,6 +424,10 @@ subroutine register_restart_field_3d(f_ptr, name, mandatory, CS, longname, units type(vardesc) :: vd + if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart: " // & + "register_restart_field_3d: Module must be initialized before "//& + "it is used to register "//trim(name)) + call lock_check(CS, name=name) vd = var_desc(name, units=units, longname=longname, hor_grid=hor_grid, & @@ -431,6 +455,10 @@ subroutine register_restart_field_2d(f_ptr, name, mandatory, CS, longname, units type(vardesc) :: vd character(len=8) :: Zgrid + if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart: " // & + "register_restart_field_2d: Module must be initialized before "//& + "it is used to register "//trim(name)) + zgrid = '1' ; if (present(z_grid)) zgrid = z_grid call lock_check(CS, name=name) @@ -459,6 +487,10 @@ subroutine register_restart_field_1d(f_ptr, name, mandatory, CS, longname, units type(vardesc) :: vd character(len=8) :: hgrid + if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart: " // & + "register_restart_field_3d: Module must be initialized before "//& + "it is used to register "//trim(name)) + hgrid = '1' ; if (present(hor_grid)) hgrid = hor_grid call lock_check(CS, name=name) @@ -484,6 +516,10 @@ subroutine register_restart_field_0d(f_ptr, name, mandatory, CS, longname, units type(vardesc) :: vd + if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart: " // & + "register_restart_field_0d: Module must be initialized before "//& + "it is used to register "//trim(name)) + call lock_check(CS, name=name) vd = var_desc(name, units=units, longname=longname, hor_grid='1', & @@ -503,6 +539,9 @@ function query_initialized_name(name, CS) result(query_initialized) integer :: m, n + if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // & + "query_initialized: Module must be initialized before it is used.") + if (CS%novars > CS%max_fields) call restart_error(CS) query_initialized = .false. @@ -533,6 +572,9 @@ function query_initialized_0d(f_ptr, CS) result(query_initialized) integer :: m, n + if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // & + "query_initialized: Module must be initialized before it is used.") + if (CS%novars > CS%max_fields) call restart_error(CS) query_initialized = .false. @@ -556,6 +598,9 @@ function query_initialized_1d(f_ptr, CS) result(query_initialized) integer :: m, n + if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // & + "query_initialized: Module must be initialized before it is used.") + if (CS%novars > CS%max_fields) call restart_error(CS) query_initialized = .false. @@ -580,6 +625,9 @@ function query_initialized_2d(f_ptr, CS) result(query_initialized) integer :: m, n + if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // & + "query_initialized: Module must be initialized before it is used.") + if (CS%novars > CS%max_fields) call restart_error(CS) query_initialized = .false. @@ -604,6 +652,9 @@ function query_initialized_3d(f_ptr, CS) result(query_initialized) integer :: m, n + if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // & + "query_initialized: Module must be initialized before it is used.") + if (CS%novars > CS%max_fields) call restart_error(CS) query_initialized = .false. @@ -628,6 +679,9 @@ function query_initialized_4d(f_ptr, CS) result(query_initialized) integer :: m, n + if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // & + "query_initialized: Module must be initialized before it is used.") + if (CS%novars > CS%max_fields) call restart_error(CS) query_initialized = .false. @@ -653,6 +707,9 @@ function query_initialized_0d_name(f_ptr, name, CS) result(query_initialized) integer :: m, n + if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // & + "query_initialized: Module must be initialized before it is used.") + if (CS%novars > CS%max_fields) call restart_error(CS) query_initialized = .false. @@ -685,6 +742,9 @@ function query_initialized_1d_name(f_ptr, name, CS) result(query_initialized) integer :: m, n + if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // & + "query_initialized: Module must be initialized before it is used.") + if (CS%novars > CS%max_fields) call restart_error(CS) query_initialized = .false. @@ -717,6 +777,9 @@ function query_initialized_2d_name(f_ptr, name, CS) result(query_initialized) integer :: m, n + if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // & + "query_initialized: Module must be initialized before it is used.") + if (CS%novars > CS%max_fields) call restart_error(CS) query_initialized = .false. @@ -749,6 +812,9 @@ function query_initialized_3d_name(f_ptr, name, CS) result(query_initialized) integer :: m, n + if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // & + "query_initialized: Module must be initialized before it is used.") + if (CS%novars > CS%max_fields) call restart_error(CS) query_initialized = .false. @@ -781,6 +847,9 @@ function query_initialized_4d_name(f_ptr, name, CS) result(query_initialized) integer :: m, n + if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // & + "query_initialized: Module must be initialized before it is used.") + if (CS%novars > CS%max_fields) call restart_error(CS) query_initialized = .false. @@ -850,6 +919,9 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV, num_ turns = CS%turns + if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // & + "save_restart: Module must be initialized before it is used.") + if (CS%novars > CS%max_fields) call restart_error(CS) ! With parallel read & write, it is possible to disable the following... @@ -1038,6 +1110,9 @@ subroutine restore_state(filename, directory, day, G, CS) integer(kind=8) :: checksum_file ! The checksum value recorded in the input file. integer(kind=8) :: checksum_data ! The checksum value for the data that was read in. + if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // & + "restore_state: Module must be initialized before it is used.") + if (CS%novars > CS%max_fields) call restart_error(CS) ! Get NetCDF ids for all of the restart files. @@ -1230,6 +1305,9 @@ function restart_files_exist(filename, directory, G, CS) !! restart files exist in directory integer :: num_files + if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // & + "restart_files_exist: Module must be initialized before it is used.") + if ((LEN_TRIM(filename) == 1) .and. (filename(1:1) == 'F')) then num_files = get_num_restart_files('r', directory, G, CS) else @@ -1251,6 +1329,9 @@ function determine_is_new_run(filename, directory, G, CS) result(is_new_run) !! this is a new run, based on the value of !! filename and whether restart files exist + if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // & + "determine_is_new_run: Module must be initialized before it is used.") + if (LEN_TRIM(filename) > 1) then CS%new_run = .false. elseif (LEN_TRIM(filename) == 0) then @@ -1275,6 +1356,9 @@ function is_new_run(CS) logical :: is_new_run !< The function result, which had been stored in CS during !! a previous call to determine_is_new_run + if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // & + "is_new_run: Module must be initialized before it is used.") + if (.not.CS%new_run_set) call MOM_error(FATAL, "MOM_restart " // & "determine_is_new_run must be called for a restart file before is_new_run.") @@ -1317,6 +1401,9 @@ function open_restart_units(filename, directory, G, CS, IO_handles, file_paths, character(len=32) :: filename_appendix = '' ! Filename appendix for ensemble runs character(len=80) :: restartname + if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // & + "open_restart_units: Module must be initialized before it is used.") + ! Get NetCDF ids for all of the restart files. num_restart = 0 ; nf = 0 ; start_char = 1 do while (start_char <= len_trim(filename) ) @@ -1429,6 +1516,9 @@ function get_num_restart_files(filenames, directory, G, CS, file_paths) result(n integer :: num_files !< The function result, the number of files (both automatically named !! restart files and others explicitly in filename) that have been opened + if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // & + "get_num_restart_files: Module must be initialized before it is used.") + ! This call uses open_restart_units without the optional arguments needed to actually ! open the files to determine the number of restart files. num_files = open_restart_units(filenames, directory, G, CS, file_paths=file_paths) @@ -1458,6 +1548,8 @@ subroutine restart_init(param_file, CS, restart_root) endif allocate(CS) + CS%initialized = .true. + ! Determine whether all paramters are set to their default values. call get_param(param_file, mdl, "PARALLEL_RESTARTFILES", CS%parallel_restartfiles, & default=.false., do_not_log=.true.) diff --git a/src/framework/MOM_write_cputime.F90 b/src/framework/MOM_write_cputime.F90 index c9200cf41c..9df994448b 100644 --- a/src/framework/MOM_write_cputime.F90 +++ b/src/framework/MOM_write_cputime.F90 @@ -20,6 +20,7 @@ module MOM_write_cputime !> A control structure that regulates the writing of CPU time type, public :: write_cputime_CS ; private + logical :: initialized = .false. !< True if this control structure has been initialized. real :: maxcpu !< The maximum amount of cpu time per processor !! for which MOM should run before saving a restart !! file and quiting with a return value that @@ -71,6 +72,8 @@ subroutine MOM_write_cputime_init(param_file, directory, Input_start_time, CS) CS%prev_cputime = new_cputime endif + CS%initialized = .true. + ! Read all relevant parameters and write them to the model log. ! Determine whether all paramters are set to their default values. @@ -141,6 +144,9 @@ subroutine write_cputime(day, n, CS, nmax, call_end) if (.not.associated(CS)) call MOM_error(FATAL, & "write_energy: Module must be initialized before it is used.") + if (.not.CS%initialized) call MOM_error(FATAL, & + "write_cputime: Module must be initialized before it is used.") + call SYSTEM_CLOCK(new_cputime, CLOCKS_PER_SEC, MAX_TICKS) ! The following lines extract useful information even if the clock has rolled ! over, assuming a 32-bit SYSTEM_CLOCK. With more bits, rollover is essentially diff --git a/src/initialization/MOM_coord_initialization.F90 b/src/initialization/MOM_coord_initialization.F90 index eb0db68726..286dfa7d95 100644 --- a/src/initialization/MOM_coord_initialization.F90 +++ b/src/initialization/MOM_coord_initialization.F90 @@ -309,13 +309,13 @@ subroutine set_coord_from_TS_range(Rlay, g_prime, GV, US, param_file, eqn_of_sta ! Local variables real, dimension(GV%ke) :: T0, S0, Pref real :: S_Ref, S_Light, S_Dense ! Salinity range parameters [ppt]. - real :: T_Ref, T_Light, T_Dense ! Temperature range parameters [decC]. + real :: T_Ref, T_Light, T_Dense ! Temperature range parameters [degC]. real :: res_rat ! The ratio of density space resolution in the denser part ! of the range to that in the lighter part of the range. ! Setting this greater than 1 increases the resolution for - ! the denser water. + ! the denser water [nondim]. real :: g_fs ! Reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2]. - real :: a1, frac_dense, k_frac + real :: a1, frac_dense, k_frac ! Nondimensional temporary variables [nondim] integer :: k, nz, k_light character(len=40) :: mdl = "set_coord_from_TS_range" ! This subroutine's name. character(len=200) :: filename, coord_file, inputdir ! Strings for file/path diff --git a/src/initialization/MOM_fixed_initialization.F90 b/src/initialization/MOM_fixed_initialization.F90 index 069d576b2c..b67d21ebcb 100644 --- a/src/initialization/MOM_fixed_initialization.F90 +++ b/src/initialization/MOM_fixed_initialization.F90 @@ -144,7 +144,7 @@ subroutine MOM_initialize_fixed(G, US, OBC, PF, write_geom, output_dir) endif ! Calculate the value of the Coriolis parameter at the latitude ! -! of the q grid points [s-1]. +! of the q grid points [T-1 ~> s-1]. call MOM_initialize_rotation(G%CoriolisBu, G, PF, US=US) ! Calculate the components of grad f (beta) call MOM_calculate_grad_Coriolis(G%dF_dx, G%dF_dy, G, US=US) @@ -167,13 +167,13 @@ subroutine MOM_initialize_fixed(G, US, OBC, PF, write_geom, output_dir) end subroutine MOM_initialize_fixed !> MOM_initialize_topography makes the appropriate call to set up the bathymetry. At this -!! point the topography is in units of [m], but this can be changed later. +!! point the topography is in units of [Z ~> m] or [m], depending on the presence of US. subroutine MOM_initialize_topography(D, max_depth, G, PF, 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 [m] + 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 [m] + 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 ! This subroutine makes the appropriate call to set up the bottom depth. diff --git a/src/initialization/MOM_grid_initialize.F90 b/src/initialization/MOM_grid_initialize.F90 index 81e7b66d7a..f67a977d27 100644 --- a/src/initialization/MOM_grid_initialize.F90 +++ b/src/initialization/MOM_grid_initialize.F90 @@ -26,14 +26,14 @@ module MOM_grid_initialize ! vary with the Boussinesq approximation, the Boussinesq variant is given first. !> Global positioning system (aka container for information to describe the grid) -type, public :: GPS ; private +type, private :: GPS ; private real :: len_lon !< The longitudinal or x-direction length of the domain. real :: len_lat !< The latitudinal or y-direction length of the domain. real :: west_lon !< The western longitude of the domain or the equivalent !! starting value for the x-axis. real :: south_lat !< The southern latitude of the domain or the equivalent !! starting value for the y-axis. - real :: Rad_Earth !< The radius of the Earth [m]. + real :: Rad_Earth_L !< The radius of the Earth in rescaled units [L ~> m] real :: Lat_enhance_factor !< The amount by which the meridional resolution !! is enhanced within LAT_EQ_ENHANCE of the equator. real :: Lat_eq_enhance !< The latitude range to the north and south of the equator @@ -59,12 +59,18 @@ 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 + ! Local variables -! This include declares and sets the variable "version". -#include "version_variable.h" + 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, & @@ -82,6 +88,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 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) @@ -93,6 +100,15 @@ subroutine set_grid_metrics(G, param_file, US) case default ; call MOM_error(FATAL, "MOM_grid_init: set_grid_metrics "//& "Unrecognized grid configuration "//trim(config)) end select + if (G%Rad_Earth_L <= 0.0) then + ! 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) + ! 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 + endif + G%Rad_Earth = 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") @@ -113,8 +129,8 @@ subroutine grid_metrics_chksum(parent, G, US) 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 - real :: m_to_L ! A unit conversion factor [L m-1 ~> nondim] - real :: L_to_m ! A unit conversion factor [m L-1 ~> nondim] + 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 @@ -181,7 +197,7 @@ 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 ~> nondim] + 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 @@ -386,11 +402,10 @@ subroutine set_grid_metrics_cartesian(G, param_file, US) integer :: niglobal, njglobal 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 :: dx_everywhere, dy_everywhere ! Grid spacings [m]. - real :: I_dx, I_dy ! Inverse grid spacings [m-1]. + 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 ~> nondim] - real :: L_to_m ! A unit conversion factor [m L-1 ~> nondim] + 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" @@ -402,7 +417,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 - L_to_m = 1.0 ; if (present(US)) L_to_m = US%L_to_m PI = 4.0*atan(1.0) call get_param(param_file, mdl, "AXIS_UNITS", units_temp, & @@ -423,8 +437,8 @@ subroutine set_grid_metrics_cartesian(G, param_file, US) call get_param(param_file, mdl, "LENLON", G%len_lon, & "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, & - "The radius of the Earth.", units="m", default=6.378e6) + 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) if (units_temp(1:1) == 'k') then G%x_axis_units = "kilometers" ; G%y_axis_units = "kilometers" @@ -462,14 +476,14 @@ 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 * G%len_lon / (REAL(niglobal)) - dy_everywhere = 1000.0 * G%len_lat / (REAL(njglobal)) + dx_everywhere = 1000.0*m_to_L * G%len_lon / (REAL(niglobal)) + dy_everywhere = 1000.0*m_to_L * G%len_lat / (REAL(njglobal)) elseif (units_temp(1:1) == 'm') then ! Axes are measured in m. - dx_everywhere = G%len_lon / (REAL(niglobal)) - dy_everywhere = G%len_lat / (REAL(njglobal)) + dx_everywhere = m_to_L*G%len_lon / (REAL(niglobal)) + dy_everywhere = m_to_L*G%len_lat / (REAL(njglobal)) else ! Axes are measured in degrees of latitude and longitude. - dx_everywhere = G%Rad_Earth * G%len_lon * PI / (180.0 * niglobal) - dy_everywhere = G%Rad_Earth * G%len_lat * PI / (180.0 * njglobal) + 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) endif I_dx = 1.0 / dx_everywhere ; I_dy = 1.0 / dy_everywhere @@ -477,30 +491,30 @@ subroutine set_grid_metrics_cartesian(G, param_file, US) do J=JsdB,JedB ; do I=IsdB,IedB G%geoLonBu(I,J) = grid_lonB(I) ; G%geoLatBu(I,J) = grid_latB(J) - G%dxBu(I,J) = m_to_L*dx_everywhere ; G%IdxBu(I,J) = L_to_m*I_dx - G%dyBu(I,J) = m_to_L*dy_everywhere ; G%IdyBu(I,J) = L_to_m*I_dy - G%areaBu(I,J) = m_to_L**2*dx_everywhere * dy_everywhere ; G%IareaBu(I,J) = L_to_m**2*I_dx * I_dy + G%dxBu(I,J) = dx_everywhere ; G%IdxBu(I,J) = I_dx + G%dyBu(I,J) = dy_everywhere ; G%IdyBu(I,J) = I_dy + G%areaBu(I,J) = dx_everywhere * dy_everywhere ; G%IareaBu(I,J) = I_dx * I_dy enddo ; enddo do j=jsd,jed ; do i=isd,ied G%geoLonT(i,j) = grid_lonT(i) ; G%geoLatT(i,j) = grid_LatT(j) - G%dxT(i,j) = m_to_L*dx_everywhere ; G%IdxT(i,j) = L_to_m*I_dx - G%dyT(i,j) = m_to_L*dy_everywhere ; G%IdyT(i,j) = L_to_m*I_dy - G%areaT(i,j) = m_to_L**2*dx_everywhere * dy_everywhere ; G%IareaT(i,j) = L_to_m**2*I_dx * I_dy + G%dxT(i,j) = dx_everywhere ; G%IdxT(i,j) = I_dx + G%dyT(i,j) = dy_everywhere ; G%IdyT(i,j) = I_dy + G%areaT(i,j) = dx_everywhere * dy_everywhere ; G%IareaT(i,j) = I_dx * I_dy enddo ; enddo do j=jsd,jed ; do I=IsdB,IedB G%geoLonCu(I,j) = grid_lonB(I) ; G%geoLatCu(I,j) = grid_LatT(j) - G%dxCu(I,j) = m_to_L*dx_everywhere ; G%IdxCu(I,j) = L_to_m*I_dx - G%dyCu(I,j) = m_to_L*dy_everywhere ; G%IdyCu(I,j) = L_to_m*I_dy + G%dxCu(I,j) = dx_everywhere ; G%IdxCu(I,j) = I_dx + G%dyCu(I,j) = dy_everywhere ; G%IdyCu(I,j) = I_dy enddo ; enddo do J=JsdB,JedB ; do i=isd,ied G%geoLonCv(i,J) = grid_lonT(i) ; G%geoLatCv(i,J) = grid_latB(J) - G%dxCv(i,J) = m_to_L*dx_everywhere ; G%IdxCv(i,J) = L_to_m*I_dx - G%dyCv(i,J) = m_to_L*dy_everywhere ; G%IdyCv(i,J) = L_to_m*I_dy + G%dxCv(i,J) = dx_everywhere ; G%IdxCv(i,J) = I_dx + G%dyCv(i,J) = dy_everywhere ; G%IdyCv(i,J) = I_dy enddo ; enddo call callTree_leave("set_grid_metrics_cartesian()") @@ -527,7 +541,7 @@ 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 ~> nondim] + 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 @@ -555,8 +569,8 @@ subroutine set_grid_metrics_spherical(G, param_file, US) call get_param(param_file, mdl, "LENLON", G%len_lon, & "The longitudinal length of the domain.", units="degrees", & fail_if_missing=.true.) - call get_param(param_file, mdl, "RAD_EARTH", G%Rad_Earth, & - "The radius of the Earth.", units="m", default=6.378e6) + 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) dLon = G%len_lon/G%Domain%niglobal dLat = G%len_lat/G%Domain%njglobal @@ -600,9 +614,9 @@ subroutine set_grid_metrics_spherical(G, param_file, US) ! The following line is needed to reproduce the solution from ! set_grid_metrics_mercator when used to generate a simple spherical grid. - G%dxBu(I,J) = m_to_L*G%Rad_Earth * COS( G%geoLatBu(I,J)*PI_180 ) * dL_di -! G%dxBu(I,J) = m_to_L*G%Rad_Earth * dLon*PI_180 * COS( G%geoLatBu(I,J)*PI_180 ) - G%dyBu(I,J) = m_to_L*G%Rad_Earth * dLat*PI_180 + G%dxBu(I,J) = G%Rad_Earth_L * COS( G%geoLatBu(I,J)*PI_180 ) * dL_di +! G%dxBu(I,J) = G%Rad_Earth_L * dLon*PI_180 * COS( G%geoLatBu(I,J)*PI_180 ) + G%dyBu(I,J) = G%Rad_Earth_L * dLat*PI_180 G%areaBu(I,J) = G%dxBu(I,J) * G%dyBu(I,J) enddo ; enddo @@ -612,9 +626,9 @@ subroutine set_grid_metrics_spherical(G, param_file, US) ! The following line is needed to reproduce the solution from ! set_grid_metrics_mercator when used to generate a simple spherical grid. - G%dxCv(i,J) = m_to_L*G%Rad_Earth * COS( G%geoLatCv(i,J)*PI_180 ) * dL_di -! G%dxCv(i,J) = m_to_L*G%Rad_Earth * (dLon*PI_180) * COS( G%geoLatCv(i,J)*PI_180 ) - G%dyCv(i,J) = m_to_L*G%Rad_Earth * dLat*PI_180 + G%dxCv(i,J) = G%Rad_Earth_L * COS( G%geoLatCv(i,J)*PI_180 ) * dL_di +! G%dxCv(i,J) = G%Rad_Earth_L * (dLon*PI_180) * COS( G%geoLatCv(i,J)*PI_180 ) + G%dyCv(i,J) = G%Rad_Earth_L * dLat*PI_180 enddo ; enddo do j=jsd,jed ; do I=IsdB,IedB @@ -623,9 +637,9 @@ subroutine set_grid_metrics_spherical(G, param_file, US) ! The following line is needed to reproduce the solution from ! set_grid_metrics_mercator when used to generate a simple spherical grid. - G%dxCu(I,j) = m_to_L*G%Rad_Earth * COS( G%geoLatCu(I,j)*PI_180 ) * dL_di -! G%dxCu(I,j) = m_to_L*G%Rad_Earth * dLon*PI_180 * COS( latitude ) - G%dyCu(I,j) = m_to_L*G%Rad_Earth * dLat*PI_180 + G%dxCu(I,j) = G%Rad_Earth_L * COS( G%geoLatCu(I,j)*PI_180 ) * dL_di +! G%dxCu(I,j) = G%Rad_Earth_L * dLon*PI_180 * COS( latitude ) + G%dyCu(I,j) = G%Rad_Earth_L * dLat*PI_180 enddo ; enddo do j=jsd,jed ; do i=isd,ied @@ -634,13 +648,13 @@ subroutine set_grid_metrics_spherical(G, param_file, US) ! The following line is needed to reproduce the solution from ! set_grid_metrics_mercator when used to generate a simple spherical grid. - G%dxT(i,j) = m_to_L*G%Rad_Earth * COS( G%geoLatT(i,j)*PI_180 ) * dL_di -! G%dxT(i,j) = G%Rad_Earth * dLon*PI_180 * COS( latitude ) - G%dyT(i,j) = m_to_L*G%Rad_Earth * dLat*PI_180 + G%dxT(i,j) = G%Rad_Earth_L * COS( G%geoLatT(i,j)*PI_180 ) * dL_di +! G%dxT(i,j) = G%Rad_Earth_L * dLon*PI_180 * COS( latitude ) + G%dyT(i,j) = G%Rad_Earth_L * dLat*PI_180 ! latitude = G%geoLatCv(i,J)*PI_180 ! In radians ! dL_di = G%geoLatCv(i,max(jsd,J-1))*PI_180 ! In radians -! G%areaT(i,j) = m_to_L**2 * Rad_Earth**2*dLon*dLat*ABS(SIN(latitude)-SIN(dL_di)) +! G%areaT(i,j) = Rad_Earth_L**2*dLon*dLat*ABS(SIN(latitude)-SIN(dL_di)) G%areaT(i,j) = G%dxT(i,j) * G%dyT(i,j) enddo ; enddo @@ -677,7 +691,7 @@ 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 ~> nondim] + 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 @@ -713,11 +727,12 @@ subroutine set_grid_metrics_mercator(G, param_file, US) call get_param(param_file, mdl, "LENLON", GP%len_lon, & "The longitudinal length of the domain.", units="degrees", & fail_if_missing=.true.) - call get_param(param_file, mdl, "RAD_EARTH", GP%Rad_Earth, & - "The radius of the Earth.", units="m", default=6.378e6) + 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) 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 = GP%Rad_Earth + G%Rad_Earth_L = GP%Rad_Earth_L + call get_param(param_file, mdl, "ISOTROPIC", GP%isotropic, & "If true, an isotropic grid on a sphere (also known as "//& "a Mercator grid) is used. With an isotropic grid, the "//& @@ -826,8 +841,8 @@ subroutine set_grid_metrics_mercator(G, param_file, US) do J=JsdB,JedB ; do I=IsdB,IedB G%geoLonBu(I,J) = xq(I,J)*180.0/PI G%geoLatBu(I,J) = yq(I,J)*180.0/PI - G%dxBu(I,J) = m_to_L*ds_di(xq(I,J), yq(I,J), GP) - G%dyBu(I,J) = m_to_L*ds_dj(xq(I,J), yq(I,J), GP) + G%dxBu(I,J) = ds_di(xq(I,J), yq(I,J), GP) + G%dyBu(I,J) = ds_dj(xq(I,J), yq(I,J), GP) G%areaBu(I,J) = G%dxBu(I,J) * G%dyBu(I,J) G%IareaBu(I,J) = 1.0 / (G%areaBu(I,J)) @@ -836,8 +851,8 @@ subroutine set_grid_metrics_mercator(G, param_file, US) do j=jsd,jed ; do i=isd,ied G%geoLonT(i,j) = xh(i,j)*180.0/PI G%geoLatT(i,j) = yh(i,j)*180.0/PI - G%dxT(i,j) = m_to_L*ds_di(xh(i,j), yh(i,j), GP) - G%dyT(i,j) = m_to_L*ds_dj(xh(i,j), yh(i,j), GP) + G%dxT(i,j) = ds_di(xh(i,j), yh(i,j), GP) + G%dyT(i,j) = ds_dj(xh(i,j), yh(i,j), GP) G%areaT(i,j) = G%dxT(i,j)*G%dyT(i,j) G%IareaT(i,j) = 1.0 / (G%areaT(i,j)) @@ -846,20 +861,20 @@ subroutine set_grid_metrics_mercator(G, param_file, US) do j=jsd,jed ; do I=IsdB,IedB G%geoLonCu(I,j) = xu(I,j)*180.0/PI G%geoLatCu(I,j) = yu(I,j)*180.0/PI - G%dxCu(I,j) = m_to_L*ds_di(xu(I,j), yu(I,j), GP) - G%dyCu(I,j) = m_to_L*ds_dj(xu(I,j), yu(I,j), GP) + G%dxCu(I,j) = ds_di(xu(I,j), yu(I,j), GP) + G%dyCu(I,j) = ds_dj(xu(I,j), yu(I,j), GP) enddo ; enddo do J=JsdB,JedB ; do i=isd,ied G%geoLonCv(i,J) = xv(i,J)*180.0/PI G%geoLatCv(i,J) = yv(i,J)*180.0/PI - G%dxCv(i,J) = m_to_L*ds_di(xv(i,J), yv(i,J), GP) - G%dyCv(i,J) = m_to_L*ds_dj(xv(i,J), yv(i,J), GP) + G%dxCv(i,J) = ds_di(xv(i,J), yv(i,J), GP) + G%dyCv(i,J) = ds_dj(xv(i,J), yv(i,J), GP) enddo ; enddo if (.not.simple_area) then do j=JsdB+1,jed ; do i=IsdB+1,ied - G%areaT(I,J) = m_to_L**2*GP%Rad_Earth**2 * & + G%areaT(I,J) = GP%Rad_Earth_L**2 * & (dL(xq(I-1,J-1),xq(I-1,J),yq(I-1,J-1),yq(I-1,J)) + & (dL(xq(I-1,J),xq(I,J),yq(I-1,J),yq(I,J)) + & (dL(xq(I,J),xq(I,J-1),yq(I,J),yq(I,J-1)) + & @@ -884,31 +899,31 @@ subroutine set_grid_metrics_mercator(G, param_file, US) end subroutine set_grid_metrics_mercator -!> This function returns the grid spacing in the logical x direction. +!> This function returns the grid spacing in the logical x direction in [L ~> m]. function ds_di(x, y, GP) real, intent(in) :: x !< The longitude in question real, intent(in) :: y !< The latitude in question type(GPS), intent(in) :: GP !< A structure of grid parameters - real :: ds_di - ! Local variables - ds_di = GP%Rad_Earth * cos(y) * dx_di(x,GP) + real :: ds_di ! The returned grid spacing [L ~> m] + + ds_di = GP%Rad_Earth_L * cos(y) * dx_di(x,GP) ! In general, this might be... - ! ds_di = GP%Rad_Earth * sqrt( cos(y)*cos(y) * dx_di(x,y,GP)*dx_di(x,y,GP) + & + ! ds_di = GP%Rad_Earth_L * sqrt( cos(y)*cos(y) * dx_di(x,y,GP)*dx_di(x,y,GP) + & ! dy_di(x,y,GP)*dy_di(x,y,GP)) end function ds_di -!> This function returns the grid spacing in the logical y direction. +!> This function returns the grid spacing in the logical y direction in [L ~> m]. function ds_dj(x, y, GP) real, intent(in) :: x !< The longitude in question real, intent(in) :: y !< The latitude in question type(GPS), intent(in) :: GP !< A structure of grid parameters - ! Local variables - real :: ds_dj - ds_dj = GP%Rad_Earth * dy_dj(y,GP) + real :: ds_dj ! The returned grid spacing [L ~> m] + + ds_dj = GP%Rad_Earth_L * dy_dj(y,GP) ! In general, this might be... - ! ds_dj = GP%Rad_Earth * sqrt( cos(y)*cos(y) * dx_dj(x,y,GP)*dx_dj(x,y,GP) + & + ! ds_dj = GP%Rad_Earth_L * sqrt( cos(y)*cos(y) * dx_dj(x,y,GP)*dx_dj(x,y,GP) + & ! dy_dj(x,y,GP)*dy_dj(x,y,GP)) end function ds_dj @@ -1199,8 +1214,7 @@ subroutine initialize_masks(G, PF, US) type(param_file_type), intent(in) :: PF !< Parameter file structure type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type ! Local variables - real :: m_to_Z_scale ! A unit conversion factor from m to Z. - real :: m_to_L ! A unit conversion factor [L m-1 ~> nondim] + 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]. @@ -1209,7 +1223,6 @@ subroutine initialize_masks(G, PF, US) 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 - m_to_L = 1.0 ; if (present(US)) m_to_L = US%m_to_L call get_param(PF, mdl, "MINIMUM_DEPTH", min_depth, & "If MASKING_DEPTH is unspecified, then anything shallower than "//& diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90 index c252e296a5..49a9bfa6f2 100644 --- a/src/initialization/MOM_shared_initialization.F90 +++ b/src/initialization/MOM_shared_initialization.F90 @@ -95,7 +95,7 @@ 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 ~> nondim] + 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 @@ -184,7 +184,8 @@ end subroutine initialize_topography_from_file subroutine apply_topography_edits_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(inout) :: D !< Ocean bottom depth in m or Z if US is present + 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 @@ -298,13 +299,13 @@ subroutine initialize_topography_named(D, G, param_file, topog_config, max_depth ! 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. + 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. - real :: expdecay ! A decay scale of associated with the sloping boundaries [m]. - real :: Dedge ! The depth [Z ~> m], at the basin edge -! real :: south_lat, west_lon, len_lon, len_lat, Rad_earth + real :: D0 ! A constant to make the maximum basin depth MAXIMUM_DEPTH [Z ~> m] + real :: expdecay ! A decay scale of associated with the sloping boundaries [L ~> m] + real :: Dedge ! The depth at the basin edge [Z ~> m] integer :: i, j, is, ie, js, je, isd, ied, jsd, jed character(len=40) :: mdl = "initialize_topography_named" ! This subroutine's name. is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec @@ -315,6 +316,7 @@ subroutine initialize_topography_named(D, G, param_file, topog_config, max_depth "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) @@ -325,23 +327,9 @@ subroutine initialize_topography_named(D, G, param_file, topog_config, max_depth 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) -! call get_param(param_file, mdl, "SOUTHLAT", south_lat, & -! "The southern latitude of the domain.", units="degrees", & -! fail_if_missing=.true.) -! call get_param(param_file, mdl, "LENLAT", len_lat, & -! "The latitudinal length of the domain.", units="degrees", & -! fail_if_missing=.true.) -! call get_param(param_file, mdl, "WESTLON", west_lon, & -! "The western longitude of the domain.", units="degrees", & -! default=0.0) -! call get_param(param_file, mdl, "LENLON", len_lon, & -! "The longitudinal length of the domain.", units="degrees", & -! fail_if_missing=.true.) -! call get_param(param_file, mdl, "RAD_EARTH", Rad_Earth, & -! "The radius of the Earth.", units="m", default=6.378e6) 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) + "the named topographies.", units="m", default=400000.0, scale=m_to_L) endif @@ -351,30 +339,30 @@ subroutine initialize_topography_named(D, G, param_file, topog_config, max_depth do i=is,ie ; do j=js,je ; D(i,j) = max_depth ; enddo ; enddo elseif (trim(topog_config) == "spoon") then D0 = (max_depth - Dedge) / & - ((1.0 - exp(-0.5*G%len_lat*G%Rad_earth*PI/(180.0 *expdecay))) * & - (1.0 - exp(-0.5*G%len_lat*G%Rad_earth*PI/(180.0 *expdecay)))) + ((1.0 - exp(-0.5*G%len_lat*G%Rad_Earth_L*PI/(180.0 *expdecay))) * & + (1.0 - exp(-0.5*G%len_lat*G%Rad_Earth_L*PI/(180.0 *expdecay)))) do i=is,ie ; do j=js,je ! This sets a bowl shaped (sort of) bottom topography, with a ! ! maximum depth of max_depth. ! D(i,j) = Dedge + D0 * & (sin(PI * (G%geoLonT(i,j) - (G%west_lon)) / G%len_lon) * & - (1.0 - exp((G%geoLatT(i,j) - (G%south_lat+G%len_lat))*G%Rad_earth*PI / & + (1.0 - exp((G%geoLatT(i,j) - (G%south_lat+G%len_lat))*G%Rad_Earth_L*PI / & (180.0*expdecay)) )) enddo ; enddo elseif (trim(topog_config) == "bowl") then D0 = (max_depth - Dedge) / & - ((1.0 - exp(-0.5*G%len_lat*G%Rad_earth*PI/(180.0 *expdecay))) * & - (1.0 - exp(-0.5*G%len_lat*G%Rad_earth*PI/(180.0 *expdecay)))) + ((1.0 - exp(-0.5*G%len_lat*G%Rad_Earth_L*PI/(180.0 *expdecay))) * & + (1.0 - exp(-0.5*G%len_lat*G%Rad_Earth_L*PI/(180.0 *expdecay)))) ! This sets a bowl shaped (sort of) bottom topography, with a ! maximum depth of max_depth. do i=is,ie ; do j=js,je D(i,j) = Dedge + D0 * & (sin(PI * (G%geoLonT(i,j) - G%west_lon) / G%len_lon) * & - ((1.0 - exp(-(G%geoLatT(i,j) - G%south_lat)*G%Rad_Earth*PI/ & + ((1.0 - exp(-(G%geoLatT(i,j) - G%south_lat)*G%Rad_Earth_L*PI/ & (180.0*expdecay))) * & (1.0 - exp((G%geoLatT(i,j) - (G%south_lat+G%len_lat))* & - G%Rad_Earth*PI/(180.0*expdecay))))) + G%Rad_Earth_L*PI/(180.0*expdecay))))) enddo ; enddo elseif (trim(topog_config) == "halfpipe") then D0 = max_depth - Dedge @@ -510,10 +498,13 @@ subroutine set_rotation_beta_plane(f, G, param_file, US) ! This subroutine sets up the Coriolis parameter for a beta-plane integer :: I, J 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 m-1 ~> s-1 m-1] - real :: beta_lat_ref ! The reference latitude for the beta plane [degrees/km/m/cm] - real :: y_scl, Rad_Earth - real :: T_to_s ! A time unit conversion factor + 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/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 @@ -522,28 +513,30 @@ 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) 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) + "the betaplane option.", units="m-1 s-1", default=0.0, scale=T_to_s*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, & - "The radius of the Earth.", units="m", default=6.378e6) + 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) beta_lat_ref_units = "degrees" - y_scl = PI * Rad_Earth/ 180. + y_scl = PI * Rad_Earth_L / 180. case ("k") beta_lat_ref_units = "kilometers" - y_scl = 1.E3 + y_scl = 1.E3 * m_to_L case ("m") beta_lat_ref_units = "meters" - y_scl = 1. + y_scl = 1. * m_to_L case default ; call MOM_error(FATAL, & " set_rotation_beta_plane: unknown AXIS_UNITS = "//trim(axis_units)) end select @@ -644,8 +637,8 @@ subroutine reset_face_lengths_named(G, param_file, name, US) ! Local variables character(len=256) :: mesg ! Message for error messages. - real :: m_to_L ! A unit conversion factor [L m-1 ~> nondim] - real :: L_to_m ! A unit conversion factor [m L-1 ~> nondim] + 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 :: pi_180 integer :: option = -1 @@ -772,8 +765,8 @@ subroutine reset_face_lengths_file(G, param_file, US) 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 ~> nondim] - real :: L_to_m ! A unit conversion factor [m L-1 ~> nondim] + 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 @@ -850,8 +843,8 @@ subroutine reset_face_lengths_list(G, param_file, US) integer, allocatable, dimension(:) :: & u_line_no, v_line_no, & ! The line numbers in lines of u- and v-face lines u_line_used, v_line_used ! The number of times each u- and v-line is used. - real :: m_to_L ! A unit conversion factor [L m-1 ~> nondim] - real :: L_to_m ! A unit conversion factor [m L-1 ~> nondim] + 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 :: 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. @@ -1022,7 +1015,7 @@ 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) * m_to_L*min(L_to_m*G%dyCu(I,j), max(u_width(npt), 0.0)) + G%dy_Cu(I,j) = G%mask2dCu(I,j) * min(G%dyCu(I,j), max(m_to_L*u_width(npt), 0.0)) 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 write(stdout,'(A,2F8.2,A,4F8.2,A)') "read_face_lengths_list : G%mask2dCu=0 at ",lat,lon," (",& @@ -1052,7 +1045,7 @@ 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) * m_to_L*min(L_to_m*G%dxCv(i,J), max(v_width(npt), 0.0)) + G%dx_Cv(i,J) = G%mask2dCv(i,J) * min(G%dxCv(i,J), max(m_to_L*v_width(npt), 0.0)) 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 write(stdout,'(A,2F8.2,A,4F8.2,A)') "read_face_lengths_list : G%mask2dCv=0 at ",lat,lon," (",& diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index dfcd097be0..ce1f9ad92f 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -1763,10 +1763,10 @@ subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, depth_t real, allocatable, dimension(:,:,:) :: tmp_tr ! A temporary array for reading sponge fields real, allocatable, dimension(:,:,:) :: tmp_u,tmp_v ! A temporary array for reading sponge fields - real :: Idamp(SZI_(G),SZJ_(G)) ! The inverse damping rate [T-1 ~> s-1]. - real :: Idamp_u(SZIB_(G),SZJ_(G)) ! The inverse damping rate for velocity fields [T-1 ~> s-1]. - real :: Idamp_v(SZI_(G),SZJB_(G)) ! The inverse damping rate for velocity fields [T-1 ~> s-1]. - real :: pres(SZI_(G)) ! An array of the reference pressure [R L2 T-2 ~> Pa]. + real :: Idamp(SZI_(G),SZJ_(G)) ! The sponge damping rate [T-1 ~> s-1] + real :: Idamp_u(SZIB_(G),SZJ_(G)) ! The sponge damping rate for velocity fields [T-1 ~> s-1] + real :: Idamp_v(SZI_(G),SZJB_(G)) ! The sponge damping rate for velocity fields [T-1 ~> s-1] + real :: pres(SZI_(G)) ! An array of the reference pressure [R L2 T-2 ~> Pa] integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state integer :: i, j, k, is, ie, js, je, nz @@ -1852,7 +1852,7 @@ subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, depth_t "of sponge restoring data.", default=time_space_interp_sponge) - ! Read in inverse damping rate for tracers + ! Read in sponge damping rate for tracers filename = trim(inputdir)//trim(damping_file) call log_param(param_file, mdl, "INPUTDIR/SPONGE_DAMPING_FILE", filename) if (.not.file_exists(filename, G%Domain)) & @@ -1863,7 +1863,7 @@ subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, depth_t call MOM_read_data(filename, Idamp_var, Idamp(:,:), G%Domain, scale=US%T_to_s) - ! Read in inverse damping rate for velocities + ! Read in sponge damping rate for velocities if (sponge_uv) then if (separate_idamp_for_uv()) then filename = trim(inputdir)//trim(uv_damping_file) @@ -1911,7 +1911,7 @@ subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, depth_t if (eta(i,j,K) < (eta(i,j,K+1) + GV%Angstrom_Z)) & eta(i,j,K) = eta(i,j,K+1) + GV%Angstrom_Z enddo ; enddo ; enddo - ! Set the inverse damping rates so that the model will know where to + ! Set the sponge damping rates so that the model will know where to ! apply the sponges, along with the interface heights. call initialize_sponge(Idamp, eta, G, param_file, Layer_CSp, GV) deallocate(eta) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index e4f18e75d7..c786f395cf 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -29,6 +29,7 @@ module MOM_MEKE !> Control structure that contains MEKE parameters and diagnostics handles type, public :: MEKE_CS ; private + logical :: initialized = .false. !< True if this control structure has been initialized. ! Parameters real :: MEKE_FrCoeff !< Efficiency of conversion of ME into MEKE [nondim] real :: MEKE_GMcoeff !< Efficiency of conversion of PE into MEKE [nondim] @@ -79,7 +80,8 @@ module MOM_MEKE real :: MEKE_advection_factor !< A scaling in front of the advection of MEKE [nondim] real :: MEKE_topographic_beta !< Weight for how much topographic beta is considered !! when computing beta in Rhines scale [nondim] - real :: MEKE_restoring_rate !< Inverse of the timescale used to nudge MEKE toward its equilibrium value [s-1]. + real :: MEKE_restoring_rate !< Inverse of the timescale used to nudge MEKE toward its + !! equilibrium value [T-1 ~> s-1]. logical :: MEKE_advection_bug !< If true, recover a bug in the calculation of the barotropic !! transport for the advection of MEKE, wherein only the transports in the !! deepest layer are used. @@ -175,6 +177,9 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB + if (.not.CS%initialized) call MOM_error(FATAL, & + "MOM_MEKE: Module must be initialized before it is used.") + if ((CS%MEKE_Cd_scale > 0.0) .or. (CS%MEKE_Cb>0.) .or. CS%visc_drag) then use_drag_rate = .true. else @@ -1048,6 +1053,7 @@ logical function MEKE_init(Time, G, US, param_file, diag, CS, MEKE, restart_CS) "a sub-grid mesoscale eddy kinetic energy budget.", & default=.false.) if (.not. MEKE_init) return + CS%initialized = .true. call MOM_mesg("MEKE_init: reading parameters ", 5) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 1a417ea650..4381af9d84 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -30,6 +30,7 @@ module MOM_hor_visc !> Control structure for horizontal viscosity type, public :: hor_visc_CS ; private + logical :: initialized = .false. !< True if this control structure has been initialized. logical :: Laplacian !< Use a Laplacian horizontal viscosity if true. logical :: biharmonic !< Use a biharmonic horizontal viscosity if true. logical :: debug !< If true, write verbose checksums for debugging purposes. @@ -415,6 +416,9 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, apply_OBC = .true. endif ; endif ; endif + if (.not.CS%initialized) call MOM_error(FATAL, & + "MOM_hor_visc: Module must be initialized before it is used.") + if (.not.(CS%Laplacian .or. CS%biharmonic)) return find_FrictWork = (CS%id_FrictWork > 0) @@ -1815,6 +1819,8 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB + CS%initialized = .true. + CS%diag => diag ! Read parameters and write them to the model log. call log_version(param_file, mdl, version, "") diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 5902f98b56..da8e936642 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -25,6 +25,7 @@ module MOM_lateral_mixing_coeffs !> Variable mixing coefficients type, public :: VarMix_CS + logical :: initialized = .false. !< True if this control structure has been initialized. logical :: use_variable_mixing !< If true, use the variable mixing. logical :: Resoln_scaling_used !< If true, a resolution function is used somewhere to scale !! away one of the viscosities or diffusivities when the @@ -175,6 +176,9 @@ subroutine calc_depth_function(G, CS) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB + if (.not. CS%initialized) call MOM_error(FATAL, "calc_depth_function: "// & + "Module must be initialized before it is used.") + if (.not. CS%calculate_depth_fns) return if (.not. allocated(CS%Depth_fn_u)) call MOM_error(FATAL, & "calc_depth_function: %Depth_fn_u is not associated with Depth_scaled_KhTh.") @@ -216,6 +220,9 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB + if (.not. CS%initialized) call MOM_error(FATAL, "calc_resoln_function: "// & + "Module must be initialized before it is used.") + if (CS%calculate_cg1) then if (.not. allocated(CS%cg1)) call MOM_error(FATAL, & "calc_resoln_function: %cg1 is not associated with Resoln_scaled_Kh.") @@ -458,6 +465,9 @@ subroutine calc_slope_functions(h, tv, dt, G, GV, US, CS, OBC) real, dimension(SZIB_(G), SZJ_(G),SZK_(GV)+1) :: dzSxN ! |Sx| N times dz at u-points [Z T-1 ~> m s-1] real, dimension(SZI_(G), SZJB_(G),SZK_(GV)+1) :: dzSyN ! |Sy| N times dz at v-points [Z T-1 ~> m s-1] + if (.not. CS%initialized) call MOM_error(FATAL, "MOM_lateral_mixing_coeffs.F90, calc_slope_functions: "//& + "Module must be initialized before it is used.") + if (CS%calculate_Eady_growth_rate) then if (CS%use_simpler_Eady_growth_rate) then call find_eta(h, tv, G, GV, US, e, halo_size=2) @@ -524,6 +534,9 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C real :: S2_v(SZI_(G), SZJB_(G)) logical :: local_open_u_BC, local_open_v_BC + if (.not. CS%initialized) call MOM_error(FATAL, "calc_Visbeck_coeffs_old: "// & + "Module must be initialized before it is used.") + if (.not. CS%calculate_Eady_growth_rate) return if (.not. allocated(CS%SN_u)) call MOM_error(FATAL, "calc_slope_function:"// & "%SN_u is not associated with use_variable_mixing.") @@ -682,7 +695,7 @@ subroutine calc_Eady_growth_rate_2D(CS, G, GV, US, OBC, h, e, dzu, dzv, dzSxN, d real :: vint_SN(SZIB_(G)) ! Cumulative integral of SN [Z T-1 ~> m s-1] real, dimension(SZIB_(G),SZJ_(G)) :: SN_cpy !< SN at u-points [T-1 ~> s-1] real :: dz_neglect ! An incy wincy distance to avoid division by zero [Z ~> m] - real :: r_crp_dist ! The inverse of the distance over which to scale the cropping [Z-1 ~. m-1] + real :: r_crp_dist ! The inverse of the distance over which to scale the cropping [Z-1 ~> m-1] real :: dB, dT ! Elevation variables used when cropping [Z ~> m] integer :: i, j, k, l_seg logical :: local_open_u_BC, local_open_v_BC, crop @@ -869,6 +882,9 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slop real :: S2N2_v_local(SZI_(G), SZJB_(G),SZK_(GV)) logical :: local_open_u_BC, local_open_v_BC + if (.not. CS%initialized) call MOM_error(FATAL, "calc_slope_functions_using_just_e: "// & + "Module must be initialized before it is used.") + if (.not. CS%calculate_Eady_growth_rate) return if (.not. allocated(CS%SN_u)) call MOM_error(FATAL, "calc_slope_function:"// & "%SN_u is not associated with use_variable_mixing.") @@ -1168,6 +1184,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, 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 + CS%initialized = .true. in_use = .false. ! Set to true to avoid deallocating CS%diag => diag ! Diagnostics pointer CS%calculate_cg1 = .false. diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 index bb37245c5b..04982d7171 100644 --- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 +++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 @@ -36,6 +36,7 @@ module MOM_mixed_layer_restrat !> Control structure for mom_mixed_layer_restrat type, public :: mixedlayer_restrat_CS ; private + logical :: initialized = .false. !< True if this control structure has been initialized. real :: ml_restrat_coef !< A non-dimensional factor by which the instability is enhanced !! over what would be predicted based on the resolved gradients !! [nondim]. This increases with grid spacing^2, up to something @@ -104,6 +105,9 @@ subroutine mixedlayer_restrat(h, uhtr, vhtr, tv, forces, dt, MLD, VarMix, G, GV, type(VarMix_CS), intent(in) :: VarMix !< Variable mixing control struct type(mixedlayer_restrat_CS), intent(inout) :: CS !< Module control structure + if (.not. CS%initialized) call MOM_error(FATAL, "MOM_mixedlayer_restrat: "// & + "Module must be initialized before it is used.") + if (GV%nkml>0) then call mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS) else @@ -610,6 +614,9 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nkml = GV%nkml + if (.not. CS%initialized) call MOM_error(FATAL, "MOM_mixedlayer_restrat: "// & + "Module must be initialized before it is used.") + if ((nkml<2) .or. (CS%ml_restrat_coef<=0.0)) return uDml(:) = 0.0 ; vDml(:) = 0.0 @@ -816,6 +823,8 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS, "BULKMIXEDLAYER is true.", default=.false.) if (.not. mixedlayer_restrat_init) return + CS%initialized = .true. + ! Nonsense values to cause problems when these parameters are not used CS%MLE_MLD_decay_time = -9.e9*US%s_to_T CS%MLE_density_diff = -9.e9*US%kg_m3_to_R diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index 7382459c16..8303d30621 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -35,6 +35,7 @@ module MOM_thickness_diffuse !> Control structure for thickness diffusion type, public :: thickness_diffuse_CS ; private + logical :: initialized = .false. !< True if this control structure has been initialized. real :: Khth !< Background interface depth diffusivity [L2 T-1 ~> m2 s-1] real :: Khth_Slope_Cff !< Slope dependence coefficient of Khth [nondim] real :: max_Khth_CFL !< Maximum value of the diffusive CFL for thickness diffusion @@ -160,6 +161,9 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp real :: KH_u_lay(SZI_(G), SZJ_(G)) ! layer ave thickness diffusivities [L2 T-1 ~> m2 s-1] real :: KH_v_lay(SZI_(G), SZJ_(G)) ! layer ave thickness diffusivities [L2 T-1 ~> m2 s-1] + if (.not. CS%initialized) call MOM_error(FATAL, "MOM_thickness_diffuse: "//& + "Module must be initialized before it is used.") + if ((.not.CS%thickness_diffuse) & .or. .not. (CS%Khth > 0.0 .or. VarMix%use_variable_mixing)) return @@ -1897,6 +1901,7 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS) ! rotation [nondim]. logical :: default_2018_answers ! The default setting for the various 2018_ANSWERS flags. + CS%initialized = .true. CS%diag => diag ! Read all relevant parameters and write them to the model log. diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index b9ceb85cc5..5342c621dd 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -641,7 +641,7 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, buoyFlux, Kt, Ks, Kv, & real :: surfFricVel, surfBuoyFlux real :: sigma, sigmaRatio - real :: buoy_scale ! A unit conversion factor for buoyancy fluxes [m2 T3 L-2 s-3 ~> nondim] + real :: buoy_scale ! A unit conversion factor for buoyancy fluxes [m2 T3 L-2 s-3 ~> 1] real :: dh ! The local thickness used for calculating interface positions [m] real :: hcorr ! A cumulative correction arising from inflation of vanished layers [m] @@ -948,7 +948,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl real :: zBottomMinusOffset ! Height of bottom plus a little bit [m] real :: SLdepth_0d ! Surface layer depth = surf_layer_ext*OBLdepth. real :: hTot ! Running sum of thickness used in the surface layer average [m] - real :: buoy_scale ! A unit conversion factor for buoyancy fluxes [m2 T3 L-2 s-3 ~> nondim] + real :: buoy_scale ! A unit conversion factor for buoyancy fluxes [m2 T3 L-2 s-3 ~> 1] real :: delH ! Thickness of a layer [m] real :: surfHtemp, surfTemp ! Integral and average of temp over the surface layer real :: surfHsalt, surfSalt ! Integral and average of saln over the surface layer diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 index ca545c14ad..c80ee0ea61 100644 --- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 +++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 @@ -30,6 +30,7 @@ module MOM_bulk_mixed_layer !> The control structure with parameters for the MOM_bulk_mixed_layer module type, public :: bulkmixedlayer_CS ; private + logical :: initialized = .false. !< True if this control structure has been initialized. integer :: nkml !< The number of layers in the mixed layer. integer :: nkbl !< The number of buffer layers. integer :: nsw !< The number of bands of penetrating shortwave radiation. @@ -186,9 +187,9 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C !! layer; this should be increased due to !! mixed layer entrainment [H ~> m or kg m-2]. type(bulkmixedlayer_CS), intent(inout) :: CS !< Bulk mixed layer control struct - type(optics_type), pointer :: optics !< The structure containing the inverse of the - !! vertical absorption decay scale for - !! penetrating shortwave radiation [m-1]. + type(optics_type), pointer :: optics !< The structure that can be queried for the + !! inverse of the vertical absorption decay + !! scale for penetrating shortwave radiation. real, dimension(:,:), pointer :: Hml !< Active mixed layer depth [Z ~> m]. logical, intent(in) :: aggregate_FW_forcing !< If true, the net incoming and !! outgoing surface freshwater fluxes are @@ -328,6 +329,8 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + if (.not. CS%initialized) call MOM_error(FATAL, "MOM_bulk_mixed_layer: "//& + "Module must be initialized before it is used.") if (GV%nkml < 1) return if (.not. associated(tv%eqn_of_state)) call MOM_error(FATAL, & @@ -474,9 +477,9 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C ! Here we add an additional source of TKE to the mixed layer where river ! is present to simulate unresolved estuaries. The TKE input is diagnosed ! as follows: - ! TKE_river[m3 s-3] = 0.5*rivermix_depth*g*Irho0*drho_ds* + ! TKE_river[Z L2 T-3 ~> m3 s-3] = 0.5*rivermix_depth * g * Irho0**2 * drho_ds * ! River*(Samb - Sriver) = CS%mstar*U_star^3 - ! where River is in units of [m s-1]. + ! where River is in units of [R Z T-1 ~> kg m-2 s-1]. ! Samb = Ambient salinity at the mouth of the estuary ! rivermix_depth = The prescribed depth over which to mix river inflow ! drho_ds = The gradient of density wrt salt at the ambient surface salinity. @@ -3362,6 +3365,7 @@ subroutine bulkmixedlayer_init(Time, G, GV, US, param_file, diag, CS) logical :: use_temperature, use_omega isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + CS%initialized = .true. CS%diag => diag CS%Time => Time diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 9aa953ed06..bac311bd6d 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -88,6 +88,7 @@ module MOM_diabatic_driver !> Control structure for this module type, public :: diabatic_CS ; private + logical :: initialized = .false. !< True if this control structure has been initialized. logical :: use_legacy_diabatic !< If true (default), use a legacy version of the diabatic !! algorithm. This is temporary and is needed to avoid change @@ -305,6 +306,10 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & if (.not. associated(CS)) call MOM_error(FATAL, "MOM_diabatic_driver: "// & "Module must be initialized before it is used.") + + if (.not. CS%initialized) call MOM_error(FATAL, "MOM_diabatic_driver: "// & + "Module must be initialized before it is used.") + if (dt == 0.0) call MOM_error(FATAL, "MOM_diabatic_driver: "// & "diabatic was called with a zero length timestep.") if (dt < 0.0) call MOM_error(FATAL, "MOM_diabatic_driver: "// & @@ -2908,6 +2913,8 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di allocate(CS) endif + CS%initialized = .true. + CS%diag => diag CS%Time => Time diff --git a/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 b/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 index 2cca587e9e..379275196e 100644 --- a/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 +++ b/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 @@ -68,7 +68,7 @@ subroutine diapyc_energy_req_test(h_3d, dt, tv, G, GV, US, CS, Kd_int) Kd, & ! A column of diapycnal diffusivities at interfaces [Z2 T-1 ~> m2 s-1]. h_top, h_bot ! Distances from the top or bottom [H ~> m or kg m-2]. real :: ustar, absf, htot - real :: energy_Kd ! The energy used by diapycnal mixing [W m-2]. + real :: energy_Kd ! The energy used by diapycnal mixing [R Z L2 T-3 ~> W m-2]. real :: tmp1 ! A temporary array. integer :: i, j, k, is, ie, js, je, nz, itt logical :: may_print @@ -77,6 +77,9 @@ subroutine diapyc_energy_req_test(h_3d, dt, tv, G, GV, US, CS, Kd_int) if (.not. associated(CS)) call MOM_error(FATAL, "diapyc_energy_req_test: "// & "Module must be initialized before it is used.") + if (.not. CS%initialized) call MOM_error(FATAL, "diapyc_energy_req_test: "// & + "Module must be initialized before it is used.") + !$OMP do do j=js,je ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then if (present(Kd_int) .and. .not.CS%use_test_Kh_profile) then @@ -1193,7 +1196,7 @@ subroutine find_PE_chg_orig(Kddt_h, h_k, b_den_1, dTe_term, dSe_term, & real :: b1Kd ! Temporary array [nondim] real :: ColHt_chg ! The change in column thickness [Z ~> m]. real :: dColHt_max ! The change in column thickness for infinite diffusivity [Z ~> m]. - real :: dColHt_dKd ! The partial derivative of column thickness with Kddt_h [Z H-1 ~> 1 or m3 kg-1]. + real :: dColHt_dKd ! The partial derivative of column thickness with Kddt_h [Z H-1 ~> nondim or m3 kg-1] real :: dT_k, dT_km1 ! Temperature changes in layers k and k-1 [degC] real :: dS_k, dS_km1 ! Salinity changes in layers k and k-1 [ppt] real :: I_Kr_denom ! Temporary array [H-2 ~> m-2 or m4 kg-2] diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 4a762cd34c..25b06d51fa 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -32,6 +32,7 @@ module MOM_energetic_PBL !> This control structure holds parameters for the MOM_energetic_PBL module type, public :: energetic_PBL_CS ; private + logical :: initialized = .false. !< True if this control structure has been initialized. !/ Constants real :: VonKar = 0.41 !< The von Karman coefficient. This should be a runtime parameter, @@ -344,6 +345,8 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + if (.not. CS%initialized) call MOM_error(FATAL, "energetic_PBL: "//& + "Module must be initialized before it is used.") if (.not. associated(tv%eqn_of_state)) call MOM_error(FATAL, & "energetic_PBL: Temperature, salinity and an equation of state "//& "must now be used.") @@ -636,7 +639,7 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs real :: LA ! The value of the Langmuir number [nondim] real :: LAmod ! The modified Langmuir number by convection [nondim] real :: hbs_here ! The local minimum of hb_hs and MixLen_shape, times a - ! conversion factor from H to Z [Z H-1 ~> 1 or m3 kg-1]. + ! conversion factor from H to Z [Z H-1 ~> nondim or m3 kg-1]. real :: nstar_FC ! The fraction of conv_PErel that can be converted to mixing [nondim]. real :: TKE_reduc ! The fraction by which TKE and other energy fields are ! reduced to support mixing [nondim]. between 0 and 1. @@ -1634,7 +1637,7 @@ subroutine find_PE_chg_orig(Kddt_h, h_k, b_den_1, dTe_term, dSe_term, & real :: b1Kd ! Temporary array [nondim] real :: ColHt_chg ! The change in column thickness [Z ~> m]. real :: dColHt_max ! The change in column thickness for infinite diffusivity [Z ~> m]. - real :: dColHt_dKd ! The partial derivative of column thickness with Kddt_h [Z H-1 ~> 1 or m3 kg-1]. + real :: dColHt_dKd ! The partial derivative of column thickness with Kddt_h [Z H-1 ~> nondim or m3 kg-1] real :: dT_k, dT_km1 ! Temperature changes in layers k and k-1 [degC] real :: dS_k, dS_km1 ! Salinity changes in layers k and k-1 [ppt] real :: I_Kr_denom ! Temporary array [H-2 ~> m-2 or m4 kg-2] @@ -1924,6 +1927,7 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) logical :: use_la_windsea isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + CS%initialized = .true. CS%diag => diag CS%Time => Time diff --git a/src/parameterizations/vertical/MOM_entrain_diffusive.F90 b/src/parameterizations/vertical/MOM_entrain_diffusive.F90 index e279092fef..fdc38ebf1e 100644 --- a/src/parameterizations/vertical/MOM_entrain_diffusive.F90 +++ b/src/parameterizations/vertical/MOM_entrain_diffusive.F90 @@ -27,6 +27,7 @@ module MOM_entrain_diffusive !> The control structure holding parametes for the MOM_entrain_diffusive module type, public :: entrain_diffusive_CS ; private + logical :: initialized = .false. !< True if this control structure has been initialized. logical :: bulkmixedlayer !< If true, a refined bulk mixed layer is used with !! GV%nk_rho_varies variable density mixed & buffer layers. integer :: max_ent_it !< The maximum number of iterations that may be used to @@ -207,6 +208,9 @@ subroutine entrainment_diffusive(h, tv, fluxes, dt, G, GV, US, CS, ea, eb, & Angstrom = GV%Angstrom_H h_neglect = GV%H_subroundoff + if (.not. CS%initialized) call MOM_error(FATAL, & + "MOM_entrain_diffusive: Module must be initialized before it is used.") + if (.not.(present(Kd_Lay) .or. present(Kd_int))) call MOM_error(FATAL, & "MOM_entrain_diffusive: Either Kd_Lay or Kd_int must be present in call.") @@ -2076,6 +2080,7 @@ subroutine entrain_diffusive_init(Time, G, GV, US, param_file, diag, CS, just_re # include "version_variable.h" character(len=40) :: mdl = "MOM_entrain_diffusive" ! This module's name. + CS%initialized = .true. CS%diag => diag CS%bulkmixedlayer = (GV%nkml > 0) diff --git a/src/parameterizations/vertical/MOM_geothermal.F90 b/src/parameterizations/vertical/MOM_geothermal.F90 index 7fdaa6abda..9884eec7a2 100644 --- a/src/parameterizations/vertical/MOM_geothermal.F90 +++ b/src/parameterizations/vertical/MOM_geothermal.F90 @@ -23,6 +23,7 @@ module MOM_geothermal !> Control structure for geothermal heating type, public :: geothermal_CS ; private + logical :: initialized = .false. !< True if this control structure has been initialized. real :: dRcv_dT_inplace !< The value of dRcv_dT above which (dRcv_dT is negative) the !! water is heated in place instead of moving upward between !! layers in non-ALE layered mode [R degC-1 ~> kg m-3 degC-1] @@ -119,6 +120,8 @@ subroutine geothermal_entraining(h, tv, dt, ea, eb, G, GV, US, CS, halo) is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo endif + if (.not. CS%initialized) call MOM_error(FATAL, "MOM_geothermal: "//& + "Module must be initialized before it is used.") if (.not.CS%apply_geothermal) return nkmb = GV%nk_rho_varies @@ -392,6 +395,8 @@ subroutine geothermal_in_place(h, tv, dt, G, GV, US, CS, halo) is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo endif + if (.not. CS%initialized) call MOM_error(FATAL, "MOM_geothermal: "//& + "Module must be initialized before it is used.") if (.not.CS%apply_geothermal) return Irho_cp = 1.0 / (GV%H_to_RZ * tv%C_p) @@ -502,10 +507,11 @@ subroutine geothermal_init(Time, G, GV, US, param_file, diag, CS, useALEalgorith ! Local variables character(len=200) :: inputdir, geo_file, filename, geotherm_var real :: geo_scale ! A constant heat flux or dimensionally rescaled geothermal flux scaling factor - ! [Q R Z T-1 ~> W m-2] or [Q R Z m2 s J-1 T-1 ~> 1] + ! [Q R Z T-1 ~> W m-2] or [Q R Z m2 s J-1 T-1 ~> nondim] integer :: i, j, isd, ied, jsd, jed, id isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + CS%initialized = .true. CS%diag => diag CS%Time => Time diff --git a/src/parameterizations/vertical/MOM_internal_tide_input.F90 b/src/parameterizations/vertical/MOM_internal_tide_input.F90 index df24d3f4e9..300cdcbe1e 100644 --- a/src/parameterizations/vertical/MOM_internal_tide_input.F90 +++ b/src/parameterizations/vertical/MOM_internal_tide_input.F90 @@ -34,6 +34,7 @@ module MOM_int_tide_input !> This control structure holds parameters that regulate internal tide energy inputs. type, public :: int_tide_input_CS ; private + logical :: initialized = .false. !< True if this control structure has been initialized. logical :: debug !< If true, write verbose checksums for debugging. type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to !! regulate the timing of diagnostic output. @@ -111,6 +112,9 @@ subroutine set_int_tide_input(u, v, h, tv, fluxes, itide, dt, G, GV, US, CS) if (.not.associated(CS)) call MOM_error(FATAL,"set_diffusivity: "//& "Module must be initialized before it is used.") + if (.not.CS%initialized) call MOM_error(FATAL,"set_diffusivity: "//& + "Module must be initialized before it is used.") + use_EOS = associated(tv%eqn_of_state) ! Smooth the properties through massless layers. @@ -330,6 +334,8 @@ subroutine int_tide_input_init(Time, G, GV, US, param_file, diag, CS, itide) allocate(CS) allocate(itide) + CS%initialized = .true. + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed diff --git a/src/parameterizations/vertical/MOM_opacity.F90 b/src/parameterizations/vertical/MOM_opacity.F90 index e61cc3736b..5dec767f5b 100644 --- a/src/parameterizations/vertical/MOM_opacity.F90 +++ b/src/parameterizations/vertical/MOM_opacity.F90 @@ -552,8 +552,8 @@ subroutine absorbRemainingSW(G, GV, US, h, opacity_band, nsw, optics, j, dt, H_l real, dimension(SZI_(G)), optional, intent(in) :: htot !< Total mixed layer thickness [H ~> m or kg m-2]. real, dimension(SZI_(G)), optional, intent(inout) :: Ttot !< Depth integrated mixed layer !! temperature [degC H ~> degC m or degC kg m-2] - real, dimension(SZI_(G),SZK_(GV)), optional, intent(in) :: dSV_dT !< The partial derivative of specific - !! volume with temperature [R-1 degC-1]. + real, dimension(SZI_(G),SZK_(GV)), optional, intent(in) :: dSV_dT !< The partial derivative of specific volume + !! with temperature [R-1 degC-1 ~> m3 kg-1 degC-1] real, dimension(SZI_(G),SZK_(GV)), optional, intent(inout) :: TKE !< The TKE sink from mixing the heating !! throughout a layer [R Z3 T-2 ~> J m-2]. @@ -935,7 +935,7 @@ subroutine opacity_init(Time, G, GV, US, param_file, diag, CS, optics) ! This include declares and sets the variable "version". # include "version_variable.h" real :: PenSW_absorb_minthick ! A thickness that is used to absorb the remaining shortwave heat - ! flux when that flux drops below PEN_SW_FLUX_ABSORB [m]. + ! flux when that flux drops below PEN_SW_FLUX_ABSORB [H ~> m or kg m-2] real :: PenSW_minthick_dflt ! The default for PenSW_absorb_minthick [m] logical :: default_2018_answers logical :: use_scheme diff --git a/src/parameterizations/vertical/MOM_regularize_layers.F90 b/src/parameterizations/vertical/MOM_regularize_layers.F90 index f42b1ae7ee..868ff6a832 100644 --- a/src/parameterizations/vertical/MOM_regularize_layers.F90 +++ b/src/parameterizations/vertical/MOM_regularize_layers.F90 @@ -23,6 +23,7 @@ module MOM_regularize_layers !> This control structure holds parameters used by the MOM_regularize_layers module type, public :: regularize_layers_CS ; private + logical :: initialized = .false. !< True if this control structure has been initialized. logical :: regularize_surface_layers !< If true, vertically restructure the !! near-surface layers when they have too much !! lateral variations to allow for sensible lateral @@ -93,6 +94,9 @@ subroutine regularize_layers(h, tv, dt, ea, eb, G, GV, US, CS) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + if (.not. CS%initialized) call MOM_error(FATAL, "MOM_regularize_layers: "//& + "Module must be initialized before it is used.") + if (CS%regularize_surface_layers) then call pass_var(h, G%Domain, clock=id_clock_pass) call regularize_surface(h, tv, dt, ea, eb, G, GV, US, CS) @@ -191,6 +195,9 @@ subroutine regularize_surface(h, tv, dt, ea, eb, G, GV, US, CS) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + if (.not. CS%initialized) call MOM_error(FATAL, "MOM_regularize_layers: "//& + "Module must be initialized before it is used.") + if (GV%nkml<1) return nkmb = GV%nk_rho_varies ; nkml = GV%nkml if (.not.associated(tv%eqn_of_state)) call MOM_error(FATAL, & @@ -722,6 +729,8 @@ subroutine regularize_layers_init(Time, G, GV, param_file, diag, CS) integer :: isd, ied, jsd, jed isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + CS%initialized = .true. + CS%diag => diag CS%Time => Time diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index dce302372e..061fe776e1 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -56,6 +56,7 @@ module MOM_set_diffusivity !> This control structure contains parameters for MOM_set_diffusivity. type, public :: set_diffusivity_CS ; private + logical :: initialized = .false. !< True if this control structure has been initialized. logical :: debug !< If true, write verbose checksums for debugging. logical :: bulkmixedlayer !< If true, a refined bulk mixed layer is used with @@ -282,6 +283,9 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i if (.not.associated(CS)) call MOM_error(FATAL,"set_diffusivity: "//& "Module must be initialized before it is used.") + if (.not.CS%initialized) call MOM_error(FATAL,"set_diffusivity: "//& + "Module must be initialized before it is used.") + if (CS%answers_2018) then ! These hard-coded dimensional parameters are being replaced. kappa_dt_fill = US%m_to_Z**2 * 1.e-3 * 7200. @@ -1709,6 +1713,9 @@ subroutine set_BBL_TKE(u, v, h, fluxes, visc, G, GV, US, CS, OBC) if (.not.associated(CS)) call MOM_error(FATAL,"set_BBL_TKE: "//& "Module must be initialized before it is used.") + if (.not.CS%initialized) call MOM_error(FATAL,"set_BBL_TKE: "//& + "Module must be initialized before it is used.") + if (.not.CS%bottomdraglaw .or. (CS%BBL_effic<=0.0)) then if (associated(visc%ustar_BBL)) then do j=js,je ; do i=is,ie ; visc%ustar_BBL(i,j) = 0.0 ; enddo ; enddo @@ -1987,6 +1994,8 @@ subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, int_tide_ endif allocate(CS) + CS%initialized = .true. + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 75fcb04831..83558b9a08 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -43,6 +43,7 @@ module MOM_set_visc !> Control structure for MOM_set_visc type, public :: set_visc_CS ; private + logical :: initialized = .false. !< True if this control structure has been initialized. real :: Hbbl !< The static bottom boundary layer thickness [H ~> m or kg m-2]. !! Runtime parameter `HBBL`. real :: cdrag !< The quadratic drag coefficient. @@ -177,7 +178,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS) ! Rho0 divided by G_Earth and the conversion ! from m to thickness units [H R ~> kg m-2 or kg2 m-5]. real :: cdrag_sqrt_Z ! Square root of the drag coefficient, times a unit conversion - ! factor from lateral lengths to vertical depths [Z L-1 ~> 1]. + ! factor from lateral lengths to vertical depths [Z L-1 ~> nondim]. real :: cdrag_sqrt ! Square root of the drag coefficient [nondim]. real :: oldfn ! The integrated energy required to ! entrain up to the bottom of the layer, @@ -256,7 +257,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS) real :: Cell_width ! The transverse width of the velocity cell [L ~> m]. real :: Rayleigh ! A nondimensional value that is multiplied by the layer's ! velocity magnitude to give the Rayleigh drag velocity, times - ! a lateral to vertical distance conversion factor [Z L-1 ~> 1]. + ! a lateral to vertical distance conversion factor [Z L-1 ~> nondim]. real :: gam ! The ratio of the change in the open interface width ! to the open interface width atop a cell [nondim]. real :: BBL_frac ! The fraction of a layer's drag that goes into the @@ -284,6 +285,9 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS) Vol_quit = 0.9*GV%Angstrom_H + h_neglect C2pi_3 = 8.0*atan(1.0)/3.0 + if (.not.CS%initialized) call MOM_error(FATAL,"MOM_set_viscosity(BBL): "//& + "Module must be initialized before it is used.") + if (.not.CS%bottomdraglaw) return if (CS%debug) then @@ -592,7 +596,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS) ! When stratification dominates h_N< kg m-2 or kg m-5] htot = 0.0 ! Calculate the thickness of a stratification limited BBL ignoring rotation: @@ -947,10 +951,10 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS) ! set kv_bbl to the bound and recompute bbl_thick to be consistent ! but with a ridiculously large upper bound on thickness (for Cd u*=0) kv_bbl = CS%Kv_BBL_min - if (cdrag_sqrt*ustar(i)*BBL_visc_frac*G%Rad_Earth*US%m_to_Z > kv_bbl) then + if (cdrag_sqrt*ustar(i)*BBL_visc_frac*G%Rad_Earth_L*US%L_to_Z > kv_bbl) then bbl_thick_Z = kv_bbl / ( cdrag_sqrt*ustar(i)*BBL_visc_frac ) else - bbl_thick_Z = G%Rad_Earth * US%m_to_Z + bbl_thick_Z = G%Rad_Earth_L * US%L_to_Z endif else kv_bbl = cdrag_sqrt*ustar(i)*bbl_thick_Z*BBL_visc_frac @@ -979,10 +983,10 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS) ! set kv_bbl to the bound and recompute bbl_thick to be consistent ! but with a ridiculously large upper bound on thickness (for Cd u*=0) kv_bbl = CS%Kv_BBL_min - if (cdrag_sqrt*ustar(i)*G%Rad_Earth*US%m_to_Z > kv_bbl) then + if (cdrag_sqrt*ustar(i)*G%Rad_Earth_L*US%L_to_Z > kv_bbl) then bbl_thick_Z = kv_bbl / ( cdrag_sqrt*ustar(i) ) else - bbl_thick_Z = G%Rad_Earth * US%m_to_Z + bbl_thick_Z = G%Rad_Earth_L * US%L_to_Z endif else kv_bbl = cdrag_sqrt*ustar(i)*bbl_thick_Z @@ -1214,7 +1218,7 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS) ! Rho0 divided by G_Earth and the conversion ! from m to thickness units [H R ~> kg m-2 or kg2 m-5]. real :: cdrag_sqrt_Z ! Square root of the drag coefficient, times a unit conversion - ! factor from lateral lengths to vertical depths [Z L-1 ~> 1]. + ! factor from lateral lengths to vertical depths [Z L-1 ~> nondim] real :: cdrag_sqrt ! Square root of the drag coefficient [nondim]. real :: oldfn ! The integrated energy required to ! entrain up to the bottom of the layer, @@ -1245,6 +1249,9 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS) Isq = G%isc-1 ; Ieq = G%IecB ; Jsq = G%jsc-1 ; Jeq = G%JecB nkmb = GV%nk_rho_varies ; nkml = GV%nkml + if (.not.CS%initialized) call MOM_error(FATAL,"MOM_set_viscosity(visc_ML): "//& + "Module must be initialized before it is used.") + if (.not.(CS%dynamic_viscous_ML .or. associated(forces%frac_shelf_u) .or. & associated(forces%frac_shelf_v)) ) return @@ -1927,6 +1934,7 @@ subroutine set_visc_init(Time, G, GV, US, param_file, diag, visc, CS, restart_CS # include "version_variable.h" character(len=40) :: mdl = "MOM_set_visc" ! This module's name. + CS%initialized = .true. CS%OBC => OBC is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index b11bb2d8b2..c8166c47b8 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -747,7 +747,7 @@ subroutine calculate_CVMix_tidal(h, j, N2_int, G, GV, US, CS, Kv, Kd_lay, Kd_int integer :: i, k, is, ie real :: dh, hcorr, Simmons_coeff - real, parameter :: rho_fw = 1000.0 ! fresh water density [kg/m^3] + real, parameter :: rho_fw = 1000.0 ! fresh water density [kg m-3] ! TODO: when coupled, get this from CESM (SHR_CONST_RHOFW) is = G%isc ; ie = G%iec @@ -1001,7 +1001,7 @@ subroutine add_int_tide_diffusivity(h, j, N2_bot, N2_lay, TKE_to_Kd, max_TKE, & ! multiplied by N2_bot/N2_meanz to be coherent with the WKB scaled z ! z*=int(N2/N2_bot) * N2_bot/N2_meanz = int(N2/N2_meanz) ! z0_Polzin_scaled = z0_Polzin * N2_bot/N2_meanz - N2_meanz, & ! vertically averaged squared buoyancy frequency [T-2] for WKB scaling + N2_meanz, & ! vertically averaged squared buoyancy frequency [T-2 ~> s-2] for WKB scaling TKE_itidal_rem, & ! remaining internal tide TKE (from barotropic source) [Z3 T-3 ~> m3 s-3] TKE_Niku_rem, & ! remaining lee-wave TKE [Z3 T-3 ~> m3 s-3] TKE_lowmode_rem, & ! remaining internal tide TKE (from propagating low mode source) [Z3 T-3 ~> m3 s-3] (BDM) @@ -1578,13 +1578,13 @@ subroutine read_tidal_constituents(G, US, tidal_energy_file, CS) tidal_qk1, & ! qk1 coefficient used in Schmittner & Egbert tidal_qo1 ! qo1 coefficient used in Schmittner & Egbert real, allocatable, dimension(:) :: & - z_t, & ! depth from surface to midpoint of input layer [Z] - z_w ! depth from surface to top of input layer [Z] + z_t, & ! depth from surface to midpoint of input layer [Z ~> m] + z_w ! depth from surface to top of input layer [Z ~> m] real, allocatable, dimension(:,:,:) :: & - tc_m2, & ! input lunar semidiurnal tidal energy flux [W/m^2] - tc_s2, & ! input solar semidiurnal tidal energy flux [W/m^2] - tc_k1, & ! input lunar diurnal tidal energy flux [W/m^2] - tc_o1 ! input lunar diurnal tidal energy flux [W/m^2] + tc_m2, & ! input lunar semidiurnal tidal energy flux [W m-2] + tc_s2, & ! input solar semidiurnal tidal energy flux [W m-2] + tc_k1, & ! input lunar diurnal tidal energy flux [W m-2] + tc_o1 ! input lunar diurnal tidal energy flux [W m-2] integer, dimension(4) :: nz_in integer :: k, is, ie, js, je, isd, ied, jsd, jed, i, j diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 878ee771cc..b332346c6c 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -37,6 +37,7 @@ module MOM_vert_friction !> The control structure with parameters and memory for the MOM_vert_friction module type, public :: vertvisc_CS ; private + logical :: initialized = .false. !< True if this control structure has been initialized. real :: Hmix !< The mixed layer thickness in thickness units [H ~> m or kg m-2]. real :: Hmix_stress !< The mixed layer thickness over which the wind !! stress is applied with direct_stress [H ~> m or kg m-2]. @@ -236,6 +237,9 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & if (.not.associated(CS)) call MOM_error(FATAL,"MOM_vert_friction(visc): "// & "Module must be initialized before it is used.") + if (.not.CS%initialized) call MOM_error(FATAL,"MOM_vert_friction(visc): "// & + "Module must be initialized before it is used.") + if (CS%direct_stress) then Hmix = CS%Hmix_stress I_Hmix = 1.0 / Hmix @@ -648,6 +652,9 @@ subroutine vertvisc_remnant(visc, visc_rem_u, visc_rem_v, dt, G, GV, US, CS) if (.not.associated(CS)) call MOM_error(FATAL,"MOM_vert_friction(visc): "// & "Module must be initialized before it is used.") + if (.not.CS%initialized) call MOM_error(FATAL,"MOM_vert_friction(remant): "// & + "Module must be initialized before it is used.") + dt_Z_to_H = dt*GV%Z_to_H do k=1,nz ; do i=Isq,Ieq ; Ray(i,k) = 0.0 ; enddo ; enddo @@ -800,6 +807,9 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC) if (.not.associated(CS)) call MOM_error(FATAL,"MOM_vert_friction(coef): "// & "Module must be initialized before it is used.") + if (.not.CS%initialized) call MOM_error(FATAL,"MOM_vert_friction(coef): "// & + "Module must be initialized before it is used.") + h_neglect = GV%H_subroundoff a_cpl_max = 1.0e37 * US%m_to_Z * US%T_to_s I_Hbbl(:) = 1.0 / (CS%Hbbl + h_neglect) @@ -1702,6 +1712,8 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & endif allocate(CS) + CS%initialized = .true. + if (GV%Boussinesq) then; thickness_units = "m" else; thickness_units = "kg m-2"; endif diff --git a/src/tracer/MOM_tracer_flow_control.F90 b/src/tracer/MOM_tracer_flow_control.F90 index a44750c1fc..4278594913 100644 --- a/src/tracer/MOM_tracer_flow_control.F90 +++ b/src/tracer/MOM_tracer_flow_control.F90 @@ -310,7 +310,7 @@ subroutine tracer_flow_control_init(restart, day, G, GV, US, h, param_file, diag ! Add other user-provided calls here. if (CS%use_USER_tracer_example) & - call USER_initialize_tracer(restart, day, G, GV, h, diag, OBC, CS%USER_tracer_example_CSp, & + call USER_initialize_tracer(restart, day, G, GV, US, h, diag, OBC, CS%USER_tracer_example_CSp, & sponge_CSp) if (CS%use_DOME_tracer) & call initialize_DOME_tracer(restart, day, G, GV, US, h, diag, OBC, CS%DOME_tracer_CSp, & diff --git a/src/tracer/tracer_example.F90 b/src/tracer/tracer_example.F90 index b58e45b366..10551ea247 100644 --- a/src/tracer/tracer_example.F90 +++ b/src/tracer/tracer_example.F90 @@ -134,13 +134,14 @@ end function USER_register_tracer_example !> This subroutine initializes the NTR tracer fields in tr(:,:,:,:) !! and it sets up the tracer output. -subroutine USER_initialize_tracer(restart, day, G, GV, h, diag, OBC, CS, & +subroutine USER_initialize_tracer(restart, day, G, GV, US, h, diag, OBC, CS, & sponge_CSp) logical, intent(in) :: restart !< .true. if the fields have already !! been read from a restart file. type(time_type), target, intent(in) :: day !< Time of the start of the run. type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate @@ -163,7 +164,7 @@ subroutine USER_initialize_tracer(restart, day, G, GV, h, diag, OBC, CS, & real, pointer :: tr_ptr(:,:,:) => NULL() real :: PI ! 3.1415926... calculated as 4*atan(1) real :: tr_y ! Initial zonally uniform tracer concentrations. - real :: dist2 ! The distance squared from a line [m2]. + real :: dist2 ! The distance squared from a line [L2 ~> m2]. integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz, m integer :: IsdB, IedB, JsdB, JedB, lntr @@ -196,9 +197,9 @@ subroutine USER_initialize_tracer(restart, day, G, GV, h, diag, OBC, CS, & ! This sets a stripe of tracer across the basin. PI = 4.0*atan(1.0) do j=js,je - dist2 = (G%Rad_Earth * PI / 180.0)**2 * & + dist2 = (G%Rad_Earth_L * PI / 180.0)**2 * & (G%geoLatT(i,j) - 40.0) * (G%geoLatT(i,j) - 40.0) - tr_y = 0.5*exp(-dist2/(1.0e5*1.0e5)) + tr_y = 0.5 * exp( -dist2 / (1.0e5*US%m_to_L)**2 ) do k=1,nz ; do i=is,ie ! This adds the stripes of tracer to every layer. @@ -282,10 +283,10 @@ subroutine tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, ! Local variables real :: hold0(SZI_(G)) ! The original topmost layer thickness, - ! with surface mass fluxes added back, m. - real :: b1(SZI_(G)) ! b1 and c1 are variables used by the - real :: c1(SZI_(G),SZK_(GV)) ! tridiagonal solver. - real :: d1(SZI_(G)) ! d1=1-c1 is used by the tridiagonal solver. + ! with surface mass fluxes added back [H ~> m or kg m-2]. + real :: b1(SZI_(G)) ! b1 is a variable used by the tridiagonal solver [H ~> m or kg m-2]. + real :: c1(SZI_(G),SZK_(GV)) ! c1 is a variable used by the tridiagonal solver [nondim]. + real :: d1(SZI_(G)) ! d1=1-c1 is used by the tridiagonal solver [nondim]. real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m or kg m-2]. real :: b_denom_1 ! The first term in the denominator of b1 [H ~> m or kg m-2]. @@ -374,7 +375,7 @@ function USER_tracer_stock(h, stocks, G, GV, CS, names, units, stock_index) !! stocks calculated here. ! Local variables - real :: stock_scale ! The dimensional scaling factor to convert stocks to kg [kg H-1 L-2 ~> kg m-3 or nondim] + real :: stock_scale ! The dimensional scaling factor to convert stocks to kg [kg H-1 L-2 ~> kg m-3 or 1] integer :: i, j, k, is, ie, js, je, nz, m is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke diff --git a/src/user/BFB_initialization.F90 b/src/user/BFB_initialization.F90 index 922ae60fc5..8ef21d190f 100644 --- a/src/user/BFB_initialization.F90 +++ b/src/user/BFB_initialization.F90 @@ -91,7 +91,7 @@ subroutine BFB_initialize_sponges_southonly(G, GV, US, use_temperature, tv, dept ! Local variables real :: eta(SZI_(G),SZJ_(G),SZK_(GV)+1) ! A temporary array for eta, in depth units [Z ~> m]. - real :: Idamp(SZI_(G),SZJ_(G)) ! The inverse damping rate [T-1 ~> s-1]. + real :: Idamp(SZI_(G),SZJ_(G)) ! The sponge damping rate [T-1 ~> s-1] real :: H0(SZK_(GV)) ! Resting layer thicknesses in depth units [Z ~> m]. real :: min_depth ! The minimum ocean depth in depth units [Z ~> m]. real :: slat, wlon, lenlat, lenlon, nlat diff --git a/src/user/BFB_surface_forcing.F90 b/src/user/BFB_surface_forcing.F90 index 3963d4d90d..6b17d64697 100644 --- a/src/user/BFB_surface_forcing.F90 +++ b/src/user/BFB_surface_forcing.F90 @@ -29,8 +29,6 @@ module BFB_surface_forcing real :: Rho0 !< The density used in the Boussinesq approximation [R ~> kg m-3]. real :: G_Earth !< The gravitational acceleration [L2 Z-1 T-2 ~> m s-2] real :: Flux_const !< The restoring rate at the surface [Z T-1 ~> m s-1]. - real :: gust_const !< A constant unresolved background gustiness - !! that contributes to ustar [Pa]. real :: SST_s !< SST at the southern edge of the linear forcing ramp [degC] real :: SST_n !< SST at the northern edge of the linear forcing ramp [degC] real :: lfrslat !< Southern latitude where the linear forcing ramp begins [degLat] @@ -220,8 +218,6 @@ subroutine BFB_surface_forcing_init(Time, G, US, param_file, diag, CS) call get_param(param_file, mdl, "DRHO_DT", CS%drho_dt, & "The rate of change of density with temperature.", & units="kg m-3 K-1", default=-0.2, scale=US%kg_m3_to_R) - call get_param(param_file, mdl, "GUST_CONST", CS%gust_const, & - "The background gustiness in the winds.", units="Pa", default=0.0) call get_param(param_file, mdl, "RESTOREBUOY", CS%restorebuoy, & "If true, the buoyancy fluxes drive the model back "//& diff --git a/src/user/DOME2d_initialization.F90 b/src/user/DOME2d_initialization.F90 index 42279be8e3..70b4bbc27d 100644 --- a/src/user/DOME2d_initialization.F90 +++ b/src/user/DOME2d_initialization.F90 @@ -362,7 +362,7 @@ subroutine DOME2d_initialize_sponges(G, GV, US, tv, depth_tot, param_file, use_A real :: S(SZI_(G),SZJ_(G),SZK_(GV)) ! A temporary array for salt [ppt] real :: h(SZI_(G),SZJ_(G),SZK_(GV)) ! A temporary array for thickness [H ~> m or kg m-2]. real :: eta(SZI_(G),SZJ_(G),SZK_(GV)+1) ! A temporary array for interface heights [Z ~> m] - real :: Idamp(SZI_(G),SZJ_(G)) ! The inverse damping rate [T-1 ~> s-1]. + real :: Idamp(SZI_(G),SZJ_(G)) ! The sponge damping rate [T-1 ~> s-1] real :: S_ref ! Reference salinity within the surface layer [ppt] real :: T_ref ! Reference temerature within the surface layer [degC] real :: S_range ! Range of salinities in the vertical [ppt] @@ -418,7 +418,7 @@ subroutine DOME2d_initialize_sponges(G, GV, US, tv, depth_tot, param_file, use_A call get_param(param_file, mdl, "T_RANGE", T_range, default=0.0) - ! Set the inverse damping rate as a function of position + ! Set the sponge damping rate as a function of position Idamp(:,:) = 0.0 do j=js,je ; do i=is,ie if (G%mask2dT(i,j) > 0.) then ! Only set damping rate for wet points diff --git a/src/user/DOME_initialization.F90 b/src/user/DOME_initialization.F90 index ee4491799a..b5c14517c2 100644 --- a/src/user/DOME_initialization.F90 +++ b/src/user/DOME_initialization.F90 @@ -160,7 +160,7 @@ subroutine DOME_initialize_sponges(G, GV, US, tv, depth_tot, PF, CSp) real :: eta(SZI_(G),SZJ_(G),SZK_(GV)+1) ! A temporary array for interface heights [Z ~> m]. real :: temp(SZI_(G),SZJ_(G),SZK_(GV)) ! A temporary array for other variables. ! - real :: Idamp(SZI_(G),SZJ_(G)) ! The inverse damping rate [T-1 ~> s-1]. + real :: Idamp(SZI_(G),SZJ_(G)) ! The sponge damping rate [T-1 ~> s-1]. real :: H0(SZK_(GV)) ! Interface heights [Z ~> m]. real :: min_depth ! The minimum depth at which to apply damping [Z ~> m] @@ -174,7 +174,7 @@ subroutine DOME_initialize_sponges(G, GV, US, tv, depth_tot, PF, CSp) eta(:,:,:) = 0.0 ; temp(:,:,:) = 0.0 ; Idamp(:,:) = 0.0 -! Here the inverse damping time [s-1], is set. Set Idamp to 0 ! +! 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 set up the sponges only where Idamp is positive! ! and mask2dT is 1. ! diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index 5fe228e278..617ac0da3d 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -440,7 +440,7 @@ subroutine ISOMIP_initialize_sponges(G, GV, US, tv, depth_tot, PF, use_ALE, CSp, real :: S(SZI_(G),SZJ_(G),SZK_(GV)) ! A temporary array for salt [ppt] ! real :: RHO(SZI_(G),SZJ_(G),SZK_(GV)) ! A temporary array for RHO [R ~> kg m-3] real :: h(SZI_(G),SZJ_(G),SZK_(GV)) ! A temporary array for thickness [H ~> m or kg m-2] - real :: Idamp(SZI_(G),SZJ_(G)) ! The inverse damping rate [T-1 ~> s-1]. + real :: Idamp(SZI_(G),SZJ_(G)) ! The sponge damping rate [T-1 ~> s-1] real :: TNUDG ! Nudging time scale [T ~> s] real :: S_sur, T_sur ! Surface salinity and temerature in sponge real :: S_bot, T_bot ! Bottom salinity and temerature in sponge @@ -660,7 +660,7 @@ subroutine ISOMIP_initialize_sponges(G, GV, US, tv, depth_tot, PF, use_ALE, CSp, ! call MOM_mesg(mesg,5) !enddo - ! Set the inverse damping rates so that the model will know where to + ! Set the sponge damping rates so that the model will know where to ! apply the sponges, along with the interface heights. call initialize_sponge(Idamp, eta, G, PF, CSp, GV) ! Apply sponge in tracer fields diff --git a/src/user/Phillips_initialization.F90 b/src/user/Phillips_initialization.F90 index a2dd76519d..ed7bc07ba3 100644 --- a/src/user/Phillips_initialization.F90 +++ b/src/user/Phillips_initialization.F90 @@ -237,7 +237,7 @@ subroutine Phillips_initialize_sponges(G, GV, US, tv, param_file, CSp, h) real :: eta0(SZK_(GV)+1) ! The 1-d nominal positions of the interfaces. real :: eta(SZI_(G),SZJ_(G),SZK_(GV)+1) ! A temporary array for interface heights [Z ~> m]. real :: temp(SZI_(G),SZJ_(G),SZK_(GV)) ! A temporary array for other variables. - real :: Idamp(SZI_(G),SZJ_(G)) ! The inverse damping rate [T-1 ~> s-1]. + real :: Idamp(SZI_(G),SZJ_(G)) ! The sponge damping rate [T-1 ~> s-1] real :: eta_im(SZJ_(G),SZK_(GV)+1) ! A temporary array for zonal-mean eta [Z ~> m]. real :: Idamp_im(SZJ_(G)) ! The inverse zonal-mean damping rate [T-1 ~> s-1]. real :: damp_rate ! The inverse zonal-mean damping rate [T-1 ~> s-1]. diff --git a/src/user/RGC_initialization.F90 b/src/user/RGC_initialization.F90 index d051bccc6c..b8eae3c704 100644 --- a/src/user/RGC_initialization.F90 +++ b/src/user/RGC_initialization.F90 @@ -76,7 +76,7 @@ subroutine RGC_initialize_sponges(G, GV, US, tv, u, v, depth_tot, PF, use_ALE, C real :: RHO(SZI_(G),SZJ_(G),SZK_(GV)) ! A temporary array for RHO real :: tmp(SZI_(G),SZJ_(G)) ! A temporary array for tracers. real :: h(SZI_(G),SZJ_(G),SZK_(GV)) ! A temporary array for thickness at h points - real :: Idamp(SZI_(G),SZJ_(G)) ! The inverse damping rate at h points [T-1 ~> s-1]. + real :: Idamp(SZI_(G),SZJ_(G)) ! The sponge damping rate at h points [T-1 ~> s-1] real :: TNUDG ! Nudging time scale [T ~> s] real :: pres(SZI_(G)) ! An array of the reference pressure [R L2 T-2 ~> Pa] real :: eta(SZI_(G),SZJ_(G),SZK_(GV)+1) ! A temporary array for eta, positive upward [m]. @@ -196,7 +196,7 @@ subroutine RGC_initialize_sponges(G, GV, US, tv, u, v, depth_tot, PF, use_ALE, C !read eta call MOM_read_data(filename, eta_var, eta(:,:,:), G%Domain) - ! Set the inverse damping rates so that the model will know where to + ! Set the sponge damping rates so that the model will know where to ! apply the sponges, along with the interface heights. call initialize_sponge(Idamp, eta, G, PF, CSp, GV) diff --git a/src/user/user_change_diffusivity.F90 b/src/user/user_change_diffusivity.F90 index 92570e3caa..0308a3b008 100644 --- a/src/user/user_change_diffusivity.F90 +++ b/src/user/user_change_diffusivity.F90 @@ -25,6 +25,7 @@ module user_change_diffusivity !> Control structure for user_change_diffusivity type, public :: user_change_diff_CS ; private + logical :: initialized = .false. !< True if this control structure has been initialized. real :: Kd_add !< The scale of a diffusivity that is added everywhere !! without any filtering or scaling [Z2 T-1 ~> m2 s-1]. real :: lat_range(4) !< 4 values that define the latitude range over which @@ -86,6 +87,9 @@ subroutine user_change_diff(h, tv, G, GV, US, CS, Kd_lay, Kd_int, T_f, S_f, Kd_i if (.not.associated(CS)) call MOM_error(FATAL,"user_set_diffusivity: "//& "Module must be initialized before it is used.") + if (.not.CS%initialized) call MOM_error(FATAL,"user_set_diffusivity: "//& + "Module must be initialized before it is used.") + use_EOS = associated(tv%eqn_of_state) if (.not.use_EOS) return store_Kd_add = .false. @@ -214,6 +218,8 @@ subroutine user_change_diff_init(Time, G, GV, US, param_file, diag, CS) endif allocate(CS) + CS%initialized = .true. + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec CS%diag => diag diff --git a/src/user/user_initialization.F90 b/src/user/user_initialization.F90 index d719e5867c..18b1fa5225 100644 --- a/src/user/user_initialization.F90 +++ b/src/user/user_initialization.F90 @@ -239,7 +239,7 @@ end subroutine write_user_log !! - GV%Rlay - Layer potential density (coordinate variable) [R ~> kg m-3]. !! If ENABLE_THERMODYNAMICS is defined: !! - T - Temperature [degC]. -!! - S - Salinity [psu]. +!! - S - Salinity [ppt]. !! If BULKMIXEDLAYER is defined: !! - Rml - Mixed layer and buffer layer potential densities [R ~> kg m-3]. !! If SPONGE is defined: