diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90 index 1666b4a97e..30b2e90d1a 100644 --- a/src/core/MOM_PressureForce_FV.F90 +++ b/src/core/MOM_PressureForce_FV.F90 @@ -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) @@ -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. diff --git a/src/core/MOM_PressureForce_Montgomery.F90 b/src/core/MOM_PressureForce_Montgomery.F90 index a827fb12d0..18ea07b313 100644 --- a/src/core/MOM_PressureForce_Montgomery.F90 +++ b/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. diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index f1a2b2469d..68b844562f 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -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 diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index 88a11e071c..fcc4c3d49b 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/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) diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index 26bd00aaf5..694d88f2ea 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/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) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 6d8696216a..2c3f016005 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/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 diff --git a/src/parameterizations/lateral/MOM_tidal_forcing.F90 b/src/parameterizations/lateral/MOM_tidal_forcing.F90 index b8d5c44098..cc4517a473 100644 --- a/src/parameterizations/lateral/MOM_tidal_forcing.F90 +++ b/src/parameterizations/lateral/MOM_tidal_forcing.F90 @@ -11,6 +11,7 @@ module MOM_tidal_forcing use MOM_grid, only : ocean_grid_type use MOM_io, only : field_exists, file_exists, MOM_read_data use MOM_time_manager, only : set_date, time_type, time_type_to_real, operator(-) +use MOM_unit_scaling, only : unit_scale_type implicit none ; private @@ -47,12 +48,12 @@ module MOM_tidal_forcing !! astronomical/equilibrium argument. real :: sal_scalar !< The constant of proportionality between sea surface !! height (really it should be bottom pressure) anomalies - !! and bottom geopotential anomalies. + !! and bottom geopotential anomalies [nondim]. integer :: nc !< The number of tidal constituents in use. real, dimension(MAX_CONSTITUENTS) :: & - freq, & !< The frequency of a tidal constituent [s-1]. - phase0, & !< The phase of a tidal constituent at time 0, in radians. - amp, & !< The amplitude of a tidal constituent at time 0 [m]. + freq, & !< The frequency of a tidal constituent [T-1 ~> s-1]. + phase0, & !< The phase of a tidal constituent at time 0 [rad]. + amp, & !< The amplitude of a tidal constituent at time 0 [Z ~> m]. love_no !< The Love number of a tidal constituent at time 0 [nondim]. integer :: struct(MAX_CONSTITUENTS) !< An encoded spatial structure for each constituent character (len=16) :: const_name(MAX_CONSTITUENTS) !< The name of each constituent @@ -62,13 +63,13 @@ module MOM_tidal_forcing !! tidal phases at t = 0. real, allocatable :: & sin_struct(:,:,:), & !< The sine and cosine based structures that can - cos_struct(:,:,:), & !< be associated with the astronomical forcing. + cos_struct(:,:,:), & !< be associated with the astronomical forcing [nondim]. cosphasesal(:,:,:), & !< The cosine and sine of the phase of the sinphasesal(:,:,:), & !< self-attraction and loading amphidromes. - ampsal(:,:,:), & !< The amplitude of the SAL [m]. + ampsal(:,:,:), & !< The amplitude of the SAL [Z ~> m]. cosphase_prev(:,:,:), & !< The cosine and sine of the phase of the sinphase_prev(:,:,:), & !< amphidromes in the previous tidal solutions. - amp_prev(:,:,:) !< The amplitude of the previous tidal solution [m]. + amp_prev(:,:,:) !< The amplitude of the previous tidal solution [Z ~> m]. end type tidal_forcing_CS integer :: id_clock_tides !< CPU clock for tides @@ -87,8 +88,9 @@ module MOM_tidal_forcing subroutine astro_longitudes_init(time_ref, longitudes) type(time_type), intent(in) :: time_ref !> Time to calculate longitudes for. type(astro_longitudes), intent(out) :: longitudes !> Lunar and solar longitudes at time_ref. - real :: D, T !> Date offsets - real, parameter :: PI = 4.0 * atan(1.0) !> 3.14159... + real :: D !> Time since the reference date [days] + real :: T !> Time in Julian centuries [centuries] + real, parameter :: PI = 4.0 * atan(1.0) !> 3.14159... [nondim] ! Find date at time_ref in days since 1900-01-01 D = time_type_to_real(time_ref - set_date(1900, 1, 1)) / (24.0 * 3600.0) ! Time since 1900-01-01 in Julian centuries @@ -176,44 +178,45 @@ end function tidal_frequency !> Find amplitude (f) and phase (u) modulation of tidal constituents by the 18.6 !! year nodal cycle. Values here follow Table I.6 in Kowalik and Luick, !! "Modern Theory and Practice of Tide Analysis and Tidal Power", 2019. -subroutine nodal_fu(constit, N, fn, un) - character (len=2), intent(in) :: constit !> Tidal constituent to find modulation for. - real, intent(in) :: N !> Longitude of ascending node [rad]. - !! Calculate using astro_longitudes_init. - real, parameter :: RADIANS = 4.0 * atan(1.0) / 180.0 !> Converts degrees to radians. - real, intent(out) :: & - fn, & !> Amplitude modulation [nondim] - un !> Phase modulation [rad] +subroutine nodal_fu(constit, nodelon, fn, un) + character (len=2), intent(in) :: constit !> Tidal constituent to find modulation for. + real, intent(in) :: nodelon !> Longitude of ascending node [rad], which + !! can be calculated using astro_longitudes_init. + real, intent(out) :: fn !> Amplitude modulation [nondim] + real, intent(out) :: un !> Phase modulation [rad] + + real, parameter :: RADIANS = 4.0 * atan(1.0) / 180.0 !> Converts degrees to radians [nondim] + select case (constit) case ("M2") - fn = 1.0 - 0.037 * cos(N) - un = -2.1 * RADIANS * sin(N) + fn = 1.0 - 0.037 * cos(nodelon) + un = -2.1 * RADIANS * sin(nodelon) case ("S2") fn = 1.0 ! Solar S2 has no amplitude modulation. un = 0.0 ! S2 has no phase modulation. case ("N2") - fn = 1.0 - 0.037 * cos(N) - un = -2.1 * RADIANS * sin(N) + fn = 1.0 - 0.037 * cos(nodelon) + un = -2.1 * RADIANS * sin(nodelon) case ("K2") - fn = 1.024 + 0.286 * cos(N) - un = -17.7 * RADIANS * sin(N) + fn = 1.024 + 0.286 * cos(nodelon) + un = -17.7 * RADIANS * sin(nodelon) case ("K1") - fn = 1.006 + 0.115 * cos(N) - un = -8.9 * RADIANS * sin(N) + fn = 1.006 + 0.115 * cos(nodelon) + un = -8.9 * RADIANS * sin(nodelon) case ("O1") - fn = 1.009 + 0.187 * cos(N) - un = 10.8 * RADIANS * sin(N) + fn = 1.009 + 0.187 * cos(nodelon) + un = 10.8 * RADIANS * sin(nodelon) case ("P1") fn = 1.0 ! P1 has no amplitude modulation. un = 0.0 ! P1 has no phase modulation. case ("Q1") - fn = 1.009 + 0.187 * cos(N) - un = 10.8 * RADIANS * sin(N) + fn = 1.009 + 0.187 * cos(nodelon) + un = 10.8 * RADIANS * sin(nodelon) case ("MF") - fn = 1.043 + 0.414 * cos(N) - un = -23.7 * RADIANS * sin(N) + fn = 1.043 + 0.414 * cos(nodelon) + un = -23.7 * RADIANS * sin(nodelon) case ("MM") - fn = 1.0 - 0.130 * cos(N) + fn = 1.0 - 0.130 * cos(nodelon) un = 0.0 ! MM has no phase modulation. case default call MOM_error(FATAL, "nodal_fu: unrecognized constituent") @@ -226,10 +229,11 @@ end subroutine nodal_fu !! while fields like the background viscosities are 2-D arrays. !! ALLOC is a macro defined in MOM_memory.h for allocate or nothing with !! static memory. -subroutine tidal_forcing_init(Time, G, param_file, CS) - type(time_type), intent(in) :: Time !< The current model time. - type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. - type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters. +subroutine tidal_forcing_init(Time, G, US, param_file, CS) + type(time_type), intent(in) :: Time !< The current model time. + 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(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters. type(tidal_forcing_CS), intent(inout) :: CS !< Tidal forcing control struct ! Local variables @@ -237,15 +241,18 @@ subroutine tidal_forcing_init(Time, G, param_file, CS) phase, & ! The phase of some tidal constituent. lat_rad, lon_rad ! Latitudes and longitudes of h-points in radians. real :: deg_to_rad - real, dimension(MAX_CONSTITUENTS) :: freq_def, phase0_def, amp_def, love_def + real, dimension(MAX_CONSTITUENTS) :: freq_def ! Default frequency for each tidal constituent [s-1] + real, dimension(MAX_CONSTITUENTS) :: phase0_def ! Default reference phase for each tidal constituent [rad] + real, dimension(MAX_CONSTITUENTS) :: amp_def ! Default amplitude for each tidal constituent [m] + real, dimension(MAX_CONSTITUENTS) :: love_def ! Default love number for each constituent [nondim] integer, dimension(3) :: tide_ref_date !< Reference date (t = 0) for tidal forcing. logical :: use_const ! True if a constituent is being used. logical :: use_M2, use_S2, use_N2, use_K2, use_K1, use_O1, use_P1, use_Q1 logical :: use_MF, use_MM logical :: tides ! True if a tidal forcing is to be used. logical :: FAIL_IF_MISSING = .true. -! 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_tidal_forcing" ! This module's name. character(len=128) :: mesg character(len=200) :: tidal_input_files(4*MAX_CONSTITUENTS) @@ -389,68 +396,68 @@ subroutine tidal_forcing_init(Time, G, param_file, CS) endif CS%time_ref = set_date(tide_ref_date(1), tide_ref_date(2), tide_ref_date(3)) endif - ! Set the parameters for all components that are in use. - ! Initialize reference time for tides and - ! find relevant lunar and solar longitudes at the reference time. + + ! Initialize reference time for tides and find relevant lunar and solar + ! longitudes at the reference time. if (CS%use_eq_phase) call astro_longitudes_init(CS%time_ref, CS%tidal_longitudes) + + ! Set the parameters for all components that are in use. c=0 if (use_M2) then c=c+1 ; CS%const_name(c) = "M2" ; CS%struct(c) = 2 - CS%love_no(c) = 0.693 ; CS%amp(c) = 0.242334 + CS%love_no(c) = 0.693 ; amp_def(c) = 0.242334 ! Default amplitude in m. endif if (use_S2) then c=c+1 ; CS%const_name(c) = "S2" ; CS%struct(c) = 2 - CS%love_no(c) = 0.693 ; CS%amp(c) = 0.112743 + CS%love_no(c) = 0.693 ; amp_def(c) = 0.112743 ! Default amplitude in m. endif if (use_N2) then c=c+1 ; CS%const_name(c) = "N2" ; CS%struct(c) = 2 - CS%love_no(c) = 0.693 ; CS%amp(c) = 0.046397 + CS%love_no(c) = 0.693 ; amp_def(c) = 0.046397 ! Default amplitude in m. endif if (use_K2) then c=c+1 ; CS%const_name(c) = "K2" ; CS%struct(c) = 2 - CS%love_no(c) = 0.693 ; CS%amp(c) = 0.030684 + CS%love_no(c) = 0.693 ; amp_def(c) = 0.030684 ! Default amplitude in m. endif if (use_K1) then c=c+1 ; CS%const_name(c) = "K1" ; CS%struct(c) = 1 - CS%love_no(c) = 0.736 ; CS%amp(c) = 0.141565 + CS%love_no(c) = 0.736 ; amp_def(c) = 0.141565 ! Default amplitude in m. endif if (use_O1) then c=c+1 ; CS%const_name(c) = "O1" ; CS%struct(c) = 1 - CS%love_no(c) = 0.695 ; CS%amp(c) = 0.100661 + CS%love_no(c) = 0.695 ; amp_def(c) = 0.100661 ! Default amplitude in m. endif if (use_P1) then c=c+1 ; CS%const_name(c) = "P1" ; CS%struct(c) = 1 - CS%love_no(c) = 0.706 ; CS%amp(c) = 0.046848 + CS%love_no(c) = 0.706 ; amp_def(c) = 0.046848 ! Default amplitude in m. endif if (use_Q1) then c=c+1 ; CS%const_name(c) = "Q1" ; CS%struct(c) = 1 - CS%love_no(c) = 0.695 ; CS%amp(c) = 0.019273 + CS%love_no(c) = 0.695 ; amp_def(c) = 0.019273 ! Default amplitude in m. endif if (use_MF) then c=c+1 ; CS%const_name(c) = "MF" ; CS%struct(c) = 3 - CS%love_no(c) = 0.693 ; CS%amp(c) = 0.042041 + CS%love_no(c) = 0.693 ; amp_def(c) = 0.042041 ! Default amplitude in m. endif if (use_MM) then c=c+1 ; CS%const_name(c) = "MM" ; CS%struct(c) = 3 - CS%love_no(c) = 0.693 ; CS%amp(c) = 0.022191 + CS%love_no(c) = 0.693 ; amp_def(c) = 0.022191 ! Default amplitude in m. endif ! Set defaults for all included constituents ! and things that can be set by functions do c=1,nc - CS%freq(c) = tidal_frequency(CS%const_name(c)) - freq_def(c) = CS%freq(c) + freq_def(c) = tidal_frequency(CS%const_name(c)) love_def(c) = CS%love_no(c) - amp_def(c) = CS%amp(c) CS%phase0(c) = 0.0 if (CS%use_eq_phase) then phase0_def(c) = eq_phase(CS%const_name(c), CS%tidal_longitudes) @@ -467,11 +474,11 @@ subroutine tidal_forcing_init(Time, G, param_file, CS) "Frequency of the "//trim(CS%const_name(c))//" tidal constituent. "//& "This is only used if TIDES and TIDE_"//trim(CS%const_name(c))// & " are true, or if OBC_TIDE_N_CONSTITUENTS > 0 and "//trim(CS%const_name(c))// & - " is in OBC_TIDE_CONSTITUENTS.", units="s-1", default=freq_def(c)) + " is in OBC_TIDE_CONSTITUENTS.", units="s-1", default=freq_def(c), scale=US%T_to_s) call get_param(param_file, mdl, "TIDE_"//trim(CS%const_name(c))//"_AMP", CS%amp(c), & "Amplitude of the "//trim(CS%const_name(c))//" tidal constituent. "//& "This is only used if TIDES and TIDE_"//trim(CS%const_name(c))// & - " are true.", units="m", default=amp_def(c)) + " are true.", units="m", default=amp_def(c), scale=US%m_to_Z) call get_param(param_file, mdl, "TIDE_"//trim(CS%const_name(c))//"_PHASE_T0", CS%phase0(c), & "Phase of the "//trim(CS%const_name(c))//" tidal constituent at time 0. "//& "This is only used if TIDES and TIDE_"//trim(CS%const_name(c))// & @@ -484,8 +491,9 @@ subroutine tidal_forcing_init(Time, G, param_file, CS) allocate(CS%ampsal(isd:ied,jsd:jed,nc)) do c=1,nc ! Read variables with names like PHASE_SAL_M2 and AMP_SAL_M2. - call find_in_files(tidal_input_files,"PHASE_SAL_"//trim(CS%const_name(c)),phase,G) - call find_in_files(tidal_input_files,"AMP_SAL_"//trim(CS%const_name(c)),CS%ampsal(:,:,c),G) + call find_in_files(tidal_input_files, "PHASE_SAL_"//trim(CS%const_name(c)), phase, G) + call find_in_files(tidal_input_files, "AMP_SAL_"//trim(CS%const_name(c)), CS%ampsal(:,:,c), & + G, scale=US%m_to_Z) call pass_var(phase, G%domain,complete=.false.) call pass_var(CS%ampsal(:,:,c),G%domain,complete=.true.) do j=js-1,je+1 ; do i=is-1,ie+1 @@ -501,8 +509,9 @@ subroutine tidal_forcing_init(Time, G, param_file, CS) allocate(CS%amp_prev(isd:ied,jsd:jed,nc)) do c=1,nc ! Read variables with names like PHASE_PREV_M2 and AMP_PREV_M2. - call find_in_files(tidal_input_files,"PHASE_PREV_"//trim(CS%const_name(c)),phase,G) - call find_in_files(tidal_input_files,"AMP_PREV_"//trim(CS%const_name(c)),CS%amp_prev(:,:,c),G) + call find_in_files(tidal_input_files, "PHASE_PREV_"//trim(CS%const_name(c)), phase, G) + call find_in_files(tidal_input_files, "AMP_PREV_"//trim(CS%const_name(c)), CS%amp_prev(:,:,c), & + G, scale=US%m_to_Z) call pass_var(phase, G%domain,complete=.false.) call pass_var(CS%amp_prev(:,:,c),G%domain,complete=.true.) do j=js-1,je+1 ; do i=is-1,ie+1 @@ -518,18 +527,19 @@ end subroutine tidal_forcing_init !> This subroutine finds a named variable in a list of files and reads its !! values into a domain-decomposed 2-d array -subroutine find_in_files(filenames, varname, array, G) +subroutine find_in_files(filenames, varname, array, G, scale) character(len=*), dimension(:), intent(in) :: filenames !< The names of the files to search for the named variable character(len=*), intent(in) :: varname !< The name of the variable to read type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure real, dimension(SZI_(G),SZJ_(G)), intent(out) :: array !< The array to fill with the data + real, optional, intent(in) :: scale !< A factor by which to rescale the array. ! Local variables integer :: nf do nf=1,size(filenames) if (LEN_TRIM(filenames(nf)) == 0) cycle if (field_exists(filenames(nf), varname, MOM_domain=G%Domain)) then - call MOM_read_data(filenames(nf), varname, array, G%Domain) + call MOM_read_data(filenames(nf), varname, array, G%Domain, scale=scale) return endif enddo @@ -571,22 +581,22 @@ end subroutine tidal_forcing_sensitivity !! height. For now, eta and eta_tidal are both geopotential heights in depth !! units, but probably the input for eta should really be replaced with the !! column mass anomalies. -subroutine calc_tidal_forcing(Time, eta, eta_tidal, G, CS, m_to_Z) +subroutine calc_tidal_forcing(Time, eta, eta_tidal, G, US, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(time_type), intent(in) :: Time !< The time for the caluculation. real, dimension(SZI_(G),SZJ_(G)), intent(in) :: eta !< The sea surface height anomaly from !! a time-mean geoid [Z ~> m]. real, dimension(SZI_(G),SZJ_(G)), intent(out) :: eta_tidal !< The tidal forcing geopotential height !! anomalies [Z ~> m]. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(tidal_forcing_CS), intent(in) :: CS !< The control structure returned by a !! previous call to tidal_forcing_init. - real, intent(in) :: m_to_Z !< A scaling factor from m to the units of eta. ! Local variables - real :: now ! The relative time in seconds. - real :: amp_cosomegat, amp_sinomegat - real :: cosomegat, sinomegat - real :: eta_prop ! The nondimenional constant of proportionality beteen eta and eta_tidal. + real :: now ! The relative time compared with the tidal reference [T ~> s] + real :: amp_cosomegat, amp_sinomegat ! The tidal amplitudes times the components of phase [Z ~> m] + real :: cosomegat, sinomegat ! The components of the phase [nondim] + real :: eta_prop ! The nondimenional constant of proportionality beteen eta and eta_tidal [nondim] integer :: i, j, c, m, is, ie, js, je, Isq, Ieq, Jsq, Jeq is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB @@ -598,7 +608,7 @@ subroutine calc_tidal_forcing(Time, eta, eta_tidal, G, CS, m_to_Z) return endif - now = time_type_to_real(Time - cs%time_ref) + now = US%s_to_T * time_type_to_real(Time - cs%time_ref) if (CS%USE_SAL_SCALAR .and. CS%USE_PREV_TIDES) then eta_prop = 2.0*CS%SAL_SCALAR @@ -614,8 +624,8 @@ subroutine calc_tidal_forcing(Time, eta, eta_tidal, G, CS, m_to_Z) do c=1,CS%nc m = CS%struct(c) - amp_cosomegat = m_to_Z*CS%amp(c)*CS%love_no(c) * cos(CS%freq(c)*now + CS%phase0(c)) - amp_sinomegat = m_to_Z*CS%amp(c)*CS%love_no(c) * sin(CS%freq(c)*now + CS%phase0(c)) + amp_cosomegat = CS%amp(c)*CS%love_no(c) * cos(CS%freq(c)*now + CS%phase0(c)) + amp_sinomegat = CS%amp(c)*CS%love_no(c) * sin(CS%freq(c)*now + CS%phase0(c)) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 eta_tidal(i,j) = eta_tidal(i,j) + (amp_cosomegat*CS%cos_struct(i,j,m) + & amp_sinomegat*CS%sin_struct(i,j,m)) @@ -626,7 +636,7 @@ subroutine calc_tidal_forcing(Time, eta, eta_tidal, G, CS, m_to_Z) cosomegat = cos(CS%freq(c)*now) sinomegat = sin(CS%freq(c)*now) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - eta_tidal(i,j) = eta_tidal(i,j) + m_to_Z*CS%ampsal(i,j,c) * & + eta_tidal(i,j) = eta_tidal(i,j) + CS%ampsal(i,j,c) * & (cosomegat*CS%cosphasesal(i,j,c) + sinomegat*CS%sinphasesal(i,j,c)) enddo ; enddo enddo ; endif @@ -635,7 +645,7 @@ subroutine calc_tidal_forcing(Time, eta, eta_tidal, G, CS, m_to_Z) cosomegat = cos(CS%freq(c)*now) sinomegat = sin(CS%freq(c)*now) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - eta_tidal(i,j) = eta_tidal(i,j) - m_to_Z*CS%SAL_SCALAR*CS%amp_prev(i,j,c) * & + eta_tidal(i,j) = eta_tidal(i,j) - CS%SAL_SCALAR*CS%amp_prev(i,j,c) * & (cosomegat*CS%cosphase_prev(i,j,c) + sinomegat*CS%sinphase_prev(i,j,c)) enddo ; enddo enddo ; endif diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index adac9e83f4..d384500c3d 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -56,12 +56,12 @@ module MOM_vert_friction !! absolute velocities. real :: CFL_trunc !< Velocity components will be truncated when they !! are large enough that the corresponding CFL number - !! exceeds this value, nondim. + !! exceeds this value [nondim]. real :: CFL_report !< The value of the CFL number that will cause the - !! accelerations to be reported, nondim. CFL_report + !! accelerations to be reported [nondim]. CFL_report !! will often equal CFL_trunc. real :: truncRampTime !< The time-scale over which to ramp up the value of - !! CFL_trunc from CFL_truncS to CFL_truncE + !! CFL_trunc from CFL_truncS to CFL_truncE [T ~> s] real :: CFL_truncS !< The start value of CFL_trunc real :: CFL_truncE !< The end/target value of CFL_trunc logical :: CFLrampingIsActivated = .false. !< True if the ramping has been initialized @@ -105,7 +105,7 @@ module MOM_vert_friction !! thickness for viscosity. logical :: answers_2018 !< If true, use the order of arithmetic and expressions that recover the !! answers from the end of 2018. Otherwise, use expressions that do not - !! use an arbitary and hard-coded maximum viscous coupling coefficient + !! use an arbitrary and hard-coded maximum viscous coupling coefficient !! between layers. logical :: debug !< If true, write verbose checksums for debugging purposes. integer :: nkml !< The number of layers in the mixed layer. @@ -533,7 +533,7 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & endif if (associated(ADp%du_dt_str) .and. associated(ADp%dv_dt_str)) then - ! Diagnostics for thickness x wind stress acclerations + ! Diagnostics for thickness x wind stress accelerations if (CS%id_h_du_dt_str > 0) call post_product_u(CS%id_h_du_dt_str, ADp%du_dt_str, ADp%diag_hu, G, nz, CS%diag) if (CS%id_h_dv_dt_str > 0) call post_product_v(CS%id_h_dv_dt_str, ADp%dv_dt_str, ADp%diag_hv, G, nz, CS%diag) @@ -555,11 +555,11 @@ subroutine vertvisc_remnant(visc, visc_rem_u, visc_rem_v, dt, G, GV, US, CS) type(vertvisc_type), intent(in) :: visc !< Viscosities and bottom drag real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & intent(inout) :: visc_rem_u !< Fraction of a time-step's worth of a - !! barotopic acceleration that a layer experiences after + !! barotropic acceleration that a layer experiences after !! viscosity is applied in the zonal direction [nondim] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & intent(inout) :: visc_rem_v !< Fraction of a time-step's worth of a - !! barotopic acceleration that a layer experiences after + !! barotropic acceleration that a layer experiences after !! viscosity is applied in the meridional direction [nondim] real, intent(in) :: dt !< Time increment [T ~> s] type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -692,7 +692,7 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC) a_shelf, & ! The drag coefficients across interfaces in water columns under ! ice shelves [Z T-1 ~> m s-1]. z_i ! An estimate of each interface's height above the bottom, - ! normalized by the bottom boundary layer thickness, nondim. + ! normalized by the bottom boundary layer thickness [nondim] real, dimension(SZIB_(G)) :: & kv_bbl, & ! The bottom boundary layer viscosity [Z2 T-1 ~> m2 s-1]. bbl_thick, & ! The bottom boundary layer thickness [H ~> m or kg m-2]. @@ -715,10 +715,10 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC) ! than Hbbl into the interior. real :: topfn ! A function which goes from 1 at the top to 0 much more ! than Htbl into the interior. - real :: z2 ! The distance from the bottom, normalized by Hbbl, nondim. + real :: z2 ! The distance from the bottom, normalized by Hbbl [nondim] real :: z2_wt ! A nondimensional (0-1) weight used when calculating z2. real :: z_clear ! The clearance of an interface above the surrounding topography [H ~> m or kg m-2]. - real :: a_cpl_max ! The maximum drag doefficient across interfaces, set so that it will be + real :: a_cpl_max ! The maximum drag coefficient across interfaces, set so that it will be ! representable as a 32-bit float in MKS units [Z T-1 ~> m s-1] 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]. @@ -1193,7 +1193,7 @@ subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, h_neglect = GV%H_subroundoff if (CS%answers_2018) then - ! The maximum coupling coefficent was originally introduced to avoid + ! The maximum coupling coefficient was originally introduced to avoid ! truncation error problems in the tridiagonal solver. Effectively, the 1e-10 ! sets the maximum coupling coefficient increment to 1e10 m per timestep. I_amax = (1.0e-10*US%Z_to_m) * dt @@ -1759,7 +1759,7 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & call get_param(param_file, mdl, "CFL_TRUNCATE_RAMP_TIME", CS%truncRampTime, & "The time over which the CFL truncation value is ramped "//& "up at the beginning of the run.", & - units="s", default=0.) + units="s", default=0., scale=US%s_to_T) CS%CFL_truncE = CS%CFL_trunc call get_param(param_file, mdl, "CFL_TRUNCATE_START", CS%CFL_truncS, & "The start value of the truncation CFL number used when "//& @@ -1937,14 +1937,16 @@ end subroutine vertvisc_init !> Update the CFL truncation 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 updateCFLtruncationValue(Time, CS, activate) +subroutine updateCFLtruncationValue(Time, CS, US, activate) type(time_type), target, intent(in) :: Time !< Current model time type(vertvisc_CS), pointer :: CS !< Vertical viscosity control 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, wghtA + real :: deltaTime ! The time since CS%rampStartTime [T ~> s], which may be negative. + real :: wghtA ! The relative weight of the final value [nondim] character(len=12) :: msg if (CS%truncRampTime==0.) return ! This indicates to ramping is turned off @@ -1958,7 +1960,7 @@ subroutine updateCFLtruncationValue(Time, CS, activate) endif endif if (.not.CS%CFLrampingIsActivated) return - deltaTime = max( 0., time_type_to_real( Time - CS%rampStartTime ) ) + deltaTime = max( 0., US%s_to_T*time_type_to_real( Time - CS%rampStartTime ) ) if (deltaTime >= CS%truncRampTime) then CS%CFL_trunc = CS%CFL_truncE CS%truncRampTime = 0. ! This turns off ramping after this call @@ -1966,7 +1968,7 @@ subroutine updateCFLtruncationValue(Time, CS, activate) wghtA = min( 1., deltaTime / CS%truncRampTime ) ! Linear profile in time !wghtA = wghtA*wghtA ! Convert linear profile to parabolic profile in time !wghtA = wghtA*wghtA*(3. - 2.*wghtA) ! Convert linear profile to cosine profile - wghtA = 1. - ( (1. - wghtA)**2 ) ! Convert linear profiel to nverted parabolic profile + wghtA = 1. - ( (1. - wghtA)**2 ) ! Convert linear profile to inverted parabolic profile CS%CFL_trunc = CS%CFL_truncS + wghtA * ( CS%CFL_truncE - CS%CFL_truncS ) endif write(msg(1:12),'(es12.3)') CS%CFL_trunc