diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90 index 5ead019717..30b2e90d1a 100644 --- a/src/core/MOM_PressureForce_FV.F90 +++ b/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. 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_continuity_PPM.F90 b/src/core/MOM_continuity_PPM.F90 index 17d2f830c0..7bc67e2fdf 100644 --- a/src/core/MOM_continuity_PPM.F90 +++ b/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 diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index a762da7f33..68b844562f 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/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 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/diagnostics/MOM_wave_structure.F90 b/src/diagnostics/MOM_wave_structure.F90 index 2dd272d409..c833e973c5 100644 --- a/src/diagnostics/MOM_wave_structure.F90 +++ b/src/diagnostics/MOM_wave_structure.F90 @@ -171,11 +171,11 @@ 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 @@ -183,8 +183,8 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo 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 diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90 index bb5a84033b..fc5ceaf3e4 100644 --- a/src/initialization/MOM_shared_initialization.F90 +++ b/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. 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_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 index 046329523d..dd160c300c 100644 --- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 +++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 @@ -489,7 +489,7 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C ! netMassInOut = water [H ~> m or kg m-2] added/removed via surface fluxes ! netMassOut = water [H ~> m or kg m-2] removed via evaporating surface fluxes ! net_heat = heat via surface fluxes [degC H ~> degC m or degC kg m-2] - ! net_salt = salt via surface fluxes [ppt H ~> dppt m or gSalt m-2] + ! net_salt = salt via surface fluxes [ppt H ~> ppt m or gSalt m-2] ! Pen_SW_bnd = components to penetrative shortwave radiation call extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & CS%H_limit_fluxes, CS%use_river_heat_content, CS%use_calving_heat_content, & @@ -1527,7 +1527,7 @@ subroutine mechanical_entrainment(h, d_eb, htot, Ttot, Stot, uhtot, vhtot, & real :: TKE_full_ent ! The TKE remaining if a layer is fully entrained ! [Z L2 T-2 ~> m3 s-2]. real :: dRL ! Work required to mix water from the next layer - ! across the mixed layer [L2 T-2 ~> L2 s-2]. + ! across the mixed layer [L2 T-2 ~> m2 s-2]. real :: Pen_En_Contrib ! Penetrating SW contributions to the changes in ! TKE, divided by layer thickness in m [L2 T-2 ~> m2 s-2]. real :: Cpen1 ! A temporary variable [L2 T-2 ~> m2 s-2]. diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 312d114dde..c421b3a0f7 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -1171,7 +1171,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t ! netMassOut < 0 means mass leaves ocean. ! netHeat = heat via surface fluxes [degC H ~> degC m or degC kg m-2], excluding the part ! contained in Pen_SW_bnd; and excluding heat_content of netMassOut < 0. - ! netSalt = surface salt fluxes [ppt H ~> dppt m or gSalt m-2] + ! netSalt = surface salt fluxes [ppt H ~> ppt m or gSalt m-2] ! Pen_SW_bnd = components to penetrative shortwave radiation split according to bands. ! This field provides that portion of SW from atmosphere that in fact ! enters to the ocean and participates in pentrative SW heating. diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 861a8957c1..f7d4b0cc0d 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -200,7 +200,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) real :: hwtot ! Sum of the thicknesses used to calculate ! the near-bottom velocity magnitude [H ~> m or kg m-2]. real :: hutot ! Running sum of thicknesses times the - ! velocity magnitudes [H T T-1 ~> m2 s-1 or kg m-1 s-1]. + ! velocity magnitudes [H L T-1 ~> m2 s-1 or kg m-1 s-1]. real :: Thtot ! Running sum of thickness times temperature [degC H ~> degC m or degC kg m-2]. real :: Shtot ! Running sum of thickness times salinity [ppt H ~> ppt m or ppt kg m-2]. real :: hweight ! The thickness of a layer that is within Hbbl @@ -597,7 +597,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) ! When stratification dominates h_N< kg m-2 or kg m-5] + ustarsq = Rho0x400_G * ustar(i)**2 ! Note not in units of u*^2 but [H R ~> kg m-2 or kg2 m-5] htot = 0.0 ! Calculate the thickness of a stratification limited BBL ignoring rotation: 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 diff --git a/src/tracer/MOM_lateral_boundary_diffusion.F90 b/src/tracer/MOM_lateral_boundary_diffusion.F90 index 19d40f2db1..c11bc9856c 100644 --- a/src/tracer/MOM_lateral_boundary_diffusion.F90 +++ b/src/tracer/MOM_lateral_boundary_diffusion.F90 @@ -596,7 +596,7 @@ subroutine fluxes_layer_method(boundary, ke, hbl_L, hbl_R, h_L, h_R, phi_L, phi_ real, allocatable :: F_layer_z(:) !< Diffusive flux at U/V-point in the ztop grid [H L2 conc ~> m3 conc] real :: h_vel(ke) !< Thicknesses at u- and v-points in the native grid !! The harmonic mean is used to avoid zero values [H ~> m or kg m-2] - real :: khtr_avg !< Thickness-weighted diffusivity at the velocity-point [L2 T-1 ~> m s-1] + real :: khtr_avg !< Thickness-weighted diffusivity at the velocity-point [L2 T-1 ~> m2 s-1] !! This is just to remind developers that khtr_avg should be !! computed once khtr is 3D. real :: htot !< Total column thickness [H ~> m or kg m-2] diff --git a/src/tracer/MOM_tracer_registry.F90 b/src/tracer/MOM_tracer_registry.F90 index c729231927..2c77df3e74 100644 --- a/src/tracer/MOM_tracer_registry.F90 +++ b/src/tracer/MOM_tracer_registry.F90 @@ -84,7 +84,8 @@ module MOM_tracer_registry ! real, dimension(:,:,:), pointer :: diff_cont_xy => NULL() !< convergence of lateral diffusive tracer fluxes ! !! [conc H T-1 ~> conc m s-1 or conc kg m-2 s-1] ! real, dimension(:,:,:), pointer :: diff_conc_xy => NULL() !< convergence of lateral diffusive tracer fluxes -! !! expressed as a change in concentration [conc T-1] +! !! expressed as a change in concentration +! !! [conc T-1 ~> conc s-1] real, dimension(:,:,:), pointer :: t_prev => NULL() !< tracer concentration array at a previous !! timestep used for diagnostics [conc] real, dimension(:,:,:), pointer :: Trxh_prev => NULL() !< layer integrated tracer concentration array diff --git a/src/user/Idealized_Hurricane.F90 b/src/user/Idealized_Hurricane.F90 index 7182fc364a..707a0972f9 100644 --- a/src/user/Idealized_Hurricane.F90 +++ b/src/user/Idealized_Hurricane.F90 @@ -595,7 +595,7 @@ subroutine SCM_idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, C U_TS = CS%hurr_translation_spd*0.5*cos(transdir) V_TS = CS%hurr_translation_spd*0.5*sin(transdir) - ! Set the surface wind stresses, in [Pa]. A positive taux + ! Set the surface wind stresses, in [R L Z T-2 ~> Pa]. A positive taux ! accelerates the ocean to the (pseudo-)east. ! The i-loop extends to is-1 so that taux can be used later in the ! calculation of ustar - otherwise the lower bound would be Isq. diff --git a/src/user/MOM_controlled_forcing.F90 b/src/user/MOM_controlled_forcing.F90 index f783a271a6..d136d58a19 100644 --- a/src/user/MOM_controlled_forcing.F90 +++ b/src/user/MOM_controlled_forcing.F90 @@ -46,7 +46,7 @@ module MOM_controlled_forcing real :: lam_prec !< A constant of proportionality between SSS anomalies !! (normalised by mean SSS) and precipitation [R Z T-1 ~> kg m-2 s-1] real :: lam_cyc_heat !< A constant of proportionality between cyclical SST - !! anomalies and corrective heat fluxes [W m-2 degC-1] + !! anomalies and corrective heat fluxes [Q R Z T-1 degC-1 ~> W m-2 degC-1] real :: lam_cyc_prec !< A constant of proportionality between cyclical SSS !! anomalies (normalised by mean SSS) and corrective !! precipitation [R Z T-1 ~> kg m-2 s-1] @@ -270,7 +270,8 @@ subroutine apply_ctrl_forcing(SST_anom, SSS_anom, SSS_mean, virt_heat, virt_prec ! Accumulate the average anomalies for this period. dt_wt = wt_per1 * dt CS%avg_time(m_mid) = CS%avg_time(m_mid) + dt_wt - ! These loops temporarily change the units of the CS%avg_ variables to [degC s] or [ppt s]. + ! These loops temporarily change the units of the CS%avg_ variables to [degC T ~> degC s] + ! or [ppt T ~> ppt s]. do j=js,je ; do i=is,ie CS%avg_SST_anom(i,j,m_mid) = CS%avg_SST_anom(i,j,m_mid) + & dt_wt * G%mask2dT(i,j) * SST_anom(i,j) diff --git a/src/user/SCM_CVMix_tests.F90 b/src/user/SCM_CVMix_tests.F90 index 1fbc7a2b62..5bbe65b8d8 100644 --- a/src/user/SCM_CVMix_tests.F90 +++ b/src/user/SCM_CVMix_tests.F90 @@ -36,8 +36,8 @@ module SCM_CVMix_tests logical :: UseHeatFlux !< True to use heat flux logical :: UseEvaporation !< True to use evaporation logical :: UseDiurnalSW !< True to use diurnal sw radiation - real :: tau_x !< (Constant) Wind stress, X [Pa] - real :: tau_y !< (Constant) Wind stress, Y [Pa] + real :: tau_x !< (Constant) Wind stress, X [R L Z T-2 ~> Pa] + real :: tau_y !< (Constant) Wind stress, Y [R L Z T-2 ~> Pa] real :: surf_HF !< (Constant) Heat flux [degC Z T-1 ~> m degC s-1] real :: surf_evap !< (Constant) Evaporation rate [Z T-1 ~> m s-1] real :: Max_sw !< maximum of diurnal sw radiation [degC Z T-1 ~> degC m s-1] @@ -56,7 +56,7 @@ subroutine SCM_CVMix_tests_TS_init(T, S, h, G, GV, US, param_file, just_read) type(ocean_grid_type), intent(in) :: G !< Grid structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T !< Potential temperature [degC] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S !< Salinity [psu] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S !< Salinity [ppt] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2] type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(param_file_type), intent(in) :: param_file !< Input parameter structure diff --git a/src/user/baroclinic_zone_initialization.F90 b/src/user/baroclinic_zone_initialization.F90 index 1555f4ecad..a214012541 100644 --- a/src/user/baroclinic_zone_initialization.F90 +++ b/src/user/baroclinic_zone_initialization.F90 @@ -36,7 +36,7 @@ subroutine bcz_params(G, GV, US, param_file, S_ref, dSdz, delta_S, dSdx, T_ref, real, intent(out) :: S_ref !< Reference salinity [ppt] real, intent(out) :: dSdz !< Salinity stratification [ppt Z-1 ~> ppt m-1] real, intent(out) :: delta_S !< Salinity difference across baroclinic zone [ppt] - real, intent(out) :: dSdx !< Linear salinity gradient [ppt m-1] + real, intent(out) :: dSdx !< Linear salinity gradient [ppt G%xaxis_units-1] real, intent(out) :: T_ref !< Reference temperature [degC] real, intent(out) :: dTdz !< Temperature stratification [degC Z-1 ~> degC m-1] real, intent(out) :: delta_T !< Temperature difference across baroclinic zone [degC]