18 use regrid_consts,
only : coordinatemode, default_coordinate_mode
21 use regrid_consts,
only : regridding_arbitrary, regridding_sigma_shelf_zstar
22 use regrid_consts,
only : regridding_hycom1, regridding_slight, regridding_adaptive
25 use coord_zlike,
only : init_coord_zlike,
zlike_cs, set_zlike_params, build_zstar_column, end_coord_zlike
26 use coord_sigma,
only : init_coord_sigma,
sigma_cs, set_sigma_params, build_sigma_column, end_coord_sigma
27 use coord_rho,
only : init_coord_rho,
rho_cs, set_rho_params, build_rho_column, end_coord_rho
28 use coord_rho,
only : old_inflate_layers_1d
29 use coord_hycom,
only : init_coord_hycom,
hycom_cs, set_hycom_params, build_hycom1_column, end_coord_hycom
30 use coord_slight,
only : init_coord_slight,
slight_cs, set_slight_params, build_slight_column, end_coord_slight
31 use coord_adapt,
only : init_coord_adapt,
adapt_cs, set_adapt_params, build_adapt_column, end_coord_adapt
35 implicit none ;
private
37 #include <MOM_memory.h>
51 real,
dimension(:),
allocatable :: coordinateresolution
55 real :: coord_scale = 1.0
62 real,
dimension(:),
allocatable :: target_density
65 logical :: target_density_set = .false.
69 real,
dimension(:),
allocatable :: max_interface_depths
73 real,
dimension(:),
allocatable :: max_layer_thickness
79 integer :: regridding_scheme
88 real :: ref_pressure = 2.e7
94 real :: old_grid_weight = 0.
97 real :: depth_of_time_filter_shallow = 0.
100 real :: depth_of_time_filter_deep = 0.
104 real :: compressibility_fraction = 0.
108 logical :: set_maximum_depths = .false.
112 real :: max_depth_index_scale = 2.0
116 logical :: integrate_downward_for_e = .true.
128 public initialize_regridding, end_regridding, regridding_main
129 public inflate_vanished_layers_old, check_remapping_grid, check_grid_column
130 public set_regrid_params, get_regrid_size
131 public uniformresolution, setcoordinateresolution
132 public build_rho_column
133 public set_target_densities_from_gv, set_target_densities
134 public set_regrid_max_depths, set_regrid_max_thickness
135 public getcoordinateresolution, getcoordinateinterfaces
136 public getcoordinateunits, getcoordinateshortname, getstaticthickness
137 public default_coordinate_mode
138 public get_zlike_cs, get_sigma_cs, get_rho_cs
141 character(len=*),
parameter,
public :: regriddingcoordinatemodedoc = &
142 " LAYER - Isopycnal or stacked shallow water layers\n"//&
143 " ZSTAR, Z* - stretched geopotential z*\n"//&
144 " SIGMA_SHELF_ZSTAR - stretched geopotential z* ignoring shelf\n"//&
145 " SIGMA - terrain following coordinates\n"//&
146 " RHO - continuous isopycnal\n"//&
147 " HYCOM1 - HyCOM-like hybrid coordinate\n"//&
148 " SLIGHT - stretched coordinates above continuous isopycnal\n"//&
149 " ADAPTIVE - optimize for smooth neutral density surfaces"
152 character(len=*),
parameter,
public :: regriddinginterpschemedoc = &
153 " P1M_H2 (2nd-order accurate)\n"//&
154 " P1M_H4 (2nd-order accurate)\n"//&
155 " P1M_IH4 (2nd-order accurate)\n"//&
156 " PLM (2nd-order accurate)\n"//&
157 " PPM_H4 (3rd-order accurate)\n"//&
158 " PPM_IH4 (3rd-order accurate)\n"//&
159 " P3M_IH4IH3 (4th-order accurate)\n"//&
160 " P3M_IH6IH5 (4th-order accurate)\n"//&
161 " PQM_IH4IH3 (4th-order accurate)\n"//&
162 " PQM_IH6IH5 (5th-order accurate)"
165 character(len=*),
parameter,
public :: regriddingdefaultinterpscheme =
"P1M_H2"
167 logical,
parameter,
public :: regriddingdefaultboundaryextrapolation = .false.
169 real,
parameter,
public :: regriddingdefaultminthickness = 1.e-3
171 #undef __DO_SAFETY_CHECKS__
176 subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_mode, param_prefix, param_suffix)
180 real,
intent(in) :: max_depth
182 character(len=*),
intent(in) :: mdl
183 character(len=*),
intent(in) :: coord_mode
184 character(len=*),
intent(in) :: param_prefix
186 character(len=*),
intent(in) :: param_suffix
190 character(len=80) :: string, string2, varname
191 character(len=40) :: coord_units, param_name, coord_res_param
192 character(len=200) :: inputdir, filename
193 character(len=320) :: message
194 character(len=12) :: expected_units
195 logical :: tmplogical, fix_haloclines, set_max, do_sum, main_parameters
196 logical :: coord_is_state_dependent, ierr
197 real :: filt_len, strat_tol, index_scale, tmpreal
198 real :: maximum_depth
199 real :: dz_fixed_sfc, rho_avg_depth, nlay_sfc_int
200 real :: adapttimeratio, adaptzoom, adaptzoomcoeff, adaptbuoycoeff, adaptalpha
201 integer :: nz_fixed_sfc, k, nzf(4)
202 real,
dimension(:),
allocatable :: dz
204 real,
dimension(:),
allocatable :: h_max
205 real,
dimension(:),
allocatable :: z_max
207 real,
dimension(:),
allocatable :: dz_max
209 real,
dimension(:),
allocatable :: rho_target
211 real,
dimension(40) :: woa09_dz = (/ 5., 10., 10., 15., 22.5, 25., 25., 25., &
212 37.5, 50., 50., 75., 100., 100., 100., 100., &
213 100., 100., 100., 100., 100., 100., 100., 175., &
214 250., 375., 500., 500., 500., 500., 500., 500., &
215 500., 500., 500., 500., 500., 500., 500., 500. /)
217 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
218 inputdir = slasher(inputdir)
220 main_parameters=.false.
221 if (len_trim(param_prefix)==0) main_parameters=.true.
222 if (main_parameters .and. len_trim(param_suffix)>0)
call mom_error(fatal,trim(mdl)//
', initialize_regridding: '// &
223 'Suffix provided without prefix for parameter names!')
226 cs%regridding_scheme = coordinatemode(coord_mode)
228 maximum_depth = us%Z_to_m*max_depth
230 if (main_parameters)
then
232 call get_param(param_file, mdl,
"REGRIDDING_COORDINATE_UNITS", coord_units, &
233 "Units of the regridding coordinate.", default=
coordinateunits(coord_mode))
238 if (coord_is_state_dependent)
then
239 if (main_parameters)
then
240 param_name =
"INTERPOLATION_SCHEME"
241 string2 = regriddingdefaultinterpscheme
243 param_name = trim(param_prefix)//
"_INTERP_SCHEME_"//trim(param_suffix)
246 call get_param(param_file, mdl,
"INTERPOLATION_SCHEME", string, &
247 "This sets the interpolation scheme to use to "//&
248 "determine the new grid. These parameters are "//&
249 "only relevant when REGRIDDING_COORDINATE_MODE is "//&
250 "set to a function of state. Otherwise, it is not "//&
251 "used. It can be one of the following schemes: "//&
252 trim(regriddinginterpschemedoc), default=trim(string2))
253 call set_regrid_params(cs, interp_scheme=string)
256 if (main_parameters .and. coord_is_state_dependent)
then
257 call get_param(param_file, mdl,
"BOUNDARY_EXTRAPOLATION", tmplogical, &
258 "When defined, a proper high-order reconstruction "//&
259 "scheme is used within boundary cells rather "//&
260 "than PCM. E.g., if PPM is used for remapping, a "//&
261 "PPM reconstruction will also be used within "//&
262 "boundary cells.", default=regriddingdefaultboundaryextrapolation)
263 call set_regrid_params(cs, boundary_extrapolation=tmplogical)
265 call set_regrid_params(cs, boundary_extrapolation=.false.)
269 if (main_parameters)
then
270 param_name =
"ALE_COORDINATE_CONFIG"
271 coord_res_param =
"ALE_RESOLUTION"
274 param_name = trim(param_prefix)//
"_DEF_"//trim(param_suffix)
275 coord_res_param = trim(param_prefix)//
"_RES_"//trim(param_suffix)
277 if (maximum_depth>3000.) string2=
'WOA09'
279 call get_param(param_file, mdl, param_name, string, &
280 "Determines how to specify the coordinate "//&
281 "resolution. Valid options are:\n"//&
282 " PARAM - use the vector-parameter "//trim(coord_res_param)//
"\n"//&
283 " UNIFORM[:N] - uniformly distributed\n"//&
284 " FILE:string - read from a file. The string specifies\n"//&
285 " the filename and variable name, separated\n"//&
286 " by a comma or space, e.g. FILE:lev.nc,dz\n"//&
287 " or FILE:lev.nc,interfaces=zw\n"//&
288 " WOA09[:N] - the WOA09 vertical grid (approximately)\n"//&
289 " FNC1:string - FNC1:dz_min,H_total,power,precision\n"//&
290 " HYBRID:string - read from a file. The string specifies\n"//&
291 " the filename and two variable names, separated\n"//&
292 " by a comma or space, for sigma-2 and dz. e.g.\n"//&
293 " HYBRID:vgrid.nc,sigma2,dz",&
294 default=trim(string2))
295 message =
"The distribution of vertical resolution for the target\n"//&
296 "grid used for Eulerian-like coordinates. For example,\n"//&
297 "in z-coordinate mode, the parameter is a list of level\n"//&
298 "thicknesses (in m). In sigma-coordinate mode, the list\n"//&
299 "is of non-dimensional fractions of the water column."
300 if (index(trim(string),
'UNIFORM')==1)
then
301 if (len_trim(string)==7)
then
303 tmpreal = maximum_depth
304 elseif (index(trim(string),
'UNIFORM:')==1 .and. len_trim(string)>8)
then
306 ke = extract_integer(string(9:len_trim(string)),
'',1)
307 tmpreal = extract_real(string(9:len_trim(string)),
',',2,missing_value=maximum_depth)
309 call mom_error(fatal,trim(mdl)//
', initialize_regridding: '// &
310 'Unable to interpret "'//trim(string)//
'".')
314 dz(:) = uniformresolution(ke, coord_mode, tmpreal, gv%Rlay(1), gv%Rlay(1))
316 dz(:) = uniformresolution(ke, coord_mode, tmpreal, &
317 gv%Rlay(1)+0.5*(gv%Rlay(1)-gv%Rlay(2)), &
318 gv%Rlay(ke)+0.5*(gv%Rlay(ke)-gv%Rlay(ke-1)) )
320 if (main_parameters)
call log_param(param_file, mdl,
"!"//coord_res_param, dz, &
321 trim(message), units=trim(coord_units))
322 elseif (trim(string)==
'PARAM')
then
326 call get_param(param_file, mdl, coord_res_param, dz, &
327 trim(message), units=trim(coord_units), fail_if_missing=.true.)
328 elseif (index(trim(string),
'FILE:')==1)
then
331 if (string(6:6)==
'.' .or. string(6:6)==
'/')
then
333 filename = trim( extractword(trim(string(6:80)), 1) )
336 filename = trim(inputdir) // trim( extractword(trim(string(6:80)), 1) )
338 if (.not.
file_exists(filename))
call mom_error(fatal,trim(mdl)//
", initialize_regridding: "// &
339 "Specified file not found: Looking for '"//trim(filename)//
"' ("//trim(string)//
")")
341 varname = trim( extractword(trim(string(6:)), 2) )
342 if (len_trim(varname)==0)
then
343 if (field_exists(filename,
'dz')) then; varname =
'dz'
344 elseif (field_exists(filename,
'dsigma')) then; varname =
'dsigma'
345 elseif (field_exists(filename,
'ztest')) then; varname =
'ztest'
346 else ;
call mom_error(fatal,trim(mdl)//
", initialize_regridding: "// &
347 "Coordinate variable not specified and none could be guessed.")
353 if (cs%regridding_scheme == regridding_sigma)
then
354 expected_units =
'nondim'
355 elseif (cs%regridding_scheme == regridding_rho)
then
356 expected_units =
'kg m-3'
358 expected_units =
'meters'
360 if (index(trim(varname),
'interfaces=')==1)
then
361 varname=trim(varname(12:))
362 call check_grid_def(filename, varname, expected_units, message, ierr)
363 if (ierr)
call mom_error(fatal,trim(mdl)//
", initialize_regridding: "//&
364 "Unsupported format in grid definition '"//trim(filename)//
"'. Error message "//trim(message))
365 call field_size(trim(filename), trim(varname), nzf)
367 if (cs%regridding_scheme == regridding_rho)
then
368 allocate(rho_target(ke+1))
369 call mom_read_data(trim(filename), trim(varname), rho_target)
372 allocate(z_max(ke+1))
374 dz(:) = abs(z_max(1:ke) - z_max(2:ke+1))
379 call field_size(trim(filename), trim(varname), nzf)
384 if (main_parameters .and. ke/=gv%ke)
then
385 call mom_error(fatal,trim(mdl)//
', initialize_regridding: '// &
386 'Mismatch in number of model levels and "'//trim(string)//
'".')
388 if (main_parameters)
call log_param(param_file, mdl,
"!"//coord_res_param, dz, &
390 elseif (index(trim(string),
'FNC1:')==1)
then
391 ke = gv%ke;
allocate(dz(ke))
392 call dz_function1( trim(string(6:)), dz )
393 if (main_parameters)
call log_param(param_file, mdl,
"!"//coord_res_param, dz, &
395 elseif (index(trim(string),
'RFNC1:')==1)
then
397 ke = rho_function1( trim(string(7:)), rho_target )
398 elseif (index(trim(string),
'HYBRID:')==1)
then
399 ke = gv%ke;
allocate(dz(ke))
401 allocate(rho_target(ke+1))
402 filename = trim( extractword(trim(string(8:)), 1) )
403 if (filename(1:1)/=
'.' .and. filename(1:1)/=
'/') filename = trim(inputdir) // trim( filename )
404 if (.not.
file_exists(filename))
call mom_error(fatal,trim(mdl)//
", initialize_regridding: HYBRID "// &
405 "Specified file not found: Looking for '"//trim(filename)//
"' ("//trim(string)//
")")
406 varname = trim( extractword(trim(string(8:)), 2) )
407 if (.not. field_exists(filename,varname))
call mom_error(fatal,trim(mdl)//
", initialize_regridding: HYBRID "// &
408 "Specified field not found: Looking for '"//trim(varname)//
"' ("//trim(string)//
")")
409 call mom_read_data(trim(filename), trim(varname), rho_target)
410 varname = trim( extractword(trim(string(8:)), 3) )
411 if (varname(1:5) ==
'FNC1:')
then
412 call dz_function1( trim(string((index(trim(string),
'FNC1:')+5):)), dz )
414 if (.not. field_exists(filename,varname))
call mom_error(fatal,trim(mdl)//
", initialize_regridding: HYBRID "// &
415 "Specified field not found: Looking for '"//trim(varname)//
"' ("//trim(string)//
")")
418 if (main_parameters)
then
419 call log_param(param_file, mdl,
"!"//coord_res_param, dz, &
421 call log_param(param_file, mdl,
"!TARGET_DENSITIES", rho_target, &
422 'HYBRID target densities for interfaces', units=
coordinateunits(coord_mode))
424 elseif (index(trim(string),
'WOA09')==1)
then
425 if (len_trim(string)==5)
then
426 tmpreal = 0. ; ke = 0
427 do while (tmpreal<maximum_depth)
429 tmpreal = tmpreal + woa09_dz(ke)
431 elseif (index(trim(string),
'WOA09:')==1)
then
432 if (len_trim(string)==6)
call mom_error(fatal,trim(mdl)//
', initialize_regridding: '// &
433 'Expected string of form "WOA09:N" but got "'//trim(string)//
'".')
434 ke = extract_integer(string(7:len_trim(string)),
'',1)
436 if (ke>40 .or. ke<1)
call mom_error(fatal,trim(mdl)//
', initialize_regridding: '// &
437 'For "WOA05:N" N must 0<N<41 but got "'//trim(string)//
'".')
439 dz(1:ke) = woa09_dz(1:ke)
441 call mom_error(fatal,trim(mdl)//
", initialize_regridding: "// &
442 "Unrecognized coordinate configuration"//trim(string))
445 if (main_parameters)
then
447 if (coordinatemode(coord_mode) == regridding_zstar .or. &
448 coordinatemode(coord_mode) == regridding_hycom1 .or. &
449 coordinatemode(coord_mode) == regridding_slight .or. &
450 coordinatemode(coord_mode) == regridding_adaptive)
then
452 tmpreal = sum( dz(:) )
453 if (tmpreal < maximum_depth)
then
454 dz(ke) = dz(ke) + ( maximum_depth - tmpreal )
455 elseif (tmpreal > maximum_depth)
then
456 if ( dz(ke) + ( maximum_depth - tmpreal ) > 0. )
then
457 dz(ke) = dz(ke) + ( maximum_depth - tmpreal )
459 call mom_error(fatal,trim(mdl)//
", initialize_regridding: "// &
460 "MAXIMUM_DEPTH was too shallow to adjust bottom layer of DZ!"//trim(string))
469 allocate( cs%coordinateResolution(cs%nk) ); cs%coordinateResolution(:) = -1.e30
472 allocate( cs%target_density(cs%nk+1) ); cs%target_density(:) = -1.e30
475 if (
allocated(dz))
then
476 if ((coordinatemode(coord_mode) == regridding_sigma) .or. &
477 (coordinatemode(coord_mode) == regridding_rho))
then
478 call setcoordinateresolution(dz, cs, scale=1.0)
479 elseif (coordinatemode(coord_mode) == regridding_adaptive)
then
480 call setcoordinateresolution(dz, cs, scale=gv%m_to_H)
481 cs%coord_scale = gv%H_to_m
483 call setcoordinateresolution(dz, cs, scale=us%m_to_Z)
484 cs%coord_scale = us%Z_to_m
488 if (
allocated(rho_target))
then
489 call set_target_densities(cs, rho_target)
490 deallocate(rho_target)
493 elseif (coordinatemode(coord_mode) == regridding_rho)
then
494 call set_target_densities_from_gv(gv, cs)
495 call log_param(param_file, mdl,
"!TARGET_DENSITIES", cs%target_density, &
496 'RHO target densities for interfaces', units=
coordinateunits(coord_mode))
500 call initcoord(cs, gv, coord_mode)
502 if (main_parameters .and. coord_is_state_dependent)
then
503 call get_param(param_file, mdl,
"REGRID_COMPRESSIBILITY_FRACTION", tmpreal, &
504 "When interpolating potential density profiles we can add "//&
505 "some artificial compressibility solely to make homogeneous "//&
506 "regions appear stratified.", default=0.)
507 call set_regrid_params(cs, compress_fraction=tmpreal)
510 if (main_parameters)
then
511 call get_param(param_file, mdl,
"MIN_THICKNESS", tmpreal, &
512 "When regridding, this is the minimum layer "//&
513 "thickness allowed.", units=
"m", scale=gv%m_to_H, &
514 default=regriddingdefaultminthickness )
515 call set_regrid_params(cs, min_thickness=tmpreal)
517 call set_regrid_params(cs, min_thickness=0.)
520 if (coordinatemode(coord_mode) == regridding_slight)
then
522 call get_param(param_file, mdl,
"SLIGHT_DZ_SURFACE", dz_fixed_sfc, &
523 "The nominal thickness of fixed thickness near-surface "//&
524 "layers with the SLight coordinate.", units=
"m", default=1.0, scale=gv%m_to_H)
525 call get_param(param_file, mdl,
"SLIGHT_NZ_SURFACE_FIXED", nz_fixed_sfc, &
526 "The number of fixed-depth surface layers with the SLight "//&
527 "coordinate.", units=
"nondimensional", default=2)
528 call get_param(param_file, mdl,
"SLIGHT_SURFACE_AVG_DEPTH", rho_avg_depth, &
529 "The thickness of the surface region over which to average "//&
530 "when calculating the density to use to define the interior "//&
531 "with the SLight coordinate.", units=
"m", default=1.0, scale=gv%m_to_H)
532 call get_param(param_file, mdl,
"SLIGHT_NLAY_TO_INTERIOR", nlay_sfc_int, &
533 "The number of layers to offset the surface density when "//&
534 "defining where the interior ocean starts with SLight.", &
535 units=
"nondimensional", default=2.0)
536 call get_param(param_file, mdl,
"SLIGHT_FIX_HALOCLINES", fix_haloclines, &
537 "If true, identify regions above the reference pressure "//&
538 "where the reference pressure systematically underestimates "//&
539 "the stratification and use this in the definition of the "//&
540 "interior with the SLight coordinate.", default=.false.)
542 call set_regrid_params(cs, dz_min_surface=dz_fixed_sfc, &
543 nz_fixed_surface=nz_fixed_sfc, rho_ml_avg_depth=rho_avg_depth, &
544 nlay_ml_to_interior=nlay_sfc_int, fix_haloclines=fix_haloclines)
545 if (fix_haloclines)
then
547 call get_param(param_file, mdl,
"HALOCLINE_FILTER_LENGTH", filt_len, &
548 "A length scale over which to smooth the temperature and "//&
549 "salinity before identifying erroneously unstable haloclines.", &
550 units=
"m", default=2.0)
551 call get_param(param_file, mdl,
"HALOCLINE_STRAT_TOL", strat_tol, &
552 "A tolerance for the ratio of the stratification of the "//&
553 "apparent coordinate stratification to the actual value "//&
554 "that is used to identify erroneously unstable haloclines. "//&
555 "This ratio is 1 when they are equal, and sensible values "//&
556 "are between 0 and 0.5.", units=
"nondimensional", default=0.2)
557 call set_regrid_params(cs, halocline_filt_len=filt_len, &
558 halocline_strat_tol=strat_tol)
563 if (coordinatemode(coord_mode) == regridding_adaptive)
then
564 call get_param(param_file, mdl,
"ADAPT_TIME_RATIO", adapttimeratio, &
565 "Ratio of ALE timestep to grid timescale.", units=
"s", default=1e-1)
566 call get_param(param_file, mdl,
"ADAPT_ZOOM_DEPTH", adaptzoom, &
567 "Depth of near-surface zooming region.", units=
"m", default=200.0, scale=gv%m_to_H)
568 call get_param(param_file, mdl,
"ADAPT_ZOOM_COEFF", adaptzoomcoeff, &
569 "Coefficient of near-surface zooming diffusivity.", &
570 units=
"nondim", default=0.2)
571 call get_param(param_file, mdl,
"ADAPT_BUOY_COEFF", adaptbuoycoeff, &
572 "Coefficient of buoyancy diffusivity.", &
573 units=
"nondim", default=0.8)
574 call get_param(param_file, mdl,
"ADAPT_ALPHA", adaptalpha, &
575 "Scaling on optimization tendency.", &
576 units=
"nondim", default=1.0)
577 call get_param(param_file, mdl,
"ADAPT_DO_MIN_DEPTH", tmplogical, &
578 "If true, make a HyCOM-like mixed layer by preventing interfaces "//&
579 "from being shallower than the depths specified by the regridding coordinate.", &
582 call set_regrid_params(cs, adapttimeratio=adapttimeratio, adaptzoom=adaptzoom, &
583 adaptzoomcoeff=adaptzoomcoeff, adaptbuoycoeff=adaptbuoycoeff, adaptalpha=adaptalpha, &
584 adaptdomin=tmplogical)
587 if (main_parameters .and. coord_is_state_dependent)
then
588 call get_param(param_file, mdl,
"MAXIMUM_INT_DEPTH_CONFIG", string, &
589 "Determines how to specify the maximum interface depths.\n"//&
590 "Valid options are:\n"//&
591 " NONE - there are no maximum interface depths\n"//&
592 " PARAM - use the vector-parameter MAXIMUM_INTERFACE_DEPTHS\n"//&
593 " FILE:string - read from a file. The string specifies\n"//&
594 " the filename and variable name, separated\n"//&
595 " by a comma or space, e.g. FILE:lev.nc,Z\n"//&
596 " FNC1:string - FNC1:dz_min,H_total,power,precision",&
598 message =
"The list of maximum depths for each interface."
599 allocate(z_max(ke+1))
601 if ( trim(string) ==
"NONE")
then
603 elseif ( trim(string) ==
"PARAM")
then
604 call get_param(param_file, mdl,
"MAXIMUM_INTERFACE_DEPTHS", z_max, &
605 trim(message), units=
"m", scale=gv%m_to_H, fail_if_missing=.true.)
606 call set_regrid_max_depths(cs, z_max)
607 elseif (index(trim(string),
'FILE:')==1)
then
608 if (string(6:6)==
'.' .or. string(6:6)==
'/')
then
610 filename = trim( extractword(trim(string(6:80)), 1) )
613 filename = trim(inputdir) // trim( extractword(trim(string(6:80)), 1) )
615 if (.not.
file_exists(filename))
call mom_error(fatal,trim(mdl)//
", initialize_regridding: "// &
616 "Specified file not found: Looking for '"//trim(filename)//
"' ("//trim(string)//
")")
619 varname = trim( extractword(trim(string(6:)), 2) )
620 if (.not. field_exists(filename,varname))
call mom_error(fatal,trim(mdl)//
", initialize_regridding: "// &
621 "Specified field not found: Looking for '"//trim(varname)//
"' ("//trim(string)//
")")
622 if (len_trim(varname)==0)
then
623 if (field_exists(filename,
'z_max')) then; varname =
'z_max'
624 elseif (field_exists(filename,
'dz')) then; varname =
'dz' ; do_sum = .true.
625 elseif (field_exists(filename,
'dz_max')) then; varname =
'dz_max' ; do_sum = .true.
626 else ;
call mom_error(fatal,trim(mdl)//
", initialize_regridding: "// &
627 "MAXIMUM_INT_DEPTHS variable not specified and none could be guessed.")
632 z_max(1) = 0.0 ;
do k=1,ke ; z_max(k+1) = z_max(k) + dz_max(k) ;
enddo
636 call log_param(param_file, mdl,
"!MAXIMUM_INT_DEPTHS", z_max, &
638 call set_regrid_max_depths(cs, z_max, gv%m_to_H)
639 elseif (index(trim(string),
'FNC1:')==1)
then
640 call dz_function1( trim(string(6:)), dz_max )
641 if ((coordinatemode(coord_mode) == regridding_slight) .and. &
642 (dz_fixed_sfc > 0.0))
then
643 do k=1,nz_fixed_sfc ; dz_max(k) = dz_fixed_sfc ;
enddo
645 z_max(1) = 0.0 ;
do k=1,ke ; z_max(k+1) = z_max(k) + dz_max(k) ;
enddo
646 call log_param(param_file, mdl,
"!MAXIMUM_INT_DEPTHS", z_max, &
648 call set_regrid_max_depths(cs, z_max, gv%m_to_H)
650 call mom_error(fatal,trim(mdl)//
", initialize_regridding: "// &
651 "Unrecognized MAXIMUM_INT_DEPTH_CONFIG "//trim(string))
658 call get_param(param_file, mdl,
"MAX_LAYER_THICKNESS_CONFIG", string, &
659 "Determines how to specify the maximum layer thicknesses.\n"//&
660 "Valid options are:\n"//&
661 " NONE - there are no maximum layer thicknesses\n"//&
662 " PARAM - use the vector-parameter MAX_LAYER_THICKNESS\n"//&
663 " FILE:string - read from a file. The string specifies\n"//&
664 " the filename and variable name, separated\n"//&
665 " by a comma or space, e.g. FILE:lev.nc,Z\n"//&
666 " FNC1:string - FNC1:dz_min,H_total,power,precision",&
668 message =
"The list of maximum thickness for each layer."
670 if ( trim(string) ==
"NONE")
then
672 elseif ( trim(string) ==
"PARAM")
then
673 call get_param(param_file, mdl,
"MAX_LAYER_THICKNESS", h_max, &
674 trim(message), units=
"m", fail_if_missing=.true., scale=gv%m_to_H)
675 call set_regrid_max_thickness(cs, h_max)
676 elseif (index(trim(string),
'FILE:')==1)
then
677 if (string(6:6)==
'.' .or. string(6:6)==
'/')
then
679 filename = trim( extractword(trim(string(6:80)), 1) )
682 filename = trim(inputdir) // trim( extractword(trim(string(6:80)), 1) )
684 if (.not.
file_exists(filename))
call mom_error(fatal,trim(mdl)//
", initialize_regridding: "// &
685 "Specified file not found: Looking for '"//trim(filename)//
"' ("//trim(string)//
")")
687 varname = trim( extractword(trim(string(6:)), 2) )
688 if (.not. field_exists(filename,varname))
call mom_error(fatal,trim(mdl)//
", initialize_regridding: "// &
689 "Specified field not found: Looking for '"//trim(varname)//
"' ("//trim(string)//
")")
690 if (len_trim(varname)==0)
then
691 if (field_exists(filename,
'h_max')) then; varname =
'h_max'
692 elseif (field_exists(filename,
'dz_max')) then; varname =
'dz_max'
693 else ;
call mom_error(fatal,trim(mdl)//
", initialize_regridding: "// &
694 "MAXIMUM_INT_DEPTHS variable not specified and none could be guessed.")
698 call log_param(param_file, mdl,
"!MAX_LAYER_THICKNESS", h_max, &
700 call set_regrid_max_thickness(cs, h_max, gv%m_to_H)
701 elseif (index(trim(string),
'FNC1:')==1)
then
702 call dz_function1( trim(string(6:)), h_max )
703 call log_param(param_file, mdl,
"!MAX_LAYER_THICKNESS", h_max, &
705 call set_regrid_max_thickness(cs, h_max, gv%m_to_H)
707 call mom_error(fatal,trim(mdl)//
", initialize_regridding: "// &
708 "Unrecognized MAX_LAYER_THICKNESS_CONFIG "//trim(string))
713 if (
allocated(dz))
deallocate(dz)
714 end subroutine initialize_regridding
717 subroutine check_grid_def(filename, varname, expected_units, msg, ierr)
718 character(len=*),
intent(in) :: filename
719 character(len=*),
intent(in) :: varname
720 character(len=*),
intent(in) :: expected_units
721 character(len=*),
intent(inout) :: msg
722 logical,
intent(out) :: ierr
724 character (len=200) :: units, long_name
725 integer :: ncid, status, intid, vid
729 status = nf90_open(trim(filename), nf90_nowrite, ncid)
730 if (status /= nf90_noerr)
then
732 msg =
'File not found: '//trim(filename)
736 status = nf90_inq_varid(ncid, trim(varname), vid)
737 if (status /= nf90_noerr)
then
739 msg =
'Var not found: '//trim(varname)
743 status = nf90_get_att(ncid, vid,
"units", units)
744 if (status /= nf90_noerr)
then
746 msg =
'Attribute not found: units'
752 do i=1,len_trim(units)
753 if (units(i:i) == char(0)) units(i:i) =
" "
756 if (trim(units) /= trim(expected_units))
then
757 if (trim(expected_units) ==
"meters")
then
758 if (trim(units) /=
"m")
then
767 msg =
'Units incorrect: '//trim(units)//
' /= '//trim(expected_units)
770 end subroutine check_grid_def
773 subroutine end_regridding(CS)
776 if (
associated(cs%zlike_CS))
call end_coord_zlike(cs%zlike_CS)
777 if (
associated(cs%sigma_CS))
call end_coord_sigma(cs%sigma_CS)
778 if (
associated(cs%rho_CS))
call end_coord_rho(cs%rho_CS)
779 if (
associated(cs%hycom_CS))
call end_coord_hycom(cs%hycom_CS)
780 if (
associated(cs%slight_CS))
call end_coord_slight(cs%slight_CS)
781 if (
associated(cs%adapt_CS))
call end_coord_adapt(cs%adapt_CS)
783 deallocate( cs%coordinateResolution )
784 if (
allocated(cs%target_density))
deallocate( cs%target_density )
785 if (
allocated(cs%max_interface_depths) )
deallocate( cs%max_interface_depths )
786 if (
allocated(cs%max_layer_thickness) )
deallocate( cs%max_layer_thickness )
788 end subroutine end_regridding
792 subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, frac_shelf_h, conv_adjust)
813 type(ocean_grid_type),
intent(in) :: g
815 real,
dimension(SZI_(G),SZJ_(G), SZK_(GV)),
intent(inout) :: h
818 real,
dimension(SZI_(G),SZJ_(G), CS%nk),
intent(inout) :: h_new
819 real,
dimension(SZI_(G),SZJ_(G), CS%nk+1),
intent(inout) :: dzinterface
820 real,
dimension(:,:),
optional,
pointer :: frac_shelf_h
821 logical,
optional,
intent(in ) :: conv_adjust
823 real :: trickgnucompiler
824 logical :: use_ice_shelf
825 logical :: do_convective_adjustment
827 do_convective_adjustment = .true.
828 if (
present(conv_adjust)) do_convective_adjustment = conv_adjust
830 use_ice_shelf = .false.
831 if (
present(frac_shelf_h))
then
832 if (
associated(frac_shelf_h)) use_ice_shelf = .true.
835 select case ( cs%regridding_scheme )
837 case ( regridding_zstar )
838 if (use_ice_shelf)
then
839 call build_zstar_grid( cs, g, gv, h, dzinterface, frac_shelf_h )
841 call build_zstar_grid( cs, g, gv, h, dzinterface )
843 call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
845 case ( regridding_sigma_shelf_zstar)
846 call build_zstar_grid( cs, g, gv, h, dzinterface )
847 call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
849 case ( regridding_sigma )
850 call build_sigma_grid( cs, g, gv, h, dzinterface )
851 call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
853 case ( regridding_rho )
854 if (do_convective_adjustment)
call convective_adjustment(g, gv, h, tv)
855 call build_rho_grid( g, gv, h, tv, dzinterface, remapcs, cs )
856 call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
858 case ( regridding_arbitrary )
859 call build_grid_arbitrary( g, gv, h, dzinterface, trickgnucompiler, cs )
860 call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
862 case ( regridding_hycom1 )
863 call build_grid_hycom1( g, gv, h, tv, h_new, dzinterface, cs )
865 case ( regridding_slight )
866 call build_grid_slight( g, gv, h, tv, dzinterface, cs )
867 call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
869 case ( regridding_adaptive )
870 call build_grid_adaptive(g, gv, h, tv, dzinterface, remapcs, cs)
871 call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
874 call mom_error(fatal,
'MOM_regridding, regridding_main: '//&
875 'Unknown regridding scheme selected!')
879 #ifdef __DO_SAFETY_CHECKS__
880 call check_remapping_grid(g, gv, h, dzinterface,
'in regridding_main')
883 end subroutine regridding_main
886 subroutine calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new)
888 type(ocean_grid_type),
intent(in) :: G
890 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(in) :: h
891 real,
dimension(SZI_(G),SZJ_(G),CS%nk+1),
intent(in) :: dzInterface
892 real,
dimension(SZI_(G),SZJ_(G),CS%nk),
intent(inout) :: h_new
894 integer :: i, j, k, nki
896 nki = min(cs%nk, gv%ke)
899 do j = g%jsc-1,g%jec+1
900 do i = g%isc-1,g%iec+1
901 if (g%mask2dT(i,j)>0.)
then
903 h_new(i,j,k) = max( 0., h(i,j,k) + ( dzinterface(i,j,k) - dzinterface(i,j,k+1) ) )
905 if (cs%nk > gv%ke)
then
907 h_new(i,j,k) = max( 0., dzinterface(i,j,k) - dzinterface(i,j,k+1) )
911 h_new(i,j,1:nki) = h(i,j,1:nki)
912 if (cs%nk > gv%ke) h_new(i,j,nki+1:cs%nk) = 0.
918 end subroutine calc_h_new_by_dz
921 subroutine check_remapping_grid( G, GV, h, dzInterface, msg )
922 type(ocean_grid_type),
intent(in) :: g
924 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(in) :: h
925 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1),
intent(in) :: dzinterface
927 character(len=*),
intent(in) :: msg
932 do j = g%jsc-1,g%jec+1 ;
do i = g%isc-1,g%iec+1
933 if (g%mask2dT(i,j)>0.)
call check_grid_column( gv%ke, gv%Z_to_H*g%bathyT(i,j), h(i,j,:), dzinterface(i,j,:), msg )
936 end subroutine check_remapping_grid
939 subroutine check_grid_column( nk, depth, h, dzInterface, msg )
940 integer,
intent(in) :: nk
941 real,
intent(in) :: depth
942 real,
dimension(nk),
intent(in) :: h
943 real,
dimension(nk+1),
intent(in) :: dzinterface
944 character(len=*),
intent(in) :: msg
947 real :: eps, total_h_old, total_h_new, h_new, z_old, z_new
949 eps =1. ; eps = epsilon(eps)
954 total_h_old = total_h_old + h(k)
959 if (depth == 0.) z_old = - total_h_old
963 z_new = z_old + dzinterface(k)
964 h_new = h(k) + ( dzinterface(k) - dzinterface(k+1) )
966 write(0,*)
'k,h,hnew=',k,h(k),h_new
967 write(0,*)
'dzI(k+1),dzI(k)=',dzinterface(k+1),dzinterface(k)
968 call mom_error( fatal,
'MOM_regridding, check_grid_column: '//&
969 'Negative layer thickness implied by re-gridding, '//trim(msg))
971 total_h_new = total_h_new + h_new
976 if (abs(total_h_new-total_h_old)>real(nk-1)*0.5*(total_h_old+total_h_new)*eps)
then
979 write(0,*)
'k,h,hnew=',k,h(k),h(k)+(dzinterface(k)-dzinterface(k+1))
981 write(0,*)
'Hold,Hnew,Hnew-Hold=',total_h_old,total_h_new,total_h_new-total_h_old
982 write(0,*)
'eps,(n)/2*eps*H=',eps,real(nk-1)*0.5*(total_h_old+total_h_new)*eps
983 call mom_error( fatal,
'MOM_regridding, check_grid_column: '//&
984 'Re-gridding did NOT conserve total thickness to within roundoff '//trim(msg))
988 if (dzinterface(1) /= 0.)
call mom_error( fatal, &
989 'MOM_regridding, check_grid_column: Non-zero dzInterface at surface! '//trim(msg))
990 if (dzinterface(nk+1) /= 0.)
call mom_error( fatal, &
991 'MOM_regridding, check_grid_column: Non-zero dzInterface at bottom! '//trim(msg))
993 end subroutine check_grid_column
1000 subroutine filtered_grid_motion( CS, nk, z_old, z_new, dz_g )
1002 integer,
intent(in) :: nk
1003 real,
dimension(nk+1),
intent(in) :: z_old
1004 real,
dimension(CS%nk+1),
intent(in) :: z_new
1005 real,
dimension(CS%nk+1),
intent(inout) :: dz_g
1008 real :: dz_tgt, zr1, z_old_k
1009 real :: Aq, Bq, dz0, z0, F0
1010 real :: zs, zd, dzwt, Idzwt
1012 real :: Int_zs, Int_zd, dInt_zs_zd
1014 real,
dimension(nk+1) :: z_act
1016 logical :: debug = .false.
1019 if ((z_old(nk+1) - z_old(1)) * (z_new(cs%nk+1) - z_new(1)) < 0.0)
then
1020 call mom_error(fatal,
"filtered_grid_motion: z_old and z_new use different sign conventions.")
1021 elseif ((z_old(nk+1) - z_old(1)) * (z_new(cs%nk+1) - z_new(1)) == 0.0)
then
1023 do k=1,cs%nk+1 ; dz_g(k) = 0.0 ;
enddo ;
return
1024 elseif ((z_old(nk+1) - z_old(1)) + (z_new(cs%nk+1) - z_new(1)) > 0.0)
then
1032 if (sgn*(z_new(k)-z_new(k-1)) < -5e-16*(abs(z_new(k))+abs(z_new(k-1))) ) &
1033 call mom_error(fatal,
"filtered_grid_motion: z_new is tangled.")
1036 if (sgn*(z_old(k)-z_old(k-1)) < -5e-16*(abs(z_old(k))+abs(z_old(k-1))) ) &
1037 call mom_error(fatal,
"filtered_grid_motion: z_old is tangled.")
1042 zs = cs%depth_of_time_filter_shallow
1043 zd = cs%depth_of_time_filter_deep
1044 wtd = 1.0 - cs%old_grid_weight
1048 idzwt = 0.0 ;
if (abs(zd - zs) > 0.0) idzwt = 1.0 / (zd - zs)
1049 dint_zs_zd = 0.5*(1.0 + iwtd) * (zd - zs)
1050 aq = 0.5*(iwtd - 1.0)
1055 if (k<=nk+1) z_old_k = z_old(k)
1057 dz_tgt = sgn*(z_new(k) - z_old_k)
1058 zr1 = sgn*(z_old_k - z_old(1))
1062 if ((zr1 > zd) .and. (zr1 + wtd * dz_tgt > zd))
then
1063 dz_g(k) = sgn * wtd * dz_tgt
1064 elseif ((zr1 < zs) .and. (zr1 + dz_tgt < zs))
then
1065 dz_g(k) = sgn * dz_tgt
1074 int_zd = iwtd*(zd - zr1)
1075 int_zs = int_zd - dint_zs_zd
1076 elseif (zr1 <= zs)
then
1078 int_zd = dint_zs_zd + (zs - zr1)
1081 int_zd = (zd - zr1) * (iwtd*(0.5*(zd+zr1) - zs) + 0.5*(zd - zr1)) * idzwt
1082 int_zs = (zs - zr1) * (0.5*iwtd * ((zr1 - zs)) + (zd - 0.5*(zr1+zs))) * idzwt
1086 if (dz_tgt >= int_zd)
then
1087 dz_g(k) = sgn * ((zd-zr1) + wtd*(dz_tgt - int_zd))
1088 elseif (dz_tgt <= int_zs)
then
1089 dz_g(k) = sgn * ((zs-zr1) + (dz_tgt - int_zs))
1095 if (zr1 <= zs)
then ; dz0 = zs-zr1 ; z0 = zs ; f0 = dz_tgt - int_zs
1096 elseif (zr1 >= zd)
then ; dz0 = zd-zr1 ; z0 = zd ; f0 = dz_tgt - int_zd
1097 else ; dz0 = 0.0 ; z0 = zr1 ; f0 = dz_tgt ;
endif
1099 bq = (dzwt + 2.0*aq*(z0-zs))
1102 dz_g(k) = sgn * (dz0 + 2.0*f0*dzwt / (bq + sqrt(bq**2 + 4.0*aq*f0*dzwt) ))
1124 if (k<=nk+1) z_old_k = z_old(k)
1125 z_act(k) = z_old_k + dz_g(k)
1128 if (sgn*((z_act(k))-z_act(k-1)) < -1e-15*(abs(z_act(k))+abs(z_act(k-1))) ) &
1129 call mom_error(fatal,
"filtered_grid_motion: z_output is tangled.")
1133 end subroutine filtered_grid_motion
1138 subroutine build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h)
1142 type(ocean_grid_type),
intent(in) :: G
1144 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(in) :: h
1145 real,
dimension(SZI_(G),SZJ_(G), CS%nk+1),
intent(inout) :: dzInterface
1147 real,
dimension(:,:),
optional,
pointer :: frac_shelf_h
1151 real :: nominalDepth, totalThickness, dh
1152 real,
dimension(SZK_(GV)+1) :: zOld, zNew
1153 real :: minThickness
1154 logical :: ice_shelf
1157 minthickness = cs%min_thickness
1159 if (
present(frac_shelf_h))
then
1160 if (
associated(frac_shelf_h)) ice_shelf = .true.
1167 do j = g%jsc-1,g%jec+1
1168 do i = g%isc-1,g%iec+1
1170 if (g%mask2dT(i,j)==0.)
then
1171 dzinterface(i,j,:) = 0.
1176 nominaldepth = g%bathyT(i,j)*gv%Z_to_H
1179 totalthickness = 0.0
1181 totalthickness = totalthickness + h(i,j,k)
1184 zold(nz+1) = - nominaldepth
1186 zold(k) = zold(k+1) + h(i,j,k)
1190 if (frac_shelf_h(i,j) > 0.)
then
1191 call build_zstar_column(cs%zlike_CS, nominaldepth, totalthickness, znew, &
1192 z_rigid_top = totalthickness-nominaldepth, &
1193 eta_orig=zold(1), zscale=gv%Z_to_H)
1195 call build_zstar_column(cs%zlike_CS, nominaldepth, totalthickness, &
1196 znew, zscale=gv%Z_to_H)
1199 call build_zstar_column(cs%zlike_CS, nominaldepth, totalthickness, &
1200 znew, zscale=gv%Z_to_H)
1204 call filtered_grid_motion( cs, nz, zold, znew, dzinterface(i,j,:) )
1206 #ifdef __DO_SAFETY_CHECKS__
1207 dh=max(nominaldepth,totalthickness)
1208 if (abs(znew(1)-zold(1))>(nz-1)*0.5*epsilon(dh)*dh)
then
1209 write(0,*)
'min_thickness=',minthickness
1210 write(0,*)
'nominalDepth=',nominaldepth,
'totalThickness=',totalthickness
1211 write(0,*)
'dzInterface(1) = ',dzinterface(i,j,1),epsilon(dh),nz
1213 write(0,*) k,zold(k),znew(k)
1216 write(0,*) k,h(i,j,k),znew(k)-znew(k+1),cs%coordinateResolution(k)
1218 call mom_error( fatal, &
1219 'MOM_regridding, build_zstar_grid(): top surface has moved!!!' )
1223 call adjust_interface_motion( cs, nz, h(i,j,:), dzinterface(i,j,:) )
1228 end subroutine build_zstar_grid
1233 subroutine build_sigma_grid( CS, G, GV, h, dzInterface )
1243 type(ocean_grid_type),
intent(in) :: G
1245 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(in) :: h
1246 real,
dimension(SZI_(G),SZJ_(G), CS%nk+1),
intent(inout) :: dzInterface
1252 real :: nominalDepth, totalThickness, dh
1253 real,
dimension(SZK_(GV)+1) :: zOld, zNew
1257 do i = g%isc-1,g%iec+1
1258 do j = g%jsc-1,g%jec+1
1260 if (g%mask2dT(i,j)==0.)
then
1261 dzinterface(i,j,:) = 0.
1266 nominaldepth = g%bathyT(i,j)*gv%Z_to_H
1269 totalthickness = 0.0
1271 totalthickness = totalthickness + h(i,j,k)
1274 call build_sigma_column(cs%sigma_CS, nominaldepth, totalthickness, znew)
1277 zold(nz+1) = -nominaldepth
1279 zold(k) = zold(k+1) + h(i, j, k)
1282 call filtered_grid_motion( cs, nz, zold, znew, dzinterface(i,j,:) )
1284 #ifdef __DO_SAFETY_CHECKS__
1285 dh=max(nominaldepth,totalthickness)
1286 if (abs(znew(1)-zold(1))>(cs%nk-1)*0.5*epsilon(dh)*dh)
then
1287 write(0,*)
'min_thickness=',cs%min_thickness
1288 write(0,*)
'nominalDepth=',nominaldepth,
'totalThickness=',totalthickness
1289 write(0,*)
'dzInterface(1) = ',dzinterface(i,j,1),epsilon(dh),nz,cs%nk
1291 write(0,*) k,zold(k),znew(k)
1294 write(0,*) k,h(i,j,k),znew(k)-znew(k+1),totalthickness*cs%coordinateResolution(k),cs%coordinateResolution(k)
1296 call mom_error( fatal, &
1297 'MOM_regridding, build_sigma_grid: top surface has moved!!!' )
1299 dzinterface(i,j,1) = 0.
1300 dzinterface(i,j,cs%nk+1) = 0.
1306 end subroutine build_sigma_grid
1312 subroutine build_rho_grid( G, GV, h, tv, dzInterface, remapCS, CS )
1329 type(ocean_grid_type),
intent(in) :: G
1331 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(in) :: h
1333 real,
dimension(SZI_(G),SZJ_(G), SZK_(GV)+1),
intent(inout) :: dzInterface
1341 real :: nominalDepth, totalThickness
1342 real,
dimension(SZK_(GV)+1) :: zOld, zNew
1343 real :: h_neglect, h_neglect_edge
1344 #ifdef __DO_SAFETY_CHECKS__
1349 if (gv%Boussinesq)
then
1350 h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
1352 h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
1357 if (.not.cs%target_density_set)
call mom_error(fatal,
"build_rho_grid: "//&
1358 "Target densities must be set before build_rho_grid is called.")
1361 do j = g%jsc-1,g%jec+1
1362 do i = g%isc-1,g%iec+1
1364 if (g%mask2dT(i,j)==0.)
then
1365 dzinterface(i,j,:) = 0.
1371 nominaldepth = g%bathyT(i,j)*gv%Z_to_H
1373 call build_rho_column(cs%rho_CS, nz, nominaldepth, h(i, j, :), &
1374 tv%T(i, j, :), tv%S(i, j, :), tv%eqn_of_state, znew, &
1375 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
1377 if (cs%integrate_downward_for_e)
then
1380 zold(k+1) = zold(k) - h(i,j,k)
1384 zold(nz+1) = - nominaldepth
1386 zold(k) = zold(k+1) + h(i,j,k)
1391 call filtered_grid_motion( cs, nz, zold, znew, dzinterface(i,j,:) )
1393 #ifdef __DO_SAFETY_CHECKS__
1395 if (znew(k) > zold(1))
then
1396 write(0,*)
'zOld=',zold
1397 write(0,*)
'zNew=',znew
1398 call mom_error( fatal,
'MOM_regridding, build_rho_grid: '//&
1399 'interior interface above surface!' )
1401 if (znew(k) > znew(k-1))
then
1402 write(0,*)
'zOld=',zold
1403 write(0,*)
'zNew=',znew
1404 call mom_error( fatal,
'MOM_regridding, build_rho_grid: '//&
1405 'interior interfaces cross!' )
1409 totalthickness = 0.0
1411 totalthickness = totalthickness + h(i,j,k)
1414 dh=max(nominaldepth,totalthickness)
1415 if (abs(znew(1)-zold(1))>(nz-1)*0.5*epsilon(dh)*dh)
then
1416 write(0,*)
'min_thickness=',cs%min_thickness
1417 write(0,*)
'nominalDepth=',nominaldepth,
'totalThickness=',totalthickness
1418 write(0,*)
'zNew(1)-zOld(1) = ',znew(1)-zold(1),epsilon(dh),nz
1420 write(0,*) k,zold(k),znew(k)
1423 write(0,*) k,h(i,j,k),znew(k)-znew(k+1)
1425 call mom_error( fatal, &
1426 'MOM_regridding, build_rho_grid: top surface has moved!!!' )
1433 end subroutine build_rho_grid
1442 subroutine build_grid_hycom1( G, GV, h, tv, h_new, dzInterface, CS )
1443 type(ocean_grid_type),
intent(in) :: G
1445 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(in) :: h
1448 real,
dimension(SZI_(G),SZJ_(G),CS%nk),
intent(inout) :: h_new
1449 real,
dimension(SZI_(G),SZJ_(G),CS%nk+1),
intent(inout) :: dzInterface
1452 real,
dimension(SZK_(GV)+1) :: z_col
1453 real,
dimension(CS%nk+1) :: z_col_new
1454 real,
dimension(SZK_(GV)+1) :: dz_col
1455 real,
dimension(SZK_(GV)) :: p_col
1456 integer :: i, j, k, nki
1458 real :: h_neglect, h_neglect_edge
1461 if (gv%Boussinesq)
then
1462 h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
1464 h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
1467 if (.not.cs%target_density_set)
call mom_error(fatal,
"build_grid_HyCOM1 : "//&
1468 "Target densities must be set before build_grid_HyCOM1 is called.")
1470 nki = min(gv%ke, cs%nk)
1473 do j = g%jsc-1,g%jec+1 ;
do i = g%isc-1,g%iec+1
1474 if (g%mask2dT(i,j)>0.)
then
1476 depth = g%bathyT(i,j) * gv%Z_to_H
1480 z_col(k+1) = z_col(k) + h(i,j,k)
1481 p_col(k) = cs%ref_pressure + cs%compressibility_fraction * &
1482 ( 0.5 * ( z_col(k) + z_col(k+1) ) * gv%H_to_Pa - cs%ref_pressure )
1485 call build_hycom1_column(cs%hycom_CS, tv%eqn_of_state, gv%ke, depth, &
1486 h(i, j, :), tv%T(i, j, :), tv%S(i, j, :), p_col, &
1487 z_col, z_col_new, zscale=gv%Z_to_H, &
1488 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
1491 call filtered_grid_motion( cs, gv%ke, z_col, z_col_new, dz_col )
1494 dz_col(:) = -dz_col(:)
1495 call adjust_interface_motion( cs, gv%ke, h(i,j,:), dz_col(:) )
1497 dzinterface(i,j,1:nki+1) = dz_col(1:nki+1)
1498 if (nki<cs%nk) dzinterface(i,j,nki+2:cs%nk+1) = 0.
1501 dzinterface(i,j,:) = 0.
1505 call calc_h_new_by_dz(cs, g, gv, h, dzinterface, h_new)
1507 end subroutine build_grid_hycom1
1511 subroutine build_grid_adaptive(G, GV, h, tv, dzInterface, remapCS, CS)
1512 type(ocean_grid_type),
intent(in) :: G
1514 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(in) :: h
1517 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1),
intent(inout) :: dzInterface
1523 integer :: i, j, k, nz
1525 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: tInt, sInt
1528 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: zInt
1529 real,
dimension(SZK_(GV)+1) :: zNext
1537 do k = 2, nz ;
do j = g%jsc-2,g%jec+2 ;
do i = g%isc-2,g%iec+2
1538 tint(i,j,k) = 0.5 * (tv%T(i,j,k-1) + tv%T(i,j,k))
1539 sint(i,j,k) = 0.5 * (tv%S(i,j,k-1) + tv%S(i,j,k))
1540 zint(i,j,k) = zint(i,j,k-1) + h(i,j,k-1)
1541 enddo ;
enddo ;
enddo
1545 tint(:,:,1) = tv%T(:,:,1) ; tint(:,:,nz+1) = tv%T(:,:,nz)
1546 sint(:,:,1) = tv%S(:,:,1) ; sint(:,:,nz+1) = tv%S(:,:,nz)
1549 zint(:,:,nz+1) = zint(:,:,nz) + h(:,:,nz)
1553 do j = g%jsc-1,g%jec+1 ;
do i = g%isc-1,g%iec+1
1554 if (g%mask2dT(i,j) < 0.5)
then
1555 dzinterface(i,j,:) = 0.
1559 call build_adapt_column(cs%adapt_CS, g, gv, tv, i, j, zint, tint, sint, h, znext)
1561 call filtered_grid_motion(cs, nz, zint(i,j,:), znext, dzinterface(i,j,:))
1563 do k = 1, nz+1 ; dzinterface(i,j,k) = -dzinterface(i,j,k) ;
enddo
1564 call adjust_interface_motion(cs, nz, h(i,j,:), dzinterface(i,j,:))
1566 end subroutine build_grid_adaptive
1577 subroutine build_grid_slight(G, GV, h, tv, dzInterface, CS)
1578 type(ocean_grid_type),
intent(in) :: G
1580 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(in) :: h
1582 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1),
intent(inout) :: dzInterface
1585 real,
dimension(SZK_(GV)+1) :: z_col
1586 real,
dimension(SZK_(GV)+1) :: z_col_new
1587 real,
dimension(SZK_(GV)+1) :: dz_col
1588 real,
dimension(SZK_(GV)) :: p_col
1590 integer :: i, j, k, nz
1591 real :: h_neglect, h_neglect_edge
1594 if (gv%Boussinesq)
then
1595 h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
1597 h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
1602 if (.not.cs%target_density_set)
call mom_error(fatal,
"build_grid_SLight : "//&
1603 "Target densities must be set before build_grid_SLight is called.")
1606 do j = g%jsc-1,g%jec+1 ;
do i = g%isc-1,g%iec+1
1607 if (g%mask2dT(i,j)>0.)
then
1609 depth = g%bathyT(i,j) * gv%Z_to_H
1612 z_col(k+1) = z_col(k) + h(i,j,k)
1613 p_col(k) = cs%ref_pressure + cs%compressibility_fraction * &
1614 ( 0.5 * ( z_col(k) + z_col(k+1) ) * gv%H_to_Pa - cs%ref_pressure )
1617 call build_slight_column(cs%slight_CS, tv%eqn_of_state, gv%H_to_Pa, &
1618 gv%H_subroundoff, nz, depth, h(i, j, :), &
1619 tv%T(i, j, :), tv%S(i, j, :), p_col, z_col, z_col_new, &
1620 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
1623 call filtered_grid_motion( cs, nz, z_col, z_col_new, dz_col )
1624 do k=1,nz+1 ; dzinterface(i,j,k) = -dz_col(k) ;
enddo
1625 #ifdef __DO_SAFETY_CHECKS__
1626 if (dzinterface(i,j,1) /= 0.) stop
'build_grid_SLight: Surface moved?!'
1627 if (dzinterface(i,j,nz+1) /= 0.) stop
'build_grid_SLight: Bottom moved?!'
1631 call adjust_interface_motion( cs, nz, h(i,j,:), dzinterface(i,j,:) )
1634 dzinterface(i,j,:) = 0.
1638 end subroutine build_grid_slight
1641 subroutine adjust_interface_motion( CS, nk, h_old, dz_int )
1643 integer,
intent(in) :: nk
1644 real,
dimension(nk),
intent(in) :: h_old
1645 real,
dimension(CS%nk+1),
intent(inout) :: dz_int
1648 real :: h_new, eps, h_total, h_err
1650 eps = 1. ; eps = epsilon(eps)
1652 h_total = 0. ; h_err = 0.
1653 do k = 1, min(cs%nk,nk)
1654 h_total = h_total + h_old(k)
1655 h_err = h_err + max( h_old(k), abs(dz_int(k)), abs(dz_int(k+1)) )*eps
1656 h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
1657 if (h_new < -3.0*h_err)
then
1658 write(0,*)
'h<0 at k=',k,
'h_old=',h_old(k), &
1659 'wup=',dz_int(k),
'wdn=',dz_int(k+1),
'dw_dz=',dz_int(k) - dz_int(k+1), &
1660 'h_new=',h_new,
'h_err=',h_err
1661 call mom_error( fatal,
'MOM_regridding: adjust_interface_motion() - '//&
1662 'implied h<0 is larger than roundoff!')
1667 h_err = h_err + max( abs(dz_int(k)), abs(dz_int(k+1)) )*eps
1668 h_new = ( dz_int(k) - dz_int(k+1) )
1669 if (h_new < -3.0*h_err)
then
1670 write(0,*)
'h<0 at k=',k,
'h_old was empty',&
1671 'wup=',dz_int(k),
'wdn=',dz_int(k+1),
'dw_dz=',dz_int(k) - dz_int(k+1), &
1672 'h_new=',h_new,
'h_err=',h_err
1673 call mom_error( fatal,
'MOM_regridding: adjust_interface_motion() - '//&
1674 'implied h<0 is larger than roundoff!')
1678 do k = min(cs%nk,nk),2,-1
1679 h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
1680 if (h_new<cs%min_thickness) &
1681 dz_int(k) = ( dz_int(k+1) - h_old(k) ) + cs%min_thickness
1682 h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
1684 dz_int(k) = ( 1. - eps ) * ( dz_int(k+1) - h_old(k) )
1685 h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
1687 write(0,*)
'h<0 at k=',k,
'h_old=',h_old(k), &
1688 'wup=',dz_int(k),
'wdn=',dz_int(k+1),
'dw_dz=',dz_int(k) - dz_int(k+1), &
1690 stop
'Still did not work!'
1691 call mom_error( fatal,
'MOM_regridding: adjust_interface_motion() - '//&
1692 'Repeated adjustment for roundoff h<0 failed!')
1697 end subroutine adjust_interface_motion
1702 subroutine build_grid_arbitrary( G, GV, h, dzInterface, h_new, CS )
1708 type(ocean_grid_type),
intent(in) :: G
1710 real,
dimension(SZI_(G),SZJ_(G), SZK_(GV)),
intent(in) :: h
1711 real,
dimension(SZI_(G),SZJ_(G), SZK_(GV)+1),
intent(inout) :: dzInterface
1713 real,
intent(inout) :: h_new
1719 real :: z_inter(SZK_(GV)+1)
1720 real :: total_height
1725 real :: x1, y1, x2, y2
1729 max_depth = g%max_depth*gv%Z_to_H
1731 do j = g%jsc-1,g%jec+1
1732 do i = g%isc-1,g%iec+1
1735 local_depth = g%bathyT(i,j)*gv%Z_to_H
1740 total_height = total_height + h(i,j,k)
1743 eta = total_height - local_depth
1746 delta_h = (max_depth + eta) / nz
1751 z_inter(k+1) = z_inter(k) - delta_h
1756 x1 = 0.35; y1 = 0.45; x2 = 0.65; y2 = 0.55
1758 x = - ( z_inter(k) - eta ) / max_depth
1762 elseif ( (x > x1 ) .and. ( x < x2 ))
then
1763 t = y1 + (y2-y1) * (x-x1) / (x2-x1)
1765 t = y2 + (1.0-y2) * (x-x2) / (1.0-x2)
1768 z_inter(k) = -t * max_depth + eta
1773 z_inter(nz+1) = - local_depth
1777 if ( z_inter(k) < (z_inter(k+1) + cs%min_thickness) )
then
1778 z_inter(k) = z_inter(k+1) + cs%min_thickness
1784 dzinterface(i,j,1) = 0.
1787 dzinterface(i,j,k) = z_inter(k) - x
1789 dzinterface(i,j,nz+1) = 0.
1797 end subroutine build_grid_arbitrary
1804 subroutine inflate_vanished_layers_old( CS, G, GV, h )
1815 type(ocean_grid_type),
intent(in) :: g
1817 real,
dimension(SZI_(G),SZJ_(G), SZK_(GV)),
intent(inout) :: h
1823 do i = g%isc-1,g%iec+1
1824 do j = g%jsc-1,g%jec+1
1831 call old_inflate_layers_1d( cs%min_thickness, gv%ke, htmp )
1841 end subroutine inflate_vanished_layers_old
1845 subroutine convective_adjustment(G, GV, h, tv)
1846 type(ocean_grid_type),
intent(in) :: G
1848 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1862 logical :: stratified
1863 real,
dimension(GV%ke) :: p_col, densities
1868 do j = g%jsc-1,g%jec+1 ;
do i = g%isc-1,g%iec+1
1872 densities, 1, gv%ke, tv%eqn_of_state )
1879 t0 = tv%T(i,j,k) ; t1 = tv%T(i,j,k+1)
1880 s0 = tv%S(i,j,k) ; s1 = tv%S(i,j,k+1)
1881 r0 = densities(k) ; r1 = densities(k+1)
1882 h0 = h(i,j,k) ; h1 = h(i,j,k+1)
1887 tv%T(i,j,k) = t1 ; tv%T(i,j,k+1) = t0
1888 tv%S(i,j,k) = s1 ; tv%S(i,j,k+1) = s0
1889 h(i,j,k) = h1 ; h(i,j,k+1) = h0
1892 densities(k), tv%eqn_of_state )
1894 densities(k+1), tv%eqn_of_state )
1895 stratified = .false.
1899 if ( stratified )
exit
1904 end subroutine convective_adjustment
1909 function uniformresolution(nk,coordMode,maxDepth,rhoLight,rhoHeavy)
1914 integer,
intent(in) :: nk
1915 character(len=*),
intent(in) :: coordmode
1918 real,
intent(in) :: maxdepth
1919 real,
intent(in) :: rholight
1920 real,
intent(in) :: rhoheavy
1922 real :: uniformresolution(nk)
1927 scheme = coordinatemode(coordmode)
1928 select case ( scheme )
1930 case ( regridding_zstar, regridding_hycom1, regridding_slight, regridding_sigma_shelf_zstar, &
1931 regridding_adaptive )
1932 uniformresolution(:) = maxdepth / real(nk)
1934 case ( regridding_rho )
1935 uniformresolution(:) = (rhoheavy - rholight) / real(nk)
1937 case ( regridding_sigma )
1938 uniformresolution(:) = 1. / real(nk)
1941 call mom_error(fatal,
"MOM_regridding, uniformResolution: "//&
1942 "Unrecognized choice for coordinate mode ("//trim(coordmode)//
").")
1946 end function uniformresolution
1950 subroutine initcoord(CS, GV, coord_mode)
1952 character(len=*),
intent(in) :: coord_mode
1957 select case (coordinatemode(coord_mode))
1958 case (regridding_zstar)
1959 call init_coord_zlike(cs%zlike_CS, cs%nk, cs%coordinateResolution)
1960 case (regridding_sigma_shelf_zstar)
1961 call init_coord_zlike(cs%zlike_CS, cs%nk, cs%coordinateResolution)
1962 case (regridding_sigma)
1963 call init_coord_sigma(cs%sigma_CS, cs%nk, cs%coordinateResolution)
1964 case (regridding_rho)
1965 call init_coord_rho(cs%rho_CS, cs%nk, cs%ref_pressure, cs%target_density, cs%interp_CS)
1966 case (regridding_hycom1)
1967 call init_coord_hycom(cs%hycom_CS, cs%nk, cs%coordinateResolution, cs%target_density, cs%interp_CS)
1968 case (regridding_slight)
1969 call init_coord_slight(cs%slight_CS, cs%nk, cs%ref_pressure, cs%target_density, cs%interp_CS, gv%m_to_H)
1970 case (regridding_adaptive)
1971 call init_coord_adapt(cs%adapt_CS, cs%nk, cs%coordinateResolution, gv%m_to_H)
1973 end subroutine initcoord
1977 subroutine setcoordinateresolution( dz, CS, scale )
1978 real,
dimension(:),
intent(in) :: dz
1980 real,
optional,
intent(in) :: scale
1982 if (
size(dz)/=cs%nk)
call mom_error( fatal, &
1983 'setCoordinateResolution: inconsistent number of levels' )
1985 if (
present(scale))
then
1986 cs%coordinateResolution(:) = scale*dz(:)
1988 cs%coordinateResolution(:) = dz(:)
1991 end subroutine setcoordinateresolution
1994 subroutine set_target_densities_from_gv( GV, CS )
2001 cs%target_density(1) = gv%Rlay(1)+0.5*(gv%Rlay(1)-gv%Rlay(2))
2002 cs%target_density(nz+1) = gv%Rlay(nz)+0.5*(gv%Rlay(nz)-gv%Rlay(nz-1))
2004 cs%target_density(k) = cs%target_density(k-1) + cs%coordinateResolution(k)
2006 cs%target_density_set = .true.
2008 end subroutine set_target_densities_from_gv
2011 subroutine set_target_densities( CS, rho_int )
2013 real,
dimension(CS%nk+1),
intent(in) :: rho_int
2015 if (
size(cs%target_density)/=
size(rho_int))
then
2016 call mom_error(fatal,
"set_target_densities inconsistent args!")
2019 cs%target_density(:) = rho_int(:)
2020 cs%target_density_set = .true.
2022 end subroutine set_target_densities
2025 subroutine set_regrid_max_depths( CS, max_depths, units_to_H )
2027 real,
dimension(CS%nk+1),
intent(in) :: max_depths
2028 real,
optional,
intent(in) :: units_to_h
2033 if (.not.
allocated(cs%max_interface_depths))
allocate(cs%max_interface_depths(1:cs%nk+1))
2035 val_to_h = 1.0 ;
if (
present(units_to_h)) val_to_h = units_to_h
2036 if (max_depths(cs%nk+1) < max_depths(1)) val_to_h = -1.0*val_to_h
2039 if (max_depths(cs%nk+1) < max_depths(1))
then
2040 do k=1,cs%nk ;
if (max_depths(k+1) > max_depths(k)) &
2041 call mom_error(fatal,
"Unordered list of maximum depths sent to set_regrid_max_depths!")
2044 do k=1,cs%nk ;
if (max_depths(k+1) < max_depths(k)) &
2045 call mom_error(fatal,
"Unordered list of maximum depths sent to set_regrid_max_depths.")
2050 cs%max_interface_depths(k) = val_to_h * max_depths(k)
2054 select case (cs%regridding_scheme)
2055 case (regridding_hycom1)
2056 call set_hycom_params(cs%hycom_CS, max_interface_depths=cs%max_interface_depths)
2057 case (regridding_slight)
2058 call set_slight_params(cs%slight_CS, max_interface_depths=cs%max_interface_depths)
2060 end subroutine set_regrid_max_depths
2063 subroutine set_regrid_max_thickness( CS, max_h, units_to_H )
2065 real,
dimension(CS%nk+1),
intent(in) :: max_h
2066 real,
optional,
intent(in) :: units_to_h
2071 if (.not.
allocated(cs%max_layer_thickness))
allocate(cs%max_layer_thickness(1:cs%nk))
2073 val_to_h = 1.0 ;
if (
present( units_to_h)) val_to_h = units_to_h
2076 cs%max_layer_thickness(k) = val_to_h * max_h(k)
2080 select case (cs%regridding_scheme)
2081 case (regridding_hycom1)
2082 call set_hycom_params(cs%hycom_CS, max_layer_thickness=cs%max_layer_thickness)
2083 case (regridding_slight)
2084 call set_slight_params(cs%slight_CS, max_layer_thickness=cs%max_layer_thickness)
2086 end subroutine set_regrid_max_thickness
2091 function getcoordinateresolution( CS, undo_scaling )
2093 logical,
optional,
intent(in) :: undo_scaling
2095 real,
dimension(CS%nk) :: getcoordinateresolution
2098 unscale = .false. ;
if (
present(undo_scaling)) unscale = undo_scaling
2101 getcoordinateresolution(:) = cs%coord_scale * cs%coordinateResolution(:)
2103 getcoordinateresolution(:) = cs%coordinateResolution(:)
2106 end function getcoordinateresolution
2109 function getcoordinateinterfaces( CS, undo_scaling )
2111 logical,
optional,
intent(in) :: undo_scaling
2113 real,
dimension(CS%nk+1) :: getcoordinateinterfaces
2117 unscale = .false. ;
if (
present(undo_scaling)) unscale = undo_scaling
2121 if (cs%regridding_scheme == regridding_rho)
then
2122 if (.not. cs%target_density_set) &
2123 call mom_error(fatal,
'MOM_regridding, getCoordinateInterfaces: '//&
2124 'target densities not set!')
2126 getcoordinateinterfaces(:) = cs%target_density(:)
2129 getcoordinateinterfaces(1) = 0.
2131 getcoordinateinterfaces(k+1) = getcoordinateinterfaces(k) - &
2132 cs%coord_scale * cs%coordinateResolution(k)
2135 getcoordinateinterfaces(1) = 0.
2137 getcoordinateinterfaces(k+1) = getcoordinateinterfaces(k) - &
2138 cs%coordinateResolution(k)
2143 getcoordinateinterfaces(:) = abs( getcoordinateinterfaces(:) )
2146 end function getcoordinateinterfaces
2150 function getcoordinateunits( CS )
2152 character(len=20) :: getcoordinateunits
2154 select case ( cs%regridding_scheme )
2155 case ( regridding_zstar, regridding_hycom1, regridding_slight, regridding_adaptive )
2156 getcoordinateunits =
'meter'
2157 case ( regridding_sigma_shelf_zstar )
2158 getcoordinateunits =
'meter/fraction'
2159 case ( regridding_sigma )
2160 getcoordinateunits =
'fraction'
2161 case ( regridding_rho )
2162 getcoordinateunits =
'kg/m3'
2163 case ( regridding_arbitrary )
2164 getcoordinateunits =
'unknown'
2166 call mom_error(fatal,
'MOM_regridding, getCoordinateUnits: '//&
2167 'Unknown regridding scheme selected!')
2170 end function getcoordinateunits
2174 function getcoordinateshortname( CS )
2176 character(len=20) :: getcoordinateshortname
2178 select case ( cs%regridding_scheme )
2179 case ( regridding_zstar )
2182 getcoordinateshortname =
'pseudo-depth, -z*'
2183 case ( regridding_sigma_shelf_zstar )
2184 getcoordinateshortname =
'pseudo-depth, -z*/sigma'
2185 case ( regridding_sigma )
2186 getcoordinateshortname =
'sigma'
2187 case ( regridding_rho )
2188 getcoordinateshortname =
'rho'
2189 case ( regridding_arbitrary )
2190 getcoordinateshortname =
'coordinate'
2191 case ( regridding_hycom1 )
2192 getcoordinateshortname =
'z-rho'
2193 case ( regridding_slight )
2194 getcoordinateshortname =
's-rho'
2195 case ( regridding_adaptive )
2196 getcoordinateshortname =
'adaptive'
2198 call mom_error(fatal,
'MOM_regridding, getCoordinateShortName: '//&
2199 'Unknown regridding scheme selected!')
2202 end function getcoordinateshortname
2205 subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_grid_weight, &
2206 interp_scheme, depth_of_time_filter_shallow, depth_of_time_filter_deep, &
2207 compress_fraction, dz_min_surface, nz_fixed_surface, Rho_ML_avg_depth, &
2208 nlay_ML_to_interior, fix_haloclines, halocline_filt_len, &
2209 halocline_strat_tol, integrate_downward_for_e, &
2210 adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, adaptAlpha, adaptDoMin)
2212 logical,
optional,
intent(in) :: boundary_extrapolation
2213 real,
optional,
intent(in) :: min_thickness
2215 real,
optional,
intent(in) :: old_grid_weight
2216 character(len=*),
optional,
intent(in) :: interp_scheme
2217 real,
optional,
intent(in) :: depth_of_time_filter_shallow
2218 real,
optional,
intent(in) :: depth_of_time_filter_deep
2219 real,
optional,
intent(in) :: compress_fraction
2220 real,
optional,
intent(in) :: dz_min_surface
2222 integer,
optional,
intent(in) :: nz_fixed_surface
2223 real,
optional,
intent(in) :: rho_ml_avg_depth
2225 real,
optional,
intent(in) :: nlay_ml_to_interior
2227 logical,
optional,
intent(in) :: fix_haloclines
2228 real,
optional,
intent(in) :: halocline_filt_len
2230 real,
optional,
intent(in) :: halocline_strat_tol
2232 logical,
optional,
intent(in) :: integrate_downward_for_e
2234 real,
optional,
intent(in) :: adapttimeratio
2235 real,
optional,
intent(in) :: adaptzoom
2236 real,
optional,
intent(in) :: adaptzoomcoeff
2237 real,
optional,
intent(in) :: adaptbuoycoeff
2238 real,
optional,
intent(in) :: adaptalpha
2239 logical,
optional,
intent(in) :: adaptdomin
2243 if (
present(interp_scheme))
call set_interp_scheme(cs%interp_CS, interp_scheme)
2244 if (
present(boundary_extrapolation))
call set_interp_extrap(cs%interp_CS, boundary_extrapolation)
2246 if (
present(old_grid_weight))
then
2247 if (old_grid_weight<0. .or. old_grid_weight>1.) &
2248 call mom_error(fatal,
'MOM_regridding, set_regrid_params: Weight is out side the range 0..1!')
2249 cs%old_grid_weight = old_grid_weight
2251 if (
present(depth_of_time_filter_shallow)) cs%depth_of_time_filter_shallow = depth_of_time_filter_shallow
2252 if (
present(depth_of_time_filter_deep)) cs%depth_of_time_filter_deep = depth_of_time_filter_deep
2253 if (
present(depth_of_time_filter_shallow) .or.
present(depth_of_time_filter_deep))
then
2254 if (cs%depth_of_time_filter_deep<cs%depth_of_time_filter_shallow)
call mom_error(fatal,
'MOM_regridding, '//&
2255 'set_regrid_params: depth_of_time_filter_deep<depth_of_time_filter_shallow!')
2258 if (
present(min_thickness)) cs%min_thickness = min_thickness
2259 if (
present(compress_fraction)) cs%compressibility_fraction = compress_fraction
2260 if (
present(integrate_downward_for_e)) cs%integrate_downward_for_e = integrate_downward_for_e
2262 select case (cs%regridding_scheme)
2263 case (regridding_zstar)
2264 if (
present(min_thickness))
call set_zlike_params(cs%zlike_CS, min_thickness=min_thickness)
2265 case (regridding_sigma_shelf_zstar)
2266 if (
present(min_thickness))
call set_zlike_params(cs%zlike_CS, min_thickness=min_thickness)
2267 case (regridding_sigma)
2268 if (
present(min_thickness))
call set_sigma_params(cs%sigma_CS, min_thickness=min_thickness)
2269 case (regridding_rho)
2270 if (
present(min_thickness))
call set_rho_params(cs%rho_CS, min_thickness=min_thickness)
2271 if (
present(integrate_downward_for_e)) &
2272 call set_rho_params(cs%rho_CS, integrate_downward_for_e=integrate_downward_for_e)
2273 if (
associated(cs%rho_CS) .and. (
present(interp_scheme) .or.
present(boundary_extrapolation))) &
2274 call set_rho_params(cs%rho_CS, interp_cs=cs%interp_CS)
2275 case (regridding_hycom1)
2276 if (
associated(cs%hycom_CS) .and. (
present(interp_scheme) .or.
present(boundary_extrapolation))) &
2277 call set_hycom_params(cs%hycom_CS, interp_cs=cs%interp_CS)
2278 case (regridding_slight)
2279 if (
present(min_thickness))
call set_slight_params(cs%slight_CS, min_thickness=min_thickness)
2280 if (
present(dz_min_surface))
call set_slight_params(cs%slight_CS, dz_ml_min=dz_min_surface)
2281 if (
present(nz_fixed_surface))
call set_slight_params(cs%slight_CS, nz_fixed_surface=nz_fixed_surface)
2282 if (
present(rho_ml_avg_depth))
call set_slight_params(cs%slight_CS, rho_ml_avg_depth=rho_ml_avg_depth)
2283 if (
present(nlay_ml_to_interior))
call set_slight_params(cs%slight_CS, nlay_ml_offset=nlay_ml_to_interior)
2284 if (
present(fix_haloclines))
call set_slight_params(cs%slight_CS, fix_haloclines=fix_haloclines)
2285 if (
present(halocline_filt_len))
call set_slight_params(cs%slight_CS, halocline_filter_length=halocline_filt_len)
2286 if (
present(halocline_strat_tol))
call set_slight_params(cs%slight_CS, halocline_strat_tol=halocline_strat_tol)
2287 if (
present(compress_fraction))
call set_slight_params(cs%slight_CS, compressibility_fraction=compress_fraction)
2288 if (
associated(cs%slight_CS) .and. (
present(interp_scheme) .or.
present(boundary_extrapolation))) &
2289 call set_slight_params(cs%slight_CS, interp_cs=cs%interp_CS)
2290 case (regridding_adaptive)
2291 if (
present(adapttimeratio))
call set_adapt_params(cs%adapt_CS, adapttimeratio=adapttimeratio)
2292 if (
present(adaptzoom))
call set_adapt_params(cs%adapt_CS, adaptzoom=adaptzoom)
2293 if (
present(adaptzoomcoeff))
call set_adapt_params(cs%adapt_CS, adaptzoomcoeff=adaptzoomcoeff)
2294 if (
present(adaptbuoycoeff))
call set_adapt_params(cs%adapt_CS, adaptbuoycoeff=adaptbuoycoeff)
2295 if (
present(adaptalpha))
call set_adapt_params(cs%adapt_CS, adaptalpha=adaptalpha)
2296 if (
present(adaptdomin))
call set_adapt_params(cs%adapt_CS, adaptdomin=adaptdomin)
2299 end subroutine set_regrid_params
2302 integer function get_regrid_size(CS)
2305 get_regrid_size = cs%nk
2307 end function get_regrid_size
2310 function get_zlike_cs(CS)
2314 get_zlike_cs = cs%zlike_CS
2315 end function get_zlike_cs
2318 function get_sigma_cs(CS)
2322 get_sigma_cs = cs%sigma_CS
2323 end function get_sigma_cs
2326 function get_rho_cs(CS)
2328 type(
rho_cs) :: get_rho_cs
2330 get_rho_cs = cs%rho_CS
2331 end function get_rho_cs
2335 function getstaticthickness( CS, SSH, depth )
2337 real,
intent(in) :: ssh
2338 real,
intent(in) :: depth
2339 real,
dimension(CS%nk) :: getstaticthickness
2344 select case ( cs%regridding_scheme )
2345 case ( regridding_zstar, regridding_sigma_shelf_zstar, regridding_hycom1, regridding_slight, regridding_adaptive )
2349 dz = cs%coordinateResolution(k) * ( 1. + ssh/depth )
2351 dz = min(dz, depth - z)
2353 getstaticthickness(k) = dz
2356 getstaticthickness(:) = 0.
2358 case ( regridding_sigma )
2359 getstaticthickness(:) = cs%coordinateResolution(:) * ( depth + ssh )
2360 case ( regridding_rho )
2361 getstaticthickness(:) = 0.
2362 case ( regridding_arbitrary )
2363 getstaticthickness(:) = 0.
2365 call mom_error(fatal,
'MOM_regridding, getStaticThickness: '//&
2366 'Unknown regridding scheme selected!')
2369 end function getstaticthickness
2372 subroutine dz_function1( string, dz )
2373 character(len=*),
intent(in) :: string
2375 real,
dimension(:),
intent(inout) :: dz
2378 real :: dz_min, power, prec, H_total
2382 read( string, *) dz_min, h_total, power, prec
2383 if (prec == -1024.)
call mom_error(fatal,
"dz_function1: "// &
2384 "Problem reading FNC1: string ="//trim(string))
2387 dz(k) = (real(k-1)/real(nk-1))**power
2389 dz(:) = ( h_total - real(nk) * dz_min ) * ( dz(:) / sum(dz) )
2390 dz(:) = anint( dz(:) / prec ) * prec
2391 dz(:) = ( h_total - real(nk) * dz_min ) * ( dz(:) / sum(dz) )
2392 dz(:) = anint( dz(:) / prec ) * prec
2393 dz(nk) = dz(nk) + ( h_total - sum( dz(:) + dz_min ) )
2394 dz(:) = anint( dz(:) / prec ) * prec
2395 dz(:) = dz(:) + dz_min
2397 end subroutine dz_function1
2401 integer function rho_function1( string, rho_target )
2402 character(len=*),
intent(in) :: string
2404 real,
dimension(:),
allocatable,
intent(inout) :: rho_target
2406 integer :: nki, k, nk
2407 real :: ddx, dx, rho_1, rho_2, rho_3, drho, rho_4, drho_min
2409 read( string, *) nk, rho_1, rho_2, rho_3, drho, rho_4, drho_min
2410 allocate(rho_target(nk+1))
2412 rho_target(1) = rho_1
2413 rho_target(2) = rho_2
2416 ddx = max( drho_min, real(nki-k)/real(nki*nki) )
2418 rho_target(3+k) = rho_3 + (2. * drho) * dx
2420 rho_target(nki+4) = rho_4
2424 end function rho_function1