MOM6
adjustment_initialization.F90
Go to the documentation of this file.
1 !> Configures the model for the geostrophic adjustment test case.
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
8 use mom_get_input, only : directories
9 use mom_grid, only : ocean_grid_type
17 
18 implicit none ; private
19 
20 character(len=40) :: mdl = "adjustment_initialization" !< This module's name.
21 
22 #include <MOM_memory.h>
23 
26 
27 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
28 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
29 ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
30 ! vary with the Boussinesq approximation, the Boussinesq variant is given first.
31 
32 contains
33 
34 !> Initializes the layer thicknesses in the adjustment test case
35 subroutine adjustment_initialize_thickness ( h, G, GV, US, param_file, just_read_params)
36  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
37  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
38  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
39  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
40  intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2].
41  type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
42  !! to parse for model parameter values.
43  logical, optional, intent(in) :: just_read_params !< If present and true, this call will
44  !! only read parameters without changing h.
45  ! Local variables
46  real :: e0(szk_(g)+1) ! The resting interface heights, in depth units [Z ~> m], usually
47  ! negative because it is positive upward.
48  real :: eta1d(szk_(g)+1)! Interface height relative to the sea surface
49  ! positive upward, in depth units [Z ~> m].
50  real :: drho_ds ! The partial derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1].
51  ! In this subroutine it is hard coded at 1.0 kg m-3 ppt-1.
52  real :: x, y, yy
53  real :: delta_s_strat, dsdz, delta_s, s_ref
54  real :: min_thickness, adjustment_width, adjustment_delta
55  real :: adjustment_deltas
56  real :: front_wave_amp, front_wave_length, front_wave_asym
57  real :: target_values(szk_(g)+1) ! Target densities or density anomalies [R ~> kg m-3]
58  logical :: just_read ! If true, just read parameters but set nothing.
59  character(len=20) :: verticalcoordinate
60 ! This include declares and sets the variable "version".
61 #include "version_variable.h"
62  integer :: i, j, k, is, ie, js, je, nz
63 
64  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
65 
66  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params
67 
68  if (.not.just_read) &
69  call mom_mesg("initialize_thickness_uniform: setting thickness")
70 
71  ! Parameters used by main model initialization
72  if (.not.just_read) call log_version(param_file, mdl, version, "")
73  call get_param(param_file, mdl,"S_REF",s_ref,fail_if_missing=.true.,do_not_log=.true.)
74  call get_param(param_file, mdl,"MIN_THICKNESS",min_thickness,'Minimum layer thickness', &
75  units='m', default=1.0e-3, do_not_log=just_read, scale=us%m_to_Z)
76 
77  ! Parameters specific to this experiment configuration
78  call get_param(param_file, mdl,"REGRIDDING_COORDINATE_MODE",verticalcoordinate, &
79  default=default_coordinate_mode, do_not_log=just_read)
80  call get_param(param_file, mdl,"ADJUSTMENT_WIDTH",adjustment_width, &
81  "Width of frontal zone", &
82  units="same as x,y", fail_if_missing=.not.just_read, do_not_log=just_read)
83  call get_param(param_file, mdl,"DELTA_S_STRAT",delta_s_strat, &
84  "Top-to-bottom salinity difference of stratification", &
85  units="1e-3", fail_if_missing=.not.just_read, do_not_log=just_read)
86  call get_param(param_file, mdl,"ADJUSTMENT_DELTAS",adjustment_deltas, &
87  "Salinity difference across front", &
88  units="1e-3", fail_if_missing=.not.just_read, do_not_log=just_read)
89  call get_param(param_file, mdl,"FRONT_WAVE_AMP",front_wave_amp, &
90  "Amplitude of trans-frontal wave perturbation", &
91  units="same as x,y",default=0., do_not_log=just_read)
92  call get_param(param_file, mdl,"FRONT_WAVE_LENGTH",front_wave_length, &
93  "Wave-length of trans-frontal wave perturbation", &
94  units="same as x,y",default=0., do_not_log=just_read)
95  call get_param(param_file, mdl,"FRONT_WAVE_ASYM",front_wave_asym, &
96  "Amplitude of frontal asymmetric perturbation", &
97  default=0., do_not_log=just_read)
98 
99  if (just_read) return ! All run-time parameters have been read, so return.
100 
101  ! WARNING: this routine specifies the interface heights so that the last layer
102  ! is vanished, even at maximum depth. In order to have a uniform
103  ! layer distribution, use this line of code within the loop:
104  ! e0(k) = -G%max_depth * real(k-1) / real(nz)
105  ! To obtain a thickness distribution where the last layer is
106  ! vanished and the other thicknesses uniformly distributed, use:
107  ! e0(k) = -G%max_depth * real(k-1) / real(nz-1)
108 
109  dsdz = -delta_s_strat / g%max_depth
110 
111  select case ( coordinatemode(verticalcoordinate) )
112 
113  case ( regridding_layer, regridding_rho )
114  drho_ds = 1.0 * us%kg_m3_to_R
115  if (delta_s_strat /= 0.) then
116  ! This was previously coded ambiguously.
117  adjustment_delta = (adjustment_deltas / delta_s_strat) * g%max_depth
118  do k=1,nz+1
119  e0(k) = adjustment_delta - (g%max_depth + 2*adjustment_delta) * (real(k-1) / real(nz))
120  enddo
121  else
122  adjustment_delta = 2.*g%max_depth
123  do k=1,nz+1
124  e0(k) = -g%max_depth * (real(k-1) / real(nz))
125  enddo
126  endif
127  target_values(1) = ( gv%Rlay(1) + 0.5*(gv%Rlay(1)-gv%Rlay(2)) )
128  target_values(nz+1) = ( gv%Rlay(nz) + 0.5*(gv%Rlay(nz)-gv%Rlay(nz-1)) )
129  do k = 2,nz
130  target_values(k) = target_values(k-1) + ( gv%Rlay(nz) - gv%Rlay(1) ) / (nz-1)
131  enddo
132  target_values(:) = target_values(:) - 1000.*us%kg_m3_to_R
133  do j=js,je ; do i=is,ie
134  if (front_wave_length /= 0.) then
135  y = ( 0.125 + g%geoLatT(i,j) / front_wave_length ) * ( 4. * acos(0.) )
136  yy = 2. * ( g%geoLatT(i,j) - 0.5 * g%len_lat ) / adjustment_width
137  yy = min(1.0, yy); yy = max(-1.0, yy)
138  yy = yy * 2. * acos( 0. )
139  y = front_wave_amp*sin(y) + front_wave_asym*sin(yy)
140  else
141  y = 0.
142  endif
143  x = ( ( g%geoLonT(i,j) - 0.5 * g%len_lon ) + y ) / adjustment_width
144  x = min(1.0, x); x = max(-1.0, x)
145  x = x * acos( 0. )
146  delta_s = adjustment_deltas * 0.5 * (1. - sin( x ) )
147  do k=2,nz
148  if (drho_ds*dsdz /= 0.) then
149  eta1d(k) = ( target_values(k) - drho_ds*( s_ref + delta_s ) ) / (drho_ds*dsdz)
150  else
151  eta1d(k) = e0(k) - (0.5*adjustment_delta) * sin( x )
152  endif
153  eta1d(k) = max( eta1d(k), -g%max_depth )
154  eta1d(k) = min( eta1d(k), 0. )
155  enddo
156  eta1d(1) = 0.; eta1d(nz+1) = -g%max_depth
157  do k=nz,1,-1
158  if (eta1d(k) > 0.) then
159  eta1d(k) = max( eta1d(k+1) + min_thickness, 0. )
160  h(i,j,k) = gv%Z_to_H * max( eta1d(k) - eta1d(k+1), min_thickness )
161  elseif (eta1d(k) <= (eta1d(k+1) + min_thickness)) then
162  eta1d(k) = eta1d(k+1) + min_thickness
163  h(i,j,k) = gv%Z_to_H * min_thickness
164  else
165  h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
166  endif
167  enddo
168  enddo ; enddo
169 
170  case ( regridding_zstar, regridding_sigma )
171  do k=1,nz+1
172  eta1d(k) = -g%max_depth * (real(k-1) / real(nz))
173  eta1d(k) = max(min(eta1d(k), 0.), -g%max_depth)
174  enddo
175  do j=js,je ; do i=is,ie
176  do k=nz,1,-1
177  h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
178  enddo
179  enddo ; enddo
180 
181  case default
182  call mom_error(fatal,"adjustment_initialize_thickness: "// &
183  "Unrecognized i.c. setup - set ADJUSTMENT_IC")
184 
185  end select
186 
187 end subroutine adjustment_initialize_thickness
188 
189 !> Initialization of temperature and salinity in the adjustment test case
190 subroutine adjustment_initialize_temperature_salinity(T, S, h, G, GV, param_file, &
191  eqn_of_state, just_read_params)
192  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
193  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
194  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: t !< The temperature that is being initialized.
195  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: s !< The salinity that is being initialized.
196  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(in) :: h !< The model thicknesses [H ~> m or kg m-2].
197  type(param_file_type), intent(in) :: param_file !< A structure indicating the open file to
198  !! parse for model parameter values.
199  type(eos_type), pointer :: eqn_of_state !< Equation of state.
200  logical, optional, intent(in) :: just_read_params !< If present and true, this call will
201  !! only read parameters without changing T & S.
202 
203  integer :: i, j, k, is, ie, js, je, nz
204  real :: x, y, yy
205  integer :: index_bay_z
206  real :: s_ref, t_ref ! Reference salinity and temerature within
207  ! surface layer
208  real :: s_range, t_range ! Range of salinities and temperatures over the
209  ! vertical
210  real :: xi0, xi1, dsdz, delta_s, delta_s_strat
211  real :: adjustment_width, adjustment_deltas
212  real :: front_wave_amp, front_wave_length, front_wave_asym
213  real :: eta1d(szk_(g)+1)
214  logical :: just_read ! If true, just read parameters but set nothing.
215  character(len=20) :: verticalcoordinate
216 
217  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
218 
219  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params
220 
221  ! Parameters used by main model initialization
222  call get_param(param_file, mdl,"S_REF",s_ref,'Reference salinity', units='1e-3', &
223  fail_if_missing=.not.just_read, do_not_log=just_read)
224  call get_param(param_file, mdl,"T_REF",t_ref,'Reference temperature', units='C', &
225  fail_if_missing=.not.just_read, do_not_log=just_read)
226  call get_param(param_file, mdl,"S_RANGE",s_range,'Initial salinity range', units='1e-3', &
227  default=2.0, do_not_log=just_read)
228  call get_param(param_file, mdl,"T_RANGE",t_range,'Initial temperature range', units='C', &
229  default=0.0, do_not_log=just_read)
230  ! Parameters specific to this experiment configuration BUT logged in previous s/r
231  call get_param(param_file, mdl,"REGRIDDING_COORDINATE_MODE",verticalcoordinate, &
232  default=default_coordinate_mode, do_not_log=just_read)
233  call get_param(param_file, mdl,"ADJUSTMENT_WIDTH", adjustment_width, &
234  fail_if_missing=.not.just_read, do_not_log=.true.)
235  call get_param(param_file, mdl,"ADJUSTMENT_DELTAS", adjustment_deltas, &
236  fail_if_missing=.not.just_read, do_not_log=.true.)
237  call get_param(param_file, mdl,"DELTA_S_STRAT", delta_s_strat, &
238  fail_if_missing=.not.just_read, do_not_log=.true.)
239  call get_param(param_file, mdl,"FRONT_WAVE_AMP", front_wave_amp, default=0., &
240  do_not_log=.true.)
241  call get_param(param_file, mdl,"FRONT_WAVE_LENGTH",front_wave_length, &
242  default=0., do_not_log=.true.)
243  call get_param(param_file, mdl,"FRONT_WAVE_ASYM", front_wave_asym, default=0., &
244  do_not_log=.true.)
245 
246  if (just_read) return ! All run-time parameters have been read, so return.
247 
248  t(:,:,:) = 0.0
249  s(:,:,:) = 0.0
250 
251  ! Linear salinity profile
252  select case ( coordinatemode(verticalcoordinate) )
253 
254  case ( regridding_zstar, regridding_sigma )
255  dsdz = -delta_s_strat / g%max_depth
256  do j=js,je ; do i=is,ie
257  eta1d(nz+1) = -g%bathyT(i,j)
258  do k=nz,1,-1
259  eta1d(k) = eta1d(k+1) + h(i,j,k)*gv%H_to_Z
260  enddo
261  if (front_wave_length /= 0.) then
262  y = ( 0.125 + g%geoLatT(i,j) / front_wave_length ) * ( 4. * acos(0.) )
263  yy = 2. * ( g%geoLatT(i,j) - 0.5 * g%len_lat ) / front_wave_length
264  yy = min(1.0, yy); yy = max(-1.0, yy)
265  yy = yy * 2. * acos( 0. )
266  y = front_wave_amp*sin(y) + front_wave_asym*sin(yy)
267  else
268  y = 0.
269  endif
270  x = ( ( g%geoLonT(i,j) - 0.5 * g%len_lon ) + y ) / adjustment_width
271  x = min(1.0, x); x = max(-1.0, x)
272  x = x * acos( 0. )
273  delta_s = adjustment_deltas * 0.5 * (1. - sin( x ) )
274  do k=1,nz
275  s(i,j,k) = s_ref + delta_s + 0.5 * ( eta1d(k)+eta1d(k+1) ) * dsdz
276  x = abs(s(i,j,k) - 0.5*real(nz-1)/real(nz)*s_range)/s_range*real(2*nz)
277  x = 1. - min(1., x)
278  t(i,j,k) = x
279  enddo
280  ! x = GV%H_to_Z*sum(T(i,j,:)*h(i,j,:))
281  ! T(i,j,:) = (T(i,j,:) / x) * (G%max_depth*1.5/real(nz))
282  enddo ; enddo
283 
284  case ( regridding_layer, regridding_rho )
285  do k = 1,nz
286  s(:,:,k) = s_ref + s_range * ( (real(k)-0.5) / real( nz ) )
287  ! x = abs(S(1,1,k) - 0.5*real(nz-1)/real(nz)*S_range)/S_range*real(2*nz)
288  ! x = 1.-min(1., x)
289  ! T(:,:,k) = x
290  enddo
291 
292  case default
293  call mom_error(fatal,"adjustment_initialize_temperature_salinity: "// &
294  "Unrecognized i.c. setup - set ADJUSTMENT_IC")
295 
296  end select
297 
299 
300 end module adjustment_initialization
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_variables::thermo_var_ptrs
Pointers to an assortment of thermodynamic fields that may be available, including potential temperat...
Definition: MOM_variables.F90:78
mom_get_input::directories
Container for paths and parameter file names.
Definition: MOM_get_input.F90:20
mom_error_handler::mom_mesg
subroutine, public mom_mesg(message, verb, all_print)
This provides a convenient interface for writing an informative comment.
Definition: MOM_error_handler.F90:53
mom_eos
Provides subroutines for quantities specific to the equation of state.
Definition: MOM_EOS.F90:2
mom_file_parser::param_file_type
A structure that can be parsed to read and document run-time parameters.
Definition: MOM_file_parser.F90:54
mom_file_parser::get_param
An overloaded interface to read and log the values of various types of parameters.
Definition: MOM_file_parser.F90:102
adjustment_initialization::adjustment_initialize_thickness
subroutine, public adjustment_initialize_thickness(h, G, GV, US, param_file, just_read_params)
Initializes the layer thicknesses in the adjustment test case.
Definition: adjustment_initialization.F90:36
regrid_consts::regridding_layer
integer, parameter regridding_layer
Layer mode identifier.
Definition: regrid_consts.F90:13
mom_get_input
Reads the only Fortran name list needed to boot-strap the model.
Definition: MOM_get_input.F90:6
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
regrid_consts
Contains constants for interpreting input parameters that control regridding.
Definition: regrid_consts.F90:2
mom_eos::eos_type
A control structure for the equation of state.
Definition: MOM_EOS.F90:91
regrid_consts::regridding_sigma
integer, parameter regridding_sigma
Sigma coordinates identifier.
Definition: regrid_consts.F90:16
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
regrid_consts::regridding_zstar
integer, parameter regridding_zstar
z* coordinates identifier
Definition: regrid_consts.F90:14
mom_eos::calculate_density_derivs
Calculate the derivatives of density with temperature and salinity from T, S, and P.
Definition: MOM_EOS.F90:70
adjustment_initialization::mdl
character(len=40) mdl
This module's name.
Definition: adjustment_initialization.F90:20
adjustment_initialization::adjustment_initialize_temperature_salinity
subroutine, public adjustment_initialize_temperature_salinity(T, S, h, G, GV, param_file, eqn_of_state, just_read_params)
Initialization of temperature and salinity in the adjustment test case.
Definition: adjustment_initialization.F90:192
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_error_handler::is_root_pe
logical function, public is_root_pe()
This returns .true. if the current PE is the root PE.
Definition: MOM_error_handler.F90:44
regrid_consts::coordinatemode
integer function coordinatemode(string)
Parse a string parameter specifying the coordinate mode and return the appropriate enumerated integer...
Definition: regrid_consts.F90:54
regrid_consts::default_coordinate_mode
character(len= *), parameter default_coordinate_mode
Default coordinate mode.
Definition: regrid_consts.F90:35
mom_error_handler::mom_error
subroutine, public mom_error(level, message, all_print)
This provides a convenient interface for writing an mpp_error message with run-time filter based on a...
Definition: MOM_error_handler.F90:72
regrid_consts::regridding_rho
integer, parameter regridding_rho
Density coordinates identifier.
Definition: regrid_consts.F90:15
mom_file_parser::log_param
An overloaded interface to log the values of various types of parameters.
Definition: MOM_file_parser.F90:96
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:26
adjustment_initialization
Configures the model for the geostrophic adjustment test case.
Definition: adjustment_initialization.F90:2
mom_eos::calculate_density
Calculates density of sea water from T, S and P.
Definition: MOM_EOS.F90:60