Skip to content

Commit

Permalink
(*)Use dz_to_thickness in 4 user modules
Browse files Browse the repository at this point in the history 
  Use dz_to_thickness to convert vertical distances to layer thicknesses in the
sponge initialization routines in the DOME2d_initialization,
ISOMIP_initialization, dumbbell_initialization and dense_water_initialization
modules, and also in MOM_initialize_tracer_from_Z.  For the user modules,
the presence or absence of an equation of state is known and handled properly,
but MOM_initialize_tracer_from_Z works with the generic tracer code and it
it outside of the scope of MOM6 code to provide any information about the
equation of state or the state variables that would be needed to initialize
a non-Boussinesq model properly from a depth-space input file.  For now we
are doing the best we can, but this should be revisited.  All examples in
existing test cases are bitwise identical, but answers could change (and be
improved) in any non-Boussinesq variants of the relevant test cases.
  • Loading branch information
@Hallberg-NOAA @marshallward
Hallberg-NOAA authored and marshallward committed May 5, 2023
1 parent fb5f4d7 commit debe45e
Show file tree
Hide file tree
Showing 5 changed files with 105 additions and 52 deletions.
14 changes: 12 additions & 2 deletions src/initialization/MOM_tracer_initialization_from_Z.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@ module MOM_tracer_initialization_from_Z
use MOM_file_parser, only : get_param, param_file_type, log_version
use MOM_grid, only : ocean_grid_type
use MOM_horizontal_regridding, only : myStats, horiz_interp_and_extrap_tracer
use MOM_interface_heights, only : dz_to_thickness_simple
use MOM_remapping, only : remapping_CS, initialize_remapping
use MOM_unit_scaling, only : unit_scale_type
use MOM_verticalGrid, only : verticalGrid_type
Expand Down Expand Up @@ -75,10 +76,12 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_
real, allocatable, dimension(:), target :: z_in ! Cell center depths for input data [Z ~> m]

! Local variables for ALE remapping
real, dimension(:,:,:), allocatable :: hSrc ! Source thicknesses [H ~> m or kg m-2].
real, dimension(:,:,:), allocatable :: dzSrc ! Source thicknesses in height units [Z ~> m]
real, dimension(:,:,:), allocatable :: hSrc ! Source thicknesses [H ~> m or kg m-2]
real, dimension(:), allocatable :: h1 ! A 1-d column of source thicknesses [Z ~> m].
real :: zTopOfCell, zBottomOfCell, z_bathy ! Heights [Z ~> m].
type(remapping_CS) :: remapCS ! Remapping parameters and work arrays
type(verticalGrid_type) :: GV_loc ! A temporary vertical grid structure

real :: missing_value ! A value indicating that there is no valid input data at this point [CU ~> conc]
integer :: nPoints ! The number of valid input data points in a column
Expand Down Expand Up @@ -180,6 +183,7 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_
call cpu_clock_begin(id_clock_ALE)
! First we reserve a work space for reconstructions of the source data
allocate( h1(kd) )
allocate( dzSrc(isd:ied,jsd:jed,kd) )
allocate( hSrc(isd:ied,jsd:jed,kd) )
! Set parameters for reconstructions
call initialize_remapping( remapCS, remapScheme, boundary_extrapolation=.false., answer_date=remap_answer_date )
Expand All @@ -204,12 +208,18 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_
else
tr(i,j,:) = 0.
endif ! mask2dT
hSrc(i,j,:) = GV%Z_to_H * h1(:)
dzSrc(i,j,:) = h1(:)
enddo ; enddo

! Equation of state data is not available, so a simpler rescaling will have to suffice,
! but it might be problematic in non-Boussinesq mode.
GV_loc = GV ; GV_loc%ke = kd
call dz_to_thickness_simple(dzSrc, hSrc, G, GV_loc, US)

call ALE_remap_scalar(remapCS, G, GV, kd, hSrc, tr_z, h, tr, all_cells=.false., answer_date=remap_answer_date )

deallocate( hSrc )
deallocate( dzSrc )
deallocate( h1 )

do k=1,nz
Expand Down
24 changes: 17 additions & 7 deletions src/user/DOME2d_initialization.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@ module DOME2d_initialization
use MOM_file_parser, only : get_param, log_version, param_file_type
use MOM_get_input, only : directories
use MOM_grid, only : ocean_grid_type
use MOM_interface_heights, only : dz_to_thickness, dz_to_thickness_simple
use MOM_sponge, only : sponge_CS, set_up_sponge_field, initialize_sponge
use MOM_unit_scaling, only : unit_scale_type
use MOM_variables, only : thermo_var_ptrs
Expand Down Expand Up @@ -373,7 +374,8 @@ subroutine DOME2d_initialize_sponges(G, GV, US, tv, depth_tot, param_file, use_A
! Local variables
real :: T(SZI_(G),SZJ_(G),SZK_(GV)) ! A temporary array for temp [C ~> degC]
real :: S(SZI_(G),SZJ_(G),SZK_(GV)) ! A temporary array for salt [S ~> ppt]
real :: h(SZI_(G),SZJ_(G),SZK_(GV)) ! A temporary array for thickness [H ~> m or kg m-2].
real :: dz(SZI_(G),SZJ_(G),SZK_(GV)) ! A temporary array for thickness in height units [Z ~> m]
real :: h(SZI_(G),SZJ_(G),SZK_(GV)) ! A temporary array for thickness [H ~> m or kg m-2]
real :: eta(SZI_(G),SZJ_(G),SZK_(GV)+1) ! A temporary array for interface heights [Z ~> m]
real :: Idamp(SZI_(G),SZJ_(G)) ! The sponge damping rate [T-1 ~> s-1]
real :: S_ref ! Reference salinity within the surface layer [S ~> ppt]
Expand Down Expand Up @@ -478,30 +480,38 @@ subroutine DOME2d_initialize_sponges(G, GV, US, tv, depth_tot, param_file, use_A
eta1D(k) = e0(k)
if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_Z)) then
eta1D(k) = eta1D(k+1) + GV%Angstrom_Z
h(i,j,k) = GV%Angstrom_H
dz(i,j,k) = GV%Angstrom_Z
else
h(i,j,k) = GV%Z_to_H * (eta1D(k) - eta1D(k+1))
dz(i,j,k) = eta1D(k) - eta1D(k+1)
endif
enddo
enddo ; enddo
! Store the grid on which the T/S sponge data will reside
call initialize_ALE_sponge(Idamp, G, GV, param_file, ACSp, h, nz)

! Construct temperature and salinity on the arbitrary grid
T(:,:,:) = 0.0 ; S(:,:,:) = 0.0
do j=js,je ; do i=is,ie
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
z = z + 0.5 * dz(i,j,k) ! Position of the center of layer k
! Use salinity stratification in the eastern sponge.
S(i,j,k) = S_surf - S_range_sponge * (z / G%max_depth)
! Use a constant salinity in the western sponge.
if ( ( G%geoLonT(i,j) - G%west_lon ) / G%len_lon < dome2d_west_sponge_width ) &
S(i,j,k) = S_ref + S_range
z = z + 0.5 * GV%H_to_Z * h(i,j,k) ! Position of the interface k
z = z + 0.5 * dz(i,j,k) ! Position of the interface k
enddo
enddo ; enddo

! Convert thicknesses from height units to thickness units
if (associated(tv%eqn_of_state)) then
call dz_to_thickness(dz, T, S, tv%eqn_of_state, h, G, GV, US)
else
call dz_to_thickness_simple(dz, h, G, GV, US, layer_mode=.true.)
endif

! Store damping rates and the grid on which the T/S sponge data will reside
call initialize_ALE_sponge(Idamp, G, GV, param_file, ACSp, h, nz)

if ( associated(tv%T) ) call set_up_ALE_sponge_field(T, G, GV, tv%T, ACSp, 'temp', &
sp_long_name='temperature', sp_unit='degC s-1')
if ( associated(tv%S) ) call set_up_ALE_sponge_field(S, G, GV, tv%S, ACSp, 'salt', &
Expand Down
63 changes: 34 additions & 29 deletions src/user/ISOMIP_initialization.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@ module ISOMIP_initialization
use MOM_file_parser, only : get_param, log_version, param_file_type
use MOM_get_input, only : directories
use MOM_grid, only : ocean_grid_type
use MOM_interface_heights, only : dz_to_thickness
use MOM_io, only : file_exists, MOM_read_data, slasher
use MOM_unit_scaling, only : unit_scale_type
use MOM_variables, only : thermo_var_ptrs
Expand Down Expand Up @@ -146,8 +147,7 @@ subroutine ISOMIP_initialize_thickness ( h, depth_tot, G, GV, US, param_file, tv
intent(out) :: h !< The thickness that is being initialized [Z ~> m]
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(param_file_type), intent(in) :: param_file !< A structure to parse for model parameter values
type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any
!! available thermodynamic fields, including
!! the eqn. of state.
Expand Down Expand Up @@ -440,27 +440,25 @@ end subroutine ISOMIP_initialize_temperature_salinity
! 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, 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
type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers
!! to any available thermodynamic
!! fields, potential temperature and
!! salinity or mixed layer density.
!! Absent fields have NULL ptrs.
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, 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.
logical, intent(in) :: use_ALE !< If true, indicates model is in ALE mode
type(sponge_CS), pointer :: CSp !< Layer-mode sponge structure
type(ALE_sponge_CS), pointer :: ACSp !< ALE-mode sponge structure
intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
type(param_file_type), intent(in) :: PF !< A structure to parse for model parameter values
logical, intent(in) :: use_ALE !< If true, indicates model is in ALE mode
type(sponge_CS), pointer :: CSp !< Layer-mode sponge structure
type(ALE_sponge_CS), pointer :: ACSp !< ALE-mode sponge structure
! Local variables
real :: T(SZI_(G),SZJ_(G),SZK_(GV)) ! A temporary array for temp [C ~> degC]
real :: S(SZI_(G),SZJ_(G),SZK_(GV)) ! A temporary array for salt [S ~> ppt]
! real :: RHO(SZI_(G),SZJ_(G),SZK_(GV)) ! A temporary array for RHO [R ~> kg m-3]
real :: h(SZI_(G),SZJ_(G),SZK_(GV)) ! A temporary array for thickness [H ~> m or kg m-2]
real :: dz(SZI_(G),SZJ_(G),SZK_(GV)) ! Sponge layer thicknesses in height units [Z ~> m]
real :: h(SZI_(G),SZJ_(G),SZK_(GV)) ! Sponge layer thicknesses [H ~> m or kg m-2]
real :: Idamp(SZI_(G),SZJ_(G)) ! The sponge damping rate [T-1 ~> s-1]
real :: TNUDG ! Nudging time scale [T ~> s]
real :: S_sur, S_bot ! Surface and bottom salinities in the sponge region [S ~> ppt]
Expand Down Expand Up @@ -582,9 +580,9 @@ subroutine ISOMIP_initialize_sponges(G, GV, US, tv, depth_tot, PF, use_ALE, CSp,
eta1D(k) = e0(k)
if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_Z)) then
eta1D(k) = eta1D(k+1) + GV%Angstrom_Z
h(i,j,k) = GV%Angstrom_H
dz(i,j,k) = GV%Angstrom_Z
else
h(i,j,k) = GV%Z_to_H*(eta1D(k) - eta1D(k+1))
dz(i,j,k) = eta1D(k) - eta1D(k+1)
endif
enddo
enddo ; enddo
Expand All @@ -596,16 +594,16 @@ subroutine ISOMIP_initialize_sponges(G, GV, US, tv, depth_tot, PF, use_ALE, CSp,
eta1D(k) = -G%max_depth * real(k-1) / real(nz)
if (eta1D(k) < (eta1D(k+1) + min_thickness)) then
eta1D(k) = eta1D(k+1) + min_thickness
h(i,j,k) = min_thickness * GV%Z_to_H
dz(i,j,k) = min_thickness
else
h(i,j,k) = GV%Z_to_H*(eta1D(k) - eta1D(k+1))
dz(i,j,k) = eta1D(k) - eta1D(k+1)
endif
enddo
enddo ; enddo

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

case default
Expand All @@ -614,21 +612,25 @@ subroutine ISOMIP_initialize_sponges(G, GV, US, tv, depth_tot, PF, use_ALE, CSp,

end select

! This call sets up the damping rates and interface heights.
! This sets the inverse damping timescale fields in the sponges.
call initialize_ALE_sponge(Idamp, G, GV, PF, ACSp, h, nz)

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 = -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
xi0 = xi0 + 0.5 * dz(i,j,k) ! Depth in middle of layer
S(i,j,k) = S_sur + dS_dz * xi0
T(i,j,k) = T_sur + dT_dz * xi0
xi0 = xi0 + 0.5 * h(i,j,k) * GV%H_to_Z ! Depth at top of layer
xi0 = xi0 + 0.5 * dz(i,j,k) ! Depth at top of layer
enddo
enddo ; enddo

! Convert thicknesses from height units to thickness units
if (associated(tv%eqn_of_state)) then
call dz_to_thickness(dz, T, S, tv%eqn_of_state, h, G, GV, US)
else
call MOM_error(FATAL, "The ISOMIP test case requires an equation of state.")
endif

! for debugging
!i=G%iec; j=G%jec
!do k = 1,nz
Expand All @@ -637,6 +639,9 @@ subroutine ISOMIP_initialize_sponges(G, GV, US, tv, depth_tot, PF, use_ALE, CSp,
! call MOM_mesg(mesg,5)
!enddo

! This call sets up the damping rates and interface heights in the sponges.
call initialize_ALE_sponge(Idamp, G, GV, PF, ACSp, h, nz)

! Now register all of the fields which are damped in the sponge. !
! By default, momentum is advected vertically within the sponge, but !
! momentum is typically not damped within the sponge. !
Expand Down
24 changes: 17 additions & 7 deletions src/user/dense_water_initialization.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@ module dense_water_initialization
use MOM_EOS, only : EOS_type
use MOM_error_handler, only : MOM_error, FATAL
use MOM_file_parser, only : get_param, param_file_type
use MOM_interface_heights, only : dz_to_thickness, dz_to_thickness_simple
use MOM_grid, only : ocean_grid_type
use MOM_sponge, only : sponge_CS
use MOM_unit_scaling, only : unit_scale_type
Expand Down Expand Up @@ -172,7 +173,8 @@ subroutine dense_water_initialize_sponges(G, GV, US, tv, depth_tot, param_file,
real :: east_sponge_width ! The fraction of the domain in which the eastern (outflow) sponge is active [nondim]

real, dimension(SZI_(G),SZJ_(G)) :: Idamp ! inverse damping timescale [T-1 ~> s-1]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h ! sponge thicknesses [H ~> m or kg m-2]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dz ! sponge layer thicknesses in height units [Z ~> m]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h ! sponge layer thicknesses [H ~> m or kg m-2]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: T ! sponge temperature [C ~> degC]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: S ! sponge salinity [S ~> ppt]
real, dimension(SZK_(GV)+1) :: e0, eta1D ! interface positions for ALE sponge [Z ~> m]
Expand Down Expand Up @@ -256,16 +258,14 @@ subroutine dense_water_initialize_sponges(G, GV, US, tv, depth_tot, param_file,
if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_Z)) then
! is this layer vanished?
eta1D(k) = eta1D(k+1) + GV%Angstrom_Z
h(i,j,k) = GV%Angstrom_H
dz(i,j,k) = GV%Angstrom_Z
else
h(i,j,k) = GV%Z_to_H * (eta1D(k) - eta1D(k+1))
dz(i,j,k) = eta1D(k) - eta1D(k+1)
endif
enddo
enddo
enddo

call initialize_ALE_sponge(Idamp, G, GV, param_file, ACSp, h, nz)

! construct temperature and salinity for the sponge
! start with initial condition
T(:,:,:) = T_ref
Expand All @@ -277,7 +277,7 @@ subroutine dense_water_initialize_sponges(G, GV, US, tv, depth_tot, param_file,
x = (G%geoLonT(i,j) - G%west_lon) / G%len_lon
do k = 1,nz
! nondimensional middle of layer
zmid = zi + 0.5 * h(i,j,k) / (GV%Z_to_H * G%max_depth)
zmid = zi + 0.5 * dz(i,j,k) / G%max_depth

if (x > (1. - east_sponge_width)) then
!if (zmid >= 0.9 * sill_frac) &
Expand All @@ -288,11 +288,21 @@ subroutine dense_water_initialize_sponges(G, GV, US, tv, depth_tot, param_file,
S(i,j,k) = S_ref + S_range * (zmid - mld) / (1.0 - mld)
endif

zi = zi + h(i,j,k) / (GV%Z_to_H * G%max_depth)
zi = zi + dz(i,j,k) / G%max_depth
enddo
enddo
enddo

! Convert thicknesses from height units to thickness units
if (associated(tv%eqn_of_state)) then
call dz_to_thickness(dz, T, S, tv%eqn_of_state, h, G, GV, US)
else
call dz_to_thickness_simple(dz, h, G, GV, US, layer_mode=.true.)
endif

! This call sets up the damping rates and interface heights in the sponges.
call initialize_ALE_sponge(Idamp, G, GV, param_file, ACSp, h, nz)

if ( associated(tv%T) ) call set_up_ALE_sponge_field(T, G, GV, tv%T, ACSp, 'temp', &
sp_long_name='temperature', sp_unit='degC s-1')
if ( associated(tv%S) ) call set_up_ALE_sponge_field(S, G, GV, tv%S, ACSp, 'salt', &
Expand Down
Loading

0 comments on commit debe45e

Please sign in to comment.