sloshing_initialization.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Initialization for the "sloshing" internal waves configuration.
7
8use mom_domains, only : sum_across_pes
10use mom_error_handler, only : mom_mesg, mom_error, fatal, is_root_pe
11use mom_file_parser, only : get_param, log_version, param_file_type
12use mom_get_input, only : directories
13use mom_grid, only : ocean_grid_type
15use mom_tracer_registry, only : tracer_registry_type
19
20implicit none ; private
21
22#include <MOM_memory.h>
23
24! The following routines are visible to the outside world
28
29contains
30
31!> Initialization of topography.
32subroutine sloshing_initialize_topography( D, G, param_file, max_depth )
33 type(dyn_horgrid_type), intent(in) :: g !< The dynamic horizontal grid type
34 real, dimension(G%isd:G%ied,G%jsd:G%jed), &
35 intent(out) :: d !< Ocean bottom depth [Z ~> m]
36 type(param_file_type), intent(in) :: param_file !< Parameter file structure
37 real, intent(in) :: max_depth !< Maximum ocean depth [Z ~> m]
38
39 ! Local variables
40 integer :: i, j
41
42 do i=g%isc,g%iec ; do j=g%jsc,g%jec
43 d(i,j) = max_depth
44 enddo ; enddo
45
47
48
49!> Initialization of thicknesses
50!! This routine is called when THICKNESS_CONFIG is set to 'sloshing'
51!!
52!! This routine initializes layer positions to set off a sloshing motion in
53!! the zonal direction in a rectangular basin. All layers have initially the
54!! same thickness but all interfaces (except bottom and sea surface) are
55!! displaced according to a half-period cosine, with maximum value on the
56!! left and minimum value on the right. This sets off a regular sloshing motion.
57subroutine sloshing_initialize_thickness ( h, depth_tot, G, GV, US, param_file, just_read)
58 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
59 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
60 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
61 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
62 intent(out) :: h !< The thickness that is being initialized [Z ~> m]
63 real, dimension(SZI_(G),SZJ_(G)), &
64 intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
65 type(param_file_type), intent(in) :: param_file !< A structure to parse for model parameter values.
66 logical, intent(in) :: just_read !< If true, this call will only read
67 !! parameters without changing h.
68
69 ! Local variables
70 real :: displ(szk_(gv)+1) ! The interface displacement [Z ~> m].
71 real :: z_unif(szk_(gv)+1) ! Fractional uniform interface heights [nondim].
72 real :: z_inter(szk_(gv)+1) ! Interface heights [Z ~> m]
73 real :: a0 ! The displacement amplitude [Z ~> m].
74 real :: weight_z ! A depth-space weighting [nondim].
75 real :: x1, y1, x2, y2 ! Dimensonless parameters specifying the depth profile [nondim]
76 real :: x, t ! Dimensionless depth coordinates scales [nondim]
77 logical :: use_ic_bug ! If true, set the initial conditions retaining an old bug.
78 ! This include declares and sets the variable "version".
79# include "version_variable.h"
80 character(len=40) :: mdl = "sloshing_initialization" !< This module's name.
81 integer :: i, j, k, is, ie, js, je, nz
82
83 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
84
85 if (.not.just_read) call log_version(param_file, mdl, version, "")
86 call get_param(param_file, mdl, "SLOSHING_IC_AMPLITUDE", a0, &
87 "Initial amplitude of sloshing internal interface height "//&
88 "displacements it the sloshing test case.", &
89 units='m', default=75.0, scale=us%m_to_Z, do_not_log=just_read)
90 call get_param(param_file, mdl, "SLOSHING_IC_BUG", use_ic_bug, &
91 "If true, use code with a bug to set the sloshing initial conditions.", &
92 default=.false., do_not_log=just_read)
93
94 if (just_read) return ! All run-time parameters have been read, so return.
95
96 ! Define thicknesses
97 do j=g%jsc,g%jec ; do i=g%isc,g%iec
98
99 ! Define uniform interfaces
100 do k = 0,nz
101 z_unif(k+1) = -real(k)/real(nz)
102 enddo
103
104 ! 1. Define stratification
105 do k = 1,nz+1
106
107 ! Thin pycnocline in the middle
108 !z_inter(k) = (2.0**(n-1)) * (z_unif(k) + 0.5)**n - 0.5
109
110 ! Thin pycnocline in the middle (piecewise linear profile)
111 x1 = 0.30 ; y1 = 0.48 ; x2 = 0.70 ; y2 = 0.52
112
113 x = -z_unif(k)
114
115 if ( x <= x1 ) then
116 t = y1*x/x1
117 elseif ( (x > x1 ) .and. ( x < x2 )) then
118 t = y1 + (y2-y1) * (x-x1) / (x2-x1)
119 else
120 t = y2 + (1.0-y2) * (x-x2) / (1.0-x2)
121 endif
122
123 t = - z_unif(k)
124
125 z_inter(k) = -t * g%max_depth
126
127 enddo
128
129 ! 2. Define displacement
130 ! a0 is set via get_param; by default a0 is a 75m Displacement amplitude in depth units.
131 do k = 1,nz+1
132
133 weight_z = - 4.0 * ( z_unif(k) + 0.5 )**2 + 1.0
134
135 x = g%geoLonT(i,j) / g%len_lon
136 if (use_ic_bug) then
137 displ(k) = a0 * cos(acos(-1.0)*x) + weight_z * us%m_to_Z ! There is a flag to fix this bug.
138 else
139 displ(k) = a0 * cos(acos(-1.0)*x) * weight_z
140 endif
141
142 if ( k == 1 ) then
143 displ(k) = 0.0
144 endif
145
146 if ( k == nz+1 ) then
147 displ(k) = 0.0
148 endif
149
150 z_inter(k) = z_inter(k) + displ(k)
151
152 enddo
153
154 ! 3. The last interface must coincide with the seabed
155 z_inter(nz+1) = -depth_tot(i,j)
156 ! Modify interface heights to make sure all thicknesses are strictly positive
157 do k = nz,1,-1
158 if ( z_inter(k) < (z_inter(k+1) + gv%Angstrom_Z) ) then
159 z_inter(k) = z_inter(k+1) + gv%Angstrom_Z
160 endif
161 enddo
162
163 ! 4. Define layers
164 do k = 1,nz
165 h(i,j,k) = z_inter(k) - z_inter(k+1)
166 enddo
167
168 enddo ; enddo
169
171
172
173!> Initialization of temperature and salinity
174!!
175!! This subroutine initializes linear profiles for T and S according to
176!! reference surface layer salinity and temperature and a specified range.
177!! Note that the linear distribution is set up with respect to the layer
178!! number, not the physical position).
179subroutine sloshing_initialize_temperature_salinity ( T, S, h, G, GV, US, param_file, just_read)
180 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure.
181 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
182 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: t !< Potential temperature [C ~> degC].
183 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: s !< Salinity [S ~> ppt].
184 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [Z ~> m].
185 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
186 type(param_file_type), intent(in) :: param_file !< A structure to parse
187 !! for model parameter values.
188 logical, intent(in) :: just_read !< If true, this call will only read
189 !! parameters without changing T & S.
190
191 ! Local variables
192 real :: delta_t ! Temperature difference between layers [C ~> degC]
193 real :: s_ref, t_ref ! Reference salinity [S ~> ppt] and temperature [C ~> degC] within surface layer
194 real :: s_range, t_range ! Range of salinities [S ~> ppt] and temperatures [C ~> degC] over the vertical
195 real :: s_surf ! Initial surface salinity [S ~> ppt]
196 real :: t_pert ! A perturbed temperature [C ~> degC]
197 integer :: kdelta ! Half the number of layers with the temperature perturbation
198 real :: deltah ! Thickness of each layer [Z ~> m]
199 real :: xi0, xi1 ! Fractional vertical positions [nondim]
200 character(len=40) :: mdl = "sloshing_initialization" ! This module's name.
201 integer :: i, j, k, is, ie, js, je, nz
202
203 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
204
205 call get_param(param_file, mdl, "S_REF", s_ref, 'Reference value for salinity', &
206 default=35.0, units="ppt", scale=us%ppt_to_S, do_not_log=just_read)
207 call get_param(param_file, mdl, "T_REF", t_ref, 'Reference value for temperature', &
208 units='degC', scale=us%degC_to_C, fail_if_missing=.not.just_read, do_not_log=just_read)
209
210 ! The default is to assume an increase by 2 ppt for the salinity and a uniform temperature.
211 call get_param(param_file, mdl, "S_RANGE", s_range, 'Initial salinity range.', &
212 units="ppt", default=2.0, scale=us%ppt_to_S, do_not_log=just_read)
213 call get_param(param_file, mdl, "T_RANGE", t_range, 'Initial temperature range', &
214 units='degC', default=0.0, scale=us%degC_to_C, do_not_log=just_read)
215 call get_param(param_file, mdl, "INITIAL_SSS", s_surf, "Initial surface salinity", &
216 units="ppt", default=34.0, scale=us%ppt_to_S, do_not_log=just_read)
217 call get_param(param_file, mdl, "SLOSHING_T_PERT", t_pert, &
218 'A mid-column temperature perturbation in the sloshing test case', &
219 units='degC', default=1.0, scale=us%degC_to_C, do_not_log=just_read)
220
221 if (just_read) return ! All run-time parameters have been read, so return.
222
223 ! Prescribe salinity
224 !delta_S = S_range / ( GV%ke - 1.0 )
225
226 !S(:,:,1) = S_ref
227 !do k = 2,GV%ke
228 ! S(:,:,k) = S(:,:,k-1) + delta_S
229 !enddo
230
231 deltah = g%max_depth / nz
232 do j=js,je ; do i=is,ie
233 xi0 = 0.0
234 do k = 1,nz
235 xi1 = xi0 + deltah / g%max_depth ! = xi0 + 1.0 / real(nz)
236 s(i,j,k) = s_surf + 0.5 * s_range * (xi0 + xi1)
237 xi0 = xi1
238 enddo
239 enddo ; enddo
240
241 ! Prescribe temperature
242 delta_t = t_range / ( gv%ke - 1.0 )
243
244 t(:,:,1) = t_ref
245 do k = 2,gv%ke
246 t(:,:,k) = t(:,:,k-1) + delta_t
247 enddo
248 kdelta = 2
249 ! Perhaps the following lines should instead assign T() = T_pert + T_ref
250 t(:,:,gv%ke/2 - (kdelta-1):gv%ke/2 + kdelta) = t_pert
251
253
254!> \namespace sloshing_initialization
255!!
256!! The module configures the model for the non-rotating sloshing test case.