Skip to content

Commit

Permalink
*+Revise non-Boussinesq gprime expressions
Browse files Browse the repository at this point in the history 
  Revised the calculation of gprime and the coordinate densities (GV%Rlay) in
fully non-Boussinesq mode to use the arithmetic mean of adjacent coordinate
densities in the denominator of the expression for g_prime in place of RHO_0.
Also use LIGHTEST_DENSITY in place of RHO_0 to specify the top-level coordinate
density in certain coordinate modes.  Also made corresponding changes to the
fully non-Boussinesq APE calculation when CALCULATE_APE is true, and eliminated
an incorrect calculation of the layer volumes in non-Boussinesq mode using the
Boussinesq reference density that was never actually being used when
CALCULATE_APE is false.

  This commit will change answers in some fully non-Boussinesq calculations, and
an existing runtime parameter is used and logged in some new cases, changing
the MOM_parameter_doc file in those cases.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Jul 3, 2023
1 parent d44c228 commit 9120755
Show file tree
Hide file tree
Showing 2 changed files with 101 additions and 31 deletions.
38 changes: 23 additions & 15 deletions src/diagnostics/MOM_sum_output.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -510,24 +510,18 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci
do k=1,nz ; vol_lay(k) = (US%m_to_L**2*GV%H_to_Z/GV%H_to_kg_m2)*mass_lay(k) ; enddo
else
tmp1(:,:,:) = 0.0
if (CS%do_APE_calc) then
do k=1,nz ; do j=js,je ; do i=is,ie
tmp1(i,j,k) = HL2_to_kg * h(i,j,k) * areaTm(i,j)
enddo ; enddo ; enddo
mass_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, sums=mass_lay, EFP_sum=mass_EFP)
do k=1,nz ; do j=js,je ; do i=is,ie
tmp1(i,j,k) = HL2_to_kg * h(i,j,k) * areaTm(i,j)
enddo ; enddo ; enddo
mass_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, sums=mass_lay, EFP_sum=mass_EFP)

if (CS%do_APE_calc) then
call find_eta(h, tv, G, GV, US, eta, dZref=G%Z_ref)
do k=1,nz ; do j=js,je ; do i=is,ie
tmp1(i,j,k) = US%Z_to_m*US%L_to_m**2*(eta(i,j,K)-eta(i,j,K+1)) * areaTm(i,j)
enddo ; enddo ; enddo
vol_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, sums=vol_lay)
do k=1,nz ; vol_lay(k) = US%m_to_Z*US%m_to_L**2 * vol_lay(k) ; enddo
else
do k=1,nz ; do j=js,je ; do i=is,ie
tmp1(i,j,k) = HL2_to_kg * h(i,j,k) * areaTm(i,j)
enddo ; enddo ; enddo
mass_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, sums=mass_lay, EFP_sum=mass_EFP)
do k=1,nz ; vol_lay(k) = US%m_to_Z*US%m_to_L**2*US%kg_m3_to_R * (mass_lay(k) / GV%Rho0) ; enddo
endif
endif ! Boussinesq

Expand Down Expand Up @@ -643,7 +637,7 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci
if (GV%Boussinesq) then
do j=js,je ; do i=is,ie
hbelow = 0.0
do k=nz,1,-1
do K=nz,1,-1
hbelow = hbelow + h(i,j,k) * GV%H_to_Z
hint = Z_0APE(K) + (hbelow - (G%bathyT(i,j) + G%Z_ref))
hbot = Z_0APE(K) - (G%bathyT(i,j) + G%Z_ref)
Expand All @@ -652,14 +646,28 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci
(hint * hint - hbot * hbot)
enddo
enddo ; enddo
else
elseif (GV%semi_Boussinesq) then
do j=js,je ; do i=is,ie
do k=nz,1,-1
do K=nz,1,-1
hint = Z_0APE(K) + eta(i,j,K) ! eta and H_0 have opposite signs.
hbot = max(Z_0APE(K) - (G%bathyT(i,j) + G%Z_ref), 0.0)
PE_pt(i,j,K) = (0.5 * PE_scale_factor * areaTm(i,j) * (GV%Rho0*GV%g_prime(K))) * &
(hint * hint - hbot * hbot)
(hint * hint - hbot * hbot)
enddo
enddo ; enddo
else
do j=js,je ; do i=is,ie
do K=nz,2,-1
hint = Z_0APE(K) + eta(i,j,K) ! eta and H_0 have opposite signs.
hbot = max(Z_0APE(K) - (G%bathyT(i,j) + G%Z_ref), 0.0)
PE_pt(i,j,K) = (0.25 * PE_scale_factor * areaTm(i,j) * &
((GV%Rlay(k)+GV%Rlay(k-1))*GV%g_prime(K))) * &
(hint * hint - hbot * hbot)
enddo
hint = Z_0APE(1) + eta(i,j,1) ! eta and H_0 have opposite signs.
hbot = max(Z_0APE(1) - (G%bathyT(i,j) + G%Z_ref), 0.0)
PE_pt(i,j,1) = (0.5 * PE_scale_factor * areaTm(i,j) * (GV%Rlay(1)*GV%g_prime(1))) * &
(hint * hint - hbot * hbot)
enddo ; enddo
endif

Expand Down
94 changes: 78 additions & 16 deletions src/initialization/MOM_coord_initialization.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -126,6 +126,7 @@ subroutine set_coord_from_gprime(Rlay, g_prime, GV, US, param_file)
! Local variables
real :: g_int ! Reduced gravities across the internal interfaces [L2 Z-1 T-2 ~> m s-2].
real :: g_fs ! Reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2].
real :: Rlay_Ref ! The target density of the surface layer [R ~> kg m-3].
character(len=40) :: mdl = "set_coord_from_gprime" ! This subroutine's name.
integer :: k, nz
nz = GV%ke
Expand All @@ -138,11 +139,20 @@ subroutine set_coord_from_gprime(Rlay, g_prime, GV, US, param_file)
call get_param(param_file, mdl, "GINT", g_int, &
"The reduced gravity across internal interfaces.", &
units="m s-2", fail_if_missing=.true., scale=US%m_s_to_L_T**2*US%Z_to_m)
call get_param(param_file, mdl, "LIGHTEST_DENSITY", Rlay_Ref, &
"The reference potential density used for layer 1.", &
units="kg m-3", default=US%R_to_kg_m3*GV%Rho0, scale=US%kg_m3_to_R)

g_prime(1) = g_fs
do k=2,nz ; g_prime(k) = g_int ; enddo
Rlay(1) = GV%Rho0
do k=2,nz ; Rlay(k) = Rlay(k-1) + g_prime(k)*(GV%Rho0/GV%g_Earth) ; enddo
Rlay(1) = Rlay_Ref
if (GV%Boussinesq .or. GV%semi_Boussinesq) then
do k=2,nz ; Rlay(k) = Rlay(k-1) + g_prime(k)*(GV%Rho0/GV%g_Earth) ; enddo
else
do k=2,nz
Rlay(k) = Rlay(k-1) * ((GV%g_Earth + 0.5*g_prime(k)) / (GV%g_Earth - 0.5*g_prime(k)))
enddo
endif

call callTree_leave(trim(mdl)//'()')

Expand Down Expand Up @@ -184,9 +194,15 @@ subroutine set_coord_from_layer_density(Rlay, g_prime, GV, US, param_file)
enddo
! These statements set the interface reduced gravities. !
g_prime(1) = g_fs
do k=2,nz
g_prime(k) = (GV%g_Earth/(GV%Rho0)) * (Rlay(k) - Rlay(k-1))
enddo
if (GV%Boussinesq .or. GV%semi_Boussinesq) then
do k=2,nz
g_prime(k) = (GV%g_Earth/GV%Rho0) * (Rlay(k) - Rlay(k-1))
enddo
else
do k=2,nz
g_prime(k) = 2.0*GV%g_Earth * (Rlay(k) - Rlay(k-1)) / (Rlay(k) + Rlay(k-1))
enddo
endif

call callTree_leave(trim(mdl)//'()')
end subroutine set_coord_from_layer_density
Expand Down Expand Up @@ -237,7 +253,13 @@ subroutine set_coord_from_TS_ref(Rlay, g_prime, GV, US, param_file, eqn_of_state
call calculate_density(T_ref, S_ref, P_ref, Rlay(1), eqn_of_state)

! These statements set the layer densities. !
do k=2,nz ; Rlay(k) = Rlay(k-1) + g_prime(k)*(GV%Rho0/GV%g_Earth) ; enddo
if (GV%Boussinesq .or. GV%semi_Boussinesq) then
do k=2,nz ; Rlay(k) = Rlay(k-1) + g_prime(k)*(GV%Rho0/GV%g_Earth) ; enddo
else
do k=2,nz
Rlay(k) = Rlay(k-1) * ((GV%g_Earth + 0.5*g_prime(k)) / (GV%g_Earth - 0.5*g_prime(k)))
enddo
endif

call callTree_leave(trim(mdl)//'()')
end subroutine set_coord_from_TS_ref
Expand Down Expand Up @@ -294,7 +316,15 @@ subroutine set_coord_from_TS_profile(Rlay, g_prime, GV, US, param_file, eqn_of_s
g_prime(1) = g_fs
do k=1,nz ; Pref(k) = P_Ref ; enddo
call calculate_density(T0, S0, Pref, Rlay, eqn_of_state, (/1,nz/) )
do k=2,nz; g_prime(k) = (GV%g_Earth/(GV%Rho0)) * (Rlay(k) - Rlay(k-1)) ; enddo
if (GV%Boussinesq .or. GV%semi_Boussinesq) then
do k=2,nz
g_prime(k) = (GV%g_Earth/GV%Rho0) * (Rlay(k) - Rlay(k-1))
enddo
else
do k=2,nz
g_prime(k) = 2.0*GV%g_Earth * (Rlay(k) - Rlay(k-1)) / (Rlay(k) + Rlay(k-1))
enddo
endif

call callTree_leave(trim(mdl)//'()')
end subroutine set_coord_from_TS_profile
Expand Down Expand Up @@ -387,7 +417,15 @@ subroutine set_coord_from_TS_range(Rlay, g_prime, GV, US, param_file, eqn_of_sta
do k=k_light-1,1,-1
Rlay(k) = 2.0*Rlay(k+1) - Rlay(k+2)
enddo
do k=2,nz ; g_prime(k) = (GV%g_Earth/GV%Rho0) * (Rlay(k) - Rlay(k-1)) ; enddo
if (GV%Boussinesq .or. GV%semi_Boussinesq) then
do k=2,nz
g_prime(k) = (GV%g_Earth/GV%Rho0) * (Rlay(k) - Rlay(k-1))
enddo
else
do k=2,nz
g_prime(k) = 2.0*GV%g_Earth * (Rlay(k) - Rlay(k-1)) / (Rlay(k) + Rlay(k-1))
enddo
endif

call callTree_leave(trim(mdl)//'()')
end subroutine set_coord_from_TS_range
Expand Down Expand Up @@ -429,7 +467,15 @@ subroutine set_coord_from_file(Rlay, g_prime, GV, US, param_file)

call MOM_read_data(filename, coord_var, Rlay, scale=US%kg_m3_to_R)
g_prime(1) = g_fs
do k=2,nz ; g_prime(k) = (GV%g_Earth/GV%Rho0) * (Rlay(k) - Rlay(k-1)) ; enddo
if (GV%Boussinesq .or. GV%semi_Boussinesq) then
do k=2,nz
g_prime(k) = (GV%g_Earth/GV%Rho0) * (Rlay(k) - Rlay(k-1))
enddo
else
do k=2,nz
g_prime(k) = 2.0*GV%g_Earth * (Rlay(k) - Rlay(k-1)) / (Rlay(k) + Rlay(k-1))
enddo
endif
do k=1,nz ; if (g_prime(k) <= 0.0) then
call MOM_error(FATAL, "MOM_initialization set_coord_from_file: "//&
"Zero or negative g_primes read from variable "//"Layer"//" in file "//&
Expand Down Expand Up @@ -479,9 +525,15 @@ subroutine set_coord_linear(Rlay, g_prime, GV, US, param_file)
enddo
! These statements set the interface reduced gravities.
g_prime(1) = g_fs
do k=2,nz
g_prime(k) = (GV%g_Earth/(GV%Rho0)) * (Rlay(k) - Rlay(k-1))
enddo
if (GV%Boussinesq .or. GV%semi_Boussinesq) then
do k=2,nz
g_prime(k) = (GV%g_Earth/GV%Rho0) * (Rlay(k) - Rlay(k-1))
enddo
else
do k=2,nz
g_prime(k) = 2.0*GV%g_Earth * (Rlay(k) - Rlay(k-1)) / (Rlay(k) + Rlay(k-1))
enddo
endif

call callTree_leave(trim(mdl)//'()')
end subroutine set_coord_linear
Expand All @@ -498,6 +550,7 @@ subroutine set_coord_to_none(Rlay, g_prime, GV, US, param_file)
type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
! Local variables
real :: g_fs ! Reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2].
real :: Rlay_Ref ! The target density of the surface layer [R ~> kg m-3].
character(len=40) :: mdl = "set_coord_to_none" ! This subroutine's name.
integer :: k, nz
nz = GV%ke
Expand All @@ -507,11 +560,20 @@ subroutine set_coord_to_none(Rlay, g_prime, GV, US, param_file)
call get_param(param_file, mdl, "GFS" , g_fs, &
"The reduced gravity at the free surface.", units="m s-2", &
default=GV%g_Earth*US%L_T_to_m_s**2*US%m_to_Z, scale=US%m_s_to_L_T**2*US%Z_to_m)
call get_param(param_file, mdl, "LIGHTEST_DENSITY", Rlay_Ref, &
"The reference potential density used for layer 1.", &
units="kg m-3", default=US%R_to_kg_m3*GV%Rho0, scale=US%kg_m3_to_R)

g_prime(1) = g_fs
do k=2,nz ; g_prime(k) = 0. ; enddo
Rlay(1) = GV%Rho0
do k=2,nz ; Rlay(k) = Rlay(k-1) + g_prime(k)*(GV%Rho0/GV%g_Earth) ; enddo
Rlay(1) = Rlay_Ref
if (GV%Boussinesq .or. GV%semi_Boussinesq) then
do k=2,nz ; Rlay(k) = Rlay(k-1) + g_prime(k)*(GV%Rho0/GV%g_Earth) ; enddo
else
do k=2,nz
Rlay(k) = Rlay(k-1) * ((GV%g_Earth + 0.5*g_prime(k)) / (GV%g_Earth - 0.5*g_prime(k)))
enddo
endif

call callTree_leave(trim(mdl)//'()')

Expand All @@ -522,8 +584,8 @@ end subroutine set_coord_to_none
subroutine write_vertgrid_file(GV, US, param_file, directory)
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
character(len=*), intent(in) :: directory !< The directory into which to place the file.
type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
character(len=*), intent(in) :: directory !< The directory into which to place the file.
! Local variables
character(len=240) :: filepath
type(vardesc) :: vars(2)
Expand Down

0 comments on commit 9120755

Please sign in to comment.