Skip to content

Commit

Permalink
+Rescale tides and ramp-up times
Browse files Browse the repository at this point in the history
  Rescaled the dimensions of the tidal amplitudes and frequencies used
internally in calc_tidal_forcing() and ramp-up times used by update_OBC_ramp()
and updateCFLtruncationValue() for nearly complete dimensional consistency
testing.  New unit_scale_type arguments were added to 5 routines, in the case of
calc_tidal_forcing() replacing a previous optional argument that was always
being used.  One overly short internal variable, "N", was renamed "nodelon" to
make its purpose clearer and easier to search for.  All answers are bitwise
identical, but there are changes to the argument lists of 5 routines.
  • Loading branch information
Hallberg-NOAA committed Dec 14, 2021
1 parent dab2eb5 commit 7b81cf2
Show file tree
Hide file tree
Showing 8 changed files with 127 additions and 112 deletions.
4 changes: 2 additions & 2 deletions src/core/MOM_PressureForce_FV.F90
Original file line number Diff line number Diff line change
Expand Up @@ -306,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)
Expand Down Expand Up @@ -574,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.
Expand Down
4 changes: 2 additions & 2 deletions src/core/MOM_PressureForce_Montgomery.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down Expand Up @@ -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.
Expand Down
10 changes: 5 additions & 5 deletions src/core/MOM_dynamics_split_RK2.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion src/core/MOM_dynamics_unsplit.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
2 changes: 1 addition & 1 deletion src/core/MOM_dynamics_unsplit_RK2.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
25 changes: 14 additions & 11 deletions src/core/MOM_open_boundary.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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].
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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).
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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

Expand All @@ -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
Expand Down
Loading

0 comments on commit 7b81cf2

Please sign in to comment.