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

+Add 7 runtime parameters for user init modules #292

Merged
merged 5 commits into from
Jan 5, 2023
Merged

+Add 7 runtime parameters for user init modules #292

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
fcfa2f2
+Add runtime parameters for benchmark_initialization
Hallberg-NOAA Dec 23, 2022
dd7c787
+Add runtime parameters for DOME2d_initialization
Hallberg-NOAA Dec 23, 2022
5831dc6
+Add runtime parameters for sloshing_initialization
Hallberg-NOAA Dec 23, 2022
157b806
+Add runtime parameters for DOME_initialization
Hallberg-NOAA Dec 23, 2022
13fd0e1
Merge branch 'dev/gfdl' into explicit_user_params
marshallward Jan 5, 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
43 changes: 29 additions & 14 deletions src/user/DOME2d_initialization.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,9 +40,9 @@ module DOME2d_initialization
subroutine DOME2d_initialize_topography( D, G, param_file, max_depth )
type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type
real, dimension(G%isd:G%ied,G%jsd:G%jed), &
intent(out) :: D !< Ocean bottom depth in the units of depth_max
intent(out) :: D !< Ocean bottom depth [Z ~> m]
type(param_file_type), intent(in) :: param_file !< Parameter file structure
real, intent(in) :: max_depth !< Maximum ocean depth in arbitrary units
real, intent(in) :: max_depth !< Maximum ocean depth [Z ~> m]

! Local variables
real :: bay_depth ! Depth of shelf, as fraction of basin depth [nondim]
Expand Down Expand Up @@ -239,6 +239,8 @@ subroutine DOME2d_initialize_temperature_salinity ( T, S, h, G, GV, US, param_fi
real :: delta_S ! Change in salinity between layers [S ~> ppt]
real :: S_ref, T_ref ! Reference salinity [S ~> ppt] and temperature [C ~> degC] within surface layer
real :: S_range, T_range ! Range of salinities [S ~> ppt] and temperatures [C ~> degC] over the vertical
real :: S_surf ! Initial surface salinity [S ~> ppt]
real :: T_bay ! Temperature in the inflow embayment [C ~> degC]
real :: xi0, xi1 ! Fractional vertical positions [nondim]
real :: dome2d_width_bay ! Width of shelf, as fraction of domain [nondim]
real :: dome2d_width_bottom ! Width of deep ocean basin, as fraction of domain [nondim]
Expand All @@ -262,9 +264,14 @@ subroutine DOME2d_initialize_temperature_salinity ( T, S, h, G, GV, US, param_fi
call get_param(param_file, mdl, "T_REF", T_ref, 'Reference temperature', &
units='degC', scale=US%degC_to_C, fail_if_missing=.not.just_read, do_not_log=just_read)
call get_param(param_file, mdl, "S_RANGE", S_range,' Initial salinity range', &
units='1e-3', default=2.0, scale=US%ppt_to_S, do_not_log=just_read)
units='1e-3', default=2.0, scale=US%ppt_to_S, do_not_log=just_read)
call get_param(param_file, mdl, "T_RANGE", T_range, 'Initial temperature range', &
units='degC', default=0.0, scale=US%degC_to_C, do_not_log=just_read)
units='degC', default=0.0, scale=US%degC_to_C, do_not_log=just_read)
call get_param(param_file, mdl, "INITIAL_SSS", S_surf, "Initial surface salinity", &
units="1e-3", default=34.0, scale=US%ppt_to_S, do_not_log=just_read)
call get_param(param_file, mdl, "DOME2D_T_BAY", T_bay, &
"Temperature in the inflow embayment in the DOME2d test case", &
units="degC", default=1.0, scale=US%degC_to_C, do_not_log=just_read)

if (just_read) return ! All run-time parameters have been read, so return.

Expand All @@ -281,7 +288,7 @@ subroutine DOME2d_initialize_temperature_salinity ( T, S, h, G, GV, US, param_fi
xi0 = 0.0
do k = 1,nz
xi1 = xi0 + (GV%H_to_Z * h(i,j,k)) / G%max_depth
S(i,j,k) = 34.0*US%ppt_to_S + 0.5 * S_range * (xi0 + xi1)
S(i,j,k) = S_surf + 0.5 * S_range * (xi0 + xi1)
xi0 = xi1
enddo
enddo ; enddo
Expand All @@ -292,12 +299,12 @@ subroutine DOME2d_initialize_temperature_salinity ( T, S, h, G, GV, US, param_fi
xi0 = 0.0
do k = 1,nz
xi1 = xi0 + (GV%H_to_Z * h(i,j,k)) / G%max_depth
S(i,j,k) = 34.0*US%ppt_to_S + 0.5 * S_range * (xi0 + xi1)
S(i,j,k) = S_surf + 0.5 * S_range * (xi0 + xi1)
xi0 = xi1
enddo
x = ( G%geoLonT(i,j) - G%west_lon ) / G%len_lon
if ( x <= dome2d_width_bay ) then
S(i,j,nz) = 34.0*US%ppt_to_S + S_range
S(i,j,nz) = S_surf + S_range
endif
enddo ; enddo

Expand All @@ -322,7 +329,7 @@ subroutine DOME2d_initialize_temperature_salinity ( T, S, h, G, GV, US, param_fi
x = ( G%geoLonT(i,j) - G%west_lon ) / G%len_lon
if ( x <= dome2d_width_bay ) then
S(i,j,1:index_bay_z) = S_ref + S_range ! Use for z coordinates
T(i,j,1:index_bay_z) = 1.0*US%degC_to_C ! Use for z coordinates
T(i,j,1:index_bay_z) = T_bay ! Use for z coordinates
endif
enddo ; enddo ! i and j loops
endif ! Z initial conditions
Expand All @@ -332,8 +339,8 @@ subroutine DOME2d_initialize_temperature_salinity ( T, S, h, G, GV, US, param_fi
do i = G%isc,G%iec ; do j = G%jsc,G%jec
x = ( G%geoLonT(i,j) - G%west_lon ) / G%len_lon
if ( x <= dome2d_width_bay ) then
S(i,j,1:GV%ke) = S_ref + S_range ! Use for sigma coordinates
T(i,j,1:GV%ke) = 1.0*US%degC_to_C ! Use for sigma coordinates
S(i,j,1:GV%ke) = S_ref + S_range ! Use for sigma coordinates
T(i,j,1:GV%ke) = T_bay ! Use for sigma coordinates
endif
enddo ; enddo
endif
Expand All @@ -344,7 +351,7 @@ subroutine DOME2d_initialize_temperature_salinity ( T, S, h, G, GV, US, param_fi
do i = G%isc,G%iec ; do j = G%jsc,G%jec
x = ( G%geoLonT(i,j) - G%west_lon ) / G%len_lon
if ( x <= dome2d_width_bay ) then
T(i,j,GV%ke) = 1.0*US%degC_to_C
T(i,j,GV%ke) = T_bay
endif
enddo ; enddo
endif
Expand Down Expand Up @@ -373,6 +380,8 @@ subroutine DOME2d_initialize_sponges(G, GV, US, tv, depth_tot, param_file, use_A
real :: T_ref ! Reference temperature within the surface layer [C ~> degC]
real :: S_range ! Range of salinities in the vertical [S ~> ppt]
real :: T_range ! Range of temperatures in the vertical [C ~> degC]
real :: S_range_sponge ! Range of salinities in the vertical in the east sponge [S ~> ppt]
real :: S_surf ! Initial surface salinity [S ~> ppt]
real :: e0(SZK_(GV)+1) ! The resting interface heights [Z ~> m],
! usually negative because it is positive upward.
real :: eta1D(SZK_(GV)+1) ! Interface height relative to the sea surface
Expand Down Expand Up @@ -428,7 +437,11 @@ subroutine DOME2d_initialize_sponges(G, GV, US, tv, depth_tot, param_file, use_A
call get_param(param_file, mdl, "T_REF", T_ref, units="degC", scale=US%degC_to_C, fail_if_missing=.false.)
call get_param(param_file, mdl, "S_RANGE", S_range, units="ppt", default=2.0, scale=US%ppt_to_S)
call get_param(param_file, mdl, "T_RANGE", T_range, units="degC", default=0.0, scale=US%degC_to_C)

call get_param(param_file, mdl, "INITIAL_SSS", S_surf, "Initial surface salinity", &
units="1e-3", default=34.0, scale=US%ppt_to_S, do_not_log=.true.)
call get_param(param_file, mdl, "DOME2D_EAST_SPONGE_S_RANGE", S_range_sponge, &
"Range of salinities in the eastern sponge region in the DOME2D configuration", &
units="1e-3", default=1.0, scale=US%ppt_to_S)

! Set the sponge damping rate as a function of position
Idamp(:,:) = 0.0
Expand All @@ -454,7 +467,7 @@ subroutine DOME2d_initialize_sponges(G, GV, US, tv, depth_tot, param_file, use_A

if (use_ALE) then

! Construct a grid (somewhat arbitrarily) to describe the sponge T/S on
! Construct a grid (somewhat arbitrarily) to describe the sponge T/S on
do k=1,nz
e0(k) = -G%max_depth * ( real(k-1) / real(nz) )
enddo
Expand All @@ -480,7 +493,9 @@ subroutine DOME2d_initialize_sponges(G, GV, US, tv, depth_tot, param_file, use_A
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*US%ppt_to_S - 1.0*US%ppt_to_S * (z / (G%max_depth))
! 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
Expand Down
17 changes: 13 additions & 4 deletions src/user/DOME_initialization.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -306,6 +306,8 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, US, PF, tr_Reg)
real :: drho_dT(SZK_(GV)) ! Derivative of density with temperature [R C-1 ~> kg m-3 degC-1].
real :: drho_dS(SZK_(GV)) ! Derivative of density with salinity [R S-1 ~> kg m-3 ppt-1].
real :: rho_guess(SZK_(GV)) ! Potential density at T0 & S0 [R ~> kg m-3].
real :: S_ref ! A default value for salinities [S ~> ppt]
real :: T_light ! A first guess at the temperature of the lightest layer [C ~> degC]
! The following variables are used to set up the transport in the DOME example.
real :: tr_0 ! The total integrated inflow transport [H L2 T-1 ~> m3 s-1 or kg s-1]
real :: tr_k ! The integrated inflow transport of a layer [H L2 T-1 ~> m3 s-1 or kg s-1]
Expand Down Expand Up @@ -357,6 +359,13 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, US, PF, tr_Reg)
"inflow properties.", units="s-1", default=f_0*US%s_to_T, scale=US%T_to_s)
call get_param(PF, mdl, "DOME_INFLOW_LON", inflow_lon, &
"The edge longitude of the DOME inflow.", units="km", default=1000.0)
if (associated(tv%S) .or. associated(tv%T)) then
call get_param(PF, mdl, "S_REF", S_ref, &
units="ppt", default=35.0, scale=US%ppt_to_S, do_not_log=.true.)
call get_param(PF, mdl, "DOME_T_LIGHT", T_light, &
"A first guess at the temperature of the lightest layer in the DOME test case.", &
units="degC", default=25.0, scale=US%degC_to_C)
endif

if (.not.associated(OBC)) return

Expand Down Expand Up @@ -413,16 +422,16 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, US, PF, tr_Reg)
! The inflow values of temperature and salinity also need to be set here if
! these variables are used. The following code is just a naive example.
if (associated(tv%S)) then
! In this example, all S inflows have values of 35 psu.
! In this example, all S inflows have values given by S_ref.
name = 'salt'
call tracer_name_lookup(tr_Reg, tr_ptr, name)
call register_segment_tracer(tr_ptr, PF, GV, segment, OBC_scalar=35.0*US%ppt_to_S, scale=US%ppt_to_S)
call register_segment_tracer(tr_ptr, PF, GV, segment, OBC_scalar=S_ref, scale=US%ppt_to_S)
endif
if (associated(tv%T)) then
! In this example, the T values are set to be consistent with the layer
! target density and a salinity of 35 psu. This code is taken from
! target density and a salinity of S_ref. This code is taken from
! USER_initialize_temp_sal.
pres(:) = tv%P_Ref ; S0(:) = 35.0*US%ppt_to_S ; T0(1) = 25.0*US%degC_to_C
pres(:) = tv%P_Ref ; S0(:) = S_ref ; T0(1) = T_light
call calculate_density(T0(1), S0(1), pres(1), rho_guess(1), tv%eqn_of_state)
call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, tv%eqn_of_state, (/1,1/) )

Expand Down
30 changes: 23 additions & 7 deletions src/user/benchmark_initialization.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -101,9 +101,11 @@ subroutine benchmark_initialize_thickness(h, depth_tot, G, GV, US, param_file, e
! in depth units [Z ~> m].
real :: eta1D(SZK_(GV)+1) ! Interface height relative to the sea surface
! positive upward, in depth units [Z ~> m].
real :: SST ! The initial sea surface temperature [C ~> degC].
real :: T_int ! The initial temperature of an interface [C ~> degC].
real :: ML_depth ! The specified initial mixed layer depth, in depth units [Z ~> m].
real :: SST ! The initial sea surface temperature [C ~> degC].
real :: S_ref ! A default value for salinities [S ~> ppt]
real :: T_light ! A first guess at the temperature of the lightest layer [C ~> degC]
real :: T_int ! The initial temperature of an interface [C ~> degC].
real :: ML_depth ! The specified initial mixed layer depth, in depth units [Z ~> m].
real :: thermocline_scale ! The e-folding scale of the thermocline, in depth units [Z ~> m].
real, dimension(SZK_(GV)) :: &
T0, S0, & ! Profiles of temperature [C ~> degC] and salinity [S ~> ppt]
Expand Down Expand Up @@ -135,6 +137,12 @@ subroutine benchmark_initialize_thickness(h, depth_tot, G, GV, US, param_file, e
call get_param(param_file, mdl, "BENCHMARK_THERMOCLINE_SCALE", thermocline_scale, &
"Initial thermocline depth scale in the benchmark test case.", &
default=500.0, units="m", scale=US%m_to_Z, do_not_log=just_read)
call get_param(param_file, mdl, "BENCHMARK_T_LIGHT", T_light, &
"A first guess at the temperature of the lightest layer in the benchmark test case.", &
units="degC", default=29.0, scale=US%degC_to_C, do_not_log=just_read)
call get_param(param_file, mdl, "S_REF", S_ref, &
"The uniform salinities used to initialize the benchmark test case.", &
units="ppt", default=35.0, scale=US%ppt_to_S, do_not_log=just_read)

if (just_read) return ! This subroutine has no run-time parameters.

Expand All @@ -147,9 +155,9 @@ subroutine benchmark_initialize_thickness(h, depth_tot, G, GV, US, param_file, e
! This block calculates T0(k) for the purpose of diagnosing where the
! interfaces will be found.
do k=1,nz
pres(k) = P_Ref ; S0(k) = 35.0*US%ppt_to_S
pres(k) = P_Ref ; S0(k) = S_ref
enddo
T0(k1) = 29.0*US%degC_to_C
T0(k1) = T_light
call calculate_density(T0(k1), S0(k1), pres(k1), rho_guess(k1), eqn_of_state)
call calculate_density_derivs(T0(k1), S0(k1), pres(k1), drho_dT(k1), drho_dS(k1), eqn_of_state)

Expand Down Expand Up @@ -232,25 +240,33 @@ subroutine benchmark_init_temperature_salinity(T, S, G, GV, US, param_file, &
! Local variables
real :: T0(SZK_(GV)) ! A profile of temperatures [C ~> degC]
real :: S0(SZK_(GV)) ! A profile of salinities [S ~> ppt]
real :: S_ref ! A default value for salinities [S ~> ppt]
real :: T_light ! A first guess at the temperature of the lightest layer [C ~> degC]
real :: pres(SZK_(GV)) ! Reference pressure [R L2 T-2 ~> Pa]
real :: drho_dT(SZK_(GV)) ! Derivative of density with temperature [R C-1 ~> kg m-3 degC-1]
real :: drho_dS(SZK_(GV)) ! Derivative of density with salinity [R S-1 ~> kg m-3 ppt-1]
real :: rho_guess(SZK_(GV)) ! Potential density at T0 & S0 [R ~> kg m-3]
real :: PI ! 3.1415926... calculated as 4*atan(1)
real :: SST ! The initial sea surface temperature [C ~> degC]
character(len=40) :: mdl = "benchmark_init_temperature_salinity" ! This subroutine's name.
integer :: i, j, k, k1, is, ie, js, je, nz, itt

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

call get_param(param_file, mdl, "S_REF", S_ref, &
units="ppt", default=35.0, scale=US%ppt_to_S, do_not_log=.true.)
call get_param(param_file, mdl, "BENCHMARK_T_LIGHT", T_light, &
units="degC", default=29.0, scale=US%degC_to_C, do_not_log=.true.)

if (just_read) return ! All run-time parameters have been read, so return.

k1 = GV%nk_rho_varies + 1

do k=1,nz
pres(k) = P_Ref ; S0(k) = 35.0*US%ppt_to_S
pres(k) = P_Ref ; S0(k) = S_ref
enddo

T0(k1) = 29.0*US%degC_to_C
T0(k1) = T_light
call calculate_density(T0(k1), S0(k1), pres(k1), rho_guess(k1), eqn_of_state)
call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, eqn_of_state, (/k1,k1/) )

Expand Down
Loading