diff --git a/src/diagnostics/MOM_wave_speed.F90 b/src/diagnostics/MOM_wave_speed.F90 index c2b671f1c6..59dbfc184e 100644 --- a/src/diagnostics/MOM_wave_speed.F90 +++ b/src/diagnostics/MOM_wave_speed.F90 @@ -7,11 +7,12 @@ module MOM_wave_speed use MOM_error_handler, only : MOM_error, FATAL, WARNING use MOM_file_parser, only : log_version use MOM_grid, only : ocean_grid_type +use MOM_interface_heights, only : thickness_to_dz use MOM_remapping, only : remapping_CS, initialize_remapping, remapping_core_h, interpolate_column use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type -use MOM_EOS, only : calculate_density_derivs +use MOM_EOS, only : calculate_density_derivs, calculate_specific_vol_derivs implicit none ; private @@ -62,7 +63,7 @@ module MOM_wave_speed contains !> Calculates the wave speed of the first baroclinic mode. -subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_N2_column_fraction, & +subroutine wave_speed(h, tv, G, GV, US, cg1, CS, halo_size, use_ebt_mode, mono_N2_column_fraction, & mono_N2_depth, modal_structure) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure @@ -72,8 +73,8 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_ type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables real, dimension(SZI_(G),SZJ_(G)), intent(out) :: cg1 !< First mode internal wave speed [L T-1 ~> m s-1] type(wave_speed_CS), intent(in) :: CS !< Wave speed control struct - logical, optional, intent(in) :: full_halos !< If true, do the calculation - !! over the entire computational domain. + integer, optional, intent(in) :: halo_size !< Width of halo within which to + !! calculate wave speeds logical, optional, intent(in) :: use_ebt_mode !< If true, use the equivalent !! barotropic mode instead of the first baroclinic mode. real, optional, intent(in) :: mono_N2_column_fraction !< The lower fraction @@ -89,27 +90,28 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_ real, dimension(SZK_(GV)+1) :: & dRho_dT, & ! Partial derivative of density with temperature [R C-1 ~> kg m-3 degC-1] dRho_dS, & ! Partial derivative of density with salinity [R S-1 ~> kg m-3 ppt-1] + dSpV_dT, & ! Partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1] + dSpV_dS, & ! Partial derivative of specific volume with salinity [R-1 S-1 ~> m3 kg-1 ppt-1] pres, & ! Interface pressure [R L2 T-2 ~> Pa] T_int, & ! Temperature interpolated to interfaces [C ~> degC] S_int, & ! Salinity interpolated to interfaces [S ~> ppt] - H_top, & ! The distance of each filtered interface from the ocean surface [Z ~> m] - H_bot, & ! The distance of each filtered interface from the bottom [Z ~> m] - gprime ! The reduced gravity across each interface [L2 Z-1 T-2 ~> m s-2]. + H_top, & ! The distance of each filtered interface from the ocean surface [H ~> m or kg m-2] + H_bot, & ! The distance of each filtered interface from the bottom [H ~> m or kg m-2] + gprime ! The reduced gravity across each interface [L2 H-1 T-2 ~> m s-2 or m4 s-1 kg-1]. real, dimension(SZK_(GV)) :: & Igl, Igu ! The inverse of the reduced gravity across an interface times ! the thickness of the layer below (Igl) or above (Igu) it, in [T2 L-2 ~> s2 m-2]. real, dimension(SZK_(GV),SZI_(G)) :: & - Hf, & ! Layer thicknesses after very thin layers are combined [Z ~> m] + Hf, & ! Layer thicknesses after very thin layers are combined [H ~> m or kg m-2] Tf, & ! Layer temperatures after very thin layers are combined [C ~> degC] Sf, & ! Layer salinities after very thin layers are combined [S ~> ppt] Rf ! Layer densities after very thin layers are combined [R ~> kg m-3] real, dimension(SZK_(GV)) :: & - Hc, & ! A column of layer thicknesses after convective instabilities are removed [Z ~> m] + Hc, & ! A column of layer thicknesses after convective instabilities are removed [H ~> m or kg m-2] Tc, & ! A column of layer temperatures after convective instabilities are removed [C ~> degC] Sc, & ! A column of layer salinities after convective instabilities are removed [S ~> ppt] - Rc, & ! A column of layer densities after convective instabilities are removed [R ~> kg m-3] - Hc_H ! Hc(:) rescaled from Z to thickness units [H ~> m or kg m-2] - real :: I_Htot ! The inverse of the total filtered thicknesses [Z ~> m] + Rc ! A column of layer densities after convective instabilities are removed [R ~> kg m-3] + real :: I_Htot ! The inverse of the total filtered thicknesses [H-1 ~> m-1 or m2 kg-1] real :: det, ddet ! Determinant of the eigen system and its derivative with lam. Because the ! units of the eigenvalue change with the number of layers and because of the ! dynamic rescaling that is used to keep det in a numerically representable range, @@ -118,18 +120,21 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_ real :: lam ! The eigenvalue [T2 L-2 ~> s2 m-2] real :: dlam ! The change in estimates of the eigenvalue [T2 L-2 ~> s2 m-2] real :: lam0 ! The first guess of the eigenvalue [T2 L-2 ~> s2 m-2] - real :: Z_to_pres ! A conversion factor from thicknesses to pressure [R L2 T-2 Z-1 ~> Pa m-1] + real :: H_to_pres ! A conversion factor from thicknesses to pressure [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1] real, dimension(SZI_(G)) :: & - htot, hmin, & ! Thicknesses [Z ~> m] - H_here, & ! A thickness [Z ~> m] - HxT_here, & ! A layer integrated temperature [C Z ~> degC m] - HxS_here, & ! A layer integrated salinity [S Z ~> ppt m] - HxR_here ! A layer integrated density [R Z ~> kg m-2] + htot, hmin, & ! Thicknesses [H ~> m or kg m-2] + H_here, & ! A thickness [H ~> m or kg m-2] + HxT_here, & ! A layer integrated temperature [C H ~> degC m or degC kg m-2] + HxS_here, & ! A layer integrated salinity [S H ~> ppt m or ppt kg m-2] + HxR_here ! A layer integrated density [R H ~> kg m-2 or kg2 m-5] real :: speed2_tot ! overestimate of the mode-1 speed squared [L2 T-2 ~> m2 s-2] real :: cg1_min2 ! A floor in the squared first mode speed below which 0 is returned [L2 T-2 ~> m2 s-2] - real :: I_Hnew ! The inverse of a new layer thickness [Z-1 ~> m-1] - real :: drxh_sum ! The sum of density differences across interfaces times thicknesses [R Z ~> kg m-2] - real :: g_Rho0 ! G_Earth/Rho0 [L2 T-2 Z-1 R-1 ~> m4 s-2 kg-1]. + real :: cg1_est ! An initial estimate of the squared first mode speed [L2 T-2 ~> m2 s-2] + real :: I_Hnew ! The inverse of a new layer thickness [H-1 ~> m-1 or m2 kg-1] + real :: drxh_sum ! The sum of density differences across interfaces times thicknesses [R H ~> kg m-2 or kg2 m-5] + real :: dSpVxh_sum ! The sum of specific volume differences across interfaces times + ! thicknesses [R-1 H ~> m4 kg-1 or m], negative for stable stratification. + real :: g_Rho0 ! G_Earth/Rho0 [L2 T-2 H-1 R-1 ~> m4 s-2 kg-1 or m7 s-2 kg-2]. real :: c2_scale ! A scaling factor for wave speeds to help control the growth of the determinant and ! its derivative with lam between rows of the Thomas algorithm solver [L2 s2 T-2 m-2 ~> nondim]. ! The exact value should not matter for the final result if it is an even power of 2. @@ -148,15 +153,16 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_ ! with each iteration. Because of all of the dynamic rescaling of the determinant ! between rows, its units are not easily interpretable, but the ratio of det/ddet ! always has units of [T2 L-2 ~> s2 m-2] - logical :: use_EOS ! If true, density is calculated from T & S using an equation of state. + logical :: use_EOS ! If true, density or specific volume is calculated from T & S using an equation of state. + logical :: nonBous ! If true, do not make the Boussinesq approximation. logical :: better_est ! If true, use an improved estimate of the first mode internal wave speed. logical :: merge ! If true, merge the current layer with the one above. integer :: kc ! The number of layers in the column after merging - integer :: i, j, k, k2, itt, is, ie, js, je, nz - real :: hw ! The mean of the adjacent layer thicknesses [Z ~> m] - real :: sum_hc ! The sum of the layer thicknesses [Z ~> m] - real :: gp ! A limited local copy of gprime [L2 Z-1 T-2 ~> m s-2] - real :: N2min ! A minimum buoyancy frequency, including a slope rescaling factor [L2 Z-2 T-2 ~> s-2] + integer :: i, j, k, k2, itt, is, ie, js, je, nz, halo + real :: hw ! The mean of the adjacent layer thicknesses [H ~> m or kg m-2] + real :: sum_hc ! The sum of the layer thicknesses [H ~> m or kg m-2] + real :: gp ! A limited local copy of gprime [L2 H-1 T-2 ~> m s-2 or m4 s-1 kg-1] + real :: N2min ! A minimum buoyancy frequency, including a slope rescaling factor [L2 H-2 T-2 ~> s-2 or m6 kg-2 s-2] logical :: below_mono_N2_frac ! True if an interface is below the fractional depth where N2 should not increase. logical :: below_mono_N2_depth ! True if an interface is below the absolute depth where N2 should not increase. logical :: l_use_ebt_mode, calc_modal_structure @@ -167,14 +173,15 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_ real :: ms_min, ms_max ! The minimum and maximum mode structure values returned from tdma6 [L2 T-2 ~> m2 s-2] real :: ms_sq ! The sum of the square of the values returned from tdma6 [L4 T-4 ~> m4 s-4] - is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke ; halo = 0 if (.not. CS%initialized) call MOM_error(FATAL, "MOM_wave_speed / wave_speed: "// & "Module must be initialized before it is used.") - if (present(full_halos)) then ; if (full_halos) then - is = G%isd ; ie = G%ied ; js = G%jsd ; je = G%jed - endif ; endif + if (present(halo_size)) then + halo = halo_size + is = G%isc - halo ; ie = G%iec + halo ; js = G%jsc - halo ; je = G%jec + halo + endif l_use_ebt_mode = CS%use_ebt_mode if (present(use_ebt_mode)) l_use_ebt_mode = use_ebt_mode @@ -190,9 +197,10 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_ enddo ; enddo ; enddo endif - g_Rho0 = GV%g_Earth / GV%Rho0 - ! Simplifying the following could change answers at roundoff. - Z_to_pres = GV%Z_to_H * (GV%H_to_RZ * GV%g_Earth) + g_Rho0 = GV%g_Earth*GV%H_to_Z / GV%Rho0 + H_to_pres = GV%H_to_RZ * GV%g_Earth + ! Note that g_Rho0 = H_to_pres / GV%Rho0**2 + nonBous = .not.(GV%Boussinesq .or. GV%semi_Boussinesq) use_EOS = associated(tv%eqn_of_state) better_est = CS%better_cg1_est @@ -216,17 +224,17 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_ c2_scale = US%m_s_to_L_T**2 / 4096.0**2 ! Other powers of 2 give identical results. min_h_frac = tol_Hfrac / real(nz) - !$OMP parallel do default(private) shared(is,ie,js,je,nz,h,G,GV,US,tv,use_EOS, & + !$OMP parallel do default(private) shared(is,ie,js,je,nz,h,G,GV,US,tv,use_EOS,nonBous, & !$OMP CS,min_h_frac,calc_modal_structure,l_use_ebt_mode, & !$OMP modal_structure,l_mono_N2_column_fraction,l_mono_N2_depth, & - !$OMP Z_to_pres,cg1,g_Rho0,rescale,I_rescale,cg1_min2, & + !$OMP H_to_pres,cg1,g_Rho0,rescale,I_rescale,cg1_min2, & !$OMP better_est,tol_solve,tol_merge,c2_scale) do j=js,je ! First merge very thin layers with the one above (or below if they are ! at the top). This also transposes the row order so that columns can ! be worked upon one at a time. do i=is,ie ; htot(i) = 0.0 ; enddo - do k=1,nz ; do i=is,ie ; htot(i) = htot(i) + h(i,j,k)*GV%H_to_Z ; enddo ; enddo + do k=1,nz ; do i=is,ie ; htot(i) = htot(i) + h(i,j,k) ; enddo ; enddo do i=is,ie hmin(i) = htot(i)*min_h_frac ; kf(i) = 1 ; H_here(i) = 0.0 @@ -234,20 +242,20 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_ enddo if (use_EOS) then do k=1,nz ; do i=is,ie - if ((H_here(i) > hmin(i)) .and. (h(i,j,k)*GV%H_to_Z > hmin(i))) then + if ((H_here(i) > hmin(i)) .and. (h(i,j,k) > hmin(i))) then Hf(kf(i),i) = H_here(i) Tf(kf(i),i) = HxT_here(i) / H_here(i) Sf(kf(i),i) = HxS_here(i) / H_here(i) kf(i) = kf(i) + 1 ! Start a new layer - H_here(i) = h(i,j,k)*GV%H_to_Z - HxT_here(i) = (h(i,j,k) * GV%H_to_Z) * tv%T(i,j,k) - HxS_here(i) = (h(i,j,k) * GV%H_to_Z) * tv%S(i,j,k) + H_here(i) = h(i,j,k) + HxT_here(i) = h(i,j,k) * tv%T(i,j,k) + HxS_here(i) = h(i,j,k) * tv%S(i,j,k) else - H_here(i) = H_here(i) + h(i,j,k)*GV%H_to_Z - HxT_here(i) = HxT_here(i) + (h(i,j,k) * GV%H_to_Z) * tv%T(i,j,k) - HxS_here(i) = HxS_here(i) + (h(i,j,k) * GV%H_to_Z) * tv%S(i,j,k) + H_here(i) = H_here(i) + h(i,j,k) + HxT_here(i) = HxT_here(i) + h(i,j,k) * tv%T(i,j,k) + HxS_here(i) = HxS_here(i) + h(i,j,k) * tv%S(i,j,k) endif enddo ; enddo do i=is,ie ; if (H_here(i) > 0.0) then @@ -255,18 +263,18 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_ Tf(kf(i),i) = HxT_here(i) / H_here(i) Sf(kf(i),i) = HxS_here(i) / H_here(i) endif ; enddo - else + else ! .not. (use_EOS) do k=1,nz ; do i=is,ie - if ((H_here(i) > hmin(i)) .and. (h(i,j,k)*GV%H_to_Z > hmin(i))) then + if ((H_here(i) > hmin(i)) .and. (h(i,j,k) > hmin(i))) then Hf(kf(i),i) = H_here(i) ; Rf(kf(i),i) = HxR_here(i) / H_here(i) kf(i) = kf(i) + 1 ! Start a new layer - H_here(i) = h(i,j,k)*GV%H_to_Z - HxR_here(i) = (h(i,j,k)*GV%H_to_Z)*GV%Rlay(k) + H_here(i) = h(i,j,k) + HxR_here(i) = h(i,j,k)*GV%Rlay(k) else - H_here(i) = H_here(i) + h(i,j,k)*GV%H_to_Z - HxR_here(i) = HxR_here(i) + (h(i,j,k)*GV%H_to_Z)*GV%Rlay(k) + H_here(i) = H_here(i) + h(i,j,k) + HxR_here(i) = HxR_here(i) + h(i,j,k)*GV%Rlay(k) endif enddo ; enddo do i=is,ie ; if (H_here(i) > 0.0) then @@ -279,16 +287,21 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_ if (use_EOS) then pres(1) = 0.0 ; H_top(1) = 0.0 do K=2,kf(i) - pres(K) = pres(K-1) + Z_to_pres*Hf(k-1,i) + pres(K) = pres(K-1) + H_to_pres*Hf(k-1,i) T_int(K) = 0.5*(Tf(k,i)+Tf(k-1,i)) S_int(K) = 0.5*(Sf(k,i)+Sf(k-1,i)) H_top(K) = H_top(K-1) + Hf(k-1,i) enddo - call calculate_density_derivs(T_int, S_int, pres, drho_dT, drho_dS, & - tv%eqn_of_state, (/2,kf(i)/) ) + if (nonBous) then + call calculate_specific_vol_derivs(T_int, S_int, pres, dSpV_dT, dSpV_dS, & + tv%eqn_of_state, (/2,kf(i)/) ) + else + call calculate_density_derivs(T_int, S_int, pres, drho_dT, drho_dS, & + tv%eqn_of_state, (/2,kf(i)/) ) + endif ! Sum the reduced gravities to find out how small a density difference is negligibly small. - drxh_sum = 0.0 + drxh_sum = 0.0 ; dSpVxh_sum = 0.0 if (better_est) then ! This is an estimate that is correct for the non-EBT mode for 2 or 3 layers, or for ! clusters of massless layers at interfaces that can be grouped into 2 or 3 layers. @@ -297,44 +310,81 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_ if (H_top(kf(i)) > 0.0) then I_Htot = 1.0 / (H_top(kf(i)) + Hf(kf(i),i)) ! = 1.0 / (H_top(K) + H_bot(K)) for all K. H_bot(kf(i)+1) = 0.0 - do K=kf(i),2,-1 - H_bot(K) = H_bot(K+1) + Hf(k,i) - drxh_sum = drxh_sum + ((H_top(K) * H_bot(K)) * I_Htot) * & - max(0.0, drho_dT(K)*(Tf(k,i)-Tf(k-1,i)) + drho_dS(K)*(Sf(k,i)-Sf(k-1,i))) - enddo + if (nonBous) then + do K=kf(i),2,-1 + H_bot(K) = H_bot(K+1) + Hf(k,i) + dSpVxh_sum = dSpVxh_sum + ((H_top(K) * H_bot(K)) * I_Htot) * & + min(0.0, dSpV_dT(K)*(Tf(k,i)-Tf(k-1,i)) + dSpV_dS(K)*(Sf(k,i)-Sf(k-1,i))) + enddo + else + do K=kf(i),2,-1 + H_bot(K) = H_bot(K+1) + Hf(k,i) + drxh_sum = drxh_sum + ((H_top(K) * H_bot(K)) * I_Htot) * & + max(0.0, drho_dT(K)*(Tf(k,i)-Tf(k-1,i)) + drho_dS(K)*(Sf(k,i)-Sf(k-1,i))) + enddo + endif endif else ! This estimate is problematic in that it goes like 1/nz for a large number of layers, ! but it is an overestimate (as desired) for a small number of layers, by at a factor ! of (H1+H2)**2/(H1*H2) >= 4 for two thick layers. - do K=2,kf(i) - drxh_sum = drxh_sum + 0.5*(Hf(k-1,i)+Hf(k,i)) * & - max(0.0, drho_dT(K)*(Tf(k,i)-Tf(k-1,i)) + drho_dS(K)*(Sf(k,i)-Sf(k-1,i))) - enddo + if (nonBous) then + do K=2,kf(i) + dSpVxh_sum = dSpVxh_sum + 0.5*(Hf(k-1,i)+Hf(k,i)) * & + min(0.0, dSpV_dT(K)*(Tf(k,i)-Tf(k-1,i)) + dSpV_dS(K)*(Sf(k,i)-Sf(k-1,i))) + enddo + else + do K=2,kf(i) + drxh_sum = drxh_sum + 0.5*(Hf(k-1,i)+Hf(k,i)) * & + max(0.0, drho_dT(K)*(Tf(k,i)-Tf(k-1,i)) + drho_dS(K)*(Sf(k,i)-Sf(k-1,i))) + enddo + endif endif - else - drxh_sum = 0.0 + else ! .not. (use_EOS) + drxh_sum = 0.0 ; dSpVxh_sum = 0.0 if (better_est) then H_top(1) = 0.0 do K=2,kf(i) ; H_top(K) = H_top(K-1) + Hf(k-1,i) ; enddo if (H_top(kf(i)) > 0.0) then I_Htot = 1.0 / (H_top(kf(i)) + Hf(kf(i),i)) ! = 1.0 / (H_top(K) + H_bot(K)) for all K. H_bot(kf(i)+1) = 0.0 - do K=kf(i),2,-1 - H_bot(K) = H_bot(K+1) + Hf(k,i) - drxh_sum = drxh_sum + ((H_top(K) * H_bot(K)) * I_Htot) * max(0.0,Rf(k,i)-Rf(k-1,i)) - enddo + if (nonBous) then + do K=kf(i),2,-1 + H_bot(K) = H_bot(K+1) + Hf(k,i) + dSpVxh_sum = dSpVxh_sum + ((H_top(K) * H_bot(K)) * I_Htot) * & + min(0.0, (Rf(k-1,i)-Rf(k,i)) / (Rf(k,i)*Rf(k-1,i))) + enddo + else + do K=kf(i),2,-1 + H_bot(K) = H_bot(K+1) + Hf(k,i) + drxh_sum = drxh_sum + ((H_top(K) * H_bot(K)) * I_Htot) * max(0.0,Rf(k,i)-Rf(k-1,i)) + enddo + endif endif else - do K=2,kf(i) - drxh_sum = drxh_sum + 0.5*(Hf(k-1,i)+Hf(k,i)) * max(0.0,Rf(k,i)-Rf(k-1,i)) - enddo + if (nonBous) then + do K=2,kf(i) + dSpVxh_sum = dSpVxh_sum + 0.5*(Hf(k-1,i)+Hf(k,i)) * & + min(0.0, (Rf(k-1,i)-Rf(k,i)) / (Rf(k,i)*Rf(k-1,i))) + enddo + else + do K=2,kf(i) + drxh_sum = drxh_sum + 0.5*(Hf(k-1,i)+Hf(k,i)) * max(0.0,Rf(k,i)-Rf(k-1,i)) + enddo + endif endif endif ! use_EOS + if (nonBous) then + ! Note that dSpVxh_sum is negative for stable stratification. + cg1_est = H_to_pres * abs(dSpVxh_sum) + else + cg1_est = g_Rho0 * drxh_sum + endif + ! Find gprime across each internal interface, taking care of convective instabilities by ! merging layers. If the estimated wave speed is too small, simply return zero. - if (g_Rho0 * drxh_sum <= cg1_min2) then + if (cg1_est <= cg1_min2) then cg1(i,j) = 0.0 if (present(modal_structure)) modal_structure(i,j,:) = 0. else @@ -345,9 +395,15 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_ kc = 1 Hc(1) = Hf(1,i) ; Tc(1) = Tf(1,i) ; Sc(1) = Sf(1,i) do k=2,kf(i) - if (better_est) then + if (better_est .and. nonBous) then + merge = ((dSpV_dT(K)*(Tc(kc)-Tf(k,i)) + dSpV_dS(K)*(Sc(kc)-Sf(k,i))) * & + ((Hc(kc) * Hf(k,i))*I_Htot) < abs(2.0 * tol_merge * dSpVxh_sum)) + elseif (better_est) then merge = ((drho_dT(K)*(Tf(k,i)-Tc(kc)) + drho_dS(K)*(Sf(k,i)-Sc(kc))) * & ((Hc(kc) * Hf(k,i))*I_Htot) < 2.0 * tol_merge*drxh_sum) + elseif (nonBous) then + merge = ((dSpV_dT(K)*(Tc(kc)-Tf(k,i)) + dSpV_dS(K)*(Sc(kc)-Sf(k,i))) * & + (Hc(kc) + Hf(k,i)) < abs(2.0 * tol_merge * dSpVxh_sum)) else merge = ((drho_dT(K)*(Tf(k,i)-Tc(kc)) + drho_dS(K)*(Sf(k,i)-Sc(kc))) * & (Hc(kc) + Hf(k,i)) < 2.0 * tol_merge*drxh_sum) @@ -362,9 +418,15 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_ ! that the tolerance is a factor of two larger, to avoid limit how ! far back we go. do K2=kc,2,-1 - if (better_est) then + if (better_est .and. nonBous) then + merge = ( (dSpV_dT(K2)*(Tc(k2-1)-Tc(k2)) + dSpV_dS(K2)*(Sc(k2-1)-Sc(k2))) * & + ((Hc(k2) * Hc(k2-1))*I_Htot) < abs(tol_merge * dSpVxh_sum) ) + elseif (better_est) then merge = ((drho_dT(K2)*(Tc(k2)-Tc(k2-1)) + drho_dS(K2)*(Sc(k2)-Sc(k2-1))) * & ((Hc(k2) * Hc(k2-1))*I_Htot) < tol_merge*drxh_sum) + elseif (nonBous) then + merge = ( (dSpV_dT(K2)*(Tc(k2-1)-Tc(k2)) + dSpV_dS(K2)*(Sc(k2-1)-Sc(k2))) * & + (Hc(k2) + Hc(k2-1)) < abs(tol_merge * dSpVxh_sum) ) else merge = ((drho_dT(K2)*(Tc(k2)-Tc(k2-1)) + drho_dS(K2)*(Sc(k2)-Sc(k2-1))) * & (Hc(k2) + Hc(k2-1)) < tol_merge*drxh_sum) @@ -381,20 +443,36 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_ else ! Add a new layer to the column. kc = kc + 1 - drho_dS(Kc) = drho_dS(K) ; drho_dT(Kc) = drho_dT(K) + if (nonBous) then + dSpV_dS(Kc) = dSpV_dS(K) ; dSpV_dT(Kc) = dSpV_dT(K) + else + drho_dS(Kc) = drho_dS(K) ; drho_dT(Kc) = drho_dT(K) + endif Tc(kc) = Tf(k,i) ; Sc(kc) = Sf(k,i) ; Hc(kc) = Hf(k,i) endif enddo ! At this point there are kc layers and the gprimes should be positive. - do K=2,kc ! Revisit this if non-Boussinesq. - gprime(K) = g_Rho0 * (drho_dT(K)*(Tc(k)-Tc(k-1)) + drho_dS(K)*(Sc(k)-Sc(k-1))) - enddo - else ! .not.use_EOS + if (nonBous) then + do K=2,kc + gprime(K) = H_to_pres * (dSpV_dT(K)*(Tc(k-1)-Tc(k)) + dSpV_dS(K)*(Sc(k-1)-Sc(k))) + enddo + else + do K=2,kc + gprime(K) = g_Rho0 * (drho_dT(K)*(Tc(k)-Tc(k-1)) + drho_dS(K)*(Sc(k)-Sc(k-1))) + enddo + endif + else ! .not. (use_EOS) ! Do the same with density directly... kc = 1 Hc(1) = Hf(1,i) ; Rc(1) = Rf(1,i) do k=2,kf(i) - if (better_est) then + if (nonBous .and. better_est) then + merge = ((Rf(k,i) - Rc(kc)) * ((Hc(kc) * Hf(k,i))*I_Htot) < & + (Rc(kc)*Rf(k,i)) * abs(2.0 * tol_merge * dSpVxh_sum)) + elseif (nonBous) then + merge = ((Rf(k,i) - Rc(kc)) * (Hc(kc) + Hf(k,i)) < & + (Rc(kc)*Rf(k,i)) * abs(2.0 * tol_merge * dSpVxh_sum)) + elseif (better_est) then merge = ((Rf(k,i) - Rc(kc)) * ((Hc(kc) * Hf(k,i))*I_Htot) < 2.0*tol_merge*drxh_sum) else merge = ((Rf(k,i) - Rc(kc)) * (Hc(kc) + Hf(k,i)) < 2.0*tol_merge*drxh_sum) @@ -407,7 +485,13 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_ ! that the tolerance is a factor of two larger, to avoid limit how ! far back we go. do k2=kc,2,-1 - if (better_est) then + if (nonBous .and. better_est) then + merge = ((Rc(k2) - Rc(k2-1)) * ((Hc(kc) * Hf(k,i))*I_Htot) < & + (Rc(k2-1)*Rc(k2)) * abs(2.0 * tol_merge * dSpVxh_sum)) + elseif (nonBous) then + merge = ((Rc(k2) - Rc(k2-1)) * (Hc(kc) + Hf(k,i)) < & + (Rc(k2-1)*Rc(k2)) * abs(2.0 * tol_merge * dSpVxh_sum)) + elseif (better_est) then merge = ((Rc(k2)-Rc(k2-1)) * ((Hc(k2) * Hc(k2-1))*I_Htot) < tol_merge*drxh_sum) else merge = ((Rc(k2)-Rc(k2-1)) * (Hc(k2)+Hc(k2-1)) < tol_merge*drxh_sum) @@ -426,9 +510,15 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_ endif enddo ! At this point there are kc layers and the gprimes should be positive. - do K=2,kc ! Revisit this if non-Boussinesq. - gprime(K) = g_Rho0 * (Rc(k)-Rc(k-1)) - enddo + if (nonBous) then + do K=2,kc + gprime(K) = H_to_pres * (Rc(k) - Rc(k-1)) / (Rc(k) * Rc(k-1)) + enddo + else + do K=2,kc + gprime(K) = g_Rho0 * (Rc(k)-Rc(k-1)) + enddo + endif endif ! use_EOS ! Sum the contributions from all of the interfaces to give an over-estimate @@ -453,14 +543,18 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_ do k=2,kc hw = 0.5*(Hc(k-1)+Hc(k)) gp = gprime(K) + if (l_mono_N2_column_fraction>0. .or. l_mono_N2_depth>=0.) then - ! Determine whether N2 estimates should not be allowed to increase with depth. + ! Determine whether N2 estimates should not be allowed to increase with depth. if (l_mono_N2_column_fraction>0.) then - !### Change to: (htot(i) - sum_hc < l_mono_N2_column_fraction*htot(i)) - below_mono_N2_frac = ((G%bathyT(i,j)+G%Z_ref) - GV%H_to_Z*sum_hc < & - l_mono_N2_column_fraction*(G%bathyT(i,j)+G%Z_ref)) + if (GV%Boussinesq .or. GV%semi_Boussinesq) then + below_mono_N2_frac = ((G%bathyT(i,j)+G%Z_ref) - GV%H_to_Z*sum_hc < & + l_mono_N2_column_fraction*(G%bathyT(i,j)+G%Z_ref)) + else + below_mono_N2_frac = (htot(i) - sum_hc < l_mono_N2_column_fraction*htot(i)) + endif endif - if (l_mono_N2_depth >= 0.) below_mono_N2_depth = (sum_hc > GV%H_to_Z*l_mono_N2_depth) + if (l_mono_N2_depth >= 0.) below_mono_N2_depth = (sum_hc > l_mono_N2_depth) if ( (gp > N2min*hw) .and. (below_mono_N2_frac .or. below_mono_N2_depth) ) then ! Filters out regions where N2 increases with depth, but only in a lower fraction @@ -578,17 +672,13 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_ else mode_struct(1:kc)=0. endif - ! Note that remapping_core_h requires that the same units be used - ! for both the source and target grid thicknesses, here [H ~> m or kg m-2]. - do k = 1,kc - Hc_H(k) = GV%Z_to_H * Hc(k) - enddo + if (CS%remap_answer_date < 20190101) then - call remapping_core_h(CS%remapping_CS, kc, Hc_H(:), mode_struct, & + call remapping_core_h(CS%remapping_CS, kc, Hc(:), mode_struct, & nz, h(i,j,:), modal_structure(i,j,:), & 1.0e-30*GV%m_to_H, 1.0e-10*GV%m_to_H) else - call remapping_core_h(CS%remapping_CS, kc, Hc_H(:), mode_struct, & + call remapping_core_h(CS%remapping_CS, kc, Hc(:), mode_struct, & nz, h(i,j,:), modal_structure(i,j,:), & GV%H_subroundoff, GV%H_subroundoff) endif @@ -658,7 +748,7 @@ end subroutine tdma6 !> Calculates the wave speeds for the first few barolinic modes. subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_struct_max, u_struct_bot, Nb, int_w2, & - int_U2, int_N2w2, full_halos) + int_U2, int_N2w2, halo_size) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -666,50 +756,56 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables integer, intent(in) :: nmodes !< Number of modes type(wave_speed_CS), intent(in) :: CS !< Wave speed control struct - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1,nmodes),intent(out) :: w_struct !< Wave Vertical profile [nondim] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV),nmodes),intent(out) :: u_struct !< Wave Horizontal profile [Z-1 ~> m-1] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1,nmodes),intent(out) :: w_struct !< Wave vertical velocity profile [nondim] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV),nmodes),intent(out) :: u_struct !< Wave horizontal velocity profile + !! [Z-1 ~> m-1] real, dimension(SZI_(G),SZJ_(G),nmodes), intent(out) :: cn !< Waves speeds [L T-1 ~> m s-1] - real, dimension(SZI_(G),SZJ_(G),nmodes), intent(out) :: u_struct_max !< Maximum of wave horizontal profile - !! [Z-1 ~> m-1] + real, dimension(SZI_(G),SZJ_(G),nmodes), intent(out) :: u_struct_max !< Maximum of wave horizontal velocity + !! profile [Z-1 ~> m-1] real, dimension(SZI_(G),SZJ_(G),nmodes), intent(out) :: u_struct_bot !< Bottom value of wave horizontal - !! profile [Z-1 ~> m-1] - real, dimension(SZI_(G),SZJ_(G)), intent(out) :: Nb !< Bottom value of Brunt Vaissalla freqency + !! velocity profile [Z-1 ~> m-1] + real, dimension(SZI_(G),SZJ_(G)), intent(out) :: Nb !< Bottom value of buoyancy freqency !! [T-1 ~> s-1] - real, dimension(SZI_(G),SZJ_(G),nmodes), intent(out) :: int_w2 !< depth-integrated - !! vertical profile squared [Z ~> m] - real, dimension(SZI_(G),SZJ_(G),nmodes), intent(out) :: int_U2 !< depth-integrated - !! horizontal profile squared [Z-1 ~> m-1] - real, dimension(SZI_(G),SZJ_(G),nmodes), intent(out) :: int_N2w2 !< depth-integrated Brunt Vaissalla - !! frequency times vertical - !! profile squared [Z T-2 ~> m s-2] - logical, optional, intent(in) :: full_halos !< If true, do the calculation - !! over the entire data domain. + real, dimension(SZI_(G),SZJ_(G),nmodes), intent(out) :: int_w2 !< depth-integrated vertical velocity + !! profile squared [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),nmodes), intent(out) :: int_U2 !< depth-integrated horizontal velocity + !! profile squared [H Z-2 ~> m-1 or kg m-4] + real, dimension(SZI_(G),SZJ_(G),nmodes), intent(out) :: int_N2w2 !< depth-integrated buoyancy frequency + !! times vertical velocity profile + !! squared [H T-2 ~> m s-2 or kg m-2 s-2] + integer, optional, intent(in) :: halo_size !< Width of halo within which to + !! calculate wave speeds ! Local variables real, dimension(SZK_(GV)+1) :: & dRho_dT, & ! Partial derivative of density with temperature [R C-1 ~> kg m-3 degC-1] dRho_dS, & ! Partial derivative of density with salinity [R S-1 ~> kg m-3 ppt-1] + dSpV_dT, & ! Partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1] + dSpV_dS, & ! Partial derivative of specific volume with salinity [R-1 S-1 ~> m3 kg-1 ppt-1] pres, & ! Interface pressure [R L2 T-2 ~> Pa] T_int, & ! Temperature interpolated to interfaces [C ~> degC] S_int, & ! Salinity interpolated to interfaces [S ~> ppt] - H_top, & ! The distance of each filtered interface from the ocean surface [Z ~> m] - H_bot, & ! The distance of each filtered interface from the bottom [Z ~> m] - gprime, & ! The reduced gravity across each interface [L2 Z-1 T-2 ~> m s-2]. + H_top, & ! The distance of each filtered interface from the ocean surface [H ~> m or kg m-2] + H_bot, & ! The distance of each filtered interface from the bottom [H ~> m or kg m-2] + gprime, & ! The reduced gravity across each interface [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2]. N2 ! The buoyancy freqency squared [T-2 ~> s-2] real, dimension(SZK_(GV),SZI_(G)) :: & - Hf, & ! Layer thicknesses after very thin layers are combined [Z ~> m] + Hf, & ! Layer thicknesses after very thin layers are combined [H ~> m or kg m-2] + dzf, & ! Layer vertical extents after very thin layers are combined [Z ~> m] Tf, & ! Layer temperatures after very thin layers are combined [C ~> degC] Sf, & ! Layer salinities after very thin layers are combined [S ~> ppt] Rf ! Layer densities after very thin layers are combined [R ~> kg m-3] + real, dimension(SZI_(G),SZK_(GV)) :: & + dz_2d ! Height change across layers [Z ~> m] real, dimension(SZK_(GV)) :: & Igl, Igu, & ! The inverse of the reduced gravity across an interface times ! the thickness of the layer below (Igl) or above (Igu) it, in [T2 L-2 ~> s2 m-2]. - Hc, & ! A column of layer thicknesses after convective instabilities are removed [Z ~> m] + Hc, & ! A column of layer thicknesses after convective instabilities are removed [H ~> m or kg m-2] + dzc, & ! A column of layer vertical extents after convective instabilities are removed [Z ~> m] Tc, & ! A column of layer temperatures after convective instabilities are removed [C ~> degC] Sc, & ! A column of layer salinities after convective instabilities are removed [S ~> ppt] - Rc, & ! A column of layer densities after convective instabilities are removed [R ~> kg m-3] - Hc_H ! Hc(:) rescaled from Z to thickness units [H ~> m or kg m-2] - real :: I_Htot ! The inverse of the total filtered thicknesses [Z-1 ~> m-1] + Rc ! A column of layer densities after convective instabilities are removed [R ~> kg m-3] + real :: I_Htot ! The inverse of the total filtered thicknesses [H-1 ~> m-1 or m2 kg-1] real :: c2_scale ! A scaling factor for wave speeds to help control the growth of the determinant and its ! derivative with lam between rows of the Thomas algorithm solver [L2 s2 T-2 m-2 ~> nondim]. ! The exact value should not matter for the final result if it is an even power of 2. @@ -740,20 +836,24 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s xbl, xbr ! lam guesses bracketing a zero-crossing (root) [T2 L-2 ~> s2 m-2] integer :: numint ! number of widows (intervals) in root searching range integer :: nrootsfound ! number of extra roots found (not including 1st root) - real :: Z_to_pres ! A conversion factor from thicknesses to pressure [R L2 T-2 Z-1 ~> Pa m-1] + real :: H_to_pres ! A conversion factor from thicknesses to pressure [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1] real, dimension(SZI_(G)) :: & - htot, hmin, & ! Thicknesses [Z ~> m] - H_here, & ! A thickness [Z ~> m] - HxT_here, & ! A layer integrated temperature [C Z ~> degC m] - HxS_here, & ! A layer integrated salinity [S Z ~> ppt m] - HxR_here ! A layer integrated density [R Z ~> kg m-2] + htot, hmin, & ! Thicknesses [H ~> m or kg m-2] + H_here, & ! A layer thickness [H ~> m or kg m-2] + dz_here, & ! A layer vertical extent [Z ~> m] + HxT_here, & ! A layer integrated temperature [C H ~> degC m or degC kg m-2] + HxS_here, & ! A layer integrated salinity [S H ~> ppt m or ppt kg m-2] + HxR_here ! A layer integrated density [R H ~> kg m-2 or kg2 m-5] real :: speed2_tot ! overestimate of the mode-1 speed squared [L2 T-2 ~> m2 s-2] real :: speed2_min ! minimum mode speed (squared) to consider in root searching [L2 T-2 ~> m2 s-2] real :: cg1_min2 ! A floor in the squared first mode speed below which 0 is returned [L2 T-2 ~> m2 s-2] + real :: cg1_est ! An initial estimate of the squared first mode speed [L2 T-2 ~> m2 s-2] real, parameter :: reduct_factor = 0.5 ! A factor used in setting speed2_min [nondim] - real :: I_Hnew ! The inverse of a new layer thickness [Z-1 ~> m-1] - real :: drxh_sum ! The sum of density differences across interfaces times thicknesses [R Z ~> kg m-2] - real :: g_Rho0 ! G_Earth/Rho0 [L2 T-2 Z-1 R-1 ~> m4 s-2 kg-1]. + real :: I_Hnew ! The inverse of a new layer thickness [H-1 ~> m-1 or m2 kg-1] + real :: drxh_sum ! The sum of density differences across interfaces times thicknesses [R H ~> kg m-2 or kg2 m-5] + real :: dSpVxh_sum ! The sum of specific volume differences across interfaces times + ! thicknesses [R-1 H ~> m4 kg-1 or m], negative for stable stratification. + real :: g_Rho0 ! G_Earth/Rho0 [L2 T-2 H-1 R-1 ~> m4 s-2 kg-1 pr m7 s-2 kg-1]. real :: tol_Hfrac ! Layers that together are smaller than this fraction of ! the total water column can be merged for efficiency [nondim]. real :: min_h_frac ! tol_Hfrac divided by the total number of layers [nondim]. @@ -762,7 +862,8 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s ! when deciding to merge layers in the calculation [nondim] integer :: kf(SZI_(G)) ! The number of active layers after filtering. integer, parameter :: max_itt = 30 - logical :: use_EOS ! If true, density is calculated from T & S using the equation of state. + logical :: use_EOS ! If true, density or specific volume is calculated from T & S using the equation of state. + logical :: nonBous ! If true, do not make the Boussinesq approximation. logical :: better_est ! If true, use an improved estimate of the first mode internal wave speed. logical :: merge ! If true, merge the current layer with the one above. integer :: nsub ! number of subintervals used for root finding @@ -772,35 +873,36 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s logical :: sub_rootfound ! if true, subdivision has located root integer :: kc ! The number of layers in the column after merging integer :: sub, sub_it - integer :: i, j, k, k2, itt, is, ie, js, je, nz, iint, m + integer :: i, j, k, k2, itt, is, ie, js, je, nz, iint, m, halo real, dimension(SZK_(GV)+1) :: modal_structure !< Normalized model structure [nondim] real, dimension(SZK_(GV)) :: modal_structure_fder !< Normalized model structure [Z-1 ~> m-1] real :: mode_struct(SZK_(GV)+1) ! The mode structure [nondim], but it is also temporarily ! in units of [L2 T-2 ~> m2 s-2] after it is modified inside of tdma6. - real :: mode_struct_fder(SZK_(GV)) ! The mode structure 1st derivative [nondim], but it is also temporarily - ! in units of [L2 T-2 ~> m2 s-2] after it is modified inside of tdma6. + real :: mode_struct_fder(SZK_(GV)) ! The mode structure 1st derivative [Z-1 ~> m-1], but it is also temporarily + ! in units of [Z-1 L2 T-2 ~> m s-2] after it is modified inside of tdma6. real :: mode_struct_sq(SZK_(GV)+1) ! The square of mode structure [nondim] real :: mode_struct_fder_sq(SZK_(GV)) ! The square of mode structure 1st derivative [Z-2 ~> m-2] real :: ms_min, ms_max ! The minimum and maximum mode structure values returned from tdma6 [L2 T-2 ~> m2 s-2] real :: ms_sq ! The sum of the square of the values returned from tdma6 [L4 T-4 ~> m4 s-4] - real :: w2avg ! A total for renormalization - real, parameter :: a_int = 0.5 ! Integral total for normalization - real :: renorm ! Normalization factor + real :: w2avg ! A total for renormalization [H L4 T-4 ~> m5 s-4 or kg m2 s-4] + real, parameter :: a_int = 0.5 ! Integral total for normalization [nondim] + real :: renorm ! Normalization factor [T2 L-2 ~> s2 m-2] - is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke ; halo = 0 if (.not. CS%initialized) call MOM_error(FATAL, "MOM_wave_speed / wave_speeds: "// & "Module must be initialized before it is used.") - if (present(full_halos)) then ; if (full_halos) then - is = G%isd ; ie = G%ied ; js = G%jsd ; je = G%jed - endif ; endif + if (present(halo_size)) then + halo = halo_size + is = G%isc - halo ; ie = G%iec + halo ; js = G%jsc - halo ; je = G%jec + halo + endif - g_Rho0 = GV%g_Earth / GV%Rho0 - ! Simplifying the following could change answers at roundoff. - Z_to_pres = GV%Z_to_H * (GV%H_to_RZ * GV%g_Earth) + g_Rho0 = GV%g_Earth * GV%H_to_Z / GV%Rho0 + H_to_pres = GV%H_to_RZ * GV%g_Earth + nonBous = .not.(GV%Boussinesq .or. GV%semi_Boussinesq) use_EOS = associated(tv%eqn_of_state) if (CS%c1_thresh < 0.0) & @@ -830,59 +932,69 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s w_struct(:,:,:,:) = 0.0 min_h_frac = tol_Hfrac / real(nz) - !$OMP parallel do default(private) shared(is,ie,js,je,nz,h,G,GV,US,CS,min_h_frac,use_EOS, & - !$OMP Z_to_pres,tv,cn,g_Rho0,nmodes,cg1_min2,better_est, & - !$OMP tol_solve,tol_merge,c2_scale) + !$OMP parallel do default(private) shared(is,ie,js,je,nz,h,G,GV,US,CS,use_EOS,nonBous, & + !$OMP min_h_frac,H_to_pres,tv,cn,g_Rho0,nmodes,cg1_min2, & + !$OMP better_est,tol_solve,tol_merge,c2_scale) do j=js,je ! First merge very thin layers with the one above (or below if they are ! at the top). This also transposes the row order so that columns can ! be worked upon one at a time. do i=is,ie ; htot(i) = 0.0 ; enddo - do k=1,nz ; do i=is,ie ; htot(i) = htot(i) + h(i,j,k)*GV%H_to_Z ; enddo ; enddo + do k=1,nz ; do i=is,ie ; htot(i) = htot(i) + h(i,j,k) ; enddo ; enddo + + call thickness_to_dz(h, tv, dz_2d, j, G, GV, halo_size=halo) do i=is,ie - hmin(i) = htot(i)*min_h_frac ; kf(i) = 1 ; H_here(i) = 0.0 + hmin(i) = htot(i)*min_h_frac ; kf(i) = 1 ; H_here(i) = 0.0 ; dz_here(i) = 0.0 HxT_here(i) = 0.0 ; HxS_here(i) = 0.0 ; HxR_here(i) = 0.0 enddo if (use_EOS) then do k=1,nz ; do i=is,ie - if ((H_here(i) > hmin(i)) .and. (h(i,j,k)*GV%H_to_Z > hmin(i))) then + if ((H_here(i) > hmin(i)) .and. (h(i,j,k) > hmin(i))) then Hf(kf(i),i) = H_here(i) + dzf(kf(i),i) = dz_here(i) Tf(kf(i),i) = HxT_here(i) / H_here(i) Sf(kf(i),i) = HxS_here(i) / H_here(i) kf(i) = kf(i) + 1 ! Start a new layer - H_here(i) = h(i,j,k)*GV%H_to_Z - HxT_here(i) = (h(i,j,k)*GV%H_to_Z)*tv%T(i,j,k) - HxS_here(i) = (h(i,j,k)*GV%H_to_Z)*tv%S(i,j,k) + H_here(i) = h(i,j,k) + dz_here(i) = dz_2d(i,k) + HxT_here(i) = h(i,j,k)*tv%T(i,j,k) + HxS_here(i) = h(i,j,k)*tv%S(i,j,k) else - H_here(i) = H_here(i) + h(i,j,k)*GV%H_to_Z - HxT_here(i) = HxT_here(i) + (h(i,j,k)*GV%H_to_Z)*tv%T(i,j,k) - HxS_here(i) = HxS_here(i) + (h(i,j,k)*GV%H_to_Z)*tv%S(i,j,k) + H_here(i) = H_here(i) + h(i,j,k) + dz_here(i) = dz_here(i) + dz_2d(i,k) + HxT_here(i) = HxT_here(i) + h(i,j,k)*tv%T(i,j,k) + HxS_here(i) = HxS_here(i) + h(i,j,k)*tv%S(i,j,k) endif enddo ; enddo do i=is,ie ; if (H_here(i) > 0.0) then Hf(kf(i),i) = H_here(i) + dzf(kf(i),i) = dz_here(i) Tf(kf(i),i) = HxT_here(i) / H_here(i) Sf(kf(i),i) = HxS_here(i) / H_here(i) endif ; enddo - else + else ! .not. (use_EOS) do k=1,nz ; do i=is,ie - if ((H_here(i) > hmin(i)) .and. (h(i,j,k)*GV%H_to_Z > hmin(i))) then + if ((H_here(i) > hmin(i)) .and. (h(i,j,k) > hmin(i))) then Hf(kf(i),i) = H_here(i) ; Rf(kf(i),i) = HxR_here(i) / H_here(i) + dzf(kf(i),i) = dz_here(i) kf(i) = kf(i) + 1 ! Start a new layer - H_here(i) = h(i,j,k)*GV%H_to_Z - HxR_here(i) = (h(i,j,k)*GV%H_to_Z)*GV%Rlay(k) + H_here(i) = h(i,j,k) + dz_here(i) = dz_2d(i,k) + HxR_here(i) = h(i,j,k)*GV%Rlay(k) else - H_here(i) = H_here(i) + h(i,j,k)*GV%H_to_Z - HxR_here(i) = HxR_here(i) + (h(i,j,k)*GV%H_to_Z)*GV%Rlay(k) + H_here(i) = H_here(i) + h(i,j,k) + dz_here(i) = dz_here(i) + dz_2d(i,k) + HxR_here(i) = HxR_here(i) + h(i,j,k)*GV%Rlay(k) endif enddo ; enddo do i=is,ie ; if (H_here(i) > 0.0) then Hf(kf(i),i) = H_here(i) ; Rf(kf(i),i) = HxR_here(i) / H_here(i) + dzf(kf(i),i) = dz_here(i) endif ; enddo endif @@ -892,16 +1004,21 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s if (use_EOS) then pres(1) = 0.0 ; H_top(1) = 0.0 do K=2,kf(i) - pres(K) = pres(K-1) + Z_to_pres*Hf(k-1,i) + pres(K) = pres(K-1) + H_to_pres*Hf(k-1,i) T_int(K) = 0.5*(Tf(k,i)+Tf(k-1,i)) S_int(K) = 0.5*(Sf(k,i)+Sf(k-1,i)) H_top(K) = H_top(K-1) + Hf(k-1,i) enddo - call calculate_density_derivs(T_int, S_int, pres, drho_dT, drho_dS, & - tv%eqn_of_state, (/2,kf(i)/) ) + if (nonBous) then + call calculate_specific_vol_derivs(T_int, S_int, pres, dSpV_dT, dSpV_dS, & + tv%eqn_of_state, (/2,kf(i)/) ) + else + call calculate_density_derivs(T_int, S_int, pres, drho_dT, drho_dS, & + tv%eqn_of_state, (/2,kf(i)/) ) + endif ! Sum the reduced gravities to find out how small a density difference is negligibly small. - drxh_sum = 0.0 + drxh_sum = 0.0 ; dSpVxh_sum = 0.0 if (better_est) then ! This is an estimate that is correct for the non-EBT mode for 2 or 3 layers, or for ! clusters of massless layers at interfaces that can be grouped into 2 or 3 layers. @@ -910,33 +1027,57 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s if (H_top(kf(i)) > 0.0) then I_Htot = 1.0 / (H_top(kf(i)) + Hf(kf(i),i)) ! = 1.0 / (H_top(K) + H_bot(K)) for all K. H_bot(kf(i)+1) = 0.0 - do K=kf(i),2,-1 - H_bot(K) = H_bot(K+1) + Hf(k,i) - drxh_sum = drxh_sum + ((H_top(K) * H_bot(K)) * I_Htot) * & - max(0.0, drho_dT(K)*(Tf(k,i)-Tf(k-1,i)) + drho_dS(K)*(Sf(k,i)-Sf(k-1,i))) - enddo + if (nonBous) then + do K=kf(i),2,-1 + H_bot(K) = H_bot(K+1) + Hf(k,i) + dSpVxh_sum = dSpVxh_sum + ((H_top(K) * H_bot(K)) * I_Htot) * & + min(0.0, dSpV_dT(K)*(Tf(k,i)-Tf(k-1,i)) + dSpV_dS(K)*(Sf(k,i)-Sf(k-1,i))) + enddo + else + do K=kf(i),2,-1 + H_bot(K) = H_bot(K+1) + Hf(k,i) + drxh_sum = drxh_sum + ((H_top(K) * H_bot(K)) * I_Htot) * & + max(0.0, drho_dT(K)*(Tf(k,i)-Tf(k-1,i)) + drho_dS(K)*(Sf(k,i)-Sf(k-1,i))) + enddo + endif endif else ! This estimate is problematic in that it goes like 1/nz for a large number of layers, ! but it is an overestimate (as desired) for a small number of layers, by at a factor ! of (H1+H2)**2/(H1*H2) >= 4 for two thick layers. - do K=2,kf(i) - drxh_sum = drxh_sum + 0.5*(Hf(k-1,i)+Hf(k,i)) * & - max(0.0, drho_dT(K)*(Tf(k,i)-Tf(k-1,i)) + drho_dS(K)*(Sf(k,i)-Sf(k-1,i))) - enddo + if (nonBous) then + do K=2,kf(i) + dSpVxh_sum = dSpVxh_sum + 0.5*(Hf(k-1,i)+Hf(k,i)) * & + min(0.0, dSpV_dT(K)*(Tf(k,i)-Tf(k-1,i)) + dSpV_dS(K)*(Sf(k,i)-Sf(k-1,i))) + enddo + else + do K=2,kf(i) + drxh_sum = drxh_sum + 0.5*(Hf(k-1,i)+Hf(k,i)) * & + max(0.0, drho_dT(K)*(Tf(k,i)-Tf(k-1,i)) + drho_dS(K)*(Sf(k,i)-Sf(k-1,i))) + enddo + endif endif - else - drxh_sum = 0.0 + cg1_est = g_Rho0 * drxh_sum + else ! Not use_EOS + drxh_sum = 0.0 ; dSpVxh_sum = 0.0 if (better_est) then H_top(1) = 0.0 do K=2,kf(i) ; H_top(K) = H_top(K-1) + Hf(k-1,i) ; enddo if (H_top(kf(i)) > 0.0) then I_Htot = 1.0 / (H_top(kf(i)) + Hf(kf(i),i)) ! = 1.0 / (H_top(K) + H_bot(K)) for all K. H_bot(kf(i)+1) = 0.0 - do K=kf(i),2,-1 - H_bot(K) = H_bot(K+1) + Hf(k,i) - drxh_sum = drxh_sum + ((H_top(K) * H_bot(K)) * I_Htot) * max(0.0,Rf(k,i)-Rf(k-1,i)) - enddo + if (nonBous) then + do K=kf(i),2,-1 + H_bot(K) = H_bot(K+1) + Hf(k,i) + dSpVxh_sum = dSpVxh_sum + ((H_top(K) * H_bot(K)) * I_Htot) * & + min(0.0, (Rf(k-1,i)-Rf(k,i)) / (Rf(k,i)*Rf(k-1,i))) + enddo + else + do K=kf(i),2,-1 + H_bot(K) = H_bot(K+1) + Hf(k,i) + drxh_sum = drxh_sum + ((H_top(K) * H_bot(K)) * I_Htot) * max(0.0,Rf(k,i)-Rf(k-1,i)) + enddo + endif endif else do K=2,kf(i) @@ -945,19 +1086,32 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s endif endif + if (nonBous) then + ! Note that dSpVxh_sum is negative for stable stratification. + cg1_est = H_to_pres * abs(dSpVxh_sum) + else + cg1_est = g_Rho0 * drxh_sum + endif + ! Find gprime across each internal interface, taking care of convective ! instabilities by merging layers. - if (g_Rho0 * drxh_sum > cg1_min2) then + if (cg1_est > cg1_min2) then ! Merge layers to eliminate convective instabilities or exceedingly ! small reduced gravities. Merging layers reduces the estimated wave speed by ! (rho(2)-rho(1))*h(1)*h(2) / H_tot. if (use_EOS) then kc = 1 - Hc(1) = Hf(1,i) ; Tc(1) = Tf(1,i) ; Sc(1) = Sf(1,i) + Hc(1) = Hf(1,i) ; dzc(1) = dzf(1,i) ; Tc(1) = Tf(1,i) ; Sc(1) = Sf(1,i) do k=2,kf(i) - if (better_est) then + if (better_est .and. nonBous) then + merge = ((dSpV_dT(K)*(Tc(kc)-Tf(k,i)) + dSpV_dS(K)*(Sc(kc)-Sf(k,i))) * & + ((Hc(kc) * Hf(k,i))*I_Htot) < abs(2.0 * tol_merge * dSpVxh_sum)) + elseif (better_est) then merge = ((drho_dT(K)*(Tf(k,i)-Tc(kc)) + drho_dS(K)*(Sf(k,i)-Sc(kc))) * & ((Hc(kc) * Hf(k,i))*I_Htot) < 2.0 * tol_merge*drxh_sum) + elseif (nonBous) then + merge = ((dSpV_dT(K)*(Tc(kc)-Tf(k,i)) + dSpV_dS(K)*(Sc(kc)-Sf(k,i))) * & + (Hc(kc) + Hf(k,i)) < abs(2.0 * tol_merge * dSpVxh_sum)) else merge = ((drho_dT(K)*(Tf(k,i)-Tc(kc)) + drho_dS(K)*(Sf(k,i)-Sc(kc))) * & (Hc(kc) + Hf(k,i)) < 2.0 * tol_merge*drxh_sum) @@ -967,14 +1121,21 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s I_Hnew = 1.0 / (Hc(kc) + Hf(k,i)) Tc(kc) = (Hc(kc)*Tc(kc) + Hf(k,i)*Tf(k,i)) * I_Hnew Sc(kc) = (Hc(kc)*Sc(kc) + Hf(k,i)*Sf(k,i)) * I_Hnew - Hc(kc) = (Hc(kc) + Hf(k,i)) + Hc(kc) = Hc(kc) + Hf(k,i) + dzc(kc) = dzc(kc) + dzf(k,i) ! Backtrack to remove any convective instabilities above... Note ! that the tolerance is a factor of two larger, to avoid limit how ! far back we go. do K2=kc,2,-1 - if (better_est) then + if (better_est .and. nonBous) then + merge = ( (dSpV_dT(K2)*(Tc(k2-1)-Tc(k2)) + dSpV_dS(K2)*(Sc(k2-1)-Sc(k2))) * & + ((Hc(k2) * Hc(k2-1))*I_Htot) < abs(tol_merge * dSpVxh_sum) ) + elseif (better_est) then merge = ((drho_dT(K2)*(Tc(k2)-Tc(k2-1)) + drho_dS(K2)*(Sc(k2)-Sc(k2-1))) * & ((Hc(k2) * Hc(k2-1))*I_Htot) < tol_merge*drxh_sum) + elseif (nonBous) then + merge = ( (dSpV_dT(K2)*(Tc(k2-1)-Tc(k2)) + dSpV_dS(K2)*(Sc(k2-1)-Sc(k2))) * & + (Hc(k2) + Hc(k2-1)) < abs(tol_merge * dSpVxh_sum) ) else merge = ((drho_dT(K2)*(Tc(k2)-Tc(k2-1)) + drho_dS(K2)*(Sc(k2)-Sc(k2-1))) * & (Hc(k2) + Hc(k2-1)) < tol_merge*drxh_sum) @@ -985,27 +1146,44 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s Tc(kc-1) = (Hc(kc)*Tc(kc) + Hc(kc-1)*Tc(kc-1)) * I_Hnew Sc(kc-1) = (Hc(kc)*Sc(kc) + Hc(kc-1)*Sc(kc-1)) * I_Hnew Hc(kc-1) = Hc(kc) + Hc(kc-1) + dzc(kc-1) = dzc(kc) + dzc(kc-1) kc = kc - 1 else ; exit ; endif enddo else ! Add a new layer to the column. kc = kc + 1 - drho_dS(Kc) = drho_dS(K) ; drho_dT(Kc) = drho_dT(K) - Tc(kc) = Tf(k,i) ; Sc(kc) = Sf(k,i) ; Hc(kc) = Hf(k,i) + if (nonBous) then + dSpV_dS(Kc) = dSpV_dS(K) ; dSpV_dT(Kc) = dSpV_dT(K) + else + drho_dS(Kc) = drho_dS(K) ; drho_dT(Kc) = drho_dT(K) + endif + Tc(kc) = Tf(k,i) ; Sc(kc) = Sf(k,i) ; Hc(kc) = Hf(k,i) ; dzc(kc) = dzf(k,i) endif enddo ! At this point there are kc layers and the gprimes should be positive. - do K=2,kc ! Revisit this if non-Boussinesq. - gprime(K) = g_Rho0 * (drho_dT(K)*(Tc(k)-Tc(k-1)) + drho_dS(K)*(Sc(k)-Sc(k-1))) - enddo - else ! .not.use_EOS + if (nonBous) then + do K=2,kc + gprime(K) = H_to_pres * (dSpV_dT(K)*(Tc(k-1)-Tc(k)) + dSpV_dS(K)*(Sc(k-1)-Sc(k))) + enddo + else + do K=2,kc + gprime(K) = g_Rho0 * (drho_dT(K)*(Tc(k)-Tc(k-1)) + drho_dS(K)*(Sc(k)-Sc(k-1))) + enddo + endif + else ! .not. (use_EOS) ! Do the same with density directly... kc = 1 - Hc(1) = Hf(1,i) ; Rc(1) = Rf(1,i) + Hc(1) = Hf(1,i) ; dzc(1) = dzf(1,i) ; Rc(1) = Rf(1,i) do k=2,kf(i) - if (better_est) then - merge = ((Rf(k,i) - Rc(kc)) * ((Hc(kc) * Hf(k,i))*I_Htot) < 2.0 * tol_merge*drxh_sum) + if (nonBous .and. better_est) then + merge = ((Rf(k,i) - Rc(kc)) * ((Hc(kc) * Hf(k,i))*I_Htot) < & + (Rc(kc)*Rf(k,i)) * abs(2.0 * tol_merge * dSpVxh_sum)) + elseif (nonBous) then + merge = ((Rf(k,i) - Rc(kc)) * (Hc(kc) + Hf(k,i)) < & + (Rc(kc)*Rf(k,i)) * abs(2.0 * tol_merge * dSpVxh_sum)) + elseif (better_est) then + merge = ((Rf(k,i) - Rc(kc)) * ((Hc(kc) * Hf(k,i))*I_Htot) < 2.0*tol_merge*drxh_sum) else merge = ((Rf(k,i) - Rc(kc)) * (Hc(kc) + Hf(k,i)) < 2.0*tol_merge*drxh_sum) endif @@ -1013,6 +1191,7 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s ! Merge this layer with the one above and backtrack. Rc(kc) = (Hc(kc)*Rc(kc) + Hf(k,i)*Rf(k,i)) / (Hc(kc) + Hf(k,i)) Hc(kc) = Hc(kc) + Hf(k,i) + dzc(kc) = dzc(kc) + dzf(k,i) ! Backtrack to remove any convective instabilities above... Note ! that the tolerance is a factor of two larger, to avoid limit how ! far back we go. @@ -1026,19 +1205,26 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s ! Merge the two bottommost layers. At this point kc = k2. Rc(kc-1) = (Hc(kc)*Rc(kc) + Hc(kc-1)*Rc(kc-1)) / (Hc(kc) + Hc(kc-1)) Hc(kc-1) = Hc(kc) + Hc(kc-1) + dzc(kc-1) = dzc(kc) + dzc(kc-1) kc = kc - 1 else ; exit ; endif enddo else ! Add a new layer to the column. kc = kc + 1 - Rc(kc) = Rf(k,i) ; Hc(kc) = Hf(k,i) + Rc(kc) = Rf(k,i) ; Hc(kc) = Hf(k,i) ; dzc(kc) = dzf(k,i) endif enddo ! At this point there are kc layers and the gprimes should be positive. - do K=2,kc ! Revisit this if non-Boussinesq. - gprime(K) = g_Rho0 * (Rc(k)-Rc(k-1)) - enddo + if (nonBous) then + do K=2,kc + gprime(K) = H_to_pres * (Rc(k) - Rc(k-1)) / (Rc(k) * Rc(k-1)) + enddo + else + do K=2,kc + gprime(K) = g_Rho0 * (Rc(k)-Rc(k-1)) + enddo + endif endif ! use_EOS !-----------------NOW FIND WAVE SPEEDS--------------------------------------- @@ -1063,8 +1249,13 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s N2(:) = 0. do K=2,kc - Igl(K) = 1.0/(gprime(K)*Hc(k)) ; Igu(K) = 1.0/(gprime(K)*Hc(k-1)) - N2(K) = US%L_to_Z**2*gprime(K)/(0.5*(Hc(k)+Hc(k-1))) + Igl(K) = 1.0 / (gprime(K)*Hc(k)) ; Igu(K) = 1.0 / (gprime(K)*Hc(k-1)) + if (nonBous) then + N2(K) = 2.0*US%L_to_Z**2*gprime(K) * (Hc(k) + Hc(k-1)) / & ! Units are [T-2 ~> s-2] + (dzc(k) + dzc(k-1))**2 + else + N2(K) = 2.0*US%L_to_Z**2*GV%Z_to_H*gprime(K) / (dzc(k) + dzc(k-1)) ! Units are [T-2 ~> s-2] + endif if (better_est) then speed2_tot = speed2_tot + gprime(K)*((H_top(K) * H_bot(K)) * I_Htot) else @@ -1113,12 +1304,11 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s ! renormalization of the integral of the profile w2avg = 0.0 do k=1,kc - w2avg = w2avg + 0.5*(mode_struct(K)**2+mode_struct(K+1)**2)*Hc(k) ![Z L4 T-4] + w2avg = w2avg + 0.5*(mode_struct(K)**2+mode_struct(K+1)**2)*Hc(k) ! [H L4 T-4] enddo renorm = sqrt(htot(i)*a_int/w2avg) ! [T2 L-2] do K=1,kc+1 ; mode_struct(K) = renorm * mode_struct(K) ; enddo ! after renorm, mode_struct is again [nondim] - if (abs(dlam) < tol_solve*lam_1) exit enddo @@ -1131,7 +1321,7 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s ! vertical derivative of w at interfaces lives on the layer points do k=1,kc - mode_struct_fder(k) = (mode_struct(k) - mode_struct(k+1)) / Hc(k) + mode_struct_fder(k) = (mode_struct(k) - mode_struct(k+1)) / dzc(k) enddo ! boundary condition for derivative is no-gradient @@ -1163,18 +1353,12 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s mode_struct_sq(K+1)*N2(K+1)) * Hc(k) enddo - ! Note that remapping_core_h requires that the same units be used - ! for both the source and target grid thicknesses, here [H ~> m or kg m-2]. - do k = 1,kc - Hc_H(k) = GV%Z_to_H * Hc(k) - enddo - ! for w (diag) interpolate onto all interfaces - call interpolate_column(kc, Hc_H(1:kc), mode_struct(1:kc+1), & + call interpolate_column(kc, Hc(1:kc), mode_struct(1:kc+1), & nz, h(i,j,:), modal_structure(:), .false.) ! for u (remap) onto all layers - call remapping_core_h(CS%remapping_CS, kc, Hc_H(1:kc), mode_struct_fder(1:kc), & + call remapping_core_h(CS%remapping_CS, kc, Hc(1:kc), mode_struct_fder(1:kc), & nz, h(i,j,:), modal_structure_fder(:), & GV%H_subroundoff, GV%H_subroundoff) @@ -1313,7 +1497,7 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s ! derivative of vertical profile (i.e. dw/dz) is evaluated at the layer point do k=1,kc - mode_struct_fder(k) = (mode_struct(k) - mode_struct(k+1)) / Hc(k) + mode_struct_fder(k) = (mode_struct(k) - mode_struct(k+1)) / dzc(k) enddo ! boundary condition for 1st derivative is no-gradient @@ -1345,18 +1529,12 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, w_struct, u_struct, u_s mode_struct_sq(K+1)*N2(K+1)) * Hc(k) enddo - ! Note that remapping_core_h requires that the same units be used - ! for both the source and target grid thicknesses, here [H ~> m or kg m-2]. - do k = 1,kc - Hc_H(k) = GV%Z_to_H * Hc(k) - enddo - ! for w (diag) interpolate onto all interfaces - call interpolate_column(kc, Hc_H(1:kc), mode_struct(1:kc+1), & + call interpolate_column(kc, Hc(1:kc), mode_struct(1:kc+1), & nz, h(i,j,:), modal_structure(:), .false.) ! for u (remap) onto all layers - call remapping_core_h(CS%remapping_CS, kc, Hc_H(1:kc), mode_struct_fder(1:kc), & + call remapping_core_h(CS%remapping_CS, kc, Hc(1:kc), mode_struct_fder(1:kc), & nz, h(i,j,:), modal_structure_fder(:), & GV%H_subroundoff, GV%H_subroundoff) diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index 83910e6690..172d2459d5 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -110,11 +110,11 @@ module MOM_internal_tides real, allocatable, dimension(:,:,:) :: u_struct_bot !< Bottom value of u_struct, !! for each mode [Z-1 ~> m-1] real, allocatable, dimension(:,:,:) :: int_w2 !< Vertical integral of w_struct squared, - !! for each mode [Z ~> m] + !! for each mode [H ~> m or kg m-2] real, allocatable, dimension(:,:,:) :: int_U2 !< Vertical integral of u_struct squared, - !! for each mode [Z-1 ~> m-1] + !! for each mode [H Z-2 ~> m-1 or kg m-4] real, allocatable, dimension(:,:,:) :: int_N2w2 !< Depth-integrated Brunt Vaissalla freqency times - !! vertical profile squared, for each mode [Z T-2 ~> m s-2] + !! vertical profile squared, for each mode [H T-2 ~> m s-2 or kg m-2 s-2] real :: q_itides !< fraction of local dissipation [nondim] real :: En_sum !< global sum of energy for use in debugging, in MKS units [J] type(time_type), pointer :: Time => NULL() !< A pointer to the model's clock. @@ -259,6 +259,7 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, real :: En_initial, Delta_E_check ! Energies for debugging [R Z3 T-2 ~> J m-2] real :: TKE_Froude_loss_check, TKE_Froude_loss_tot ! Energy losses for debugging [R Z3 T-3 ~> W m-2] character(len=160) :: mesg ! The text of an error message + integer :: En_halo_ij_stencil ! The halo size needed for energy advection integer :: a, m, fr, i, j, k, is, ie, js, je, isd, ied, jsd, jed, nAngle integer :: id_g, jd_g ! global (decomp-invar) indices (for debugging) type(group_pass_type), save :: pass_test, pass_En @@ -314,13 +315,16 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, else call wave_speeds(h, tv, G, GV, US, CS%nMode, cn, CS%wave_speed, & CS%w_struct, CS%u_struct, CS%u_struct_max, CS%u_struct_bot, & - Nb, CS%int_w2, CS%int_U2, CS%int_N2w2, full_halos=.true.) + Nb, CS%int_w2, CS%int_U2, CS%int_N2w2, halo_size=2) + ! The value of halo_size above would have to be larger if there were + ! not a halo update between the calls to propagate_x and propagate_y. + ! It can be 1 point smaller if teleport is not used. endif ! Set the wave speeds for the modes, using cg(n) ~ cg(1)/n.********************** ! This is wrong, of course, but it works reasonably in some cases. ! Uncomment if wave_speed is not used to calculate the true values (BDM). - !do m=1,CS%nMode ; do j=jsd,jed ; do i=isd,ied + !do m=1,CS%nMode ; do j=js-2,je+2 ; do i=is-2,ie+2 ! cn(i,j,m) = cn(i,j,1) / real(m) !enddo ; enddo ; enddo @@ -362,6 +366,7 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, call refract(CS%En(:,:,:,fr,m), cn(:,:,m), CS%frequency(fr), 0.5*dt, & G, US, CS%nAngle, CS%use_PPMang) enddo ; enddo + ! A this point, CS%En is only valid on the computational domain. ! Check for En<0 - for debugging, delete later do m=1,CS%nMode ; do fr=1,CS%Nfreq ; do a=1,CS%nAngle @@ -381,8 +386,13 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, call complete_group_pass(pass_test, G%domain) + ! Set the halo size to work on, using similar logic to that used in propagate. This may need + ! to be adjusted depending on the advection scheme and whether teleport is used. + if (CS%upwind_1st) then ; En_halo_ij_stencil = 2 + else ; En_halo_ij_stencil = 3 ; endif + ! Rotate points in the halos as necessary. - call correct_halo_rotation(CS%En, test, G, CS%nAngle) + call correct_halo_rotation(CS%En, test, G, CS%nAngle, halo=En_halo_ij_stencil) ! Propagate the waves. do m=1,CS%nMode ; do fr=1,CS%Nfreq @@ -414,6 +424,7 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, call refract(CS%En(:,:,:,fr,m), cn(:,:,m), CS%frequency(fr), 0.5*dt, & G, US, CS%NAngle, CS%use_PPMang) enddo ; enddo + ! A this point, CS%En is only valid on the computational domain. ! Check for En<0 - for debugging, delete later do m=1,CS%nMode ; do fr=1,CS%Nfreq ; do a=1,CS%nAngle @@ -421,7 +432,7 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, if (CS%En(i,j,a,fr,m)<0.0) then id_g = i + G%idg_offset ; jd_g = j + G%jdg_offset ! for debugging write(mesg,*) 'After second refraction: En<0.0 at ig=', id_g, ', jg=', jd_g, & - 'En=',CS%En(i,j,a,fr,m) + 'En=', CS%En(i,j,a,fr,m) call MOM_error(WARNING, "propagate_int_tide: "//trim(mesg)) CS%En(i,j,a,fr,m) = 0.0 ! call MOM_error(FATAL, "propagate_int_tide: stopped due to negative energy.") @@ -436,7 +447,7 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, tot_En(:,:) = 0.0 tot_En_mode(:,:,:,:) = 0.0 do m=1,CS%nMode ; do fr=1,CS%Nfreq - do j=jsd,jed ; do i=isd,ied ; do a=1,CS%nAngle + do j=js,je ; do i=is,ie ; do a=1,CS%nAngle tot_En(i,j) = tot_En(i,j) + CS%En(i,j,a,fr,m) tot_En_mode(i,j,fr,m) = tot_En_mode(i,j,fr,m) + CS%En(i,j,a,fr,m) enddo ; enddo ; enddo @@ -445,7 +456,7 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, ! Extract the energy for mixing due to misc. processes (background leakage)------ if (CS%apply_background_drag) then - do m=1,CS%nMode ; do fr=1,CS%nFreq ; do a=1,CS%nAngle ; do j=jsd,jed ; do i=isd,ied + do m=1,CS%nMode ; do fr=1,CS%nFreq ; do a=1,CS%nAngle ; do j=js,je ; do i=is,ie ! Calculate loss rate and apply loss over the time step ; apply the same drag timescale ! to each En component (technically not correct; fix later) CS%TKE_leak_loss(i,j,a,fr,m) = CS%En(i,j,a,fr,m) * CS%decay_rate ! loss rate [R Z3 T-3 ~> W m-2] @@ -468,25 +479,25 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, ! Extract the energy for mixing due to bottom drag------------------------------- if (CS%apply_bottom_drag) then - do j=jsd,jed ; do i=isd,ied ; htot(i,j) = 0.0 ; enddo ; enddo - do k=1,GV%ke ; do j=jsd,jed ; do i=isd,ied + do j=js,je ; do i=is,ie ; htot(i,j) = 0.0 ; enddo ; enddo + do k=1,GV%ke ; do j=js,je ; do i=is,ie htot(i,j) = htot(i,j) + h(i,j,k) enddo ; enddo ; enddo if (GV%Boussinesq) then ! This is mathematically equivalent to the form in the option below, but they differ at roundoff. - do j=jsd,jed ; do i=isd,ied + do j=js,je ; do i=is,ie I_D_here = 1.0 / (max(htot(i,j), CS%drag_min_depth)) drag_scale(i,j) = CS%cdrag * sqrt(max(0.0, US%L_to_Z**2*vel_btTide(i,j)**2 + & tot_En(i,j) * GV%RZ_to_H * I_D_here)) * GV%Z_to_H*I_D_here enddo ; enddo else - do j=jsd,jed ; do i=isd,ied + do j=js,je ; do i=is,ie I_mass = GV%RZ_to_H / (max(htot(i,j), CS%drag_min_depth)) drag_scale(i,j) = (CS%cdrag * (Rho_bot(i,j)*I_mass)) * & sqrt(max(0.0, US%L_to_Z**2*vel_btTide(i,j)**2 + tot_En(i,j) * I_mass)) enddo ; enddo endif - do m=1,CS%nMode ; do fr=1,CS%nFreq ; do a=1,CS%nAngle ; do j=jsd,jed ; do i=isd,ied + do m=1,CS%nMode ; do fr=1,CS%nFreq ; do a=1,CS%nAngle ; do j=js,je ; do i=is,ie ! Calculate loss rate and apply loss over the time step ; apply the same drag timescale ! to each En component (technically not correct; fix later) CS%TKE_quad_loss(i,j,a,fr,m) = CS%En(i,j,a,fr,m) * drag_scale(i,j) ! loss rate @@ -515,7 +526,7 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, do m=1,CS%nMode ; do fr=1,CS%Nfreq ! compute near-bottom and max horizontal baroclinic velocity values at each point - do j=jsd,jed ; do i=isd,ied + do j=js,je ; do i=is,ie id_g = i + G%idg_offset ; jd_g = j + G%jdg_offset ! for debugging ! Calculate wavenumber magnitude @@ -529,9 +540,9 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, ! Back-calculate amplitude from energy equation if ( (G%mask2dT(i,j) > 0.5) .and. (freq2*Kmag2 > 0.0)) then ! Units here are [R Z ~> kg m-2] - KE_term = 0.25*GV%Rho0*( ((freq2 + f2) / (freq2*Kmag2))*US%L_to_Z**2*CS%int_U2(i,j,m) + & + KE_term = 0.25*GV%H_to_RZ*( ((freq2 + f2) / (freq2*Kmag2))*US%L_to_Z**2*CS%int_U2(i,j,m) + & CS%int_w2(i,j,m) ) - PE_term = 0.25*GV%Rho0*( CS%int_N2w2(i,j,m) / freq2 ) + PE_term = 0.25*GV%H_to_RZ*( CS%int_N2w2(i,j,m) / freq2 ) if (KE_term + PE_term > 0.0) then W0 = sqrt( tot_En_mode(i,j,fr,m) / (KE_term + PE_term) ) @@ -557,7 +568,7 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, if (CS%apply_wave_drag) then ! Calculate loss rate and apply loss over the time step call itidal_lowmode_loss(G, GV, US, CS, Nb, Rho_bot, Ub, CS%En, CS%TKE_itidal_loss_fixed, & - CS%TKE_itidal_loss, dt, full_halos=.false.) + CS%TKE_itidal_loss, dt, halo_size=0) endif ! Check for En<0 - for debugging, delete later do m=1,CS%nMode ; do fr=1,CS%Nfreq ; do a=1,CS%nAngle @@ -644,7 +655,7 @@ subroutine propagate_int_tide(h, tv, TKE_itidal_input, vel_btTide, Nb, Rho_bot, ! loss from residual of reflection/transmission coefficients if (CS%apply_residual_drag) then - do m=1,CS%nMode ; do fr=1,CS%nFreq ; do a=1,CS%nAngle ; do j=jsd,jed ; do i=isd,ied + do m=1,CS%nMode ; do fr=1,CS%nFreq ; do a=1,CS%nAngle ; do j=js,je ; do i=is,ie ! implicit form !CS%En(i,j,a,fr,m) = CS%En(i,j,a,fr,m) / (1.0 + dt * CS%TKE_residual_loss(i,j,a,fr,m) / & ! (CS%En(i,j,a,fr,m) + en_subRO)) @@ -863,7 +874,7 @@ end subroutine sum_En !> Calculates the energy lost from the propagating internal tide due to !! scattering over small-scale roughness along the lines of Jayne & St. Laurent (2001). -subroutine itidal_lowmode_loss(G, GV, US, CS, Nb, Rho_bot, Ub, En, TKE_loss_fixed, TKE_loss, dt, full_halos) +subroutine itidal_lowmode_loss(G, GV, US, CS, Nb, Rho_bot, Ub, En, TKE_loss_fixed, TKE_loss, dt, halo_size) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -884,10 +895,9 @@ subroutine itidal_lowmode_loss(G, GV, US, CS, Nb, Rho_bot, Ub, En, TKE_loss_fixe intent(out) :: TKE_loss !< Energy loss rate [R Z3 T-3 ~> W m-2] !! (q*rho*kappa*h^2*N*U^2). real, intent(in) :: dt !< Time increment [T ~> s]. - logical,optional, intent(in) :: full_halos !< If true, do the calculation over the - !! entire computational domain. + integer, optional, intent(in) :: halo_size !< The halo size over which to do the calculations ! Local variables - integer :: j,i,m,fr,a, is, ie, js, je + integer :: j, i, m, fr, a, is, ie, js, je, halo real :: En_tot ! energy for a given mode, frequency, and point summed over angles [R Z3 T-2 ~> J m-2] real :: TKE_loss_tot ! dissipation for a given mode, frequency, and point summed over angles [R Z3 T-3 ~> W m-2] real :: frac_per_sector ! fraction of energy in each wedge [nondim] @@ -901,9 +911,10 @@ subroutine itidal_lowmode_loss(G, GV, US, CS, Nb, Rho_bot, Ub, En, TKE_loss_fixe q_itides = CS%q_itides En_negl = 1e-30*US%kg_m3_to_R*US%m_to_Z**3*US%T_to_s**2 - if (present(full_halos)) then ; if (full_halos) then - is = G%isd ; ie = G%ied ; js = G%jsd ; je = G%jed - endif ; endif + if (present(halo_size)) then + halo = halo_size + is = G%isc - halo ; ie = G%iec + halo ; js = G%jsc - halo ; je = G%jec + halo + endif do j=js,je ; do i=is,ie ; do m=1,CS%nMode ; do fr=1,CS%nFreq @@ -931,7 +942,9 @@ subroutine itidal_lowmode_loss(G, GV, US, CS, Nb, Rho_bot, Ub, En, TKE_loss_fixe enddo else ! no loss if no energy - TKE_loss(i,j,:,fr,m) = 0.0 + do a=1,CS%nAngle + TKE_loss(i,j,a,fr,m) = 0.0 + enddo endif ! Update energy remaining (this is the old explicit calc) @@ -2099,7 +2112,7 @@ end subroutine teleport !> Rotates points in the halos where required to accommodate !! changes in grid orientation, such as at the tripolar fold. -subroutine correct_halo_rotation(En, test, G, NAngle) +subroutine correct_halo_rotation(En, test, G, NAngle, halo) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure real, dimension(:,:,:,:,:), intent(inout) :: En !< The internal gravity wave energy density as a !! function of space, angular orientation, frequency, @@ -2110,18 +2123,19 @@ subroutine correct_halo_rotation(En, test, G, NAngle) !! wave energies in the halo region to be corrected [nondim]. integer, intent(in) :: NAngle !< The number of wave orientations in the !! discretized wave energy spectrum. + integer, intent(in) :: halo !< The halo size over which to do the calculations ! Local variables real, dimension(G%isd:G%ied,NAngle) :: En2d ! A zonal row of the internal gravity wave energy density ! in a frequency band and mode [R Z3 T-2 ~> J m-2]. integer, dimension(G%isd:G%ied) :: a_shift integer :: i_first, i_last, a_new - integer :: a, i, j, isd, ied, jsd, jed, m, fr + integer :: a, i, j, ish, ieh, jsh, jeh, m, fr character(len=160) :: mesg ! The text of an error message - isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + ish = G%isc-halo ; ieh = G%iec+halo ; jsh = G%jsc-halo ; jeh = G%jec+halo - do j=jsd,jed - i_first = ied+1 ; i_last = isd-1 - do i=isd,ied + do j=jsh,jeh + i_first = ieh+1 ; i_last = ish-1 + do i=ish,ieh a_shift(i) = 0 if (test(i,j,1) /= 1.0) then if (i 0.0) then - CS%refl_pref_logical(i,j) = .true. - endif - enddo - enddo + do j=jsd,jed ; do i=isd,ied + ! flag cells with partial reflection + if ((CS%refl_angle(i,j) /= CS%nullangle) .and. & + (CS%refl_pref(i,j) < 1.0) .and. (CS%refl_pref(i,j) > 0.0)) then + CS%refl_pref_logical(i,j) = .true. + endif + enddo ; enddo ! Read in double-reflective (ridge) tags from file call get_param(param_file, mdl, "REFL_DBL_FILE", refl_dbl_file, & @@ -2776,11 +2788,10 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) if (trim(refl_dbl_file) /= '' ) call MOM_error(FATAL, & "REFL_DBL_FILE: "//trim(filename)//" not found") endif - call pass_var(ridge_temp,G%domain) + call pass_var(ridge_temp, G%domain) allocate(CS%refl_dbl(isd:ied,jsd:jed), source=.false.) - do i=isd,ied ; do j=jsd,jed - if (ridge_temp(i,j) == 1) then; CS%refl_dbl(i,j) = .true. - else ; CS%refl_dbl(i,j) = .false. ; endif + do j=jsd,jed ; do i=isd,ied + CS%refl_dbl(i,j) = (ridge_temp(i,j) == 1) enddo ; enddo ! Read in the transmission coefficient and infer the residual @@ -2797,17 +2808,16 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) "TRANS_FILE: "//trim(filename)//" not found") endif - call pass_var(CS%trans,G%domain) + call pass_var(CS%trans, G%domain) + ! residual allocate(CS%residual(isd:ied,jsd:jed), source=0.0) - do j=jsd,jed - do i=isd,ied - if (CS%refl_pref_logical(i,j)) then - CS%residual(i,j) = 1. - CS%refl_pref(i,j) - CS%trans(i,j) - endif - enddo - enddo - call pass_var(CS%residual,G%domain) + do j=G%jsc,G%jec ; do i=G%isc,G%iec + if (CS%refl_pref_logical(i,j)) then + CS%residual(i,j) = 1. - CS%refl_pref(i,j) - CS%trans(i,j) + endif + enddo ; enddo + call pass_var(CS%residual, G%domain) CS%id_cg1 = register_diag_field('ocean_model', 'cn1', diag%axesT1, & Time, 'First baroclinic mode (eigen) speed', 'm s-1', conversion=US%L_T_to_m_s) @@ -2820,7 +2830,6 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) call MOM_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5) enddo - ! Register maps of reflection parameters CS%id_refl_ang = register_diag_field('ocean_model', 'refl_angle', diag%axesT1, & Time, 'Local angle of coastline/ridge/shelf with respect to equator', 'rad') @@ -2971,19 +2980,19 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) write(var_name, '("Itide_int_w2","_mode",i1)') m write(var_descript, '("integral of w2 for mode ",i1)') m CS%id_int_w2_mode(m) = register_diag_field('ocean_model', var_name, & - diag%axesT1, Time, var_descript, 'm', conversion=US%Z_to_m) + diag%axesT1, Time, var_descript, 'm', conversion=GV%H_to_m) call MOM_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5) write(var_name, '("Itide_int_U2","_mode",i1)') m write(var_descript, '("integral of U2 for mode ",i1)') m CS%id_int_U2_mode(m) = register_diag_field('ocean_model', var_name, & - diag%axesT1, Time, var_descript, 'm-1', conversion=US%m_to_L) + diag%axesT1, Time, var_descript, 'm-1', conversion=US%m_to_Z*GV%H_to_Z) call MOM_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5) write(var_name, '("Itide_int_N2w2","_mode",i1)') m write(var_descript, '("integral of N2w2 for mode ",i1)') m CS%id_int_N2w2_mode(m) = register_diag_field('ocean_model', var_name, & - diag%axesT1, Time, var_descript, 'm s-2', conversion=US%Z_to_m*US%s_to_T**2) + diag%axesT1, Time, var_descript, 'm s-2', conversion=GV%H_to_m*US%s_to_T**2) call MOM_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5) enddo diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 1ccd6a7fb2..5b89c8c726 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -159,6 +159,9 @@ module MOM_diabatic_driver !! evaporated in one time-step [nondim]. integer :: halo_TS_diff = 0 !< The temperature, salinity and thickness halo size that !! must be valid for the diffusivity calculations. + integer :: halo_diabatic = 0 !< The temperature, salinity, specific volume and thickness + !! halo size that must be valid for the diabatic calculations, + !! including vertical mixing and internal tide propagation. logical :: useKPP = .false. !< use CVMix/KPP diffusivities and non-local transport logical :: KPPisPassive !< If true, KPP is in passive mode, not changing answers. logical :: debug !< If true, write verbose checksums for debugging purposes. @@ -2661,7 +2664,7 @@ subroutine extract_diabatic_member(CS, opacity_CSp, optics_CSp, evap_CFL_limit, ! Constants within diabatic_CS if (present(evap_CFL_limit)) evap_CFL_limit = CS%evap_CFL_limit if (present(minimum_forcing_depth)) minimum_forcing_depth = CS%minimum_forcing_depth - if (present(diabatic_halo)) diabatic_halo = CS%halo_TS_diff + if (present(diabatic_halo)) diabatic_halo = CS%halo_diabatic if (present(use_KPP)) use_KPP = CS%use_KPP end subroutine extract_diabatic_member @@ -3513,6 +3516,9 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di halo_TS=CS%halo_TS_diff, double_diffuse=CS%double_diffuse, & physical_OBL_scheme=physical_OBL_scheme) + CS%halo_diabatic = CS%halo_TS_diff + if (CS%use_int_tides) CS%halo_diabatic = max(CS%halo_TS_diff, 2) + if (CS%useKPP .and. (CS%double_diffuse .and. .not.CS%use_CVMix_ddiff)) & call MOM_error(FATAL, 'diabatic_driver_init: DOUBLE_DIFFUSION (old method) does not work '//& 'with KPP. Please set DOUBLE_DIFFUSION=False and USE_CVMIX_DDIFF=True.')