From 10d9523e11546b76f2e00ea504dae312e8285768 Mon Sep 17 00:00:00 2001 From: He Wang Date: Mon, 22 Jul 2024 14:23:15 -0400 Subject: [PATCH 01/17] Refactor SAL in MOM_PressureForce_FV Self-attraction and loading calculation in Boussinesq pressure gradient force is refactored to be consistent with the algorithm in non-Boussinesq version. --- src/core/MOM_PressureForce_FV.F90 | 112 +++++++++++------------------- 1 file changed, 40 insertions(+), 72 deletions(-) diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90 index 42e6514ab9..905df97d68 100644 --- a/src/core/MOM_PressureForce_FV.F90 +++ b/src/core/MOM_PressureForce_FV.F90 @@ -1003,77 +1003,50 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm MassWt_u(:,:,:) = 0.0 ; MassWt_v(:,:,:) = 0.0 endif - if (CS%tides_answer_date>20230630) then - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - e(i,j,nz+1) = -G%bathyT(i,j) - enddo ; enddo + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + e(i,j,nz+1) = -G%bathyT(i,j) + enddo ; enddo - ! Calculate and add the self-attraction and loading geopotential anomaly. - if (CS%calculate_SAL) then - ! Determine the surface height anomaly for calculating self attraction - ! and loading. This should really be based on bottom pressure anomalies, - ! but that is not yet implemented, and the current form is correct for - ! barotropic tides. - !$OMP parallel do default(shared) - do j=Jsq,Jeq+1 - do i=Isq,Ieq+1 - SSH(i,j) = min(-G%bathyT(i,j) - G%Z_ref, 0.0) - enddo - do k=1,nz ; do i=Isq,Ieq+1 - SSH(i,j) = SSH(i,j) + h(i,j,k)*GV%H_to_Z - enddo ; enddo + ! Calculate and add the self-attraction and loading geopotential anomaly. + if (CS%calculate_SAL) then + ! Determine the surface height anomaly for calculating self attraction + ! and loading. This should really be based on bottom pressure anomalies, + ! but that is not yet implemented, and the current form is correct for + ! barotropic tides. + !$OMP parallel do default(shared) + do j=Jsq,Jeq+1 + do i=Isq,Ieq+1 + SSH(i,j) = min(-G%bathyT(i,j) - G%Z_ref, 0.0) enddo - call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m) + do k=1,nz ; do i=Isq,Ieq+1 + SSH(i,j) = SSH(i,j) + h(i,j,k)*GV%H_to_Z + enddo ; enddo + enddo + call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m) + + if (CS%tides_answer_date>20230630) then !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 e(i,j,nz+1) = e(i,j,nz+1) - e_sal(i,j) enddo ; enddo endif + endif - ! Calculate and add the tidal geopotential anomaly. - if (CS%tides) then + ! Calculate and add the tidal geopotential anomaly. + if (CS%tides) then + if (CS%tides_answer_date>20230630) then call calc_tidal_forcing(CS%Time, e_tide_eq, e_tide_sal, G, US, CS%tides_CSp) - !$OMP parallel do default(shared) + !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 e(i,j,nz+1) = e(i,j,nz+1) - (e_tide_eq(i,j) + e_tide_sal(i,j)) enddo ; enddo - endif - else ! Old answers - ! Calculate and add the self-attraction and loading geopotential anomaly. - if (CS%calculate_SAL) then - ! Determine the surface height anomaly for calculating self attraction - ! and loading. This should really be based on bottom pressure anomalies, - ! but that is not yet implemented, and the current form is correct for - ! barotropic tides. - !$OMP parallel do default(shared) - do j=Jsq,Jeq+1 - do i=Isq,Ieq+1 - SSH(i,j) = min(-G%bathyT(i,j) - G%Z_ref, 0.0) - enddo - do k=1,nz ; do i=Isq,Ieq+1 - SSH(i,j) = SSH(i,j) + h(i,j,k)*GV%H_to_Z - enddo ; enddo - enddo - call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m) - else - !$OMP parallel do default(shared) - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - e_sal(i,j) = 0.0 - enddo ; enddo - endif - - ! Calculate and add the tidal geopotential anomaly. - if (CS%tides) then + else ! Old answers + if (.not.CS%calculate_SAL) e_sal(:,:) = 0.0 call calc_tidal_forcing_legacy(CS%Time, e_sal, e_sal_tide, e_tide_eq, e_tide_sal, & G, US, CS%tides_CSp) !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - e(i,j,nz+1) = -(G%bathyT(i,j) + e_sal_tide(i,j)) - enddo ; enddo - else - !$OMP parallel do default(shared) - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - e(i,j,nz+1) = -(G%bathyT(i,j) + e_sal(i,j)) + e(i,j,nz+1) = e(i,j,nz+1) - e_sal_tide(i,j) enddo ; enddo endif endif @@ -1640,33 +1613,28 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm if (present(eta)) then ! eta is the sea surface height relative to a time-invariant geoid, for comparison with ! what is used for eta in btstep. See how e was calculated about 200 lines above. - if (CS%tides_answer_date>20230630) then - !$OMP parallel do default(shared) - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - eta(i,j) = e(i,j,1)*GV%Z_to_H - enddo ; enddo - if (CS%tides) then + !$OMP parallel do default(shared) + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + eta(i,j) = e(i,j,1)*GV%Z_to_H + enddo ; enddo + if (CS%tides) then + if (CS%tides_answer_date>20230630) then !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 eta(i,j) = eta(i,j) + (e_tide_eq(i,j)+e_tide_sal(i,j))*GV%Z_to_H enddo ; enddo - endif - if (CS%calculate_SAL) then + else !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - eta(i,j) = eta(i,j) + e_sal(i,j)*GV%Z_to_H + eta(i,j) = eta(i,j) + e_sal_tide(i,j)*GV%Z_to_H enddo ; enddo endif - else ! Old answers - if (CS%tides) then - !$OMP parallel do default(shared) - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - eta(i,j) = e(i,j,1)*GV%Z_to_H + (e_sal_tide(i,j))*GV%Z_to_H - enddo ; enddo - else + endif + if (CS%calculate_SAL) then + if (CS%tides_answer_date>20230630) then !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - eta(i,j) = (e(i,j,1) + e_sal(i,j))*GV%Z_to_H + eta(i,j) = eta(i,j) + e_sal(i,j)*GV%Z_to_H enddo ; enddo endif endif From 89fec6e7f4cb0dfb3bee966d14118a9235b4c241 Mon Sep 17 00:00:00 2001 From: He Wang Date: Sat, 27 Jul 2024 06:31:19 -0400 Subject: [PATCH 02/17] Add option for SAL to use bottom pressure anomaly Using SSH to calculate self-attraction and loading (SAL) is only accurate for barotropic flows. Bottom pressure anomaly should really be used for general purposes. * New runtime parameter A runtime parameter SAL_USE_BPA is added to use bottom pressure anomaly. The option requires an input file for long term mean reference bottom pressure. The reference bottom pressure field is stored with SAL_CS. * Refactor SAL and tides in Boussinesq mode As the total bottom pressure is needed for bottom pressure anomaly, SAL calculation in Boussinesq mode needs to be refactored. In addition, there is a longstanding bug in Boussinesq mode that interface height is modified by SAL and tides, and the modified interface height is erroneously used for density and pressure calculation later on. Therefore the SAL and tides are refactored by moving their calculations after pressure is calculated. Tide answers before 20230630 is retained. Answers after 20230630 are changed for Boussinesq mode. --- src/core/MOM_PressureForce_FV.F90 | 158 +++++++------ src/core/MOM_dynamics_split_RK2.F90 | 2 +- src/core/MOM_dynamics_split_RK2b.F90 | 2 +- src/core/MOM_dynamics_unsplit.F90 | 2 +- src/core/MOM_dynamics_unsplit_RK2.F90 | 2 +- .../lateral/MOM_self_attr_load.F90 | 211 ++++++++++++------ 6 files changed, 238 insertions(+), 139 deletions(-) diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90 index 905df97d68..8fec00fb90 100644 --- a/src/core/MOM_PressureForce_FV.F90 +++ b/src/core/MOM_PressureForce_FV.F90 @@ -39,8 +39,10 @@ module MOM_PressureForce_FV !> Finite volume pressure gradient control structure type, public :: PressureForce_FV_CS ; private logical :: initialized = .false. !< True if this control structure has been initialized. - logical :: calculate_SAL !< If true, calculate self-attraction and loading. - logical :: tides !< If true, apply tidal momentum forcing. + logical :: calculate_SAL = .false. !< If true, calculate self-attraction and loading. + logical :: sal_use_bpa = .false. !< If true, use bottom pressure anomaly instead of SSH + !! to calculate SAL. + logical :: tides = .false. !< If true, apply tidal momentum forcing. real :: Rho0 !< The density used in the Boussinesq !! approximation [R ~> kg m-3]. real :: GFS_scale !< A scaling of the surface pressure gradients to @@ -80,7 +82,7 @@ module MOM_PressureForce_FV !! equation of state is 0 to account for the displacement of the sea !! surface including adjustments for atmospheric or sea-ice pressure. logical :: use_stanley_pgf !< If true, turn on Stanley parameterization in the PGF - integer :: tides_answer_date !< Recover old answers with tides in Boussinesq mode + integer :: tides_answer_date = 99991231 !< Recover old answers with tides integer :: id_e_tide = -1 !< Diagnostic identifier integer :: id_e_tide_eq = -1 !< Diagnostic identifier integer :: id_e_tide_sal = -1 !< Diagnostic identifier @@ -217,6 +219,7 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ real :: dp_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [R L2 T-2 ~> Pa]. real :: I_gEarth ! The inverse of GV%g_Earth [T2 Z L-2 ~> s2 m-1] + real :: I_g_rho ! The inverse of the density times the gravitational acceleration [Z T2 L-2 R-1 ~> m Pa-1] real :: alpha_anom ! The in-situ specific volume, averaged over a ! layer, less alpha_ref [R-1 ~> m3 kg-1]. logical :: use_p_atm ! If true, use the atmospheric pressure. @@ -399,15 +402,24 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ enddo ; enddo enddo - ! Calculate and add the self-attraction and loading geopotential anomaly. + ! Calculate and add self-attraction and loading (SAL) geopotential height anomaly to interface height. if (CS%calculate_SAL) then - !$OMP parallel do default(shared) - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - SSH(i,j) = (za(i,j,1) - alpha_ref*p(i,j,1)) * I_gEarth - G%Z_ref & - - max(-G%bathyT(i,j)-G%Z_ref, 0.0) - enddo ; enddo + if (CS%sal_use_bpa) then + I_g_rho = 1.0 / (GV%rho0*GV%g_Earth) + !$OMP parallel do default(shared) + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + SSH(i,j) = p(i,j,nz+1) * I_g_rho + enddo ; enddo + else + !$OMP parallel do default(shared) + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + SSH(i,j) = (za(i,j,1) - alpha_ref*p(i,j,1)) * I_gEarth - G%Z_ref & + - max(-G%bathyT(i,j)-G%Z_ref, 0.0) + enddo ; enddo + endif call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m) + ! This gives new answers after the change of separating SAL from tidal forcing module. if ((CS%tides_answer_date>20230630) .or. (.not.GV%semi_Boussinesq) .or. (.not.CS%tides)) then !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 @@ -416,7 +428,7 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ endif endif - ! Calculate and add the tidal geopotential anomaly. + ! Calculate and add tidal geopotential height anomaly to interface height. if (CS%tides) then if ((CS%tides_answer_date>20230630) .or. (.not.GV%semi_Boussinesq)) then call calc_tidal_forcing(CS%Time, e_tide_eq, e_tide_sal, G, US, CS%tides_CSp) @@ -950,7 +962,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m]. real :: I_Rho0 ! The inverse of the Boussinesq reference density [R-1 ~> m3 kg-1]. - real :: G_Rho0 ! G_Earth / Rho0 in [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1]. + real :: G_Rho0 ! G_Earth / Rho_0 in [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1]. real :: I_g_rho ! The inverse of the density times the gravitational acceleration [Z T2 L-2 R-1 ~> m Pa-1] real :: rho_ref ! The reference density [R ~> kg m-3]. real :: dz_neglect ! A minimal thickness [Z ~> m], like e. @@ -998,6 +1010,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm G_Rho0 = GV%g_Earth / GV%Rho0 GxRho = GV%g_Earth * GV%Rho0 rho_ref = CS%Rho0 + I_g_rho = 1.0 / (CS%rho0*GV%g_Earth) ! I think it should be I_g_rho = 1.0 / (GV%rho0*GV%g_Earth) if ((CS%id_MassWt_u > 0) .or. (CS%id_MassWt_v > 0)) then MassWt_u(:,:,:) = 0.0 ; MassWt_v(:,:,:) = 0.0 @@ -1007,48 +1020,28 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm e(i,j,nz+1) = -G%bathyT(i,j) enddo ; enddo - ! Calculate and add the self-attraction and loading geopotential anomaly. - if (CS%calculate_SAL) then - ! Determine the surface height anomaly for calculating self attraction - ! and loading. This should really be based on bottom pressure anomalies, - ! but that is not yet implemented, and the current form is correct for - ! barotropic tides. + ! Self-attraction and loading (SAL) and tides (old answers) + ! SAL and tidal geopotential height anomalies are calculated and added as corrections to interface heights. + ! The following code is for recovering old answers only. The algorithm moves interface heights before density + ! calculations, and therefore is incorrect without SSH_IN_EOS_PRESSURE_FOR_PGF=True (added in August 2024). + ! See the code right after Pa calculation loop for the new method. + if (CS%tides_answer_date<=20230630) then !$OMP parallel do default(shared) do j=Jsq,Jeq+1 do i=Isq,Ieq+1 - SSH(i,j) = min(-G%bathyT(i,j) - G%Z_ref, 0.0) + SSH(i,j) = min(-G%bathyT(i,j) - G%Z_ref, 0.0) ! reference bottom pressure enddo do k=1,nz ; do i=Isq,Ieq+1 SSH(i,j) = SSH(i,j) + h(i,j,k)*GV%H_to_Z enddo ; enddo enddo call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m) - - if (CS%tides_answer_date>20230630) then - !$OMP parallel do default(shared) - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - e(i,j,nz+1) = e(i,j,nz+1) - e_sal(i,j) - enddo ; enddo - endif - endif - - ! Calculate and add the tidal geopotential anomaly. - if (CS%tides) then - if (CS%tides_answer_date>20230630) then - call calc_tidal_forcing(CS%Time, e_tide_eq, e_tide_sal, G, US, CS%tides_CSp) - !$OMP parallel do default(shared) - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - e(i,j,nz+1) = e(i,j,nz+1) - (e_tide_eq(i,j) + e_tide_sal(i,j)) - enddo ; enddo - else ! Old answers - if (.not.CS%calculate_SAL) e_sal(:,:) = 0.0 - call calc_tidal_forcing_legacy(CS%Time, e_sal, e_sal_tide, e_tide_eq, e_tide_sal, & - G, US, CS%tides_CSp) - !$OMP parallel do default(shared) - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - e(i,j,nz+1) = e(i,j,nz+1) - e_sal_tide(i,j) - enddo ; enddo - endif + call calc_tidal_forcing_legacy(CS%Time, e_sal, e_sal_tide, e_tide_eq, e_tide_sal, & + G, US, CS%tides_CSp) + !$OMP parallel do default(shared) + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + e(i,j,nz+1) = e(i,j,nz+1) - e_sal_tide(i,j) + enddo ; enddo endif !$OMP parallel do default(shared) @@ -1117,7 +1110,6 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm endif if (CS%use_SSH_in_Z0p .and. use_p_atm) then - I_g_rho = 1.0 / (CS%rho0*GV%g_Earth) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 Z_0p(i,j) = e(i,j,1) + p_atm(i,j) * I_g_rho enddo ; enddo @@ -1196,6 +1188,45 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm enddo ; enddo enddo + ! Self-attraction and loading (SAL) and tides (new answers) + if (CS%tides_answer_date>20230630) then + ! SAL + if (CS%calculate_SAL) then + if (CS%sal_use_bpa) then + !$OMP parallel do default(shared) + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + SSH(i,j) = pa(i,j,nz+1) * I_g_rho - e(i,j,nz+1) + enddo ; enddo + else + !$OMP parallel do default(shared) + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + SSH(i,j) = e(i,j,1) - max(-G%bathyT(i,j) - G%Z_ref, 0.0) ! Remove topogrpahy above sea level + enddo ; enddo + endif + call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m) + do K=1,nz+1 + !$OMP parallel do default(shared) + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + e(i,j,K) = e(i,j,K) - e_sal(i,j) + pa(i,j,K) = pa(i,j,K) - GxRho * e_sal(i,j) + enddo ; enddo + enddo + endif + + ! Tides + if (CS%tides) then + call calc_tidal_forcing(CS%Time, e_tide_eq, e_tide_sal, G, US, CS%tides_CSp) + !$OMP parallel do default(shared) + do K=1,nz+1 + !$OMP parallel do default(shared) + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + e(i,j,K) = e(i,j,K) - (e_tide_eq(i,j) + e_tide_sal(i,j)) + pa(i,j,K) = pa(i,j,K) - GxRho * (e_tide_eq(i,j) + e_tide_sal(i,j)) + enddo ; enddo + enddo + endif + endif + if (CS%correction_intxpa .or. CS%reset_intxpa_integral) then ! Determine surface temperature and salinity for use in the pressure gradient corrections if (use_ALE .and. (CS%Recon_Scheme > 0)) then @@ -1630,13 +1661,11 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm enddo ; enddo endif endif - if (CS%calculate_SAL) then - if (CS%tides_answer_date>20230630) then - !$OMP parallel do default(shared) - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - eta(i,j) = eta(i,j) + e_sal(i,j)*GV%Z_to_H - enddo ; enddo - endif + if (CS%calculate_SAL .and. (CS%tides_answer_date>20230630)) then + !$OMP parallel do default(shared) + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + eta(i,j) = eta(i,j) + e_sal(i,j)*GV%Z_to_H + enddo ; enddo endif endif @@ -1737,19 +1766,19 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp, units="kg m-3", default=GV%Rho0*US%R_to_kg_m3, scale=US%kg_m3_to_R) call get_param(param_file, mdl, "TIDES", CS%tides, & "If true, apply tidal momentum forcing.", default=.false.) - if (CS%tides) then - call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, & - "This sets the default value for the various _ANSWER_DATE parameters.", & - default=99991231) - call get_param(param_file, mdl, "TIDES_ANSWER_DATE", CS%tides_answer_date, & - "The vintage of self-attraction and loading (SAL) and tidal forcing calculations in "//& - "Boussinesq mode. Values below 20230701 recover the old answers in which the SAL is "//& - "part of the tidal forcing calculation. The change is due to a reordered summation "//& - "and the difference is only at bit level.", default=20230630) - endif + call get_param(param_file, '', "DEFAULT_ANSWER_DATE", default_answer_date, default=99991231) + if (CS%tides) & + call get_param(param_file, mdl, "TIDES_ANSWER_DATE", CS%tides_answer_date, "The vintage of "//& + "self-attraction and loading (SAL) and tidal forcing calculations. In both "//& + "Boussinesq and non-Boussinesq modes, values below 20230701 recover the old "//& + "answers in which the SAL is part of the tidal forcing calculation. The "//& + "change is due to a reordered summation and the difference is only at bit "//& + "level.", default=20230630, do_not_log=(.not.CS%tides)) call get_param(param_file, mdl, "CALCULATE_SAL", CS%calculate_SAL, & "If true, calculate self-attraction and loading.", default=CS%tides) - + if (CS%calculate_SAL) & + call get_param(param_file, '', "SAL_USE_BPA", CS%sal_use_bpa, default=.false., & + do_not_log=.true.) call get_param(param_file, "MOM", "ENABLE_THERMODYNAMICS", use_temperature, & "If true, Temperature and salinity are used as state variables.", & default=.true., do_not_log=.true.) @@ -1764,6 +1793,9 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp, "pressure used in the equation of state calculations for the Boussinesq pressure "//& "gradient forces, including adjustments for atmospheric or sea-ice pressure.", & default=.false., do_not_log=.not.GV%Boussinesq) + if (CS%tides .and. CS%tides_answer_date<=20230630 .and. CS%use_SSH_in_Z0p) & + call MOM_error(FATAL, trim(mdl) // ", PressureForce_FV_init: SSH_IN_EOS_PRESSURE_FOR_PGF "//& + "needs to be FALSE to recover tide answers before 20230630.") call get_param(param_file, "MOM", "USE_REGRIDDING", use_ALE, & "If True, use the ALE algorithm (regridding/remapping). "//& diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index f602b65240..2978bd46b9 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -1506,7 +1506,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, p 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 (CS%calculate_SAL) call SAL_init(G, US, param_file, CS%SAL_CSp) + if (CS%calculate_SAL) call SAL_init(G, GV, US, param_file, CS%SAL_CSp) if (CS%use_tides) then call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp, CS%HA_CSp) HA_CSp => CS%HA_CSp diff --git a/src/core/MOM_dynamics_split_RK2b.F90 b/src/core/MOM_dynamics_split_RK2b.F90 index 87e46795b5..31ee923aae 100644 --- a/src/core/MOM_dynamics_split_RK2b.F90 +++ b/src/core/MOM_dynamics_split_RK2b.F90 @@ -1419,7 +1419,7 @@ subroutine initialize_dyn_split_RK2b(u, v, h, tv, uh, vh, eta, Time, G, GV, US, 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 (CS%calculate_SAL) call SAL_init(G, US, param_file, CS%SAL_CSp) + if (CS%calculate_SAL) call SAL_init(G, GV, US, param_file, CS%SAL_CSp) if (CS%use_tides) then call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp, CS%HA_CSp) HA_CSp => CS%HA_CSp diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index 72c7dbe6cd..c4b7c89edc 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -707,7 +707,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 (CS%calculate_SAL) call SAL_init(G, US, param_file, CS%SAL_CSp) + if (CS%calculate_SAL) call SAL_init(G, GV, US, param_file, CS%SAL_CSp) if (CS%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%SAL_CSp, CS%tides_CSp) diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index cc37f1c2bc..7ca2c700f7 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -671,7 +671,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 (CS%calculate_SAL) call SAL_init(G, US, param_file, CS%SAL_CSp) + if (CS%calculate_SAL) call SAL_init(G, GV, US, param_file, CS%SAL_CSp) if (CS%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%SAL_CSp, CS%tides_CSp) diff --git a/src/parameterizations/lateral/MOM_self_attr_load.F90 b/src/parameterizations/lateral/MOM_self_attr_load.F90 index f72b6f513e..09adb06730 100644 --- a/src/parameterizations/lateral/MOM_self_attr_load.F90 +++ b/src/parameterizations/lateral/MOM_self_attr_load.F90 @@ -1,16 +1,18 @@ module MOM_self_attr_load -use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_MODULE -use MOM_domains, only : pass_var -use MOM_error_handler, only : MOM_error, FATAL, WARNING -use MOM_file_parser, only : read_param, get_param, log_version, param_file_type +use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_MODULE +use MOM_domains, only : pass_var +use MOM_error_handler, only : MOM_error, FATAL, WARNING +use MOM_file_parser, only : get_param, log_param, log_version, param_file_type +use MOM_grid, only : ocean_grid_type +use MOM_io, only : slasher, MOM_read_data +use MOM_load_love_numbers, only : Love_Data use MOM_obsolete_params, only : obsolete_logical, obsolete_int -use MOM_grid, only : ocean_grid_type -use MOM_unit_scaling, only : unit_scale_type use MOM_spherical_harmonics, only : spherical_harmonics_init, spherical_harmonics_end use MOM_spherical_harmonics, only : spherical_harmonics_forward, spherical_harmonics_inverse use MOM_spherical_harmonics, only : sht_CS, order2index, calc_lmax -use MOM_load_love_numbers, only : Love_Data +use MOM_unit_scaling, only : unit_scale_type +use MOM_verticalGrid, only : verticalGrid_type implicit none ; private @@ -27,61 +29,81 @@ module MOM_self_attr_load logical :: use_tidal_sal_prev = .false. !< If true, read the tidal SAL from the previous iteration of the tides to !! facilitate convergence. - real :: sal_scalar_value !< The constant of proportionality between sea surface height - !! (really it should be bottom pressure) anomalies and bottom - !! geopotential anomalies [nondim]. - type(sht_CS), allocatable :: sht !< Spherical harmonic transforms (SHT) control structure - integer :: sal_sht_Nd !< Maximum degree for SHT [nodim] - real, allocatable :: Love_Scaling(:) !< Love number for each SHT mode [nodim] - real, allocatable :: Snm_Re(:), & !< Real SHT coefficient for SHT SAL [Z ~> m] - Snm_Im(:) !< Imaginary SHT coefficient for SHT SAL [Z ~> m] + logical :: use_bpa = .false. + !< If true, use bottom pressure anomaly instead of SSH to calculate SAL. + real :: eta_prop + !< The partial derivative of eta_sal with the local value of eta [nondim]. + real :: linear_scaling + !< A dimensional coefficient for scalar approximation SAL, equal to eta_prop * unit_convert + !! [nondim or L2 Z-1 T-2 ~> m s-2 or R-1 ~> m3 kg-1 or L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1] + type(sht_CS), allocatable :: sht + !< Spherical harmonic transforms (SHT) control structure + integer :: sal_sht_Nd + !< Maximum degree for spherical harmonic transforms [nodim] + real, allocatable :: ebot_ref(:,:) + !< Reference bottom pressure normalized by Rho_0 and G_Earth[Z ~> m] + real, allocatable :: Love_scaling(:) + !< Dimensional coefficients for harmonic SAL, which are functions of Love numbers + !! [nondim or L2 Z-1 T-2 ~> m s-2 or R-1 ~> m3 kg-1 or L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1] + real, allocatable :: Snm_Re(:), Snm_Im(:) + !< Real and imaginary coefficients for harmonic SAL [Z ~> m or L2 T-2 ~> m2 s-2] end type SAL_CS integer :: id_clock_SAL !< CPU clock for self-attraction and loading contains -!> This subroutine calculates seawater self-attraction and loading based on sea surface height. This should -!! be changed into bottom pressure anomaly in the future. Note that the SAL calculation applies to all motions -!! across the spectrum. Tidal-specific methods that assume periodicity, i.e. iterative and read-in SAL, are -!! stored in MOM_tidal_forcing module. +!> This subroutine calculates seawater self-attraction and loading based on either sea surface height (SSH) or bottom +!! pressure anomaly. Note that the SAL calculation applies to all motions across the spectrum. Tidal-specific methods +!! that assume periodicity, i.e. iterative and read-in SAL, are stored in MOM_tidal_forcing module. +!! The input field is always assume to have the unit of [Z ~> m], which can be either SSH or total bottom pressure +!! normalized by mean seawater density Rho_0 and earth gravity G_Earth. For spherical harmonic method, the mean +!! seawater density would be cancelled by the same parameter in Love Number scalings. If total bottom pressure is used +!! as input, bottom pressure anomaly is calculated in the subroutine by subtracting a reference pressure from the +!! input bottom pressure. +!! The output field is expressed as geopotential height anomaly, and therefore has the unit of [Z ~> m]. subroutine calc_SAL(eta, eta_sal, G, CS, tmp_scale) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. 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_sal !< The sea surface height anomaly from - !! self-attraction and loading [Z ~> m]. + !! a time-mean geoid or total bottom pressure normalized by mean density [Z ~> m]. + real, dimension(SZI_(G),SZJ_(G)), intent(out) :: eta_sal !< The geopotential height anomaly from + !! self-attraction and loading [Z ~> m]. type(SAL_CS), intent(inout) :: CS !< The control structure returned by a previous call to SAL_init. - real, optional, intent(in) :: tmp_scale !< A rescaling factor to temporarily convert eta - !! to MKS units in reproducing sumes [m Z-1 ~> 1] + real, optional, intent(in) :: tmp_scale !< A rescaling factor to temporarily convert eta + !! to MKS units in reproducing sumes [m Z-1 ~> 1] ! Local variables integer :: n, m, l + real, dimension(SZI_(G),SZJ_(G)) :: bpa ! [Z ~> m] integer :: Isq, Ieq, Jsq, Jeq integer :: i, j - real :: eta_prop ! The scalar constant of proportionality between eta and eta_sal [nondim] call cpu_clock_begin(id_clock_SAL) Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB + if (CS%use_bpa) then ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + bpa(i,j) = eta(i,j) - CS%ebot_ref(i,j) + enddo ; enddo ; else ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + bpa(i,j) = eta(i,j) + enddo ; enddo ; endif + ! use the scalar approximation and/or iterative tidal SAL if (CS%use_sal_scalar .or. CS%use_tidal_sal_prev) then - call scalar_SAL_sensitivity(CS, eta_prop) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - eta_sal(i,j) = eta_prop*eta(i,j) + eta_sal(i,j) = CS%eta_prop * bpa(i,j) enddo ; enddo ! use the spherical harmonics method elseif (CS%use_sal_sht) then - call spherical_harmonics_forward(G, CS%sht, eta, CS%Snm_Re, CS%Snm_Im, CS%sal_sht_Nd, tmp_scale=tmp_scale) + call spherical_harmonics_forward(G, CS%sht, bpa, CS%Snm_Re, CS%Snm_Im, CS%sal_sht_Nd, tmp_scale=tmp_scale) ! Multiply scaling factors to each mode do m = 0,CS%sal_sht_Nd l = order2index(m, CS%sal_sht_Nd) do n = m,CS%sal_sht_Nd - CS%Snm_Re(l+n-m) = CS%Snm_Re(l+n-m) * CS%Love_Scaling(l+n-m) - CS%Snm_Im(l+n-m) = CS%Snm_Im(l+n-m) * CS%Love_Scaling(l+n-m) + CS%Snm_Re(l+n-m) = CS%Snm_Re(l+n-m) * CS%Love_scaling(l+n-m) + CS%Snm_Im(l+n-m) = CS%Snm_Im(l+n-m) * CS%Love_scaling(l+n-m) enddo enddo @@ -98,36 +120,32 @@ subroutine calc_SAL(eta, eta_sal, G, CS, tmp_scale) call cpu_clock_end(id_clock_SAL) end subroutine calc_SAL -!> This subroutine calculates the partial derivative of the local geopotential height with the input -!! sea surface height due to the scalar approximation of self-attraction and loading. +!> This subroutine returns eta_prop member of SAL_CS type, which is the non-dimensional partial +!! derivative of the local geopotential height with the input sea surface height due to the scalar +!! approximation of self-attraction and loading. subroutine scalar_SAL_sensitivity(CS, deta_sal_deta) type(SAL_CS), intent(in) :: CS !< The control structure returned by a previous call to SAL_init. real, intent(out) :: deta_sal_deta !< The partial derivative of eta_sal with !! the local value of eta [nondim]. - - if (CS%use_sal_scalar .and. CS%use_tidal_sal_prev) then - deta_sal_deta = 2.0*CS%sal_scalar_value - elseif (CS%use_sal_scalar .or. CS%use_tidal_sal_prev) then - deta_sal_deta = CS%sal_scalar_value - else - deta_sal_deta = 0.0 - endif + deta_sal_deta = CS%eta_prop end subroutine scalar_SAL_sensitivity !> This subroutine calculates coefficients of the spherical harmonic modes for self-attraction and loading. !! The algorithm is based on the SAL implementation in MPAS-ocean, which was modified by Kristin Barton from !! routine written by K. Quinn (March 2010) and modified by M. Schindelegger (May 2017). -subroutine calc_love_scaling(nlm, rhoW, rhoE, Love_Scaling) +subroutine calc_love_scaling(nlm, rhoW, rhoE, Love_scaling) integer, intent(in) :: nlm !< Maximum spherical harmonics degree [nondim] real, intent(in) :: rhoW !< The average density of sea water [R ~> kg m-3] real, intent(in) :: rhoE !< The average density of Earth [R ~> kg m-3] - real, dimension(:), intent(out) :: Love_Scaling !< Scaling factors for inverse SHT [nondim] + real, dimension(:), intent(out) :: Love_scaling !< Scaling factors for inverse spherical harmonic + !! transforms [nondim] ! Local variables real, dimension(:), allocatable :: HDat, LDat, KDat ! Love numbers converted in CF reference frames [nondim] real :: H1, L1, K1 ! Temporary variables to store degree 1 Love numbers [nondim] integer :: n_tot ! Size of the stored Love numbers integer :: n, m, l + real :: unit n_tot = size(Love_Data, dim=2) @@ -148,34 +166,40 @@ subroutine calc_love_scaling(nlm, rhoW, rhoE, Love_Scaling) do m=0,nlm ; do n=m,nlm l = order2index(m,nlm) - Love_Scaling(l+n-m) = (3.0 / real(2*n+1)) * (rhoW / rhoE) * (1.0 + KDat(n+1) - HDat(n+1)) + Love_scaling(l+n-m) = (3.0 / real(2*n+1)) * (rhoW / rhoE) * (1.0 + KDat(n+1) - HDat(n+1)) enddo ; enddo end subroutine calc_love_scaling !> This subroutine initializes the self-attraction and loading control structure. -subroutine SAL_init(G, US, param_file, CS) - 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. +subroutine SAL_init(G, GV, US, param_file, CS) + type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< Vertical 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(SAL_CS), intent(inout) :: CS !< Self-attraction and loading control structure ! Local variables # include "version_variable.h" character(len=40) :: mdl = "MOM_self_attr_load" ! This module's name. integer :: lmax ! Total modes of the real spherical harmonics [nondim] - real :: rhoW ! The average density of sea water [R ~> kg m-3]. real :: rhoE ! The average density of Earth [R ~> kg m-3]. - - logical :: calculate_sal - logical :: tides, use_tidal_sal_file - real :: tide_sal_scalar_value ! Scaling SAL factor [nondim] + character(len=200) :: filename, ebot_ref_file, inputdir ! Strings for file/path + character(len=200) :: ebot_ref_varname ! Variable name in file + real :: I_g_rho ! The inverse of the density times the gravitational acceleration [Z T2 L-2 R-1 ~> m Pa-1] + logical :: calculate_sal=.false. + logical :: tides=.false., use_tidal_sal_file=.false., bq_sal_tides_bug=.false. + integer :: tides_answer_date=99991203 ! Recover old answers with tides + real :: sal_scalar_value, tide_sal_scalar_value ! Scaling SAL factors [nondim] + + integer :: isd, ied, jsd, jed + integer :: i, j + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ! Read all relevant parameters and write them to the model log. call log_version(param_file, mdl, version, "") call get_param(param_file, '', "TIDES", tides, default=.false., do_not_log=.True.) - call get_param(param_file, mdl, "CALCULATE_SAL", calculate_sal, "If true, calculate "//& - " self-attraction and loading.", default=tides, do_not_log=.True.) + call get_param(param_file, '', "CALCULATE_SAL", calculate_sal, default=tides, do_not_log=.True.) if (.not. calculate_sal) return if (tides) then @@ -183,17 +207,51 @@ subroutine SAL_init(G, US, param_file, CS) default=.false., do_not_log=.True.) call get_param(param_file, '', "TIDAL_SAL_FROM_FILE", use_tidal_sal_file, & default=.false., do_not_log=.True.) + call get_param(param_file, '', "TIDES_ANSWER_DATE", tides_answer_date, & + default=20230630, do_not_log=.True.) endif + call get_param(param_file, mdl, "SAL_USE_BPA", CS%use_bpa, & + "If true, use bottom pressure anomaly to calculate self-attraction and "// & + "loading (SAL). Otherwise sea surface height anomaly is used, which is "// & + "only correct for homogenous flow.", default=.False.) + if (CS%use_bpa) then + call get_param(param_file, '', "INPUTDIR", inputdir, default=".", do_not_log=.True.) + inputdir = slasher(inputdir) + call get_param(param_file, mdl, "REF_BOT_PRES_FILE", ebot_ref_file, & + "Reference bottom pressure file used by self-attraction and loading (SAL).", & + default="pbot.nc") + call get_param(param_file, mdl, "REF_BOT_PRES_VARNAME", ebot_ref_varname, & + "The name of the variable in REF_BOT_PRES_FILE with reference bottom "//& + "pressure. The variable should have the unit of Pa.", & + default="pbot") + filename = trim(inputdir)//trim(ebot_ref_file) + call log_param(param_file, mdl, "INPUTDIR/REF_BOT_PRES_FILE", filename) + + allocate(CS%ebot_ref(isd:ied, jsd:jed), source=0.0) + call MOM_read_data(filename, trim(ebot_ref_varname), CS%ebot_ref, G%Domain,& + scale=US%Pa_to_RL2_T2) + I_g_rho = 1.0 / (GV%g_Earth * GV%Rho0) + do j=jsd,jed ; do i=isd,ied + CS%ebot_ref(i,j) = CS%ebot_ref(i,j) * I_g_rho + enddo ; enddo + call pass_var(CS%ebot_ref, G%Domain) + endif + if (tides_answer_date<=20230630 .and. CS%use_bpa) & + call MOM_error(FATAL, trim(mdl) // ", SAL_init: SAL_USE_BPA needs to be false to recover "//& + "tide answers before 20230630.") call get_param(param_file, mdl, "SAL_SCALAR_APPROX", CS%use_sal_scalar, & "If true, use the scalar approximation to calculate self-attraction and "//& "loading.", default=tides .and. (.not. use_tidal_sal_file)) + ! if (CS%use_sal_scalar .and. CS%use_bpa) & + ! call MOM_error(WARNING, trim(mdl) // ", SAL_init: Using bottom pressure anomaly for scalar "//& + ! "approximation SAL is unsubstantiated.") call get_param(param_file, '', "TIDE_SAL_SCALAR_VALUE", tide_sal_scalar_value, & units="m m-1", default=0.0, do_not_log=.True.) if (tide_sal_scalar_value/=0.0) & call MOM_error(WARNING, "TIDE_SAL_SCALAR_VALUE is a deprecated parameter. "//& "Use SAL_SCALAR_VALUE instead." ) - call get_param(param_file, mdl, "SAL_SCALAR_VALUE", CS%sal_scalar_value, & + call get_param(param_file, mdl, "SAL_SCALAR_VALUE", sal_scalar_value, & "The constant of proportionality between sea surface "//& "height (really it should be bottom pressure) anomalies "//& "and bottom geopotential anomalies. This is only used if "//& @@ -207,20 +265,28 @@ subroutine SAL_init(G, US, param_file, CS) "The maximum degree of the spherical harmonics transformation used for "// & "calculating the self-attraction and loading term.", & default=0, do_not_log=(.not. CS%use_sal_sht)) - call get_param(param_file, '', "RHO_0", rhoW, default=1035.0, scale=US%kg_m3_to_R, & - units="kg m-3", do_not_log=.True.) call get_param(param_file, mdl, "RHO_SOLID_EARTH", rhoE, & "The mean solid earth density. This is used for calculating the "// & "self-attraction and loading term.", units="kg m-3", & default=5517.0, scale=US%kg_m3_to_R, do_not_log=(.not. CS%use_sal_sht)) + ! Set scaling coefficients for scalar approximation + if (CS%use_sal_scalar .and. CS%use_tidal_sal_prev) then + CS%eta_prop = 2.0 * sal_scalar_value + elseif (CS%use_sal_scalar .or. CS%use_tidal_sal_prev) then + CS%eta_prop = sal_scalar_value + else + CS%eta_prop = 0.0 + endif + + ! Set scaling coefficients for spherical harmonics if (CS%use_sal_sht) then lmax = calc_lmax(CS%sal_sht_Nd) allocate(CS%Snm_Re(lmax)); CS%Snm_Re(:) = 0.0 allocate(CS%Snm_Im(lmax)); CS%Snm_Im(:) = 0.0 - allocate(CS%Love_Scaling(lmax)); CS%Love_Scaling(:) = 0.0 - call calc_love_scaling(CS%sal_sht_Nd, rhoW, rhoE, CS%Love_Scaling) + allocate(CS%Love_scaling(lmax)); CS%Love_scaling(:) = 0.0 + call calc_love_scaling(CS%sal_sht_Nd, GV%Rho0, rhoE, CS%Love_scaling) allocate(CS%sht) call spherical_harmonics_init(G, param_file, CS%sht) @@ -234,8 +300,11 @@ end subroutine SAL_init subroutine SAL_end(CS) type(SAL_CS), intent(inout) :: CS !< The control structure returned by a previous call !! to SAL_init; it is deallocated here. + + if (allocated(CS%ebot_ref)) deallocate(CS%ebot_ref) + if (CS%use_sal_sht) then - if (allocated(CS%Love_Scaling)) deallocate(CS%Love_Scaling) + if (allocated(CS%Love_scaling)) deallocate(CS%Love_scaling) if (allocated(CS%Snm_Re)) deallocate(CS%Snm_Re) if (allocated(CS%Snm_Im)) deallocate(CS%Snm_Im) call spherical_harmonics_end(CS%sht) @@ -245,22 +314,20 @@ end subroutine SAL_end !> \namespace self_attr_load !! -!! \section section_SAL Self attraction and loading -!! -!! This module contains methods to calculate self-attraction and loading (SAL) as a function of sea surface height (SSH) -!! (rather, it should be bottom pressure anomaly). SAL is primarily used for fast evolving processes like tides or -!! storm surges, but the effect applies to all motions. +!! This module contains methods to calculate self-attraction and loading (SAL) as a function of sea surface height or +!! bottom pressure anomaly. SAL is primarily used for fast evolving processes like tides or storm surges, but the +!! effect applies to all motions. !! !! If SAL_SCALAR_APPROX is true, a scalar approximation is applied (\cite Accad1978) and the SAL is simply -!! a fraction (set by SAL_SCALAR_VALUE, usually around 10% for global tides) of local SSH. -!! For tides, the scalar approximation can also be used to iterate the SAL to convergence [see -!! USE_PREVIOUS_TIDES in MOM_tidal_forcing, \cite Arbic2004]. +!! a fraction (set by SAL_SCALAR_VALUE, usually around 10% for global tides) of local SSH. For tides, the +!! scalar approximation can also be used to iterate the SAL to convergence [see USE_PREVIOUS_TIDES in +!! MOM_tidal_forcing, \cite Arbic2004]. !! !! If SAL_HARMONICS is true, a more accurate online spherical harmonic transforms are used to calculate -!! SAL. Subroutines in module MOM_spherical_harmonics are called and the degree of spherical harmonic transforms is -!! set by SAL_HARMONICS_DEGREE. The algorithm is based on SAL calculation in Model for Prediction Across -!! Scales (MPAS)-Ocean -!! developed by Los Alamos National Laboratory and University of Michigan [\cite Barton2022 and \cite Brus2023]. +!! SAL. Subroutines in module MOM_spherical_harmonics are called and the degree of spherical harmonic transforms is set +!! by SAL_HARMONICS_DEGREE. The algorithm is based on SAL calculation in Model for Prediction Across +!! Scales (MPAS)-Ocean developed by Los Alamos National Laboratory and University of Michigan +!! [\cite Barton2022 and \cite Brus2023]. !! !! References: !! From 402949e6fb95b7eb87206d28dd820bccf12d15e8 Mon Sep 17 00:00:00 2001 From: He Wang Date: Fri, 27 Sep 2024 11:03:27 -0400 Subject: [PATCH 03/17] Rename variables with tides in Pressure Force * Local variable SSH is renamed to pbot_anom. * Tides related local variable names are changed to be less confusing. Notably, `e_sal_tide` is renamed to `e_sal_and_tide` (the summation of sal and tide), not to be confused with `e_tide_sal`, which is renamed to `e_tidal_sal` (sal from tides). * Fix e_tidal diagnostics in Boussinesq mode. The term was not assigned if tide answer date is >20230630. --- src/core/MOM_PressureForce_FV.F90 | 79 ++++++++++--------- .../lateral/MOM_self_attr_load.F90 | 15 ++-- 2 files changed, 47 insertions(+), 47 deletions(-) diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90 index 8fec00fb90..f86ce5e57f 100644 --- a/src/core/MOM_PressureForce_FV.F90 +++ b/src/core/MOM_PressureForce_FV.F90 @@ -84,8 +84,8 @@ module MOM_PressureForce_FV logical :: use_stanley_pgf !< If true, turn on Stanley parameterization in the PGF integer :: tides_answer_date = 99991231 !< Recover old answers with tides integer :: id_e_tide = -1 !< Diagnostic identifier - integer :: id_e_tide_eq = -1 !< Diagnostic identifier - integer :: id_e_tide_sal = -1 !< Diagnostic identifier + integer :: id_e_tidal_eq = -1 !< Diagnostic identifier + integer :: id_e_tidal_sal = -1 !< Diagnostic identifier integer :: id_e_sal = -1 !< Diagnostic identifier integer :: id_rho_pgf = -1 !< Diagnostic identifier integer :: id_rho_stanley_pgf = -1 !< Diagnostic identifier @@ -143,12 +143,13 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ ! the pressure anomaly at the top of the layer [R L4 T-4 ~> Pa m2 s-2]. real, dimension(SZI_(G),SZJ_(G)) :: & dp, & ! The (positive) change in pressure across a layer [R L2 T-2 ~> Pa]. - SSH, & ! The sea surface height anomaly, in depth units [Z ~> m]. + pbot_anom, & ! Bottom pressure (anomaly) in depth units, used for self-attraction and loading [Z ~> m]. e_sal, & ! The bottom geopotential anomaly due to self-attraction and loading [Z ~> m]. - e_tide_eq, & ! The bottom geopotential anomaly due to tidal forces from astronomical sources [Z ~> m]. - e_tide_sal, & ! The bottom geopotential anomaly due to harmonic self-attraction and loading + e_tidal_eq, & ! The bottom geopotential anomaly due to tidal forces from astronomical sources [Z ~> m]. + e_tidal_sal, & ! The bottom geopotential anomaly due to harmonic self-attraction and loading ! specific to tides [Z ~> m]. - e_sal_tide, & ! The summation of self-attraction and loading and tidal forcing [Z ~> m]. + e_sal_and_tide, & ! The summation of self-attraction and loading and tidal forcing, used for recovering + ! old answers only [Z ~> m]. dM ! The barotropic adjustment to the Montgomery potential to ! account for a reduced gravity model [L2 T-2 ~> m2 s-2]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: & @@ -408,16 +409,16 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ I_g_rho = 1.0 / (GV%rho0*GV%g_Earth) !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - SSH(i,j) = p(i,j,nz+1) * I_g_rho + pbot_anom(i,j) = p(i,j,nz+1) * I_g_rho enddo ; enddo else !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - SSH(i,j) = (za(i,j,1) - alpha_ref*p(i,j,1)) * I_gEarth - G%Z_ref & + pbot_anom(i,j) = (za(i,j,1) - alpha_ref*p(i,j,1)) * I_gEarth - G%Z_ref & - max(-G%bathyT(i,j)-G%Z_ref, 0.0) enddo ; enddo endif - call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m) + call calc_SAL(pbot_anom, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m) ! This gives new answers after the change of separating SAL from tidal forcing module. if ((CS%tides_answer_date>20230630) .or. (.not.GV%semi_Boussinesq) .or. (.not.CS%tides)) then @@ -431,18 +432,18 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ ! Calculate and add tidal geopotential height anomaly to interface height. if (CS%tides) then if ((CS%tides_answer_date>20230630) .or. (.not.GV%semi_Boussinesq)) then - call calc_tidal_forcing(CS%Time, e_tide_eq, e_tide_sal, G, US, CS%tides_CSp) + call calc_tidal_forcing(CS%Time, e_tidal_eq, e_tidal_sal, G, US, CS%tides_CSp) !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - za(i,j,1) = za(i,j,1) - GV%g_Earth * (e_tide_eq(i,j) + e_tide_sal(i,j)) + za(i,j,1) = za(i,j,1) - GV%g_Earth * (e_tidal_eq(i,j) + e_tidal_sal(i,j)) enddo ; enddo else ! This block recreates older answers with tides. if (.not.CS%calculate_SAL) e_sal(:,:) = 0.0 - call calc_tidal_forcing_legacy(CS%Time, e_sal, e_sal_tide, e_tide_eq, e_tide_sal, & + call calc_tidal_forcing_legacy(CS%Time, e_sal, e_sal_and_tide, e_tidal_eq, e_tidal_sal, & G, US, CS%tides_CSp) !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - za(i,j,1) = za(i,j,1) - GV%g_Earth * e_sal_tide(i,j) + za(i,j,1) = za(i,j,1) - GV%g_Earth * e_sal_and_tide(i,j) enddo ; enddo endif endif @@ -820,10 +821,10 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ ! To be consistent with old runs, tidal forcing diagnostic also includes total SAL. ! New diagnostics are given for each individual field. - if (CS%id_e_tide>0) call post_data(CS%id_e_tide, e_sal+e_tide_eq+e_tide_sal, CS%diag) + if (CS%id_e_tide>0) call post_data(CS%id_e_tide, e_sal+e_tidal_eq+e_tidal_sal, CS%diag) if (CS%id_e_sal>0) call post_data(CS%id_e_sal, e_sal, CS%diag) - if (CS%id_e_tide_eq>0) call post_data(CS%id_e_tide_eq, e_tide_eq, CS%diag) - if (CS%id_e_tide_sal>0) call post_data(CS%id_e_tide_sal, e_tide_sal, CS%diag) + if (CS%id_e_tidal_eq>0) call post_data(CS%id_e_tidal_eq, e_tidal_eq, CS%diag) + if (CS%id_e_tidal_sal>0) call post_data(CS%id_e_tidal_sal, e_tidal_sal, CS%diag) if (CS%id_MassWt_u>0) call post_data(CS%id_MassWt_u, MassWt_u, CS%diag) if (CS%id_MassWt_v>0) call post_data(CS%id_MassWt_v, MassWt_v, CS%diag) @@ -858,14 +859,14 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm ! Local variables real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: e ! Interface height in depth units [Z ~> m]. real, dimension(SZI_(G),SZJ_(G)) :: & - e_sal_tide, & ! The summation of self-attraction and loading and tidal forcing [Z ~> m]. + e_sal_and_tide, & ! The summation of self-attraction and loading and tidal forcing [Z ~> m]. e_sal, & ! The bottom geopotential anomaly due to self-attraction and loading [Z ~> m]. - e_tide_eq, & ! The bottom geopotential anomaly due to tidal forces from astronomical sources + e_tidal_eq, & ! The bottom geopotential anomaly due to tidal forces from astronomical sources ! [Z ~> m]. - e_tide_sal, & ! The bottom geopotential anomaly due to harmonic self-attraction and loading + e_tidal_sal, & ! The bottom geopotential anomaly due to harmonic self-attraction and loading ! specific to tides [Z ~> m]. Z_0p, & ! The height at which the pressure used in the equation of state is 0 [Z ~> m] - SSH, & ! The sea surface height anomaly, in depth units [Z ~> m]. + pbot_anom, & ! Bottom pressure (anomaly) in depth units, used for self-attraction and loading [Z ~> m]. dM ! The barotropic adjustment to the Montgomery potential to ! account for a reduced gravity model [L2 T-2 ~> m2 s-2]. real, dimension(SZI_(G)) :: & @@ -1029,18 +1030,18 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm !$OMP parallel do default(shared) do j=Jsq,Jeq+1 do i=Isq,Ieq+1 - SSH(i,j) = min(-G%bathyT(i,j) - G%Z_ref, 0.0) ! reference bottom pressure + pbot_anom(i,j) = min(-G%bathyT(i,j) - G%Z_ref, 0.0) ! reference bottom pressure enddo do k=1,nz ; do i=Isq,Ieq+1 - SSH(i,j) = SSH(i,j) + h(i,j,k)*GV%H_to_Z + pbot_anom(i,j) = pbot_anom(i,j) + h(i,j,k)*GV%H_to_Z enddo ; enddo enddo - call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m) - call calc_tidal_forcing_legacy(CS%Time, e_sal, e_sal_tide, e_tide_eq, e_tide_sal, & + call calc_SAL(pbot_anom, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m) + call calc_tidal_forcing_legacy(CS%Time, e_sal, e_sal_and_tide, e_tidal_eq, e_tidal_sal, & G, US, CS%tides_CSp) !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - e(i,j,nz+1) = e(i,j,nz+1) - e_sal_tide(i,j) + e(i,j,nz+1) = e(i,j,nz+1) - e_sal_and_tide(i,j) enddo ; enddo endif @@ -1195,15 +1196,15 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm if (CS%sal_use_bpa) then !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - SSH(i,j) = pa(i,j,nz+1) * I_g_rho - e(i,j,nz+1) + pbot_anom(i,j) = pa(i,j,nz+1) * I_g_rho - e(i,j,nz+1) enddo ; enddo else !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - SSH(i,j) = e(i,j,1) - max(-G%bathyT(i,j) - G%Z_ref, 0.0) ! Remove topogrpahy above sea level + pbot_anom(i,j) = e(i,j,1) - max(-G%bathyT(i,j) - G%Z_ref, 0.0) ! Remove topogrpahy above sea level enddo ; enddo endif - call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m) + call calc_SAL(pbot_anom, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m) do K=1,nz+1 !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 @@ -1215,13 +1216,13 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm ! Tides if (CS%tides) then - call calc_tidal_forcing(CS%Time, e_tide_eq, e_tide_sal, G, US, CS%tides_CSp) + call calc_tidal_forcing(CS%Time, e_tidal_eq, e_tidal_sal, G, US, CS%tides_CSp) !$OMP parallel do default(shared) do K=1,nz+1 !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - e(i,j,K) = e(i,j,K) - (e_tide_eq(i,j) + e_tide_sal(i,j)) - pa(i,j,K) = pa(i,j,K) - GxRho * (e_tide_eq(i,j) + e_tide_sal(i,j)) + e(i,j,K) = e(i,j,K) - (e_tidal_eq(i,j) + e_tidal_sal(i,j)) + pa(i,j,K) = pa(i,j,K) - GxRho * (e_tidal_eq(i,j) + e_tidal_sal(i,j)) enddo ; enddo enddo endif @@ -1652,12 +1653,12 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm if (CS%tides_answer_date>20230630) then !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - eta(i,j) = eta(i,j) + (e_tide_eq(i,j)+e_tide_sal(i,j))*GV%Z_to_H + eta(i,j) = eta(i,j) + (e_tidal_eq(i,j)+e_tidal_sal(i,j))*GV%Z_to_H enddo ; enddo else !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - eta(i,j) = eta(i,j) + e_sal_tide(i,j)*GV%Z_to_H + eta(i,j) = eta(i,j) + e_sal_and_tide(i,j)*GV%Z_to_H enddo ; enddo endif endif @@ -1708,10 +1709,10 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm ! To be consistent with old runs, tidal forcing diagnostic also includes total SAL. ! New diagnostics are given for each individual field. - if (CS%id_e_tide>0) call post_data(CS%id_e_tide, e_sal_tide, CS%diag) + if (CS%id_e_tide>0) call post_data(CS%id_e_tide, e_sal+e_tidal_eq+e_tidal_sal, CS%diag) if (CS%id_e_sal>0) call post_data(CS%id_e_sal, e_sal, CS%diag) - if (CS%id_e_tide_eq>0) call post_data(CS%id_e_tide_eq, e_tide_eq, CS%diag) - if (CS%id_e_tide_sal>0) call post_data(CS%id_e_tide_sal, e_tide_sal, CS%diag) + if (CS%id_e_tidal_eq>0) call post_data(CS%id_e_tidal_eq, e_tidal_eq, CS%diag) + if (CS%id_e_tidal_sal>0) call post_data(CS%id_e_tidal_sal, e_tidal_sal, CS%diag) if (CS%id_MassWt_u>0) call post_data(CS%id_MassWt_u, MassWt_u, CS%diag) if (CS%id_MassWt_v>0) call post_data(CS%id_MassWt_v, MassWt_v, CS%diag) @@ -1795,7 +1796,7 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp, default=.false., do_not_log=.not.GV%Boussinesq) if (CS%tides .and. CS%tides_answer_date<=20230630 .and. CS%use_SSH_in_Z0p) & call MOM_error(FATAL, trim(mdl) // ", PressureForce_FV_init: SSH_IN_EOS_PRESSURE_FOR_PGF "//& - "needs to be FALSE to recover tide answers before 20230630.") + "needs to be FALSE to recover tide answers before 20230701.") call get_param(param_file, "MOM", "USE_REGRIDDING", use_ALE, & "If True, use the ALE algorithm (regridding/remapping). "//& @@ -1880,9 +1881,9 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp, if (CS%tides) then CS%id_e_tide = register_diag_field('ocean_model', 'e_tidal', diag%axesT1, Time, & 'Tidal Forcing Astronomical and SAL Height Anomaly', 'meter', conversion=US%Z_to_m) - CS%id_e_tide_eq = register_diag_field('ocean_model', 'e_tide_eq', diag%axesT1, Time, & + CS%id_e_tidal_eq = register_diag_field('ocean_model', 'e_tide_eq', diag%axesT1, Time, & 'Equilibrium tides height anomaly', 'meter', conversion=US%Z_to_m) - CS%id_e_tide_sal = register_diag_field('ocean_model', 'e_tide_sal', diag%axesT1, Time, & + CS%id_e_tidal_sal = register_diag_field('ocean_model', 'e_tide_sal', diag%axesT1, Time, & 'Read-in tidal self-attraction and loading height anomaly', 'meter', conversion=US%Z_to_m) endif diff --git a/src/parameterizations/lateral/MOM_self_attr_load.F90 b/src/parameterizations/lateral/MOM_self_attr_load.F90 index 09adb06730..fe4149d602 100644 --- a/src/parameterizations/lateral/MOM_self_attr_load.F90 +++ b/src/parameterizations/lateral/MOM_self_attr_load.F90 @@ -41,7 +41,7 @@ module MOM_self_attr_load integer :: sal_sht_Nd !< Maximum degree for spherical harmonic transforms [nodim] real, allocatable :: ebot_ref(:,:) - !< Reference bottom pressure normalized by Rho_0 and G_Earth[Z ~> m] + !< Reference bottom pressure scaled by Rho_0 and G_Earth[Z ~> m] real, allocatable :: Love_scaling(:) !< Dimensional coefficients for harmonic SAL, which are functions of Love numbers !! [nondim or L2 Z-1 T-2 ~> m s-2 or R-1 ~> m3 kg-1 or L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1] @@ -56,16 +56,15 @@ module MOM_self_attr_load !> This subroutine calculates seawater self-attraction and loading based on either sea surface height (SSH) or bottom !! pressure anomaly. Note that the SAL calculation applies to all motions across the spectrum. Tidal-specific methods !! that assume periodicity, i.e. iterative and read-in SAL, are stored in MOM_tidal_forcing module. -!! The input field is always assume to have the unit of [Z ~> m], which can be either SSH or total bottom pressure -!! normalized by mean seawater density Rho_0 and earth gravity G_Earth. For spherical harmonic method, the mean -!! seawater density would be cancelled by the same parameter in Love Number scalings. If total bottom pressure is used -!! as input, bottom pressure anomaly is calculated in the subroutine by subtracting a reference pressure from the -!! input bottom pressure. -!! The output field is expressed as geopotential height anomaly, and therefore has the unit of [Z ~> m]. +!! The input field is always assume to have the unit of [Z ~> m], which can be either SSH or total bottom pressure +!! scaled by mean seawater density Rho_0 and earth gravity G_Earth. For spherical harmonics, the mean seawater density +!! would be cancelled by the same parameter in Love number scalings. If total bottom pressure is used as input, bottom +!! pressure anomaly is calculated in the subroutine by subtracting a reference pressure from the input bottom pressure. +!! The output field is expressed as geopotential height anomaly, and therefore has the unit of [Z ~> m]. subroutine calc_SAL(eta, eta_sal, G, CS, tmp_scale) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. real, dimension(SZI_(G),SZJ_(G)), intent(in) :: eta !< The sea surface height anomaly from - !! a time-mean geoid or total bottom pressure normalized by mean density [Z ~> m]. + !! a time-mean geoid or total bottom pressure scaled by mean density and earth gravity [Z ~> m]. real, dimension(SZI_(G),SZJ_(G)), intent(out) :: eta_sal !< The geopotential height anomaly from !! self-attraction and loading [Z ~> m]. type(SAL_CS), intent(inout) :: CS !< The control structure returned by a previous call to SAL_init. From 8039c33655831238200a62e6bef33e8dc2affb82 Mon Sep 17 00:00:00 2001 From: He Wang Date: Fri, 27 Sep 2024 11:55:51 -0400 Subject: [PATCH 04/17] Alternative method for SAL and tides in Boussinesq Add an alternative method for SAL and tides in Boussinesq mode. The current method adjusts interface heights with geopotential height anomaly for SAL and tides. For non-Boussinesq mode, the current method is algebraically the same as taking the gradient of SAL and tide geopotential (body forcing). For Boussinesq mode, there is a baroclinic component of tidal forcing and SAL. The alternative method is added to calculate the gradient of tidal forcing and SAL directly at the cost of additional multiplications. The new method is controlled by runtime parameter BOUSSINESQ_SAL_TIDES. --- src/core/MOM_PressureForce_FV.F90 | 120 ++++++++++++------ .../lateral/MOM_self_attr_load.F90 | 4 +- 2 files changed, 84 insertions(+), 40 deletions(-) diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90 index f86ce5e57f..7a96de15ac 100644 --- a/src/core/MOM_PressureForce_FV.F90 +++ b/src/core/MOM_PressureForce_FV.F90 @@ -82,6 +82,8 @@ module MOM_PressureForce_FV !! equation of state is 0 to account for the displacement of the sea !! surface including adjustments for atmospheric or sea-ice pressure. logical :: use_stanley_pgf !< If true, turn on Stanley parameterization in the PGF + logical :: bq_sal_tides = .false. !< If true, use an alternative method for SAL and tides + !! in Boussinesq mode integer :: tides_answer_date = 99991231 !< Recover old answers with tides integer :: id_e_tide = -1 !< Diagnostic identifier integer :: id_e_tidal_eq = -1 !< Diagnostic identifier @@ -821,7 +823,12 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ ! To be consistent with old runs, tidal forcing diagnostic also includes total SAL. ! New diagnostics are given for each individual field. - if (CS%id_e_tide>0) call post_data(CS%id_e_tide, e_sal+e_tidal_eq+e_tidal_sal, CS%diag) + if (CS%id_e_tide>0) then + if (CS%tides_answer_date>20230630) then ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + e_sal_and_tide(i,j) = e_sal(i,j) + e_tidal_eq(i,j) + e_tidal_sal(i,j) + enddo ; enddo ; endif + call post_data(CS%id_e_tide, e_sal_and_tide, CS%diag) + endif if (CS%id_e_sal>0) call post_data(CS%id_e_sal, e_sal, CS%diag) if (CS%id_e_tidal_eq>0) call post_data(CS%id_e_tidal_eq, e_tidal_eq, CS%diag) if (CS%id_e_tidal_sal>0) call post_data(CS%id_e_tidal_sal, e_tidal_sal, CS%diag) @@ -1189,43 +1196,39 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm enddo ; enddo enddo - ! Self-attraction and loading (SAL) and tides (new answers) - if (CS%tides_answer_date>20230630) then - ! SAL - if (CS%calculate_SAL) then - if (CS%sal_use_bpa) then - !$OMP parallel do default(shared) - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - pbot_anom(i,j) = pa(i,j,nz+1) * I_g_rho - e(i,j,nz+1) - enddo ; enddo - else - !$OMP parallel do default(shared) - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - pbot_anom(i,j) = e(i,j,1) - max(-G%bathyT(i,j) - G%Z_ref, 0.0) ! Remove topogrpahy above sea level - enddo ; enddo - endif - call calc_SAL(pbot_anom, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m) - do K=1,nz+1 - !$OMP parallel do default(shared) - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - e(i,j,K) = e(i,j,K) - e_sal(i,j) - pa(i,j,K) = pa(i,j,K) - GxRho * e_sal(i,j) - enddo ; enddo - enddo + ! Self-attraction and loading (SAL) (new answers) + if (CS%calculate_SAL .and. CS%tides_answer_date>20230630) then + if (CS%sal_use_bpa) then + !$OMP parallel do default(shared) + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + pbot_anom(i,j) = pa(i,j,nz+1) * I_g_rho - e(i,j,nz+1) + enddo ; enddo + else + !$OMP parallel do default(shared) + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + pbot_anom(i,j) = e(i,j,1) - max(-G%bathyT(i,j) - G%Z_ref, 0.0) ! Remove topography above sea level + enddo ; enddo endif + call calc_SAL(pbot_anom, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m) + if (.not.CS%bq_sal_tides) then ; do K=1,nz+1 + !$OMP parallel do default(shared) + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + e(i,j,K) = e(i,j,K) - e_sal(i,j) + pa(i,j,K) = pa(i,j,K) - GxRho * e_sal(i,j) + enddo ; enddo + enddo ; endif + endif - ! Tides - if (CS%tides) then - call calc_tidal_forcing(CS%Time, e_tidal_eq, e_tidal_sal, G, US, CS%tides_CSp) + ! Tides (new answers) + if (CS%tides .and. CS%tides_answer_date>20230630) then + call calc_tidal_forcing(CS%Time, e_tidal_eq, e_tidal_sal, G, US, CS%tides_CSp) + if (.not.CS%bq_sal_tides) then ; do K=1,nz+1 !$OMP parallel do default(shared) - do K=1,nz+1 - !$OMP parallel do default(shared) - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - e(i,j,K) = e(i,j,K) - (e_tidal_eq(i,j) + e_tidal_sal(i,j)) - pa(i,j,K) = pa(i,j,K) - GxRho * (e_tidal_eq(i,j) + e_tidal_sal(i,j)) - enddo ; enddo - enddo - endif + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + e(i,j,K) = e(i,j,K) - (e_tidal_eq(i,j) + e_tidal_sal(i,j)) + pa(i,j,K) = pa(i,j,K) - GxRho * (e_tidal_eq(i,j) + e_tidal_sal(i,j)) + enddo ; enddo + enddo ; endif endif if (CS%correction_intxpa .or. CS%reset_intxpa_integral) then @@ -1604,6 +1607,34 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm ((h(i,j,k) + h(i,j+1,k)) + h_neglect)) enddo ; enddo ; enddo + ! Calculate SAL geopotential anomaly and add its gradient to pressure gradient force + if (CS%calculate_SAL .and. CS%tides_answer_date>20230630 .and. CS%bq_sal_tides) then + !$OMP parallel do default(shared) + do k=1,nz + do j=js,je ; do I=Isq,Ieq + PFu(I,j,k) = PFu(I,j,k) + (e_sal(i+1,j) - e_sal(i,j)) * GV%g_Earth * G%IdxCu(I,j) + enddo ; enddo + do J=Jsq,Jeq ; do i=is,ie + PFv(i,J,k) = PFv(i,J,k) + (e_sal(i,j+1) - e_sal(i,j)) * GV%g_Earth * G%IdyCv(i,J) + enddo ; enddo + enddo + endif + + ! Calculate tidal geopotential anomaly and add its gradient to pressure gradient force + if (CS%tides .and. CS%tides_answer_date>20230630 .and. CS%bq_sal_tides) then + !$OMP parallel do default(shared) + do k=1,nz + do j=js,je ; do I=Isq,Ieq + PFu(I,j,k) = PFu(I,j,k) + ((e_tidal_eq(i+1,j) + e_tidal_sal(i+1,j)) & + - (e_tidal_eq(i,j) + e_tidal_sal(i,j))) * GV%g_Earth * G%IdxCu(I,j) + enddo ; enddo + do J=Jsq,Jeq ; do i=is,ie + PFv(i,J,k) = PFv(i,J,k) + ((e_tidal_eq(i,j+1) + e_tidal_sal(i,j+1)) & + - (e_tidal_eq(i,j) + e_tidal_sal(i,j))) * GV%g_Earth * G%IdyCv(i,J) + enddo ; enddo + enddo + endif + if (CS%GFS_scale < 1.0) then ! Adjust the Montgomery potential to make this a reduced gravity model. if (use_EOS) then @@ -1649,7 +1680,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 eta(i,j) = e(i,j,1)*GV%Z_to_H enddo ; enddo - if (CS%tides) then + if (CS%tides .and. (.not.CS%bq_sal_tides)) then if (CS%tides_answer_date>20230630) then !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 @@ -1662,7 +1693,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm enddo ; enddo endif endif - if (CS%calculate_SAL .and. (CS%tides_answer_date>20230630)) then + if (CS%calculate_SAL .and. (CS%tides_answer_date>20230630) .and. (.not.CS%bq_sal_tides)) then !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 eta(i,j) = eta(i,j) + e_sal(i,j)*GV%Z_to_H @@ -1709,7 +1740,12 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm ! To be consistent with old runs, tidal forcing diagnostic also includes total SAL. ! New diagnostics are given for each individual field. - if (CS%id_e_tide>0) call post_data(CS%id_e_tide, e_sal+e_tidal_eq+e_tidal_sal, CS%diag) + if (CS%id_e_tide>0) then + if (CS%tides_answer_date>20230630) then ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + e_sal_and_tide(i,j) = e_sal(i,j) + e_tidal_eq(i,j) + e_tidal_sal(i,j) + enddo ; enddo ; endif + call post_data(CS%id_e_tide, e_sal_and_tide, CS%diag) + endif if (CS%id_e_sal>0) call post_data(CS%id_e_sal, e_sal, CS%diag) if (CS%id_e_tidal_eq>0) call post_data(CS%id_e_tidal_eq, e_tidal_eq, CS%diag) if (CS%id_e_tidal_sal>0) call post_data(CS%id_e_tidal_sal, e_tidal_sal, CS%diag) @@ -1780,6 +1816,14 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp, if (CS%calculate_SAL) & call get_param(param_file, '', "SAL_USE_BPA", CS%sal_use_bpa, default=.false., & do_not_log=.true.) + if ((CS%tides .or. CS%calculate_SAL) .and. GV%Boussinesq) & + call get_param(param_file, mdl, "BOUSSINESQ_SAL_TIDES", CS%bq_sal_tides, "If true, "//& + "in Boussinesq mode, use an alternative method to include self-attraction "//& + "and loading (SAL) and tidal forcings in pressure gradient, in which their "//& + "gradients are calculated separately, instead of adding geopotential "//& + "anomalies as corrections to the interface height. This alternative method "//& + "elimates a baroclinic component of the SAL and tidal forcings.", & + default=.false.) call get_param(param_file, "MOM", "ENABLE_THERMODYNAMICS", use_temperature, & "If true, Temperature and salinity are used as state variables.", & default=.true., do_not_log=.true.) diff --git a/src/parameterizations/lateral/MOM_self_attr_load.F90 b/src/parameterizations/lateral/MOM_self_attr_load.F90 index fe4149d602..cde06505e6 100644 --- a/src/parameterizations/lateral/MOM_self_attr_load.F90 +++ b/src/parameterizations/lateral/MOM_self_attr_load.F90 @@ -45,8 +45,8 @@ module MOM_self_attr_load real, allocatable :: Love_scaling(:) !< Dimensional coefficients for harmonic SAL, which are functions of Love numbers !! [nondim or L2 Z-1 T-2 ~> m s-2 or R-1 ~> m3 kg-1 or L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1] - real, allocatable :: Snm_Re(:), Snm_Im(:) - !< Real and imaginary coefficients for harmonic SAL [Z ~> m or L2 T-2 ~> m2 s-2] + real, allocatable :: Snm_Re(:), & !< Real SHT coefficient for SHT SAL [Z ~> m] + Snm_Im(:) !< Imaginary SHT coefficient for SHT SAL [Z ~> m] end type SAL_CS integer :: id_clock_SAL !< CPU clock for self-attraction and loading From 36f763c4f877e4ef3a5e81dd8a5bf90c07d57f61 Mon Sep 17 00:00:00 2001 From: He Wang Date: Mon, 7 Oct 2024 10:19:39 -0400 Subject: [PATCH 05/17] Remove rescaling bottom pressure in SAL * Instead of rescaling bottom pressure to height unit, calc_Loving_scaling is modified to be conditionally dimensional. When calculating self-attraction and loading, Love numbers are now dimensional when bottom pressure anomaly is used as an input. This change eliminate Love numbers' dependence on mean seawater density. * A new coefficient called linear_scaling is added to SAL CS for the same purpose, although to use bottom pressure anomaly for scalar approximation is not quite justifiable. A WARNING is given when users try to do that. * Two separate field, SSH (sea surface height anomaly) and pbot (total bottom pressure), are used for calculating self attraction and loading when SAL_USE_BPA = False and SAL_USE_BPA = True, respectively. The change is to enhance code readability. --- src/core/MOM_PressureForce_FV.F90 | 36 +++++---- .../lateral/MOM_self_attr_load.F90 | 80 ++++++++++--------- .../lateral/MOM_spherical_harmonics.F90 | 4 +- 3 files changed, 64 insertions(+), 56 deletions(-) diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90 index 7a96de15ac..7accb052a5 100644 --- a/src/core/MOM_PressureForce_FV.F90 +++ b/src/core/MOM_PressureForce_FV.F90 @@ -145,7 +145,10 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ ! the pressure anomaly at the top of the layer [R L4 T-4 ~> Pa m2 s-2]. real, dimension(SZI_(G),SZJ_(G)) :: & dp, & ! The (positive) change in pressure across a layer [R L2 T-2 ~> Pa]. - pbot_anom, & ! Bottom pressure (anomaly) in depth units, used for self-attraction and loading [Z ~> m]. + SSH, & ! Sea surfae height anomaly for self-attraction and loading. Used if + ! CALCULATE_SAL is True and SAL_USE_BPA is False [Z ~> m]. + pbot, & ! Total bottom pressure for self-attraction and loading. Used if + ! CALCULATE_SAL is True and SAL_USE_BPA is True [R L2 T-2 ~> Pa]. e_sal, & ! The bottom geopotential anomaly due to self-attraction and loading [Z ~> m]. e_tidal_eq, & ! The bottom geopotential anomaly due to tidal forces from astronomical sources [Z ~> m]. e_tidal_sal, & ! The bottom geopotential anomaly due to harmonic self-attraction and loading @@ -222,7 +225,6 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ real :: dp_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [R L2 T-2 ~> Pa]. real :: I_gEarth ! The inverse of GV%g_Earth [T2 Z L-2 ~> s2 m-1] - real :: I_g_rho ! The inverse of the density times the gravitational acceleration [Z T2 L-2 R-1 ~> m Pa-1] real :: alpha_anom ! The in-situ specific volume, averaged over a ! layer, less alpha_ref [R-1 ~> m3 kg-1]. logical :: use_p_atm ! If true, use the atmospheric pressure. @@ -408,19 +410,19 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ ! Calculate and add self-attraction and loading (SAL) geopotential height anomaly to interface height. if (CS%calculate_SAL) then if (CS%sal_use_bpa) then - I_g_rho = 1.0 / (GV%rho0*GV%g_Earth) !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - pbot_anom(i,j) = p(i,j,nz+1) * I_g_rho + pbot(i,j) = p(i,j,nz+1) enddo ; enddo + call calc_SAL(pbot, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m) else !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - pbot_anom(i,j) = (za(i,j,1) - alpha_ref*p(i,j,1)) * I_gEarth - G%Z_ref & + SSH(i,j) = (za(i,j,1) - alpha_ref*p(i,j,1)) * I_gEarth - G%Z_ref & - max(-G%bathyT(i,j)-G%Z_ref, 0.0) enddo ; enddo + call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m) endif - call calc_SAL(pbot_anom, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m) ! This gives new answers after the change of separating SAL from tidal forcing module. if ((CS%tides_answer_date>20230630) .or. (.not.GV%semi_Boussinesq) .or. (.not.CS%tides)) then @@ -873,7 +875,10 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm e_tidal_sal, & ! The bottom geopotential anomaly due to harmonic self-attraction and loading ! specific to tides [Z ~> m]. Z_0p, & ! The height at which the pressure used in the equation of state is 0 [Z ~> m] - pbot_anom, & ! Bottom pressure (anomaly) in depth units, used for self-attraction and loading [Z ~> m]. + SSH, & ! Sea surfae height anomaly for self-attraction and loading. Used if + ! CALCULATE_SAL is True and SAL_USE_BPA is False [Z ~> m]. + pbot, & ! Total bottom pressure for self-attraction and loading. Used if + ! CALCULATE_SAL is True and SAL_USE_BPA is True [R L2 T-2 ~> Pa]. dM ! The barotropic adjustment to the Montgomery potential to ! account for a reduced gravity model [L2 T-2 ~> m2 s-2]. real, dimension(SZI_(G)) :: & @@ -1018,7 +1023,6 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm G_Rho0 = GV%g_Earth / GV%Rho0 GxRho = GV%g_Earth * GV%Rho0 rho_ref = CS%Rho0 - I_g_rho = 1.0 / (CS%rho0*GV%g_Earth) ! I think it should be I_g_rho = 1.0 / (GV%rho0*GV%g_Earth) if ((CS%id_MassWt_u > 0) .or. (CS%id_MassWt_v > 0)) then MassWt_u(:,:,:) = 0.0 ; MassWt_v(:,:,:) = 0.0 @@ -1033,17 +1037,17 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm ! The following code is for recovering old answers only. The algorithm moves interface heights before density ! calculations, and therefore is incorrect without SSH_IN_EOS_PRESSURE_FOR_PGF=True (added in August 2024). ! See the code right after Pa calculation loop for the new method. - if (CS%tides_answer_date<=20230630) then + if (CS%tides .and. CS%tides_answer_date<=20230630) then !$OMP parallel do default(shared) do j=Jsq,Jeq+1 do i=Isq,Ieq+1 - pbot_anom(i,j) = min(-G%bathyT(i,j) - G%Z_ref, 0.0) ! reference bottom pressure + SSH(i,j) = min(-G%bathyT(i,j) - G%Z_ref, 0.0) ! reference bottom pressure enddo do k=1,nz ; do i=Isq,Ieq+1 - pbot_anom(i,j) = pbot_anom(i,j) + h(i,j,k)*GV%H_to_Z + SSH(i,j) = SSH(i,j) + h(i,j,k)*GV%H_to_Z enddo ; enddo enddo - call calc_SAL(pbot_anom, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m) + call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m) call calc_tidal_forcing_legacy(CS%Time, e_sal, e_sal_and_tide, e_tidal_eq, e_tidal_sal, & G, US, CS%tides_CSp) !$OMP parallel do default(shared) @@ -1118,6 +1122,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm endif if (CS%use_SSH_in_Z0p .and. use_p_atm) then + I_g_rho = 1.0 / (CS%rho0*GV%g_Earth) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 Z_0p(i,j) = e(i,j,1) + p_atm(i,j) * I_g_rho enddo ; enddo @@ -1201,15 +1206,16 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm if (CS%sal_use_bpa) then !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - pbot_anom(i,j) = pa(i,j,nz+1) * I_g_rho - e(i,j,nz+1) + pbot(i,j) = pa(i,j,nz+1) - GxRho * e(i,j,nz+1) enddo ; enddo + call calc_SAL(pbot, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m) else !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - pbot_anom(i,j) = e(i,j,1) - max(-G%bathyT(i,j) - G%Z_ref, 0.0) ! Remove topography above sea level + SSH(i,j) = e(i,j,1) - max(-G%bathyT(i,j) - G%Z_ref, 0.0) ! Remove topography above sea level enddo ; enddo + call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m) endif - call calc_SAL(pbot_anom, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m) if (.not.CS%bq_sal_tides) then ; do K=1,nz+1 !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 diff --git a/src/parameterizations/lateral/MOM_self_attr_load.F90 b/src/parameterizations/lateral/MOM_self_attr_load.F90 index cde06505e6..1344ca1c04 100644 --- a/src/parameterizations/lateral/MOM_self_attr_load.F90 +++ b/src/parameterizations/lateral/MOM_self_attr_load.F90 @@ -34,8 +34,7 @@ module MOM_self_attr_load real :: eta_prop !< The partial derivative of eta_sal with the local value of eta [nondim]. real :: linear_scaling - !< A dimensional coefficient for scalar approximation SAL, equal to eta_prop * unit_convert - !! [nondim or L2 Z-1 T-2 ~> m s-2 or R-1 ~> m3 kg-1 or L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1] + !< Dimensional coefficients for scalar SAL [nondim or Z T2 L-2 R-1 ~> m Pa-1] type(sht_CS), allocatable :: sht !< Spherical harmonic transforms (SHT) control structure integer :: sal_sht_Nd @@ -44,7 +43,7 @@ module MOM_self_attr_load !< Reference bottom pressure scaled by Rho_0 and G_Earth[Z ~> m] real, allocatable :: Love_scaling(:) !< Dimensional coefficients for harmonic SAL, which are functions of Love numbers - !! [nondim or L2 Z-1 T-2 ~> m s-2 or R-1 ~> m3 kg-1 or L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1] + !! [nondim or Z T2 L-2 R-1 ~> m Pa-1] real, allocatable :: Snm_Re(:), & !< Real SHT coefficient for SHT SAL [Z ~> m] Snm_Im(:) !< Imaginary SHT coefficient for SHT SAL [Z ~> m] end type SAL_CS @@ -56,15 +55,13 @@ module MOM_self_attr_load !> This subroutine calculates seawater self-attraction and loading based on either sea surface height (SSH) or bottom !! pressure anomaly. Note that the SAL calculation applies to all motions across the spectrum. Tidal-specific methods !! that assume periodicity, i.e. iterative and read-in SAL, are stored in MOM_tidal_forcing module. -!! The input field is always assume to have the unit of [Z ~> m], which can be either SSH or total bottom pressure -!! scaled by mean seawater density Rho_0 and earth gravity G_Earth. For spherical harmonics, the mean seawater density -!! would be cancelled by the same parameter in Love number scalings. If total bottom pressure is used as input, bottom -!! pressure anomaly is calculated in the subroutine by subtracting a reference pressure from the input bottom pressure. +!! The input field can be either SSH [Z ~> m] or total bottom pressure [R L2 T-2 ~> Pa]. If total bottom pressure is +!! used, bottom pressure anomaly is first calculated by subtracting a reference bottom pressure from an input file. !! The output field is expressed as geopotential height anomaly, and therefore has the unit of [Z ~> m]. subroutine calc_SAL(eta, eta_sal, G, CS, tmp_scale) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. real, dimension(SZI_(G),SZJ_(G)), intent(in) :: eta !< The sea surface height anomaly from - !! a time-mean geoid or total bottom pressure scaled by mean density and earth gravity [Z ~> m]. + !! a time-mean geoid or total bottom pressure [Z ~> m] or [R L2 T-2 ~> Pa]. real, dimension(SZI_(G),SZJ_(G)), intent(out) :: eta_sal !< The geopotential height anomaly from !! self-attraction and loading [Z ~> m]. type(SAL_CS), intent(inout) :: CS !< The control structure returned by a previous call to SAL_init. @@ -72,8 +69,8 @@ subroutine calc_SAL(eta, eta_sal, G, CS, tmp_scale) !! to MKS units in reproducing sumes [m Z-1 ~> 1] ! Local variables + real, dimension(SZI_(G),SZJ_(G)) :: bpa ! SSH or bottom pressure anomaly [Z ~> m or R L2 T-2 ~> Pa] integer :: n, m, l - real, dimension(SZI_(G),SZJ_(G)) :: bpa ! [Z ~> m] integer :: Isq, Ieq, Jsq, Jeq integer :: i, j @@ -90,7 +87,7 @@ subroutine calc_SAL(eta, eta_sal, G, CS, tmp_scale) ! use the scalar approximation and/or iterative tidal SAL if (CS%use_sal_scalar .or. CS%use_tidal_sal_prev) then do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - eta_sal(i,j) = CS%eta_prop * bpa(i,j) + eta_sal(i,j) = CS%linear_scaling * bpa(i,j) enddo ; enddo ! use the spherical harmonics method @@ -132,21 +129,21 @@ end subroutine scalar_SAL_sensitivity !> This subroutine calculates coefficients of the spherical harmonic modes for self-attraction and loading. !! The algorithm is based on the SAL implementation in MPAS-ocean, which was modified by Kristin Barton from !! routine written by K. Quinn (March 2010) and modified by M. Schindelegger (May 2017). -subroutine calc_love_scaling(nlm, rhoW, rhoE, Love_scaling) - integer, intent(in) :: nlm !< Maximum spherical harmonics degree [nondim] - real, intent(in) :: rhoW !< The average density of sea water [R ~> kg m-3] - real, intent(in) :: rhoE !< The average density of Earth [R ~> kg m-3] - real, dimension(:), intent(out) :: Love_scaling !< Scaling factors for inverse spherical harmonic - !! transforms [nondim] +subroutine calc_love_scaling(rhoW, rhoE, grav, CS) + real, intent(in) :: rhoW !< The average density of sea water [R ~> kg m-3] + real, intent(in) :: rhoE !< The average density of Earth [R ~> kg m-3] + real, intent(in) :: grav !< The gravitational acceleration [L2 Z-1 T-2 ~> m s-2] + type(SAL_CS), intent(inout) :: CS !< The control structure returned by a previous call to SAL_init. ! Local variables real, dimension(:), allocatable :: HDat, LDat, KDat ! Love numbers converted in CF reference frames [nondim] real :: H1, L1, K1 ! Temporary variables to store degree 1 Love numbers [nondim] - integer :: n_tot ! Size of the stored Love numbers + integer :: n_tot ! Size of the stored Love numbers [nondim] + integer :: nlm ! Maximum spherical harmonics degree [nondim] integer :: n, m, l - real :: unit n_tot = size(Love_Data, dim=2) + nlm = CS%sal_sht_Nd if (nlm+1 > n_tot) call MOM_error(FATAL, "MOM_tidal_forcing " // & "calc_love_scaling: maximum spherical harmonics degree is larger than " // & @@ -165,7 +162,11 @@ subroutine calc_love_scaling(nlm, rhoW, rhoE, Love_scaling) do m=0,nlm ; do n=m,nlm l = order2index(m,nlm) - Love_scaling(l+n-m) = (3.0 / real(2*n+1)) * (rhoW / rhoE) * (1.0 + KDat(n+1) - HDat(n+1)) + if (CS%use_bpa) then + CS%Love_scaling(l+n-m) = (3.0 / real(2*n+1)) * (1.0 / (rhoE * grav)) * (1.0 + KDat(n+1) - HDat(n+1)) + else + CS%Love_scaling(l+n-m) = (3.0 / real(2*n+1)) * (rhoW / rhoE) * (1.0 + KDat(n+1) - HDat(n+1)) + endif enddo ; enddo end subroutine calc_love_scaling @@ -184,14 +185,12 @@ subroutine SAL_init(G, GV, US, param_file, CS) real :: rhoE ! The average density of Earth [R ~> kg m-3]. character(len=200) :: filename, ebot_ref_file, inputdir ! Strings for file/path character(len=200) :: ebot_ref_varname ! Variable name in file - real :: I_g_rho ! The inverse of the density times the gravitational acceleration [Z T2 L-2 R-1 ~> m Pa-1] logical :: calculate_sal=.false. logical :: tides=.false., use_tidal_sal_file=.false., bq_sal_tides_bug=.false. integer :: tides_answer_date=99991203 ! Recover old answers with tides real :: sal_scalar_value, tide_sal_scalar_value ! Scaling SAL factors [nondim] - integer :: isd, ied, jsd, jed - integer :: i, j + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ! Read all relevant parameters and write them to the model log. @@ -230,10 +229,6 @@ subroutine SAL_init(G, GV, US, param_file, CS) allocate(CS%ebot_ref(isd:ied, jsd:jed), source=0.0) call MOM_read_data(filename, trim(ebot_ref_varname), CS%ebot_ref, G%Domain,& scale=US%Pa_to_RL2_T2) - I_g_rho = 1.0 / (GV%g_Earth * GV%Rho0) - do j=jsd,jed ; do i=isd,ied - CS%ebot_ref(i,j) = CS%ebot_ref(i,j) * I_g_rho - enddo ; enddo call pass_var(CS%ebot_ref, G%Domain) endif if (tides_answer_date<=20230630 .and. CS%use_bpa) & @@ -242,9 +237,9 @@ subroutine SAL_init(G, GV, US, param_file, CS) call get_param(param_file, mdl, "SAL_SCALAR_APPROX", CS%use_sal_scalar, & "If true, use the scalar approximation to calculate self-attraction and "//& "loading.", default=tides .and. (.not. use_tidal_sal_file)) - ! if (CS%use_sal_scalar .and. CS%use_bpa) & - ! call MOM_error(WARNING, trim(mdl) // ", SAL_init: Using bottom pressure anomaly for scalar "//& - ! "approximation SAL is unsubstantiated.") + if (CS%use_sal_scalar .and. CS%use_bpa) & + call MOM_error(WARNING, trim(mdl) // ", SAL_init: Using bottom pressure anomaly for scalar "//& + "approximation SAL is unsubstantiated.") call get_param(param_file, '', "TIDE_SAL_SCALAR_VALUE", tide_sal_scalar_value, & units="m m-1", default=0.0, do_not_log=.True.) if (tide_sal_scalar_value/=0.0) & @@ -270,22 +265,29 @@ subroutine SAL_init(G, GV, US, param_file, CS) default=5517.0, scale=US%kg_m3_to_R, do_not_log=(.not. CS%use_sal_sht)) ! Set scaling coefficients for scalar approximation - if (CS%use_sal_scalar .and. CS%use_tidal_sal_prev) then - CS%eta_prop = 2.0 * sal_scalar_value - elseif (CS%use_sal_scalar .or. CS%use_tidal_sal_prev) then - CS%eta_prop = sal_scalar_value + if (CS%use_sal_scalar .or. CS%use_tidal_sal_prev) then + if (CS%use_sal_scalar .and. CS%use_tidal_sal_prev) then + CS%eta_prop = 2.0 * sal_scalar_value + else + CS%eta_prop = sal_scalar_value + endif + if (CS%use_bpa) then + CS%linear_scaling = CS%eta_prop / (GV%Rho0 * GV%g_Earth) + else + CS%linear_scaling = CS%eta_prop + endif else - CS%eta_prop = 0.0 + CS%eta_prop = 0.0 ; CS%linear_scaling = 0.0 endif ! Set scaling coefficients for spherical harmonics if (CS%use_sal_sht) then lmax = calc_lmax(CS%sal_sht_Nd) - allocate(CS%Snm_Re(lmax)); CS%Snm_Re(:) = 0.0 - allocate(CS%Snm_Im(lmax)); CS%Snm_Im(:) = 0.0 + allocate(CS%Snm_Re(lmax), source=0.0) + allocate(CS%Snm_Im(lmax), source=0.0) - allocate(CS%Love_scaling(lmax)); CS%Love_scaling(:) = 0.0 - call calc_love_scaling(CS%sal_sht_Nd, GV%Rho0, rhoE, CS%Love_scaling) + allocate(CS%Love_scaling(lmax), source=0.0) + call calc_love_scaling(GV%Rho0, rhoE, GV%g_Earth, CS) allocate(CS%sht) call spherical_harmonics_init(G, param_file, CS%sht) diff --git a/src/parameterizations/lateral/MOM_spherical_harmonics.F90 b/src/parameterizations/lateral/MOM_spherical_harmonics.F90 index 4f9cee03aa..067e2e3e94 100644 --- a/src/parameterizations/lateral/MOM_spherical_harmonics.F90 +++ b/src/parameterizations/lateral/MOM_spherical_harmonics.F90 @@ -1,12 +1,12 @@ !> Laplace's spherical harmonic transforms (SHT) module MOM_spherical_harmonics +use MOM_coms_infra, only : sum_across_PEs +use MOM_coms, only : reproducing_sum use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, & CLOCK_MODULE, CLOCK_ROUTINE, CLOCK_LOOP use MOM_error_handler, only : MOM_error, FATAL use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type -use MOM_coms_infra, only : sum_across_PEs -use MOM_coms, only : reproducing_sum implicit none ; private From 73a625df738d2ee5aed9f2b1a5f21b8aaca1c944 Mon Sep 17 00:00:00 2001 From: He Wang Date: Tue, 7 Jan 2025 14:50:07 -0500 Subject: [PATCH 06/17] Add answer date flag for tide/SAL A new date for TIDES_ANSWER_DATE is added to recover answers for tides in Boussinesq mode before 20250201. --- src/core/MOM_PressureForce_FV.F90 | 71 ++++++++++++------- .../lateral/MOM_self_attr_load.F90 | 6 +- 2 files changed, 50 insertions(+), 27 deletions(-) diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90 index 7accb052a5..cccfa1c2e4 100644 --- a/src/core/MOM_PressureForce_FV.F90 +++ b/src/core/MOM_PressureForce_FV.F90 @@ -1032,28 +1032,49 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm e(i,j,nz+1) = -G%bathyT(i,j) enddo ; enddo - ! Self-attraction and loading (SAL) and tides (old answers) - ! SAL and tidal geopotential height anomalies are calculated and added as corrections to interface heights. - ! The following code is for recovering old answers only. The algorithm moves interface heights before density - ! calculations, and therefore is incorrect without SSH_IN_EOS_PRESSURE_FOR_PGF=True (added in August 2024). - ! See the code right after Pa calculation loop for the new method. - if (CS%tides .and. CS%tides_answer_date<=20230630) then + ! The following two if-blocks are used to recover old answers for self-attraction and loading + ! (SAL) and tides only. The old algorithm moves interface heights before density calculations, + ! and therefore is incorrect without SSH_IN_EOS_PRESSURE_FOR_PGF=True (added in August 2024). + ! See the code right after Pa calculation loop for the new algorithm. + + ! Calculate and add SAL geopotential anomaly to interface height (old answers) + if (CS%calculate_SAL .and. CS%tides_answer_date<=20250131) then !$OMP parallel do default(shared) do j=Jsq,Jeq+1 do i=Isq,Ieq+1 - SSH(i,j) = min(-G%bathyT(i,j) - G%Z_ref, 0.0) ! reference bottom pressure + SSH(i,j) = min(-G%bathyT(i,j) - G%Z_ref, 0.0) enddo do k=1,nz ; do i=Isq,Ieq+1 SSH(i,j) = SSH(i,j) + h(i,j,k)*GV%H_to_Z enddo ; enddo enddo call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m) - call calc_tidal_forcing_legacy(CS%Time, e_sal, e_sal_and_tide, e_tidal_eq, e_tidal_sal, & - G, US, CS%tides_CSp) - !$OMP parallel do default(shared) - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - e(i,j,nz+1) = e(i,j,nz+1) - e_sal_and_tide(i,j) - enddo ; enddo + + if (CS%tides_answer_date>20230630) then ! answers_date between [20230701, 20250131] + !$OMP parallel do default(shared) + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + e(i,j,nz+1) = e(i,j,nz+1) - e_sal(i,j) + enddo ; enddo + endif + endif + + ! Calculate and add tidal geopotential anomaly to interface height (old answers) + if (CS%tides .and. CS%tides_answer_date<=20250131) then + if (CS%tides_answer_date>20230630) then ! answers_date between [20230701, 20250131] + call calc_tidal_forcing(CS%Time, e_tidal_eq, e_tidal_sal, G, US, CS%tides_CSp) + !$OMP parallel do default(shared) + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + e(i,j,nz+1) = e(i,j,nz+1) - (e_tidal_eq(i,j) + e_tidal_sal(i,j)) + enddo ; enddo + else ! answers_date before 20230701 + if (.not.CS%calculate_SAL) e_sal(:,:) = 0.0 + call calc_tidal_forcing_legacy(CS%Time, e_sal, e_sal_and_tide, e_tidal_eq, e_tidal_sal, & + G, US, CS%tides_CSp) + !$OMP parallel do default(shared) + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + e(i,j,nz+1) = e(i,j,nz+1) - e_sal_and_tide(i,j) + enddo ; enddo + endif endif !$OMP parallel do default(shared) @@ -1201,8 +1222,8 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm enddo ; enddo enddo - ! Self-attraction and loading (SAL) (new answers) - if (CS%calculate_SAL .and. CS%tides_answer_date>20230630) then + ! Calculate and add SAL geopotential anomaly to interface height (new answers) + if (CS%calculate_SAL .and. CS%tides_answer_date>20250131) then if (CS%sal_use_bpa) then !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 @@ -1225,8 +1246,8 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm enddo ; endif endif - ! Tides (new answers) - if (CS%tides .and. CS%tides_answer_date>20230630) then + ! Calculate and add tidal geopotential anomaly to interface height (new answers) + if (CS%tides .and. CS%tides_answer_date>20250131) then call calc_tidal_forcing(CS%Time, e_tidal_eq, e_tidal_sal, G, US, CS%tides_CSp) if (.not.CS%bq_sal_tides) then ; do K=1,nz+1 !$OMP parallel do default(shared) @@ -1812,11 +1833,13 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp, call get_param(param_file, '', "DEFAULT_ANSWER_DATE", default_answer_date, default=99991231) if (CS%tides) & call get_param(param_file, mdl, "TIDES_ANSWER_DATE", CS%tides_answer_date, "The vintage of "//& - "self-attraction and loading (SAL) and tidal forcing calculations. In both "//& - "Boussinesq and non-Boussinesq modes, values below 20230701 recover the old "//& - "answers in which the SAL is part of the tidal forcing calculation. The "//& - "change is due to a reordered summation and the difference is only at bit "//& - "level.", default=20230630, do_not_log=(.not.CS%tides)) + "self-attraction and loading (SAL) and tidal forcing calculations. Setting "//& + "dates before 20230701 recovers old answers (Boussinesq and non-Boussinesq "//& + "modes) when SAL is part of the tidal forcing calculation. The answer "//& + "difference is only at bit level and due to a reordered summation. Setting "//& + "dates before 20250201 recovers answers (Boussinesq mode) that interface "//& + "heights are modified before pressure force integrals are calculated.", & + default=20230630, do_not_log=(.not.CS%tides)) call get_param(param_file, mdl, "CALCULATE_SAL", CS%calculate_SAL, & "If true, calculate self-attraction and loading.", default=CS%tides) if (CS%calculate_SAL) & @@ -1844,9 +1867,9 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp, "pressure used in the equation of state calculations for the Boussinesq pressure "//& "gradient forces, including adjustments for atmospheric or sea-ice pressure.", & default=.false., do_not_log=.not.GV%Boussinesq) - if (CS%tides .and. CS%tides_answer_date<=20230630 .and. CS%use_SSH_in_Z0p) & + if (CS%tides .and. CS%tides_answer_date<=20250131 .and. CS%use_SSH_in_Z0p) & call MOM_error(FATAL, trim(mdl) // ", PressureForce_FV_init: SSH_IN_EOS_PRESSURE_FOR_PGF "//& - "needs to be FALSE to recover tide answers before 20230701.") + "needs to be FALSE to recover tide answers before 20250131.") call get_param(param_file, "MOM", "USE_REGRIDDING", use_ALE, & "If True, use the ALE algorithm (regridding/remapping). "//& diff --git a/src/parameterizations/lateral/MOM_self_attr_load.F90 b/src/parameterizations/lateral/MOM_self_attr_load.F90 index 1344ca1c04..fa5a46110c 100644 --- a/src/parameterizations/lateral/MOM_self_attr_load.F90 +++ b/src/parameterizations/lateral/MOM_self_attr_load.F90 @@ -187,7 +187,7 @@ subroutine SAL_init(G, GV, US, param_file, CS) character(len=200) :: ebot_ref_varname ! Variable name in file logical :: calculate_sal=.false. logical :: tides=.false., use_tidal_sal_file=.false., bq_sal_tides_bug=.false. - integer :: tides_answer_date=99991203 ! Recover old answers with tides + integer :: tides_answer_date=99991231 ! Recover old answers with tides real :: sal_scalar_value, tide_sal_scalar_value ! Scaling SAL factors [nondim] integer :: isd, ied, jsd, jed @@ -231,9 +231,9 @@ subroutine SAL_init(G, GV, US, param_file, CS) scale=US%Pa_to_RL2_T2) call pass_var(CS%ebot_ref, G%Domain) endif - if (tides_answer_date<=20230630 .and. CS%use_bpa) & + if (tides_answer_date<=20250131 .and. CS%use_bpa) & call MOM_error(FATAL, trim(mdl) // ", SAL_init: SAL_USE_BPA needs to be false to recover "//& - "tide answers before 20230630.") + "tide answers before 20250131.") call get_param(param_file, mdl, "SAL_SCALAR_APPROX", CS%use_sal_scalar, & "If true, use the scalar approximation to calculate self-attraction and "//& "loading.", default=tides .and. (.not. use_tidal_sal_file)) From 8e0b83cdaac726b5f77cfb6cdfee9c30c9cf8531 Mon Sep 17 00:00:00 2001 From: He Wang Date: Mon, 13 Jan 2025 09:40:14 -0500 Subject: [PATCH 07/17] Refactor Love_scaling calculation in SAL module * Precalcualte a local field `coef_rhoE` to avoid in-loop division and if-blocks. The unit of coef_rhoE depends on use_bpa. * Fix a few unit description typos in SAL module and two other files. --- src/core/MOM_open_boundary.F90 | 2 +- .../lateral/MOM_self_attr_load.F90 | 25 ++++++++++++------- .../vertical/MOM_diapyc_energy_req.F90 | 2 +- 3 files changed, 18 insertions(+), 11 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index c15b6bd54b..fb833189ec 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -5470,7 +5470,7 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg) endif I_scale = 1.0 ; if (segment%tr_Reg%Tr(m)%scale /= 0.0) I_scale = 1.0 / segment%tr_Reg%Tr(m)%scale if (allocated(segment%tr_Reg%Tr(m)%tres)) then ; do k=1,nz - ! Calculate weights. Both a and u_L are nodim. Adding them together has no meaning. + ! Calculate weights. Both a and u_L are nondim. Adding them together has no meaning. ! However, since they cannot be both non-zero, adding them works like a switch. ! When InvLscale_out is 0 and outflow, only interior data is applied to reservoirs ! When InvLscale_in is 0 and inflow, only nudged data is applied to reservoirs diff --git a/src/parameterizations/lateral/MOM_self_attr_load.F90 b/src/parameterizations/lateral/MOM_self_attr_load.F90 index fa5a46110c..9aa325cfb8 100644 --- a/src/parameterizations/lateral/MOM_self_attr_load.F90 +++ b/src/parameterizations/lateral/MOM_self_attr_load.F90 @@ -38,12 +38,12 @@ module MOM_self_attr_load type(sht_CS), allocatable :: sht !< Spherical harmonic transforms (SHT) control structure integer :: sal_sht_Nd - !< Maximum degree for spherical harmonic transforms [nodim] + !< Maximum degree for spherical harmonic transforms [nondim] real, allocatable :: ebot_ref(:,:) !< Reference bottom pressure scaled by Rho_0 and G_Earth[Z ~> m] real, allocatable :: Love_scaling(:) !< Dimensional coefficients for harmonic SAL, which are functions of Love numbers - !! [nondim or Z T2 L-2 R-1 ~> m Pa-1] + !! [nondim] or [Z T2 L-2 R-1 ~> m Pa-1], depending on the value of use_ppa. real, allocatable :: Snm_Re(:), & !< Real SHT coefficient for SHT SAL [Z ~> m] Snm_Im(:) !< Imaginary SHT coefficient for SHT SAL [Z ~> m] end type SAL_CS @@ -69,7 +69,7 @@ subroutine calc_SAL(eta, eta_sal, G, CS, tmp_scale) !! to MKS units in reproducing sumes [m Z-1 ~> 1] ! Local variables - real, dimension(SZI_(G),SZJ_(G)) :: bpa ! SSH or bottom pressure anomaly [Z ~> m or R L2 T-2 ~> Pa] + real, dimension(SZI_(G),SZJ_(G)) :: bpa ! SSH or bottom pressure anomaly [Z ~> m] or [R L2 T-2 ~> Pa] integer :: n, m, l integer :: Isq, Ieq, Jsq, Jeq integer :: i, j @@ -136,6 +136,8 @@ subroutine calc_love_scaling(rhoW, rhoE, grav, CS) type(SAL_CS), intent(inout) :: CS !< The control structure returned by a previous call to SAL_init. ! Local variables + real :: coef_rhoE ! A scaling coefficient of solid Earth density. coef_rhoE = rhoW / rhoE with USE_BPA=False + ! and coef_rhoE = 1.0 / (rhoE * grav) with USE_BPA=True. [nondim] or [Z T2 L-2 R-1 ~> m Pa-1] real, dimension(:), allocatable :: HDat, LDat, KDat ! Love numbers converted in CF reference frames [nondim] real :: H1, L1, K1 ! Temporary variables to store degree 1 Love numbers [nondim] integer :: n_tot ! Size of the stored Love numbers [nondim] @@ -160,13 +162,16 @@ subroutine calc_love_scaling(rhoW, rhoE, grav, CS) KDat(2) = (-1.0 / 3.0) * H1 - (2.0 / 3.0) * L1 - 1.0 endif + if (CS%use_bpa) then + coef_rhoE = 1.0 / (rhoE * grav) ! [Z T2 L-2 R-1 ~> m Pa-1] + else + coef_rhoE = rhoW / rhoE ! [nondim] + endif + do m=0,nlm ; do n=m,nlm - l = order2index(m,nlm) - if (CS%use_bpa) then - CS%Love_scaling(l+n-m) = (3.0 / real(2*n+1)) * (1.0 / (rhoE * grav)) * (1.0 + KDat(n+1) - HDat(n+1)) - else - CS%Love_scaling(l+n-m) = (3.0 / real(2*n+1)) * (rhoW / rhoE) * (1.0 + KDat(n+1) - HDat(n+1)) - endif + l = order2index(m, nlm) + ! Love_scaling has the same as coef_rhoE. + CS%Love_scaling(l+n-m) = (3.0 / real(2*n+1)) * coef_rhoE * (1.0 + KDat(n+1) - HDat(n+1)) enddo ; enddo end subroutine calc_love_scaling @@ -315,6 +320,8 @@ end subroutine SAL_end !> \namespace self_attr_load !! +!! \section section_SAL Self attraction and loading +!! !! This module contains methods to calculate self-attraction and loading (SAL) as a function of sea surface height or !! bottom pressure anomaly. SAL is primarily used for fast evolving processes like tides or storm surges, but the !! effect applies to all motions. diff --git a/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 b/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 index 7ca432fea4..8ed75a9f44 100644 --- a/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 +++ b/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 @@ -73,7 +73,7 @@ subroutine diapyc_energy_req_test(h_3d, dt, tv, G, GV, US, CS, Kd_int) Kd, & ! A column of diapycnal diffusivities at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]. h_top, h_bot ! Distances from the top or bottom [H ~> m or kg m-2]. real :: dz_h_int ! The ratio of the vertical distances across the layers surrounding an interface - ! over the layer thicknesses [H Z-1 ~> nonodim or kg m-3] + ! over the layer thicknesses [H Z-1 ~> nondim or kg m-3] real :: ustar ! The local friction velocity [Z T-1 ~> m s-1] real :: absf ! The absolute value of the Coriolis parameter [T-1 ~> s-1] real :: htot ! The sum of the thicknesses [H ~> m or kg m-2]. From 8249510a23e9815046a867f47d21b000dac23671 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 11 Dec 2024 09:36:26 -0500 Subject: [PATCH 08/17] +Refactor MOM_opacity to replace hard-coded params Refactored MOM_opacity to replace hard-coded dimensional parameters for the Manizza and Morel opacity fits with run-time parameters, and also added the runtime parameter OPACITY_BAND_WAVELENGTHS to provide the ability to set the wavelengths of the bands, even though these are not actually used in MOM6. By default, these parameters are all set to the previous hard-coded values, using the recently added defaults argument to get_param_real_array(). The bounds of the frequency band label arrays with the MANIZZA_05 opacity scheme were also corrected when PEN_SW_NBANDS > 3, but it would not be typical to use so many bands for no purpose and these labeling arrays (optics%min_wavelength_band and optics%max_wavelength_band) do not appear to be used anywhere. In addition, the unused publicly visible routines opacity_manizza and opacity_morel were eliminated or made private. All answers are bitwise identical, but there are new entries in some MOM_parameter_doc files. --- .../vertical/MOM_opacity.F90 | 196 +++++++++++++----- 1 file changed, 149 insertions(+), 47 deletions(-) diff --git a/src/parameterizations/vertical/MOM_opacity.F90 b/src/parameterizations/vertical/MOM_opacity.F90 index 61a7a0c7d0..abd46d5485 100644 --- a/src/parameterizations/vertical/MOM_opacity.F90 +++ b/src/parameterizations/vertical/MOM_opacity.F90 @@ -17,7 +17,7 @@ module MOM_opacity #include -public set_opacity, opacity_init, opacity_end, opacity_manizza, opacity_morel +public set_opacity, opacity_init, opacity_end public extract_optics_slice, extract_optics_fields, optics_nbands public absorbRemainingSW, sumSWoverBands @@ -66,8 +66,21 @@ module MOM_opacity !! radiation that is in the blue band [nondim]. real :: opacity_land_value !< The value to use for opacity over land [Z-1 ~> m-1]. !! The default is 10 m-1 - a value for muddy water. + real, allocatable, dimension(:,:) & + :: opacity_coef !< Groups of coefficients, in [Z-1 ~> m-1] or [Z ~> m] depending on the + !! scheme, in expressions for opacity, with the second index being the + !! wavelength band. For example, when OPACITY_SCHEME = MANIZZA_05, + !! these are coef_1 and coef_2 in the + !! expression opacity = coef_1 + coef_2 * chl**pow. + real, allocatable, dimension(:) & + :: sw_pen_frac_coef !< Coefficients in the expression for the penetrating shortwave + !! fracetion [nondim] + real, allocatable, dimension(:) & + :: chl_power !< Powers of chlorophyll [nondim] for each band for expressions for + !! opacity of the form opacity = coef_1 + coef_2 * chl**pow. type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to !! regulate the timing of diagnostic output. + integer :: chl_dep_bands !< The number of bands that depend on the Chlorophyll concentrations. logical :: warning_issued !< A flag that is used to avoid repetitive warnings. !>@{ Diagnostic IDs @@ -344,9 +357,9 @@ subroutine opacity_from_chl(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir SW_pen_tot = 0.0 if (G%mask2dT(i,j) > 0.0) then if (multiband_vis_input) then - SW_pen_tot = SW_pen_frac_morel(chl_data(i,j)) * (sw_vis_dir(i,j) + sw_vis_dif(i,j)) + SW_pen_tot = SW_pen_frac_morel(chl_data(i,j), CS) * (sw_vis_dir(i,j) + sw_vis_dif(i,j)) elseif (total_sw_input) then - SW_pen_tot = SW_pen_frac_morel(chl_data(i,j)) * 0.5*sw_total(i,j) + SW_pen_tot = SW_pen_frac_morel(chl_data(i,j), CS) * 0.5*sw_total(i,j) endif endif @@ -372,19 +385,21 @@ subroutine opacity_from_chl(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir optics%opacity_band(n,i,j,k) = CS%opacity_land_value enddo else - ! Band 1 is Manizza blue. - optics%opacity_band(1,i,j,k) = (0.0232 + 0.074*chl_data(i,j)**0.674) * US%Z_to_m - if (nbands >= 2) & ! Band 2 is Manizza red. - optics%opacity_band(2,i,j,k) = (0.225 + 0.037*chl_data(i,j)**0.629) * US%Z_to_m - ! All remaining bands are NIR, for lack of something better to do. - do n=3,nbands ; optics%opacity_band(n,i,j,k) = 2.86*US%Z_to_m ; enddo + do n=1,CS%chl_dep_bands + optics%opacity_band(n,i,j,k) = CS%opacity_coef(1,n) + & + CS%opacity_coef(2,n) * chl_data(i,j)**CS%chl_power(n) + enddo + do n=CS%chl_dep_bands+1,optics%nbands ! These bands do not depend on the chlorophyll. + ! Any nonzero values that were in opacity_coef(2,n) have been added to opacity_coef(1,n). + optics%opacity_band(n,i,j,k) = CS%opacity_coef(1,n) + enddo endif enddo ; enddo case (MOREL_88) do j=js,je ; do i=is,ie optics%opacity_band(1,i,j,k) = CS%opacity_land_value if (G%mask2dT(i,j) > 0.0) & - optics%opacity_band(1,i,j,k) = US%Z_to_m * opacity_morel(chl_data(i,j)) + optics%opacity_band(1,i,j,k) = opacity_morel(chl_data(i,j), CS) do n=2,optics%nbands optics%opacity_band(n,i,j,k) = optics%opacity_band(1,i,j,k) @@ -401,28 +416,25 @@ end subroutine opacity_from_chl !> This sets the blue-wavelength opacity according to the scheme proposed by !! Morel and Antoine (1994). -function opacity_morel(chl_data) +function opacity_morel(chl_data, CS) real, intent(in) :: chl_data !< The chlorophyll-A concentration in [mg m-3] - real :: opacity_morel !< The returned opacity [m-1] + type(opacity_CS) :: CS !< Opacity control structure + real :: opacity_morel !< The returned opacity [Z-1 ~> m-1] - ! The following are coefficients for the optical model taken from Morel and - ! Antoine (1994). These coefficients represent a non uniform distribution of - ! chlorophyll-a through the water column. Other approaches may be more - ! appropriate when using an interactive ecosystem model that predicts - ! three-dimensional chl-a values. - real, dimension(6), parameter :: & - Z2_coef = (/7.925, -6.644, 3.662, -1.815, -0.218, 0.502/) ! Extinction length coefficients [m] real :: Chl, Chl2 ! The log10 of chl_data (in mg m-3), and Chl^2 [nondim] Chl = log10(min(max(chl_data,0.02),60.0)) ; Chl2 = Chl*Chl - opacity_morel = 1.0 / ( (Z2_coef(1) + Z2_coef(2)*Chl) + Chl2 * & - ((Z2_coef(3) + Chl*Z2_coef(4)) + Chl2*(Z2_coef(5) + Chl*Z2_coef(6))) ) -end function + ! All frequency bands currently use the same opacities. + opacity_morel = 1.0 / ( (CS%opacity_coef(1,1) + CS%opacity_coef(2,1)*Chl) + Chl2 * & + ((CS%opacity_coef(3,1) + Chl*CS%opacity_coef(4,1)) + & + Chl2*(CS%opacity_coef(5,1) + Chl*CS%opacity_coef(6,1))) ) +end function opacity_morel !> This sets the penetrating shortwave fraction according to the scheme proposed by !! Morel and Antoine (1994). -function SW_pen_frac_morel(chl_data) +function SW_pen_frac_morel(chl_data, CS) real, intent(in) :: chl_data !< The chlorophyll-A concentration [mg m-3] + type(opacity_CS) :: CS !< Opacity control structure real :: SW_pen_frac_morel !< The returned penetrating shortwave fraction [nondim] ! The following are coefficients for the optical model taken from Morel and @@ -431,24 +443,13 @@ function SW_pen_frac_morel(chl_data) ! appropriate when using an interactive ecosystem model that predicts ! three-dimensional chl-a values. real :: Chl, Chl2 ! The log10 of chl_data in mg m-3, and Chl^2 [nondim] - real, dimension(6), parameter :: & - V1_coef = (/0.321, 0.008, 0.132, 0.038, -0.017, -0.007/) ! Penetrating fraction coefficients [nondim] Chl = log10(min(max(chl_data,0.02),60.0)) ; Chl2 = Chl*Chl - SW_pen_frac_morel = 1.0 - ( (V1_coef(1) + V1_coef(2)*Chl) + Chl2 * & - ((V1_coef(3) + Chl*V1_coef(4)) + Chl2*(V1_coef(5) + Chl*V1_coef(6))) ) + SW_pen_frac_morel = 1.0 - ( (CS%SW_pen_frac_coef(1) + CS%SW_pen_frac_coef(2)*Chl) + Chl2 * & + ((CS%SW_pen_frac_coef(3) + Chl*CS%SW_pen_frac_coef(4)) + & + Chl2*(CS%SW_pen_frac_coef(5) + Chl*CS%SW_pen_frac_coef(6))) ) end function SW_pen_frac_morel -!> This sets the blue-wavelength opacity according to the scheme proposed by -!! Manizza, M. et al, 2005. -function opacity_manizza(chl_data) - real, intent(in) :: chl_data !< The chlorophyll-A concentration [mg m-3] - real :: opacity_manizza !< The returned opacity [m-1] -! This sets the blue-wavelength opacity according to the scheme proposed by Manizza, M. et al, 2005. - - opacity_manizza = 0.0232 + 0.074*chl_data**0.674 -end function - !> This subroutine returns a 2-d slice at constant j of fields from an optics_type, with the potential !! for rescaling these fields. subroutine extract_optics_slice(optics, j, G, GV, opacity, opacity_scale, penSW_top, penSW_scale, SpV_avg) @@ -973,9 +974,25 @@ subroutine opacity_init(Time, G, GV, US, param_file, diag, CS, optics) character(len=40) :: scheme_string ! This include declares and sets the variable "version". # include "version_variable.h" + real :: opacity_coefs(6) ! Pairs of opacity coefficients [Z-1 ~> m-1] for blue, red and + ! near-infrared radiation with parameterizations following the + ! functional form from Manizza et al., GRL 2005, namely in the form + ! opacity = coef_1 + coef_2 * chl**pow for each band. + real :: opacity_powers(3) ! Powers of chlorophyll [nondim] for blue, red and near-infrared + ! radiation bands, in expressions for opacity of the form + ! opacity = coef_1 + coef_2 * chl**pow. + real :: extinction_coefs(6) ! Extinction length coefficients [Z ~> m] for penetrating shortwave + ! radiation in the form proposed by Morel and Antoine (1994), namely + ! opacity = 1 / (sum(n=1:6, Coef(n) * log10(Chl)**(n-1))) + real :: sw_pen_frac_coefs(6) ! Coefficients for the shortwave radiation fraction [nondim] in a + ! fifth order polynomial fit as a funciton of log10(Chlorophyll). real :: PenSW_absorb_minthick ! A thickness that is used to absorb the remaining shortwave heat ! flux when that flux drops below PEN_SW_FLUX_ABSORB [H ~> m or kg m-2] - real :: PenSW_minthick_dflt ! The default for PenSW_absorb_minthick [m] + real :: PenSW_minthick_dflt ! The default for PenSW_absorb_minthick [m] + real :: I_NIR_bands ! The inverse of the number of near-infrared bands being used [nondim] + real, allocatable :: band_wavelengths(:) ! The bounding wavelengths for the penetrating shortwave + ! radiation bands [nm] + real, allocatable :: band_wavelen_default(:) ! The defaults for band_wavelengths [nm] integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags integer :: isd, ied, jsd, jed, nz, n isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = GV%ke @@ -1098,25 +1115,104 @@ subroutine opacity_init(Time, G, GV, US, param_file, diag, CS, optics) default=PenSW_minthick_dflt, units="m", scale=GV%m_to_H) optics%PenSW_absorb_Invlen = 1.0 / (PenSW_absorb_minthick + GV%H_subroundoff) + ! The defaults for the following coefficients are taken from Manizza et al., GRL, 2005. + call get_param(param_file, mdl, "OPACITY_VALUES_MANIZZA", opacity_coefs, & + "Pairs of opacity coefficients for blue, red and near-infrared radiation with "//& + "parameterizations following the functional form from Manizza et al., GRL 2005, "//& + "namely in the form opacity = coef_1 + coef_2 * chl**pow for each band. Although "//& + "coefficients are set for 3 bands, more or less bands may actually be used, with "//& + "extra bands following the same properties as band 3.", & + units="m-1", scale=US%Z_to_m, defaults=(/0.0232, 0.074, 0.225, 0.037, 2.86, 0.0/), & + do_not_log=(CS%opacity_scheme/=MANIZZA_05)) + call get_param(param_file, mdl, "CHOROPHYLL_POWER_MANIZZA", opacity_powers, & + "Powers of chlorophyll for blue, red and near-infrared radiation bands in "//& + "expressions for opacity of the form opacity = coef_1 + coef_2 * chl**pow.", & + units="nondim", defaults=(/0.674, 0.629, 0.0/), & + do_not_log=(CS%opacity_scheme/=MANIZZA_05)) + + ! The defaults for the following coefficients are taken from Morel and Antoine (1994). + call get_param(param_file, mdl, "OPACITY_VALUES_MOREL", extinction_coefs, & + "Shortwave extinction length coefficients for shortwave radiation in the form "//& + "proposed by Morel (1988), opacity = 1 / (sum(Coef(n) * log10(Chl)**(n-1))).", & + units="m", scale=US%m_to_Z, defaults=(/7.925, -6.644, 3.662, -1.815, -0.218, 0.502/), & + do_not_log=(CS%opacity_scheme/=MOREL_88)) + call get_param(param_file, mdl, "SW_PEN_FRAC_COEFS_MOREL", sw_pen_frac_coefs, & + "Coefficients for the shortwave radiation fraction in a fifth order polynomial "//& + "fit as a function of log10(Chlorophyll).", & + units="nondim", defaults=(/0.321, 0.008, 0.132, 0.038, -0.017, -0.007/), & + do_not_log=(CS%opacity_scheme/=MOREL_88)) + if (.not.allocated(optics%min_wavelength_band)) & allocate(optics%min_wavelength_band(optics%nbands)) if (.not.allocated(optics%max_wavelength_band)) & allocate(optics%max_wavelength_band(optics%nbands)) + ! Set the wavelengths of the opacity bands + allocate(band_wavelengths(optics%nbands+1), source=0.0) + allocate(band_wavelen_default(optics%nbands+1), source=0.0) if (CS%opacity_scheme == MANIZZA_05) then - optics%min_wavelength_band(1) =0 - optics%max_wavelength_band(1) =550 - if (optics%nbands >= 2) then - optics%min_wavelength_band(2)=550 - optics%max_wavelength_band(2)=700 - endif - if (optics%nbands > 2) then + if (optics%nbands >= 1) band_wavelen_default(2) = 550.0 + if (optics%nbands >= 2) band_wavelen_default(3) = 700.0 + if (optics%nbands >= 3) then + I_NIR_bands = 1.0 / real(optics%nbands - 2) do n=3,optics%nbands - optics%min_wavelength_band(n) =700 - optics%max_wavelength_band(n) =2800 + band_wavelen_default(n+1) = 2800. - (optics%nbands-n)*2100.0*I_NIR_bands enddo endif endif + call get_param(param_file, mdl, "OPACITY_BAND_WAVELENGTHS", band_wavelengths, & + "The bounding wavelengths for the various bands of shortwave radiation, with "//& + "defaults that depend on the setting for OPACITY_SCHEME.", & + units="nm", defaults=band_wavelen_default, do_not_log=(optics%nbands<2)) + do n=1,optics%nbands + optics%min_wavelength_band(n) = band_wavelengths(n) + optics%max_wavelength_band(n) = band_wavelengths(n+1) + enddo + deallocate(band_wavelengths, band_wavelen_default) + + ! Set opacity scheme dependent parameters. + + if (CS%opacity_scheme == MANIZZA_05) then + allocate(CS%opacity_coef(2,optics%nbands)) + allocate(CS%chl_power(optics%nbands)) + do n=1,min(3,optics%nbands) + CS%opacity_coef(1,n) = opacity_coefs(2*n-1) ; CS%opacity_coef(2,n) = opacity_coefs(2*n) + CS%chl_power(n) = opacity_powers(n) + enddo + ! All remaining bands use the same properties as NIR, for lack of something better to do. + do n=4,optics%nbands + CS%opacity_coef(1,n) = CS%opacity_coef(1,n-1) ; CS%opacity_coef(2,n) = CS%opacity_coef(2,n-1) + CS%chl_power(n) = CS%chl_power(n-1) + enddo + ! Determine the last band that is dependent on chlorophyll. + CS%chl_dep_bands = optics%nbands + do n=optics%nbands,1,-1 + if (CS%chl_power(n) /= 0.0) exit + CS%chl_dep_bands = n - 1 + enddo + do n=CS%chl_dep_bands+1,optics%nbands + if (CS%opacity_coef(2,n) /= 0.0) then + call MOM_error(WARNING, "set_opacity: A non-zero value of the chlorophyll dependence in "//& + "OPACITY_VALUES_MANIZZA was set for a band with zero power in its chlorophyll dependence "//& + "as set by CHOROPHYLL_POWER_MANIZZA.") + CS%opacity_coef(1,n) = CS%opacity_coef(1,n) + CS%opacity_coef(2,n) + CS%opacity_coef(2,n) = 0.0 + endif + enddo + + elseif (CS%opacity_scheme == MOREL_88) then + ! The Morel opacity scheme represents a non uniform distribution of chlorophyll-a through the + ! water column. Other approaches may be more appropriate when using an interactive ecosystem + ! model that predicts three-dimensional chl-a values. + allocate(CS%opacity_coef(6, optics%nbands)) + allocate(CS%sw_pen_frac_coef(6)) + + ! As presently implemented, all frequency bands use the same opacities. + do n=1,optics%nbands + CS%opacity_coef(1:6,n) = extinction_coefs(1:6) + enddo + CS%sw_pen_frac_coef(:) = sw_pen_frac_coefs(:) + endif call get_param(param_file, mdl, "OPACITY_LAND_VALUE", CS%opacity_land_value, & "The value to use for opacity over land. The default is "//& @@ -1152,6 +1248,12 @@ subroutine opacity_end(CS, optics) if (allocated(CS%id_opacity)) & deallocate(CS%id_opacity) + if (allocated(CS%opacity_coef)) & + deallocate(CS%opacity_coef) + if (allocated(CS%sw_pen_frac_coef)) & + deallocate(CS%sw_pen_frac_coef) + if (allocated(CS%chl_power)) & + deallocate(CS%chl_power) if (allocated(optics%sw_pen_band)) & deallocate(optics%sw_pen_band) if (allocated(optics%opacity_band)) & From 585d24df2a7efd778ae240abb091f5e6f25b7af9 Mon Sep 17 00:00:00 2001 From: raphael dussin Date: Tue, 14 Jan 2025 12:09:27 -0500 Subject: [PATCH 09/17] implement spatially varying decay rates for internal tides leakage (#794) * add option to specify horizontally varying decay * extend to vertical modes * fix comments * corrected description of decay_rate array --------- Co-authored-by: Raphael Dussin --- .../lateral/MOM_internal_tides.F90 | 46 +++++++++++++++++-- 1 file changed, 41 insertions(+), 5 deletions(-) diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index 899dcbbbf0..65e9bf55f1 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -152,8 +152,10 @@ module MOM_internal_tides integer :: En_restart_power !< A power factor of 2 by which to multiply the energy in restart [nondim] type(time_type), pointer :: Time => NULL() !< A pointer to the model's clock. character(len=200) :: inputdir !< directory to look for coastline angle file - real :: decay_rate !< A constant rate at which internal tide energy is - !! lost to the interior ocean internal wave field [T-1 ~> s-1]. + real, allocatable, dimension(:,:,:,:) :: decay_rate_2d !< rate at which internal tide energy is + !! lost to the interior ocean internal wave field + !! as a function of longitude, latitude, frequency + !! and vertical mode [T-1 ~> s-1]. real :: cdrag !< The bottom drag coefficient [nondim]. real :: drag_min_depth !< The minimum total ocean thickness that will be used in the denominator !! of the quadratic drag terms for internal tides when @@ -174,6 +176,8 @@ module MOM_internal_tides !! internal tide energy [H Z2 T-2 ~> m3 s-2 or J m-2] logical :: apply_residual_drag !< If true, apply sink from residual term of reflection/transmission. + logical :: use_2d_decay_rate + !< If true, use a spatially varying decay rate for each harmonic. real, allocatable :: En(:,:,:,:,:) !< The internal wave energy density as a function of (i,j,angle,frequency,mode) !! integrated within an angular and frequency band [H Z2 T-2 ~> m3 s-2 or J m-2] @@ -677,7 +681,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C ! Calculate loss rate and apply loss over the time step ; apply the same drag timescale ! to each En component (technically not correct; fix later) En_b = CS%En(i,j,a,fr,m) ! save previous value - En_a = CS%En(i,j,a,fr,m) / (1.0 + (dt * CS%decay_rate)) ! implicit update + En_a = CS%En(i,j,a,fr,m) / (1.0 + (dt * CS%decay_rate_2d(i,j,fr,m))) ! implicit update CS%TKE_leak_loss(i,j,a,fr,m) = (En_b - En_a) * I_dt ! compute exact loss rate [H Z2 T-3 ~> m3 s-3 or W m-2] CS%En(i,j,a,fr,m) = En_a ! update value enddo ; enddo ; enddo ; enddo ; enddo @@ -3386,6 +3390,9 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) real :: kappa_itides ! characteristic topographic wave number [L-1 ~> m-1] real, dimension(:,:), allocatable :: ridge_temp ! array for temporary storage of flags ! of cells with double-reflecting ridges [nondim] + real, dimension(:,:), allocatable :: tmp_decay ! a temp array to store decay rates [T-1 ~> s-1] + real :: decay_rate ! A constant rate at which internal tide energy is + ! lost to the interior ocean internal wave field [T-1 ~> s-1]. logical :: use_int_tides, use_temperature logical :: om4_remap_via_sub_cells ! Use the OM4-era ramap_via_sub_cells for calculating the EBT structure real :: IGW_c1_thresh ! A threshold first mode internal wave speed below which all higher @@ -3411,7 +3418,7 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) character(len=200) :: filename character(len=200) :: refl_angle_file character(len=200) :: refl_pref_file, refl_dbl_file, trans_file - character(len=200) :: h2_file + character(len=200) :: h2_file, decay_file character(len=80) :: rough_var ! Input file variable names character(len=240), dimension(:), allocatable :: energy_fractions @@ -3546,10 +3553,13 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) call get_param(param_file, mdl, "MAX_TKE_TO_KD", CS%max_TKE_to_Kd, & "Limiter for TKE_to_Kd.", & units="", default=1e9, scale=US%Z_to_m*US%s_to_T**2) - call get_param(param_file, mdl, "INTERNAL_TIDE_DECAY_RATE", CS%decay_rate, & + call get_param(param_file, mdl, "INTERNAL_TIDE_DECAY_RATE", decay_rate, & "The rate at which internal tide energy is lost to the "//& "interior ocean internal wave field.", & units="s-1", default=0.0, scale=US%T_to_s) + call get_param(param_file, mdl, "USE_2D_INTERNAL_TIDE_DECAY_RATE", CS%use_2d_decay_rate, & + "If true, use a spatially varying decay rate for leakage loss in the "// & + "internal tide code.", default=.false.) call get_param(param_file, mdl, "INTERNAL_TIDE_VOLUME_BASED_CFL", CS%vol_CFL, & "If true, use the ratio of the open face lengths to the "//& "tracer cell areas when estimating CFL numbers in the "//& @@ -3679,6 +3689,32 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) allocate(CS%TKE_itidal_loss_glo_dt(num_freq,num_mode), source=0.0) allocate(CS%TKE_residual_loss_glo_dt(num_freq,num_mode), source=0.0) allocate(CS%TKE_input_glo_dt(num_freq,num_mode), source=0.0) + allocate(CS%decay_rate_2d(isd:ied,jsd:jed,num_freq,num_mode), source=0.0) + allocate(tmp_decay(isd:ied,jsd:jed), source=0.0) + + if (CS%use_2d_decay_rate) then + call get_param(param_file, mdl, "ITIDES_DECAY_FILE", decay_file, & + "The path to the file containing the decay rates "//& + "for internal tides with USE_2D_INTERNAL_TIDE_DECAY_RATE.", & + fail_if_missing=.true.) + do m=1,num_mode ; do fr=1,num_freq + ! read 2d field for each harmonic + filename = trim(CS%inputdir) // trim(decay_file) + write(var_name, '("decay_rate_freq",i1,"_mode",i1)') fr, m + call MOM_read_data(filename, var_name, tmp_decay, G%domain, scale=US%T_to_s) + do j=G%jsc,G%jec ; do i=G%isc,G%iec + CS%decay_rate_2d(i,j,fr,m) = tmp_decay(i,j) + enddo ; enddo + enddo ; enddo + else + do m=1,num_mode ; do fr=1,num_freq ; do j=G%jsc,G%jec ; do i=G%isc,G%iec + CS%decay_rate_2d(i,j,fr,m) = decay_rate + enddo ; enddo ; enddo ; enddo + endif + + do m=1,num_mode + call pass_var(CS%decay_rate_2d(:,:,:,m), G%domain) + enddo ! Compute the fixed part of the bottom drag loss from baroclinic modes call get_param(param_file, mdl, "H2_FILE", h2_file, & From dd1f3f52f5fef0e78ae7216e64d3489f28123cc9 Mon Sep 17 00:00:00 2001 From: claireyung <61528379+claireyung@users.noreply.github.com> Date: Wed, 15 Jan 2025 09:06:36 +1100 Subject: [PATCH 10/17] Add ice shelf pressure initialisation bug fix (#800) * Add ice shelf pressure initialisation bug fix This commit adds a new parameter FRAC_DP_AT_POS_NEGATIVE_P_BUGFIX and associated bug fix. This bug occurs when TRIM_IC_FOR_P_SURF pressure initialisation is used for ice shelf, whilst in Boussinesq mode and a nonlinear equation of state. The subroutine trim_for_ice calls cut_off_column_top. This routine in Boussinesq mode calls find_depth_of_pressure_in_cell in MOM_density_integrals.F90. This subsequently calls the function frac_dp_at_pos which calls the density equation of state based on given salinity, temperature and depth, which is incorrectly converted into pressure by z position (which is negative) multiplied by g and rho0. The bug results in incorrect densities from the equation of state and therefore an imperfect initialisation of layer thicknesses and sea surface height due to the displacement of water by the ice. The bug is fixed by multiplying the pressure by -1. This commit introduces parameter FRAC_DP_AT_POS_NEGATIVE_P_BUGFIX with default value False to preserve answers. If true, the ice shelf initialisation is accurate to a high degree with a nonlinear equation of state. Answers should not change, except for the added parameter in MOM_parameter_doc. The same functions are called by a unit test in MOM_state_initialization; in this commit I set the MOM_state_initialization unit test to use the existing bug (with a false flag). * Resolve indenting and white space inconsitencies with MOM6 style for ice shelf pressure initialisation bug fix FRAC_DP_AT_POS_NEGATIVE_P_BUGFIX --- src/core/MOM_density_integrals.F90 | 16 +++++++++----- .../MOM_state_initialization.F90 | 21 +++++++++++++------ 2 files changed, 26 insertions(+), 11 deletions(-) diff --git a/src/core/MOM_density_integrals.F90 b/src/core/MOM_density_integrals.F90 index 90994dd073..151bc4ef3c 100644 --- a/src/core/MOM_density_integrals.F90 +++ b/src/core/MOM_density_integrals.F90 @@ -2001,7 +2001,7 @@ end subroutine diagnose_mass_weight_p !> Find the depth at which the reconstructed pressure matches P_tgt subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_tgt, & - rho_ref, G_e, EOS, US, P_b, z_out, z_tol) + rho_ref, G_e, EOS, US, P_b, z_out, z_tol, frac_dp_bugfix) real, intent(in) :: T_t !< Potential temperature at the cell top [C ~> degC] real, intent(in) :: T_b !< Potential temperature at the cell bottom [C ~> degC] real, intent(in) :: S_t !< Salinity at the cell top [S ~> ppt] @@ -2020,6 +2020,7 @@ subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_t real, intent(out) :: P_b !< Pressure at the bottom of the cell [R L2 T-2 ~> Pa] real, intent(out) :: z_out !< Absolute depth at which anomalous pressure = p_tgt [Z ~> m] real, intent(in) :: z_tol !< The tolerance in finding z_out [Z ~> m] + logical, intent(in) :: frac_dp_bugfix !< If true, use bugfix in frac_dp_at_pos ! Local variables real :: dp ! Pressure thickness of the layer [R L2 T-2 ~> Pa] @@ -2032,7 +2033,7 @@ subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_t GxRho = G_e * rho_ref ! Anomalous pressure difference across whole cell - dp = frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, 1.0, EOS) + dp = frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, 1.0, EOS, frac_dp_bugfix) P_b = P_t + dp ! Anomalous pressure at bottom of cell @@ -2063,7 +2064,7 @@ subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_t call MOM_error(FATAL, 'find_depth_of_pressure_in_cell completes too many iterations: '//msg) endif z_out = z_t + ( z_b - z_t ) * F_guess - Pa = frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, F_guess, EOS) - ( P_tgt - P_t ) + Pa = frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, F_guess, EOS, frac_dp_bugfix) - ( P_tgt - P_t ) if (Pa Returns change in anomalous pressure change from top to non-dimensional !! position pos between z_t and z_b [R L2 T-2 ~> Pa] -real function frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, pos, EOS) +real function frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, pos, EOS, frac_dp_bugfix) real, intent(in) :: T_t !< Potential temperature at the cell top [C ~> degC] real, intent(in) :: T_b !< Potential temperature at the cell bottom [C ~> degC] real, intent(in) :: S_t !< Salinity at the cell top [S ~> ppt] @@ -2131,6 +2132,7 @@ real function frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, pos, EO real, intent(in) :: G_e !< The Earth's gravitational acceleration [L2 Z-1 T-2 ~> m s-2] real, intent(in) :: pos !< The fractional vertical position, 0 to 1 [nondim] type(EOS_type), intent(in) :: EOS !< Equation of state structure + logical, intent(in) :: frac_dp_bugfix !< If true, use bugfix in frac_dp_at_pos ! Local variables real, parameter :: C1_90 = 1.0/90.0 ! A rational constant [nondim] @@ -2150,7 +2152,11 @@ real function frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, pos, EO ! Salinity and temperature points are linearly interpolated S5(n) = top_weight * S_t + bottom_weight * S_b T5(n) = top_weight * T_t + bottom_weight * T_b - p5(n) = ( top_weight * z_t + bottom_weight * z_b ) * ( G_e * rho_ref ) + if (frac_dp_bugfix) then + p5(n) = (-1) * ( top_weight * z_t + bottom_weight * z_b ) * ( G_e * rho_ref ) + else + p5(n) = ( top_weight * z_t + bottom_weight * z_b ) * ( G_e * rho_ref ) + endif !bugfix enddo call calculate_density(T5, S5, p5, rho5, EOS) rho5(:) = rho5(:) !- rho_ref ! Work with anomalies relative to rho_ref diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index d908ec23a0..7909bddb80 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -1205,6 +1205,7 @@ subroutine trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read) ! answers from 2018, while higher values use more robust ! forms of the same remapping expressions. logical :: use_remapping ! If true, remap the initial conditions. + logical :: use_frac_dp_bugfix ! If true, use bugfix. Otherwise, pressure input to EOS is negative. type(remapping_CS), pointer :: remap_CS => NULL() call get_param(PF, mdl, "SURFACE_PRESSURE_FILE", p_surf_file, & @@ -1227,7 +1228,10 @@ subroutine trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read) "The tolerance with which to find the depth matching the specified "//& "surface pressure with TRIM_IC_FOR_P_SURF.", & units="m", default=1.0e-5, scale=US%m_to_Z, do_not_log=just_read) - + call get_param(PF, mdl, "FRAC_DP_AT_POS_NEGATIVE_P_BUGFIX", use_frac_dp_bugfix, & + "If true, use bugfix in ice shelf TRIM_IC initialization. "//& + "Otherwise, pressure input to density EOS is negative.", & + default=.false., do_not_log=just_read) call get_param(PF, mdl, "TRIMMING_USES_REMAPPING", use_remapping, & 'When trimming the column, also remap T and S.', & default=.false., do_not_log=just_read) @@ -1277,7 +1281,8 @@ subroutine trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read) do j=G%jsc,G%jec ; do i=G%isc,G%iec call cut_off_column_top(GV%ke, tv, GV, US, GV%g_Earth, G%bathyT(i,j)+G%Z_ref, min_thickness, & tv%T(i,j,:), T_t(i,j,:), T_b(i,j,:), tv%S(i,j,:), S_t(i,j,:), S_b(i,j,:), & - p_surf(i,j), h(i,j,:), remap_CS, z_tol=z_tolerance) + p_surf(i,j), h(i,j,:), remap_CS, z_tol=z_tolerance, & + frac_dp_bugfix=use_frac_dp_bugfix) enddo ; enddo end subroutine trim_for_ice @@ -1368,7 +1373,7 @@ end subroutine calc_sfc_displacement !> Adjust the layer thicknesses by removing the top of the water column above the !! depth where the hydrostatic pressure matches p_surf subroutine cut_off_column_top(nk, tv, GV, US, G_earth, depth, min_thickness, T, T_t, T_b, & - S, S_t, S_b, p_surf, h, remap_CS, z_tol) + S, S_t, S_b, p_surf, h, remap_CS, z_tol, frac_dp_bugfix) integer, intent(in) :: nk !< Number of layers type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. @@ -1388,6 +1393,7 @@ subroutine cut_off_column_top(nk, tv, GV, US, G_earth, depth, min_thickness, T, !! if associated real, intent(in) :: z_tol !< The tolerance with which to find the depth !! matching the specified pressure [Z ~> m]. + logical, intent(in) :: frac_dp_bugfix !< If true, use bugfix in frac_dp_at_pos ! Local variables real, dimension(nk+1) :: e ! Top and bottom edge positions for reconstructions [Z ~> m] @@ -1416,7 +1422,8 @@ subroutine cut_off_column_top(nk, tv, GV, US, G_earth, depth, min_thickness, T, do k=1,nk call find_depth_of_pressure_in_cell(T_t(k), T_b(k), S_t(k), S_b(k), e(K), e(K+1), & P_t, p_surf, GV%Rho0, G_earth, tv%eqn_of_state, & - US, P_b, z_out, z_tol=z_tol) + US, P_b, z_out, z_tol=z_tol, & + frac_dp_bugfix=frac_dp_bugfix) if (z_out>=e(K)) then ! Imposed pressure was less that pressure at top of cell exit @@ -3139,7 +3146,8 @@ subroutine MOM_state_init_tests(G, GV, US, tv) P_t = 0. do k = 1, nk call find_depth_of_pressure_in_cell(T_t(k), T_b(k), S_t(k), S_b(k), e(K), e(K+1), P_t, 0.5*P_tot, & - GV%Rho0, GV%g_Earth, tv%eqn_of_state, US, P_b, z_out, z_tol=z_tol) + GV%Rho0, GV%g_Earth, tv%eqn_of_state, US, P_b, z_out, z_tol=z_tol, & + frac_dp_bugfix=.false.) write(0,*) k, US%RL2_T2_to_Pa*P_t, US%RL2_T2_to_Pa*P_b, 0.5*US%RL2_T2_to_Pa*P_tot, & US%Z_to_m*e(K), US%Z_to_m*e(K+1), US%Z_to_m*z_out P_t = P_b @@ -3158,7 +3166,8 @@ subroutine MOM_state_init_tests(G, GV, US, tv) ! h_neglect=GV%H_subroundoff, h_neglect_edge=GV%H_subroundoff) ! endif call cut_off_column_top(nk, tv, GV, US, GV%g_Earth, -e(nk+1), GV%Angstrom_H, & - T, T_t, T_b, S, S_t, S_b, 0.5*P_tot, h, remap_CS, z_tol=z_tol) + T, T_t, T_b, S, S_t, S_b, 0.5*P_tot, h, remap_CS, z_tol=z_tol, & + frac_dp_bugfix=.false.) write(0,*) GV%H_to_m*h(:) if (associated(remap_CS)) deallocate(remap_CS) From d8da512c48a4445114d4b371144f625feace89cd Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 5 Jan 2025 06:22:40 -0500 Subject: [PATCH 11/17] +Add the optional argument old_name to get_param Added the new optional argument old_name to the 8 get_param() routines. This new capability allows for an archaic parameter name to be specified and for appropriate warnings encouraging the user to migrate to using the new name while still setting the parameter as intended, or error messages in the case of inconsistent setting via the archaic name and the correct name. The logging inside of the MOM_parameter_doc files only uses the correct parameter name. Also added the new optional argument set to the 8 read_param routines, to indicate whether a parameter has been found and successfully set. The new set argument is now being used in read_param() calls in obsolete_int(), obsolete_real(), obsolete_char() and obsolete_logical(). Obsolete_logical() in particular was substantially simplified by the use of this new argument, and is now only about half as long as it was. The read_param() set argument is also used in all of the get_param() routines when they are given an old_name argument. The new old_name argument to get_param() is not yet being used in the version of the MOM6 code that is being checked in, but it has been tested extensively by adding or modifying get_param calls in a variant of the initialization code, and it will be used in an updated version of github.com/NOAA-GFDL/MOM6/pull/725 to gracefully handle the deprecation of 4 parameter names. All answers are bitwise identical, but there are new optional arguments to two widely used interfaces. --- src/diagnostics/MOM_obsolete_params.F90 | 37 ++- src/framework/MOM_file_parser.F90 | 288 +++++++++++++++++++++--- 2 files changed, 274 insertions(+), 51 deletions(-) diff --git a/src/diagnostics/MOM_obsolete_params.F90 b/src/diagnostics/MOM_obsolete_params.F90 index a590ae3893..a0401d7199 100644 --- a/src/diagnostics/MOM_obsolete_params.F90 +++ b/src/diagnostics/MOM_obsolete_params.F90 @@ -168,29 +168,15 @@ subroutine obsolete_logical(param_file, varname, warning_val, hint) character(len=*), optional, intent(in) :: hint !< A hint to the user about what to do. ! Local variables logical :: test_logic, fatal_err + logical :: var_is_set ! True if this value was read by read_param. character(len=128) :: hint_msg - test_logic = .false. ; call read_param(param_file, varname, test_logic) + test_logic = .false. ; call read_param(param_file, varname, test_logic, set=var_is_set) fatal_err = .true. - if (present(warning_val)) fatal_err = (warning_val .neqv. .true.) + if (var_is_set .and. present(warning_val)) fatal_err = (warning_val .neqv. test_logic) hint_msg = " " ; if (present(hint)) hint_msg = hint - if (test_logic) then - if (fatal_err) then - call MOM_ERROR(FATAL, "MOM_obsolete_params: "//trim(varname)// & - " is an obsolete run-time flag, and should not be used. "// & - trim(hint_msg)) - else - call MOM_ERROR(WARNING, "MOM_obsolete_params: "//trim(varname)// & - " is an obsolete run-time flag. "//trim(hint_msg)) - endif - endif - - test_logic = .true. ; call read_param(param_file, varname, test_logic) - fatal_err = .true. - if (present(warning_val)) fatal_err = (warning_val .neqv. .false.) - - if (.not.test_logic) then + if (var_is_set) then if (fatal_err) then call MOM_ERROR(FATAL, "MOM_obsolete_params: "//trim(varname)// & " is an obsolete run-time flag, and should not be used. "// & @@ -211,12 +197,13 @@ subroutine obsolete_char(param_file, varname, warning_val, hint) character(len=*), optional, intent(in) :: hint !< A hint to the user about what to do. ! Local variables character(len=200) :: test_string, hint_msg + logical :: var_is_set ! True if this value was read by read_param. logical :: only_warn - test_string = ''; call read_param(param_file, varname, test_string) + test_string = ''; call read_param(param_file, varname, test_string, set=var_is_set) hint_msg = " " ; if (present(hint)) hint_msg = hint - if (len_trim(test_string) > 0) then + if (var_is_set) then only_warn = .false. if (present(warning_val)) then ! Check if test_string and warning_val are the same. if (len_trim(warning_val) == len_trim(test_string)) then @@ -246,15 +233,16 @@ subroutine obsolete_real(param_file, varname, warning_val, hint, only_warn) ! Local variables real :: test_val, warn_val + logical :: var_is_set ! True if this value was read by read_param. logical :: issue_warning character(len=128) :: hint_msg - test_val = -9e35; call read_param(param_file, varname, test_val) + test_val = -9e35; call read_param(param_file, varname, test_val, set=var_is_set) warn_val = -9e35; if (present(warning_val)) warn_val = warning_val hint_msg = " " ; if (present(hint)) hint_msg = hint issue_warning = .false. ; if (present(only_warn)) issue_warning = only_warn - if (test_val /= -9e35) then + if (var_is_set) then if ((test_val == warn_val) .or. issue_warning) then call MOM_ERROR(WARNING, "MOM_obsolete_params: "//trim(varname)// & " is an obsolete run-time flag. "//trim(hint_msg)) @@ -273,14 +261,15 @@ subroutine obsolete_int(param_file, varname, warning_val, hint) integer, optional, intent(in) :: warning_val !< An allowed value that causes a warning instead of an error. character(len=*), optional, intent(in) :: hint !< A hint to the user about what to do. ! Local variables + logical :: var_is_set ! True if this value was read by read_param. integer :: test_val, warn_val character(len=128) :: hint_msg - test_val = -123456788; call read_param(param_file, varname, test_val) + test_val = -123456788; call read_param(param_file, varname, test_val, set=var_is_set) warn_val = -123456788; if (present(warning_val)) warn_val = warning_val hint_msg = " " ; if (present(hint)) hint_msg = hint - if (test_val /= -123456788) then + if (var_is_set) then if (test_val == warn_val) then call MOM_ERROR(WARNING, "MOM_obsolete_params: "//trim(varname)// & " is an obsolete run-time flag. "//trim(hint_msg)) diff --git a/src/framework/MOM_file_parser.F90 b/src/framework/MOM_file_parser.F90 index 7d3337ea24..291d44492d 100644 --- a/src/framework/MOM_file_parser.F90 +++ b/src/framework/MOM_file_parser.F90 @@ -8,7 +8,7 @@ module MOM_file_parser use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, assert use MOM_error_handler, only : is_root_pe, stdlog, stdout use MOM_time_manager, only : get_time, time_type, get_ticks_per_second -use MOM_time_manager, only : set_date, get_date, real_to_time, operator(-), set_time +use MOM_time_manager, only : set_date, get_date, real_to_time, operator(-), operator(==), set_time use MOM_document, only : doc_param, doc_module, doc_init, doc_end, doc_type use MOM_document, only : doc_openBlock, doc_closeBlock use MOM_string_functions, only : left_int, left_ints, slasher @@ -618,7 +618,7 @@ function simplifyWhiteSpace(string) end function simplifyWhiteSpace !> This subroutine reads the value of an integer model parameter from a parameter file. -subroutine read_param_int(CS, varname, value, fail_if_missing) +subroutine read_param_int(CS, varname, value, fail_if_missing, set) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: varname !< The case-sensitive name of the parameter to read @@ -626,6 +626,8 @@ subroutine read_param_int(CS, varname, value, fail_if_missing) !! read from the parameter file logical, optional, intent(in) :: fail_if_missing !< If present and true, a fatal error occurs !! if this variable is not found in the parameter file + logical, optional, intent(out) :: set !< If present, this indicates whether this parameter + !! has been found and successfully set in the input files. ! Local variables character(len=CS%max_line_len) :: value_string(1) logical :: found, defined @@ -633,6 +635,7 @@ subroutine read_param_int(CS, varname, value, fail_if_missing) call get_variable_line(CS, varname, found, defined, value_string) if (found .and. defined .and. (LEN_TRIM(value_string(1)) > 0)) then read(value_string(1),*,err = 1001) value + if (present(set)) set = .true. else if (present(fail_if_missing)) then ; if (fail_if_missing) then if (.not.found) then @@ -643,6 +646,7 @@ subroutine read_param_int(CS, varname, value, fail_if_missing) ' found but not set in input files.') endif endif ; endif + if (present(set)) set = .false. endif return 1001 call MOM_error(FATAL,'read_param_int: read error for integer variable '//trim(varname)// & @@ -650,7 +654,7 @@ subroutine read_param_int(CS, varname, value, fail_if_missing) end subroutine read_param_int !> This subroutine reads the values of an array of integer model parameters from a parameter file. -subroutine read_param_int_array(CS, varname, value, fail_if_missing) +subroutine read_param_int_array(CS, varname, value, fail_if_missing, set) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: varname !< The case-sensitive name of the parameter to read @@ -658,12 +662,15 @@ subroutine read_param_int_array(CS, varname, value, fail_if_missing) !! read from the parameter file logical, optional, intent(in) :: fail_if_missing !< If present and true, a fatal error occurs !! if this variable is not found in the parameter file + logical, optional, intent(out) :: set !< If present, this indicates whether this parameter + !! has been found and successfully set in the input files. ! Local variables character(len=CS%max_line_len) :: value_string(1) logical :: found, defined call get_variable_line(CS, varname, found, defined, value_string) if (found .and. defined .and. (LEN_TRIM(value_string(1)) > 0)) then + if (present(set)) set = .true. read(value_string(1),*,end=991,err=1002) value 991 return else @@ -676,6 +683,7 @@ subroutine read_param_int_array(CS, varname, value, fail_if_missing) ' found but not set in input files.') endif endif ; endif + if (present(set)) set = .false. endif return 1002 call MOM_error(FATAL,'read_param_int_array: read error for integer array '//trim(varname)// & @@ -683,7 +691,7 @@ subroutine read_param_int_array(CS, varname, value, fail_if_missing) end subroutine read_param_int_array !> This subroutine reads the value of a real model parameter from a parameter file. -subroutine read_param_real(CS, varname, value, fail_if_missing, scale) +subroutine read_param_real(CS, varname, value, fail_if_missing, scale, set) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: varname !< The case-sensitive name of the parameter to read @@ -693,6 +701,8 @@ subroutine read_param_real(CS, varname, value, fail_if_missing, scale) !! if this variable is not found in the parameter file real, optional, intent(in) :: scale !< A scaling factor that the parameter is multiplied !! by before it is returned. + logical, optional, intent(out) :: set !< If present, this indicates whether this parameter + !! has been found and successfully set in the input files. ! Local variables character(len=CS%max_line_len) :: value_string(1) @@ -702,6 +712,7 @@ subroutine read_param_real(CS, varname, value, fail_if_missing, scale) if (found .and. defined .and. (LEN_TRIM(value_string(1)) > 0)) then read(value_string(1),*,err=1003) value if (present(scale)) value = scale*value + if (present(set)) set = .true. else if (present(fail_if_missing)) then ; if (fail_if_missing) then if (.not.found) then @@ -712,6 +723,7 @@ subroutine read_param_real(CS, varname, value, fail_if_missing, scale) ' found but not set in input files.') endif endif ; endif + if (present(set)) set = .false. endif return 1003 call MOM_error(FATAL,'read_param_real: read error for real variable '//trim(varname)// & @@ -719,7 +731,7 @@ subroutine read_param_real(CS, varname, value, fail_if_missing, scale) end subroutine read_param_real !> This subroutine reads the values of an array of real model parameters from a parameter file. -subroutine read_param_real_array(CS, varname, value, fail_if_missing, scale) +subroutine read_param_real_array(CS, varname, value, fail_if_missing, scale, set) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: varname !< The case-sensitive name of the parameter to read @@ -729,6 +741,8 @@ subroutine read_param_real_array(CS, varname, value, fail_if_missing, scale) !! if this variable is not found in the parameter file real, optional, intent(in) :: scale !< A scaling factor that the parameter is multiplied !! by before it is returned. + logical, optional, intent(out) :: set !< If present, this indicates whether this parameter + !! has been found and successfully set in the input files. ! Local variables character(len=CS%max_line_len) :: value_string(1) @@ -739,7 +753,7 @@ subroutine read_param_real_array(CS, varname, value, fail_if_missing, scale) read(value_string(1),*,end=991,err=1004) value 991 continue if (present(scale)) value(:) = scale*value(:) - return + if (present(set)) set = .true. else if (present(fail_if_missing)) then ; if (fail_if_missing) then if (.not.found) then @@ -750,6 +764,7 @@ subroutine read_param_real_array(CS, varname, value, fail_if_missing, scale) ' found but not set in input files.') endif endif ; endif + if (present(set)) set = .false. endif return 1004 call MOM_error(FATAL,'read_param_real_array: read error for real array '//trim(varname)// & @@ -757,7 +772,7 @@ subroutine read_param_real_array(CS, varname, value, fail_if_missing, scale) end subroutine read_param_real_array !> This subroutine reads the value of a character string model parameter from a parameter file. -subroutine read_param_char(CS, varname, value, fail_if_missing) +subroutine read_param_char(CS, varname, value, fail_if_missing, set) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: varname !< The case-sensitive name of the parameter to read @@ -765,6 +780,8 @@ subroutine read_param_char(CS, varname, value, fail_if_missing) !! read from the parameter file logical, optional, intent(in) :: fail_if_missing !< If present and true, a fatal error occurs !! if this variable is not found in the parameter file + logical, optional, intent(out) :: set !< If present, this indicates whether this parameter + !! has been found and successfully set in the input files. ! Local variables character(len=CS%max_line_len) :: value_string(1) logical :: found, defined @@ -776,10 +793,12 @@ subroutine read_param_char(CS, varname, value, fail_if_missing) call MOM_error(FATAL, 'Unable to find variable '//trim(varname)//' in any input files.') endif ; endif + if (present(set)) set = found + end subroutine read_param_char !> This subroutine reads the values of an array of character string model parameters from a parameter file. -subroutine read_param_char_array(CS, varname, value, fail_if_missing) +subroutine read_param_char_array(CS, varname, value, fail_if_missing, set) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: varname !< The case-sensitive name of the parameter to read @@ -787,6 +806,8 @@ subroutine read_param_char_array(CS, varname, value, fail_if_missing) !! read from the parameter file logical, optional, intent(in) :: fail_if_missing !< If present and true, a fatal error occurs !! if this variable is not found in the parameter file + logical, optional, intent(out) :: set !< If present, this indicates whether this parameter + !! has been found and successfully set in the input files. ! Local variables character(len=CS%max_line_len) :: value_string(1), loc_string @@ -813,10 +834,12 @@ subroutine read_param_char_array(CS, varname, value, fail_if_missing) call MOM_error(FATAL, 'Unable to find variable '//trim(varname)//' in any input files.') endif ; endif + if (present(set)) set = found + end subroutine read_param_char_array !> This subroutine reads the value of a logical model parameter from a parameter file. -subroutine read_param_logical(CS, varname, value, fail_if_missing) +subroutine read_param_logical(CS, varname, value, fail_if_missing, set) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: varname !< The case-sensitive name of the parameter to read @@ -824,6 +847,8 @@ subroutine read_param_logical(CS, varname, value, fail_if_missing) !! read from the parameter file logical, optional, intent(in) :: fail_if_missing !< If present and true, a fatal error occurs !! if this variable is not found in the parameter file + logical, optional, intent(out) :: set !< If present, this indicates whether this parameter + !! has been found and successfully set in the input files. ! Local variables character(len=CS%max_line_len) :: value_string(1) @@ -835,10 +860,13 @@ subroutine read_param_logical(CS, varname, value, fail_if_missing) elseif (present(fail_if_missing)) then ; if (fail_if_missing) then call MOM_error(FATAL, 'Unable to find variable '//trim(varname)//' in any input files.') endif ; endif + + if (present(set)) set = found + end subroutine read_param_logical !> This subroutine reads the value of a time_type model parameter from a parameter file. -subroutine read_param_time(CS, varname, value, timeunit, fail_if_missing, date_format) +subroutine read_param_time(CS, varname, value, timeunit, fail_if_missing, date_format, set) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: varname !< The case-sensitive name of the parameter to read @@ -850,6 +878,8 @@ subroutine read_param_time(CS, varname, value, timeunit, fail_if_missing, date_f logical, optional, intent(out) :: date_format !< If present, this indicates whether this !! parameter was read in a date format, so that it can !! later be logged in the same format. + logical, optional, intent(out) :: set !< If present, this indicates whether this parameter + !! has been found and successfully set in the input files. ! Local variables character(len=CS%max_line_len) :: value_string(1) @@ -891,6 +921,7 @@ subroutine read_param_time(CS, varname, value, timeunit, fail_if_missing, date_f read( value_string(1), *) real_time value = real_to_time(real_time*time_unit) endif + if (present(set)) set = .true. else if (present(fail_if_missing)) then ; if (fail_if_missing) then if (.not.found) then @@ -899,6 +930,7 @@ subroutine read_param_time(CS, varname, value, timeunit, fail_if_missing, date_f call MOM_error(FATAL, 'Variable '//trim(varname)//' found but not set in input files.') endif endif ; endif + if (present(set)) set = .false. endif return @@ -1704,7 +1736,7 @@ end function convert_date_to_string !! and logs it in documentation files. subroutine get_param_int(CS, modulename, varname, value, desc, units, & default, fail_if_missing, do_not_read, do_not_log, & - layoutParam, debuggingParam) + layoutParam, debuggingParam, old_name) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: modulename !< The name of the calling module @@ -1725,15 +1757,37 @@ subroutine get_param_int(CS, modulename, varname, value, desc, units, & !! logged in the layout parameter file logical, optional, intent(in) :: debuggingParam !< If present and true, this parameter is !! logged in the debugging parameter file + character(len=*), optional, intent(in) :: old_name !< A case-sensitive archaic name of the parameter + !! to read. Errors or warnings are issued if the old name + !! is being used. + ! Local variables logical :: do_read, do_log + logical :: new_name_used, old_name_used, same_value + integer :: new_name_value ! The value that is set when the standard name is used. do_read = .true. ; if (present(do_not_read)) do_read = .not.do_not_read do_log = .true. ; if (present(do_not_log)) do_log = .not.do_not_log if (do_read) then if (present(default)) value = default - call read_param_int(CS, varname, value, fail_if_missing) + + old_name_used = .false. + if (present(old_name)) then + new_name_value = value + call read_param_int(CS, old_name, value, set=old_name_used) + if (old_name_used) then + call read_param_int(CS, varname, new_name_value, set=new_name_used) + + ! Issue appropriate warnings or error messages. + same_value = (value == new_name_value) + call archaic_param_name_message(varname, old_name, new_name_used, same_value) + endif + endif + + if (.not.old_name_used) then ! Old name is either not present or not set. + call read_param_int(CS, varname, value, fail_if_missing) + endif endif if (do_log) then @@ -1747,7 +1801,7 @@ end subroutine get_param_int !! and logs them in documentation files. subroutine get_param_int_array(CS, modulename, varname, value, desc, units, & default, defaults, fail_if_missing, do_not_read, do_not_log, & - layoutParam, debuggingParam) + layoutParam, debuggingParam, old_name) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: modulename !< The name of the calling module @@ -1769,8 +1823,15 @@ subroutine get_param_int_array(CS, modulename, varname, value, desc, units, & !! logged in the layout parameter file logical, optional, intent(in) :: debuggingParam !< If present and true, this parameter is !! logged in the debugging parameter file + character(len=*), optional, intent(in) :: old_name !< A case-sensitive archaic name of the parameter + !! to read. Errors or warnings are issued if the old name + !! is being used. + ! Local variables logical :: do_read, do_log + logical :: new_name_used, old_name_used, same_value + integer :: new_name_value(size(value)) ! The values that are set when the old name is used. + integer :: m do_read = .true. ; if (present(do_not_read)) do_read = .not.do_not_read do_log = .true. ; if (present(do_not_log)) do_log = .not.do_not_log @@ -1785,7 +1846,24 @@ subroutine get_param_int_array(CS, modulename, varname, value, desc, units, & if (do_read) then if (present(default)) value(:) = default if (present(defaults)) value(:) = defaults(:) - call read_param_int_array(CS, varname, value, fail_if_missing) + + old_name_used = .false. + if (present(old_name)) then + new_name_value(:) = value(:) + call read_param_int_array(CS, old_name, value, set=old_name_used) + if (old_name_used) then + call read_param_int_array(CS, varname, new_name_value, set=new_name_used) + + ! Issue appropriate warnings or error messages. + same_value = .true. + do m=1,size(value) ; if (value(m) /= new_name_value(m)) same_value = .false. ; enddo + call archaic_param_name_message(varname, old_name, new_name_used, same_value) + endif + endif + + if (.not.old_name_used) then ! Old name is either not present or not set. + call read_param_int_array(CS, varname, value, fail_if_missing) + endif endif if (do_log) then @@ -1799,7 +1877,7 @@ end subroutine get_param_int_array !! and logs it in documentation files. subroutine get_param_real(CS, modulename, varname, value, desc, units, & default, fail_if_missing, do_not_read, do_not_log, & - debuggingParam, scale, unscaled) + debuggingParam, scale, unscaled, old_name) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: modulename !< The name of the calling module @@ -1822,15 +1900,37 @@ subroutine get_param_real(CS, modulename, varname, value, desc, units, & !! multiplied by before it is returned. real, optional, intent(out) :: unscaled !< The value of the parameter that would be !! returned without any multiplication by a scaling factor. + character(len=*), optional, intent(in) :: old_name !< A case-sensitive archaic name of the parameter + !! to read. Errors or warnings are issued if the old name + !! is being used. + ! Local variables logical :: do_read, do_log + logical :: new_name_used, old_name_used, same_value + real :: new_name_value ! The value that is set when the old name is used. do_read = .true. ; if (present(do_not_read)) do_read = .not.do_not_read do_log = .true. ; if (present(do_not_log)) do_log = .not.do_not_log if (do_read) then if (present(default)) value = default - call read_param_real(CS, varname, value, fail_if_missing) + + old_name_used = .false. + if (present(old_name)) then + new_name_value = value + call read_param_real(CS, old_name, value, set=old_name_used) + if (old_name_used) then + call read_param_real(CS, varname, new_name_value, set=new_name_used) + + ! Issue appropriate warnings or error messages. + same_value = (new_name_used .and. old_name_used .and. (value == new_name_value)) + call archaic_param_name_message(varname, old_name, new_name_used, same_value) + endif + endif + + if (.not.old_name_used) then ! Old name is either not present or not set. + call read_param_real(CS, varname, value, fail_if_missing) + endif endif if (do_log) then @@ -1847,7 +1947,7 @@ end subroutine get_param_real !! and logs them in documentation files. subroutine get_param_real_array(CS, modulename, varname, value, desc, units, & default, defaults, fail_if_missing, do_not_read, do_not_log, debuggingParam, & - scale, unscaled) + scale, unscaled, old_name) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: modulename !< The name of the calling module @@ -1871,8 +1971,15 @@ subroutine get_param_real_array(CS, modulename, varname, value, desc, units, & !! multiplied by before it is returned. real, dimension(:), optional, intent(out) :: unscaled !< The value of the parameter that would be !! returned without any multiplication by a scaling factor. + character(len=*), optional, intent(in) :: old_name !< A case-sensitive archaic name of the parameter + !! to read. Errors or warnings are issued if the old name + !! is being used. + ! Local variables logical :: do_read, do_log + logical :: new_name_used, old_name_used, same_value + real :: new_name_value(size(value)) ! The values that are set when the standard name is used. + integer :: m do_read = .true. ; if (present(do_not_read)) do_read = .not.do_not_read do_log = .true. ; if (present(do_not_log)) do_log = .not.do_not_log @@ -1887,7 +1994,24 @@ subroutine get_param_real_array(CS, modulename, varname, value, desc, units, & if (do_read) then if (present(default)) value(:) = default if (present(defaults)) value(:) = defaults(:) - call read_param_real_array(CS, varname, value, fail_if_missing) + + old_name_used = .false. + if (present(old_name)) then + new_name_value(:) = value(:) + call read_param_real_array(CS, old_name, value, set=old_name_used) + if (old_name_used) then + call read_param_real_array(CS, varname, new_name_value, set=new_name_used) + + ! Issue appropriate warnings or error messages. + same_value = .true. + do m=1,size(value) ; if (value(m) /= new_name_value(m)) same_value = .false. ; enddo + call archaic_param_name_message(varname, old_name, new_name_used, same_value) + endif + endif + + if (.not.old_name_used) then ! Old name is either not present or not set. + call read_param_real_array(CS, varname, value, fail_if_missing) + endif endif if (do_log) then @@ -1904,7 +2028,7 @@ end subroutine get_param_real_array !! and logs it in documentation files. subroutine get_param_char(CS, modulename, varname, value, desc, units, & default, fail_if_missing, do_not_read, do_not_log, & - layoutParam, debuggingParam) + layoutParam, debuggingParam, old_name) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: modulename !< The name of the calling module @@ -1925,15 +2049,37 @@ subroutine get_param_char(CS, modulename, varname, value, desc, units, & !! logged in the layout parameter file logical, optional, intent(in) :: debuggingParam !< If present and true, this parameter is !! logged in the debugging parameter file + character(len=*), optional, intent(in) :: old_name !< A case-sensitive archaic name of the parameter + !! to read. Errors or warnings are issued if the old name + !! is being used. + ! Local variables logical :: do_read, do_log + logical :: new_name_used, old_name_used, same_value + character(len=:), allocatable :: new_name_value ! The value that is set when the standard name is used. do_read = .true. ; if (present(do_not_read)) do_read = .not.do_not_read do_log = .true. ; if (present(do_not_log)) do_log = .not.do_not_log if (do_read) then if (present(default)) value = default - call read_param_char(CS, varname, value, fail_if_missing) + + old_name_used = .false. + if (present(old_name)) then + new_name_value = value + call read_param_char(CS, old_name, value, set=old_name_used) + if (old_name_used) then + call read_param_char(CS, varname, new_name_value, set=new_name_used) + + ! Issue appropriate warnings or error messages. + same_value = (trim(value) == trim(new_name_value)) + call archaic_param_name_message(varname, old_name, new_name_used, same_value) + endif + endif + + if (.not.old_name_used) then ! Old name is either not present or not set. + call read_param_char(CS, varname, value, fail_if_missing) + endif endif if (do_log) then @@ -1946,7 +2092,7 @@ end subroutine get_param_char !> This subroutine reads the values of an array of character string model parameters !! from a parameter file and logs them in documentation files. subroutine get_param_char_array(CS, modulename, varname, value, desc, units, & - default, fail_if_missing, do_not_read, do_not_log) + default, fail_if_missing, do_not_read, do_not_log, old_name) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: modulename !< The name of the calling module @@ -1963,18 +2109,40 @@ subroutine get_param_char_array(CS, modulename, varname, value, desc, units, & !! value for this parameter, although it might be logged. logical, optional, intent(in) :: do_not_log !< If present and true, do not log this !! parameter to the documentation files + character(len=*), optional, intent(in) :: old_name !< A case-sensitive archaic name of the parameter + !! to read. Errors or warnings are issued if the old name + !! is being used. ! Local variables logical :: do_read, do_log - integer :: i, len_tot, len_val + logical :: new_name_used, old_name_used, same_value + integer :: i, m, len_tot, len_val character(len=:), allocatable :: cat_val + character(len=:), allocatable :: new_name_value(:) ! The value that is set when the standard name is used. do_read = .true. ; if (present(do_not_read)) do_read = .not.do_not_read do_log = .true. ; if (present(do_not_log)) do_log = .not.do_not_log if (do_read) then if (present(default)) value(:) = default - call read_param_char_array(CS, varname, value, fail_if_missing) + + old_name_used = .false. + if (present(old_name)) then + new_name_value(:) = value(:) + call read_param_char_array(CS, old_name, value, set=old_name_used) + if (old_name_used) then + call read_param_char_array(CS, varname, new_name_value, set=new_name_used) + + ! Issue appropriate warnings or error messages. + same_value = .true. + do m=1,size(value) ; if (trim(value(m)) /= trim(new_name_value(m))) same_value = .false. ; enddo + call archaic_param_name_message(varname, old_name, new_name_used, same_value) + endif + endif + + if (.not.old_name_used) then ! Old name is either not present or not set. + call read_param_char_array(CS, varname, value, fail_if_missing) + endif endif if (do_log) then @@ -1996,7 +2164,7 @@ end subroutine get_param_char_array !! and logs it in documentation files. subroutine get_param_logical(CS, modulename, varname, value, desc, units, & default, fail_if_missing, do_not_read, do_not_log, & - layoutParam, debuggingParam) + layoutParam, debuggingParam, old_name) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: modulename !< The name of the calling module @@ -2017,15 +2185,37 @@ subroutine get_param_logical(CS, modulename, varname, value, desc, units, & !! logged in the layout parameter file logical, optional, intent(in) :: debuggingParam !< If present and true, this parameter is !! logged in the debugging parameter file + character(len=*), optional, intent(in) :: old_name !< A case-sensitive archaic name of the parameter + !! to read. Errors or warnings are issued if the old name + !! is being used. + ! Local variables logical :: do_read, do_log + logical :: new_name_used, old_name_used, same_value + logical :: new_name_value ! The value that is set when the standard name is used. do_read = .true. ; if (present(do_not_read)) do_read = .not.do_not_read do_log = .true. ; if (present(do_not_log)) do_log = .not.do_not_log if (do_read) then if (present(default)) value = default - call read_param_logical(CS, varname, value, fail_if_missing) + + old_name_used = .false. + if (present(old_name)) then + new_name_value = value + call read_param_logical(CS, old_name, value, set=old_name_used) + if (old_name_used) then + call read_param_logical(CS, varname, new_name_value, set=new_name_used) + + ! Issue appropriate warnings or error messages. + same_value = (value .eqv. new_name_value) + call archaic_param_name_message(varname, old_name, new_name_used, same_value) + endif + endif + + if (.not.old_name_used) then ! Old name is either not present or not set. + call read_param_logical(CS, varname, value, fail_if_missing) + endif endif if (do_log) then @@ -2040,7 +2230,7 @@ end subroutine get_param_logical subroutine get_param_time(CS, modulename, varname, value, desc, units, & default, fail_if_missing, do_not_read, do_not_log, & timeunit, layoutParam, debuggingParam, & - log_as_date) + log_as_date, old_name) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: modulename !< The name of the calling module @@ -2065,8 +2255,14 @@ subroutine get_param_time(CS, modulename, varname, value, desc, units, & !! logged in the debugging parameter file logical, optional, intent(in) :: log_as_date !< If true, log the time_type in date !! format. The default is false. + character(len=*), optional, intent(in) :: old_name !< A case-sensitive archaic name of the parameter + !! to read. Errors or warnings are issued if the old name + !! is being used. + ! Local variables logical :: do_read, do_log, log_date + logical :: new_name_used, old_name_used, same_value + type(time_type) :: new_name_value ! The value that is set when the standard name is used. do_read = .true. ; if (present(do_not_read)) do_read = .not.do_not_read do_log = .true. ; if (present(do_not_log)) do_log = .not.do_not_log @@ -2074,7 +2270,23 @@ subroutine get_param_time(CS, modulename, varname, value, desc, units, & if (do_read) then if (present(default)) value = default - call read_param_time(CS, varname, value, timeunit, fail_if_missing, date_format=log_date) + + old_name_used = .false. + if (present(old_name)) then + new_name_value = value + call read_param_time(CS, old_name, value, timeunit, date_format=log_date, set=old_name_used) + if (old_name_used) then + call read_param_time(CS, varname, new_name_value, timeunit, date_format=log_date, set=new_name_used) + + ! Issue appropriate warnings or error messages. + same_value = (value == new_name_value) + call archaic_param_name_message(varname, old_name, new_name_used, same_value) + endif + endif + + if (.not.old_name_used) then ! Old name is either not present or not set. + call read_param_time(CS, varname, value, timeunit, fail_if_missing, date_format=log_date) + endif endif if (do_log) then @@ -2086,6 +2298,28 @@ subroutine get_param_time(CS, modulename, varname, value, desc, units, & end subroutine get_param_time +!> Issue error messages or warnings about the use of an archaic parameter name. +subroutine archaic_param_name_message(varname, old_name, new_name_used, same_value) + character(len=*), intent(in) :: varname !< The case-sensitive name of the parameter to read + character(len=*), intent(in) :: old_name !< The case-sensitive archaic name of the parameter + logical, intent(in) :: new_name_used !< True if varname is used in the parameter file. + logical, intent(in) :: same_value !< True if varname and old_name give the same values. + + if (new_name_used .and. same_value) then + call MOM_error(WARNING, "The runtime parameter "//trim(varname)//& + " is also being set consistently via its older name of "//trim(old_name)//& + ". Please migrate to only using "//trim(varname)//".") + elseif (new_name_used .and. .not.same_value) then + call MOM_error(FATAL, "The runtime parameter "//trim(varname)//& + " is also being set inconsistently via its older name of "//trim(old_name)//& + ". Only use "//trim(varname)//".") + else + call MOM_error(WARNING, "The runtime parameter "//trim(varname)//& + " is being set via its soon to be obsolete name of "//trim(old_name)//& + ". Please migrate to using "//trim(varname)//".") + endif +end subroutine archaic_param_name_message + ! ----------------------------------------------------------------------------- !> Resets the parameter block name to blank From 400c982b003d49e156951ffc2d97f364352322cd Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 4 Dec 2024 05:39:07 -0500 Subject: [PATCH 12/17] +Add grid_unit_to_L to the ocean_grid_type Add the new element grid_unit_to_L to the ocean_grid_type and the dyn_horgrid_type, which can be used to convert the units of the geoLat and geoLon variables to rescaled horizontal distance units ([L ~> m]) when they are Cartesian coordinates. When Cartesian coordinates are not in use, G%grid_unit_to_L is set to 0. This new element of the grid type is used to test for inconsistent grids or to eliminate rescaling variables in set_rotation_beta_plane(), initialize_velocity_circular(), DOME_initialize_topography(), DOME_initialize_sponges(), DOME_set_OBC_data(), ISOMIP_initialize_topography(), idealized_hurricane_wind_forcing(), Kelvin_set_OBC_data(), Rossby_front_initialize_velocity(), soliton_initialize_thickness(), and soliton_initialize_velocity(). These are the instances where this new variable could be used and bitwise identical answers are recovered. There are a few other places where they should be used, but where answers would change, and these will be deferred to a subsequent commit. All answers are bitwise identical, but there are new elements in two transparent and widely used types. --- src/core/MOM_grid.F90 | 4 ++ src/core/MOM_transcribe_grid.F90 | 2 + src/framework/MOM_dyn_horgrid.F90 | 5 ++ src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 | 3 + src/initialization/MOM_grid_initialize.F90 | 9 ++- .../MOM_shared_initialization.F90 | 9 +-- .../MOM_state_initialization.F90 | 7 ++- src/user/DOME_initialization.F90 | 60 +++++++++++++------ src/user/ISOMIP_initialization.F90 | 18 +++--- src/user/Idealized_Hurricane.F90 | 13 ++-- src/user/Kelvin_initialization.F90 | 21 ++++--- src/user/Rossby_front_2d_initialization.F90 | 8 +-- src/user/soliton_initialization.F90 | 8 ++- 13 files changed, 110 insertions(+), 57 deletions(-) diff --git a/src/core/MOM_grid.F90 b/src/core/MOM_grid.F90 index 6fb8426395..5e5b307103 100644 --- a/src/core/MOM_grid.F90 +++ b/src/core/MOM_grid.F90 @@ -190,6 +190,10 @@ module MOM_grid ! These parameters are run-time parameters that are used during some ! initialization routines (but not all) + real :: grid_unit_to_L !< A factor that converts a the geoLat and geoLon variables and related + !! variables like len_lat and len_lon into rescaled horizontal distance + !! units on a Cartesian grid, in [L km ~> 1000] or [L m-1 ~> 1] or + !! is 0 for a non-Cartesian grid. real :: south_lat !< The latitude (or y-coordinate) of the first v-line [degrees_N] or [km] or [m] real :: west_lon !< The longitude (or x-coordinate) of the first u-line [degrees_E] or [km] or [m] real :: len_lat !< The latitudinal (or y-coord) extent of physical domain [degrees_N] or [km] or [m] diff --git a/src/core/MOM_transcribe_grid.F90 b/src/core/MOM_transcribe_grid.F90 index f8ae58d9e1..4c1b87f37e 100644 --- a/src/core/MOM_transcribe_grid.F90 +++ b/src/core/MOM_transcribe_grid.F90 @@ -135,6 +135,7 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG, US) ! Copy various scalar variables and strings. oG%x_axis_units = dG%x_axis_units ; oG%y_axis_units = dG%y_axis_units oG%x_ax_unit_short = dG%x_ax_unit_short ; oG%y_ax_unit_short = dG%y_ax_unit_short + oG%grid_unit_to_L = dG%grid_unit_to_L oG%areaT_global = dG%areaT_global ; oG%IareaT_global = dG%IareaT_global oG%south_lat = dG%south_lat ; oG%west_lon = dG%west_lon oG%len_lat = dG%len_lat ; oG%len_lon = dG%len_lon @@ -296,6 +297,7 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG, US) ! Copy various scalar variables and strings. dG%x_axis_units = oG%x_axis_units ; dG%y_axis_units = oG%y_axis_units dG%x_ax_unit_short = oG%x_ax_unit_short ; dG%y_ax_unit_short = oG%y_ax_unit_short + dG%grid_unit_to_L = oG%grid_unit_to_L dG%areaT_global = oG%areaT_global ; dG%IareaT_global = oG%IareaT_global dG%south_lat = oG%south_lat ; dG%west_lon = oG%west_lon dG%len_lat = oG%len_lat ; dG%len_lon = oG%len_lon diff --git a/src/framework/MOM_dyn_horgrid.F90 b/src/framework/MOM_dyn_horgrid.F90 index 987d5bf502..440ea351b9 100644 --- a/src/framework/MOM_dyn_horgrid.F90 +++ b/src/framework/MOM_dyn_horgrid.F90 @@ -181,6 +181,10 @@ module MOM_dyn_horgrid ! These parameters are run-time parameters that are used during some ! initialization routines (but not all) + real :: grid_unit_to_L !< A factor that converts a the geoLat and geoLon variables and related + !! variables like len_lat and len_lon into rescaled horizontal distance + !! units on a Cartesian grid, in [L km ~> 1000] or [L m-1 ~> 1] or + !! is 0 for a non-Cartesian grid. real :: south_lat !< The latitude (or y-coordinate) of the first v-line [degrees_N] or [km] or [m] real :: west_lon !< The longitude (or x-coordinate) of the first u-line [degrees_E] or [km] or [m] real :: len_lat !< The latitudinal (or y-coord) extent of physical domain [degrees_N] or [km] or [m] @@ -400,6 +404,7 @@ subroutine rotate_dyn_horgrid(G_in, G, US, turns) G%len_lon = G_in%len_lon ! Rotation-invariant fields + G%grid_unit_to_L = G_in%grid_unit_to_L G%areaT_global = G_in%areaT_global G%IareaT_global = G_in%IareaT_global G%Rad_Earth = G_in%Rad_Earth diff --git a/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 b/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 index 5ecfb9e788..fe54dd6533 100644 --- a/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 +++ b/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 @@ -132,6 +132,7 @@ subroutine set_IS_axes_info(G, param_file, diag_cs, axes_set_name) G%x_axis_units = "degrees_E" ; G%y_axis_units = "degrees_N" G%x_ax_unit_short = "degrees_E" ; G%y_ax_unit_short = "degrees_N" + G%grid_unit_to_L = 0.0 if (index(lowercase(trim(grid_config)),"cartesian") > 0) then ! This is a cartesian grid, and may have different axis units. @@ -145,9 +146,11 @@ subroutine set_IS_axes_info(G, param_file, diag_cs, axes_set_name) if (units_temp(1:1) == 'k') then G%x_axis_units = "kilometers" ; G%y_axis_units = "kilometers" G%x_ax_unit_short = "km" ; G%y_ax_unit_short = "km" + G%grid_unit_to_L = 1000.0*diag_cs%US%m_to_L elseif (units_temp(1:1) == 'm') then G%x_axis_units = "meters" ; G%y_axis_units = "meters" G%x_ax_unit_short = "m" ; G%y_ax_unit_short = "m" + G%grid_unit_to_L = diag_cs%US%m_to_L endif call log_param(param_file, mdl, "explicit AXIS_UNITS", G%x_axis_units) else diff --git a/src/initialization/MOM_grid_initialize.F90 b/src/initialization/MOM_grid_initialize.F90 index ef78a896c3..429ce24d79 100644 --- a/src/initialization/MOM_grid_initialize.F90 +++ b/src/initialization/MOM_grid_initialize.F90 @@ -84,7 +84,7 @@ subroutine set_grid_metrics(G, param_file, US) ! These are defaults that may be changed in the next select block. G%x_axis_units = "degrees_east" ; G%y_axis_units = "degrees_north" G%x_ax_unit_short = "degrees_E" ; G%y_ax_unit_short = "degrees_N" - + G%grid_unit_to_L = 0.0 G%Rad_Earth_L = -1.0*US%m_to_L ; G%len_lat = 0.0 ; G%len_lon = 0.0 select case (trim(config)) case ("mosaic"); call set_grid_metrics_from_mosaic(G, param_file, US) @@ -251,6 +251,11 @@ subroutine set_grid_metrics_from_mosaic(G, param_file, US) G%geoLatCv(i,J) = tmpZ(i2-1,j2) enddo ; enddo + ! This routine could be modified to support the use of a mosaic using Cartesian grid coordinates, + ! in which case the values of G%x_axis_units, G%y_axis_units and G%grid_unit_to_L would need to be + ! reset appropriately here, but this option has not yet been implemented, and the grid coordinates + ! are assumed to be degrees of longitude and latitude. + ! Read DX,DY from the supergrid tmpU(:,:) = 0. ; tmpV(:,:) = 0. call MOM_read_data(filename, 'dx', tmpV, SGdom, position=NORTH_FACE, scale=US%m_to_L) @@ -440,9 +445,11 @@ subroutine set_grid_metrics_cartesian(G, param_file, US) enddo if (units_temp(1:1) == 'k') then ! Axes are measured in km. + G%grid_unit_to_L = 1000.0*US%m_to_L dx_everywhere = 1000.0*US%m_to_L * G%len_lon / (REAL(niglobal)) dy_everywhere = 1000.0*US%m_to_L * G%len_lat / (REAL(njglobal)) elseif (units_temp(1:1) == 'm') then ! Axes are measured in m. + G%grid_unit_to_L = US%m_to_L dx_everywhere = US%m_to_L*G%len_lon / (REAL(niglobal)) dy_everywhere = US%m_to_L*G%len_lat / (REAL(njglobal)) else ! Axes are measured in degrees of latitude and longitude. diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90 index 0a257b6ca1..1f7a519d9e 100644 --- a/src/initialization/MOM_shared_initialization.F90 +++ b/src/initialization/MOM_shared_initialization.F90 @@ -485,7 +485,6 @@ subroutine set_rotation_beta_plane(f, G, param_file, US) real :: f_0 ! The reference value of the Coriolis parameter [T-1 ~> s-1] real :: beta ! The meridional gradient of the Coriolis parameter [T-1 L-1 ~> s-1 m-1] real :: beta_lat_ref ! The reference latitude for the beta plane [degrees_N] or [km] or [m] - real :: Rad_Earth_L ! The radius of the planet in rescaled units [L ~> m] real :: y_scl ! A scaling factor from the units of latitude [L lat-1 ~> m lat-1] real :: PI ! The ratio of the circumference of a circle to its diameter [nondim] character(len=40) :: mdl = "set_rotation_beta_plane" ! This subroutine's name. @@ -503,18 +502,16 @@ subroutine set_rotation_beta_plane(f, G, param_file, US) call get_param(param_file, mdl, "AXIS_UNITS", axis_units, default="degrees") PI = 4.0*atan(1.0) + y_scl = G%grid_unit_to_L + if (G%grid_unit_to_L <= 0.0) y_scl = PI * G%Rad_Earth_L / 180. + select case (axis_units(1:1)) case ("d") - call get_param(param_file, mdl, "RAD_EARTH", Rad_Earth_L, & - "The radius of the Earth.", units="m", default=6.378e6, scale=US%m_to_L) beta_lat_ref_units = "degrees" - y_scl = PI * Rad_Earth_L / 180. case ("k") beta_lat_ref_units = "kilometers" - y_scl = 1.0e3 * US%m_to_L case ("m") beta_lat_ref_units = "meters" - y_scl = 1.0 * US%m_to_L case default ; call MOM_error(FATAL, & " set_rotation_beta_plane: unknown AXIS_UNITS = "//trim(axis_units)) end select diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 7909bddb80..0f5d06aa4e 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -1636,7 +1636,10 @@ subroutine initialize_velocity_circular(u, v, G, GV, US, param_file, just_read) if (just_read) return ! All run-time parameters have been read, so return. - dpi=acos(0.0)*2.0 ! pi + if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, "MOM_state_initialization.F90: "//& + "initialize_velocity_circular() is only set to work with Cartesian axis units.") + + dpi = acos(0.0)*2.0 ! pi do k=1,nz ; do j=js,je ; do I=Isq,Ieq psi1 = my_psi(I,j) @@ -1663,7 +1666,7 @@ real function my_psi(ig,jg) r = sqrt( (x**2) + (y**2) ) ! Circular stream function is a function of radius only r = min(1.0, r) ! Flatten stream function in corners of box my_psi = 0.5*(1.0 - cos(dpi*r)) - my_psi = my_psi * (circular_max_u * G%US%m_to_L*G%len_lon*1e3 / dpi) ! len_lon is in km + my_psi = my_psi * (circular_max_u * G%len_lon * G%grid_unit_to_L / dpi) ! len_lon is in km end function my_psi end subroutine initialize_velocity_circular diff --git a/src/user/DOME_initialization.F90 b/src/user/DOME_initialization.F90 index 858ca32f93..e608dbd1c2 100644 --- a/src/user/DOME_initialization.F90 +++ b/src/user/DOME_initialization.F90 @@ -49,10 +49,11 @@ subroutine DOME_initialize_topography(D, G, param_file, max_depth, US) real :: min_depth ! The minimum ocean depth [Z ~> m] real :: shelf_depth ! The ocean depth on the shelf in the DOME configuration [Z ~> m] real :: slope ! The bottom slope in the DOME configuration [Z L-1 ~> nondim] - real :: shelf_edge_lat ! The latitude of the edge of the topographic shelf [km] - real :: inflow_lon ! The edge longitude of the DOME inflow [km] - real :: inflow_width ! The longitudinal width of the DOME inflow channel [km] - real :: km_to_L ! The conversion factor from the units of latitude to L [L km-1 ~> 1e3] + real :: shelf_edge_lat ! The latitude of the edge of the topographic shelf in the same units as geolat, often [km] + real :: inflow_lon ! The edge longitude of the DOME inflow in the same units as geolon, often [km] + real :: inflow_width ! The longitudinal width of the DOME inflow channel in the same units as geolat, often [km] + real :: km_to_grid_unit ! The conversion factor from km to the units of latitude often 1 [nondim], + ! but this could be 1000 [m km-1] ! This include declares and sets the variable "version". # include "version_variable.h" character(len=40) :: mdl = "DOME_initialize_topography" ! This subroutine's name. @@ -60,7 +61,16 @@ subroutine DOME_initialize_topography(D, G, param_file, max_depth, US) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed - km_to_L = 1.0e3*US%m_to_L + if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, "DOME_initialization: "//& + "DOME_initialize_topography is only set to work with Cartesian axis units.") + if (abs(G%grid_unit_to_L*US%L_to_m - 1000.0) < 1.0e-3) then ! The grid latitudes are in km. + km_to_grid_unit = 1.0 + elseif (abs(G%grid_unit_to_L*US%L_to_m - 1.0) < 1.0e-6) then ! The grid latitudes are in m. + km_to_grid_unit = 1000.0 + else + call MOM_error(FATAL, "DOME_initialization: "//& + "DOME_initialize_topography is not recognizing the value of G%grid_unit_to_L.") + endif call MOM_mesg(" DOME_initialization.F90, DOME_initialize_topography: setting topography", 5) @@ -75,15 +85,16 @@ subroutine DOME_initialize_topography(D, G, param_file, max_depth, US) default=600.0, units="m", scale=US%m_to_Z) call get_param(param_file, mdl, "DOME_SHELF_EDGE_LAT", shelf_edge_lat, & "The latitude of the shelf edge in the DOME configuration.", & - default=600.0, units="km") + default=600.0, units="km", scale=km_to_grid_unit) call get_param(param_file, mdl, "DOME_INFLOW_LON", inflow_lon, & - "The edge longitude of the DOME inflow.", units="km", default=1000.0) + "The edge longitude of the DOME inflow.", units="km", default=1000.0, scale=km_to_grid_unit) call get_param(param_file, mdl, "DOME_INFLOW_WIDTH", inflow_width, & - "The longitudinal width of the DOME inflow channel.", units="km", default=100.0) + "The longitudinal width of the DOME inflow channel.", & + units="km", default=100.0, scale=km_to_grid_unit) do j=js,je ; do i=is,ie if (G%geoLatT(i,j) < shelf_edge_lat) then - D(i,j) = min(shelf_depth - slope * (G%geoLatT(i,j)-shelf_edge_lat)*km_to_L, max_depth) + D(i,j) = min(shelf_depth - slope * (G%geoLatT(i,j)-shelf_edge_lat)*G%grid_unit_to_L, max_depth) else if ((G%geoLonT(i,j) > inflow_lon) .AND. (G%geoLonT(i,j) < inflow_lon+inflow_width)) then D(i,j) = shelf_depth @@ -177,7 +188,6 @@ subroutine DOME_initialize_sponges(G, GV, US, tv, depth_tot, PF, CSp) real :: min_depth ! The minimum depth at which to apply damping [Z ~> m] real :: damp_W, damp_E ! Damping rates in the western and eastern sponges [T-1 ~> s-1] real :: peak_damping ! The maximum sponge damping rates as the edges [T-1 ~> s-1] - real :: km_to_L ! The conversion factor from the units of longitude to L [L km-1 ~> 1e3] real :: edge_dist ! The distance to an edge [L ~> m] real :: sponge_width ! The width of the sponges [L ~> m] character(len=40) :: mdl = "DOME_initialize_sponges" ! This subroutine's name. @@ -186,7 +196,8 @@ subroutine DOME_initialize_sponges(G, GV, US, tv, depth_tot, PF, CSp) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed - km_to_L = 1.0e3*US%m_to_L + if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, "DOME_initialization: "//& + "DOME_initialize_sponges is only set to work with Cartesian axis units.") ! Set up sponges for the DOME configuration call get_param(PF, mdl, "MINIMUM_DEPTH", min_depth, & @@ -196,7 +207,7 @@ subroutine DOME_initialize_sponges(G, GV, US, tv, depth_tot, PF, CSp) default=10.0, units="day-1", scale=1.0/(86400.0*US%s_to_T)) call get_param(PF, mdl, "DOME_SPONGE_WIDTH", sponge_width, & "The width of the the DOME sponges.", & - default=200.0, units="km", scale=km_to_L) + default=200.0, units="km", scale=1.0e3*US%m_to_L) ! Here the inverse damping time [T-1 ~> s-1], is set. Set Idamp to 0 wherever ! there is no sponge, and the subroutines that are called will automatically @@ -204,7 +215,7 @@ subroutine DOME_initialize_sponges(G, GV, US, tv, depth_tot, PF, CSp) Idamp(:,:) = 0.0 do j=js,je ; do i=is,ie ; if (depth_tot(i,j) > min_depth) then - edge_dist = (G%geoLonT(i,j) - G%west_lon) * km_to_L + edge_dist = (G%geoLonT(i,j) - G%west_lon) * G%grid_unit_to_L if (edge_dist < 0.5*sponge_width) then damp_W = peak_damping elseif (edge_dist < sponge_width) then @@ -213,7 +224,7 @@ subroutine DOME_initialize_sponges(G, GV, US, tv, depth_tot, PF, CSp) damp_W = 0.0 endif - edge_dist = ((G%len_lon + G%west_lon) - G%geoLonT(i,j)) * km_to_L + edge_dist = ((G%len_lon + G%west_lon) - G%geoLonT(i,j)) * G%grid_unit_to_L if (edge_dist < 0.5*sponge_width) then damp_E = peak_damping elseif (edge_dist < sponge_width) then @@ -328,10 +339,12 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, US, PF, tr_Reg) ! properties [T-1 ~> s-1] real :: g_prime_tot ! The reduced gravity across all layers [L2 Z-1 T-2 ~> m s-2] real :: Def_Rad ! The deformation radius, based on fluid of thickness D_edge [L ~> m] - real :: inflow_lon ! The edge longitude of the DOME inflow [km] + real :: inflow_lon ! The edge longitude of the DOME inflow in the same units as geolon, often [km] real :: I_Def_Rad ! The inverse of the deformation radius in the same units as G%geoLon [km-1] real :: Ri_trans ! The shear Richardson number in the transition ! region of the specified shear profile [nondim] + real :: km_to_grid_unit ! The conversion factor from km to the units of latitude often 1 [nondim], + ! but this could be 1000 [m km-1] character(len=32) :: name ! The name of a tracer field. character(len=40) :: mdl = "DOME_set_OBC_data" ! This subroutine's name. integer :: i, j, k, itt, is, ie, js, je, isd, ied, jsd, jed, m, nz, ntherm, ntr_id @@ -343,6 +356,17 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, US, PF, tr_Reg) isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB + if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, "DOME_initialization: "//& + "DOME_initialize_topography is only set to work with Cartesian axis units.") + if (abs(G%grid_unit_to_L*US%L_to_m - 1000.0) < 1.0e-3) then ! The grid latitudes are in km. + km_to_grid_unit = 1.0 + elseif (abs(G%grid_unit_to_L*US%L_to_m - 1.0) < 1.0e-6) then ! The grid latitudes are in m. + km_to_grid_unit = 1000.0 + else + call MOM_error(FATAL, "DOME_initialization: "//& + "DOME_initialize_topography is not recognizing the value of G%grid_unit_to_L.") + endif + call get_param(PF, mdl, "DOME_INFLOW_THICKNESS", D_edge, & "The thickness of the dense DOME inflow at the inner edge.", & default=300.0, units="m", scale=US%m_to_Z) @@ -362,7 +386,7 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, US, PF, tr_Reg) "The value of the Coriolis parameter that is used to determine the DOME "//& "inflow properties.", units="s-1", default=f_0*US%s_to_T, scale=US%T_to_s) call get_param(PF, mdl, "DOME_INFLOW_LON", inflow_lon, & - "The edge longitude of the DOME inflow.", units="km", default=1000.0) + "The edge longitude of the DOME inflow.", units="km", default=1000.0, scale=km_to_grid_unit) if (associated(tv%S) .or. associated(tv%T)) then call get_param(PF, mdl, "S_REF", S_ref, & units="ppt", default=35.0, scale=US%ppt_to_S, do_not_log=.true.) @@ -383,7 +407,9 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, US, PF, tr_Reg) tr_0 = (-D_edge*sqrt(D_edge*g_prime_tot)*0.5*Def_Rad) * (Rlay_Ref + 0.5*Rlay_range) * GV%RZ_to_H endif - I_Def_Rad = 1.0 / (1.0e-3*US%L_to_m*Def_Rad) + I_Def_Rad = 1.0 / ((1.0e-3*US%L_to_m*km_to_grid_unit) * Def_Rad) + ! This is mathematically equivalent to + ! I_Def_Rad = G%grid_unit_to_L / Def_Rad if (OBC%number_of_segments /= 1) then call MOM_error(WARNING, 'Error in DOME OBC segment setup', .true.) diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index d03a07e313..4bf7931856 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -59,7 +59,6 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth, US) real :: ly ! domain width (across ice flow) [L ~> m] real :: bx, by ! The x- and y- contributions to the bathymetric profiles at a point [Z ~> m] real :: xtil ! x-positon normalized by the characteristic along-flow length scale [nondim] - real :: km_to_L ! The conversion factor from the axis units to L [L km-1 ~> 1e3] logical :: is_2D ! If true, use a 2D setup ! This include declares and sets the variable "version". # include "version_variable.h" @@ -93,7 +92,8 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth, US) "Characteristic width of the side walls of the channel in the ISOMIP configuration.", & units="m", default=4.0e3, scale=US%m_to_L) - km_to_L = 1.0e3*US%m_to_L + if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, "ISOMIP_initialization.F90: " //& + "ISOMIP_initialize_topography is only set to work with Cartesian axis units.") ! The following variables should be transformed into runtime parameters. b0 = -150.0*US%m_to_Z ; b2 = -728.8*US%m_to_Z ; b4 = 343.91*US%m_to_Z ; b6 = -50.57*US%m_to_Z @@ -101,8 +101,8 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth, US) if (is_2D) then do j=js,je ; do i=is,ie ! For the 2D setup take a slice through the middle of the domain - xtil = G%geoLonT(i,j)*km_to_L / xbar - !xtil = 450.*km_to_L / xbar + xtil = G%geoLonT(i,j)*G%grid_unit_to_L / xbar + ! xtil = 450.0e3*US%m_to_L / xbar bx = b0 + b2*xtil**2 + b4*xtil**4 + b6*xtil**6 by = 2.0 * dc / (1.0 + exp(2.0*wc / fc)) @@ -117,17 +117,17 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth, US) ! 3D setup ! ===== TEST ===== !if (G%geoLonT(i,j)<500.) then - ! xtil = 500.*km_to_L / xbar + ! xtil = 500.0e3*US%m_to_L / xbar !else - ! xtil = G%geoLonT(i,j)*km_to_L / xbar + ! xtil = G%geoLonT(i,j)*G%grid_unit_to_L / xbar !endif ! ===== TEST ===== - xtil = G%geoLonT(i,j)*km_to_L / xbar + xtil = G%geoLonT(i,j)*G%grid_unit_to_L / xbar bx = b0 + b2*xtil**2 + b4*xtil**4 + b6*xtil**6 - by = (dc / (1.0 + exp(-2.*(G%geoLatT(i,j)*km_to_L - 0.5*ly - wc) / fc))) + & - (dc / (1.0 + exp(2.*(G%geoLatT(i,j)*km_to_L - 0.5*ly + wc) / fc))) + by = (dc / (1.0 + exp(-2.*(G%geoLatT(i,j)*G%grid_unit_to_L - 0.5*ly - wc) / fc))) + & + (dc / (1.0 + exp(2.*(G%geoLatT(i,j)*G%grid_unit_to_L - 0.5*ly + wc) / fc))) D(i,j) = -max(bx+by, -bmax) if (D(i,j) > max_depth) D(i,j) = max_depth diff --git a/src/user/Idealized_Hurricane.F90 b/src/user/Idealized_Hurricane.F90 index aca19e86d0..99cc7b88c5 100644 --- a/src/user/Idealized_Hurricane.F90 +++ b/src/user/Idealized_Hurricane.F90 @@ -349,7 +349,6 @@ subroutine idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, CS) real :: fbench !< The benchmark 'f' value [T-1 ~> s-1] real :: fbench_fac !< A factor that is set to 0 to use the !! benchmark 'f' value [nondim] - real :: km_to_L !< The conversion factor from the units of latitude to L [L km-1 ~> 1e3] real :: rel_tau_fac !< A factor that is set to 0 to disable !! current relative stress calculation [nondim] @@ -359,7 +358,8 @@ subroutine idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, CS) isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB - km_to_L = 1.0e3*US%m_to_L + if ((G%grid_unit_to_L <= 0.) .and. (.not.CS%SCM_mode)) call MOM_error(FATAL, "Idealized_Hurricane.F90: " //& + "idealized_hurricane_wind_forcing is only set to work with Cartesian axis units.") ! Allocate the forcing arrays, if necessary. call allocate_mech_forcing(G, forces, stress=.true., ustar=.true., tau_mag=.true.) @@ -376,7 +376,6 @@ subroutine idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, CS) YC = CS%Hurr_cen_Y0 + (time_type_to_real(day)*US%s_to_T * CS%hurr_translation_spd * & sin(CS%hurr_translation_dir)) - if (CS%BR_Bench) then ! f reset to value used in generated wind for benchmark test fbench = CS%f_column @@ -403,8 +402,8 @@ subroutine idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, CS) YY = CS%dy_from_center - YC XX = -XC else - YY = G%geoLatCu(I,j)*km_to_L - YC - XX = G%geoLonCu(I,j)*km_to_L - XC + YY = G%geoLatCu(I,j) * G%grid_unit_to_L - YC + XX = G%geoLonCu(I,j) * G%grid_unit_to_L - XC endif call idealized_hurricane_wind_profile(CS, US, f_local, YY, XX, Uocn, Vocn, TX, TY) forces%taux(I,j) = G%mask2dCu(I,j) * TX @@ -427,8 +426,8 @@ subroutine idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, CS) YY = CS%dy_from_center - YC XX = -XC else - YY = G%geoLatCv(i,J)*km_to_L - YC - XX = G%geoLonCv(i,J)*km_to_L - XC + YY = G%geoLatCv(i,J) * G%grid_unit_to_L - YC + XX = G%geoLonCv(i,J) * G%grid_unit_to_L - XC endif call idealized_hurricane_wind_profile(CS, US, f_local, YY, XX, Uocn, Vocn, TX, TY) forces%tauy(i,J) = G%mask2dCv(i,J) * TY diff --git a/src/user/Kelvin_initialization.F90 b/src/user/Kelvin_initialization.F90 index 1fc8a2f564..118d76ac93 100644 --- a/src/user/Kelvin_initialization.F90 +++ b/src/user/Kelvin_initialization.F90 @@ -215,10 +215,11 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time) if (.not.associated(OBC)) call MOM_error(FATAL, 'Kelvin_initialization.F90: '// & 'Kelvin_set_OBC_data() was called but OBC type was not initialized!') + if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, 'Kelvin_initialization.F90: '// & + "Kelvin_set_OBC_data() is only set to work with Cartesian axis units.") time_sec = US%s_to_T*time_type_to_real(Time) PI = 4.0*atan(1.0) - km_to_L_scale = 1000.0*US%m_to_L do j=jsd,jed ; do i=isd,ied depth_tot(i,j) = 0.0 @@ -237,6 +238,8 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time) lambda = PI * CS%mode * CS%F_0 / (CS%H0 * N0) ! Two wavelengths in domain omega = (4.0 * CS%H0 * N0) / (CS%mode * US%m_to_L*G%len_lon) + !### There is a bug here when len_lon is in km. This should be + ! omega = (4.0 * CS%H0 * N0) / (CS%mode * G%grid_unit_to_L*G%len_lon) endif sina = sin(CS%coast_angle) @@ -257,8 +260,8 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time) jsd = segment%HI%jsd ; jed = segment%HI%jed JsdB = segment%HI%JsdB ; JedB = segment%HI%JedB do j=jsd,jed ; do I=IsdB,IedB - x1 = km_to_L_scale * G%geoLonCu(I,j) - y1 = km_to_L_scale * G%geoLatCu(I,j) + x1 = G%grid_unit_to_L * G%geoLonCu(I,j) + y1 = G%grid_unit_to_L * G%geoLatCu(I,j) x = (x1 - CS%coast_offset1) * cosa + y1 * sina y = -(x1 - CS%coast_offset1) * sina + y1 * cosa if (CS%mode == 0) then @@ -299,8 +302,8 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time) enddo ; enddo if (allocated(segment%tangential_vel)) then do J=JsdB+1,JedB-1 ; do I=IsdB,IedB - x1 = km_to_L_scale * G%geoLonBu(I,J) - y1 = km_to_L_scale * G%geoLatBu(I,J) + x1 = G%grid_unit_to_L * G%geoLonBu(I,J) + y1 = G%grid_unit_to_L * G%geoLatBu(I,J) x = (x1 - CS%coast_offset1) * cosa + y1 * sina y = - (x1 - CS%coast_offset1) * sina + y1 * cosa cff = sqrt(GV%g_Earth * depth_tot(i+1,j) ) @@ -316,8 +319,8 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time) isd = segment%HI%isd ; ied = segment%HI%ied JsdB = segment%HI%JsdB ; JedB = segment%HI%JedB do J=JsdB,JedB ; do i=isd,ied - x1 = km_to_L_scale * G%geoLonCv(i,J) - y1 = km_to_L_scale * G%geoLatCv(i,J) + x1 = G%grid_unit_to_L * G%geoLonCv(i,J) + y1 = G%grid_unit_to_L * G%geoLatCv(i,J) x = (x1 - CS%coast_offset1) * cosa + y1 * sina y = - (x1 - CS%coast_offset1) * sina + y1 * cosa if (CS%mode == 0) then @@ -355,8 +358,8 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time) enddo ; enddo if (allocated(segment%tangential_vel)) then do J=JsdB,JedB ; do I=IsdB+1,IedB-1 - x1 = km_to_L_scale * G%geoLonBu(I,J) - y1 = km_to_L_scale * G%geoLatBu(I,J) + x1 = G%grid_unit_to_L * G%geoLonBu(I,J) + y1 = G%grid_unit_to_L * G%geoLatBu(I,J) x = (x1 - CS%coast_offset1) * cosa + y1 * sina y = - (x1 - CS%coast_offset1) * sina + y1 * cosa cff = sqrt(GV%g_Earth * depth_tot(i,j+1) ) diff --git a/src/user/Rossby_front_2d_initialization.F90 b/src/user/Rossby_front_2d_initialization.F90 index 33c7641a00..09219eb845 100644 --- a/src/user/Rossby_front_2d_initialization.F90 +++ b/src/user/Rossby_front_2d_initialization.F90 @@ -250,6 +250,8 @@ subroutine Rossby_front_initialize_velocity(u, v, h, G, GV, US, param_file, just if (max_depth <= 0.0) call MOM_error(FATAL, & "Rossby_front_initialize_thickness, Rossby_front_initialize_velocity: "//& "This module requires a positive value of MAXIMUM_DEPTH.") + if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, 'Rossby_front_2d_initialization.F90: '// & + "dTdy() is only set to work with Cartesian axis units.") v(:,:,:) = 0.0 u(:,:,:) = 0.0 @@ -335,18 +337,16 @@ end function Hml real function dTdy( G, dT, lat, US ) type(ocean_grid_type), intent(in) :: G !< Grid structure real, intent(in) :: dT !< Top to bottom temperature difference [C ~> degC] - real, intent(in) :: lat !< Latitude in [km] + real, intent(in) :: lat !< Latitude in the same units as geoLat, often [km] type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! Local real :: PI ! The ratio of the circumference of a circle to its diameter [nondim] real :: dHML ! The range of the mixed layer depths [Z ~> m] real :: dHdy ! The mixed layer depth gradient [Z L-1 ~> m m-1] - real :: km_to_L ! Horizontal axis unit conversion factor when AXIS_UNITS = 'k' (1000 m) [L km-1 ~> 1000] PI = 4.0 * atan(1.0) - km_to_L = 1.0e3*US%m_to_L dHML = 0.5 * ( HMLmax - HMLmin ) * G%max_depth - dHdy = dHML * ( PI / ( frontFractionalWidth * G%len_lat * km_to_L ) ) * cos( yPseudo(G, lat) ) + dHdy = dHML * ( PI / ( frontFractionalWidth * G%len_lat * G%grid_unit_to_L ) ) * cos( yPseudo(G, lat) ) dTdy = -( dT / G%max_depth ) * dHdy end function dTdy diff --git a/src/user/soliton_initialization.F90 b/src/user/soliton_initialization.F90 index 722a41b7e5..a734574995 100644 --- a/src/user/soliton_initialization.F90 +++ b/src/user/soliton_initialization.F90 @@ -79,10 +79,12 @@ subroutine soliton_initialize_thickness(h, depth_tot, G, GV, US, param_file, jus if (abs(beta) <= 0.0) call MOM_error(FATAL, & "soliton_initialization, soliton_initialize_thickness: "//& "This module requires a non-zero value of BETA.") + if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, "soliton_initialization.F90: "//& + "soliton_initialize_thickness() is only set to work with Cartesian axis units.") cg_max = sqrt(GV%g_Earth * max_depth) L_eq = sqrt(cg_max / abs(beta)) - scale_pos = US%m_to_L / L_eq + scale_pos = G%grid_unit_to_L / L_eq I_nz = 1.0 / real(nz) x0 = 2.0*G%len_lon/3.0 @@ -150,10 +152,12 @@ subroutine soliton_initialize_velocity(u, v, G, GV, US, param_file, just_read) if (abs(beta) <= 0.0) call MOM_error(FATAL, & "soliton_initialization, soliton_initialize_velocity: "//& "This module requires a non-zero value of BETA.") + if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, "soliton_initialization.F90: "//& + "soliton_initialize_velocity() is only set to work with Cartesian axis units.") cg_max = sqrt(GV%g_Earth * max_depth) L_eq = sqrt(cg_max / abs(beta)) - scale_pos = US%m_to_L / L_eq + scale_pos = G%grid_unit_to_L / L_eq x0 = 2.0*G%len_lon/3.0 y0 = 0.0 From 14f2c97aa9ee2f5d508aa281c8b2fc4a0747bcbd Mon Sep 17 00:00:00 2001 From: Cory Spencer Jones Date: Mon, 23 Dec 2024 15:59:40 -0600 Subject: [PATCH 13/17] Update particles_run call to allow accurate particle advection using u, v or uh, vh The drifters package is able to run in 2 different modes: one where the particles are advected using u and v, and one where they are advected using uh and vh. When u and v are used (i.e. the drifters are advected using the resolved velocities only, with no sub-grd component), the particles package needs the layer thickness that is consistent with these velocities, so particles_run needs to be called before any subgrid scale parameterizations are applied. This was not implemented correctly in my last PR, and is corrected here. The particles package also needs the correct timestep for particle advection. This is added here. --- .../external/drifters/MOM_particles.F90 | 4 +-- src/core/MOM.F90 | 26 +++++++++++++------ 2 files changed, 20 insertions(+), 10 deletions(-) diff --git a/config_src/external/drifters/MOM_particles.F90 b/config_src/external/drifters/MOM_particles.F90 index 95470e6510..1c41170582 100644 --- a/config_src/external/drifters/MOM_particles.F90 +++ b/config_src/external/drifters/MOM_particles.F90 @@ -28,7 +28,7 @@ subroutine particles_init(parts, Grid, Time, dt, u, v, h) end subroutine particles_init !> The main driver the steps updates particles -subroutine particles_run(parts, time, uo, vo, ho, tv, use_uh, stagger) +subroutine particles_run(parts, time, uo, vo, ho, tv, dt_adv, use_uh) ! Arguments type(particles), pointer :: parts !< Container for all types and memory type(time_type), intent(in) :: time !< Model time @@ -40,8 +40,8 @@ subroutine particles_run(parts, time, uo, vo, ho, tv, use_uh, stagger) !! that are used to advect tracers [H L2 ~> m3 or kg] real, dimension(:,:,:), intent(in) :: ho !< Ocean layer thickness [H ~> m or kg m-2] type(thermo_var_ptrs), intent(in) :: tv !< structure containing pointers to available thermodynamic fields + real, intent(in) :: dt_adv !< timestep for advecting particles [s] logical :: use_uh !< Flag for whether u and v are weighted by thickness - integer, optional, intent(in) :: stagger !< Flag for whether velocities are staggered end subroutine particles_run diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index e5d1f46086..8cfab91cfc 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1283,6 +1283,18 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & endif ! -------------------------------------------------- end SPLIT + if (CS%use_particles .and. CS%do_dynamics .and. (.not. CS%use_uh_particles)) then + if (CS%thickness_diffuse_first) call MOM_error(WARNING,"particles_run: "//& + "Thickness_diffuse_first is true and use_uh_particles is false. "//& + "This is usually a bad combination.") + !Run particles using unweighted velocity + call particles_run(CS%particles, Time_local, CS%u, CS%v, CS%h, & + CS%tv, dt, CS%use_uh_particles) + call particles_to_z_space(CS%particles, h) + endif + + + ! Update the model's current to reflect wind-wave growth if (Waves%Stokes_DDT .and. (.not.Waves%Passive_Stokes_DDT)) then do J=jsq,jeq ; do i=is,ie @@ -1368,23 +1380,21 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & endif call disable_averaging(CS%diag) + ! Advance the dynamics time by dt. + CS%t_dyn_rel_adv = CS%t_dyn_rel_adv + dt + if (CS%use_particles .and. CS%do_dynamics .and. CS%use_uh_particles) then !Run particles using thickness-weighted velocity call particles_run(CS%particles, Time_local, CS%uhtr, CS%vhtr, CS%h, & - CS%tv, CS%use_uh_particles) - elseif (CS%use_particles .and. CS%do_dynamics) then - !Run particles using unweighted velocity - call particles_run(CS%particles, Time_local, CS%u, CS%v, CS%h, & - CS%tv, CS%use_uh_particles) + CS%tv, CS%t_dyn_rel_adv, CS%use_uh_particles) endif - - ! Advance the dynamics time by dt. - CS%t_dyn_rel_adv = CS%t_dyn_rel_adv + dt CS%n_dyn_steps_in_adv = CS%n_dyn_steps_in_adv + 1 if (CS%alternate_first_direction) then call set_first_direction(G, MODULO(G%first_direction+1,2)) CS%first_dir_restart = real(G%first_direction) + elseif (CS%use_particles .and. CS%do_dynamics .and. (.not.CS%use_uh_particles)) then + call particles_to_k_space(CS%particles, h) endif CS%t_dyn_rel_thermo = CS%t_dyn_rel_thermo + dt if (abs(CS%t_dyn_rel_thermo) < 1e-6*dt) CS%t_dyn_rel_thermo = 0.0 From 40a59f749f0ffc01c5d82aef06901197a13f2b65 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 8 Dec 2024 10:19:40 -0500 Subject: [PATCH 14/17] +Remove dyn_horgrid_type%Rad_Earth Remove the unused element Rad_Earth from ocean_grid_type and dyn_horgrid_type. The dimensionally rescaled equivalent element G%Rad_Earth_L is extensively used, and it will continue to exist. G%Rad_Earth_L was introduced in November 27, 2021 as a dimensionally rescaled version of G%Rad_Earth, while the unscaled version was retained because at the time, the Rad_Earth element of the dyn_horgrid_type is also used in SIS2. However, SIS2 has not used G%Rad_Earth since December 23, 2021, so after 3 years we can now safely remove this unused element. Any cases on other branches that might be impacted by this change will not compile. All answers are bitwise identical, but there is one less element in two transparent types. --- src/core/MOM_grid.F90 | 1 - src/core/MOM_transcribe_grid.F90 | 4 ++-- src/framework/MOM_dyn_horgrid.F90 | 2 -- src/initialization/MOM_grid_initialize.F90 | 1 - 4 files changed, 2 insertions(+), 6 deletions(-) diff --git a/src/core/MOM_grid.F90 b/src/core/MOM_grid.F90 index 5e5b307103..10335de0a6 100644 --- a/src/core/MOM_grid.F90 +++ b/src/core/MOM_grid.F90 @@ -198,7 +198,6 @@ module MOM_grid real :: west_lon !< The longitude (or x-coordinate) of the first u-line [degrees_E] or [km] or [m] real :: len_lat !< The latitudinal (or y-coord) extent of physical domain [degrees_N] or [km] or [m] real :: len_lon !< The longitudinal (or x-coord) extent of physical domain [degrees_E] or [km] or [m] - real :: Rad_Earth !< The radius of the planet [m] real :: Rad_Earth_L !< The radius of the planet in rescaled units [L ~> m] real :: max_depth !< The maximum depth of the ocean in depth units [Z ~> m] end type ocean_grid_type diff --git a/src/core/MOM_transcribe_grid.F90 b/src/core/MOM_transcribe_grid.F90 index 4c1b87f37e..d9ca19985f 100644 --- a/src/core/MOM_transcribe_grid.F90 +++ b/src/core/MOM_transcribe_grid.F90 @@ -139,7 +139,7 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG, US) oG%areaT_global = dG%areaT_global ; oG%IareaT_global = dG%IareaT_global oG%south_lat = dG%south_lat ; oG%west_lon = dG%west_lon oG%len_lat = dG%len_lat ; oG%len_lon = dG%len_lon - oG%Rad_Earth = dG%Rad_Earth ; oG%Rad_Earth_L = dG%Rad_Earth_L + oG%Rad_Earth_L = dG%Rad_Earth_L oG%max_depth = dG%max_depth ! Update the halos in case the dynamic grid has smaller halos than the ocean grid. @@ -301,7 +301,7 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG, US) dG%areaT_global = oG%areaT_global ; dG%IareaT_global = oG%IareaT_global dG%south_lat = oG%south_lat ; dG%west_lon = oG%west_lon dG%len_lat = oG%len_lat ; dG%len_lon = oG%len_lon - dG%Rad_Earth = oG%Rad_Earth ; dG%Rad_Earth_L = oG%Rad_Earth_L + dG%Rad_Earth_L = oG%Rad_Earth_L dG%max_depth = oG%max_depth ! Update the halos in case the dynamic grid has smaller halos than the ocean grid. diff --git a/src/framework/MOM_dyn_horgrid.F90 b/src/framework/MOM_dyn_horgrid.F90 index 440ea351b9..3b5aeca214 100644 --- a/src/framework/MOM_dyn_horgrid.F90 +++ b/src/framework/MOM_dyn_horgrid.F90 @@ -189,7 +189,6 @@ module MOM_dyn_horgrid real :: west_lon !< The longitude (or x-coordinate) of the first u-line [degrees_E] or [km] or [m] real :: len_lat !< The latitudinal (or y-coord) extent of physical domain [degrees_N] or [km] or [m] real :: len_lon !< The longitudinal (or x-coord) extent of physical domain [degrees_E] or [km] or [m] - real :: Rad_Earth !< The radius of the planet [m] real :: Rad_Earth_L !< The radius of the planet in rescaled units [L ~> m] real :: max_depth !< The maximum depth of the ocean [Z ~> m] end type dyn_horgrid_type @@ -407,7 +406,6 @@ subroutine rotate_dyn_horgrid(G_in, G, US, turns) G%grid_unit_to_L = G_in%grid_unit_to_L G%areaT_global = G_in%areaT_global G%IareaT_global = G_in%IareaT_global - G%Rad_Earth = G_in%Rad_Earth G%Rad_Earth_L = G_in%Rad_Earth_L G%max_depth = G_in%max_depth diff --git a/src/initialization/MOM_grid_initialize.F90 b/src/initialization/MOM_grid_initialize.F90 index 429ce24d79..8b4d1eed0c 100644 --- a/src/initialization/MOM_grid_initialize.F90 +++ b/src/initialization/MOM_grid_initialize.F90 @@ -102,7 +102,6 @@ subroutine set_grid_metrics(G, param_file, US) call get_param(param_file, "MOM_grid_init", "RAD_EARTH", G%Rad_Earth_L, & "The radius of the Earth.", units="m", default=6.378e6, scale=US%m_to_L) endif - G%Rad_Earth = US%L_to_m*G%Rad_Earth_L ! Calculate derived metrics (i.e. reciprocals and products) call callTree_enter("set_derived_metrics(), MOM_grid_initialize.F90") From 73514f2dabfbf257b157eee182d52b27980942f4 Mon Sep 17 00:00:00 2001 From: William Xu Date: Fri, 6 Sep 2024 19:35:15 -0300 Subject: [PATCH 15/17] Nodal modulation This commit fixes a few (potential) inconsistencies between open boundary tidal forcing and astronomical tidal forcing. 1. There was an inconsistency in the code that the nodal modulation can be calculated in OBC tidal forcing but not in the astronomical tidal forcing. This commit fixes this bug so that nodal modulation is applied to both forcings. 2. In the previous version of MOM_open_boundary.F90, a different set of tidal parameters can be set for open boundary tidal forcing, leading to potential inconsistency with astronomical tidal forcing. This commit obsoletes those for open boundary tidal forcing. 3. Another important bug fix is that the equilibrium phase is added to the SAL term, which was missing in the original code. --- src/core/MOM_open_boundary.F90 | 20 +++-- src/diagnostics/MOM_harmonic_analysis.F90 | 26 ++++--- .../lateral/MOM_tidal_forcing.F90 | 77 +++++++++++++++---- 3 files changed, 89 insertions(+), 34 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index fb833189ec..797f60bd9b 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -1169,26 +1169,30 @@ subroutine initialize_obc_tides(OBC, US, param_file) type(astro_longitudes) :: nodal_longitudes !< Solar and lunar longitudes for tidal forcing type(time_type) :: nodal_time !< Model time to calculate nodal modulation for. integer :: c !< Index to tidal constituent. + logical :: tides !< True if astronomical tides are also used. call get_param(param_file, mdl, "OBC_TIDE_CONSTITUENTS", tide_constituent_str, & "Names of tidal constituents being added to the open boundaries.", & fail_if_missing=.true.) - call get_param(param_file, mdl, "OBC_TIDE_ADD_EQ_PHASE", OBC%add_eq_phase, & + call get_param(param_file, mdl, "TIDES", tides, & + "If true, apply tidal momentum forcing.", default=.false., do_not_log=.true.) + + call get_param(param_file, mdl, "TIDE_USE_EQ_PHASE", OBC%add_eq_phase, & "If true, add the equilibrium phase argument to the specified tidal phases.", & - default=.false., fail_if_missing=.false.) + old_name="OBC_TIDE_ADD_EQ_PHASE", default=.false., do_not_log=tides) - call get_param(param_file, mdl, "OBC_TIDE_ADD_NODAL", OBC%add_nodal_terms, & + call get_param(param_file, mdl, "TIDE_ADD_NODAL", OBC%add_nodal_terms, & "If true, include 18.6 year nodal modulation in the boundary tidal forcing.", & - default=.false.) + old_name="OBC_TIDE_ADD_NODAL", default=.false., do_not_log=tides) - call get_param(param_file, mdl, "OBC_TIDE_REF_DATE", tide_ref_date, & + call get_param(param_file, mdl, "TIDE_REF_DATE", tide_ref_date, & "Reference date to use for tidal calculations and equilibrium phase.", & - fail_if_missing=.true.) + old_name="OBC_TIDE_REF_DATE", defaults=(/0, 0, 0/), do_not_log=tides) - call get_param(param_file, mdl, "OBC_TIDE_NODAL_REF_DATE", nodal_ref_date, & + call get_param(param_file, mdl, "TIDE_NODAL_REF_DATE", nodal_ref_date, & "Fixed reference date to use for nodal modulation of boundary tides.", & - fail_if_missing=.false., defaults=(/0, 0, 0/)) + old_name="OBC_TIDE_NODAL_REF_DATE", defaults=(/0, 0, 0/), do_not_log=tides) if (.not. OBC%add_eq_phase) then ! If equilibrium phase argument is not added, the input phases diff --git a/src/diagnostics/MOM_harmonic_analysis.F90 b/src/diagnostics/MOM_harmonic_analysis.F90 index 0013a8ab83..f2585d510a 100644 --- a/src/diagnostics/MOM_harmonic_analysis.F90 +++ b/src/diagnostics/MOM_harmonic_analysis.F90 @@ -46,7 +46,9 @@ module MOM_harmonic_analysis time_ref !< Reference time (t = 0) used to calculate tidal forcing real, dimension(MAX_CONSTITUENTS) :: & freq, & !< The frequency of a tidal constituent [T-1 ~> s-1] - phase0 !< The phase of a tidal constituent at time 0 [rad] + phase0, & !< The phase of a tidal constituent at time 0 [rad] + tide_fn, & !< Amplitude modulation of tides by nodal cycle [nondim]. + tide_un !< Phase modulation of tides by nodal cycle [rad]. real, allocatable :: FtF(:,:) !< Accumulator of (F' * F) for all fields [nondim] integer :: nc !< The number of tidal constituents in use integer :: length !< Number of fields of which harmonic analysis is to be performed @@ -60,13 +62,15 @@ module MOM_harmonic_analysis !> This subroutine sets static variables used by this module and initializes CS%list. !! THIS MUST BE CALLED AT THE END OF tidal_forcing_init. -subroutine HA_init(Time, US, param_file, time_ref, nc, freq, phase0, const_name, CS) +subroutine HA_init(Time, US, param_file, time_ref, nc, freq, phase0, const_name, tide_fn, tide_un, CS) type(time_type), intent(in) :: Time !< The current model time type(time_type), intent(in) :: time_ref !< Reference time (t = 0) used to calculate tidal forcing 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 - real, dimension(MAX_CONSTITUENTS), intent(in) :: freq !< The frequency of a tidal constituent [T-1 ~> s-1] - real, dimension(MAX_CONSTITUENTS), intent(in) :: phase0 !< The phase of a tidal constituent at time 0 [rad] + real, intent(in) :: freq(MAX_CONSTITUENTS) !< The frequency of a tidal constituent [T-1 ~> s-1] + real, intent(in) :: phase0(MAX_CONSTITUENTS) !< The phase of a tidal constituent at time 0 [rad] + real, intent(in) :: tide_fn(MAX_CONSTITUENTS) !< Amplitude modulation of tides by nodal cycle [nondim]. + real, intent(in) :: tide_un(MAX_CONSTITUENTS) !< Phase modulation of tides by nodal cycle [rad]. integer, intent(in) :: nc !< The number of tidal constituents in use character(len=16), intent(in) :: const_name(MAX_CONSTITUENTS) !< The name of each constituent type(harmonic_analysis_CS), intent(out) :: CS !< Control structure of the MOM_harmonic_analysis module @@ -135,6 +139,8 @@ subroutine HA_init(Time, US, param_file, time_ref, nc, freq, phase0, const_name, CS%time_ref = time_ref CS%freq = freq CS%phase0 = phase0 + CS%tide_fn = tide_fn + CS%tide_un = tide_un CS%nc = nc CS%const_name = const_name CS%length = 0 @@ -198,8 +204,8 @@ subroutine HA_accum_FtF(Time, CS) do c=1,nc icos = 2*c isin = 2*c+1 - cosomegat = cos(CS%freq(c) * now + CS%phase0(c)) - sinomegat = sin(CS%freq(c) * now + CS%phase0(c)) + cosomegat = CS%tide_fn(c) * cos(CS%freq(c) * now + (CS%phase0(c) + CS%tide_un(c))) + sinomegat = CS%tide_fn(c) * sin(CS%freq(c) * now + (CS%phase0(c) + CS%tide_un(c))) ! First column, corresponding to the zero frequency constituent (mean) CS%FtF(icos,1) = CS%FtF(icos,1) + cosomegat @@ -208,8 +214,8 @@ subroutine HA_accum_FtF(Time, CS) do cc=1,c iccos = 2*cc issin = 2*cc+1 - ccosomegat = cos(CS%freq(cc) * now + CS%phase0(cc)) - ssinomegat = sin(CS%freq(cc) * now + CS%phase0(cc)) + ccosomegat = CS%tide_fn(cc) * cos(CS%freq(cc) * now + (CS%phase0(cc) + CS%tide_un(cc))) + ssinomegat = CS%tide_fn(cc) * sin(CS%freq(cc) * now + (CS%phase0(cc) + CS%tide_un(cc))) ! Interior of the matrix, corresponding to the products of cosine and sine terms CS%FtF(icos,iccos) = CS%FtF(icos,iccos) + cosomegat * ccosomegat @@ -290,8 +296,8 @@ subroutine HA_accum_FtSSH(key, data, Time, G, CS) do c=1,nc icos = 2*c isin = 2*c+1 - cosomegat = cos(CS%freq(c) * now + CS%phase0(c)) - sinomegat = sin(CS%freq(c) * now + CS%phase0(c)) + cosomegat = CS%tide_fn(c) * cos(CS%freq(c) * now + (CS%phase0(c) + CS%tide_un(c))) + sinomegat = CS%tide_fn(c) * sin(CS%freq(c) * now + (CS%phase0(c) + CS%tide_un(c))) do j=js,je ; do i=is,ie ha1%FtSSH(i,j,icos) = ha1%FtSSH(i,j,icos) + (data(i,j) - ha1%ref(i,j)) * cosomegat ha1%FtSSH(i,j,isin) = ha1%FtSSH(i,j,isin) + (data(i,j) - ha1%ref(i,j)) * sinomegat diff --git a/src/parameterizations/lateral/MOM_tidal_forcing.F90 b/src/parameterizations/lateral/MOM_tidal_forcing.F90 index 85c9b1ee81..5e624d4aea 100644 --- a/src/parameterizations/lateral/MOM_tidal_forcing.F90 +++ b/src/parameterizations/lateral/MOM_tidal_forcing.F90 @@ -70,7 +70,9 @@ module MOM_tidal_forcing ampsal(:,:,:), & !< The amplitude of the SAL [Z ~> m]. cosphase_prev(:,:,:), & !< The cosine of the phase of the amphidromes in the previous tidal solutions [nondim]. sinphase_prev(:,:,:), & !< The sine of the phase of the amphidromes in the previous tidal solutions [nondim]. - amp_prev(:,:,:) !< The amplitude of the previous tidal solution [Z ~> m]. + amp_prev(:,:,:), & !< The amplitude of the previous tidal solution [Z ~> m]. + tide_fn(:), & !< Amplitude modulation of tides by nodal cycle [nondim]. + tide_un(:) !< Phase modulation of tides by nodal cycle [rad]. end type tidal_forcing_CS integer :: id_clock_tides !< CPU clock for tides @@ -251,9 +253,14 @@ subroutine tidal_forcing_init(Time, G, US, param_file, CS, HA_CS) 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. + integer, dimension(3) :: nodal_ref_date !< Reference date for calculating nodal modulation for tidal forcing. 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 :: add_nodal_terms = .false. !< If true, insert terms for the 18.6 year modulation when + !! calculating tidal forcing. + type(time_type) :: nodal_time !< Model time to calculate nodal modulation for. + type(astro_longitudes) :: nodal_longitudes !< Solar and lunar longitudes for tidal forcing logical :: HA_ssh, HA_ubt, HA_vbt ! This include declares and sets the variable "version". # include "version_variable.h" @@ -385,11 +392,11 @@ subroutine tidal_forcing_init(Time, G, US, param_file, CS, HA_CS) call get_param(param_file, mdl, "TIDE_REF_DATE", tide_ref_date, & "Year,month,day to use as reference date for tidal forcing. "//& "If not specified, defaults to 0.", & - defaults=(/0, 0, 0/)) + old_name="OBC_TIDE_REF_DATE", defaults=(/0, 0, 0/)) call get_param(param_file, mdl, "TIDE_USE_EQ_PHASE", CS%use_eq_phase, & "Correct phases by calculating equilibrium phase arguments for TIDE_REF_DATE. ", & - default=.false., fail_if_missing=.false.) + old_name="OBC_TIDE_ADD_EQ_PHASE", default=.false., fail_if_missing=.false.) if (sum(tide_ref_date) == 0) then ! tide_ref_date defaults to 0. CS%time_ref = set_date(1, 1, 1, 0, 0, 0) @@ -528,8 +535,46 @@ subroutine tidal_forcing_init(Time, G, US, param_file, CS, HA_CS) enddo endif + call get_param(param_file, mdl, "TIDE_ADD_NODAL", add_nodal_terms, & + "If true, include 18.6 year nodal modulation in the astronomical tidal forcing.", & + old_name="OBC_TIDE_ADD_NODAL", default=.false.) + call get_param(param_file, mdl, "TIDE_NODAL_REF_DATE", nodal_ref_date, & + "Fixed reference date to use for nodal modulation of astronomical tidal forcing.", & + old_name="OBC_TIDE_REF_DATE", fail_if_missing=.false., defaults=(/0, 0, 0/)) + + ! If the nodal correction is based on a different time, initialize that. + ! Otherwise, it can use N from the time reference. + if (add_nodal_terms) then + if (sum(nodal_ref_date) /= 0) then + ! A reference date was provided for the nodal correction + nodal_time = set_date(nodal_ref_date(1), nodal_ref_date(2), nodal_ref_date(3)) + call astro_longitudes_init(nodal_time, nodal_longitudes) + elseif (CS%use_eq_phase) then + ! Astronomical longitudes were already calculated for use in equilibrium phases, + ! so use nodal longitude from that. + nodal_longitudes = CS%tidal_longitudes + else + ! Tidal reference time is a required parameter, so calculate the longitudes from that. + call astro_longitudes_init(CS%time_ref, nodal_longitudes) + endif + endif + + allocate(CS%tide_fn(nc)) + allocate(CS%tide_un(nc)) + + do c=1,nc + ! Find nodal corrections if needed + if (add_nodal_terms) then + call nodal_fu(trim(CS%const_name(c)), nodal_longitudes%N, CS%tide_fn(c), CS%tide_un(c)) + else + CS%tide_fn(c) = 1.0 + CS%tide_un(c) = 0.0 + endif + enddo + if (present(HA_CS)) then - call HA_init(Time, US, param_file, CS%time_ref, CS%nc, CS%freq, CS%phase0, CS%const_name, HA_CS) + call HA_init(Time, US, param_file, CS%time_ref, CS%nc, CS%freq, CS%phase0, CS%const_name, & + CS%tide_fn, CS%tide_un, HA_CS) call get_param(param_file, mdl, "HA_SSH", HA_ssh, & "If true, perform harmonic analysis of sea serface height.", default=.false.) if (HA_ssh) call HA_register('ssh', 'h', HA_CS) @@ -614,8 +659,8 @@ subroutine calc_tidal_forcing(Time, e_tide_eq, e_tide_sal, G, US, CS) do c=1,CS%nc m = CS%struct(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)) + amp_cosomegat = CS%amp(c)*CS%love_no(c)*CS%tide_fn(c) * cos(CS%freq(c)*now + (CS%phase0(c) + CS%tide_un(c))) + amp_sinomegat = CS%amp(c)*CS%love_no(c)*CS%tide_fn(c) * sin(CS%freq(c)*now + (CS%phase0(c) + CS%tide_un(c))) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 e_tide_eq(i,j) = e_tide_eq(i,j) + (amp_cosomegat*CS%cos_struct(i,j,m) + & amp_sinomegat*CS%sin_struct(i,j,m)) @@ -623,8 +668,8 @@ subroutine calc_tidal_forcing(Time, e_tide_eq, e_tide_sal, G, US, CS) enddo if (CS%use_tidal_sal_file) then ; do c=1,CS%nc - cosomegat = cos(CS%freq(c)*now) - sinomegat = sin(CS%freq(c)*now) + cosomegat = CS%tide_fn(c) * cos(CS%freq(c)*now + (CS%phase0(c) + CS%tide_un(c))) + sinomegat = CS%tide_fn(c) * sin(CS%freq(c)*now + (CS%phase0(c) + CS%tide_un(c))) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 e_tide_sal(i,j) = e_tide_sal(i,j) + CS%ampsal(i,j,c) * & (cosomegat*CS%cosphasesal(i,j,c) + sinomegat*CS%sinphasesal(i,j,c)) @@ -632,8 +677,8 @@ subroutine calc_tidal_forcing(Time, e_tide_eq, e_tide_sal, G, US, CS) enddo ; endif if (CS%use_tidal_sal_prev) then ; do c=1,CS%nc - cosomegat = cos(CS%freq(c)*now) - sinomegat = sin(CS%freq(c)*now) + cosomegat = CS%tide_fn(c) * cos(CS%freq(c)*now + (CS%phase0(c) + CS%tide_un(c))) + sinomegat = CS%tide_fn(c) * sin(CS%freq(c)*now + (CS%phase0(c) + CS%tide_un(c))) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 e_tide_sal(i,j) = e_tide_sal(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)) @@ -692,8 +737,8 @@ subroutine calc_tidal_forcing_legacy(Time, e_sal, e_sal_tide, e_tide_eq, e_tide_ do c=1,CS%nc m = CS%struct(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)) + amp_cosomegat = CS%amp(c)*CS%love_no(c)*CS%tide_fn(c) * cos(CS%freq(c)*now + (CS%phase0(c) + CS%tide_un(c))) + amp_sinomegat = CS%amp(c)*CS%love_no(c)*CS%tide_fn(c) * sin(CS%freq(c)*now + (CS%phase0(c) + CS%tide_un(c))) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 amp_cossin = (amp_cosomegat*CS%cos_struct(i,j,m) + amp_sinomegat*CS%sin_struct(i,j,m)) e_sal_tide(i,j) = e_sal_tide(i,j) + amp_cossin @@ -702,8 +747,8 @@ subroutine calc_tidal_forcing_legacy(Time, e_sal, e_sal_tide, e_tide_eq, e_tide_ enddo if (CS%use_tidal_sal_file) then ; do c=1,CS%nc - cosomegat = cos(CS%freq(c)*now) - sinomegat = sin(CS%freq(c)*now) + cosomegat = CS%tide_fn(c) * cos(CS%freq(c)*now + (CS%phase0(c) + CS%tide_un(c))) + sinomegat = CS%tide_fn(c) * sin(CS%freq(c)*now + (CS%phase0(c) + CS%tide_un(c))) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 amp_cossin = CS%ampsal(i,j,c) & * (cosomegat*CS%cosphasesal(i,j,c) + sinomegat*CS%sinphasesal(i,j,c)) @@ -713,8 +758,8 @@ subroutine calc_tidal_forcing_legacy(Time, e_sal, e_sal_tide, e_tide_eq, e_tide_ enddo ; endif if (CS%use_tidal_sal_prev) then ; do c=1,CS%nc - cosomegat = cos(CS%freq(c)*now) - sinomegat = sin(CS%freq(c)*now) + cosomegat = CS%tide_fn(c) * cos(CS%freq(c)*now + (CS%phase0(c) + CS%tide_un(c))) + sinomegat = CS%tide_fn(c) * sin(CS%freq(c)*now + (CS%phase0(c) + CS%tide_un(c))) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 amp_cossin = -CS%sal_scalar * CS%amp_prev(i,j,c) & * (cosomegat*CS%cosphase_prev(i,j,c) + sinomegat*CS%sinphase_prev(i,j,c)) From 54feb6fc40400e96d9f18e64a2d3ed08a1419826 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 12 Dec 2024 09:31:49 -0500 Subject: [PATCH 16/17] Correct unit conversion for BS_coeff_h diagnostic Added missing conversion arguments for the register_diag_field calls for the recently added diagnostics BS_coeff_h and BS_coeff_q. All answers are bitwise identical, but two diagnostics will have corrected dimensional rescaling when EY24_EBT_BS is true. --- src/parameterizations/lateral/MOM_hor_visc.F90 | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index a11b8a232e..910ce067be 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -1818,14 +1818,12 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, else if (use_kh_struct) then Kh_BS(I,J) = 0.25*( ((MEKE%Ku(i,j)*VarMix%BS_struct(i,j,k)) + & - (MEKE%Ku(i+1,j+1)*VarMix%BS_struct(i+1,j+1,k))) + & - ((MEKE%Ku(i+1,j)*VarMix%BS_struct(i+1,j,k)) + & - (MEKE%Ku(i,j+1)*VarMix%BS_struct(i,j+1,k))) ) + (MEKE%Ku(i+1,j+1)*VarMix%BS_struct(i+1,j+1,k))) + & + ((MEKE%Ku(i+1,j)*VarMix%BS_struct(i+1,j,k)) + & + (MEKE%Ku(i,j+1)*VarMix%BS_struct(i,j+1,k))) ) else - Kh_BS(I,J) = 0.25*( (MEKE%Ku(i,j) + & - MEKE%Ku(i+1,j+1)) + & - (MEKE%Ku(i+1,j) + & - MEKE%Ku(i,j+1)) ) + Kh_BS(I,J) = 0.25*( (MEKE%Ku(i,j) + MEKE%Ku(i+1,j+1)) + & + (MEKE%Ku(i+1,j) + MEKE%Ku(i,j+1)) ) endif endif enddo ; enddo @@ -3226,9 +3224,9 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) if (CS%EY24_EBT_BS) then CS%id_BS_coeff_h = register_diag_field('ocean_model', 'BS_coeff_h', diag%axesTL, Time, & - 'Backscatter coefficient at h points', 'm2 s-1') + 'Backscatter coefficient at h points', units='m2 s-1', conversion=US%L_to_m**2*US%s_to_T) CS%id_BS_coeff_q = register_diag_field('ocean_model', 'BS_coeff_q', diag%axesBL, Time, & - 'Backscatter coefficient at q points', 'm2 s-1') + 'Backscatter coefficient at q points', units='m2 s-1', conversion=US%L_to_m**2*US%s_to_T) endif CS%id_FrictWork = register_diag_field('ocean_model','FrictWork',diag%axesTL,Time,& From 83efb99bf298a09c15c5ec5bfbe88479b93db126 Mon Sep 17 00:00:00 2001 From: William Xu Date: Tue, 10 Dec 2024 09:03:13 -0400 Subject: [PATCH 17/17] Streaming Filter The filters and their target frequencies are no longer hard-coded. Instead, multiple filters with tidal or arbitrary frequencies as their target frequencies can be turned on. The filter names are specified in MOM_input and must consist of two letters/numbers. If the name of a filter is the same as the name of a tidal constituent, then the corresponding tidal frequency will be used as its target frequency. Otherwise, the user must specify the target frequency by "FILTER_${FILTER_NAME}_OMEGA" in MOM_input. The restarting capability has also been implemented. Because the filtering is a point-wise operation, all variables are considered as fields, even if they are velocity components. --- src/core/MOM_barotropic.F90 | 96 +++----- .../lateral/MOM_streaming_filter.F90 | 212 +++++++++++++----- 2 files changed, 181 insertions(+), 127 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 70e0a738c4..f3d5564ee0 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -26,8 +26,7 @@ module MOM_barotropic use MOM_restart, only : query_initialized, MOM_restart_CS use MOM_self_attr_load, only : scalar_SAL_sensitivity use MOM_self_attr_load, only : SAL_CS -use MOM_streaming_filter, only : Filt_register, Filt_accum, Filter_CS -use MOM_tidal_forcing, only : tidal_frequency +use MOM_streaming_filter, only : Filt_register, Filt_init, Filt_accum, Filter_CS use MOM_time_manager, only : time_type, real_to_time, operator(+), operator(-) use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : BT_cont_type, alloc_bt_cont_type @@ -250,10 +249,8 @@ module MOM_barotropic logical :: linearized_BT_PV !< If true, the PV and interface thicknesses used !! in the barotropic Coriolis calculation is time !! invariant and linearized. - logical :: use_filter_m2 !< If true, apply streaming band-pass filter for detecting - !! instantaneous tidal signals. - logical :: use_filter_k1 !< If true, apply streaming band-pass filter for detecting - !! instantaneous tidal signals. + logical :: use_filter !< If true, use streaming band-pass filter to detect the + !! instantaneous tidal signals in the simulation. logical :: use_wide_halos !< If true, use wide halos and march in during the !! barotropic time stepping for efficiency. logical :: clip_velocity !< If true, limit any velocity components that are @@ -297,10 +294,8 @@ module MOM_barotropic type(hor_index_type), pointer :: debug_BT_HI => NULL() !< debugging copy of horizontal index_type type(SAL_CS), pointer :: SAL_CSp => NULL() !< Control structure for SAL type(harmonic_analysis_CS), pointer :: HA_CSp => NULL() !< Control structure for harmonic analysis - type(Filter_CS) :: Filt_CS_um2, & !< Control structures for the M2 streaming filter - Filt_CS_vm2, & !< Control structures for the M2 streaming filter - Filt_CS_uk1, & !< Control structures for the K1 streaming filter - Filt_CS_vk1 !< Control structures for the K1 streaming filter + type(Filter_CS) :: Filt_CS_u, & !< Control structures for the streaming band-pass filter of ubt + Filt_CS_v !< Control structures for the streaming band-pass filter of vbt logical :: module_is_initialized = .false. !< If true, module has been initialized integer :: isdw !< The lower i-memory limit for the wide halo arrays. @@ -608,8 +603,8 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, DCor_v, & ! An averaged total thickness at v points [H ~> m or kg m-2]. Datv ! Basin depth at v-velocity grid points times the x-grid ! spacing [H L ~> m2 or kg m-1]. - real, dimension(:,:), pointer :: um2, uk1, vm2, vk1 - ! M2 and K1 velocities from the output of streaming filters [m s-1] + real, dimension(:,:,:), pointer :: ufilt, vfilt + ! Filtered velocities from the output of streaming filters [m s-1] real, target, dimension(SZIW_(CS),SZJW_(CS)) :: & eta, & ! The barotropic free surface height anomaly or column mass ! anomaly [H ~> m or kg m-2] @@ -1598,15 +1593,11 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, endif ; enddo ; enddo endif - ! Here is an example of how the filter equations are time stepped to determine the M2 and K1 velocities. - ! The filters are initialized and registered in subroutine barotropic_init. - if (CS%use_filter_m2) then - call Filt_accum(ubt, um2, CS%Time, US, CS%Filt_CS_um2) - call Filt_accum(vbt, vm2, CS%Time, US, CS%Filt_CS_vm2) - endif - if (CS%use_filter_k1) then - call Filt_accum(ubt, uk1, CS%Time, US, CS%Filt_CS_uk1) - call Filt_accum(vbt, vk1, CS%Time, US, CS%Filt_CS_vk1) + ! Note that the filtered velocities are only updated during the current predictor step, + ! and are calculated using the barotropic velocity from the previous correction step. + if (CS%use_filter) then + call Filt_accum(ubt, ufilt, CS%Time, US, CS%Filt_CS_u) + call Filt_accum(vbt, vfilt, CS%Time, US, CS%Filt_CS_v) endif ! Zero out the arrays for various time-averaged quantities. @@ -4979,6 +4970,12 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, endif ! len_trim(wave_drag_file) > 0 endif ! CS%linear_wave_drag + ! Initialize streaming band-pass filters + if (CS%use_filter) then + call Filt_init(param_file, US, CS%Filt_CS_u, restart_CS) + call Filt_init(param_file, US, CS%Filt_CS_v, restart_CS) + endif + CS%dtbt_fraction = 0.98 ; if (dtbt_input < 0.0) CS%dtbt_fraction = -dtbt_input dtbt_tmp = -1.0 @@ -5263,9 +5260,8 @@ subroutine register_barotropic_restarts(HI, GV, US, param_file, CS, restart_CS) ! Local variables type(vardesc) :: vd(3) character(len=40) :: mdl = "MOM_barotropic" ! This module's name. + integer :: n_filters !< Number of streaming band-pass filters to be used in the simulation. integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB - real :: am2, ak1 !< Bandwidth parameters of the M2 and K1 streaming filters [nondim] - real :: om2, ok1 !< Target frequencies of the M2 and K1 streaming filters [rad T-1 ~> rad s-1] isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed IsdB = HI%IsdB ; IedB = HI%IedB ; JsdB = HI%JsdB ; JedB = HI%JedB @@ -5278,33 +5274,6 @@ subroutine register_barotropic_restarts(HI, GV, US, param_file, CS, restart_CS) "sum(u dh_dt) while also correcting for truncation errors.", & default=.false., do_not_log=.true.) - call get_param(param_file, mdl, "STREAMING_FILTER_M2", CS%use_filter_m2, & - "If true, turn on streaming band-pass filter for detecting "//& - "instantaneous tidal signals.", default=.false.) - call get_param(param_file, mdl, "STREAMING_FILTER_K1", CS%use_filter_k1, & - "If true, turn on streaming band-pass filter for detecting "//& - "instantaneous tidal signals.", default=.false.) - call get_param(param_file, mdl, "FILTER_ALPHA_M2", am2, & - "Bandwidth parameter of the streaming filter targeting the M2 frequency. "//& - "Must be positive. To turn off filtering, set FILTER_ALPHA_M2 <= 0.0.", & - default=0.0, units="nondim", do_not_log=.not.CS%use_filter_m2) - call get_param(param_file, mdl, "FILTER_ALPHA_K1", ak1, & - "Bandwidth parameter of the streaming filter targeting the K1 frequency. "//& - "Must be positive. To turn off filtering, set FILTER_ALPHA_K1 <= 0.0.", & - default=0.0, units="nondim", do_not_log=.not.CS%use_filter_k1) - call get_param(param_file, mdl, "TIDE_M2_FREQ", om2, & - "Frequency of the M2 tidal constituent. "//& - "This is only used if TIDES and TIDE_M2"// & - " are true, or if OBC_TIDE_N_CONSTITUENTS > 0 and M2"// & - " is in OBC_TIDE_CONSTITUENTS.", units="rad s-1", default=tidal_frequency("M2"), & - scale=US%T_to_s, do_not_log=.true.) - call get_param(param_file, mdl, "TIDE_K1_FREQ", ok1, & - "Frequency of the K1 tidal constituent. "//& - "This is only used if TIDES and TIDE_K1"// & - " are true, or if OBC_TIDE_N_CONSTITUENTS > 0 and K1"// & - " is in OBC_TIDE_CONSTITUENTS.", units="rad s-1", default=tidal_frequency("K1"), & - scale=US%T_to_s, do_not_log=.true.) - ALLOC_(CS%ubtav(IsdB:IedB,jsd:jed)) ; CS%ubtav(:,:) = 0.0 ALLOC_(CS%vbtav(isd:ied,JsdB:JedB)) ; CS%vbtav(:,:) = 0.0 if (CS%gradual_BT_ICs) then @@ -5333,22 +5302,17 @@ subroutine register_barotropic_restarts(HI, GV, US, param_file, CS, restart_CS) call register_restart_field(CS%dtbt, "DTBT", .false., restart_CS, & longname="Barotropic timestep", units="seconds", conversion=US%T_to_s) - ! Initialize and register streaming filters - if (CS%use_filter_m2) then - if (am2 > 0.0 .and. om2 > 0.0) then - call Filt_register(am2, om2, 'u', HI, CS%Filt_CS_um2) - call Filt_register(am2, om2, 'v', HI, CS%Filt_CS_vm2) - else - CS%use_filter_m2 = .false. - endif - endif - if (CS%use_filter_k1) then - if (ak1 > 0.0 .and. ok1 > 0.0) then - call Filt_register(ak1, ok1, 'u', HI, CS%Filt_CS_uk1) - call Filt_register(ak1, ok1, 'v', HI, CS%Filt_CS_vk1) - else - CS%use_filter_k1 = .false. - endif + ! Initialize and register streaming band-pass filters + call get_param(param_file, mdl, "USE_FILTER", CS%use_filter, & + "If true, use streaming band-pass filters to detect the "//& + "instantaneous tidal signals in the simulation.", default=.false.) + call get_param(param_file, mdl, "N_FILTERS", n_filters, & + "Number of streaming band-pass filters to be used in the simulation.", & + default=0, do_not_log=.not.CS%use_filter) + if (n_filters<=0) CS%use_filter = .false. + if (CS%use_filter) then + call Filt_register(n_filters, 'ubt', 'u', HI, CS%Filt_CS_u, restart_CS) + call Filt_register(n_filters, 'vbt', 'v', HI, CS%Filt_CS_v, restart_CS) endif end subroutine register_barotropic_restarts diff --git a/src/parameterizations/lateral/MOM_streaming_filter.F90 b/src/parameterizations/lateral/MOM_streaming_filter.F90 index f7a4f736c7..a2adaa2e92 100644 --- a/src/parameterizations/lateral/MOM_streaming_filter.F90 +++ b/src/parameterizations/lateral/MOM_streaming_filter.F90 @@ -1,119 +1,209 @@ !> Streaming band-pass filter for detecting the instantaneous tidal signals in the simulation + module MOM_streaming_filter -use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL +use MOM_error_handler, only : MOM_mesg, MOM_error, NOTE, FATAL +use MOM_file_parser, only : get_param, param_file_type use MOM_hor_index, only : hor_index_type +use MOM_io, only : axis_info, set_axis_info +use MOM_restart, only : register_restart_field, query_initialized, MOM_restart_CS +use MOM_tidal_forcing, only : tidal_frequency use MOM_time_manager, only : time_type, time_type_to_real use MOM_unit_scaling, only : unit_scale_type implicit none ; private -public Filt_register, Filt_accum +public Filt_register, Filt_init, Filt_accum #include -!> The control structure for storing the filter infomation of a particular field +!> Control structure for the MOM_streaming_filter module type, public :: Filter_CS ; private - real :: a, & !< Parameter that determines the bandwidth [nondim] - om, & !< Target frequency of the filter [rad T-1 ~> rad s-1] - old_time = -1.0 !< The time of the previous accumulating step [T ~> s] - real, allocatable, dimension(:,:) :: s1, & !< Dummy variable [A] - u1 !< Filtered data [A] + integer :: nf !< Number of filters to be used in the simulation !>@{ Lower and upper bounds of input data integer :: is, ie, js, je !>@} + character(len=8) :: key !< Identifier of the variable to be filtered + character(len=2), allocatable, dimension(:) :: filter_names !< Names of filters + real, allocatable, dimension(:) :: filter_omega !< Target frequencies of filters [rad T-1 ~> rad s-1] + real, allocatable, dimension(:) :: filter_alpha !< Bandwidth parameters of filters [nondim] + real, allocatable, dimension(:,:,:) :: s1, & !< A dummy variable for solving the system of ODEs [A] + u1 !< Filtered data, representing the narrow-band signal + !< oscillating around the target frequency [A] + real :: old_time = -1.0 !< The time of the previous accumulating step [T ~> s] end type Filter_CS contains -!> This subroutine registers each of the fields to be filtered. -subroutine Filt_register(a, om, grid, HI, CS) - real, intent(in) :: a !< Parameter that determines the bandwidth [nondim] - real, intent(in) :: om !< Target frequency of the filter [rad T-1 ~> rad s-1] - character(len=*), intent(in) :: grid !< Horizontal grid location: h, u, or v - type(hor_index_type), intent(in) :: HI !< Horizontal index type structure - type(Filter_CS), intent(out) :: CS !< Control structure for the current field +!> This subroutine registers the filter variables given the number of filters and the grid +subroutine Filt_register(nf, key, grid, HI, CS, restart_CS) + integer, intent(in) :: nf !< Number of filters to be used in the simulation + character(len=*), intent(in) :: key !< Identifier of the variable to be filtered + character(len=*), intent(in) :: grid !< Horizontal grid location: "h", "u", or "v" + type(hor_index_type), intent(in) :: HI !< Horizontal index type structure + type(Filter_CS), intent(out) :: CS !< Control structure of MOM_streaming_filter + type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control structure ! Local variables - integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB - - if (a <= 0.0) call MOM_error(FATAL, "MOM_streaming_filter: bandwidth <= 0") - if (om <= 0.0) call MOM_error(FATAL, "MOM_streaming_filter: target frequency <= 0") - - CS%a = a - CS%om = om + type(axis_info) :: filter_axis(1) + real, dimension(:), allocatable :: n_filters !< Labels of filters [nondim] + integer :: c - isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed - IsdB = HI%IsdB ; IedB = HI%IedB ; JsdB = HI%JsdB ; JedB = HI%JedB + CS%nf = nf + CS%key = key select case (trim(grid)) case ('h') - allocate(CS%s1(isd:ied,jsd:jed)) ; CS%s1(:,:) = 0.0 - allocate(CS%u1(isd:ied,jsd:jed)) ; CS%u1(:,:) = 0.0 - CS%is = isd ; CS%ie = ied ; CS%js = jsd ; CS%je = jed + CS%is = HI%isd ; CS%ie = HI%ied ; CS%js = HI%jsd ; CS%je = HI%jed case ('u') - allocate(CS%s1(IsdB:IedB,jsd:jed)) ; CS%s1(:,:) = 0.0 - allocate(CS%u1(IsdB:IedB,jsd:jed)) ; CS%u1(:,:) = 0.0 - CS%is = IsdB ; CS%ie = IedB ; CS%js = jsd ; CS%je = jed + CS%is = HI%IsdB ; CS%ie = HI%IedB ; CS%js = HI%jsd ; CS%je = HI%jed case ('v') - allocate(CS%s1(isd:ied,JsdB:JedB)) ; CS%s1(:,:) = 0.0 - allocate(CS%u1(isd:ied,JsdB:JedB)) ; CS%u1(:,:) = 0.0 - CS%is = isd ; CS%ie = ied ; CS%js = JsdB ; CS%je = JedB + CS%is = HI%isd ; CS%ie = HI%ied ; CS%js = HI%JsdB ; CS%je = HI%JedB case default call MOM_error(FATAL, "MOM_streaming_filter: horizontal grid not supported") end select + allocate(CS%s1(CS%is:CS%ie, CS%js:CS%je, nf), source=0.0) + allocate(CS%u1(CS%is:CS%ie, CS%js:CS%je, nf), source=0.0) + + ! Register restarts for s1 and u1 + allocate(n_filters(nf)) + + do c=1,nf ; n_filters(c) = c ; enddo + + call set_axis_info(filter_axis(1), "n_filters", "", "number of filters", nf, n_filters, "N", 1) + + call register_restart_field(CS%s1(:,:,:), "Filter_"//trim(key)//"_s1", .false., restart_CS, & + longname="Dummy variable for streaming band-pass filter", & + hor_grid=trim(grid), z_grid="1", t_grid="s", extra_axes=filter_axis) + call register_restart_field(CS%u1(:,:,:), "Filter_"//trim(key)//"_u1", .false., restart_CS, & + longname="Output of streaming band-pass filter", & + hor_grid=trim(grid), z_grid="1", t_grid="s", extra_axes=filter_axis) + end subroutine Filt_register -!> This subroutine timesteps the filter equations. It takes model output u at the current time step as the input, -!! and returns tidal signal u1 as the output, which is the solution of a set of two ODEs (the filter equations). +!> This subroutine initializes the filters +subroutine Filt_init(param_file, US, CS, restart_CS) + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(Filter_CS), intent(inout) :: CS !< Control structure of MOM_streaming_filter + type(MOM_restart_CS), intent(in) :: restart_CS !< MOM restart control structure + + ! Local variables + character(len=40) :: mdl = "MOM_streaming_filter" !< This module's name + character(len=50) :: filter_name_str !< List of filters to be registered + character(len=200) :: mesg + integer :: c + + call get_param(param_file, mdl, "FILTER_NAMES", filter_name_str, & + "Names of streaming band-pass filters to be used in the simulation.", & + fail_if_missing=.true.) + allocate(CS%filter_names(CS%nf)) + allocate(CS%filter_omega(CS%nf)) + allocate(CS%filter_alpha(CS%nf)) + read(filter_name_str, *) CS%filter_names + + do c=1,CS%nf + ! If filter_name_str consists of tidal constituents, use tidal frequencies. + call get_param(param_file, mdl, "FILTER_"//trim(CS%filter_names(c))//"_OMEGA", & + CS%filter_omega(c), "Target frequency of the "//trim(CS%filter_names(c))//& + " filter. This is used if USE_FILTER is true and "//trim(CS%filter_names(c))//& + " is in FILTER_NAMES.", units="rad s-1", scale=US%T_to_s, default=0.0) + call get_param(param_file, mdl, "FILTER_"//trim(CS%filter_names(c))//"_ALPHA", & + CS%filter_alpha(c), "Bandwidth parameter of the "//trim(CS%filter_names(c))//& + " filter. Must be positive.", units="nondim", fail_if_missing=.true.) + + if (CS%filter_omega(c)<=0.0) CS%filter_omega(c) = tidal_frequency(trim(CS%filter_names(c))) + if (CS%filter_alpha(c)<=0.0) call MOM_error(FATAL, "MOM_streaming_filter: bandwidth <= 0") + + write(mesg,*) "MOM_streaming_filter: ", trim(CS%filter_names(c)), & + " filter registered, target frequency = ", CS%filter_omega(c), & + ", bandwidth = ", CS%filter_alpha(c) + call MOM_error(NOTE, trim(mesg)) + enddo + + if (query_initialized(CS%s1, "Filter_"//trim(CS%key)//"_s1", restart_CS)) then + write(mesg,*) "MOM_streaming_filter: Dummy variable for filter ", trim(CS%key), & + " found in restart files." + else + write(mesg,*) "MOM_streaming_filter: Dummy variable for filter ", trim(CS%key), & + " not found in restart files. The filter will spin up from zeros." + endif + call MOM_error(NOTE, trim(mesg)) + + if (query_initialized(CS%u1, "Filter_"//trim(CS%key)//"_u1", restart_CS)) then + write(mesg,*) "MOM_streaming_filter: Output of filter ", trim(CS%key), & + " found in restart files." + else + write(mesg,*) "MOM_streaming_filter: Output of filter ", trim(CS%key), & + " not found in restart files. The filter will spin up from zeros." + endif + call MOM_error(NOTE, trim(mesg)) + +end subroutine Filt_init + +!> This subroutine timesteps the filter equations. Here, u is the broadband input signal from the model, +!! and u1 is the filtered, narrowband output signal, obtained from the solution of the filter equations. subroutine Filt_accum(u, u1, Time, US, CS) - real, dimension(:,:), pointer, intent(out) :: u1 !< Output of the filter [A] - type(time_type), intent(in) :: Time !< The current model time - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - type(Filter_CS), target, intent(inout) :: CS !< Control structure of the MOM_streaming_filter module + real, dimension(:,:,:), pointer, intent(out) :: u1 !< Output of the filter [A] + type(time_type), intent(in) :: Time !< The current model time + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(Filter_CS), target, intent(inout) :: CS !< Control structure of MOM_streaming_filter real, dimension(CS%is:CS%ie,CS%js:CS%je), intent(in) :: u !< Input into the filter [A] ! Local variables real :: now, & !< The current model time [T ~> s] dt, & !< Time step size for the filter equations [T ~> s] c1, c2 !< Coefficients for the filter equations [nondim] - integer :: i, j, is, ie, js, je + integer :: i, j, k now = US%s_to_T * time_type_to_real(Time) - is = CS%is ; ie = CS%ie ; js = CS%js ; je = CS%je - ! Initialize u1 - if (CS%old_time < 0.0) then - CS%old_time = now - CS%u1(:,:) = u(:,:) - endif + ! Initialize CS%old_time at the first time step + if (CS%old_time<0.0) CS%old_time = now + + ! This condition implies Filt_accum has already been called in the predictor step + if (CS%old_time>=now) return dt = now - CS%old_time CS%old_time = now ! Timestepping - c1 = CS%om * dt - c2 = 1.0 - CS%a * c1 - - do j=js,je ; do i=is,ie - CS%s1(i,j) = c1 * CS%u1(i,j) + CS%s1(i,j) - CS%u1(i,j) = -c1 * (CS%s1(i,j) - CS%a * u(i,j)) + c2 * CS%u1(i,j) - enddo; enddo + do k=1,CS%nf + c1 = CS%filter_omega(k) * dt + c2 = 1.0 - CS%filter_alpha(k) * c1 + + do j=CS%js,CS%je ; do i=CS%is,CS%ie + CS%s1(i,j,k) = c1 * CS%u1(i,j,k) + CS%s1(i,j,k) + CS%u1(i,j,k) = -c1 * (CS%s1(i,j,k) - CS%filter_alpha(k) * u(i,j)) + c2 * CS%u1(i,j,k) + enddo; enddo + enddo u1 => CS%u1 end subroutine Filt_accum -!> \namespace streaming_filter +!> \namespace mom_streaming_filter +!! +!! By Chengzhu Xu (chengzhu.xu@oregonstate.edu) and Edward D. Zaron +!! +!! The algorithm detects the instantaneous, narrowband tidal signals (u1) from the broadband +!! model output (u) by solving a set of coupled ODEs (the filter equations) at each time step. +!! In the filter equations, u1 is approximately the part of the signal that oscillates at the +!! filter's target frequency, and s1 is approximately the imaginary complement in time of u1. +!! +!! Major revision on Dec 9, 2024: The filters are no longer hard-coded. Instead, multiple filters +!! with tidal frequencies or arbitrary frequencies as their target frequencies can be turned on. +!! The filter names are specified in MOM_input and must consist of two letters/numbers. If the +!! name of a filter is the same as the name of a tidal constituent, then the corresponding tidal +!! frequency will be used as its target frequency. Otherwise, the user must specify the target +!! frequency. In either case, the target frequency is specified by "FILTER_${FILTER_NAME}_OMEGA". !! -!! This module detects instantaneous tidal signals in the model output using a set of coupled ODEs (the filter -!! equations), given the target frequency (om) and the bandwidth parameter (a) of the filter. At each timestep, -!! the filter takes model output (u) as the input and returns a time series consisting of sinusoidal motions (u1) -!! near its target frequency. The filtered tidal signals can be used to parameterize frequency-dependent drag, or -!! to detide the model output. See Xu & Zaron (2024) for detail. +!! The restarting capability has also been implemented. Because the filtering is a point-wise +!! operation, all variables are considered as fields, even if they are velocity components. !! -!! Reference: Xu, C., & Zaron, E. D. (2024). Detecting instantaneous tidal signals in ocean models utilizing -!! streaming band-pass filters. Journal of Advances in Modeling Earth Systems. Under review. +!! Xu, C., & Zaron, E. D. (2024). Detecting instantaneous tidal signals in ocean models utilizing +!! streaming band-pass filters. Journal of Advances in Modeling Earth Systems, 16, e2024MS004319. +!! https://doi.org/10.1029/2024MS004319 end module MOM_streaming_filter