Skip to content

Commit

Permalink
+Pass depth_tot to initialization routines
Browse files Browse the repository at this point in the history 
  Added depth_tot arguments to explicitly pass the total depth to various
initialization routines for the thicknesses, sponges, or temperatures.  These
routines had previously used G%bathT for this role, but this change prepares to
allow the reference depth for the bathymetry to change without requiring these
routines to be changed to accommodate it.  All answers are bitwise identical,
but there are new arguments to about two dozen routines.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Aug 17, 2021
1 parent 5017bf6 commit 3cda6fd
Show file tree
Hide file tree
Showing 19 changed files with 212 additions and 144 deletions.
109 changes: 64 additions & 45 deletions src/initialization/MOM_state_initialization.F90
View file

Large diffs are not rendered by default.

6 changes: 4 additions & 2 deletions src/user/BFB_initialization.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,13 +76,15 @@ end subroutine BFB_set_coord

!> This subroutine sets up the sponges for the southern bouundary of the domain. Maximum damping occurs
!! within 2 degrees lat of the boundary. The damping linearly decreases northward over the next 2 degrees.
subroutine BFB_initialize_sponges_southonly(G, GV, US, use_temperature, tv, param_file, CSp, h)
subroutine BFB_initialize_sponges_southonly(G, GV, US, use_temperature, tv, depth_tot, param_file, CSp, h)
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
logical, intent(in) :: use_temperature !< If true, temperature and salinity are used as
!! state variables.
type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables
real, dimension(SZI_(G),SZJ_(G)), &
intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
type(sponge_CS), pointer :: CSp !< A pointer to the sponge control structure
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
Expand Down Expand Up @@ -129,7 +131,7 @@ subroutine BFB_initialize_sponges_southonly(G, GV, US, use_temperature, tv, para
max_damping = 1.0 / (86400.0*US%s_to_T)

do j=js,je ; do i=is,ie
if (G%bathyT(i,j) <= min_depth) then ; Idamp(i,j) = 0.0
if (depth_tot(i,j) <= min_depth) then ; Idamp(i,j) = 0.0
elseif (G%geoLatT(i,j) < slat+2.0) then ; Idamp(i,j) = max_damping
elseif (G%geoLatT(i,j) < slat+4.0) then
Idamp(i,j) = max_damping * (slat+4.0-G%geoLatT(i,j))/2.0
Expand Down
24 changes: 14 additions & 10 deletions src/user/DOME2d_initialization.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,12 +90,14 @@ subroutine DOME2d_initialize_topography( D, G, param_file, max_depth )
end subroutine DOME2d_initialize_topography

!> Initialize thicknesses according to coordinate mode
subroutine DOME2d_initialize_thickness ( h, G, GV, US, param_file, just_read_params )
subroutine DOME2d_initialize_thickness ( h, depth_tot, G, GV, US, param_file, just_read_params )
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
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2].
real, dimension(SZI_(G),SZJ_(G)), &
intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
!! to parse for model parameter values.
logical, optional, intent(in) :: just_read_params !< If present and true, this call will
Expand Down Expand Up @@ -150,7 +152,7 @@ subroutine DOME2d_initialize_thickness ( h, G, GV, US, param_file, just_read_par
case ( REGRIDDING_LAYER, REGRIDDING_RHO )

do j=js,je ; do i=is,ie
eta1D(nz+1) = -G%bathyT(i,j)
eta1D(nz+1) = -depth_tot(i,j)
do k=nz,1,-1
eta1D(k) = e0(k)
if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_Z)) then
Expand All @@ -172,7 +174,7 @@ subroutine DOME2d_initialize_thickness ( h, G, GV, US, param_file, just_read_par
! case ( IC_RHO_C )
!
! do j=js,je ; do i=is,ie
! eta1D(nz+1) = -G%bathyT(i,j)
! eta1D(nz+1) = -depth_tot(i,j)
! do k=nz,1,-1
! eta1D(k) = e0(k)
! if (eta1D(k) < (eta1D(k+1) + min_thickness)) then
Expand All @@ -194,7 +196,7 @@ subroutine DOME2d_initialize_thickness ( h, G, GV, US, param_file, just_read_par
case ( REGRIDDING_ZSTAR )

do j=js,je ; do i=is,ie
eta1D(nz+1) = -G%bathyT(i,j)
eta1D(nz+1) = -depth_tot(i,j)
do k=nz,1,-1
eta1D(k) = e0(k)
if (eta1D(k) < (eta1D(k+1) + min_thickness)) then
Expand All @@ -208,7 +210,7 @@ subroutine DOME2d_initialize_thickness ( h, G, GV, US, param_file, just_read_par

case ( REGRIDDING_SIGMA )
do j=js,je ; do i=is,ie
h(i,j,:) = GV%Z_to_H*G%bathyT(i,j) / nz
h(i,j,:) = GV%Z_to_H*depth_tot(i,j) / nz
enddo ; enddo

case default
Expand Down Expand Up @@ -353,11 +355,13 @@ subroutine DOME2d_initialize_temperature_salinity ( T, S, h, G, GV, param_file,
end subroutine DOME2d_initialize_temperature_salinity

!> Set up sponges in 2d DOME configuration
subroutine DOME2d_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CSp, ACSp)
subroutine DOME2d_initialize_sponges(G, GV, US, tv, depth_tot, param_file, use_ALE, CSp, ACSp)
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(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure
real, dimension(SZI_(G),SZJ_(G)), &
intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
type(param_file_type), intent(in) :: param_file !< Parameter file structure
logical, intent(in) :: use_ALE !< If true, indicates model is in ALE mode
type(sponge_CS), pointer :: CSp !< Layer-mode sponge structure
Expand Down Expand Up @@ -453,7 +457,7 @@ subroutine DOME2d_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CSp, AC
enddo
e0(nz+1) = -G%max_depth
do j=js,je ; do i=is,ie
eta1D(nz+1) = -G%bathyT(i,j)
eta1D(nz+1) = -depth_tot(i,j)
do k=nz,1,-1
eta1D(k) = e0(k)
if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_Z)) then
Expand All @@ -470,7 +474,7 @@ subroutine DOME2d_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CSp, AC
! Construct temperature and salinity on the arbitrary grid
T(:,:,:) = 0.0 ; S(:,:,:) = 0.0
do j=js,je ; do i=is,ie
z = -G%bathyT(i,j)
z = -depth_tot(i,j)
do k = nz,1,-1
z = z + 0.5 * GV%H_to_Z * h(i,j,k) ! Position of the center of layer k
S(i,j,k) = 34.0 - 1.0 * (z / (G%max_depth))
Expand All @@ -491,7 +495,7 @@ subroutine DOME2d_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CSp, AC

! Construct interface heights to restore toward
do j=js,je ; do i=is,ie
eta1D(nz+1) = -G%bathyT(i,j)
eta1D(nz+1) = -depth_tot(i,j)
do k=nz,1,-1
eta1D(K) = -G%max_depth * real(k-1) / real(nz)
if (eta1D(K) < (eta1D(K+1) + GV%Angstrom_Z)) then
Expand All @@ -508,7 +512,7 @@ subroutine DOME2d_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CSp, AC
d_eta(nz) = dome2d_depth_bay * G%max_depth - (nz-1) * GV%Angstrom_Z
endif

eta(i,j,nz+1) = -G%bathyT(i,j)
eta(i,j,nz+1) = -depth_tot(i,j)
do K=nz,1,-1
eta(i,j,K) = eta(i,j,K+1) + d_eta(k)
enddo
Expand Down
38 changes: 21 additions & 17 deletions src/user/DOME_initialization.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,11 +87,13 @@ end subroutine DOME_initialize_topography

! -----------------------------------------------------------------------------
!> This subroutine initializes layer thicknesses for the DOME experiment
subroutine DOME_initialize_thickness(h, G, GV, param_file, just_read_params)
subroutine DOME_initialize_thickness(h, depth_tot, G, GV, param_file, just_read_params)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2].
real, dimension(SZI_(G),SZJ_(G)), &
intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
!! to parse for model parameter values.
logical, optional, intent(in) :: just_read_params !< If present and true, this call will
Expand Down Expand Up @@ -124,7 +126,7 @@ subroutine DOME_initialize_thickness(h, G, GV, param_file, just_read_params)
! Angstrom thick, and 2. the interfaces are where they should be !
! based on the resting depths and interface height perturbations, !
! as long at this doesn't interfere with 1. !
eta1D(nz+1) = -G%bathyT(i,j)
eta1D(nz+1) = -depth_tot(i,j)
do k=nz,1,-1
eta1D(K) = e0(K)
if (eta1D(K) < (eta1D(K+1) + GV%Angstrom_Z)) then
Expand All @@ -145,17 +147,19 @@ end subroutine DOME_initialize_thickness
!! number of tracers should be restored within each sponge. The !
!! interface height is always subject to damping, and must always be !
!! the first registered field. !
subroutine DOME_initialize_sponges(G, GV, US, tv, PF, CSp)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
subroutine DOME_initialize_sponges(G, GV, US, tv, depth_tot, PF, CSp)
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
type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any available
!! thermodynamic fields, including potential temperature and
!! salinity or mixed layer density. Absent fields have NULL ptrs.
type(param_file_type), intent(in) :: PF !< A structure indicating the open file to
!! parse for model parameter values.
type(sponge_CS), pointer :: CSp !< A pointer that is set to point to the control
!! structure for this module.
type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any available
!! thermodynamic fields, including potential temperature and
!! salinity or mixed layer density. Absent fields have NULL ptrs.
real, dimension(SZI_(G),SZJ_(G)), &
intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
type(param_file_type), intent(in) :: PF !< A structure indicating the open file to
!! parse for model parameter values.
type(sponge_CS), pointer :: CSp !< A pointer that is set to point to the control
!! structure for this module.

real :: eta(SZI_(G),SZJ_(G),SZK_(GV)+1) ! A temporary array for interface heights [Z ~> m].
real :: temp(SZI_(G),SZJ_(G),SZK_(GV)) ! A temporary array for other variables. !
Expand Down Expand Up @@ -204,16 +208,16 @@ subroutine DOME_initialize_sponges(G, GV, US, tv, PF, CSp)
! depth space for Boussinesq or non-Boussinesq models.
eta(i,j,1) = 0.0
do k=2,nz
! eta(i,j,K)=max(H0(k), -G%bathyT(i,j), GV%Angstrom_Z*(nz-k+1) - G%bathyT(i,j))
e_dense = -G%bathyT(i,j)
! eta(i,j,K) = max(H0(k), -depth_tot(i,j), GV%Angstrom_Z*(nz-k+1) - depth_tot(i,j))
e_dense = -depth_tot(i,j)
if (e_dense >= H0(k)) then ; eta(i,j,K) = e_dense
else ; eta(i,j,K) = H0(k) ; endif
if (eta(i,j,K) < GV%Angstrom_Z*(nz-k+1) - G%bathyT(i,j)) &
eta(i,j,K) = GV%Angstrom_Z*(nz-k+1) - G%bathyT(i,j)
if (eta(i,j,K) < GV%Angstrom_Z*(nz-k+1) - depth_tot(i,j)) &
eta(i,j,K) = GV%Angstrom_Z*(nz-k+1) - depth_tot(i,j)
enddo
eta(i,j,nz+1) = -G%bathyT(i,j)
eta(i,j,nz+1) = -depth_tot(i,j)

if (G%bathyT(i,j) > min_depth) then
if (depth_tot(i,j) > min_depth) then
Idamp(i,j) = damp / 86400.0
else ; Idamp(i,j) = 0.0 ; endif
enddo ; enddo
Expand Down
30 changes: 18 additions & 12 deletions src/user/ISOMIP_initialization.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -128,12 +128,14 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth, US)
end subroutine ISOMIP_initialize_topography

!> Initialization of thicknesses
subroutine ISOMIP_initialize_thickness ( h, G, GV, US, param_file, tv, just_read_params)
subroutine ISOMIP_initialize_thickness ( h, depth_tot, G, GV, US, param_file, tv, just_read_params)
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
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2].
real, dimension(SZI_(G),SZJ_(G)), &
intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
!! to parse for model parameter values.
type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any
Expand Down Expand Up @@ -206,7 +208,7 @@ subroutine ISOMIP_initialize_thickness ( h, G, GV, US, param_file, tv, just_read

! Calculate thicknesses
do j=js,je ; do i=is,ie
eta1D(nz+1) = -G%bathyT(i,j)
eta1D(nz+1) = -depth_tot(i,j)
do k=nz,1,-1
eta1D(k) = e0(k)
if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_Z)) then
Expand All @@ -221,7 +223,7 @@ subroutine ISOMIP_initialize_thickness ( h, G, GV, US, param_file, tv, just_read
case ( REGRIDDING_ZSTAR, REGRIDDING_SIGMA_SHELF_ZSTAR ) ! Initial thicknesses for z coordinates
if (just_read) return ! All run-time parameters have been read, so return.
do j=js,je ; do i=is,ie
eta1D(nz+1) = -G%bathyT(i,j)
eta1D(nz+1) = -depth_tot(i,j)
do k=nz,1,-1
eta1D(k) = -G%max_depth * real(k-1) / real(nz)
if (eta1D(k) < (eta1D(k+1) + min_thickness)) then
Expand All @@ -236,7 +238,7 @@ subroutine ISOMIP_initialize_thickness ( h, G, GV, US, param_file, tv, just_read
case ( REGRIDDING_SIGMA ) ! Initial thicknesses for sigma coordinates
if (just_read) return ! All run-time parameters have been read, so return.
do j=js,je ; do i=is,ie
h(i,j,:) = GV%Z_to_H * G%bathyT(i,j) / dfloat(nz)
h(i,j,:) = GV%Z_to_H * depth_tot(i,j) / dfloat(nz)
enddo ; enddo

case default
Expand All @@ -248,14 +250,16 @@ subroutine ISOMIP_initialize_thickness ( h, G, GV, US, param_file, tv, just_read
end subroutine ISOMIP_initialize_thickness

!> Initial values for temperature and salinity
subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, US, param_file, &
subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, depth_tot, G, GV, US, param_file, &
eqn_of_state, just_read_params)
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
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T !< Potential temperature [degC]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S !< Salinity [ppt]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
real, dimension(SZI_(G),SZJ_(G)), intent(in) :: depth_tot !< The nominal total bottom-to-top
!! depth of the ocean [Z ~> m]
type(param_file_type), intent(in) :: param_file !< Parameter file structure
type(EOS_type), pointer :: eqn_of_state !< Equation of state structure
logical, optional, intent(in) :: just_read_params !< If present and true, this call will
Expand Down Expand Up @@ -315,7 +319,7 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, US, param_fi
dS_dz = (s_sur - s_bot) / G%max_depth
dT_dz = (t_sur - t_bot) / G%max_depth
do j=js,je ; do i=is,ie
xi0 = -G%bathyT(i,j)
xi0 = -depth_tot(i,j)
do k = nz,1,-1
xi0 = xi0 + 0.5 * h(i,j,k) * GV%H_to_Z ! Depth in middle of layer
S(i,j,k) = S_sur + dS_dz * xi0
Expand Down Expand Up @@ -420,7 +424,7 @@ end subroutine ISOMIP_initialize_temperature_salinity
!> Sets up the the inverse restoration time (Idamp), and
! the values towards which the interface heights and an arbitrary
! number of tracers should be restored within each sponge.
subroutine ISOMIP_initialize_sponges(G, GV, US, tv, PF, use_ALE, CSp, ACSp)
subroutine ISOMIP_initialize_sponges(G, GV, US, tv, depth_tot, PF, use_ALE, CSp, ACSp)
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
Expand All @@ -429,6 +433,8 @@ subroutine ISOMIP_initialize_sponges(G, GV, US, tv, PF, use_ALE, CSp, ACSp)
!! fields, potential temperature and
!! salinity or mixed layer density.
!! Absent fields have NULL ptrs.
real, dimension(SZI_(G),SZJ_(G)), &
intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
type(param_file_type), intent(in) :: PF !< A structure indicating the
!! open file to parse for model
!! parameter values.
Expand Down Expand Up @@ -508,7 +514,7 @@ subroutine ISOMIP_initialize_sponges(G, GV, US, tv, PF, use_ALE, CSp, ACSp)
! and mask2dT is 1.

do j=js,je ; do i=is,ie
if (G%bathyT(i,j) <= min_depth) then
if (depth_tot(i,j) <= min_depth) then
Idamp(i,j) = 0.0
elseif (G%geoLonT(i,j) >= 790.0 .AND. G%geoLonT(i,j) <= 800.0) then
dummy1 = (G%geoLonT(i,j)-790.0)/(800.0-790.0)
Expand Down Expand Up @@ -549,7 +555,7 @@ subroutine ISOMIP_initialize_sponges(G, GV, US, tv, PF, use_ALE, CSp, ACSp)

! Calculate thicknesses
do j=js,je ; do i=is,ie
eta1D(nz+1) = -G%bathyT(i,j)
eta1D(nz+1) = -depth_tot(i,j)
do k=nz,1,-1
eta1D(k) = e0(k)
if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_Z)) then
Expand All @@ -563,7 +569,7 @@ subroutine ISOMIP_initialize_sponges(G, GV, US, tv, PF, use_ALE, CSp, ACSp)

case ( REGRIDDING_ZSTAR, REGRIDDING_SIGMA_SHELF_ZSTAR ) ! Initial thicknesses for z coordinates
do j=js,je ; do i=is,ie
eta1D(nz+1) = -G%bathyT(i,j)
eta1D(nz+1) = -depth_tot(i,j)
do k=nz,1,-1
eta1D(k) = -G%max_depth * real(k-1) / real(nz)
if (eta1D(k) < (eta1D(k+1) + min_thickness)) then
Expand All @@ -577,7 +583,7 @@ subroutine ISOMIP_initialize_sponges(G, GV, US, tv, PF, use_ALE, CSp, ACSp)

case ( REGRIDDING_SIGMA ) ! Initial thicknesses for sigma coordinates
do j=js,je ; do i=is,ie
h(i,j,:) = GV%Z_to_H * (G%bathyT(i,j) / dfloat(nz))
h(i,j,:) = GV%Z_to_H * (depth_tot(i,j) / dfloat(nz))
enddo ; enddo

case default
Expand All @@ -593,7 +599,7 @@ subroutine ISOMIP_initialize_sponges(G, GV, US, tv, PF, use_ALE, CSp, ACSp)
dS_dz = (s_sur - s_bot) / G%max_depth
dT_dz = (t_sur - t_bot) / G%max_depth
do j=js,je ; do i=is,ie
xi0 = -G%bathyT(i,j)
xi0 = -depth_tot(i,j)
do k = nz,1,-1
xi0 = xi0 + 0.5 * h(i,j,k) * GV%H_to_Z ! Depth in middle of layer
S(i,j,k) = S_sur + dS_dz * xi0
Expand Down
Loading

0 comments on commit 3cda6fd

Please sign in to comment.