Skip to content

+(*)Better initialization of non-Boussinesq models #349

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

Merged
marshallward merged 5 commits into from
May 5, 2023
Merged

+(*)Better initialization of non-Boussinesq models #349

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
9919acf
+Remove optional argument eta_to_m from find_eta
Hallberg-NOAA Apr 4, 2023
aaf04fc
+Initialize thicknesses in height units
Hallberg-NOAA Apr 7, 2023
2d3c74e
+Add the new overloaded interface dz_to_thickness
Hallberg-NOAA Apr 7, 2023
6d192ca
(*)Improve non-Boussinesq initialization
Hallberg-NOAA Apr 8, 2023
916fbfa
(*)Use dz_to_thickness in 4 user modules
Hallberg-NOAA Apr 8, 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
12 changes: 9 additions & 3 deletions src/ALE/MOM_ALE.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1463,17 +1463,23 @@ subroutine ALE_writeCoordinateFile( CS, GV, directory )
end subroutine ALE_writeCoordinateFile

!> Set h to coordinate values for fixed coordinate systems
subroutine ALE_initThicknessToCoord( CS, G, GV, h )
subroutine ALE_initThicknessToCoord( CS, G, GV, h, height_units )
type(ALE_CS), intent(inout) :: CS !< module control structure
type(ocean_grid_type), intent(in) :: G !< module grid structure
type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: h !< layer thickness [H ~> m or kg m-2]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: h !< layer thickness in thickness units
!! [H ~> m or kg m-2] or height units [Z ~> m]
logical, optional, intent(in) :: height_units !< If present and true, the
!! thicknesses are in height units

! Local variables
real :: scale ! A scaling value for the thicknesses [nondim] or [H Z-1 ~> nondim or kg m-3]
integer :: i, j

scale = GV%Z_to_H
if (present(height_units)) then ; if (height_units) scale = 1.0 ; endif
do j = G%jsd,G%jed ; do i = G%isd,G%ied
h(i,j,:) = GV%Z_to_H * getStaticThickness( CS%regridCS, 0., G%bathyT(i,j)+G%Z_ref )
h(i,j,:) = scale * getStaticThickness( CS%regridCS, 0., G%bathyT(i,j)+G%Z_ref )
enddo ; enddo

end subroutine ALE_initThicknessToCoord
Expand Down
239 changes: 204 additions & 35 deletions src/core/MOM_interface_heights.F90
View file
Open in desktop

Large diffs are not rendered by default.

1 change: 1 addition & 0 deletions src/diagnostics/MOM_obsolete_params.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,7 @@ subroutine find_obsolete_params(param_file)
hint="Instead use OBC_SEGMENT_xxx_VELOCITY_NUDGING_TIMESCALES.")
enddo

call obsolete_logical(param_file, "CONVERT_THICKNESS_UNITS", .true.)
call obsolete_logical(param_file, "MASK_MASSLESS_TRACERS", .false.)

call obsolete_logical(param_file, "SALT_REJECT_BELOW_ML", .false.)
Expand Down
548 changes: 268 additions & 280 deletions src/initialization/MOM_state_initialization.F90
View file
Open in desktop

Large diffs are not rendered by default.

14 changes: 12 additions & 2 deletions src/initialization/MOM_tracer_initialization_from_Z.F90
View file
Open in desktop
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
17 changes: 5 additions & 12 deletions src/tracer/MOM_tracer_Z_init.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -556,29 +556,24 @@ end function find_limited_slope

!> This subroutine determines the potential temperature and salinity that
!! is consistent with the target density using provided initial guess
subroutine determine_temperature(temp, salt, R_tgt, EOS, p_ref, niter, h, k_start, G, GV, US, &
PF, just_read, h_massless)
subroutine determine_temperature(temp, salt, R_tgt, EOS, p_ref, niter, k_start, G, GV, US, PF, &
just_read)
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(inout) :: temp !< potential temperature [C ~> degC]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
intent(inout) :: salt !< salinity [S ~> ppt]
real, dimension(SZK_(GV)), intent(in) :: R_tgt !< desired potential density [R ~> kg m-3].
type(EOS_type), intent(in) :: EOS !< seawater equation of state control structure
type(EOS_type), intent(in) :: EOS !< seawater equation of state control structure
real, intent(in) :: p_ref !< reference pressure [R L2 T-2 ~> Pa].
integer, intent(in) :: niter !< maximum number of iterations
integer, intent(in) :: k_start !< starting index (i.e. below the buffer layer)
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
intent(in) :: h !< layer thickness, used only to avoid working on
!! massless layers [H ~> m or kg m-2]
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(param_file_type), intent(in) :: PF !< A structure indicating the open file
!! to parse for model parameter values.
logical, intent(in) :: just_read !< If true, this call will only read
!! parameters without changing T or S.
real, optional, intent(in) :: h_massless !< A threshold below which a layer is
!! determined to be massless [H ~> m or kg m-2]

! Local variables (All of which need documentation!)
real, dimension(SZI_(G),SZK_(GV)) :: &
Expand All @@ -587,7 +582,6 @@ subroutine determine_temperature(temp, salt, R_tgt, EOS, p_ref, niter, h, k_star
dT, & ! An estimated change in temperature before bounding [C ~> degC]
dS, & ! An estimated change in salinity before bounding [S ~> ppt]
rho, & ! Layer densities with the current estimate of temperature and salinity [R ~> kg m-3]
hin, & ! A 2D copy of the layer thicknesses [H ~> m or kg m-2]
drho_dT, & ! Partial derivative of density with temperature [R C-1 ~> kg m-3 degC-1]
drho_dS ! Partial derivative of density with salinity [R S-1 ~> kg m-3 ppt-1]
real, dimension(SZI_(G)) :: press ! Reference pressures [R L2 T-2 ~> Pa]
Expand Down Expand Up @@ -675,7 +669,6 @@ subroutine determine_temperature(temp, salt, R_tgt, EOS, p_ref, niter, h, k_star
dS(:,:) = 0. ! Needs to be zero everywhere since there is a maxval(abs(dS)) later...
T(:,:) = temp(:,j,:)
S(:,:) = salt(:,j,:)
hin(:,:) = h(:,j,:)
dT(:,:) = 0.0
adjust_salt = .true.
iter_loop: do itt = 1,niter
Expand All @@ -685,7 +678,7 @@ subroutine determine_temperature(temp, salt, R_tgt, EOS, p_ref, niter, h, k_star
EOS, EOSdom )
enddo
do k=k_start,nz ; do i=is,ie
! if (abs(rho(i,k)-R_tgt(k))>tol_rho .and. hin(i,k)>h_massless .and. abs(T(i,k)-land_fill) < epsln) then
! if (abs(rho(i,k)-R_tgt(k))>tol_rho .and. abs(T(i,k)-land_fill) < epsln) then
if (abs(rho(i,k)-R_tgt(k))>tol_rho) then
if (.not.fit_together) then
dT(i,k) = max(min((R_tgt(k)-rho(i,k)) / drho_dT(i,k), max_t_adj), -max_t_adj)
Expand Down Expand Up @@ -713,7 +706,7 @@ subroutine determine_temperature(temp, salt, R_tgt, EOS, p_ref, niter, h, k_star
EOS, EOSdom )
enddo
do k=k_start,nz ; do i=is,ie
! if (abs(rho(i,k)-R_tgt(k))>tol_rho .and. hin(i,k)>h_massless .and. abs(T(i,k)-land_fill) < epsln ) then
! if (abs(rho(i,k)-R_tgt(k))>tol_rho .and. abs(T(i,k)-land_fill) < epsln ) then
if (abs(rho(i,k)-R_tgt(k)) > tol_rho) then
dS(i,k) = max(min((R_tgt(k)-rho(i,k)) / drho_dS(i,k), max_s_adj), -max_s_adj)
S(i,k) = max(min(S(i,k)+dS(i,k), S_max), S_min)
Expand Down
60 changes: 35 additions & 25 deletions src/user/DOME2d_initialization.F90
View file
Open in desktop
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 @@ -98,7 +99,7 @@ subroutine DOME2d_initialize_thickness ( h, depth_tot, G, GV, US, param_file, ju
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].
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
Expand Down Expand Up @@ -158,16 +159,16 @@ subroutine DOME2d_initialize_thickness ( h, depth_tot, G, GV, US, param_file, ju
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
h(i,j,k) = GV%Angstrom_Z
else
h(i,j,k) = GV%Z_to_H * (eta1D(k) - eta1D(k+1))
h(i,j,k) = eta1D(k) - eta1D(k+1)
endif
enddo

x = ( G%geoLonT(i,j) - G%west_lon ) / G%len_lon
if ( x <= dome2d_width_bay ) then
h(i,j,1:nz-1) = GV%Angstrom_H
h(i,j,nz) = GV%Z_to_H * dome2d_depth_bay * G%max_depth - (nz-1) * GV%Angstrom_H
h(i,j,1:nz-1) = GV%Angstrom_Z
h(i,j,nz) = dome2d_depth_bay * G%max_depth - (nz-1) * GV%Angstrom_Z
endif

enddo ; enddo
Expand All @@ -180,16 +181,16 @@ subroutine DOME2d_initialize_thickness ( h, depth_tot, G, GV, US, param_file, ju
! eta1D(k) = e0(k)
! if (eta1D(k) < (eta1D(k+1) + min_thickness)) then
! eta1D(k) = eta1D(k+1) + min_thickness
! h(i,j,k) = GV%Z_to_H * min_thickness
! h(i,j,k) = min_thickness
! else
! h(i,j,k) = GV%Z_to_H * (eta1D(k) - eta1D(k+1))
! h(i,j,k) = eta1D(k) - eta1D(k+1)
! endif
! enddo
!
! x = G%geoLonT(i,j) / G%len_lon
! if ( x <= dome2d_width_bay ) then
! h(i,j,1:nz-1) = GV%Z_to_H * min_thickness
! h(i,j,nz) = GV%Z_to_H * (dome2d_depth_bay * G%max_depth - (nz-1) * min_thickness)
! h(i,j,1:nz-1) = min_thickness
! h(i,j,nz) = dome2d_depth_bay * G%max_depth - (nz-1) * min_thickness
! endif
!
! enddo ; enddo
Expand All @@ -202,16 +203,16 @@ subroutine DOME2d_initialize_thickness ( h, depth_tot, G, GV, US, param_file, ju
eta1D(k) = e0(k)
if (eta1D(k) < (eta1D(k+1) + min_thickness)) then
eta1D(k) = eta1D(k+1) + min_thickness
h(i,j,k) = GV%Z_to_H * min_thickness
h(i,j,k) = min_thickness
else
h(i,j,k) = GV%Z_to_H * (eta1D(k) - eta1D(k+1))
h(i,j,k) = eta1D(k) - eta1D(k+1)
endif
enddo
enddo ; enddo

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

case default
Expand All @@ -225,11 +226,11 @@ end subroutine DOME2d_initialize_thickness

!> Initialize temperature and salinity in the 2d DOME configuration
subroutine DOME2d_initialize_temperature_salinity ( T, S, h, G, GV, US, param_file, just_read)
type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T !< Potential temperature [C ~> degC]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S !< Salinity [S ~> 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),SZK_(GV)), intent(out) :: T !< Potential temperature [C ~> degC]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S !< Salinity [S ~> ppt]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [Z ~> m]
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(param_file_type), intent(in) :: param_file !< Parameter file structure
logical, intent(in) :: just_read !< If true, this call will
Expand Down Expand Up @@ -287,7 +288,7 @@ subroutine DOME2d_initialize_temperature_salinity ( T, S, h, G, GV, US, param_fi
do j=js,je ; do i=is,ie
xi0 = 0.0
do k = 1,nz
xi1 = xi0 + (GV%H_to_Z * h(i,j,k)) / G%max_depth
xi1 = xi0 + h(i,j,k) / G%max_depth
S(i,j,k) = S_surf + 0.5 * S_range * (xi0 + xi1)
xi0 = xi1
enddo
Expand All @@ -298,7 +299,7 @@ subroutine DOME2d_initialize_temperature_salinity ( T, S, h, G, GV, US, param_fi
do j=js,je ; do i=is,ie
xi0 = 0.0
do k = 1,nz
xi1 = xi0 + (GV%H_to_Z * h(i,j,k)) / G%max_depth
xi1 = xi0 + h(i,j,k) / G%max_depth
S(i,j,k) = S_surf + 0.5 * S_range * (xi0 + xi1)
xi0 = xi1
enddo
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
6 changes: 3 additions & 3 deletions src/user/DOME_initialization.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -105,7 +105,7 @@ subroutine DOME_initialize_thickness(h, depth_tot, G, GV, param_file, just_read)
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].
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
Expand Down Expand Up @@ -141,9 +141,9 @@ subroutine DOME_initialize_thickness(h, depth_tot, G, GV, param_file, just_read)
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
h(i,j,k) = GV%Angstrom_Z
else
h(i,j,k) = GV%Z_to_H * (eta1D(K) - eta1D(K+1))
h(i,j,k) = eta1D(K) - eta1D(K+1)
endif
enddo
enddo ; enddo
Expand Down
Loading