Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

+*Ignore SURFACE_ANSWER_DATE when non-Boussinesq #422

Merged
merged 2 commits into from
Aug 10, 2023
Merged

+*Ignore SURFACE_ANSWER_DATE when non-Boussinesq #422

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
02135ae
+*Ignore SURFACE_ANSWER_DATE when non-Boussinesq
Hallberg-NOAA Jun 30, 2023
8485b1a
Merge branch 'dev/gfdl' into nonBous_without_EOS
marshallward Aug 10, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 3 additions & 3 deletions src/core/MOM.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2103,7 +2103,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
call get_param(param_file, "MOM", "CALC_RHO_FOR_SEA_LEVEL", CS%calc_rho_for_sea_lev, &
"If true, the in-situ density is used to calculate the "//&
"effective sea level that is returned to the coupler. If false, "//&
"the Boussinesq parameter RHO_0 is used.", default=.false.)
"the Boussinesq parameter RHO_0 is used.", default=non_Bous)
call get_param(param_file, "MOM", "ENABLE_THERMODYNAMICS", use_temperature, &
"If true, Temperature and salinity are used as state "//&
"variables.", default=.true.)
Expand Down Expand Up @@ -2892,8 +2892,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
endif
endif

! Allocate any derived equation of state fields.
if (use_temperature .and. .not.(GV%Boussinesq .or. GV%semi_Boussinesq)) then
! Allocate any derived densities or other equation of state derived fields.
if (.not.(GV%Boussinesq .or. GV%semi_Boussinesq)) then
allocate(CS%tv%SpV_avg(isd:ied,jsd:jed,nz), source=0.0)
CS%tv%valid_SpV_halo = -1 ! This array does not yet have any valid data.
endif
Expand Down
44 changes: 33 additions & 11 deletions src/core/MOM_interface_heights.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -279,6 +279,8 @@ subroutine calc_derived_thermo(tv, h, G, GV, US, halo, debug)
! Local variables
real, dimension(SZI_(G),SZJ_(G)) :: p_t ! Hydrostatic pressure atop a layer [R L2 T-2 ~> Pa]
real, dimension(SZI_(G),SZJ_(G)) :: dp ! Pressure change across a layer [R L2 T-2 ~> Pa]
real, dimension(SZK_(GV)) :: SpV_lay ! The specific volume of each layer when no equation of
! state is used [R-1 ~> m3 kg-1]
logical :: do_debug ! If true, write checksums for debugging.
integer :: i, j, k, is, ie, js, je, halos, nz

Expand Down Expand Up @@ -310,6 +312,12 @@ subroutine calc_derived_thermo(tv, h, G, GV, US, halo, debug)
call hchksum(tv%T, "derived_thermo T", G%HI, haloshift=halos, scale=US%C_to_degC)
call hchksum(tv%S, "derived_thermo S", G%HI, haloshift=halos, scale=US%S_to_ppt)
endif
elseif (allocated(tv%Spv_avg)) then
do k=1,nz ; SpV_lay(k) = 1.0 / GV%Rlay(k) ; enddo
do k=1,nz ; do j=js,je ; do i=is,ie
tv%SpV_avg(i,j,k) = SpV_lay(k)
enddo ; enddo ; enddo
tv%valid_SpV_halo = halos
endif

end subroutine calc_derived_thermo
Expand Down Expand Up @@ -481,9 +489,7 @@ subroutine dz_to_thickness_tv(dz, tv, h, G, GV, US, halo_size)
endif
else
do k=1,nz ; do j=js,je ; do i=is,ie
h(i,j,k) = (GV%Z_to_H*dz(i,j,k)) * (GV%Rlay(k) / GV%Rho0)
! Consider revising this to the mathematically equivalent expression:
! h(i,j,k) = (GV%RZ_to_H * GV%Rlay(k)) * dz(i,j,k)
h(i,j,k) = (GV%RZ_to_H * GV%Rlay(k)) * dz(i,j,k)
enddo ; enddo ; enddo
endif
endif
Expand Down Expand Up @@ -551,10 +557,16 @@ subroutine dz_to_thickness_EOS(dz, Temp, Saln, EoS, h, G, GV, US, halo_size, p_s
do i=is,ie ; p_top(i,j) = p_bot(i,j) ; enddo
call calculate_density(Temp(:,j,k), Saln(:,j,k), p_top(:,j), rho, &
EoS, EOSdom)
do i=is,ie
! This could be simplified, but it would change answers at roundoff.
p_bot(i,j) = p_top(i,j) + (GV%g_Earth*GV%H_to_Z) * ((GV%Z_to_H*dz(i,j,k)) * rho(i))
enddo
! The following two expressions are mathematically equivalent.
if (GV%semi_Boussinesq) then
do i=is,ie
p_bot(i,j) = p_top(i,j) + (GV%g_Earth*GV%H_to_Z) * ((GV%Z_to_H*dz(i,j,k)) * rho(i))
enddo
else
do i=is,ie
p_bot(i,j) = p_top(i,j) + rho(i) * (GV%g_Earth * dz(i,j,k))
enddo
endif
enddo

do itt=1,max_itt
Expand All @@ -565,9 +577,15 @@ subroutine dz_to_thickness_EOS(dz, Temp, Saln, EoS, h, G, GV, US, halo_size, p_s
EoS, EOSdom)
! Use Newton's method to correct the bottom value.
! The hydrostatic equation is sufficiently linear that no bounds-checking is needed.
do i=is,ie
p_bot(i,j) = p_bot(i,j) + rho(i) * ((GV%g_Earth*GV%H_to_Z)*(GV%Z_to_H*dz(i,j,k)) - dz_geo(i,j))
enddo
if (GV%semi_Boussinesq) then
do i=is,ie
p_bot(i,j) = p_bot(i,j) + rho(i) * ((GV%g_Earth*GV%H_to_Z)*(GV%Z_to_H*dz(i,j,k)) - dz_geo(i,j))
enddo
else
do i=is,ie
p_bot(i,j) = p_bot(i,j) + rho(i) * (GV%g_Earth*dz(i,j,k) - dz_geo(i,j))
enddo
endif
enddo ; endif
enddo

Expand Down Expand Up @@ -608,14 +626,18 @@ subroutine dz_to_thickness_simple(dz, h, G, GV, US, halo_size, layer_mode)
layered = .false. ; if (present(layer_mode)) layered = layer_mode
is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo ; nz = GV%ke

if (GV%Boussinesq .or. (.not.layered)) then
if (GV%Boussinesq) then
do k=1,nz ; do j=js,je ; do i=is,ie
h(i,j,k) = GV%Z_to_H * dz(i,j,k)
enddo ; enddo ; enddo
elseif (layered) then
do k=1,nz ; do j=js,je ; do i=is,ie
h(i,j,k) = (GV%RZ_to_H * GV%Rlay(k)) * dz(i,j,k)
enddo ; enddo ; enddo
else
do k=1,nz ; do j=js,je ; do i=is,ie
h(i,j,k) = (US%Z_to_m * GV%m_to_H) * dz(i,j,k)
enddo ; enddo ; enddo
endif

end subroutine dz_to_thickness_simple
Expand Down