From a798abcd12326f39b35fe60e72fe92cd769eb758 Mon Sep 17 00:00:00 2001 From: Raphael Dussin Date: Thu, 19 Nov 2020 16:38:31 -0500 Subject: [PATCH 01/10] consistent units + fix renormalization --- src/diagnostics/MOM_wave_structure.F90 | 52 ++++++++++++++++---------- 1 file changed, 32 insertions(+), 20 deletions(-) diff --git a/src/diagnostics/MOM_wave_structure.F90 b/src/diagnostics/MOM_wave_structure.F90 index 88b062472f..fd93294f06 100644 --- a/src/diagnostics/MOM_wave_structure.F90 +++ b/src/diagnostics/MOM_wave_structure.F90 @@ -11,6 +11,7 @@ module MOM_wave_structure ! MOM_wave_speed by Hallberg, 2008. use MOM_debugging, only : isnan => is_NaN +use MOM_checksums, only : chksum0, hchksum use MOM_diag_mediator, only : post_data, query_averaging_enabled, diag_ctrl use MOM_diag_mediator, only : register_diag_field, safe_alloc_ptr, time_type use MOM_EOS, only : calculate_density_derivs @@ -195,7 +196,13 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo Pi = (4.0*atan(1.0)) S => tv%S ; T => tv%T - g_Rho0 = US%L_T_to_m_s**2 * GV%g_Earth / GV%Rho0 + g_Rho0 = GV%g_Earth / GV%Rho0 + + if (CS%debug) call chksum0(g_Rho0, "g/rho0 in wave struct", & + scale=US%Z_to_m*(US%s_to_T**2)*US%kg_m3_to_R) + + if (CS%debug) call chksum0(freq, "freq in wave_struct", scale=US%s_to_T) + cg_subRO = 1e-100*US%m_s_to_L_T ! The hard-coded value here might need to increase. use_EOS = associated(tv%eqn_of_state) @@ -267,7 +274,9 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo !----------------------------------- if (G%mask2dT(i,j) > 0.5) then - lam = 1/(US%L_T_to_m_s**2 * cn(i,j)**2) + gprime(:) = 0.0 ! init gprime + pres(:) = 0.0 ! init pres + lam = 1/(cn(i,j)**2) ! Calculate drxh_sum if (use_EOS) then @@ -386,7 +395,7 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo do K=2,kc Igl(K) = 1.0/(gprime(K)*Hc(k)) ; Igu(K) = 1.0/(gprime(K)*Hc(k-1)) z_int(K) = z_int(K-1) + Hc(k-1) - N2(K) = US%m_to_Z**2*gprime(K)/(0.5*(Hc(k)+Hc(k-1))) + N2(K) = gprime(K)/(0.5*(Hc(k)+Hc(k-1))) enddo ! Set stratification for surface and bottom (setting equal to nearest interface for now) N2(1) = N2(2) ; N2(kc+1) = N2(kc) @@ -407,7 +416,7 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo ! Frist, populate interior rows do K=3,kc-1 row = K-1 ! indexing for TD matrix rows - gp_unscaled = US%m_to_Z*gprime(K) + gp_unscaled = gprime(K) lam_z(row) = lam*gp_unscaled a_diag(row) = gp_unscaled*(-Igu(K)) b_diag(row) = gp_unscaled*(Igu(K)+Igl(K)) - lam_z(row) @@ -419,14 +428,14 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo enddo ! Populate top row of tridiagonal matrix K=2 ; row = K-1 ; - gp_unscaled = US%m_to_Z*gprime(K) + gp_unscaled = gprime(K) lam_z(row) = lam*gp_unscaled a_diag(row) = 0.0 b_diag(row) = gp_unscaled*(Igu(K)+Igl(K)) - lam_z(row) c_diag(row) = gp_unscaled*(-Igl(K)) ! Populate bottom row of tridiagonal matrix K=kc ; row = K-1 - gp_unscaled = US%m_to_Z*gprime(K) + gp_unscaled = gprime(K) lam_z(row) = lam*gp_unscaled a_diag(row) = gp_unscaled*(-Igu(K)) b_diag(row) = gp_unscaled*(Igu(K)+Igl(K)) - lam_z(row) @@ -462,12 +471,11 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo !(including surface and bottom) w2avg = 0.0 do k=1,nzm-1 - dz(k) = US%Z_to_m*Hc(k) + dz(k) = Hc(k) w2avg = w2avg + 0.5*(w_strct(K)**2+w_strct(K+1)**2)*dz(k) enddo - !### Some mathematical cancellations could occur in the next two lines. - w2avg = w2avg / htot(i,j) - w_strct(:) = w_strct(:) / sqrt(htot(i,j)*w2avg*I_a_int) + ! correct renormalization: + w_strct(:) = w_strct(:) * sqrt(htot(i,j)*a_int/w2avg) ! Calculate vertical structure function of u (i.e. dw/dz) do K=2,nzm-1 @@ -478,10 +486,9 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo u_strct(nzm) = (w_strct(nzm-1)- w_strct(nzm))/dz(nzm-1) ! Calculate wavenumber magnitude - f2 = G%CoriolisBu(I,J)**2 - !f2 = 0.25*US%s_to_T**2 *((G%CoriolisBu(I,J)**2 + G%CoriolisBu(I-1,J-1)**2) + & - ! (G%CoriolisBu(I,J-1)**2 + G%CoriolisBu(I-1,J)**2)) - Kmag2 = US%m_to_L**2 * (freq**2 - f2) / (cn(i,j)**2 + cg_subRO**2) + f2 = (0.25*(G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1) + & + G%CoriolisBu(I,J-1) + G%CoriolisBu(I-1,J)))**2 + Kmag2 = (freq**2 - f2) / (cn(i,j)**2 + cg_subRO**2) ! Calculate terms in vertically integrated energy equation int_dwdz2 = 0.0 ; int_w2 = 0.0 ; int_N2w2 = 0.0 @@ -489,16 +496,16 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo w_strct2(1:nzm) = w_strct(1:nzm)**2 ! vertical integration with Trapezoidal rule do k=1,nzm-1 - int_dwdz2 = int_dwdz2 + 0.5*(u_strct2(K)+u_strct2(K+1)) * US%m_to_Z*dz(k) - int_w2 = int_w2 + 0.5*(w_strct2(K)+w_strct2(K+1)) * US%m_to_Z*dz(k) - int_N2w2 = int_N2w2 + 0.5*(w_strct2(K)*N2(K)+w_strct2(K+1)*N2(K+1)) * US%m_to_Z*dz(k) + int_dwdz2 = int_dwdz2 + 0.5*(u_strct2(K)+u_strct2(K+1)) * dz(k) + int_w2 = int_w2 + 0.5*(w_strct2(K)+w_strct2(K+1)) * dz(k) + int_N2w2 = int_N2w2 + 0.5*(w_strct2(K)*N2(K)+w_strct2(K+1)*N2(K+1)) * dz(k) enddo ! Back-calculate amplitude from energy equation if (present(En) .and. (freq**2*Kmag2 > 0.0)) then ! Units here are [R KE_term = 0.25*GV%Rho0*( ((freq**2 + f2) / (freq**2*Kmag2))*int_dwdz2 + int_w2 ) - PE_term = 0.25*GV%Rho0*( int_N2w2 / (US%s_to_T*freq)**2 ) + PE_term = 0.25*GV%Rho0*( int_N2w2 / freq**2 ) if (En(i,j) >= 0.0) then W0 = sqrt( En(i,j) / (KE_term + PE_term) ) else @@ -510,7 +517,7 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo W_profile(:) = W0*w_strct(:) ! dWdz_profile(:) = W0*u_strct(:) ! Calculate average magnitude of actual horizontal velocity over a period - Uavg_profile(:) = US%Z_to_L*abs(W0*u_strct(:)) * sqrt((freq**2 + f2) / (2.0*freq**2*Kmag2)) + Uavg_profile(:) = abs(W0*u_strct(:)) * sqrt((freq**2 + f2) / (2.0*freq**2*Kmag2)) else W_profile(:) = 0.0 ! dWdz_profile(:) = 0.0 @@ -522,7 +529,7 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo CS%u_strct(i,j,1:nzm) = u_strct(1:nzm) CS%W_profile(i,j,1:nzm) = W_profile(1:nzm) CS%Uavg_profile(i,j,1:nzm)= Uavg_profile(1:nzm) - CS%z_depths(i,j,1:nzm) = US%Z_to_m*z_int(1:nzm) + CS%z_depths(i,j,1:nzm) = z_int(1:nzm) CS%N2(i,j,1:nzm) = N2(1:nzm) CS%num_intfaces(i,j) = nzm else @@ -554,6 +561,11 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo endif ; enddo ! if cn>0.0? ; i-loop enddo ! j-loop + if (CS%debug) call hchksum(CS%N2, 'N2 in wave_struct', G%HI, scale=US%s_to_T**2) + if (CS%debug) call hchksum(cn, 'cn in wave_struct', G%HI, scale=US%L_T_to_m_s) + if (CS%debug) call hchksum(CS%W_profile, 'Wprofile in wave_struct', G%HI, scale=US%L_T_to_m_s) + if (CS%debug) call hchksum(CS%Uavg_profile, 'Uavg_profile in wave_struct', G%HI, scale=US%L_T_to_m_s) + end subroutine wave_structure !> Solves a tri-diagonal system Ax=y using either the standard From 1f1688baad40bd967372950acf5809f7697894d9 Mon Sep 17 00:00:00 2001 From: Raphael Dussin Date: Thu, 19 Nov 2020 16:54:01 -0500 Subject: [PATCH 02/10] add debug --- src/diagnostics/MOM_wave_structure.F90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/diagnostics/MOM_wave_structure.F90 b/src/diagnostics/MOM_wave_structure.F90 index fd93294f06..1b38ce151d 100644 --- a/src/diagnostics/MOM_wave_structure.F90 +++ b/src/diagnostics/MOM_wave_structure.F90 @@ -58,6 +58,7 @@ module MOM_wave_structure !! for internal tide for testing (BDM) real :: int_tide_source_y !< Y Location of generation site !! for internal tide for testing (BDM) + logical :: debug !! debugging prints end type wave_structure_CS @@ -721,6 +722,8 @@ subroutine wave_structure_init(Time, G, param_file, diag, CS) "X Location of generation site for internal tide", default=1.) call get_param(param_file, mdl, "INTERNAL_TIDE_SOURCE_Y", CS%int_tide_source_y, & "Y Location of generation site for internal tide", default=1.) + call get_param(param_file, mdl, "DEBUG", CS%debug, & + "debugging prints", default=.false.) CS%diag => diag From 220536b7a3c30ff49d5b15ffd8556032a08a10a3 Mon Sep 17 00:00:00 2001 From: Raphael Dussin Date: Thu, 19 Nov 2020 16:57:25 -0500 Subject: [PATCH 03/10] fix doxygen error --- src/diagnostics/MOM_wave_structure.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/diagnostics/MOM_wave_structure.F90 b/src/diagnostics/MOM_wave_structure.F90 index 1b38ce151d..0a8545334e 100644 --- a/src/diagnostics/MOM_wave_structure.F90 +++ b/src/diagnostics/MOM_wave_structure.F90 @@ -58,7 +58,7 @@ module MOM_wave_structure !! for internal tide for testing (BDM) real :: int_tide_source_y !< Y Location of generation site !! for internal tide for testing (BDM) - logical :: debug !! debugging prints + logical :: debug !< debugging prints end type wave_structure_CS From 36e18819076f6716884a1d17a962367d8d623eee Mon Sep 17 00:00:00 2001 From: Raphael Dussin Date: Thu, 19 Nov 2020 17:52:50 -0500 Subject: [PATCH 04/10] remove whitespace --- src/diagnostics/MOM_wave_structure.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/diagnostics/MOM_wave_structure.F90 b/src/diagnostics/MOM_wave_structure.F90 index 0a8545334e..c290465315 100644 --- a/src/diagnostics/MOM_wave_structure.F90 +++ b/src/diagnostics/MOM_wave_structure.F90 @@ -566,7 +566,7 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo if (CS%debug) call hchksum(cn, 'cn in wave_struct', G%HI, scale=US%L_T_to_m_s) if (CS%debug) call hchksum(CS%W_profile, 'Wprofile in wave_struct', G%HI, scale=US%L_T_to_m_s) if (CS%debug) call hchksum(CS%Uavg_profile, 'Uavg_profile in wave_struct', G%HI, scale=US%L_T_to_m_s) - + end subroutine wave_structure !> Solves a tri-diagonal system Ax=y using either the standard From 83785ae706e4ad3af67133a293130f1d869ba01e Mon Sep 17 00:00:00 2001 From: Raphael Dussin Date: Mon, 23 Nov 2020 21:22:49 -0500 Subject: [PATCH 05/10] enable avg diags + more dim fixes --- src/diagnostics/MOM_wave_structure.F90 | 74 +++++++++++-------- .../lateral/MOM_internal_tides.F90 | 8 ++ .../vertical/MOM_internal_tide_input.F90 | 11 ++- 3 files changed, 61 insertions(+), 32 deletions(-) diff --git a/src/diagnostics/MOM_wave_structure.F90 b/src/diagnostics/MOM_wave_structure.F90 index c290465315..b0c543f12f 100644 --- a/src/diagnostics/MOM_wave_structure.F90 +++ b/src/diagnostics/MOM_wave_structure.F90 @@ -38,9 +38,9 @@ module MOM_wave_structure type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to !! regulate the timing of diagnostic output. real, allocatable, dimension(:,:,:) :: w_strct - !< Vertical structure of vertical velocity (normalized) [m s-1]. + !< Vertical structure of vertical velocity (normalized) [nondim]. real, allocatable, dimension(:,:,:) :: u_strct - !< Vertical structure of horizontal velocity (normalized) [m s-1]. + !< Vertical structure of horizontal velocity (normalized) [nondim]. real, allocatable, dimension(:,:,:) :: W_profile !< Vertical profile of w_hat(z), where !! w(x,y,z,t) = w_hat(z)*exp(i(kx+ly-freq*t)) is the full time- @@ -49,11 +49,11 @@ module MOM_wave_structure !< Vertical profile of the magnitude of horizontal velocity, !! (u^2+v^2)^0.5, averaged over a period [L T-1 ~> m s-1]. real, allocatable, dimension(:,:,:) :: z_depths - !< Depths of layer interfaces [m]. + !< Depths of layer interfaces [Z ~> m]. real, allocatable, dimension(:,:,:) :: N2 - !< Squared buoyancy frequency at each interface [s-2]. + !< Squared buoyancy frequency at each interface [T-2 ~> s-2]. integer, allocatable, dimension(:,:):: num_intfaces - !< Number of layer interfaces (including surface and bottom) + !< Number of layer interfaces (including surface and bottom) [nondim]. real :: int_tide_source_x !< X Location of generation site !! for internal tide for testing (BDM) real :: int_tide_source_y !< Y Location of generation site @@ -111,13 +111,13 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo real, dimension(SZK_(G)+1) :: & dRho_dT, & ! Partial derivative of density with temperature [R degC-1 ~> kg m-3 degC-1] dRho_dS, & ! Partial derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1] - pres, & ! Interface pressure [R L2 T-2 ~> Pa] + pres, & ! Interface pressure [R L H T-2 ~> Pa] T_int, & ! Temperature interpolated to interfaces [degC] S_int, & ! Salinity interpolated to interfaces [ppt] - gprime ! The reduced gravity across each interface [m2 Z-1 s-2 ~> m s-2]. + gprime ! The reduced gravity across each interface [L2 Z-1 T-2 ~> m s-2]. real, dimension(SZK_(G)) :: & Igl, Igu ! The inverse of the reduced gravity across an interface times - ! the thickness of the layer below (Igl) or above (Igu) it [s2 m-2]. + ! the thickness of the layer below (Igl) or above (Igu) it [T2 L-2 ~> s2 m-2]. real, dimension(SZK_(G),SZI_(G)) :: & Hf, & ! Layer thicknesses after very thin layers are combined [Z ~> m] Tf, & ! Layer temperatures after very thin layers are combined [degC] @@ -130,10 +130,10 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo Rc, & ! A column of layer densities after convective istabilities are removed [R ~> kg m-3] det, ddet real, dimension(SZI_(G),SZJ_(G)) :: & - htot ! The vertical sum of the thicknesses [Z ~> m] - real :: lam - real :: min_h_frac - real :: Z_to_pres ! A conversion factor from thicknesses to pressure [R L2 T-2 Z-1 ~> Pa m-1] + htot ! The vertical sum of the thicknesses [Z ~> m] + real :: lam ! inverse of wave speed squared [T2 L-2 ~> s2 m-2] + real :: min_h_frac ! fractional (per layer) minimum thickness [nondim] + real :: Z_to_pres ! A conversion factor from thicknesses to pressure [R L2 T-2 Z-1 ~> Pa m-1] real, dimension(SZI_(G)) :: & hmin, & ! Thicknesses [Z ~> m] H_here, & ! A thickness [Z ~> m] @@ -145,32 +145,38 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo real :: drxh_sum ! The sum of density diffrences across interfaces times thicknesses [R Z ~> kg m-2] real, parameter :: tol1 = 0.0001, tol2 = 0.001 real, pointer, dimension(:,:,:) :: T => NULL(), S => NULL() - real :: g_Rho0 ! G_Earth/Rho0 in [m2 s-2 Z-1 R-1 ~> m4 s-2 kg-1]. + real :: g_Rho0 ! G_Earth/Rho0 in [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1]. ! real :: rescale, I_rescale integer :: kf(SZI_(G)) integer, parameter :: max_itt = 1 ! number of times to iterate in solving for eigenvector real :: cg_subRO ! A tiny wave speed to prevent division by zero [L T-1 ~> m s-1] - real, parameter :: a_int = 0.5 ! value of normalized integral: \int(w_strct^2)dz = a_int - real :: I_a_int ! inverse of a_int + real, parameter :: a_int = 0.5 ! value of normalized integral: \int(w_strct^2)dz = a_int [nondim] + real :: I_a_int ! inverse of a_int [nondim] real :: f2 ! squared Coriolis frequency [T-2 ~> s-2] - real :: Kmag2 ! magnitude of horizontal wave number squared + real :: Kmag2 ! magnitude of horizontal wave number squared [L-2 ~> m-2] logical :: use_EOS ! If true, density is calculated from T & S using an ! equation of state. - real, dimension(SZK_(G)+1) :: w_strct, u_strct, W_profile, Uavg_profile, z_int, N2 - ! local representations of variables in CS; note, - ! not all rows will be filled if layers get merged! - real, dimension(SZK_(G)+1) :: w_strct2, u_strct2 - ! squared values - real, dimension(SZK_(G)) :: dz ! thicknesses of merged layers (same as Hc I hope) + + ! local representations of variables in CS; note, + ! not all rows will be filled if layers get merged! + real, dimension(SZK_(G)+1) :: w_strct ! Vertical structure of vertical velocity (normalized) [nondim]. + real, dimension(SZK_(G)+1) :: u_strct ! Vertical structure of horizontal velocity (normalized) [nondim]. + real, dimension(SZK_(G)+1) :: W_profile ! Vertical profile of w_hat(z) = W0*w_strct(z) [Z T-1 ~> m s-1]. + real, dimension(SZK_(G)+1) :: Uavg_profile ! Vertical profile of the magnitude of horizontal velocity [L T-1 ~> m s-1]. + real, dimension(SZK_(G)+1) :: z_int ! Integrated depth [Z ~> m] + real, dimension(SZK_(G)+1) :: N2 ! Squared buoyancy frequency at each interface [T-2 ~> s-2]. + real, dimension(SZK_(G)+1) :: w_strct2 ! squared values [nondim] + real, dimension(SZK_(G)+1) :: u_strct2 ! squared values [nondim] + real, dimension(SZK_(G)) :: dz ! thicknesses of merged layers (same as Hc I hope) [Z ~> m] ! real, dimension(SZK_(G)+1) :: dWdz_profile ! profile of dW/dz - real :: w2avg ! average of squared vertical velocity structure funtion + real :: w2avg ! average of squared vertical velocity structure funtion [Z ~> m] real :: int_dwdz2 real :: int_w2 real :: int_N2w2 real :: KE_term ! terms in vertically averaged energy equation real :: PE_term ! terms in vertically averaged energy equation real :: W0 ! A vertical velocity magnitude [Z T-1 ~> m s-1] - real :: gp_unscaled ! A version of gprime rescaled to [m s-2]. + real :: gp_unscaled ! A version of gprime rescaled to [L T-2 ~> m s-2]. real, dimension(SZK_(G)-1) :: lam_z ! product of eigen value and gprime(k); one value for each ! interface (excluding surface and bottom) real, dimension(SZK_(G)-1) :: a_diag, b_diag, c_diag @@ -199,8 +205,8 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo S => tv%S ; T => tv%T g_Rho0 = GV%g_Earth / GV%Rho0 - if (CS%debug) call chksum0(g_Rho0, "g/rho0 in wave struct", & - scale=US%Z_to_m*(US%s_to_T**2)*US%kg_m3_to_R) + !if (CS%debug) call chksum0(g_Rho0, "g/rho0 in wave struct", & + ! scale=(US%L_to_m**2)*US%m_to_Z*(US%s_to_T**2)*US%kg_m3_to_R) if (CS%debug) call chksum0(freq, "freq in wave_struct", scale=US%s_to_T) @@ -396,7 +402,7 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo do K=2,kc Igl(K) = 1.0/(gprime(K)*Hc(k)) ; Igu(K) = 1.0/(gprime(K)*Hc(k-1)) z_int(K) = z_int(K-1) + Hc(k-1) - N2(K) = gprime(K)/(0.5*(Hc(k)+Hc(k-1))) + N2(K) = US%L_to_Z**2*gprime(K)/(0.5*(Hc(k)+Hc(k-1))) enddo ! Set stratification for surface and bottom (setting equal to nearest interface for now) N2(1) = N2(2) ; N2(kc+1) = N2(kc) @@ -415,6 +421,16 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo ! [-1/H(k-1)]e(k-1) + [1/H(k-1)+1/H(k)-lam_z]e(k) + [-1/H(k)]e(k+1) = 0, ! where lam_z = lam*gprime is now a function of depth. ! Frist, populate interior rows + + ! init the values in matrix: since number of layers is variable, values need + ! to be reset + lam_z(:) = 0.0 + a_diag(:) = 0.0 + b_diag(:) = 0.0 + c_diag(:) = 0.0 + e_guess(:) = 0.0 + e_itt(:) = 0.0 + w_strct(:) = 0.0 do K=3,kc-1 row = K-1 ! indexing for TD matrix rows gp_unscaled = gprime(K) @@ -564,7 +580,7 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo if (CS%debug) call hchksum(CS%N2, 'N2 in wave_struct', G%HI, scale=US%s_to_T**2) if (CS%debug) call hchksum(cn, 'cn in wave_struct', G%HI, scale=US%L_T_to_m_s) - if (CS%debug) call hchksum(CS%W_profile, 'Wprofile in wave_struct', G%HI, scale=US%L_T_to_m_s) + if (CS%debug) call hchksum(CS%W_profile, 'Wprofile in wave_struct', G%HI, scale=US%Z_to_L*US%L_T_to_m_s) if (CS%debug) call hchksum(CS%Uavg_profile, 'Uavg_profile in wave_struct', G%HI, scale=US%L_T_to_m_s) end subroutine wave_structure @@ -653,7 +669,7 @@ subroutine tridiag_solver(a, b, c, h, y, method, x) ! Need to add a check for these conditions. do k=1,nrow-1 if (abs(a(k+1)-c(k)) > 1.e-10*(abs(a(k+1))+abs(c(k)))) then - call MOM_error(WARNING, "tridiag_solver: matrix not symmetric; need symmetry when invoking TDMA_H") + call MOM_error(FATAL, "tridiag_solver: matrix not symmetric; need symmetry when invoking TDMA_H") endif enddo alpha = -c diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index 2bb3c3b0f1..5d0cce3cd8 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -7,6 +7,7 @@ module MOM_internal_tides use MOM_debugging, only : is_NaN use MOM_diag_mediator, only : post_data, query_averaging_enabled, diag_axis_init +use MOM_diag_mediator, only : disable_averaging, enable_averages use MOM_diag_mediator, only : register_diag_field, diag_ctrl, safe_alloc_ptr use MOM_diag_mediator, only : axes_grp, define_axes_group use MOM_domains, only : AGRID, To_South, To_West, To_All @@ -202,6 +203,8 @@ subroutine propagate_int_tide(h, tv, cn, TKE_itidal_input, vel_btTide, Nb, dt, & integer :: a, m, fr, i, j, is, ie, js, je, isd, ied, jsd, jed, nAngle, nzm integer :: id_g, jd_g ! global (decomp-invar) indices (for debugging) type(group_pass_type), save :: pass_test, pass_En + type(time_type) :: time_end + logical:: avg_enabled if (.not.associated(CS)) return is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec @@ -497,6 +500,9 @@ subroutine propagate_int_tide(h, tv, cn, TKE_itidal_input, vel_btTide, Nb, dt, & enddo ; enddo ! Output diagnostics.************************************************************ + avg_enabled = query_averaging_enabled(CS%diag, time_end=time_end) + call enable_averages(dt, time_end, CS%diag) + if (query_averaging_enabled(CS%diag)) then ! Output two-dimensional diagnostistics if (CS%id_tot_En > 0) call post_data(CS%id_tot_En, tot_En, CS%diag) @@ -587,6 +593,8 @@ subroutine propagate_int_tide(h, tv, cn, TKE_itidal_input, vel_btTide, Nb, dt, & endif + call disable_averaging(CS%diag) + end subroutine propagate_int_tide !> Checks for energy conservation on computational domain diff --git a/src/parameterizations/vertical/MOM_internal_tide_input.F90 b/src/parameterizations/vertical/MOM_internal_tide_input.F90 index f5b9e7dbb7..059f5b9ab6 100644 --- a/src/parameterizations/vertical/MOM_internal_tide_input.F90 +++ b/src/parameterizations/vertical/MOM_internal_tide_input.F90 @@ -6,6 +6,7 @@ module MOM_int_tide_input use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE use MOM_diag_mediator, only : diag_ctrl, query_averaging_enabled +use MOM_diag_mediator, only : disable_averaging, enable_averages use MOM_diag_mediator, only : safe_alloc_ptr, post_data, register_diag_field use MOM_debugging, only : hchksum use MOM_error_handler, only : MOM_error, is_root_pe, FATAL, WARNING, NOTE @@ -55,7 +56,7 @@ module MOM_int_tide_input !>@{ Diagnostic IDs - integer :: id_TKE_itidal = -1, id_Nb = -1, id_N2_bot = -1 + integer :: id_TKE_itidal_itide = -1, id_Nb = -1, id_N2_bot = -1 !>@} end type int_tide_input_CS @@ -140,10 +141,14 @@ subroutine set_int_tide_input(u, v, h, tv, fluxes, itide, dt, G, GV, US, CS) scale=US%RZ3_T3_to_W_m2) endif - if (CS%id_TKE_itidal > 0) call post_data(CS%id_TKE_itidal, itide%TKE_itidal_input, CS%diag) + call enable_averages(dt, time_end, CS%diag) + + if (CS%id_TKE_itidal_itide > 0) call post_data(CS%id_TKE_itidal_itide, itide%TKE_itidal_input, CS%diag) if (CS%id_Nb > 0) call post_data(CS%id_Nb, itide%Nb, CS%diag) if (CS%id_N2_bot > 0 ) call post_data(CS%id_N2_bot, N2_bot, CS%diag) + call disable_averaging(CS%diag) + end subroutine set_int_tide_input !> Estimates the near-bottom buoyancy frequency (N^2). @@ -409,7 +414,7 @@ subroutine int_tide_input_init(Time, G, GV, US, param_file, diag, CS, itide) enddo ; enddo - CS%id_TKE_itidal = register_diag_field('ocean_model','TKE_itidal_itide',diag%axesT1,Time, & + CS%id_TKE_itidal_itide = register_diag_field('ocean_model','TKE_itidal_itide',diag%axesT1,Time, & 'Internal Tide Driven Turbulent Kinetic Energy', & 'W m-2', conversion=US%RZ3_T3_to_W_m2) From 70622895b6b71f1bea25bc0201fbae3ad38cf91b Mon Sep 17 00:00:00 2001 From: Raphael Dussin Date: Wed, 9 Dec 2020 12:59:09 -0500 Subject: [PATCH 06/10] switch to different tridiag solver --- src/ALE/regrid_edge_values.F90 | 1 + src/diagnostics/MOM_wave_structure.F90 | 10 ++++++++-- 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/src/ALE/regrid_edge_values.F90 b/src/ALE/regrid_edge_values.F90 index 46570b26b9..d9c25d8cd1 100644 --- a/src/ALE/regrid_edge_values.F90 +++ b/src/ALE/regrid_edge_values.F90 @@ -16,6 +16,7 @@ module regrid_edge_values public edge_values_explicit_h2, edge_values_explicit_h4 public edge_values_implicit_h4, edge_values_implicit_h6 public edge_slopes_implicit_h3, edge_slopes_implicit_h5 +public solve_diag_dominant_tridiag ! public solve_diag_dominant_tridiag, linear_solver ! The following parameters are used to avoid singular matrices for boundary diff --git a/src/diagnostics/MOM_wave_structure.F90 b/src/diagnostics/MOM_wave_structure.F90 index b0c543f12f..35fbe33218 100644 --- a/src/diagnostics/MOM_wave_structure.F90 +++ b/src/diagnostics/MOM_wave_structure.F90 @@ -21,6 +21,7 @@ module MOM_wave_structure use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type +use regrid_edge_values, only : solve_diag_dominant_tridiag implicit none ; private @@ -464,8 +465,13 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo ! Perform inverse iteration with tri-diag solver do itt=1,max_itt - call tridiag_solver(a_diag(1:kc-1),b_diag(1:kc-1),c_diag(1:kc-1), & - -lam_z(1:kc-1),e_guess(1:kc-1),"TDMA_H",e_itt) + ! this solver becomes unstable very quickly + !call tridiag_solver(a_diag(1:kc-1),b_diag(1:kc-1),c_diag(1:kc-1), & + ! -lam_z(1:kc-1),e_guess(1:kc-1),"TDMA_T",e_itt) + + call solve_diag_dominant_tridiag( a_diag(1:kc-1), -lam_z(1:kc-1), & + c_diag(1:kc-1), e_guess(1:kc-1), & + e_itt, kc-1 ) e_guess(1:kc-1) = e_itt(1:kc-1) / sqrt(sum(e_itt(1:kc-1)**2)) enddo ! itt-loop w_strct(2:kc) = e_guess(1:kc-1) From 723bfdbe81e1f1af167a4d1827da9ad0d79dd7e7 Mon Sep 17 00:00:00 2001 From: Raphael Dussin Date: Wed, 9 Dec 2020 13:28:59 -0500 Subject: [PATCH 07/10] fix comments doxygen --- src/diagnostics/MOM_wave_structure.F90 | 114 ++++++++++++------------- 1 file changed, 57 insertions(+), 57 deletions(-) diff --git a/src/diagnostics/MOM_wave_structure.F90 b/src/diagnostics/MOM_wave_structure.F90 index 35fbe33218..44a00095cc 100644 --- a/src/diagnostics/MOM_wave_structure.F90 +++ b/src/diagnostics/MOM_wave_structure.F90 @@ -110,81 +110,81 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo !! over the entire computational domain. ! Local variables real, dimension(SZK_(G)+1) :: & - dRho_dT, & ! Partial derivative of density with temperature [R degC-1 ~> kg m-3 degC-1] - dRho_dS, & ! Partial derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1] - pres, & ! Interface pressure [R L H T-2 ~> Pa] - T_int, & ! Temperature interpolated to interfaces [degC] - S_int, & ! Salinity interpolated to interfaces [ppt] - gprime ! The reduced gravity across each interface [L2 Z-1 T-2 ~> m s-2]. + dRho_dT, & !< Partial derivative of density with temperature [R degC-1 ~> kg m-3 degC-1] + dRho_dS, & !< Partial derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1] + pres, & !< Interface pressure [R L H T-2 ~> Pa] + T_int, & !< Temperature interpolated to interfaces [degC] + S_int, & !< Salinity interpolated to interfaces [ppt] + gprime !< The reduced gravity across each interface [L2 Z-1 T-2 ~> m s-2]. real, dimension(SZK_(G)) :: & - Igl, Igu ! The inverse of the reduced gravity across an interface times - ! the thickness of the layer below (Igl) or above (Igu) it [T2 L-2 ~> s2 m-2]. + Igl, Igu !< The inverse of the reduced gravity across an interface times + !< the thickness of the layer below (Igl) or above (Igu) it [T2 L-2 ~> s2 m-2]. real, dimension(SZK_(G),SZI_(G)) :: & - Hf, & ! Layer thicknesses after very thin layers are combined [Z ~> m] - Tf, & ! Layer temperatures after very thin layers are combined [degC] - Sf, & ! Layer salinities after very thin layers are combined [ppt] - Rf ! Layer densities after very thin layers are combined [R ~> kg m-3] + Hf, & !< Layer thicknesses after very thin layers are combined [Z ~> m] + Tf, & !< Layer temperatures after very thin layers are combined [degC] + Sf, & !< Layer salinities after very thin layers are combined [ppt] + Rf !< Layer densities after very thin layers are combined [R ~> kg m-3] real, dimension(SZK_(G)) :: & - Hc, & ! A column of layer thicknesses after convective istabilities are removed [Z ~> m] - Tc, & ! A column of layer temperatures after convective istabilities are removed [degC] - Sc, & ! A column of layer salinites after convective istabilities are removed [ppt] - Rc, & ! A column of layer densities after convective istabilities are removed [R ~> kg m-3] + Hc, & !< A column of layer thicknesses after convective istabilities are removed [Z ~> m] + Tc, & !< A column of layer temperatures after convective istabilities are removed [degC] + Sc, & !< A column of layer salinites after convective istabilities are removed [ppt] + Rc, & !< A column of layer densities after convective istabilities are removed [R ~> kg m-3] det, ddet real, dimension(SZI_(G),SZJ_(G)) :: & - htot ! The vertical sum of the thicknesses [Z ~> m] - real :: lam ! inverse of wave speed squared [T2 L-2 ~> s2 m-2] - real :: min_h_frac ! fractional (per layer) minimum thickness [nondim] - real :: Z_to_pres ! A conversion factor from thicknesses to pressure [R L2 T-2 Z-1 ~> Pa m-1] + htot !< The vertical sum of the thicknesses [Z ~> m] + real :: lam !< inverse of wave speed squared [T2 L-2 ~> s2 m-2] + real :: min_h_frac !< fractional (per layer) minimum thickness [nondim] + real :: Z_to_pres !< A conversion factor from thicknesses to pressure [R L2 T-2 Z-1 ~> Pa m-1] real, dimension(SZI_(G)) :: & - hmin, & ! Thicknesses [Z ~> m] - H_here, & ! A thickness [Z ~> m] - HxT_here, & ! A layer integrated temperature [degC Z ~> degC m] - HxS_here, & ! A layer integrated salinity [ppt Z ~> ppt m] - HxR_here ! A layer integrated density [R Z ~> kg m-2] + hmin, & !< Thicknesses [Z ~> m] + H_here, & !< A thickness [Z ~> m] + HxT_here, & !< A layer integrated temperature [degC Z ~> degC m] + HxS_here, & !< A layer integrated salinity [ppt Z ~> ppt m] + HxR_here !< A layer integrated density [R Z ~> kg m-2] real :: speed2_tot - real :: I_Hnew ! The inverse of a new layer thickness [Z-1 ~> m-1] - real :: drxh_sum ! The sum of density diffrences across interfaces times thicknesses [R Z ~> kg m-2] + real :: I_Hnew !< The inverse of a new layer thickness [Z-1 ~> m-1] + real :: drxh_sum !< The sum of density diffrences across interfaces times thicknesses [R Z ~> kg m-2] real, parameter :: tol1 = 0.0001, tol2 = 0.001 real, pointer, dimension(:,:,:) :: T => NULL(), S => NULL() - real :: g_Rho0 ! G_Earth/Rho0 in [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1]. + real :: g_Rho0 !< G_Earth/Rho0 in [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1]. ! real :: rescale, I_rescale integer :: kf(SZI_(G)) - integer, parameter :: max_itt = 1 ! number of times to iterate in solving for eigenvector - real :: cg_subRO ! A tiny wave speed to prevent division by zero [L T-1 ~> m s-1] - real, parameter :: a_int = 0.5 ! value of normalized integral: \int(w_strct^2)dz = a_int [nondim] - real :: I_a_int ! inverse of a_int [nondim] - real :: f2 ! squared Coriolis frequency [T-2 ~> s-2] - real :: Kmag2 ! magnitude of horizontal wave number squared [L-2 ~> m-2] - logical :: use_EOS ! If true, density is calculated from T & S using an - ! equation of state. + integer, parameter :: max_itt = 1 !< number of times to iterate in solving for eigenvector + real :: cg_subRO !< A tiny wave speed to prevent division by zero [L T-1 ~> m s-1] + real, parameter :: a_int = 0.5 !< value of normalized integral: \int(w_strct^2)dz = a_int [nondim] + real :: I_a_int !< inverse of a_int [nondim] + real :: f2 !< squared Coriolis frequency [T-2 ~> s-2] + real :: Kmag2 !< magnitude of horizontal wave number squared [L-2 ~> m-2] + logical :: use_EOS !< If true, density is calculated from T & S using an + !! equation of state. ! local representations of variables in CS; note, ! not all rows will be filled if layers get merged! - real, dimension(SZK_(G)+1) :: w_strct ! Vertical structure of vertical velocity (normalized) [nondim]. - real, dimension(SZK_(G)+1) :: u_strct ! Vertical structure of horizontal velocity (normalized) [nondim]. - real, dimension(SZK_(G)+1) :: W_profile ! Vertical profile of w_hat(z) = W0*w_strct(z) [Z T-1 ~> m s-1]. - real, dimension(SZK_(G)+1) :: Uavg_profile ! Vertical profile of the magnitude of horizontal velocity [L T-1 ~> m s-1]. - real, dimension(SZK_(G)+1) :: z_int ! Integrated depth [Z ~> m] - real, dimension(SZK_(G)+1) :: N2 ! Squared buoyancy frequency at each interface [T-2 ~> s-2]. - real, dimension(SZK_(G)+1) :: w_strct2 ! squared values [nondim] - real, dimension(SZK_(G)+1) :: u_strct2 ! squared values [nondim] - real, dimension(SZK_(G)) :: dz ! thicknesses of merged layers (same as Hc I hope) [Z ~> m] - ! real, dimension(SZK_(G)+1) :: dWdz_profile ! profile of dW/dz - real :: w2avg ! average of squared vertical velocity structure funtion [Z ~> m] + real, dimension(SZK_(G)+1) :: w_strct !< Vertical structure of vertical velocity (normalized) [nondim]. + real, dimension(SZK_(G)+1) :: u_strct !< Vertical structure of horizontal velocity (normalized) [nondim]. + real, dimension(SZK_(G)+1) :: W_profile !< Vertical profile of w_hat(z) = W0*w_strct(z) [Z T-1 ~> m s-1]. + real, dimension(SZK_(G)+1) :: Uavg_profile !< Vertical profile of the magnitude of horizontal velocity [L T-1 ~> m s-1]. + real, dimension(SZK_(G)+1) :: z_int !< Integrated depth [Z ~> m] + real, dimension(SZK_(G)+1) :: N2 !< Squared buoyancy frequency at each interface [T-2 ~> s-2]. + real, dimension(SZK_(G)+1) :: w_strct2 !< squared values [nondim] + real, dimension(SZK_(G)+1) :: u_strct2 !< squared values [nondim] + real, dimension(SZK_(G)) :: dz !< thicknesses of merged layers (same as Hc I hope) [Z ~> m] + ! real, dimension(SZK_(G)+1) :: dWdz_profile !< profile of dW/dz + real :: w2avg !< average of squared vertical velocity structure funtion [Z ~> m] real :: int_dwdz2 real :: int_w2 real :: int_N2w2 - real :: KE_term ! terms in vertically averaged energy equation - real :: PE_term ! terms in vertically averaged energy equation - real :: W0 ! A vertical velocity magnitude [Z T-1 ~> m s-1] - real :: gp_unscaled ! A version of gprime rescaled to [L T-2 ~> m s-2]. - real, dimension(SZK_(G)-1) :: lam_z ! product of eigen value and gprime(k); one value for each - ! interface (excluding surface and bottom) + real :: KE_term !< terms in vertically averaged energy equation + real :: PE_term !< terms in vertically averaged energy equation + real :: W0 !< A vertical velocity magnitude [Z T-1 ~> m s-1] + real :: gp_unscaled !< A version of gprime rescaled to [L T-2 ~> m s-2]. + real, dimension(SZK_(G)-1) :: lam_z !< product of eigen value and gprime(k); one value for each + !< interface (excluding surface and bottom) real, dimension(SZK_(G)-1) :: a_diag, b_diag, c_diag - ! diagonals of tridiagonal matrix; one value for each - ! interface (excluding surface and bottom) - real, dimension(SZK_(G)-1) :: e_guess ! guess at eigen vector with unit amplitde (for TDMA) - real, dimension(SZK_(G)-1) :: e_itt ! improved guess at eigen vector (from TDMA) + !< diagonals of tridiagonal matrix; one value for each + !< interface (excluding surface and bottom) + real, dimension(SZK_(G)-1) :: e_guess !< guess at eigen vector with unit amplitde (for TDMA) + real, dimension(SZK_(G)-1) :: e_itt !< improved guess at eigen vector (from TDMA) real :: Pi integer :: kc integer :: i, j, k, k2, itt, is, ie, js, je, nz, nzm, row, ig, jg, ig_stop, jg_stop From 37e799ebc17fd438522dd409f8bc32a7c5e9fd63 Mon Sep 17 00:00:00 2001 From: Raphael Dussin Date: Fri, 11 Dec 2020 11:10:23 -0500 Subject: [PATCH 08/10] split comment --- src/diagnostics/MOM_wave_structure.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/diagnostics/MOM_wave_structure.F90 b/src/diagnostics/MOM_wave_structure.F90 index 44a00095cc..b647731460 100644 --- a/src/diagnostics/MOM_wave_structure.F90 +++ b/src/diagnostics/MOM_wave_structure.F90 @@ -163,7 +163,8 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo real, dimension(SZK_(G)+1) :: w_strct !< Vertical structure of vertical velocity (normalized) [nondim]. real, dimension(SZK_(G)+1) :: u_strct !< Vertical structure of horizontal velocity (normalized) [nondim]. real, dimension(SZK_(G)+1) :: W_profile !< Vertical profile of w_hat(z) = W0*w_strct(z) [Z T-1 ~> m s-1]. - real, dimension(SZK_(G)+1) :: Uavg_profile !< Vertical profile of the magnitude of horizontal velocity [L T-1 ~> m s-1]. + real, dimension(SZK_(G)+1) :: Uavg_profile !< Vertical profile of the magnitude of + !! horizontal velocity [L T-1 ~> m s-1]. real, dimension(SZK_(G)+1) :: z_int !< Integrated depth [Z ~> m] real, dimension(SZK_(G)+1) :: N2 !< Squared buoyancy frequency at each interface [T-2 ~> s-2]. real, dimension(SZK_(G)+1) :: w_strct2 !< squared values [nondim] From 54d201fcbbf6d8dfa9ad8d16a53f46a6fb27b459 Mon Sep 17 00:00:00 2001 From: Raphael Dussin Date: Mon, 14 Dec 2020 18:02:14 -0500 Subject: [PATCH 09/10] fix out of bounds with full halos --- src/diagnostics/MOM_wave_structure.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/diagnostics/MOM_wave_structure.F90 b/src/diagnostics/MOM_wave_structure.F90 index b647731460..05494fb819 100644 --- a/src/diagnostics/MOM_wave_structure.F90 +++ b/src/diagnostics/MOM_wave_structure.F90 @@ -510,8 +510,8 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo u_strct(nzm) = (w_strct(nzm-1)- w_strct(nzm))/dz(nzm-1) ! Calculate wavenumber magnitude - f2 = (0.25*(G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1) + & - G%CoriolisBu(I,J-1) + G%CoriolisBu(I-1,J)))**2 + f2 = (0.25*(G%CoriolisBu(I,J) + G%CoriolisBu(max(I-1,1),max(J-1,1)) + & + G%CoriolisBu(I,max(J-1,1)) + G%CoriolisBu(max(I-1,1),J)))**2 Kmag2 = (freq**2 - f2) / (cn(i,j)**2 + cg_subRO**2) ! Calculate terms in vertically integrated energy equation From 8426fe7ab07f57a53e3f0e882c24a774ed8beb81 Mon Sep 17 00:00:00 2001 From: Raphael Dussin Date: Mon, 14 Dec 2020 18:04:14 -0500 Subject: [PATCH 10/10] fix diag write error --- src/parameterizations/vertical/MOM_internal_tide_input.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/parameterizations/vertical/MOM_internal_tide_input.F90 b/src/parameterizations/vertical/MOM_internal_tide_input.F90 index 059f5b9ab6..79c69be095 100644 --- a/src/parameterizations/vertical/MOM_internal_tide_input.F90 +++ b/src/parameterizations/vertical/MOM_internal_tide_input.F90 @@ -115,6 +115,8 @@ subroutine set_int_tide_input(u, v, h, tv, fluxes, itide, dt, G, GV, US, CS) call find_N2_bottom(h, tv, T_f, S_f, itide%h2, fluxes, G, GV, US, N2_bot) + avg_enabled = query_averaging_enabled(CS%diag, time_end=time_end) + !$OMP parallel do default(shared) do j=js,je ; do i=is,ie itide%Nb(i,j) = G%mask2dT(i,j) * sqrt(N2_bot(i,j)) @@ -123,7 +125,6 @@ subroutine set_int_tide_input(u, v, h, tv, fluxes, itide, dt, G, GV, US, CS) if (CS%int_tide_source_test) then itide%TKE_itidal_input(:,:) = 0.0 - avg_enabled = query_averaging_enabled(CS%diag, time_end=time_end) if (time_end <= CS%time_max_source) then do j=js,je ; do i=is,ie ! Input an arbitrary energy point source.id_