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 88b062472f..05494fb819 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 @@ -20,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 @@ -37,9 +39,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- @@ -48,15 +50,16 @@ 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 !! for internal tide for testing (BDM) + logical :: debug !< debugging prints end type wave_structure_CS @@ -107,75 +110,82 @@ 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 L2 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]. + 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 [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 - 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] - 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 [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 :: f2 ! squared Coriolis frequency [T-2 ~> s-2] - real :: Kmag2 ! magnitude of horizontal wave number squared - 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) - ! real, dimension(SZK_(G)+1) :: dWdz_profile ! profile of dW/dz - real :: w2avg ! average of squared vertical velocity structure funtion + 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 :: 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, 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 @@ -195,7 +205,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%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) + 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 +283,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 +404,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) = 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) @@ -405,9 +423,19 @@ 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 = 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 +447,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) @@ -438,8 +466,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) @@ -462,12 +495,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 +510,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(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 int_dwdz2 = 0.0 ; int_w2 = 0.0 ; int_N2w2 = 0.0 @@ -489,16 +520,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 +541,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 +553,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 +585,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%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 !> Solves a tri-diagonal system Ax=y using either the standard @@ -640,7 +676,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 @@ -709,6 +745,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 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..79c69be095 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 @@ -114,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)) @@ -122,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_ @@ -140,10 +142,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 +415,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)