Merge branch 'eliminate_Zd_to_m' of https://github.com/Hallberg-NOAA/…

…MOM6 into Hallberg-NOAA-eliminate_Zd_to_m

config_src/coupled_driver/MOM_surface_forcing.F90

@@ -794,7 +794,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
   type(time_type),         intent(in)    :: Time !< The time of the fluxes, used for interpolating the
                                                  !! salinity to the right time, when it is being restored.
   type(ocean_grid_type),   intent(inout) :: G    !< The ocean's grid structure
-  type(unit_scale_type),    intent(in)    :: US   !< A dimensional unit scaling type
+  type(unit_scale_type),   intent(in)    :: US   !< A dimensional unit scaling type
   type(surface_forcing_CS),pointer       :: CS   !< A pointer to the control structure returned by a
                                                  !! previous call to surface_forcing_init.
   real, dimension(SZIB_(G),SZJ_(G)), &
@@ -805,7 +805,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
                  optional, intent(inout) :: ustar !< The surface friction velocity, in Z s-1.
   real, dimension(SZI_(G),SZJ_(G)), &
                  optional, intent(out)   :: gustless_ustar !< The surface friction velocity without
-                                                 !! any contributions from gustiness, in m s-1.
+                                                 !! any contributions from gustiness, in Z s-1.
   integer,       optional, intent(in)    :: tau_halo !< The halo size of wind stresses to set, 0 by default.
 
   ! Local variables
@@ -817,7 +817,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
   real, dimension(SZIB_(G),SZJB_(G)) :: tauy_in_B ! Meridional wind stresses (in Pa) at q points
 
   real :: gustiness     ! unresolved gustiness that contributes to ustar (Pa)
-  real :: Irho0         ! inverse of the mean density in (m^3/kg)
+  real :: Irho0         ! Inverse of the mean density rescaled to (Z2 m / kg)
   real :: taux2, tauy2  ! squared wind stresses (Pa^2)
   real :: tau_mag       ! magnitude of the wind stress (Pa)
 
@@ -831,7 +831,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
   Isqh = G%IscB-halo ; Ieqh  = G%IecB+halo ; Jsqh = G%JscB-halo ; Jeqh = G%JecB+halo
   i0 = is - index_bounds(1) ; j0 = js - index_bounds(3)
 
-  Irho0 = 1.0/CS%Rho0
+  Irho0 = US%m_to_Z**2 / CS%Rho0
 
   do_ustar = present(ustar) ; do_gustless = present(gustless_ustar)
 
@@ -943,10 +943,10 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
                 (G%mask2dBu(I,J-1) + G%mask2dBu(I-1,J))) > 0)) ) &
             gustiness = CS%gust(i,j)
         endif
-        ustar(i,j) = US%m_to_Z * sqrt(gustiness*Irho0 + Irho0*IOB%stress_mag(i-i0,j-j0))
+        ustar(i,j) = sqrt(gustiness*Irho0 + Irho0*IOB%stress_mag(i-i0,j-j0))
       enddo ; enddo ; endif
       if (do_gustless) then ; do j=js,je ; do i=is,ie
-        gustless_ustar(i,j) = sqrt(IOB%stress_mag(i-i0,j-j0) / CS%Rho0)
+        gustless_ustar(i,j) = US%m_to_Z * sqrt(IOB%stress_mag(i-i0,j-j0) / CS%Rho0)
 !### Change to:
 !        gustless_ustar(i,j) = sqrt(Irho0 * IOB%stress_mag(i-i0,j-j0))
       enddo ; enddo ; endif
@@ -962,8 +962,8 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
             ((G%mask2dBu(I,J) + G%mask2dBu(I-1,J-1)) + (G%mask2dBu(I,J-1) + G%mask2dBu(I-1,J))) )
           if (CS%read_gust_2d) gustiness = CS%gust(i,j)
         endif
-        if (do_ustar) ustar(i,j) = US%m_to_Z * sqrt(gustiness*Irho0 + Irho0 * tau_mag)
-        if (do_gustless) gustless_ustar(i,j) = sqrt(tau_mag / CS%Rho0)
+        if (do_ustar) ustar(i,j) = sqrt(gustiness*Irho0 + Irho0 * tau_mag)
+        if (do_gustless) gustless_ustar(i,j) = US%m_to_Z * sqrt(tau_mag / CS%Rho0)
 !### Change to:
 !        if (do_gustless) gustless_ustar(i,j) = sqrt(Irho0 * tau_mag)
       enddo ; enddo
@@ -972,8 +972,8 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
         tau_mag = G%mask2dT(i,j) * sqrt(taux_in_A(i,j)**2 + tauy_in_A(i,j)**2)
         gustiness = CS%gust_const
         if (CS%read_gust_2d .and. (G%mask2dT(i,j) > 0)) gustiness = CS%gust(i,j)
-        if (do_ustar) ustar(i,j) = US%m_to_Z * sqrt(gustiness*Irho0 + Irho0 * tau_mag)
-        if (do_gustless) gustless_ustar(i,j) = sqrt(tau_mag / CS%Rho0)
+        if (do_ustar) ustar(i,j) = sqrt(gustiness*Irho0 + Irho0 * tau_mag)
+        if (do_gustless) gustless_ustar(i,j) = US%m_to_Z * sqrt(tau_mag / CS%Rho0)
 !### Change to:
 !        if (do_gustless) gustless_ustar(i,j) = sqrt(Irho0 * tau_mag)
       enddo ; enddo
@@ -991,8 +991,8 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
         gustiness = CS%gust_const
         if (CS%read_gust_2d) gustiness = CS%gust(i,j)
 
-        if (do_ustar) ustar(i,j) = US%m_to_Z * sqrt(gustiness*Irho0 + Irho0 * tau_mag)
-        if (do_gustless) gustless_ustar(i,j) = sqrt(tau_mag / CS%Rho0)
+        if (do_ustar) ustar(i,j) = sqrt(gustiness*Irho0 + Irho0 * tau_mag)
+        if (do_gustless) gustless_ustar(i,j) = US%m_to_Z * sqrt(tau_mag / CS%Rho0)
 !### Change to:
 !        if (do_gustless) gustless_ustar(i,j) = sqrt(Irho0 * tau_mag)
       enddo ; enddo

config_src/ice_solo_driver/MOM_surface_forcing.F90

@@ -167,15 +167,17 @@ module MOM_surface_forcing
 
 contains
 
-subroutine set_forcing(sfc_state, forcing, fluxes, day_start, day_interval, G, CS)
+subroutine set_forcing(sfc_state, forcing, fluxes, day_start, day_interval, G, US, CS)
   type(surface),         intent(inout) :: sfc_state !< A structure containing fields that
                                                     !! describe the surface state of the ocean.
   type(mech_forcing),    intent(inout) :: forces !< A structure with the driving mechanical forces
-  type(forcing),         intent(inout) :: fluxes
-  type(time_type),       intent(in)    :: day_start
-  type(time_type),       intent(in)    :: day_interval
+  type(forcing),         intent(inout) :: fluxes !< A structure containing thermodynamic forcing fields
+  type(time_type),       intent(in)    :: day_start !< The start time of the fluxes
+  type(time_type),       intent(in)    :: day_interval !< Length of time over which these fluxes applied
   type(ocean_grid_type), intent(inout) :: G    !< The ocean's grid structure
-  type(surface_forcing_CS), pointer    :: CS
+  type(unit_scale_type), intent(in)    :: US   !< A dimensional unit scaling type
+  type(surface_forcing_CS), pointer    :: CS   !< pointer to control struct returned by
+                                               !! a previous surface_forcing_init call
 
 ! This subroutine calls other subroutines in this file to get surface forcing fields.
 ! It also allocates and initializes the fields in the flux type.
@@ -282,7 +284,7 @@ subroutine set_forcing(sfc_state, forcing, fluxes, day_start, day_interval, G, C
   ! Fields that exist in both the forcing and mech_forcing types must be copied.
   if (CS%variable_winds .or. CS%first_call_set_forcing) then
     call copy_common_forcing_fields(forces, fluxes, G)
-    call set_derived_forcing_fields(forces, fluxes, G, CS%Rho0)
+    call set_derived_forcing_fields(forces, fluxes, G, US, CS%Rho0)
   endif
 
   if ((CS%variable_buoyforce .or. CS%first_call_set_forcing) .and. &

config_src/mct_driver/MOM_ocean_model.F90

@@ -542,7 +542,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
 #endif
   endif
 
-  call set_derived_forcing_fields(OS%forces, OS%fluxes, OS%grid, OS%GV%Rho0)
+  call set_derived_forcing_fields(OS%forces, OS%fluxes, OS%grid, OS%US, OS%GV%Rho0)
   call set_net_mass_forcing(OS%fluxes, OS%forces, OS%grid)
 
   if (OS%nstep==0) then

config_src/solo_driver/MOM_surface_forcing.F90

@@ -340,7 +340,7 @@ subroutine set_forcing(sfc_state, forces, fluxes, day_start, day_interval, G, US
   ! Fields that exist in both the forcing and mech_forcing types must be copied.
   if (CS%variable_winds .or. CS%first_call_set_forcing) then
     call copy_common_forcing_fields(forces, fluxes, G)
-    call set_derived_forcing_fields(forces, fluxes, G, CS%Rho0)
+    call set_derived_forcing_fields(forces, fluxes, G, US, CS%Rho0)
   endif
 
   if ((CS%variable_buoyforce .or. CS%first_call_set_forcing) .and. &

src/core/MOM.F90

@@ -72,7 +72,6 @@ module MOM
 use MOM_dynamics_unsplit_RK2,  only : initialize_dyn_unsplit_RK2, end_dyn_unsplit_RK2
 use MOM_dynamics_unsplit_RK2,  only : MOM_dyn_unsplit_RK2_CS
 use MOM_dyn_horgrid,           only : dyn_horgrid_type, create_dyn_horgrid, destroy_dyn_horgrid
-use MOM_dyn_horgrid,           only : rescale_dyn_horgrid_bathymetry
 use MOM_EOS,                   only : EOS_init, calculate_density, calculate_TFreeze
 use MOM_fixed_initialization,  only : MOM_initialize_fixed
 use MOM_grid,                  only : ocean_grid_type, MOM_grid_init, MOM_grid_end
@@ -1978,11 +1977,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
   call MOM_timing_init(CS)
 
   ! Allocate initialize time-invariant MOM variables.
-  call MOM_initialize_fixed(dG, CS%OBC, param_file, write_geom_files, dirs%output_directory)
+  call MOM_initialize_fixed(dG, US, CS%OBC, param_file, write_geom_files, dirs%output_directory)
   call callTree_waypoint("returned from MOM_initialize_fixed() (initialize_MOM)")
 
-  if (dG%Zd_to_m /= US%Z_to_m) call rescale_dyn_horgrid_bathymetry(dG, US%Z_to_m)
-
   if (associated(CS%OBC)) call call_OBC_register(param_file, CS%update_OBC_CSp, CS%OBC)
 
   call tracer_registry_init(param_file, CS%tracer_Reg)
@@ -2256,7 +2253,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
 
   diag => CS%diag
   ! Initialize the diag mediator.
-  call diag_mediator_init(G, GV, GV%ke, param_file, diag, doc_file_dir=dirs%output_directory)
+  call diag_mediator_init(G, GV, US, GV%ke, param_file, diag, doc_file_dir=dirs%output_directory)
   if (present(diag_ptr)) diag_ptr => CS%diag
 
   ! Initialize the diagnostics masks for native arrays.
@@ -2455,7 +2452,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
 
   CS%nstep_tot = 0
   if (present(count_calls)) CS%count_calls = count_calls
-  call MOM_sum_output_init(G, param_file, dirs%output_directory, &
+  call MOM_sum_output_init(G, US, param_file, dirs%output_directory, &
                            CS%ntrunc, Time_init, CS%sum_output_CSp)
 
   ! Flag whether to save initial conditions in finish_MOM_initialization() or not.
@@ -2982,7 +2979,7 @@ subroutine extract_surface_state(CS, sfc_state)
     numberOfErrors=0 ! count number of errors
     do j=js,je; do i=is,ie
       if (G%mask2dT(i,j)>0.) then
-        bathy_m = G%Zd_to_m*G%bathyT(i,j)
+        bathy_m = CS%US%Z_to_m * G%bathyT(i,j)
         localError = sfc_state%sea_lev(i,j)<=-bathy_m &
                 .or. sfc_state%sea_lev(i,j)>= CS%bad_val_ssh_max  &
                 .or. sfc_state%sea_lev(i,j)<=-CS%bad_val_ssh_max  &

src/core/MOM_barotropic.F90

@@ -2272,10 +2272,10 @@ subroutine set_dtbt(G, GV, US, CS, eta, pbce, BT_cont, gtot_est, SSH_add)
                                                       !! the effective open face areas as a
                                                       !! function of barotropic flow.
   real,               optional, intent(in)    :: gtot_est !< An estimate of the total gravitational
-                                                      !! acceleration, in m s-2.
+                                                      !! acceleration, in m2 Z-1 s-2.
   real,               optional, intent(in)    :: SSH_add  !< An additional contribution to SSH to
                                                       !! provide a margin of error when
-                                                      !! calculating the external wave speed, in m.
+                                                      !! calculating the external wave speed, in Z.
 
   ! Local variables
   real, dimension(SZI_(G),SZJ_(G)) :: &
@@ -2299,7 +2299,7 @@ subroutine set_dtbt(G, GV, US, CS, eta, pbce, BT_cont, gtot_est, SSH_add)
                   ! order 1.  For stability, this may be made larger
                   ! than physical problem would suggest.
   real :: add_SSH ! An additional contribution to SSH to provide a margin of error
-                  ! when calculating the external wave speed, in m.
+                  ! when calculating the external wave speed, in Z.
   real :: min_max_dt2, Idt_max2, dtbt_max
   logical :: use_BT_cont
   type(memory_size_type) :: MS
@@ -2326,7 +2326,7 @@ subroutine set_dtbt(G, GV, US, CS, eta, pbce, BT_cont, gtot_est, SSH_add)
   elseif (CS%Nonlinear_continuity .and. present(eta)) then
     call find_face_areas(Datu, Datv, G, GV, CS, MS, eta=eta, halo=0)
   else
-    call find_face_areas(Datu, Datv, G, GV, CS, MS, halo=0, add_max=add_SSH*US%m_to_Z)
+    call find_face_areas(Datu, Datv, G, GV, CS, MS, halo=0, add_max=add_SSH)
   endif
 
   det_de = 0.0
@@ -2345,8 +2345,8 @@ subroutine set_dtbt(G, GV, US, CS, eta, pbce, BT_cont, gtot_est, SSH_add)
     enddo ; enddo ; enddo
   else
     do j=js,je ; do i=is,ie
-      gtot_E(i,j) = gtot_est * GV%H_to_m ; gtot_W(i,j) = gtot_est * GV%H_to_m
-      gtot_N(i,j) = gtot_est * GV%H_to_m ; gtot_S(i,j) = gtot_est * GV%H_to_m
+      gtot_E(i,j) = gtot_est * GV%H_to_Z ; gtot_W(i,j) = gtot_est * GV%H_to_Z
+      gtot_N(i,j) = gtot_est * GV%H_to_Z ; gtot_S(i,j) = gtot_est * GV%H_to_Z
     enddo ; enddo
   endif
 
@@ -3717,9 +3717,9 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS,
   character(len=40)  :: mdl = "MOM_barotropic"  ! This module's name.
   real :: Datu(SZIBS_(G),SZJ_(G))   ! Zonal open face area in H m.
   real :: Datv(SZI_(G),SZJBS_(G))   ! Meridional open face area in H m.
-  real :: gtot_estimate ! Summing GV%g_prime gives an upper-bound estimate for pbce.
+  real :: gtot_estimate ! Summed GV%g_prime in m2 Z-1 s-2, to give an upper-bound estimate for pbce.
   real :: SSH_extra     ! An estimate of how much higher SSH might get, for use
-                        ! in calculating the safe external wave speed.
+                        ! in calculating the safe external wave speed, in Z.
   real :: dtbt_input, dtbt_tmp
   real :: wave_drag_scale ! A scaling factor for the barotropic linear wave drag
                           ! piston velocities.
@@ -3956,7 +3956,7 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS,
                  "An estimate of how much higher SSH might get, for use \n"//&
                  "in calculating the safe external wave speed. The \n"//&
                  "default is the minimum of 10 m or 5% of MAXIMUM_DEPTH.", &
-                 units="m", default=min(10.0,0.05*G%max_depth*US%Z_to_m))
+                 units="m", default=min(10.0,0.05*G%max_depth*US%Z_to_m), scale=US%m_to_Z)
 
   call get_param(param_file, mdl, "DEBUG", CS%debug, &
                  "If true, write out verbose debugging data.", &
@@ -4141,8 +4141,8 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS,
 
   ! Estimate the maximum stable barotropic time step.
   gtot_estimate = 0.0
-  do k=1,G%ke ; gtot_estimate = gtot_estimate + GV%g_prime(K)*US%m_to_Z ; enddo
-  call set_dtbt(G, GV, US, CS, gtot_est = gtot_estimate, SSH_add = SSH_extra)
+  do k=1,G%ke ; gtot_estimate = gtot_estimate + GV%g_prime(K) ; enddo
+  call set_dtbt(G, GV, US, CS, gtot_est=gtot_estimate, SSH_add=SSH_extra)
 
   if (dtbt_input > 0.0) then
     CS%dtbt = dtbt_input

src/core/MOM_forcing_type.F90

@@ -47,7 +47,7 @@ module MOM_forcing_type
   real, pointer, dimension(:,:) :: &
     ustar         => NULL(), & !< surface friction velocity scale (Z/s)
     ustar_gustless => NULL()   !< surface friction velocity scale without any
-                               !! any augmentation for gustiness (m/s)
+                               !! any augmentation for gustiness (Z/s)
 
   ! surface buoyancy force, used when temperature is not a state variable
   real, pointer, dimension(:,:) :: &
@@ -1961,19 +1961,20 @@ end subroutine copy_common_forcing_fields
 
 !> This subroutine calculates certain derived forcing fields based on information
 !! from a mech_forcing type and stores them in a (thermodynamic) forcing type.
-subroutine set_derived_forcing_fields(forces, fluxes, G, Rho0)
+subroutine set_derived_forcing_fields(forces, fluxes, G, US, Rho0)
   type(mech_forcing),      intent(in)    :: forces   !< A structure with the driving mechanical forces
   type(forcing),           intent(inout) :: fluxes   !< A structure containing thermodynamic forcing fields
   type(ocean_grid_type),   intent(in)    :: G        !< grid type
+  type(unit_scale_type),   intent(in)    :: US       !< A dimensional unit scaling type
   real,                    intent(in)    :: Rho0     !< A reference density of seawater, in kg m-3,
                                                      !! as used to calculate ustar.
 
   real :: taux2, tauy2 ! Squared wind stress components, in Pa^2.
-  real :: Irho0        ! Inverse of the mean density in (m^3/kg)
+  real :: Irho0        ! Inverse of the mean density rescaled to (Z2 m / kg)
   integer :: i, j, is, ie, js, je
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
 
-  Irho0 = 1.0/Rho0
+  Irho0 = US%m_to_Z**2 / Rho0
 
   if (associated(forces%taux) .and. associated(forces%tauy) .and. &
       associated(fluxes%ustar_gustless)) then
@@ -1989,7 +1990,7 @@ subroutine set_derived_forcing_fields(forces, fluxes, G, Rho0)
                  G%mask2dCv(i,J) * forces%tauy(i,J)**2) / &
                 (G%mask2dCv(i,J-1) + G%mask2dCv(i,J))
 
-      fluxes%ustar_gustless(i,j) = sqrt(sqrt(taux2 + tauy2) / Rho0)
+      fluxes%ustar_gustless(i,j) = US%m_to_Z * sqrt(sqrt(taux2 + tauy2) / Rho0)
 !### Change to:
 !      fluxes%ustar_gustless(i,j) = sqrt(sqrt(taux2 + tauy2) * Irho0)
     enddo ; enddo

src/core/MOM_grid.F90

@@ -132,7 +132,6 @@ module MOM_grid
 
   real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: &
     bathyT        !< Ocean bottom depth at tracer points, in depth units.
-  real :: Zd_to_m  = 1.0 !< The conversion factor between the units of bathyT and m.
 
   logical :: bathymetry_at_vel  !< If true, there are separate values for the
                   !! basin depths at velocity points.  Otherwise the effects of
@@ -165,7 +164,7 @@ module MOM_grid
   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 in meters.
-  real :: max_depth     !< The maximum depth of the ocean in depth units (scaled by Zd_to_m).
+  real :: max_depth     !< The maximum depth of the ocean in depth units (Z).
 end type ocean_grid_type
 
 contains
@@ -359,13 +358,13 @@ subroutine rescale_grid_bathymetry(G, m_in_new_units)
   isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
   IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
 
-  if (m_in_new_units == G%Zd_to_m) return
+  if (m_in_new_units == 1.0) return
   if (m_in_new_units < 0.0) &
     call MOM_error(FATAL, "rescale_grid_bathymetry: Negative depth units are not permitted.")
   if (m_in_new_units == 0.0) &
     call MOM_error(FATAL, "rescale_grid_bathymetry: Zero depth units are not permitted.")
 
-  rescale = G%Zd_to_m / m_in_new_units
+  rescale = 1.0 / m_in_new_units
   do j=jsd,jed ; do i=isd,ied
     G%bathyT(i,j) = rescale*G%bathyT(i,j)
   enddo ; enddo
@@ -376,7 +375,6 @@ subroutine rescale_grid_bathymetry(G, m_in_new_units)
     G%Dblock_v(i,J) = rescale*G%Dblock_v(i,J) ; G%Dopen_v(i,J) = rescale*G%Dopen_v(i,J)
   enddo ; enddo ; endif
   G%max_depth = rescale*G%max_depth
-  G%Zd_to_m = m_in_new_units
 
 end subroutine rescale_grid_bathymetry