adjustment_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!> Configures the model for the geostrophic adjustment test case.
7
8use mom_error_handler, only : mom_mesg, mom_error, fatal, is_root_pe
9use mom_file_parser, only : get_param, log_param, log_version, param_file_type
10use mom_get_input, only : directories
11use mom_grid, only : ocean_grid_type
15use regrid_consts, only : coordinatemode, default_coordinate_mode
16use regrid_consts, only : regridding_layer, regridding_zstar
17use regrid_consts, only : regridding_rho, regridding_sigma
18
19implicit none ; private
20
21character(len=40) :: mdl = "adjustment_initialization" !< This module's name.
22
23#include <MOM_memory.h>
24
27
28! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
29! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
30! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
31! vary with the Boussinesq approximation, the Boussinesq variant is given first.
32
33contains
34
35!> Initializes the layer thicknesses in the adjustment test case
36subroutine adjustment_initialize_thickness ( h, G, GV, US, param_file, just_read)
37 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
38 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
39 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
40 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
41 intent(out) :: h !< The thickness that is being initialized [Z ~> m]
42 type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
43 !! to parse for model parameter values.
44 logical, intent(in) :: just_read !< If true, this call will only read
45 !! parameters without changing h.
46 ! Local variables
47 real :: e0(szk_(gv)+1) ! The resting interface heights, in depth units [Z ~> m], usually
48 ! negative because it is positive upward.
49 real :: eta1d(szk_(gv)+1) ! Interface height relative to the sea surface
50 ! positive upward, in depth units [Z ~> m].
51 real :: drho_ds ! The partial derivative of density with salinity [R S-1 ~> kg m-3 ppt-1].
52 ! In this subroutine it is hard coded at 1.0 kg m-3 ppt-1.
53 real :: x, y, yy ! Fractional positions in the x- and y-directions [nondim]
54 real :: y_lat ! y-positions in the units of latitude [m] or [km] or [degrees]
55 real :: s_ref ! Reference salinity within surface layer [S ~> ppt]
56 real :: s_range ! Range of salinities in the vertical [S ~> ppt]
57 real :: dsdz ! Vertical salinity gradient [S Z-1 ~> ppt m-1]
58 real :: delta_s ! The local salinity perturbation [S ~> ppt]
59 real :: delta_s_strat ! Top-to-bottom salinity difference of stratification [S ~> ppt]
60 real :: min_thickness ! The minimum layer thickness [Z ~> m]
61 real :: adjustment_delta ! Interface height anomalies, positive downward [Z ~> m]
62 real :: adjustment_width ! Width of the frontal zone [m] or [km] or [degrees]
63 real :: adjustment_deltas ! Salinity difference across front [S ~> ppt]
64 real :: front_wave_amp ! Amplitude of trans-frontal wave perturbation [m] or [km] or [degrees]
65 real :: front_wave_length ! Wave-length of trans-frontal wave perturbation [m] or [km] or [degrees]
66 real :: front_wave_asym ! Amplitude of frontal asymmetric perturbation [m] or [km] or [degrees]
67 real :: target_values(szk_(gv)+1) ! Target densities or density anomalies [R ~> kg m-3]
68 character(len=20) :: verticalcoordinate
69 ! This include declares and sets the variable "version".
70# include "version_variable.h"
71 integer :: i, j, k, is, ie, js, je, nz
72
73 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
74
75 if (.not.just_read) &
76 call mom_mesg("adjustment_initialize_thickness: setting thickness")
77
78 ! Parameters used by main model initialization
79 if (.not.just_read) call log_version(param_file, mdl, version, "")
80 call get_param(param_file, mdl, "S_REF", s_ref, 'Reference salinity', &
81 default=35.0, units='ppt', scale=us%ppt_to_S, do_not_log=just_read)
82 call get_param(param_file, mdl, "MIN_THICKNESS", min_thickness, 'Minimum layer thickness', &
83 default=1.0e-3, units='m', scale=us%m_to_Z, do_not_log=just_read)
84 call get_param(param_file, mdl, "DRHO_DS", drho_ds, &
85 "The partial derivative of density with salinity with a linear equation of state.", &
86 units="kg m-3 ppt-1", default=0.8, scale=us%kg_m3_to_R*us%S_to_ppt)
87
88 ! Parameters specific to this experiment configuration
89 call get_param(param_file, mdl, "REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
90 default=default_coordinate_mode, do_not_log=just_read)
91 call get_param(param_file, mdl, "ADJUSTMENT_WIDTH", adjustment_width, &
92 "Width of frontal zone", &
93 units=g%x_ax_unit_short, fail_if_missing=.not.just_read, do_not_log=just_read)
94 call get_param(param_file, mdl, "DELTA_S_STRAT", delta_s_strat, &
95 "Top-to-bottom salinity difference of stratification", &
96 units="ppt", scale=us%ppt_to_S, fail_if_missing=.not.just_read, do_not_log=just_read)
97 call get_param(param_file, mdl, "ADJUSTMENT_DELTAS", adjustment_deltas, &
98 "Salinity difference across front", &
99 units="ppt", scale=us%ppt_to_S, fail_if_missing=.not.just_read, do_not_log=just_read)
100 call get_param(param_file, mdl, "FRONT_WAVE_AMP", front_wave_amp, &
101 "Amplitude of trans-frontal wave perturbation", &
102 units=g%x_ax_unit_short, default=0., do_not_log=just_read)
103 call get_param(param_file, mdl, "FRONT_WAVE_LENGTH", front_wave_length, &
104 "Wave-length of trans-frontal wave perturbation", &
105 units=g%x_ax_unit_short, default=0., do_not_log=just_read)
106 call get_param(param_file, mdl, "FRONT_WAVE_ASYM", front_wave_asym, &
107 "Amplitude of frontal asymmetric perturbation", &
108 units=g%x_ax_unit_short, default=0., do_not_log=just_read)
109
110 if (just_read) return ! All run-time parameters have been read, so return.
111
112 ! WARNING: this routine specifies the interface heights so that the last layer
113 ! is vanished, even at maximum depth. In order to have a uniform
114 ! layer distribution, use this line of code within the loop:
115 ! e0(k) = -G%max_depth * real(k-1) / real(nz)
116 ! To obtain a thickness distribution where the last layer is
117 ! vanished and the other thicknesses uniformly distributed, use:
118 ! e0(k) = -G%max_depth * real(k-1) / real(nz-1)
119
120 dsdz = -delta_s_strat / g%max_depth
121
122 select case ( coordinatemode(verticalcoordinate) )
123
124 case ( regridding_layer, regridding_rho )
125 if (delta_s_strat /= 0.) then
126 ! This was previously coded ambiguously.
127 adjustment_delta = (adjustment_deltas / delta_s_strat) * g%max_depth
128 do k=1,nz+1
129 e0(k) = adjustment_delta - (g%max_depth + 2*adjustment_delta) * (real(k-1) / real(nz))
130 enddo
131 else
132 adjustment_delta = 2.*g%max_depth
133 do k=1,nz+1
134 e0(k) = -g%max_depth * (real(k-1) / real(nz))
135 enddo
136 endif
137 if (nz > 1) then
138 target_values(1) = ( gv%Rlay(1) + 0.5*(gv%Rlay(1)-gv%Rlay(2)) )
139 target_values(nz+1) = ( gv%Rlay(nz) + 0.5*(gv%Rlay(nz)-gv%Rlay(nz-1)) )
140 else ! This might not be needed, but it avoids segmentation faults if nz=1.
141 target_values(1) = 0.0
142 target_values(nz+1) = 2.0 * gv%Rlay(1)
143 endif
144 do k = 2,nz
145 target_values(k) = target_values(k-1) + ( gv%Rlay(nz) - gv%Rlay(1) ) / (nz-1)
146 enddo
147 target_values(:) = target_values(:) - 1000.0*us%kg_m3_to_R
148 do j=js,je ; do i=is,ie
149 if (front_wave_length /= 0.) then
150 y = ( 0.125 + g%geoLatT(i,j) / front_wave_length ) * ( 4. * acos(0.) )
151 yy = 2. * ( g%geoLatT(i,j) - 0.5 * g%len_lat ) / adjustment_width
152 yy = min(1.0, yy) ; yy = max(-1.0, yy)
153 yy = yy * 2. * acos( 0. )
154 y_lat = front_wave_amp*sin(y) + front_wave_asym*sin(yy)
155 else
156 y_lat = 0.
157 endif
158 x = ( ( g%geoLonT(i,j) - 0.5 * g%len_lon ) + y_lat ) / adjustment_width
159 x = min(1.0, x) ; x = max(-1.0, x)
160 x = x * acos( 0. )
161 delta_s = adjustment_deltas * 0.5 * (1. - sin( x ) )
162 do k=2,nz
163 if (drho_ds*dsdz /= 0.) then
164 eta1d(k) = ( target_values(k) - drho_ds*( s_ref + delta_s ) ) / (drho_ds*dsdz)
165 else
166 eta1d(k) = e0(k) - (0.5*adjustment_delta) * sin( x )
167 endif
168 eta1d(k) = max( eta1d(k), -g%max_depth )
169 eta1d(k) = min( eta1d(k), 0. )
170 enddo
171 eta1d(1) = 0. ; eta1d(nz+1) = -g%max_depth
172 do k=nz,1,-1
173 if (eta1d(k) > 0.) then
174 eta1d(k) = max( eta1d(k+1) + min_thickness, 0. )
175 h(i,j,k) = max( eta1d(k) - eta1d(k+1), min_thickness )
176 elseif (eta1d(k) <= (eta1d(k+1) + min_thickness)) then
177 eta1d(k) = eta1d(k+1) + min_thickness
178 h(i,j,k) = min_thickness
179 else
180 h(i,j,k) = eta1d(k) - eta1d(k+1)
181 endif
182 enddo
183 enddo ; enddo
184
185 case ( regridding_zstar, regridding_sigma )
186 do k=1,nz+1
187 eta1d(k) = -g%max_depth * (real(k-1) / real(nz))
188 eta1d(k) = max(min(eta1d(k), 0.), -g%max_depth)
189 enddo
190 do j=js,je ; do i=is,ie
191 do k=nz,1,-1
192 h(i,j,k) = eta1d(k) - eta1d(k+1)
193 enddo
194 enddo ; enddo
195
196 case default
197 call mom_error(fatal, "adjustment_initialize_thickness: "// &
198 "Unrecognized i.c. setup - set ADJUSTMENT_IC")
199
200 end select
201
203
204!> Initialization of temperature and salinity in the adjustment test case
205subroutine adjustment_initialize_temperature_salinity(T, S, h, depth_tot, G, GV, US, param_file, just_read)
206 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
207 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
208 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
209 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
210 intent(out) :: t !< The temperature that is being initialized [C ~> degC]
211 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
212 intent(out) :: s !< The salinity that is being initialized [S ~> ppt]
213 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
214 intent(in) :: h !< The model thicknesses [Z ~> m]
215 real, dimension(SZI_(G),SZJ_(G)), &
216 intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
217 type(param_file_type), intent(in) :: param_file !< A structure indicating the open file to
218 !! parse for model parameter values.
219 logical, intent(in) :: just_read !< If true, this call will only read
220 !! parameters without changing T & S.
221
222 real :: x, y, yy ! Fractional positions in the x- and y-directions [nondim]
223 real :: y_lat ! y-position in the units of latitude [m] or [km] or [degrees]
224 real :: s_ref ! Reference salinity within surface layer [S ~> ppt]
225 real :: t_ref ! Reference temperature within surface layer [C ~> degC]
226 real :: s_range ! Range of salinities in the vertical [S ~> ppt]
227 real :: t_range ! Range of temperatures in the vertical [C ~> degC]
228 real :: dsdz ! Vertical salinity gradient [S Z-1 ~> ppt m-1]
229 real :: delta_s ! The local salinity perturbation [S ~> ppt]
230 real :: delta_s_strat ! Top-to-bottom salinity difference of stratification [S ~> ppt]
231 real :: adjustment_width ! Width of the frontal zone [m] or [km] or [degrees]
232 real :: adjustment_deltas ! Salinity difference across front [S ~> ppt]
233 real :: front_wave_amp ! Amplitude of trans-frontal wave perturbation [m] or [km] or [degrees]
234 real :: front_wave_length ! Wave-length of trans-frontal wave perturbation [m] or [km] or [degrees]
235 real :: front_wave_asym ! Amplitude of frontal asymmetric perturbation [m] or [km] or [degrees]
236 real :: eta1d(szk_(gv)+1) ! Interface heights [Z ~> m]
237 character(len=20) :: verticalcoordinate
238 integer :: i, j, k, is, ie, js, je, nz
239
240 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
241
242 ! Parameters used by main model initialization
243 call get_param(param_file, mdl, "S_REF", s_ref, 'Reference salinity', &
244 default=35.0, units="ppt", scale=us%ppt_to_S, do_not_log=just_read)
245 call get_param(param_file, mdl, "T_REF", t_ref, 'Reference temperature', &
246 units="degC", scale=us%degC_to_C, fail_if_missing=.not.just_read, do_not_log=just_read)
247 call get_param(param_file, mdl, "S_RANGE", s_range, 'Initial salinity range', &
248 default=2.0, units="ppt", scale=us%ppt_to_S, do_not_log=just_read)
249 call get_param(param_file, mdl, "T_RANGE", t_range, 'Initial temperature range', &
250 default=1.0, units='degC', scale=us%degC_to_C, do_not_log=just_read)
251 ! Parameters specific to this experiment configuration BUT logged in previous s/r
252 call get_param(param_file, mdl, "REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
253 default=default_coordinate_mode, do_not_log=just_read)
254 call get_param(param_file, mdl, "ADJUSTMENT_WIDTH", adjustment_width, &
255 units=g%x_ax_unit_short, fail_if_missing=.not.just_read, do_not_log=.true.)
256 call get_param(param_file, mdl, "ADJUSTMENT_DELTAS", adjustment_deltas, &
257 units="ppt", scale=us%ppt_to_S, fail_if_missing=.not.just_read, do_not_log=.true.)
258 call get_param(param_file, mdl, "DELTA_S_STRAT", delta_s_strat, &
259 units="ppt", scale=us%ppt_to_S, fail_if_missing=.not.just_read, do_not_log=.true.)
260 call get_param(param_file, mdl, "FRONT_WAVE_AMP", front_wave_amp, &
261 units=g%x_ax_unit_short, default=0., do_not_log=.true.)
262 call get_param(param_file, mdl, "FRONT_WAVE_LENGTH", front_wave_length, &
263 units=g%x_ax_unit_short, default=0., do_not_log=.true.)
264 call get_param(param_file, mdl, "FRONT_WAVE_ASYM", front_wave_asym, &
265 units=g%x_ax_unit_short, default=0., do_not_log=.true.)
266
267 if (just_read) return ! All run-time parameters have been read, so return.
268
269 t(:,:,:) = 0.0
270 s(:,:,:) = 0.0
271
272 ! Linear salinity profile
273 select case ( coordinatemode(verticalcoordinate) )
274
275 case ( regridding_zstar, regridding_sigma )
276 dsdz = -delta_s_strat / g%max_depth
277 do j=js,je ; do i=is,ie
278 eta1d(nz+1) = -depth_tot(i,j)
279 do k=nz,1,-1
280 eta1d(k) = eta1d(k+1) + h(i,j,k)
281 enddo
282 if (front_wave_length /= 0.) then
283 y = ( 0.125 + g%geoLatT(i,j) / front_wave_length ) * ( 4. * acos(0.) )
284 yy = 2. * ( g%geoLatT(i,j) - 0.5 * g%len_lat ) / front_wave_length
285 yy = min(1.0, yy) ; yy = max(-1.0, yy)
286 yy = yy * 2. * acos( 0. )
287 y_lat = front_wave_amp*sin(y) + front_wave_asym*sin(yy)
288 else
289 y_lat = 0.
290 endif
291 x = ( ( g%geoLonT(i,j) - 0.5 * g%len_lon ) + y_lat ) / adjustment_width
292 x = min(1.0, x) ; x = max(-1.0, x)
293 x = x * acos( 0. )
294 delta_s = adjustment_deltas * 0.5 * (1. - sin( x ) )
295 do k=1,nz
296 s(i,j,k) = s_ref + delta_s + 0.5 * ( eta1d(k)+eta1d(k+1) ) * dsdz
297 x = abs(s(i,j,k) - 0.5*real(nz-1)/real(nz)*s_range)/s_range*real(2*nz)
298 x = 1. - min(1., x)
299 t(i,j,k) = t_range * x
300 enddo
301 ! x = sum(T(i,j,:)*h(i,j,:))
302 ! T(i,j,:) = (T(i,j,:) / x) * (G%max_depth*1.5/real(nz))
303 enddo ; enddo
304
305 case ( regridding_layer, regridding_rho )
306 do k = 1,nz
307 s(:,:,k) = s_ref + s_range * ( (real(k)-0.5) / real( nz ) )
308 ! x = abs(S(1,1,k) - 0.5*real(nz-1)/real(nz)*S_range)/S_range*real(2*nz)
309 ! x = 1.-min(1., x)
310 ! T(:,:,k) = T_range * x
311 enddo
312
313 case default
314 call mom_error(fatal, "adjustment_initialize_temperature_salinity: "// &
315 "Unrecognized i.c. setup - set ADJUSTMENT_IC")
316
317 end select
318
320