+Rescale tides and ramp-up times

src/core/MOM_PressureForce_FV.F90

@@ -153,7 +153,7 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_
                         ! [H T2 R-1 L-2 ~> m Pa-1 or kg m-2 Pa-1].
   real :: H_to_RL2_T2   ! A factor to convert from thickness units (H) to pressure
                         ! units [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1].
-!  real :: oneatm = 101325.0  ! 1 atm in [Pa] = [kg m-1 s-2]
+!  real :: oneatm       ! 1 standard atmosphere of pressure in [R L2 T-2 ~> Pa]
   real, parameter :: C1_6 = 1.0/6.0
   integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb
   integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
@@ -187,6 +187,7 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_
       p(i,j,1) = p_atm(i,j)
     enddo ; enddo
   else
+    ! oneatm = 101325.0 * US%kg_m3_to_R * US%m_s_to_L_T**2 ! 1 atm scaled to [R L2 T-2 ~> Pa]
     !$OMP parallel do default(shared)
     do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
       p(i,j,1) = 0.0 ! or oneatm
@@ -305,7 +306,7 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_
     do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
       SSH(i,j) = (za(i,j) - alpha_ref*p(i,j,1)) * I_gEarth - G%Z_ref
     enddo ; enddo
-    call calc_tidal_forcing(CS%Time, SSH, e_tidal, G, CS%tides_CSp, m_to_Z=US%m_to_Z)
+    call calc_tidal_forcing(CS%Time, SSH, e_tidal, G, US, CS%tides_CSp)
     !$OMP parallel do default(shared)
     do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
       za(i,j) = za(i,j) - GV%g_Earth * e_tidal(i,j)
@@ -573,7 +574,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
         SSH(i,j) = SSH(i,j) + h(i,j,k)*GV%H_to_Z
       enddo ; enddo
     enddo
-    call calc_tidal_forcing(CS%Time, SSH, e_tidal, G, CS%tides_CSp, m_to_Z=US%m_to_Z)
+    call calc_tidal_forcing(CS%Time, SSH, e_tidal, G, US, CS%tides_CSp)
   endif
 
 !    Here layer interface heights, e, are calculated.

src/core/MOM_PressureForce_Montgomery.F90

@@ -203,7 +203,7 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb
       enddo ; enddo ; enddo
     endif
 
-    call calc_tidal_forcing(CS%Time, SSH, e_tidal, G, CS%tides_CSp, m_to_Z=US%m_to_Z)
+    call calc_tidal_forcing(CS%Time, SSH, e_tidal, G, US, CS%tides_CSp)
     !$OMP parallel do default(shared)
     do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
       geopot_bot(i,j) = -GV%g_Earth*(e_tidal(i,j) + G%bathyT(i,j))
@@ -451,7 +451,7 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce,
         SSH(i,j) = SSH(i,j) + h(i,j,k)*GV%H_to_Z
       enddo ; enddo
     enddo
-    call calc_tidal_forcing(CS%Time, SSH, e_tidal, G, CS%tides_CSp, m_to_Z=US%m_to_Z)
+    call calc_tidal_forcing(CS%Time, SSH, e_tidal, G, US, CS%tides_CSp)
   endif
 
 !    Here layer interface heights, e, are calculated.

src/core/MOM_continuity_PPM.F90

@@ -924,8 +924,8 @@ subroutine set_zonal_BT_cont(u, h_in, h_L, h_R, BT_cont, uh_tot_0, duhdu_tot_0,
     FAmt_0, &     ! test velocities [H L ~> m2 or kg m-1].
     uhtot_L, &    ! The summed transport with the westerly (uhtot_L) and
     uhtot_R       ! and easterly (uhtot_R) test velocities [H L2 T-1 ~> m3 s-1 or kg s-1].
-  real :: FA_0    ! The effective face area with 0 barotropic transport [L H ~> m2 or kg m].
-  real :: FA_avg  ! The average effective face area [L H ~> m2 or kg m], nominally given by
+  real :: FA_0    ! The effective face area with 0 barotropic transport [L H ~> m2 or kg m-1].
+  real :: FA_avg  ! The average effective face area [L H ~> m2 or kg m-1], nominally given by
                   ! the realized transport divided by the barotropic velocity.
   real :: visc_rem_lim ! The larger of visc_rem and min_visc_rem [nondim]. This
                        ! limiting is necessary to keep the inverse of visc_rem

src/core/MOM_dynamics_split_RK2.F90

@@ -337,8 +337,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
 
   real, pointer, dimension(:,:,:) :: &
     ! These pointers are used to alter which fields are passed to btstep with various options:
-    u_ptr => NULL(), &   ! A pointer to a zonal velocity [L T-1]
-    v_ptr => NULL(), &   ! A pointer to a meridional velocity [L T-1]
+    u_ptr => NULL(), &   ! A pointer to a zonal velocity [L T-1 ~> m s-1]
+    v_ptr => NULL(), &   ! A pointer to a meridional velocity [L T-1 ~> m s-1]
     uh_ptr => NULL(), &  ! A pointer to a zonal volume or mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
     vh_ptr => NULL(), &  ! A pointer to a meridional volume or mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
     ! These pointers are just used as shorthand for CS%u_av, CS%v_av, and CS%h_av.
@@ -374,7 +374,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
   enddo
 
   ! Update CFL truncation value as function of time
-  call updateCFLtruncationValue(Time_local, CS%vertvisc_CSp)
+  call updateCFLtruncationValue(Time_local, CS%vertvisc_CSp, US)
 
   if (CS%debug) then
     call MOM_state_chksum("Start predictor ", u, v, h, uh, vh, G, GV, US, symmetric=sym)
@@ -395,7 +395,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
     if (CS%debug_OBC) call open_boundary_test_extern_h(G, GV, CS%OBC, h)
 
     ! Update OBC ramp value as function of time
-    call update_OBC_ramp(Time_local, CS%OBC)
+    call update_OBC_ramp(Time_local, CS%OBC, US)
 
     do k=1,nz ; do j=G%jsd,G%jed ; do I=G%IsdB,G%IedB
       u_old_rad_OBC(I,j,k) = u_av(I,j,k)
@@ -1207,20 +1207,20 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   call continuity_init(Time, G, GV, US, param_file, diag, CS%continuity_CSp)
   cont_stencil = continuity_stencil(CS%continuity_CSp)
   call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv)
-  if (use_tides) call tidal_forcing_init(Time, G, param_file, CS%tides_CSp)
+  if (use_tides) call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp)
   call PressureForce_init(Time, G, GV, US, param_file, diag, CS%PressureForce_CSp, &
                           CS%tides_CSp)
   call hor_visc_init(Time, G, GV, US, param_file, diag, CS%hor_visc, ADp=CS%ADp)
   call vertvisc_init(MIS, Time, G, GV, US, param_file, diag, CS%ADp, dirs, &
                      ntrunc, CS%vertvisc_CSp)
   CS%set_visc_CSp => set_visc
-  call updateCFLtruncationValue(Time, CS%vertvisc_CSp, &
+  call updateCFLtruncationValue(Time, CS%vertvisc_CSp, US, &
                                 activate=is_new_run(restart_CS) )
 
   if (associated(ALE_CSp)) CS%ALE_CSp => ALE_CSp
   if (associated(OBC)) then
     CS%OBC => OBC
-    if (OBC%ramp) call update_OBC_ramp(Time, CS%OBC, &
+    if (OBC%ramp) call update_OBC_ramp(Time, CS%OBC, US, &
                                 activate=is_new_run(restart_CS) )
   endif
   if (associated(update_OBC_CSp)) CS%update_OBC_CSp => update_OBC_CSp

src/core/MOM_dynamics_unsplit.F90

@@ -666,7 +666,7 @@ subroutine initialize_dyn_unsplit(u, v, h, Time, G, GV, US, param_file, diag, CS
   call continuity_init(Time, G, GV, US, param_file, diag, CS%continuity_CSp)
   cont_stencil = continuity_stencil(CS%continuity_CSp)
   call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv)
-  if (use_tides) call tidal_forcing_init(Time, G, param_file, CS%tides_CSp)
+  if (use_tides) call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp)
   call PressureForce_init(Time, G, GV, US, param_file, diag, CS%PressureForce_CSp, &
                           CS%tides_CSp)
   call hor_visc_init(Time, G, GV, US, param_file, diag, CS%hor_visc)

src/core/MOM_dynamics_unsplit_RK2.F90

@@ -628,7 +628,7 @@ subroutine initialize_dyn_unsplit_RK2(u, v, h, Time, G, GV, US, param_file, diag
   call continuity_init(Time, G, GV, US, param_file, diag, CS%continuity_CSp)
   cont_stencil = continuity_stencil(CS%continuity_CSp)
   call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv)
-  if (use_tides) call tidal_forcing_init(Time, G, param_file, CS%tides_CSp)
+  if (use_tides) call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp)
   call PressureForce_init(Time, G, GV, US, param_file, diag, CS%PressureForce_CSp, &
                           CS%tides_CSp)
   call hor_visc_init(Time, G, GV, US, param_file, diag, CS%hor_visc)

src/core/MOM_open_boundary.F90

@@ -264,7 +264,7 @@ module MOM_open_boundary
   logical :: add_tide_constituents = .false.          !< If true, add tidal constituents to the boundary elevation
                                                       !! and velocity. Will be set to true if n_tide_constituents > 0.
   character(len=2), allocatable, dimension(:) :: tide_names  !< Names of tidal constituents to add to the boundary data.
-  real, allocatable, dimension(:) :: tide_frequencies        !< Angular frequencies of chosen tidal constituents [s-1].
+  real, allocatable, dimension(:) :: tide_frequencies !< Angular frequencies of chosen tidal constituents [T-1 ~> s-1].
   real, allocatable, dimension(:) :: tide_eq_phases   !< Equilibrium phases of chosen tidal constituents [rad].
   real, allocatable, dimension(:) :: tide_fn          !< Amplitude modulation of boundary tides by nodal cycle [nondim].
   real, allocatable, dimension(:) :: tide_un          !< Phase modulation of boundary tides by nodal cycle [rad].
@@ -305,8 +305,8 @@ module MOM_open_boundary
                    !! the independence of the OBCs to this external data [L T-1 ~> m s-1].
   logical :: ramp = .false.                 !< If True, ramp from zero to the external values for SSH.
   logical :: ramping_is_activated = .false. !< True if the ramping has been initialized
-  real :: ramp_timescale                    !< If ramp is True, use this timescale for ramping [s].
-  real :: trunc_ramp_time                   !< If ramp is True, time after which ramp is done [s].
+  real :: ramp_timescale                    !< If ramp is True, use this timescale for ramping [T ~> s].
+  real :: trunc_ramp_time                   !< If ramp is True, time after which ramp is done [T ~> s].
   real :: ramp_value                        !< If ramp is True, where we are on the ramp from
                                             !! zero to one [nondim].
   type(time_type) :: ramp_start_time        !< Time when model was started.
@@ -627,7 +627,7 @@ subroutine open_boundary_config(G, US, param_file, OBC)
        "Symmetric memory must be used when using Flather OBCs.")
   ! Need to do this last, because it depends on time_interp_external_init having already been called
   if (OBC%add_tide_constituents) then
-    call initialize_obc_tides(OBC, param_file)
+    call initialize_obc_tides(OBC, US, param_file)
     ! Tide update is done within update_OBC_segment_data, so this should be true if tides are included.
     OBC%update_OBC = .true.
   endif
@@ -948,8 +948,9 @@ subroutine initialize_segment_data(G, OBC, PF)
 
 end subroutine initialize_segment_data
 
-subroutine initialize_obc_tides(OBC, param_file)
+subroutine initialize_obc_tides(OBC, US, param_file)
   type(ocean_OBC_type), intent(inout) :: OBC  !< Open boundary control structure
+  type(unit_scale_type),   intent(in) :: US   !< A dimensional unit scaling type
   type(param_file_type), intent(in) :: param_file !< Parameter file handle
   integer, dimension(3) :: tide_ref_date      !< Reference date (t = 0) for tidal forcing (year, month, day).
   integer, dimension(3) :: nodal_ref_date     !< Date to calculate nodal modulation for (year, month, day).
@@ -1022,7 +1023,8 @@ subroutine initialize_obc_tides(OBC, param_file)
         "Frequency of the "//trim(OBC%tide_names(c))//" tidal constituent. "//&
         "This is only used if TIDES and TIDE_"//trim(OBC%tide_names(c))// &
         " are true, or if OBC_TIDE_N_CONSTITUENTS > 0 and "//trim(OBC%tide_names(c))//&
-        " is in OBC_TIDE_CONSTITUENTS.", units="s-1", default=tidal_frequency(trim(OBC%tide_names(c))))
+        " is in OBC_TIDE_CONSTITUENTS.", &
+        units="s-1", default=tidal_frequency(trim(OBC%tide_names(c))), scale=US%T_to_s)
 
     ! Find equilibrium phase if needed
     if (OBC%add_eq_phase) then
@@ -3727,7 +3729,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
   real :: tidal_elev  ! Interpolated tidal elevation at the OBC points [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 [s]
+  real :: time_delta  ! Time since tidal reference date [T ~> 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
@@ -3738,7 +3740,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
 
   if (.not. associated(OBC)) return
 
-  if (OBC%add_tide_constituents) time_delta = time_type_to_real(Time - OBC%time_ref)
+  if (OBC%add_tide_constituents) time_delta = US%s_to_T * time_type_to_real(Time - OBC%time_ref)
 
   do n = 1, OBC%number_of_segments
     segment => OBC%segment(n)
@@ -4336,14 +4338,15 @@ end subroutine update_OBC_segment_data
 !> Update the OBC ramp value as a function of time.
 !! If called with the optional argument activate=.true., record the
 !! value of Time as the beginning of the ramp period.
-subroutine update_OBC_ramp(Time, OBC, activate)
+subroutine update_OBC_ramp(Time, OBC, US, activate)
   type(time_type), target, intent(in)    :: Time     !< Current model time
   type(ocean_OBC_type),    intent(inout) :: OBC      !< Open boundary structure
+  type(unit_scale_type),   intent(in)    :: US       !< A dimensional unit scaling type
   logical, optional,       intent(in)    :: activate !< Specify whether to record the value of
                                                      !! Time as the beginning of the ramp period
 
   ! Local variables
-  real :: deltaTime ! The time since start of ramping [s]
+  real :: deltaTime ! The time since start of ramping [T ~> s]
   real :: wghtA     ! A temporary variable used to set OBC%ramp_value [nondim]
   character(len=12) :: msg
 
@@ -4359,7 +4362,7 @@ subroutine update_OBC_ramp(Time, OBC, activate)
     endif
   endif
   if (.not.OBC%ramping_is_activated) return
-  deltaTime = max( 0., time_type_to_real( Time - OBC%ramp_start_time ) )
+  deltaTime = max( 0., US%s_to_T*time_type_to_real( Time - OBC%ramp_start_time ) )
   if (deltaTime >= OBC%trunc_ramp_time) then
     OBC%ramp_value = 1.0
     OBC%ramp = .false. ! This turns off ramping after this call

src/diagnostics/MOM_wave_structure.F90

@@ -171,20 +171,20 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo
   real, dimension(SZK_(GV))   :: dz           !< thicknesses of merged layers (same as Hc I hope) [Z ~> m]
   ! real, dimension(SZK_(GV)+1) :: dWdz_profile !< profile of dW/dz
   real                        :: w2avg   !< average of squared vertical velocity structure funtion [Z ~> m]
-  real                        :: int_dwdz2
-  real                        :: int_w2
-  real                        :: int_N2w2
-  real                        :: KE_term !< terms in vertically averaged energy equation
-  real                        :: PE_term !< terms in vertically averaged energy equation
+  real                        :: int_dwdz2 !< Vertical integral of the square of u_strct [Z ~> m]
+  real                        :: int_w2   !< Vertical integral of the square of w_strct [Z ~> m]
+  real                        :: int_N2w2 !< Vertical integral of N2 [Z T-2 ~> m s-2]
+  real                        :: KE_term !< terms in vertically averaged energy equation [R Z ~> kg m-2]
+  real                        :: PE_term !< terms in vertically averaged energy equation [R Z ~> kg m-2]
   real                        :: W0      !< A vertical velocity magnitude [Z T-1 ~> m s-1]
   real                        :: gp_unscaled !< A version of gprime rescaled to [L T-2 ~> m s-2].
   real, dimension(SZK_(GV)-1) :: lam_z   !< product of eigen value and gprime(k); one value for each
                                          !< interface (excluding surface and bottom)
   real, dimension(SZK_(GV)-1) :: a_diag, b_diag, c_diag
                                          !< diagonals of tridiagonal matrix; one value for each
                                          !< interface (excluding surface and bottom)
-  real, dimension(SZK_(GV)-1) :: e_guess !< guess at eigen vector with unit amplitde (for TDMA)
-  real, dimension(SZK_(GV)-1) :: e_itt   !< improved guess at eigen vector (from TDMA)
+  real, dimension(SZK_(GV)-1) :: e_guess !< guess at eigen vector with unit amplitude (for TDMA) [nondim]
+  real, dimension(SZK_(GV)-1) :: e_itt   !< improved guess at eigen vector (from TDMA) [nondim]
   real    :: Pi
   integer :: kc
   integer :: i, j, k, k2, itt, is, ie, js, je, nz, nzm, row, ig, jg, ig_stop, jg_stop
@@ -523,7 +523,7 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo
 
             ! Back-calculate amplitude from energy equation
             if (present(En) .and. (freq**2*Kmag2 > 0.0)) then
-              ! Units here are [R
+              ! Units here are [R Z ~> kg m-2]
               KE_term = 0.25*GV%Rho0*( ((freq**2 + f2) / (freq**2*Kmag2))*int_dwdz2 + int_w2 )
               PE_term = 0.25*GV%Rho0*( int_N2w2 / freq**2 )
               if (En(i,j) >= 0.0) then

src/initialization/MOM_shared_initialization.F90

@@ -814,7 +814,7 @@ subroutine reset_face_lengths_list(G, param_file, US)
   real, allocatable, dimension(:) :: &
     Dmin_u, Dmax_u, Davg_u   ! Porous barrier monomial fit params [m]
   real, allocatable, dimension(:) :: &
-    Dmin_v, Dmax_v, Davg_v
+    Dmin_v, Dmax_v, Davg_v   ! Porous barrier monomial fit params [m]
   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.