Rossby_front_2d_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!> Initial conditions for the 2D Rossby front test
7
8use mom_error_handler, only : mom_mesg, mom_error, fatal, is_root_pe
9use mom_file_parser, only : get_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
21#include <MOM_memory.h>
22
23! Private (module-wise) parameters
24character(len=40) :: mdl = "Rossby_front_2d_initialization" !< This module's name.
25! This include declares and sets the variable "version".
26#include "version_variable.h"
27
31
32! Parameters defining the initial conditions of this test case
33real, parameter :: frontfractionalwidth = 0.5 !< Width of front as fraction of domain [nondim]
34real, parameter :: hmlmin = 0.25 !< Shallowest ML as fractional depth of ocean [nondim]
35real, parameter :: hmlmax = 0.75 !< Deepest ML as fractional depth of ocean [nondim]
36
37contains
38
39!> Initialization of thicknesses in 2D Rossby front test
40subroutine rossby_front_initialize_thickness(h, G, GV, US, param_file, just_read)
41 type(ocean_grid_type), intent(in) :: g !< Grid structure
42 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
43 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
44 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
45 intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2]
46 type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
47 !! to parse for model parameter values.
48 logical, intent(in) :: just_read !< If true, this call will only read
49 !! parameters without changing h.
50
51 ! Local variables
52 real :: tz ! Vertical temperature gradient [C H-1 ~> degC m2 kg-1]
53 real :: dml ! Mixed layer depth [H ~> m or kg m-2]
54 real :: eta ! An interface height depth [H ~> m or kg m-2]
55 real :: stretch ! A nondimensional stretching factor [nondim]
56 real :: h0 ! The stretched thickness per layer [H ~> m or kg m-2]
57 real :: t_range ! Range of temperatures over the vertical [C ~> degC]
58 real :: drho_dt ! The partial derivative of density with temperature [R C-1 ~> kg m-3 degC-1]
59 real :: max_depth ! Maximum depth of the model bathymetry [H ~> m or kg m-2]
60 character(len=40) :: verticalcoordinate
61 integer :: i, j, k, is, ie, js, je, nz
62
63 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
64
65 if (.not.just_read) &
66 call mom_mesg("Rossby_front_2d_initialization.F90, Rossby_front_initialize_thickness: setting thickness")
67
68 if (.not.just_read) call log_version(param_file, mdl, version, "")
69 ! Read parameters needed to set thickness
70 call get_param(param_file, mdl, "REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
71 default=default_coordinate_mode, do_not_log=just_read)
72 call get_param(param_file, mdl, "T_RANGE", t_range, 'Initial temperature range', &
73 units="degC", default=0.0, scale=us%degC_to_C, do_not_log=just_read)
74 call get_param(param_file, mdl, "DRHO_DT", drho_dt, &
75 units="kg m-3 degC-1", default=-0.2, scale=us%kg_m3_to_R*us%C_to_degC, do_not_log=.true.)
76 call get_param(param_file, mdl, "MAXIMUM_DEPTH", max_depth, &
77 units="m", default=-1.e9, scale=gv%m_to_H, do_not_log=.true.)
78
79 if (just_read) return ! All run-time parameters have been read, so return.
80
81 if (max_depth <= 0.0) call mom_error(fatal, &
82 "Rossby_front_initialize_thickness, Rossby_front_initialize_thickness: "//&
83 "This module requires a positive value of MAXIMUM_DEPTH.")
84
85 tz = t_range / max_depth
86
87 if (gv%Boussinesq) then
88 select case ( coordinatemode(verticalcoordinate) )
89
90 case (regridding_layer, regridding_rho)
91 ! This code is identical to the REGRIDDING_ZSTAR case but probably should not be.
92 do j = g%jsc,g%jec ; do i = g%isc,g%iec
93 dml = hml( g, g%geoLatT(i,j), max_depth )
94 eta = -( -drho_dt / gv%Rho0 ) * tz * 0.5 * ( dml * dml )
95 stretch = ( ( max_depth + eta ) / max_depth )
96 h0 = ( max_depth / real(nz) ) * stretch
97 do k = 1, nz
98 h(i,j,k) = h0
99 enddo
100 enddo ; enddo
101
102 case (regridding_zstar, regridding_sigma)
103 do j = g%jsc,g%jec ; do i = g%isc,g%iec
104 dml = hml( g, g%geoLatT(i,j), max_depth )
105 ! The free surface height is set so that the bottom pressure gradient is 0.
106 eta = -( -drho_dt / gv%Rho0 ) * tz * 0.5 * ( dml * dml )
107 stretch = ( ( max_depth + eta ) / max_depth )
108 h0 = ( max_depth / real(nz) ) * stretch
109 do k = 1, nz
110 h(i,j,k) = h0
111 enddo
112 enddo ; enddo
113
114 case default
115 call mom_error(fatal,"Rossby_front_initialize: "// &
116 "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE")
117
118 end select
119 else
120 ! In non-Boussinesq mode with a flat bottom, the only requirement for no bottom pressure
121 ! gradient and no abyssal flow is that all columns have the same mass.
122 h0 = max_depth / real(nz)
123 do k=1,nz ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
124 h(i,j,k) = h0
125 enddo ; enddo ; enddo
126 endif
127
129
130
131!> Initialization of temperature and salinity in the Rossby front test
132subroutine rossby_front_initialize_temperature_salinity(T, S, h, G, GV, US, &
133 param_file, just_read)
134 type(ocean_grid_type), intent(in) :: g !< Grid structure
135 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
136 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: t !< Potential temperature [C ~> degC]
137 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: s !< Salinity [S ~> ppt]
138 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Thickness [H ~> m or kg m-2]
139 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
140 type(param_file_type), intent(in) :: param_file !< Parameter file handle
141 logical, intent(in) :: just_read !< If true, this call will
142 !! only read parameters without changing T & S.
143 ! Local variables
144 real :: t_ref ! Reference temperature within the surface layer [C ~> degC]
145 real :: s_ref ! Reference salinity within the surface layer [S ~> ppt]
146 real :: t_range ! Range of temperatures over the vertical [C ~> degC]
147 real :: zc ! Position of the middle of the cell [H ~> m or kg m-2]
148 real :: zi ! Bottom interface position relative to the sea surface [H ~> m or kg m-2]
149 real :: dtdz ! Vertical temperature gradient [C H-1 ~> degC m-1 or degC m2 kg-1]
150 real :: dml ! Mixed layer depth [H ~> m or kg m-2]
151 real :: max_depth ! Maximum depth of the model bathymetry [H ~> m or kg m-2]
152 character(len=40) :: verticalcoordinate
153 integer :: i, j, k, is, ie, js, je, nz
154
155 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
156
157 call get_param(param_file, mdl,"REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
158 default=default_coordinate_mode, do_not_log=just_read)
159 call get_param(param_file, mdl, "S_REF", s_ref, 'Reference salinity', &
160 default=35.0, units="ppt", scale=us%ppt_to_S, do_not_log=just_read)
161 call get_param(param_file, mdl, "T_REF", t_ref, 'Reference temperature', &
162 units="degC", scale=us%degC_to_C, fail_if_missing=.not.just_read, do_not_log=just_read)
163 call get_param(param_file, mdl, "T_RANGE", t_range, 'Initial temperature range', &
164 units="degC", default=0.0, scale=us%degC_to_C, do_not_log=just_read)
165 call get_param(param_file, mdl, "MAXIMUM_DEPTH", max_depth, &
166 units="m", default=-1.e9, scale=gv%m_to_H, do_not_log=.true.)
167
168 if (just_read) return ! All run-time parameters have been read, so return.
169
170 if (max_depth <= 0.0) call mom_error(fatal, &
171 "Rossby_front_initialize_thickness, Rossby_front_initialize_temperature_salinity: "//&
172 "This module requires a positive value of MAXIMUM_DEPTH.")
173
174 t(:,:,:) = 0.0
175 s(:,:,:) = s_ref
176 dtdz = t_range / max_depth
177
178 ! This sets the temperature to the value at the base of the specified mixed layer
179 ! depth from a horizontally uniform constant thermal stratification.
180 do j = g%jsc,g%jec ; do i = g%isc,g%iec
181 zi = 0.
182 dml = hml(g, g%geoLatT(i,j), max_depth)
183 do k = 1, nz
184 zi = zi - h(i,j,k) ! Bottom interface position
185 zc = zi - 0.5*h(i,j,k) ! Position of middle of cell
186 t(i,j,k) = t_ref + dtdz * min( zc, -dml ) ! Linear temperature profile below the mixed layer
187 enddo
188 enddo ; enddo
189
191
192
193!> Initialization of u and v in the Rossby front test
194subroutine rossby_front_initialize_velocity(u, v, h, G, GV, US, param_file, just_read)
195 type(ocean_grid_type), intent(in) :: g !< Grid structure
196 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
197 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
198 intent(out) :: u !< i-component of velocity [L T-1 ~> m s-1]
199 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
200 intent(out) :: v !< j-component of velocity [L T-1 ~> m s-1]
201 real, dimension(SZI_(G),SZJ_(G), SZK_(GV)), &
202 intent(in) :: h !< Thickness [H ~> m or kg m-2]
203 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
204 type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
205 !! to parse for model parameter values.
206 logical, intent(in) :: just_read !< If present and true, this call will only
207 !! read parameters without setting u & v.
208
209 real :: t_range ! Range of temperatures over the vertical [C ~> degC]
210 real :: t_ref ! Reference temperature within the surface layer [C ~> degC]
211 real :: s_ref ! Reference salinity within the surface layer [S ~> ppt]
212 real :: dudt ! Factor to convert dT/dy into dU/dz, g*alpha/f with rescaling
213 ! [L2 H-1 T-1 C-1 ~> m s-1 degC-1 or m4 kg-1 s-1 degC-1]
214 real :: rho_t0_s0 ! The density at T=0, S=0 [R ~> kg m-3]
215 real :: drho_dt ! The partial derivative of density with temperature [R C-1 ~> kg m-3 degC-1]
216 real :: drho_ds ! The partial derivative of density with salinity [R S-1 ~> kg m-3 ppt-1]
217 real :: dspv_dt ! The partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1]
218 real :: t_here ! The temperature in the middle of a layer [C ~> degC]
219 real :: dtdz ! Vertical temperature gradient [C H-1 ~> degC m-1 or degC m2 kg-1]
220 real :: dml ! Mixed layer depth [H ~> m or kg m-2]
221 real :: zi, zc, zm ! Depths in thickness units [H ~> m or kg m-2].
222 real :: f ! The local Coriolis parameter [T-1 ~> s-1]
223 real :: i_f ! The Adcroft reciprocal of the local Coriolis parameter [T ~> s]
224 real :: ty ! The meridional temperature gradient [C L-1 ~> degC m-1]
225 real :: hatu ! Interpolated layer thickness in height units [H ~> m or kg m-2].
226 real :: u_int ! The zonal velocity at an interface [L T-1 ~> m s-1]
227 real :: max_depth ! Maximum depth of the model bathymetry [H ~> m or kg m-2]
228 integer :: i, j, k, is, ie, js, je, nz
229 character(len=40) :: verticalcoordinate
230
231 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
232
233 call get_param(param_file, mdl, "REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
234 default=default_coordinate_mode, do_not_log=just_read)
235 call get_param(param_file, mdl, "T_RANGE", t_range, 'Initial temperature range', &
236 units="degC", default=0.0, scale=us%degC_to_C, do_not_log=just_read)
237 call get_param(param_file, mdl, "S_REF", s_ref, 'Reference salinity', &
238 default=35.0, units="ppt", scale=us%ppt_to_S, do_not_log=.true.)
239 call get_param(param_file, mdl, "T_REF", t_ref, 'Reference temperature', &
240 units="degC", scale=us%degC_to_C, fail_if_missing=.not.just_read, do_not_log=.true.)
241 call get_param(param_file, mdl, "RHO_T0_S0", rho_t0_s0, &
242 units="kg m-3", default=1000.0, scale=us%kg_m3_to_R, do_not_log=.true.)
243 call get_param(param_file, mdl, "DRHO_DT", drho_dt, &
244 units="kg m-3 degC-1", default=-0.2, scale=us%kg_m3_to_R*us%C_to_degC, do_not_log=.true.)
245 call get_param(param_file, mdl, "DRHO_DS", drho_ds, &
246 units="kg m-3 ppt-1", default=0.8, scale=us%kg_m3_to_R*us%S_to_ppt, do_not_log=.true.)
247 call get_param(param_file, mdl, "MAXIMUM_DEPTH", max_depth, &
248 units="m", default=-1.e9, scale=gv%m_to_H, do_not_log=.true.)
249
250 if (just_read) return ! All run-time parameters have been read, so return.
251
252 if (max_depth <= 0.0) call mom_error(fatal, &
253 "Rossby_front_initialize_thickness, Rossby_front_initialize_velocity: "//&
254 "This module requires a positive value of MAXIMUM_DEPTH.")
255 if (g%grid_unit_to_L <= 0.) call mom_error(fatal, 'Rossby_front_2d_initialization.F90: '// &
256 "dTdy() is only set to work with Cartesian axis units.")
257
258 v(:,:,:) = 0.0
259 u(:,:,:) = 0.0
260
261 if (gv%Boussinesq) then
262 do j = g%jsc,g%jec ; do i = g%isc-1,g%iec+1
263 f = 0.5* (g%CoriolisBu(i,j) + g%CoriolisBu(i,j-1) )
264 dudt = 0.0 ; if (abs(f) > 0.0) &
265 dudt = ( gv%H_to_Z*gv%g_Earth*drho_dt ) / ( f * gv%Rho0 )
266 dml = hml( g, g%geoLatCu(i,j), max_depth )
267 ty = dtdy( g, t_range, g%geoLatCu(i,j), us )
268 zi = 0.
269 do k = 1, nz
270 hatu = 0.5 * (h(i,j,k) + h(i+1,j,k))
271 zi = zi - hatu ! Bottom interface position
272 zc = zi - 0.5*hatu ! Position of middle of cell
273 zm = max( zc + dml, 0. ) ! Height above bottom of mixed layer
274 u(i,j,k) = dudt * ty * zm ! Thermal wind starting at base of ML
275 enddo
276 enddo ; enddo
277 else
278 ! With an equation of state that is linear in density, the nonlinearies in
279 ! specific volume require that temperature be calculated for each layer.
280
281 dtdz = t_range / max_depth
282
283 do j = g%jsc,g%jec ; do i = g%isc-1,g%iec+1
284 f = 0.5* (g%CoriolisBu(i,j) + g%CoriolisBu(i,j-1) )
285 i_f = 0.0 ; if (abs(f) > 0.0) i_f = 1.0 / f
286 dml = hml( g, g%geoLatCu(i,j), max_depth )
287 ty = dtdy( g, t_range, g%geoLatCu(i,j), us )
288 zi = -max_depth
289 u_int = 0.0 ! The velocity at an interface
290 ! Work upward in non-Boussinesq mode
291 do k = nz, 1, -1
292 hatu = 0.5 * (h(i,j,k) + h(i+1,j,k))
293 zc = zi + 0.5*hatu ! Position of middle of cell
294 t_here = t_ref + dtdz * min(zc, -dml) ! Linear temperature profile below the mixed layer
295 dspv_dt = -drho_dt / (rho_t0_s0 + (drho_ds * s_ref + drho_dt * t_here) )**2
296 dudt = -( gv%H_to_RZ * gv%g_Earth * dspv_dt ) * i_f
297
298 ! There is thermal wind shear only within the mixed layer.
299 u(i,j,k) = u_int + dudt * ty * min(max((zi + dml) + 0.5*hatu, 0.0), 0.5*hatu)
300 u_int = u_int + dudt * ty * min(max((zi + dml) + hatu, 0.0), hatu)
301
302 zi = zi + hatu ! Update the layer top interface position
303 enddo
304 enddo ; enddo
305 endif
307
308!> Pseudo coordinate across domain used by Hml() and dTdy()
309!! returns a coordinate from -PI/2 .. PI/2 squashed towards the
310!! center of the domain [radians].
311real function ypseudo( G, lat )
312 type(ocean_grid_type), intent(in) :: g !< Grid structure
313 real, intent(in) :: lat !< Latitude in arbitrary units, often [km]
314 ! Local
315 real :: pi ! The ratio of the circumference of a circle to its diameter [nondim]
316
317 pi = 4.0 * atan(1.0)
318 ypseudo = ( ( lat - g%south_lat ) / g%len_lat ) - 0.5 ! -1/2 .. 1/.2
319 ypseudo = pi * max(-0.5, min(0.5, ypseudo / frontfractionalwidth))
320end function ypseudo
321
322
323!> Analytic prescription of mixed layer depth in 2d Rossby front test,
324!! in the same units as max_depth (usually [Z ~> m] or [H ~> m or kg m-2])
325real function hml( G, lat, max_depth )
326 type(ocean_grid_type), intent(in) :: g !< Grid structure
327 real, intent(in) :: lat !< Latitude in arbitrary units, often [km]
328 real, intent(in) :: max_depth !< The maximum depth of the ocean [Z ~> m] or [H ~> m or kg m-2]
329 ! Local
330 real :: dhml, hmlmean ! The range and mean of the mixed layer depths [Z ~> m] or [H ~> m or kg m-2]
331
332 dhml = 0.5 * ( hmlmax - hmlmin ) * max_depth
333 hmlmean = 0.5 * ( hmlmin + hmlmax ) * max_depth
334 hml = hmlmean + dhml * sin( ypseudo(g, lat) )
335end function hml
336
337
338!> Analytic prescription of mixed layer temperature gradient in [C L-1 ~> degC m-1] in 2d Rossby front test
339real function dtdy( G, dT, lat, US )
340 type(ocean_grid_type), intent(in) :: g !< Grid structure
341 real, intent(in) :: dt !< Top to bottom temperature difference [C ~> degC]
342 real, intent(in) :: lat !< Latitude in the same units as geoLat, often [km]
343 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
344 ! Local
345 real :: pi ! The ratio of the circumference of a circle to its diameter [nondim]
346 real :: dhml ! The range of the mixed layer depths [Z ~> m]
347 real :: dhdy ! The mixed layer depth gradient [Z L-1 ~> m m-1]
348
349 pi = 4.0 * atan(1.0)
350 dhml = 0.5 * ( hmlmax - hmlmin ) * g%max_depth
351 dhdy = dhml * ( pi / ( frontfractionalwidth * g%len_lat * g%grid_unit_to_L ) ) * cos( ypseudo(g, lat) )
352 dtdy = -( dt / g%max_depth ) * dhdy
353
354end function dtdy
355
356
357!> \namespace rossby_front_2d_initialization
358!!
359!! \section section_Rossby_front_2d Description of the 2d Rossby front initial conditions
360!!
361!! Consistent with a linear equation of state, the system has a constant stratification
362!! below the mixed layer, stratified in temperature only. Isotherms are flat below the
363!! mixed layer and vertical within. Salinity is constant. The mixed layer has a half sine
364!! form so that there are no mixed layer or temperature gradients at the side walls.
365!!
366!! Below the mixed layer the potential temperature, \f$\theta(z)\f$, is given by
367!! \f[ \theta(z) = \theta_0 - \Delta \theta \left( z + h_{ML} \right) \f]
368!! where \f$ \theta_0 \f$ and \f$ \Delta \theta \f$ are external model parameters.
369!!
370!! The depth of the mixed layer, \f$H_{ML}\f$ is
371!! \f[ h_{ML}(y) = h_{min} + \left( h_{max} - h_{min} \right) \cos{\pi y/L} \f].
372!! The temperature in mixed layer is given by the reference temperature at \f$z=h_{ML}\f$
373!! so that
374!! \f{eqnarray} \theta(y,z) =
375!! \theta_0 - \Delta \theta \left( z + h_{ML} \right) & \forall & z < h_{ML}(y) T.B.D.
376!! \f}
377