Skip to content

Commit

Permalink
Merge c5a689e into a2aaad6
Browse files Browse the repository at this point in the history
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA authored Oct 25, 2021
2 parents a2aaad6 + c5a689e commit f279c78
Show file tree
Hide file tree
Showing 5 changed files with 55 additions and 134 deletions.
17 changes: 7 additions & 10 deletions src/parameterizations/lateral/MOM_internal_tides.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1905,23 +1905,22 @@ subroutine PPM_reconstruction_x(h_in, h_l, h_r, G, LB, simple_2nd)
real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_l !< Left edge value of reconstruction (2D).
real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_r !< Right edge value of reconstruction (2D).
type(loop_bounds_type), intent(in) :: LB !< A structure with the active loop bounds.
logical, optional, intent(in) :: simple_2nd !< If true, use the arithmetic mean
logical, intent(in) :: simple_2nd !< If true, use the arithmetic mean
!! energy densities as default edge values
!! for a simple 2nd order scheme.
! Local variables
real, dimension(SZI_(G),SZJ_(G)) :: slp ! The slopes.
real, parameter :: oneSixth = 1./6.
real :: h_ip1, h_im1
real :: dMx, dMn
logical :: use_CW84, use_2nd
logical :: use_CW84
character(len=256) :: mesg ! The text of an error message
integer :: i, j, isl, iel, jsl, jel, stencil

use_2nd = .false. ; if (present(simple_2nd)) use_2nd = simple_2nd
isl = LB%ish-1 ; iel = LB%ieh+1 ; jsl = LB%jsh ; jel = LB%jeh

! This is the stencil of the reconstruction, not the scheme overall.
stencil = 2 ; if (use_2nd) stencil = 1
stencil = 2 ; if (simple_2nd) stencil = 1

if ((isl-stencil < G%isd) .or. (iel+stencil > G%ied)) then
write(mesg,'("In MOM_internal_tides, PPM_reconstruction_x called with a ", &
Expand All @@ -1936,7 +1935,7 @@ subroutine PPM_reconstruction_x(h_in, h_l, h_r, G, LB, simple_2nd)
call MOM_error(FATAL,mesg)
endif

if (use_2nd) then
if (simple_2nd) then
do j=jsl,jel ; do i=isl,iel
h_im1 = G%mask2dT(i-1,j) * h_in(i-1,j) + (1.0-G%mask2dT(i-1,j)) * h_in(i,j)
h_ip1 = G%mask2dT(i+1,j) * h_in(i+1,j) + (1.0-G%mask2dT(i+1,j)) * h_in(i,j)
Expand Down Expand Up @@ -1981,23 +1980,21 @@ subroutine PPM_reconstruction_y(h_in, h_l, h_r, G, LB, simple_2nd)
real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_l !< Left edge value of reconstruction (2D).
real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_r !< Right edge value of reconstruction (2D).
type(loop_bounds_type), intent(in) :: LB !< A structure with the active loop bounds.
logical, optional, intent(in) :: simple_2nd !< If true, use the arithmetic mean
logical, intent(in) :: simple_2nd !< If true, use the arithmetic mean
!! energy densities as default edge values
!! for a simple 2nd order scheme.
! Local variables
real, dimension(SZI_(G),SZJ_(G)) :: slp ! The slopes.
real, parameter :: oneSixth = 1./6.
real :: h_jp1, h_jm1
real :: dMx, dMn
logical :: use_2nd
character(len=256) :: mesg ! The text of an error message
integer :: i, j, isl, iel, jsl, jel, stencil

use_2nd = .false. ; if (present(simple_2nd)) use_2nd = simple_2nd
isl = LB%ish ; iel = LB%ieh ; jsl = LB%jsh-1 ; jel = LB%jeh+1

! This is the stencil of the reconstruction, not the scheme overall.
stencil = 2 ; if (use_2nd) stencil = 1
stencil = 2 ; if (simple_2nd) stencil = 1

if ((isl < G%isd) .or. (iel > G%ied)) then
write(mesg,'("In MOM_internal_tides, PPM_reconstruction_y called with a ", &
Expand All @@ -2012,7 +2009,7 @@ subroutine PPM_reconstruction_y(h_in, h_l, h_r, G, LB, simple_2nd)
call MOM_error(FATAL,mesg)
endif

if (use_2nd) then
if (simple_2nd) then
do j=jsl,jel ; do i=isl,iel
h_jm1 = G%mask2dT(i,j-1) * h_in(i,j-1) + (1.0-G%mask2dT(i,j-1)) * h_in(i,j)
h_jp1 = G%mask2dT(i,j+1) * h_in(i,j+1) + (1.0-G%mask2dT(i,j+1)) * h_in(i,j)
Expand Down
22 changes: 11 additions & 11 deletions src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -451,7 +451,7 @@ subroutine calc_slope_functions(h, tv, dt, G, GV, US, CS, OBC)
type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
real, intent(in) :: dt !< Time increment [T ~> s]
type(VarMix_CS), pointer :: CS !< Variable mixing coefficients
type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure.
type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure.
! Local variables
real, dimension(SZI_(G), SZJ_(G),SZK_(GV)+1) :: &
e ! The interface heights relative to mean sea level [Z ~> m].
Expand All @@ -477,10 +477,10 @@ subroutine calc_slope_functions(h, tv, dt, G, GV, US, CS, OBC)
if (CS%use_stored_slopes) then
call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*CS%kappa_smooth, &
CS%slope_x, CS%slope_y, N2_u=N2_u, N2_v=N2_v, halo=1, OBC=OBC)
call calc_Visbeck_coeffs_old(h, CS%slope_x, CS%slope_y, N2_u, N2_v, G, GV, US, CS, OBC=OBC)
call calc_Visbeck_coeffs_old(h, CS%slope_x, CS%slope_y, N2_u, N2_v, G, GV, US, CS, OBC)
else
!call calc_isoneutral_slopes(G, GV, h, e, tv, dt*CS%kappa_smooth, CS%slope_x, CS%slope_y)
call calc_slope_functions_using_just_e(h, G, GV, US, CS, e, .true., OBC=OBC)
call calc_slope_functions_using_just_e(h, G, GV, US, CS, e, .true., OBC)
endif
endif
endif
Expand Down Expand Up @@ -515,7 +515,7 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C
!! at v-points [L2 Z-2 T-2 ~> s-2]
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(VarMix_CS), pointer :: CS !< Variable mixing coefficients
type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure.
type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure.

! Local variables
real :: S2 ! Interface slope squared [nondim]
Expand Down Expand Up @@ -543,10 +543,10 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C

local_open_u_BC = .false.
local_open_v_BC = .false.
if (present(OBC)) then ; if (associated(OBC)) then
if (associated(OBC)) then
local_open_u_BC = OBC%open_u_BCs_exist_globally
local_open_v_BC = OBC%open_v_BCs_exist_globally
endif ; endif
endif

S2max = CS%Visbeck_S_max**2

Expand Down Expand Up @@ -673,8 +673,8 @@ subroutine calc_Eady_growth_rate_2D(CS, G, GV, US, OBC, h, e, dzu, dzv, dzSxN, d
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
type(ocean_OBC_type), pointer, intent(in) :: OBC !< Open boundaries control structure.
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Interface height [Z ~> m]
type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure.
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Interface height [Z ~> m]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: e !< Interface height [Z ~> m]
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: dzu !< dz at u-points [Z ~> m]
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: dzv !< dz at v-points [Z ~> m]
Expand Down Expand Up @@ -859,7 +859,7 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slop
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: e !< Interface position [Z ~> m]
logical, intent(in) :: calculate_slopes !< If true, calculate slopes
!! internally otherwise use slopes stored in CS
type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure.
type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure.
! Local variables
real :: E_x(SZIB_(G), SZJ_(G)) ! X-slope of interface at u points [nondim] (for diagnostics)
real :: E_y(SZI_(G), SZJB_(G)) ! Y-slope of interface at v points [nondim] (for diagnostics)
Expand Down Expand Up @@ -890,10 +890,10 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slop

local_open_u_BC = .false.
local_open_v_BC = .false.
if (present(OBC)) then ; if (associated(OBC)) then
if (associated(OBC)) then
local_open_u_BC = OBC%open_u_BCs_exist_globally
local_open_v_BC = OBC%open_v_BCs_exist_globally
endif ; endif
endif

one_meter = 1.0 * GV%m_to_H
h_neglect = GV%H_subroundoff
Expand Down
43 changes: 16 additions & 27 deletions src/parameterizations/lateral/MOM_thickness_diffuse.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -580,11 +580,11 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
real, intent(in) :: dt !< Time increment [T ~> s]
type(MEKE_type), pointer :: MEKE !< MEKE control structure
type(thickness_diffuse_CS), pointer :: CS !< Control structure for thickness diffusion
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), optional, intent(in) :: int_slope_u !< Ratio that determine how much of
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: int_slope_u !< Ratio that determine how much of
!! the isopycnal slopes are taken directly from
!! the interface slopes without consideration of
!! density gradients [nondim].
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), optional, intent(in) :: int_slope_v !< Ratio that determine how much of
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: int_slope_v !< Ratio that determine how much of
!! the isopycnal slopes are taken directly from
!! the interface slopes without consideration of
!! density gradients [nondim].
Expand Down Expand Up @@ -694,7 +694,6 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV

real, dimension(SZIB_(G), SZJ_(G),SZK_(GV)+1) :: diag_sfn_x, diag_sfn_unlim_x ! Diagnostics
real, dimension(SZI_(G), SZJB_(G),SZK_(GV)+1) :: diag_sfn_y, diag_sfn_unlim_y ! Diagnostics
logical :: present_int_slope_u, present_int_slope_v
logical :: present_slope_x, present_slope_y, calc_derivatives
integer, dimension(2) :: EOSdom_u ! The shifted i-computational domain to use for equation of
! state calculations at u-points.
Expand All @@ -715,8 +714,6 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
N2_floor = CS%N2_floor*US%Z_to_L**2

use_EOS = associated(tv%eqn_of_state)
present_int_slope_u = PRESENT(int_slope_u)
present_int_slope_v = PRESENT(int_slope_v)
present_slope_x = PRESENT(slope_x)
present_slope_y = PRESENT(slope_y)
use_Stanley = CS%Stanley_det_coeff >= 0.
Expand Down Expand Up @@ -818,12 +815,10 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV

EOSdom_u(1) = (is-1) - (G%IsdB-1) ; EOSdom_u(2) = ie - (G%IsdB-1)
!$OMP parallel do default(none) shared(nz,is,ie,js,je,find_work,use_EOS,G,GV,US,pres,T,S, &
!$OMP nk_linear,IsdB,tv,h,h_neglect,e,dz_neglect, &
!$OMP I_slope_max2,h_neglect2,present_int_slope_u, &
!$OMP int_slope_u,KH_u,uhtot,h_frac,h_avail_rsum, &
!$OMP uhD,h_avail,G_scale,Work_u,CS,slope_x,cg1, &
!$OMP diag_sfn_x, diag_sfn_unlim_x,N2_floor,EOSdom_u, &
!$OMP use_stanley, Tsgs2, &
!$OMP nk_linear,IsdB,tv,h,h_neglect,e,dz_neglect,I_slope_max2, &
!$OMP h_neglect2,int_slope_u,KH_u,uhtot,h_frac,h_avail_rsum, &
!$OMP uhD,h_avail,G_scale,Work_u,CS,slope_x,cg1,diag_sfn_x, &
!$OMP diag_sfn_unlim_x,N2_floor,EOSdom_u,use_stanley, Tsgs2, &
!$OMP present_slope_x,G_rho0,Slope_x_PE,hN2_x_PE) &
!$OMP private(drdiA,drdiB,drdkL,drdkR,pres_u,T_u,S_u, &
!$OMP drho_dT_u,drho_dS_u,hg2A,hg2B,hg2L,hg2R,haA, &
Expand Down Expand Up @@ -941,11 +936,9 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV

! Adjust real slope by weights that bias towards slope of interfaces
! that ignore density gradients along layers.
if (present_int_slope_u) then
Slope = (1.0 - int_slope_u(I,j,K)) * Slope + &
int_slope_u(I,j,K) * ((e(i+1,j,K)-e(i,j,K)) * G%IdxCu(I,j))
slope2_Ratio_u(I,K) = (1.0 - int_slope_u(I,j,K)) * slope2_Ratio_u(I,K)
endif
Slope = (1.0 - int_slope_u(I,j,K)) * Slope + &
int_slope_u(I,j,K) * ((e(i+1,j,K)-e(i,j,K)) * G%IdxCu(I,j))
slope2_Ratio_u(I,K) = (1.0 - int_slope_u(I,j,K)) * slope2_Ratio_u(I,K)

Slope_x_PE(I,j,k) = MIN(Slope,CS%slope_max)
hN2_x_PE(I,j,k) = hN2_u(I,K)
Expand Down Expand Up @@ -1090,12 +1083,10 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV
! Calculate the meridional fluxes and gradients.
EOSdom_v(:) = EOS_domain(G%HI)
!$OMP parallel do default(none) shared(nz,is,ie,js,je,find_work,use_EOS,G,GV,US,pres,T,S, &
!$OMP nk_linear,IsdB,tv,h,h_neglect,e,dz_neglect, &
!$OMP I_slope_max2,h_neglect2,present_int_slope_v, &
!$OMP int_slope_v,KH_v,vhtot,h_frac,h_avail_rsum, &
!$OMP vhD,h_avail,G_scale,Work_v,CS,slope_y,cg1, &
!$OMP diag_sfn_y,diag_sfn_unlim_y,N2_floor,EOSdom_v,&
!$OMP use_stanley, Tsgs2, &
!$OMP nk_linear,IsdB,tv,h,h_neglect,e,dz_neglect,I_slope_max2, &
!$OMP h_neglect2,int_slope_v,KH_v,vhtot,h_frac,h_avail_rsum, &
!$OMP vhD,h_avail,G_scale,Work_v,CS,slope_y,cg1,diag_sfn_y, &
!$OMP diag_sfn_unlim_y,N2_floor,EOSdom_v,use_stanley,Tsgs2, &
!$OMP present_slope_y,G_rho0,Slope_y_PE,hN2_y_PE) &
!$OMP private(drdjA,drdjB,drdkL,drdkR,pres_v,T_v,S_v, &
!$OMP drho_dT_v,drho_dS_v,hg2A,hg2B,hg2L,hg2R,haA, &
Expand Down Expand Up @@ -1211,11 +1202,9 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV

! Adjust real slope by weights that bias towards slope of interfaces
! that ignore density gradients along layers.
if (present_int_slope_v) then
Slope = (1.0 - int_slope_v(i,J,K)) * Slope + &
int_slope_v(i,J,K) * ((e(i,j+1,K)-e(i,j,K)) * G%IdyCv(i,J))
slope2_Ratio_v(i,K) = (1.0 - int_slope_v(i,J,K)) * slope2_Ratio_v(i,K)
endif
Slope = (1.0 - int_slope_v(i,J,K)) * Slope + &
int_slope_v(i,J,K) * ((e(i,j+1,K)-e(i,j,K)) * G%IdyCv(i,J))
slope2_Ratio_v(i,K) = (1.0 - int_slope_v(i,J,K)) * slope2_Ratio_v(i,K)

Slope_y_PE(i,J,k) = MIN(Slope,CS%slope_max)
hN2_y_PE(i,J,k) = hN2_v(i,K)
Expand Down
29 changes: 9 additions & 20 deletions src/parameterizations/lateral/MOM_tidal_forcing.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -579,7 +579,7 @@ end subroutine tidal_forcing_sensitivity
!! height. For now, eta and eta_tidal are both geopotential heights in depth
!! units, but probably the input for eta should really be replaced with the
!! column mass anomalies.
subroutine calc_tidal_forcing(Time, eta, eta_tidal, G, CS, deta_tidal_deta, m_to_Z)
subroutine calc_tidal_forcing(Time, eta, eta_tidal, G, CS, m_to_Z)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
type(time_type), intent(in) :: Time !< The time for the caluculation.
real, dimension(SZI_(G),SZJ_(G)), intent(in) :: eta !< The sea surface height anomaly from
Expand All @@ -588,16 +588,12 @@ subroutine calc_tidal_forcing(Time, eta, eta_tidal, G, CS, deta_tidal_deta, m_to
!! anomalies [Z ~> m].
type(tidal_forcing_CS), pointer :: CS !< The control structure returned by a
!! previous call to tidal_forcing_init.
real, optional, intent(out) :: deta_tidal_deta !< The partial derivative of
!! eta_tidal with the local value of
!! eta [nondim].
real, optional, intent(in) :: m_to_Z !< A scaling factor from m to the units of eta.
real, intent(in) :: m_to_Z !< A scaling factor from m to the units of eta.

! Local variables
real :: now ! The relative time in seconds.
real :: amp_cosomegat, amp_sinomegat
real :: cosomegat, sinomegat
real :: m_Z ! A scaling factor from m to depth units.
real :: eta_prop ! The nondimenional constant of proportionality beteen eta and eta_tidal.
integer :: i, j, c, m, is, ie, js, je, Isq, Ieq, Jsq, Jeq
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
Expand All @@ -622,21 +618,14 @@ subroutine calc_tidal_forcing(Time, eta, eta_tidal, G, CS, deta_tidal_deta, m_to
eta_prop = 0.0
endif

if (present(deta_tidal_deta)) then
deta_tidal_deta = eta_prop
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 ; eta_tidal(i,j) = 0.0 ; enddo ; enddo
else
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
eta_tidal(i,j) = eta_prop*eta(i,j)
enddo ; enddo
endif

m_Z = 1.0 ; if (present(m_to_Z)) m_Z = m_to_Z
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
eta_tidal(i,j) = eta_prop*eta(i,j)
enddo ; enddo

do c=1,CS%nc
m = CS%struct(c)
amp_cosomegat = m_Z*CS%amp(c)*CS%love_no(c) * cos(CS%freq(c)*now + CS%phase0(c))
amp_sinomegat = m_Z*CS%amp(c)*CS%love_no(c) * sin(CS%freq(c)*now + CS%phase0(c))
amp_cosomegat = m_to_Z*CS%amp(c)*CS%love_no(c) * cos(CS%freq(c)*now + CS%phase0(c))
amp_sinomegat = m_to_Z*CS%amp(c)*CS%love_no(c) * sin(CS%freq(c)*now + CS%phase0(c))
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
eta_tidal(i,j) = eta_tidal(i,j) + (amp_cosomegat*CS%cos_struct(i,j,m) + &
amp_sinomegat*CS%sin_struct(i,j,m))
Expand All @@ -647,7 +636,7 @@ subroutine calc_tidal_forcing(Time, eta, eta_tidal, G, CS, deta_tidal_deta, m_to
cosomegat = cos(CS%freq(c)*now)
sinomegat = sin(CS%freq(c)*now)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
eta_tidal(i,j) = eta_tidal(i,j) + m_Z*CS%ampsal(i,j,c) * &
eta_tidal(i,j) = eta_tidal(i,j) + m_to_Z*CS%ampsal(i,j,c) * &
(cosomegat*CS%cosphasesal(i,j,c) + sinomegat*CS%sinphasesal(i,j,c))
enddo ; enddo
enddo ; endif
Expand All @@ -656,7 +645,7 @@ subroutine calc_tidal_forcing(Time, eta, eta_tidal, G, CS, deta_tidal_deta, m_to
cosomegat = cos(CS%freq(c)*now)
sinomegat = sin(CS%freq(c)*now)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
eta_tidal(i,j) = eta_tidal(i,j) - m_Z*CS%SAL_SCALAR*CS%amp_prev(i,j,c) * &
eta_tidal(i,j) = eta_tidal(i,j) - m_to_Z*CS%SAL_SCALAR*CS%amp_prev(i,j,c) * &
(cosomegat*CS%cosphase_prev(i,j,c) + sinomegat*CS%sinphase_prev(i,j,c))
enddo ; enddo
enddo ; endif
Expand Down
Loading

0 comments on commit f279c78

Please sign in to comment.