Skip to content

Commit

Permalink
+Add hor_index_type variants of EOS routines
Browse files Browse the repository at this point in the history 
  Added variants of calculate_density routines that use a hor_index_type to
specify array extents and unit_scale_types for dimensional consistency testing,
further overloading existing interfaces.  Also replaced the recently added
rho_scale and pres_scale arguments to int_density_dz with an optional
unit_scale_type argument, and modified calls to use this new argument.  All
answers are bitwise identical, but there are changes to external interfaces.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Apr 10, 2020
1 parent 1252e3b commit ce0905e
Show file tree
Hide file tree
Showing 6 changed files with 627 additions and 184 deletions.
27 changes: 13 additions & 14 deletions src/core/MOM_PressureForce_Montgomery.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -186,8 +186,7 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb
!$OMP parallel do default(shared)
do k=1,nz
call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,k), p(:,:,k+1), &
0.0, G%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=1, &
SV_scale=US%R_to_kg_m3, pres_scale=US%RL2_T2_to_Pa)
0.0, G%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=1, US=US)
enddo
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do k=1,nz ; do i=Isq,Ieq+1
Expand Down Expand Up @@ -659,8 +658,8 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star)
Ihtot(i) = GV%H_to_Z / ((e(i,j,1)-e(i,j,nz+1)) + z_neglect)
press(i) = -Rho0xG*e(i,j,1)
enddo
call calculate_density(tv%T(:,j,1), tv%S(:,j,1), press, rho_in_situ, Isq, Ieq-Isq+2, &
tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density(tv%T(:,j,1), tv%S(:,j,1), press, rho_in_situ, G%HI, &
tv%eqn_of_state, US, halo=1)
do i=Isq,Ieq+1
pbce(i,j,1) = G_Rho0*(GFS_scale * rho_in_situ(i)) * GV%H_to_Z
enddo
Expand All @@ -670,8 +669,8 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star)
T_int(i) = 0.5*(tv%T(i,j,k-1)+tv%T(i,j,k))
S_int(i) = 0.5*(tv%S(i,j,k-1)+tv%S(i,j,k))
enddo
call calculate_density_derivs(T_int, S_int, press, dR_dT, dR_dS, Isq, Ieq-Isq+2, &
tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density_derivs(T_int, S_int, press, dR_dT, dR_dS, G%HI, &
tv%eqn_of_state, US, halo=1)
do i=Isq,Ieq+1
pbce(i,j,k) = pbce(i,j,k-1) + G_Rho0 * &
((e(i,j,K) - e(i,j,nz+1)) * Ihtot(i)) * &
Expand Down Expand Up @@ -755,9 +754,8 @@ subroutine Set_pbce_nonBouss(p, tv, G, GV, US, GFS_scale, pbce, alpha_star)
else
!$OMP parallel do default(shared) private(T_int,S_int,dR_dT,dR_dS,rho_in_situ)
do j=Jsq,Jeq+1
call calculate_density(tv%T(:,j,nz), tv%S(:,j,nz), p(:,j,nz+1), rho_in_situ, &
Isq, Ieq-Isq+2, tv%eqn_of_state, scale=US%kg_m3_to_R, &
pres_scale=US%RL2_T2_to_Pa)
call calculate_density(tv%T(:,j,nz), tv%S(:,j,nz), p(:,j,nz+1), rho_in_situ, G%HI, &
tv%eqn_of_state, US, halo=1)
do i=Isq,Ieq+1
C_htot(i,j) = dP_dH / ((p(i,j,nz+1)-p(i,j,1)) + dp_neglect)
pbce(i,j,nz) = dP_dH / (rho_in_situ(i))
Expand All @@ -767,11 +765,12 @@ subroutine Set_pbce_nonBouss(p, tv, G, GV, US, GFS_scale, pbce, alpha_star)
T_int(i) = 0.5*(tv%T(i,j,k)+tv%T(i,j,k+1))
S_int(i) = 0.5*(tv%S(i,j,k)+tv%S(i,j,k+1))
enddo
call calculate_density(T_int, S_int, p(:,j,k+1), rho_in_situ, Isq, Ieq-Isq+2, tv%eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density_derivs(T_int, S_int, p(:,j,k+1), dR_dT, dR_dS, &
Isq, Ieq-Isq+2, tv%eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density(T_int, S_int, p(:,j,k+1), rho_in_situ, G%HI, tv%eqn_of_state, US, halo=1)
! call calculate_density_derivs(T_int, S_int, p(:,j,k+1), dR_dT, dR_dS, &
! Isq, Ieq-Isq+2, tv%eqn_of_state, &
! scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density_derivs(T_int, S_int, p(:,j,k+1), dR_dT, dR_dS, G%HI, &
tv%eqn_of_state, US, halo=1)
do i=Isq,Ieq+1
pbce(i,j,k) = pbce(i,j,k+1) + ((p(i,j,K+1)-p(i,j,1))*C_htot(i,j)) * &
((dR_dT(i)*(tv%T(i,j,k+1)-tv%T(i,j,k)) + &
Expand Down
39 changes: 18 additions & 21 deletions src/core/MOM_PressureForce_analytic_FV.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -264,24 +264,22 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p
if ( CS%Recon_Scheme == 1 ) then
call int_spec_vol_dp_generic_plm( T_t(:,:,k), T_b(:,:,k), S_t(:,:,k), S_b(:,:,k), &
p(:,:,K), p(:,:,K+1), alpha_ref, dp_neglect, p(:,:,nz+1), G%HI, &
tv%eqn_of_state, dza(:,:,k), intp_dza(:,:,k), &
intx_dza(:,:,k), inty_dza(:,:,k), useMassWghtInterp=CS%useMassWghtInterp, &
SV_scale=US%R_to_kg_m3, pres_scale=US%RL2_T2_to_Pa)
tv%eqn_of_state, dza(:,:,k), intp_dza(:,:,k), intx_dza(:,:,k), inty_dza(:,:,k), &
useMassWghtInterp=CS%useMassWghtInterp, US=US)
elseif ( CS%Recon_Scheme == 2 ) then
call MOM_error(FATAL, "PressureForce_AFV_nonBouss: "//&
"int_spec_vol_dp_generic_ppm does not exist yet.")
! call int_spec_vol_dp_generic_ppm ( tv%T(:,:,k), T_t(:,:,k), T_b(:,:,k), &
! tv%S(:,:,k), S_t(:,:,k), S_b(:,:,k), p(:,:,K), p(:,:,K+1), &
! alpha_ref, G%HI, tv%eqn_of_state, dza(:,:,k), intp_dza(:,:,k), &
! intx_dza(:,:,k), inty_dza(:,:,k))
! intx_dza(:,:,k), inty_dza(:,:,k), US=US)
endif
else
call int_specific_vol_dp(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), p(:,:,K), &
p(:,:,K+1), alpha_ref, G%HI, tv%eqn_of_state, &
dza(:,:,k), intp_dza(:,:,k), intx_dza(:,:,k), &
inty_dza(:,:,k), bathyP=p(:,:,nz+1), dP_tiny=dp_neglect, &
useMassWghtInterp=CS%useMassWghtInterp, &
SV_scale=US%R_to_kg_m3, pres_scale=US%RL2_T2_to_Pa)
useMassWghtInterp=CS%useMassWghtInterp, US=US)
endif
else
alpha_anom = 1.0 / GV%Rlay(k) - alpha_ref
Expand Down Expand Up @@ -335,8 +333,8 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p
if (use_EOS) then
!$OMP parallel do default(shared) private(rho_in_situ)
do j=Jsq,Jeq+1
call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p(:,j,1), rho_in_situ, Isq, Ieq-Isq+2, &
tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p(:,j,1), rho_in_situ, G%HI, &
tv%eqn_of_state, US, halo=1)

do i=Isq,Ieq+1
dM(i,j) = (CS%GFS_scale - 1.0) * (p(i,j,1)*(1.0/rho_in_situ(i) - alpha_ref) + za(i,j))
Expand Down Expand Up @@ -601,11 +599,11 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_at
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1
if (use_p_atm) then
call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p_atm(:,j), rho_in_situ, &
Isq, Ieq-Isq+2, tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p_atm(:,j), rho_in_situ, G%HI, &
tv%eqn_of_state, US, halo=1)
else
call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p0, rho_in_situ, &
Isq, Ieq-Isq+2, tv%eqn_of_state, scale=US%kg_m3_to_R)
call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p0, rho_in_situ, G%HI, &
tv%eqn_of_state, US, halo=1)
endif
do i=Isq,Ieq+1
dM(i,j) = (CS%GFS_scale - 1.0) * (G_Rho0 * rho_in_situ(i)) * e(i,j,1)
Expand Down Expand Up @@ -668,20 +666,19 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_at
if ( use_ALE ) then
if ( CS%Recon_Scheme == 1 ) then
call int_density_dz_generic_plm( T_t(:,:,k), T_b(:,:,k), S_t(:,:,k), S_b(:,:,k),&
e(:,:,K), e(:,:,K+1), rho_ref, CS%Rho0, GV%g_Earth, dz_neglect, G%bathyT, G%HI, G%HI, &
tv%eqn_of_state, dpa, intz_dpa, intx_dpa, inty_dpa, useMassWghtInterp=CS%useMassWghtInterp, &
rho_scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
e(:,:,K), e(:,:,K+1), rho_ref, CS%Rho0, GV%g_Earth, dz_neglect, G%bathyT, &
G%HI, G%HI, tv%eqn_of_state, dpa, intz_dpa, intx_dpa, inty_dpa, &
useMassWghtInterp=CS%useMassWghtInterp, US=US)
elseif ( CS%Recon_Scheme == 2 ) then
call int_density_dz_generic_ppm( tv%T(:,:,k), T_t(:,:,k), T_b(:,:,k), &
tv%S(:,:,k), S_t(:,:,k), S_b(:,:,k), e(:,:,K), e(:,:,K+1), rho_ref, CS%Rho0, &
GV%g_Earth, G%HI, G%HI, tv%eqn_of_state, dpa, intz_dpa, intx_dpa, inty_dpa, &
rho_scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
tv%S(:,:,k), S_t(:,:,k), S_b(:,:,k), e(:,:,K), e(:,:,K+1), &
rho_ref, CS%Rho0, GV%g_Earth, G%HI, G%HI, tv%eqn_of_state, dpa, &
intz_dpa, intx_dpa, inty_dpa, US=US)
endif
else
call int_density_dz(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), e(:,:,K), e(:,:,K+1), &
rho_ref, CS%Rho0, GV%g_Earth, G%HI, G%HI, tv%eqn_of_state, &
dpa, intz_dpa, intx_dpa, inty_dpa, G%bathyT, dz_neglect, CS%useMassWghtInterp, &
rho_scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
rho_ref, CS%Rho0, GV%g_Earth, G%HI, G%HI, tv%eqn_of_state, dpa, &
intz_dpa, intx_dpa, inty_dpa, G%bathyT, dz_neglect, CS%useMassWghtInterp, US=US)
endif
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
Expand Down
14 changes: 5 additions & 9 deletions src/core/MOM_PressureForce_blocked_AFV.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -248,8 +248,7 @@ subroutine PressureForce_blk_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm,
p(:,:,K+1), alpha_ref, G%HI, tv%eqn_of_state, &
dza(:,:,k), intp_dza(:,:,k), intx_dza(:,:,k), &
inty_dza(:,:,k), bathyP=p(:,:,nz+1), dP_tiny=dp_neglect, &
useMassWghtInterp=CS%useMassWghtInterp, &
SV_scale=US%R_to_kg_m3, pres_scale=US%RL2_T2_to_Pa)
useMassWghtInterp=CS%useMassWghtInterp, US=US)
else
alpha_anom = 1.0 / GV%Rlay(k) - alpha_ref
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
Expand Down Expand Up @@ -669,21 +668,18 @@ subroutine PressureForce_blk_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp,
call int_density_dz_generic_plm( T_t(:,:,k), T_b(:,:,k), S_t(:,:,k), S_b(:,:,k), &
e(:,:,K), e(:,:,K+1), rho_ref, CS%Rho0, GV%g_Earth, dz_neglect, G%bathyT, G%HI, &
G%Block(n), tv%eqn_of_state, dpa_bk, intz_dpa_bk, intx_dpa_bk, inty_dpa_bk, &
useMassWghtInterp=CS%useMassWghtInterp, &
rho_scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
useMassWghtInterp=CS%useMassWghtInterp, US=US)
elseif ( CS%Recon_Scheme == 2 ) then
call int_density_dz_generic_ppm( tv%T(:,:,k), T_t(:,:,k), T_b(:,:,k), &
tv%S(:,:,k), S_t(:,:,k), S_b(:,:,k), e(:,:,K), e(:,:,K+1), rho_ref, CS%Rho0, &
GV%g_Earth, G%HI, G%Block(n), tv%eqn_of_state, dpa_bk, intz_dpa_bk, &
intx_dpa_bk, inty_dpa_bk, &
rho_scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
GV%g_Earth, G%HI, G%Block(n), tv%eqn_of_state, dpa_bk, intz_dpa_bk, &
intx_dpa_bk, inty_dpa_bk, US=US)
endif
else
call int_density_dz(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), e(:,:,K), e(:,:,K+1), &
rho_ref, CS%Rho0, GV%g_Earth, G%HI, G%Block(n), tv%eqn_of_state, &
dpa_bk, intz_dpa_bk, intx_dpa_bk, inty_dpa_bk, &
G%bathyT, dz_neglect, CS%useMassWghtInterp, &
rho_scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
G%bathyT, dz_neglect, CS%useMassWghtInterp, US=US)
endif
do jb=Jsq_bk,Jeq_bk+1 ; do ib=Isq_bk,Ieq_bk+1
intz_dpa_bk(ib,jb) = intz_dpa_bk(ib,jb)*GV%Z_to_H
Expand Down
6 changes: 2 additions & 4 deletions src/core/MOM_interface_heights.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -106,8 +106,7 @@ subroutine find_eta_3d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m)
!$OMP do
do k=1,nz
call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,K), p(:,:,K+1), &
0.0, G%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=halo, &
SV_scale=US%R_to_kg_m3, pres_scale=US%RL2_T2_to_Pa)
0.0, G%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=halo, US=US)
enddo
!$OMP do
do j=jsv,jev
Expand Down Expand Up @@ -209,8 +208,7 @@ subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m)
!$OMP do
do k = 1, nz
call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,k), p(:,:,k+1), 0.0, &
G%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=halo, &
SV_scale=US%R_to_kg_m3, pres_scale=US%RL2_T2_to_Pa)
G%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=halo, US=US)
enddo
!$OMP do
do j=js,je ; do k=1,nz ; do i=is,ie
Expand Down
3 changes: 1 addition & 2 deletions src/diagnostics/MOM_diagnostics.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -842,8 +842,7 @@ subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS)
z_bot(i,j) = z_top(i,j) - GV%H_to_Z*h(i,j,k)
enddo ; enddo
call int_density_dz(tv%T(:,:,k), tv%S(:,:,k), z_top, z_bot, 0.0, GV%Rho0, GV%g_Earth, &
G%HI, G%HI, tv%eqn_of_state, dpress, rho_scale=US%kg_m3_to_R, &
pres_scale=US%RL2_T2_to_Pa)
G%HI, G%HI, tv%eqn_of_state, dpress, US=US)
do j=js,je ; do i=is,ie
mass(i,j) = mass(i,j) + dpress(i,j) * IG_Earth
enddo ; enddo
Expand Down
Loading

0 comments on commit ce0905e

Please sign in to comment.