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.

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

Fix to dumbbell initialization in layer mode #160

Merged
merged 10 commits into from
Jul 18, 2022
2 changes: 1 addition & 1 deletion src/initialization/MOM_state_initialization.F90
Original file line number Diff line number Diff line change
Expand Up @@ -585,7 +585,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
case ("USER"); call user_initialize_sponges(G, GV, use_temperature, tv, PF, sponge_CSp, h)
case ("BFB"); call BFB_initialize_sponges_southonly(G, GV, US, use_temperature, tv, depth_tot, PF, &
sponge_CSp, h)
case ("DUMBBELL"); call dumbbell_initialize_sponges(G, GV, US, tv, depth_tot, PF, useALE, &
case ("DUMBBELL"); call dumbbell_initialize_sponges(G, GV, US, tv, h, depth_tot, PF, useALE, &
sponge_CSp, ALE_sponge_CSp)
case ("phillips"); call Phillips_initialize_sponges(G, GV, US, tv, PF, sponge_CSp, h)
case ("dense"); call dense_water_initialize_sponges(G, GV, US, tv, depth_tot, PF, useALE, &
Expand Down
69 changes: 64 additions & 5 deletions src/user/dumbbell_initialization.F90
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ module dumbbell_initialization
use MOM_verticalGrid, only : verticalGrid_type
use regrid_consts, only : coordinateMode, DEFAULT_COORDINATE_MODE
use regrid_consts, only : REGRIDDING_LAYER, REGRIDDING_ZSTAR
use regrid_consts, only : REGRIDDING_RHO, REGRIDDING_SIGMA
use regrid_consts, only : REGRIDDING_RHO, REGRIDDING_SIGMA, REGRIDDING_HYCOM1
use MOM_ALE_sponge, only : ALE_sponge_CS, set_up_ALE_sponge_field, initialize_ALE_sponge

implicit none ; private
Expand Down Expand Up @@ -112,10 +112,14 @@ subroutine dumbbell_initialize_thickness ( h, depth_tot, G, GV, US, param_file,
real :: S_range ! The range of salinities in this test case [S ~> ppt]
real :: S_light, S_dense ! The lightest and densest salinities in the sponges [S ~> ppt].
real :: eta_IC_quanta ! The granularity of quantization of intial interface heights [Z-1 ~> m-1].
logical :: dbrotate ! If true, rotate the domain.
logical :: use_ALE ! True if ALE is being used, False if in layered mode

! This include declares and sets the variable "version".
# include "version_variable.h"
character(len=20) :: verticalCoordinate
integer :: i, j, k, is, ie, js, je, nz
real :: x, y

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke

Expand All @@ -128,6 +132,8 @@ subroutine dumbbell_initialize_thickness ( h, depth_tot, G, GV, US, param_file,
units='m', default=1.0e-3, scale=US%m_to_Z, do_not_log=just_read)
call get_param(param_file, mdl,"REGRIDDING_COORDINATE_MODE", verticalCoordinate, &
default=DEFAULT_COORDINATE_MODE, do_not_log=just_read)
call get_param(param_file, mdl, "USE_REGRIDDING", use_ALE, do_not_log = .true.)
if (.not. use_ALE) verticalCoordinate = "LAYER"

! WARNING: this routine specifies the interface heights so that the last layer
! is vanished, even at maximum depth. In order to have a uniform
Expand All @@ -141,8 +147,36 @@ subroutine dumbbell_initialize_thickness ( h, depth_tot, G, GV, US, param_file,
!enddo

select case ( coordinateMode(verticalCoordinate) )
case ( REGRIDDING_LAYER) ! Initial thicknesses for isopycnal coordinates
call get_param(param_file, mdl, "DUMBBELL_ROTATION", dbrotate, &
'Logical for rotation of dumbbell domain.',&
units='nondim', default=.false., do_not_log=just_read)
do j=js,je
do i=is,ie
! Compute normalized zonal coordinates (x,y=0 at center of domain)
if (dbrotate) then
! This is really y in the rotated case
x = G%geoLatT(i,j)
else
x = G%geoLonT(i,j)
endif
eta1D(1) = 0.0
eta1D(nz+1) = -depth_tot(i,j)
if (x<0.0) then
do k=nz,2, -1
eta1D(k) = eta1D(k+1) + min_thickness
enddo
else
do k=2,nz
eta1D(k) = eta1D(k-1) - min_thickness
enddo
endif
do k=1,nz
h(i,j,k) = GV%Z_to_H * (eta1D(k) - eta1D(k+1))
enddo
enddo; enddo

case ( REGRIDDING_LAYER, REGRIDDING_RHO ) ! Initial thicknesses for isopycnal coordinates
case ( REGRIDDING_RHO, REGRIDDING_HYCOM1) ! Initial thicknesses for isopycnal coordinates
call get_param(param_file, mdl, "INITIAL_SSS", S_surf, &
units='1e-3', default=34., scale=US%ppt_to_S, do_not_log=.true.)
call get_param(param_file, mdl, "INITIAL_S_RANGE", S_range, &
Expand Down Expand Up @@ -231,12 +265,18 @@ subroutine dumbbell_initialize_temperature_salinity ( T, S, h, G, GV, US, param_
real :: x ! The fractional position in the domain [nondim]
real :: dblen ! The size of the dumbbell test case [axis_units]
logical :: dbrotate ! If true, rotate the domain.
logical :: use_ALE ! If false, use layer mode.
character(len=20) :: verticalCoordinate, density_profile

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke

T_surf = 20.0*US%degC_to_C

! layer mode
call get_param(param_file, mdl, "USE_REGRIDDING", use_ALE, do_not_log = .true.)
if (.not. use_ALE) call MOM_error(FATAL, "dumbbell_initialize_temperature_salinity: "//&
"Please use 'fit' for 'TS_CONFIG' in the LAYER mode.")

call get_param(param_file, mdl, "REGRIDDING_COORDINATE_MODE", verticalCoordinate, &
default=DEFAULT_COORDINATE_MODE, do_not_log=just_read)
call get_param(param_file, mdl, "INITIAL_DENSITY_PROFILE", density_profile, &
Expand Down Expand Up @@ -288,11 +328,12 @@ subroutine dumbbell_initialize_temperature_salinity ( T, S, h, G, GV, US, param_
end subroutine dumbbell_initialize_temperature_salinity

!> Initialize the restoring sponges for the dumbbell test case
subroutine dumbbell_initialize_sponges(G, GV, US, tv, depth_tot, param_file, use_ALE, CSp, ACSp)
subroutine dumbbell_initialize_sponges(G, GV, US, tv, h_in, depth_tot, param_file, use_ALE, CSp, ACSp)
type(ocean_grid_type), intent(in) :: G !< Horizontal grid control structure
type(verticalGrid_type), intent(in) :: GV !< Vertical grid control structure
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_in !< Layer thickness [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 !< Parameter file structure
Expand All @@ -306,6 +347,7 @@ subroutine dumbbell_initialize_sponges(G, GV, US, tv, depth_tot, param_file, use
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)) :: S ! sponge salinities [S ~> ppt]
real, dimension(SZK_(GV)+1) :: eta1D ! interface positions for ALE sponge
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: eta ! A temporary array for interface heights [Z ~> m].

integer :: i, j, k, nz
real :: x ! The fractional position in the domain [nondim]
Expand Down Expand Up @@ -404,9 +446,26 @@ subroutine dumbbell_initialize_sponges(G, GV, US, tv, depth_tot, param_file, use
enddo
endif
enddo ; enddo
endif

if (associated(tv%S)) call set_up_ALE_sponge_field(S, G, GV, tv%S, ACSp)
else
do j=G%jsc,G%jec ; do i=G%isc,G%iec
eta(i,j,1) = 0.0
do k=2,nz
eta(i,j,k) = eta(i,j,k-1)- GV%H_to_Z * h_in(i,j,k-1)
enddo
eta(i,j,nz+1) = -depth_tot(i,j)
do k=1,nz
S(i,j,k)= tv%S(i,j,k)
enddo
enddo ; enddo

! This call sets up the damping rates and interface heights.
! This sets the inverse damping timescale fields in the sponges. !
call initialize_sponge(Idamp, eta, G, param_file, CSp, GV)

! The remaining calls to set_up_sponge_field can be in any order. !
if ( associated(tv%S) ) call set_up_sponge_field(S, tv%S, G, GV, nz, CSp)
endif

end subroutine dumbbell_initialize_sponges

Expand Down