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.
132 public build_rho_column
137 public default_coordinate_mode
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"
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)"
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"
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: "//&
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 "//&
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
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 us%R_to_kg_m3*(gv%Rlay(1) + 0.5*(gv%Rlay(1)-gv%Rlay(min(2,ke)))), &
315 us%R_to_kg_m3*(gv%Rlay(ke) + 0.5*(gv%Rlay(ke)-gv%Rlay(max(ke-1,1)))) )
316 if (main_parameters)
call log_param(param_file, mdl,
"!"//coord_res_param, dz, &
317 trim(message), units=trim(coord_units))
318 elseif (trim(string)==
'PARAM')
then
322 call get_param(param_file, mdl, coord_res_param, dz, &
323 trim(message), units=trim(coord_units), fail_if_missing=.true.)
324 elseif (index(trim(string),
'FILE:')==1)
then
327 if (string(6:6)==
'.' .or. string(6:6)==
'/')
then
329 filename = trim(
extractword(trim(string(6:80)), 1) )
332 filename = trim(inputdir) // trim(
extractword(trim(string(6:80)), 1) )
334 if (.not.
file_exists(filename))
call mom_error(fatal,trim(mdl)//
", initialize_regridding: "// &
335 "Specified file not found: Looking for '"//trim(filename)//
"' ("//trim(string)//
")")
338 if (len_trim(varname)==0)
then
339 if (field_exists(filename,
'dz')) then; varname =
'dz'
340 elseif (field_exists(filename,
'dsigma')) then; varname =
'dsigma'
341 elseif (field_exists(filename,
'ztest')) then; varname =
'ztest'
342 else ;
call mom_error(fatal,trim(mdl)//
", initialize_regridding: "// &
343 "Coordinate variable not specified and none could be guessed.")
349 if (cs%regridding_scheme == regridding_sigma)
then
350 expected_units =
'nondim'
351 elseif (cs%regridding_scheme == regridding_rho)
then
352 expected_units =
'kg m-3'
354 expected_units =
'meters'
356 if (index(trim(varname),
'interfaces=')==1)
then
357 varname=trim(varname(12:))
358 call check_grid_def(filename, varname, expected_units, message, ierr)
359 if (ierr)
call mom_error(fatal,trim(mdl)//
", initialize_regridding: "//&
360 "Unsupported format in grid definition '"//trim(filename)//
"'. Error message "//trim(message))
361 call field_size(trim(filename), trim(varname), nzf)
363 if (cs%regridding_scheme == regridding_rho)
then
364 allocate(rho_target(ke+1))
365 call mom_read_data(trim(filename), trim(varname), rho_target)
368 allocate(z_max(ke+1))
370 dz(:) = abs(z_max(1:ke) - z_max(2:ke+1))
375 call field_size(trim(filename), trim(varname), nzf)
380 if (main_parameters .and. ke/=gv%ke)
then
381 call mom_error(fatal,trim(mdl)//
', initialize_regridding: '// &
382 'Mismatch in number of model levels and "'//trim(string)//
'".')
384 if (main_parameters)
call log_param(param_file, mdl,
"!"//coord_res_param, dz, &
386 elseif (index(trim(string),
'FNC1:')==1)
then
387 ke = gv%ke;
allocate(dz(ke))
389 if (main_parameters)
call log_param(param_file, mdl,
"!"//coord_res_param, dz, &
391 elseif (index(trim(string),
'RFNC1:')==1)
then
394 elseif (index(trim(string),
'HYBRID:')==1)
then
395 ke = gv%ke;
allocate(dz(ke))
397 allocate(rho_target(ke+1))
398 filename = trim(
extractword(trim(string(8:)), 1) )
399 if (filename(1:1)/=
'.' .and. filename(1:1)/=
'/') filename = trim(inputdir) // trim( filename )
400 if (.not.
file_exists(filename))
call mom_error(fatal,trim(mdl)//
", initialize_regridding: HYBRID "// &
401 "Specified file not found: Looking for '"//trim(filename)//
"' ("//trim(string)//
")")
403 if (.not. field_exists(filename,varname))
call mom_error(fatal,trim(mdl)//
", initialize_regridding: HYBRID "// &
404 "Specified field not found: Looking for '"//trim(varname)//
"' ("//trim(string)//
")")
405 call mom_read_data(trim(filename), trim(varname), rho_target)
407 if (varname(1:5) ==
'FNC1:')
then
408 call dz_function1( trim(string((index(trim(string),
'FNC1:')+5):)), dz )
410 if (.not. field_exists(filename,varname))
call mom_error(fatal,trim(mdl)//
", initialize_regridding: HYBRID "// &
411 "Specified field not found: Looking for '"//trim(varname)//
"' ("//trim(string)//
")")
414 if (main_parameters)
then
415 call log_param(param_file, mdl,
"!"//coord_res_param, dz, &
417 call log_param(param_file, mdl,
"!TARGET_DENSITIES", rho_target, &
418 'HYBRID target densities for interfaces', units=
coordinateunits(coord_mode))
420 elseif (index(trim(string),
'WOA09')==1)
then
421 if (len_trim(string)==5)
then
422 tmpreal = 0. ; ke = 0
423 do while (tmpreal<maximum_depth)
425 tmpreal = tmpreal + woa09_dz(ke)
427 elseif (index(trim(string),
'WOA09:')==1)
then
428 if (len_trim(string)==6)
call mom_error(fatal,trim(mdl)//
', initialize_regridding: '// &
429 'Expected string of form "WOA09:N" but got "'//trim(string)//
'".')
432 if (ke>40 .or. ke<1)
call mom_error(fatal,trim(mdl)//
', initialize_regridding: '// &
433 'For "WOA05:N" N must 0<N<41 but got "'//trim(string)//
'".')
435 dz(1:ke) = woa09_dz(1:ke)
437 call mom_error(fatal,trim(mdl)//
", initialize_regridding: "// &
438 "Unrecognized coordinate configuration"//trim(string))
441 if (main_parameters)
then
443 if (coordinatemode(coord_mode) == regridding_zstar .or. &
448 tmpreal = sum( dz(:) )
449 if (tmpreal < maximum_depth)
then
450 dz(ke) = dz(ke) + ( maximum_depth - tmpreal )
451 elseif (tmpreal > maximum_depth)
then
452 if ( dz(ke) + ( maximum_depth - tmpreal ) > 0. )
then
453 dz(ke) = dz(ke) + ( maximum_depth - tmpreal )
455 call mom_error(fatal,trim(mdl)//
", initialize_regridding: "// &
456 "MAXIMUM_DEPTH was too shallow to adjust bottom layer of DZ!"//trim(string))
465 allocate( cs%coordinateResolution(cs%nk) ); cs%coordinateResolution(:) = -1.e30
468 allocate( cs%target_density(cs%nk+1) ); cs%target_density(:) = -1.e30*us%kg_m3_to_R
471 if (
allocated(dz))
then
472 if (coordinatemode(coord_mode) == regridding_sigma)
then
474 elseif (coordinatemode(coord_mode) == regridding_rho)
then
476 cs%coord_scale = us%R_to_kg_m3
479 cs%coord_scale = gv%H_to_m
482 cs%coord_scale = us%Z_to_m
486 if (
allocated(rho_target))
then
488 deallocate(rho_target)
491 elseif (coordinatemode(coord_mode) == regridding_rho)
then
493 call log_param(param_file, mdl,
"!TARGET_DENSITIES", us%R_to_kg_m3*cs%target_density(:), &
494 'RHO target densities for interfaces', units=
coordinateunits(coord_mode))
500 if (main_parameters .and. coord_is_state_dependent)
then
501 call get_param(param_file, mdl,
"REGRID_COMPRESSIBILITY_FRACTION", tmpreal, &
502 "When interpolating potential density profiles we can add "//&
503 "some artificial compressibility solely to make homogeneous "//&
504 "regions appear stratified.", default=0.)
508 if (main_parameters)
then
509 call get_param(param_file, mdl,
"MIN_THICKNESS", tmpreal, &
510 "When regridding, this is the minimum layer "//&
511 "thickness allowed.", units=
"m", scale=gv%m_to_H, &
520 call get_param(param_file, mdl,
"SLIGHT_DZ_SURFACE", dz_fixed_sfc, &
521 "The nominal thickness of fixed thickness near-surface "//&
522 "layers with the SLight coordinate.", units=
"m", default=1.0, scale=gv%m_to_H)
523 call get_param(param_file, mdl,
"SLIGHT_NZ_SURFACE_FIXED", nz_fixed_sfc, &
524 "The number of fixed-depth surface layers with the SLight "//&
525 "coordinate.", units=
"nondimensional", default=2)
526 call get_param(param_file, mdl,
"SLIGHT_SURFACE_AVG_DEPTH", rho_avg_depth, &
527 "The thickness of the surface region over which to average "//&
528 "when calculating the density to use to define the interior "//&
529 "with the SLight coordinate.", units=
"m", default=1.0, scale=gv%m_to_H)
530 call get_param(param_file, mdl,
"SLIGHT_NLAY_TO_INTERIOR", nlay_sfc_int, &
531 "The number of layers to offset the surface density when "//&
532 "defining where the interior ocean starts with SLight.", &
533 units=
"nondimensional", default=2.0)
534 call get_param(param_file, mdl,
"SLIGHT_FIX_HALOCLINES", fix_haloclines, &
535 "If true, identify regions above the reference pressure "//&
536 "where the reference pressure systematically underestimates "//&
537 "the stratification and use this in the definition of the "//&
538 "interior with the SLight coordinate.", default=.false.)
541 nz_fixed_surface=nz_fixed_sfc, rho_ml_avg_depth=rho_avg_depth, &
542 nlay_ml_to_interior=nlay_sfc_int, fix_haloclines=fix_haloclines)
543 if (fix_haloclines)
then
545 call get_param(param_file, mdl,
"HALOCLINE_FILTER_LENGTH", filt_len, &
546 "A length scale over which to smooth the temperature and "//&
547 "salinity before identifying erroneously unstable haloclines.", &
548 units=
"m", default=2.0)
549 call get_param(param_file, mdl,
"HALOCLINE_STRAT_TOL", strat_tol, &
550 "A tolerance for the ratio of the stratification of the "//&
551 "apparent coordinate stratification to the actual value "//&
552 "that is used to identify erroneously unstable haloclines. "//&
553 "This ratio is 1 when they are equal, and sensible values "//&
554 "are between 0 and 0.5.", units=
"nondimensional", default=0.2)
556 halocline_strat_tol=strat_tol)
562 call get_param(param_file, mdl,
"ADAPT_TIME_RATIO", adapttimeratio, &
563 "Ratio of ALE timestep to grid timescale.", units=
"s", default=1e-1)
564 call get_param(param_file, mdl,
"ADAPT_ZOOM_DEPTH", adaptzoom, &
565 "Depth of near-surface zooming region.", units=
"m", default=200.0, scale=gv%m_to_H)
566 call get_param(param_file, mdl,
"ADAPT_ZOOM_COEFF", adaptzoomcoeff, &
567 "Coefficient of near-surface zooming diffusivity.", &
568 units=
"nondim", default=0.2)
569 call get_param(param_file, mdl,
"ADAPT_BUOY_COEFF", adaptbuoycoeff, &
570 "Coefficient of buoyancy diffusivity.", &
571 units=
"nondim", default=0.8)
572 call get_param(param_file, mdl,
"ADAPT_ALPHA", adaptalpha, &
573 "Scaling on optimization tendency.", &
574 units=
"nondim", default=1.0)
575 call get_param(param_file, mdl,
"ADAPT_DO_MIN_DEPTH", tmplogical, &
576 "If true, make a HyCOM-like mixed layer by preventing interfaces "//&
577 "from being shallower than the depths specified by the regridding coordinate.", &
581 adaptzoomcoeff=adaptzoomcoeff, adaptbuoycoeff=adaptbuoycoeff, adaptalpha=adaptalpha, &
582 adaptdomin=tmplogical)
585 if (main_parameters .and. coord_is_state_dependent)
then
586 call get_param(param_file, mdl,
"MAXIMUM_INT_DEPTH_CONFIG", string, &
587 "Determines how to specify the maximum interface depths.\n"//&
588 "Valid options are:\n"//&
589 " NONE - there are no maximum interface depths\n"//&
590 " PARAM - use the vector-parameter MAXIMUM_INTERFACE_DEPTHS\n"//&
591 " FILE:string - read from a file. The string specifies\n"//&
592 " the filename and variable name, separated\n"//&
593 " by a comma or space, e.g. FILE:lev.nc,Z\n"//&
594 " FNC1:string - FNC1:dz_min,H_total,power,precision",&
596 message =
"The list of maximum depths for each interface."
597 allocate(z_max(ke+1))
599 if ( trim(string) ==
"NONE")
then
601 elseif ( trim(string) ==
"PARAM")
then
602 call get_param(param_file, mdl,
"MAXIMUM_INTERFACE_DEPTHS", z_max, &
603 trim(message), units=
"m", scale=gv%m_to_H, fail_if_missing=.true.)
605 elseif (index(trim(string),
'FILE:')==1)
then
606 if (string(6:6)==
'.' .or. string(6:6)==
'/')
then
608 filename = trim(
extractword(trim(string(6:80)), 1) )
611 filename = trim(inputdir) // trim(
extractword(trim(string(6:80)), 1) )
613 if (.not.
file_exists(filename))
call mom_error(fatal,trim(mdl)//
", initialize_regridding: "// &
614 "Specified file not found: Looking for '"//trim(filename)//
"' ("//trim(string)//
")")
618 if (.not. field_exists(filename,varname))
call mom_error(fatal,trim(mdl)//
", initialize_regridding: "// &
619 "Specified field not found: Looking for '"//trim(varname)//
"' ("//trim(string)//
")")
620 if (len_trim(varname)==0)
then
621 if (field_exists(filename,
'z_max')) then; varname =
'z_max'
622 elseif (field_exists(filename,
'dz')) then; varname =
'dz' ; do_sum = .true.
623 elseif (field_exists(filename,
'dz_max')) then; varname =
'dz_max' ; do_sum = .true.
624 else ;
call mom_error(fatal,trim(mdl)//
", initialize_regridding: "// &
625 "MAXIMUM_INT_DEPTHS variable not specified and none could be guessed.")
630 z_max(1) = 0.0 ;
do k=1,ke ; z_max(k+1) = z_max(k) + dz_max(k) ;
enddo
634 call log_param(param_file, mdl,
"!MAXIMUM_INT_DEPTHS", z_max, &
637 elseif (index(trim(string),
'FNC1:')==1)
then
640 (dz_fixed_sfc > 0.0))
then
641 do k=1,nz_fixed_sfc ; dz_max(k) = dz_fixed_sfc ;
enddo
643 z_max(1) = 0.0 ;
do k=1,ke ; z_max(k+1) = z_max(k) + dz_max(k) ;
enddo
644 call log_param(param_file, mdl,
"!MAXIMUM_INT_DEPTHS", z_max, &
648 call mom_error(fatal,trim(mdl)//
", initialize_regridding: "// &
649 "Unrecognized MAXIMUM_INT_DEPTH_CONFIG "//trim(string))
656 call get_param(param_file, mdl,
"MAX_LAYER_THICKNESS_CONFIG", string, &
657 "Determines how to specify the maximum layer thicknesses.\n"//&
658 "Valid options are:\n"//&
659 " NONE - there are no maximum layer thicknesses\n"//&
660 " PARAM - use the vector-parameter MAX_LAYER_THICKNESS\n"//&
661 " FILE:string - read from a file. The string specifies\n"//&
662 " the filename and variable name, separated\n"//&
663 " by a comma or space, e.g. FILE:lev.nc,Z\n"//&
664 " FNC1:string - FNC1:dz_min,H_total,power,precision",&
666 message =
"The list of maximum thickness for each layer."
668 if ( trim(string) ==
"NONE")
then
670 elseif ( trim(string) ==
"PARAM")
then
671 call get_param(param_file, mdl,
"MAX_LAYER_THICKNESS", h_max, &
672 trim(message), units=
"m", fail_if_missing=.true., scale=gv%m_to_H)
674 elseif (index(trim(string),
'FILE:')==1)
then
675 if (string(6:6)==
'.' .or. string(6:6)==
'/')
then
677 filename = trim(
extractword(trim(string(6:80)), 1) )
680 filename = trim(inputdir) // trim(
extractword(trim(string(6:80)), 1) )
682 if (.not.
file_exists(filename))
call mom_error(fatal,trim(mdl)//
", initialize_regridding: "// &
683 "Specified file not found: Looking for '"//trim(filename)//
"' ("//trim(string)//
")")
686 if (.not. field_exists(filename,varname))
call mom_error(fatal,trim(mdl)//
", initialize_regridding: "// &
687 "Specified field not found: Looking for '"//trim(varname)//
"' ("//trim(string)//
")")
688 if (len_trim(varname)==0)
then
689 if (field_exists(filename,
'h_max')) then; varname =
'h_max'
690 elseif (field_exists(filename,
'dz_max')) then; varname =
'dz_max'
691 else ;
call mom_error(fatal,trim(mdl)//
", initialize_regridding: "// &
692 "MAXIMUM_INT_DEPTHS variable not specified and none could be guessed.")
696 call log_param(param_file, mdl,
"!MAX_LAYER_THICKNESS", h_max, &
699 elseif (index(trim(string),
'FNC1:')==1)
then
701 call log_param(param_file, mdl,
"!MAX_LAYER_THICKNESS", h_max, &
705 call mom_error(fatal,trim(mdl)//
", initialize_regridding: "// &
706 "Unrecognized MAX_LAYER_THICKNESS_CONFIG "//trim(string))
711 if (
allocated(dz))
deallocate(dz)
715 subroutine check_grid_def(filename, varname, expected_units, msg, ierr)
716 character(len=*),
intent(in) :: filename
717 character(len=*),
intent(in) :: varname
718 character(len=*),
intent(in) :: expected_units
719 character(len=*),
intent(inout) :: msg
720 logical,
intent(out) :: ierr
722 character (len=200) :: units, long_name
723 integer :: ncid, status, intid, vid
727 status = nf90_open(trim(filename), nf90_nowrite, ncid)
728 if (status /= nf90_noerr)
then
730 msg =
'File not found: '//trim(filename)
734 status = nf90_inq_varid(ncid, trim(varname), vid)
735 if (status /= nf90_noerr)
then
737 msg =
'Var not found: '//trim(varname)
741 status = nf90_get_att(ncid, vid,
"units", units)
742 if (status /= nf90_noerr)
then
744 msg =
'Attribute not found: units'
750 do i=1,len_trim(units)
751 if (units(i:i) == char(0)) units(i:i) =
" "
754 if (trim(units) /= trim(expected_units))
then
755 if (trim(expected_units) ==
"meters")
then
756 if (trim(units) /=
"m")
then
765 msg =
'Units incorrect: '//trim(units)//
' /= '//trim(expected_units)
776 if (
associated(cs%rho_CS))
call end_coord_rho(cs%rho_CS)
781 deallocate( cs%coordinateResolution )
782 if (
allocated(cs%target_density))
deallocate( cs%target_density )
783 if (
allocated(cs%max_interface_depths) )
deallocate( cs%max_interface_depths )
784 if (
allocated(cs%max_layer_thickness) )
deallocate( cs%max_layer_thickness )
790 subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, frac_shelf_h, conv_adjust)
811 type(ocean_grid_type),
intent(in) :: g
813 real,
dimension(SZI_(G),SZJ_(G), SZK_(GV)),
intent(inout) :: h
816 real,
dimension(SZI_(G),SZJ_(G), CS%nk),
intent(inout) :: h_new
817 real,
dimension(SZI_(G),SZJ_(G), CS%nk+1),
intent(inout) :: dzinterface
818 real,
dimension(:,:),
optional,
pointer :: frac_shelf_h
819 logical,
optional,
intent(in ) :: conv_adjust
821 real :: trickgnucompiler
822 logical :: use_ice_shelf
823 logical :: do_convective_adjustment
825 do_convective_adjustment = .true.
826 if (
present(conv_adjust)) do_convective_adjustment = conv_adjust
828 use_ice_shelf = .false.
829 if (
present(frac_shelf_h))
then
830 if (
associated(frac_shelf_h)) use_ice_shelf = .true.
833 select case ( cs%regridding_scheme )
835 case ( regridding_zstar )
836 if (use_ice_shelf)
then
843 case ( regridding_sigma_shelf_zstar)
847 case ( regridding_sigma )
851 case ( regridding_rho )
856 case ( regridding_arbitrary )
872 call mom_error(fatal,
'MOM_regridding, regridding_main: '//&
873 'Unknown regridding scheme selected!')
877 #ifdef __DO_SAFETY_CHECKS__
886 type(ocean_grid_type),
intent(in) :: G
888 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(in) :: h
889 real,
dimension(SZI_(G),SZJ_(G),CS%nk+1),
intent(in) :: dzInterface
890 real,
dimension(SZI_(G),SZJ_(G),CS%nk),
intent(inout) :: h_new
892 integer :: i, j, k, nki
894 nki = min(cs%nk, gv%ke)
897 do j = g%jsc-1,g%jec+1
898 do i = g%isc-1,g%iec+1
899 if (g%mask2dT(i,j)>0.)
then
901 h_new(i,j,k) = max( 0., h(i,j,k) + ( dzinterface(i,j,k) - dzinterface(i,j,k+1) ) )
903 if (cs%nk > gv%ke)
then
905 h_new(i,j,k) = max( 0., dzinterface(i,j,k) - dzinterface(i,j,k+1) )
909 h_new(i,j,1:nki) = h(i,j,1:nki)
910 if (cs%nk > gv%ke) h_new(i,j,nki+1:cs%nk) = 0.
920 type(ocean_grid_type),
intent(in) :: g
922 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(in) :: h
923 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1),
intent(in) :: dzinterface
925 character(len=*),
intent(in) :: msg
930 do j = g%jsc-1,g%jec+1 ;
do i = g%isc-1,g%iec+1
931 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 )
938 integer,
intent(in) :: nk
939 real,
intent(in) :: depth
940 real,
dimension(nk),
intent(in) :: h
941 real,
dimension(nk+1),
intent(in) :: dzinterface
942 character(len=*),
intent(in) :: msg
945 real :: eps, total_h_old, total_h_new, h_new, z_old, z_new
947 eps =1. ; eps = epsilon(eps)
952 total_h_old = total_h_old + h(k)
957 if (depth == 0.) z_old = - total_h_old
961 z_new = z_old + dzinterface(k)
962 h_new = h(k) + ( dzinterface(k) - dzinterface(k+1) )
964 write(0,*)
'k,h,hnew=',k,h(k),h_new
965 write(0,*)
'dzI(k+1),dzI(k)=',dzinterface(k+1),dzinterface(k)
966 call mom_error( fatal,
'MOM_regridding, check_grid_column: '//&
967 'Negative layer thickness implied by re-gridding, '//trim(msg))
969 total_h_new = total_h_new + h_new
974 if (abs(total_h_new-total_h_old)>real(nk-1)*0.5*(total_h_old+total_h_new)*eps)
then
977 write(0,*)
'k,h,hnew=',k,h(k),h(k)+(dzinterface(k)-dzinterface(k+1))
979 write(0,*)
'Hold,Hnew,Hnew-Hold=',total_h_old,total_h_new,total_h_new-total_h_old
980 write(0,*)
'eps,(n)/2*eps*H=',eps,real(nk-1)*0.5*(total_h_old+total_h_new)*eps
981 call mom_error( fatal,
'MOM_regridding, check_grid_column: '//&
982 'Re-gridding did NOT conserve total thickness to within roundoff '//trim(msg))
986 if (dzinterface(1) /= 0.)
call mom_error( fatal, &
987 'MOM_regridding, check_grid_column: Non-zero dzInterface at surface! '//trim(msg))
988 if (dzinterface(nk+1) /= 0.)
call mom_error( fatal, &
989 'MOM_regridding, check_grid_column: Non-zero dzInterface at bottom! '//trim(msg))
1000 integer,
intent(in) :: nk
1001 real,
dimension(nk+1),
intent(in) :: z_old
1002 real,
dimension(CS%nk+1),
intent(in) :: z_new
1003 real,
dimension(CS%nk+1),
intent(inout) :: dz_g
1006 real :: dz_tgt, zr1, z_old_k
1007 real :: Aq, Bq, dz0, z0, F0
1008 real :: zs, zd, dzwt, Idzwt
1010 real :: Int_zs, Int_zd, dInt_zs_zd
1012 real,
dimension(nk+1) :: z_act
1014 logical :: debug = .false.
1017 if ((z_old(nk+1) - z_old(1)) * (z_new(cs%nk+1) - z_new(1)) < 0.0)
then
1018 call mom_error(fatal,
"filtered_grid_motion: z_old and z_new use different sign conventions.")
1019 elseif ((z_old(nk+1) - z_old(1)) * (z_new(cs%nk+1) - z_new(1)) == 0.0)
then
1021 do k=1,cs%nk+1 ; dz_g(k) = 0.0 ;
enddo ;
return
1022 elseif ((z_old(nk+1) - z_old(1)) + (z_new(cs%nk+1) - z_new(1)) > 0.0)
then
1030 if (sgn*(z_new(k)-z_new(k-1)) < -5e-16*(abs(z_new(k))+abs(z_new(k-1))) ) &
1031 call mom_error(fatal,
"filtered_grid_motion: z_new is tangled.")
1034 if (sgn*(z_old(k)-z_old(k-1)) < -5e-16*(abs(z_old(k))+abs(z_old(k-1))) ) &
1035 call mom_error(fatal,
"filtered_grid_motion: z_old is tangled.")
1040 zs = cs%depth_of_time_filter_shallow
1041 zd = cs%depth_of_time_filter_deep
1042 wtd = 1.0 - cs%old_grid_weight
1046 idzwt = 0.0 ;
if (abs(zd - zs) > 0.0) idzwt = 1.0 / (zd - zs)
1047 dint_zs_zd = 0.5*(1.0 + iwtd) * (zd - zs)
1048 aq = 0.5*(iwtd - 1.0)
1053 if (k<=nk+1) z_old_k = z_old(k)
1055 dz_tgt = sgn*(z_new(k) - z_old_k)
1056 zr1 = sgn*(z_old_k - z_old(1))
1060 if ((zr1 > zd) .and. (zr1 + wtd * dz_tgt > zd))
then
1061 dz_g(k) = sgn * wtd * dz_tgt
1062 elseif ((zr1 < zs) .and. (zr1 + dz_tgt < zs))
then
1063 dz_g(k) = sgn * dz_tgt
1072 int_zd = iwtd*(zd - zr1)
1073 int_zs = int_zd - dint_zs_zd
1074 elseif (zr1 <= zs)
then
1076 int_zd = dint_zs_zd + (zs - zr1)
1079 int_zd = (zd - zr1) * (iwtd*(0.5*(zd+zr1) - zs) + 0.5*(zd - zr1)) * idzwt
1080 int_zs = (zs - zr1) * (0.5*iwtd * ((zr1 - zs)) + (zd - 0.5*(zr1+zs))) * idzwt
1084 if (dz_tgt >= int_zd)
then
1085 dz_g(k) = sgn * ((zd-zr1) + wtd*(dz_tgt - int_zd))
1086 elseif (dz_tgt <= int_zs)
then
1087 dz_g(k) = sgn * ((zs-zr1) + (dz_tgt - int_zs))
1093 if (zr1 <= zs)
then ; dz0 = zs-zr1 ; z0 = zs ; f0 = dz_tgt - int_zs
1094 elseif (zr1 >= zd)
then ; dz0 = zd-zr1 ; z0 = zd ; f0 = dz_tgt - int_zd
1095 else ; dz0 = 0.0 ; z0 = zr1 ; f0 = dz_tgt ;
endif
1097 bq = (dzwt + 2.0*aq*(z0-zs))
1100 dz_g(k) = sgn * (dz0 + 2.0*f0*dzwt / (bq + sqrt(bq**2 + 4.0*aq*f0*dzwt) ))
1122 if (k<=nk+1) z_old_k = z_old(k)
1123 z_act(k) = z_old_k + dz_g(k)
1126 if (sgn*((z_act(k))-z_act(k-1)) < -1e-15*(abs(z_act(k))+abs(z_act(k-1))) ) &
1127 call mom_error(fatal,
"filtered_grid_motion: z_output is tangled.")
1140 type(ocean_grid_type),
intent(in) :: G
1142 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(in) :: h
1143 real,
dimension(SZI_(G),SZJ_(G), CS%nk+1),
intent(inout) :: dzInterface
1145 real,
dimension(:,:),
optional,
pointer :: frac_shelf_h
1149 real :: nominalDepth, totalThickness, dh
1150 real,
dimension(SZK_(GV)+1) :: zOld, zNew
1151 real :: minThickness
1152 logical :: ice_shelf
1155 minthickness = cs%min_thickness
1157 if (
present(frac_shelf_h))
then
1158 if (
associated(frac_shelf_h)) ice_shelf = .true.
1165 do j = g%jsc-1,g%jec+1
1166 do i = g%isc-1,g%iec+1
1168 if (g%mask2dT(i,j)==0.)
then
1169 dzinterface(i,j,:) = 0.
1174 nominaldepth = g%bathyT(i,j)*gv%Z_to_H
1177 totalthickness = 0.0
1179 totalthickness = totalthickness + h(i,j,k)
1182 zold(nz+1) = - nominaldepth
1184 zold(k) = zold(k+1) + h(i,j,k)
1188 if (frac_shelf_h(i,j) > 0.)
then
1190 z_rigid_top = totalthickness-nominaldepth, &
1191 eta_orig=zold(1), zscale=gv%Z_to_H)
1194 znew, zscale=gv%Z_to_H)
1198 znew, zscale=gv%Z_to_H)
1204 #ifdef __DO_SAFETY_CHECKS__
1205 dh=max(nominaldepth,totalthickness)
1206 if (abs(znew(1)-zold(1))>(nz-1)*0.5*epsilon(dh)*dh)
then
1207 write(0,*)
'min_thickness=',minthickness
1208 write(0,*)
'nominalDepth=',nominaldepth,
'totalThickness=',totalthickness
1209 write(0,*)
'dzInterface(1) = ',dzinterface(i,j,1),epsilon(dh),nz
1211 write(0,*) k,zold(k),znew(k)
1214 write(0,*) k,h(i,j,k),znew(k)-znew(k+1),cs%coordinateResolution(k)
1217 'MOM_regridding, build_zstar_grid(): top surface has moved!!!' )
1241 type(ocean_grid_type),
intent(in) :: G
1243 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(in) :: h
1244 real,
dimension(SZI_(G),SZJ_(G), CS%nk+1),
intent(inout) :: dzInterface
1250 real :: nominalDepth, totalThickness, dh
1251 real,
dimension(SZK_(GV)+1) :: zOld, zNew
1255 do i = g%isc-1,g%iec+1
1256 do j = g%jsc-1,g%jec+1
1258 if (g%mask2dT(i,j)==0.)
then
1259 dzinterface(i,j,:) = 0.
1264 nominaldepth = g%bathyT(i,j)*gv%Z_to_H
1267 totalthickness = 0.0
1269 totalthickness = totalthickness + h(i,j,k)
1275 zold(nz+1) = -nominaldepth
1277 zold(k) = zold(k+1) + h(i, j, k)
1282 #ifdef __DO_SAFETY_CHECKS__
1283 dh=max(nominaldepth,totalthickness)
1284 if (abs(znew(1)-zold(1))>(cs%nk-1)*0.5*epsilon(dh)*dh)
then
1285 write(0,*)
'min_thickness=',cs%min_thickness
1286 write(0,*)
'nominalDepth=',nominaldepth,
'totalThickness=',totalthickness
1287 write(0,*)
'dzInterface(1) = ',dzinterface(i,j,1),epsilon(dh),nz,cs%nk
1289 write(0,*) k,zold(k),znew(k)
1292 write(0,*) k,h(i,j,k),znew(k)-znew(k+1),totalthickness*cs%coordinateResolution(k),cs%coordinateResolution(k)
1295 'MOM_regridding, build_sigma_grid: top surface has moved!!!' )
1297 dzinterface(i,j,1) = 0.
1298 dzinterface(i,j,cs%nk+1) = 0.
1310 subroutine build_rho_grid( G, GV, h, tv, dzInterface, remapCS, CS )
1327 type(ocean_grid_type),
intent(in) :: G
1329 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(in) :: h
1331 real,
dimension(SZI_(G),SZJ_(G), SZK_(GV)+1),
intent(inout) :: dzInterface
1339 real :: nominalDepth, totalThickness
1340 real,
dimension(SZK_(GV)+1) :: zOld, zNew
1341 real :: h_neglect, h_neglect_edge
1342 #ifdef __DO_SAFETY_CHECKS__
1347 if (gv%Boussinesq)
then
1348 h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
1350 h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
1355 if (.not.cs%target_density_set)
call mom_error(fatal,
"build_rho_grid: "//&
1356 "Target densities must be set before build_rho_grid is called.")
1359 do j = g%jsc-1,g%jec+1
1360 do i = g%isc-1,g%iec+1
1362 if (g%mask2dT(i,j)==0.)
then
1363 dzinterface(i,j,:) = 0.
1369 nominaldepth = g%bathyT(i,j)*gv%Z_to_H
1371 call build_rho_column(cs%rho_CS, nz, nominaldepth, h(i, j, :), &
1372 tv%T(i, j, :), tv%S(i, j, :), tv%eqn_of_state, znew, &
1373 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
1375 if (cs%integrate_downward_for_e)
then
1378 zold(k+1) = zold(k) - h(i,j,k)
1382 zold(nz+1) = - nominaldepth
1384 zold(k) = zold(k+1) + h(i,j,k)
1391 #ifdef __DO_SAFETY_CHECKS__
1393 if (znew(k) > zold(1))
then
1394 write(0,*)
'zOld=',zold
1395 write(0,*)
'zNew=',znew
1396 call mom_error( fatal,
'MOM_regridding, build_rho_grid: '//&
1397 'interior interface above surface!' )
1399 if (znew(k) > znew(k-1))
then
1400 write(0,*)
'zOld=',zold
1401 write(0,*)
'zNew=',znew
1402 call mom_error( fatal,
'MOM_regridding, build_rho_grid: '//&
1403 'interior interfaces cross!' )
1407 totalthickness = 0.0
1409 totalthickness = totalthickness + h(i,j,k)
1412 dh=max(nominaldepth,totalthickness)
1413 if (abs(znew(1)-zold(1))>(nz-1)*0.5*epsilon(dh)*dh)
then
1414 write(0,*)
'min_thickness=',cs%min_thickness
1415 write(0,*)
'nominalDepth=',nominaldepth,
'totalThickness=',totalthickness
1416 write(0,*)
'zNew(1)-zOld(1) = ',znew(1)-zold(1),epsilon(dh),nz
1418 write(0,*) k,zold(k),znew(k)
1421 write(0,*) k,h(i,j,k),znew(k)-znew(k+1)
1424 'MOM_regridding, build_rho_grid: top surface has moved!!!' )
1441 type(ocean_grid_type),
intent(in) :: G
1443 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(in) :: h
1446 real,
dimension(SZI_(G),SZJ_(G),CS%nk),
intent(inout) :: h_new
1447 real,
dimension(SZI_(G),SZJ_(G),CS%nk+1),
intent(inout) :: dzInterface
1450 real,
dimension(SZK_(GV)+1) :: z_col
1451 real,
dimension(CS%nk+1) :: z_col_new
1452 real,
dimension(SZK_(GV)+1) :: dz_col
1453 real,
dimension(SZK_(GV)) :: p_col
1454 integer :: i, j, k, nki
1456 real :: h_neglect, h_neglect_edge
1459 if (gv%Boussinesq)
then
1460 h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
1462 h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
1465 if (.not.cs%target_density_set)
call mom_error(fatal,
"build_grid_HyCOM1 : "//&
1466 "Target densities must be set before build_grid_HyCOM1 is called.")
1468 nki = min(gv%ke, cs%nk)
1471 do j = g%jsc-1,g%jec+1 ;
do i = g%isc-1,g%iec+1
1472 if (g%mask2dT(i,j)>0.)
then
1474 depth = g%bathyT(i,j) * gv%Z_to_H
1478 z_col(k+1) = z_col(k) + h(i,j,k)
1479 p_col(k) = cs%ref_pressure + cs%compressibility_fraction * &
1480 ( 0.5 * ( z_col(k) + z_col(k+1) ) * gv%H_to_Pa - cs%ref_pressure )
1484 h(i, j, :), tv%T(i, j, :), tv%S(i, j, :), p_col, &
1485 z_col, z_col_new, zscale=gv%Z_to_H, &
1486 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
1492 dz_col(:) = -dz_col(:)
1495 dzinterface(i,j,1:nki+1) = dz_col(1:nki+1)
1496 if (nki<cs%nk) dzinterface(i,j,nki+2:cs%nk+1) = 0.
1499 dzinterface(i,j,:) = 0.
1510 type(ocean_grid_type),
intent(in) :: G
1512 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(in) :: h
1515 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1),
intent(inout) :: dzInterface
1521 integer :: i, j, k, nz
1523 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: tInt, sInt
1526 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: zInt
1527 real,
dimension(SZK_(GV)+1) :: zNext
1535 do k = 2, nz ;
do j = g%jsc-2,g%jec+2 ;
do i = g%isc-2,g%iec+2
1536 tint(i,j,k) = 0.5 * (tv%T(i,j,k-1) + tv%T(i,j,k))
1537 sint(i,j,k) = 0.5 * (tv%S(i,j,k-1) + tv%S(i,j,k))
1538 zint(i,j,k) = zint(i,j,k-1) + h(i,j,k-1)
1539 enddo ;
enddo ;
enddo
1543 tint(:,:,1) = tv%T(:,:,1) ; tint(:,:,nz+1) = tv%T(:,:,nz)
1544 sint(:,:,1) = tv%S(:,:,1) ; sint(:,:,nz+1) = tv%S(:,:,nz)
1547 zint(:,:,nz+1) = zint(:,:,nz) + h(:,:,nz)
1551 do j = g%jsc-1,g%jec+1 ;
do i = g%isc-1,g%iec+1
1552 if (g%mask2dT(i,j) < 0.5)
then
1553 dzinterface(i,j,:) = 0.
1557 call build_adapt_column(cs%adapt_CS, g, gv, tv, i, j, zint, tint, sint, h, znext)
1561 do k = 1, nz+1 ; dzinterface(i,j,k) = -dzinterface(i,j,k) ;
enddo
1576 type(ocean_grid_type),
intent(in) :: G
1578 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(in) :: h
1580 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1),
intent(inout) :: dzInterface
1583 real,
dimension(SZK_(GV)+1) :: z_col
1584 real,
dimension(SZK_(GV)+1) :: z_col_new
1585 real,
dimension(SZK_(GV)+1) :: dz_col
1586 real,
dimension(SZK_(GV)) :: p_col
1588 integer :: i, j, k, nz
1589 real :: h_neglect, h_neglect_edge
1592 if (gv%Boussinesq)
then
1593 h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
1595 h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
1600 if (.not.cs%target_density_set)
call mom_error(fatal,
"build_grid_SLight : "//&
1601 "Target densities must be set before build_grid_SLight is called.")
1604 do j = g%jsc-1,g%jec+1 ;
do i = g%isc-1,g%iec+1
1605 if (g%mask2dT(i,j)>0.)
then
1607 depth = g%bathyT(i,j) * gv%Z_to_H
1610 z_col(k+1) = z_col(k) + h(i,j,k)
1611 p_col(k) = cs%ref_pressure + cs%compressibility_fraction * &
1612 ( 0.5 * ( z_col(k) + z_col(k+1) ) * gv%H_to_Pa - cs%ref_pressure )
1616 gv%H_subroundoff, nz, depth, h(i, j, :), &
1617 tv%T(i, j, :), tv%S(i, j, :), p_col, z_col, z_col_new, &
1618 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
1622 do k=1,nz+1 ; dzinterface(i,j,k) = -dz_col(k) ;
enddo
1623 #ifdef __DO_SAFETY_CHECKS__
1624 if (dzinterface(i,j,1) /= 0.) stop
'build_grid_SLight: Surface moved?!'
1625 if (dzinterface(i,j,nz+1) /= 0.) stop
'build_grid_SLight: Bottom moved?!'
1632 dzinterface(i,j,:) = 0.
1641 integer,
intent(in) :: nk
1642 real,
dimension(nk),
intent(in) :: h_old
1643 real,
dimension(CS%nk+1),
intent(inout) :: dz_int
1646 real :: h_new, eps, h_total, h_err
1648 eps = 1. ; eps = epsilon(eps)
1650 h_total = 0. ; h_err = 0.
1651 do k = 1, min(cs%nk,nk)
1652 h_total = h_total + h_old(k)
1653 h_err = h_err + max( h_old(k), abs(dz_int(k)), abs(dz_int(k+1)) )*eps
1654 h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
1655 if (h_new < -3.0*h_err)
then
1656 write(0,*)
'h<0 at k=',k,
'h_old=',h_old(k), &
1657 'wup=',dz_int(k),
'wdn=',dz_int(k+1),
'dw_dz=',dz_int(k) - dz_int(k+1), &
1658 'h_new=',h_new,
'h_err=',h_err
1659 call mom_error( fatal,
'MOM_regridding: adjust_interface_motion() - '//&
1660 'implied h<0 is larger than roundoff!')
1665 h_err = h_err + max( abs(dz_int(k)), abs(dz_int(k+1)) )*eps
1666 h_new = ( dz_int(k) - dz_int(k+1) )
1667 if (h_new < -3.0*h_err)
then
1668 write(0,*)
'h<0 at k=',k,
'h_old was empty',&
1669 'wup=',dz_int(k),
'wdn=',dz_int(k+1),
'dw_dz=',dz_int(k) - dz_int(k+1), &
1670 'h_new=',h_new,
'h_err=',h_err
1671 call mom_error( fatal,
'MOM_regridding: adjust_interface_motion() - '//&
1672 'implied h<0 is larger than roundoff!')
1676 do k = min(cs%nk,nk),2,-1
1677 h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
1678 if (h_new<cs%min_thickness) &
1679 dz_int(k) = ( dz_int(k+1) - h_old(k) ) + cs%min_thickness
1680 h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
1682 dz_int(k) = ( 1. - eps ) * ( dz_int(k+1) - h_old(k) )
1683 h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
1685 write(0,*)
'h<0 at k=',k,
'h_old=',h_old(k), &
1686 'wup=',dz_int(k),
'wdn=',dz_int(k+1),
'dw_dz=',dz_int(k) - dz_int(k+1), &
1688 stop
'Still did not work!'
1689 call mom_error( fatal,
'MOM_regridding: adjust_interface_motion() - '//&
1690 'Repeated adjustment for roundoff h<0 failed!')
1706 type(ocean_grid_type),
intent(in) :: G
1708 real,
dimension(SZI_(G),SZJ_(G), SZK_(GV)),
intent(in) :: h
1709 real,
dimension(SZI_(G),SZJ_(G), SZK_(GV)+1),
intent(inout) :: dzInterface
1711 real,
intent(inout) :: h_new
1717 real :: z_inter(SZK_(GV)+1)
1718 real :: total_height
1723 real :: x1, y1, x2, y2
1727 max_depth = g%max_depth*gv%Z_to_H
1729 do j = g%jsc-1,g%jec+1
1730 do i = g%isc-1,g%iec+1
1733 local_depth = g%bathyT(i,j)*gv%Z_to_H
1738 total_height = total_height + h(i,j,k)
1741 eta = total_height - local_depth
1744 delta_h = (max_depth + eta) / nz
1749 z_inter(k+1) = z_inter(k) - delta_h
1754 x1 = 0.35; y1 = 0.45; x2 = 0.65; y2 = 0.55
1756 x = - ( z_inter(k) - eta ) / max_depth
1760 elseif ( (x > x1 ) .and. ( x < x2 ))
then
1761 t = y1 + (y2-y1) * (x-x1) / (x2-x1)
1763 t = y2 + (1.0-y2) * (x-x2) / (1.0-x2)
1766 z_inter(k) = -t * max_depth + eta
1771 z_inter(nz+1) = - local_depth
1775 if ( z_inter(k) < (z_inter(k+1) + cs%min_thickness) )
then
1776 z_inter(k) = z_inter(k+1) + cs%min_thickness
1782 dzinterface(i,j,1) = 0.
1785 dzinterface(i,j,k) = z_inter(k) - x
1787 dzinterface(i,j,nz+1) = 0.
1813 type(ocean_grid_type),
intent(in) :: g
1815 real,
dimension(SZI_(G),SZJ_(G), SZK_(GV)),
intent(inout) :: h
1821 do i = g%isc-1,g%iec+1
1822 do j = g%jsc-1,g%jec+1
1844 type(ocean_grid_type),
intent(in) :: G
1846 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1860 logical :: stratified
1861 real,
dimension(GV%ke) :: p_col, densities
1866 do j = g%jsc-1,g%jec+1 ;
do i = g%isc-1,g%iec+1
1870 densities, 1, gv%ke, tv%eqn_of_state )
1877 t0 = tv%T(i,j,k) ; t1 = tv%T(i,j,k+1)
1878 s0 = tv%S(i,j,k) ; s1 = tv%S(i,j,k+1)
1879 r0 = densities(k) ; r1 = densities(k+1)
1880 h0 = h(i,j,k) ; h1 = h(i,j,k+1)
1885 tv%T(i,j,k) = t1 ; tv%T(i,j,k+1) = t0
1886 tv%S(i,j,k) = s1 ; tv%S(i,j,k+1) = s0
1887 h(i,j,k) = h1 ; h(i,j,k+1) = h0
1890 densities(k), tv%eqn_of_state )
1892 densities(k+1), tv%eqn_of_state )
1893 stratified = .false.
1897 if ( stratified )
exit
1912 integer,
intent(in) :: nk
1913 character(len=*),
intent(in) :: coordmode
1916 real,
intent(in) :: maxdepth
1917 real,
intent(in) :: rholight
1918 real,
intent(in) :: rhoheavy
1925 scheme = coordinatemode(coordmode)
1926 select case ( scheme )
1932 case ( regridding_rho )
1935 case ( regridding_sigma )
1939 call mom_error(fatal,
"MOM_regridding, uniformResolution: "//&
1940 "Unrecognized choice for coordinate mode ("//trim(coordmode)//
").")
1948 subroutine initcoord(CS, GV, US, coord_mode)
1950 character(len=*),
intent(in) :: coord_mode
1956 select case (coordinatemode(coord_mode))
1957 case (regridding_zstar)
1959 case (regridding_sigma_shelf_zstar)
1961 case (regridding_sigma)
1963 case (regridding_rho)
1964 call init_coord_rho(cs%rho_CS, cs%nk, cs%ref_pressure, cs%target_density, cs%interp_CS, &
1965 rho_scale=us%kg_m3_to_R)
1967 call init_coord_hycom(cs%hycom_CS, cs%nk, cs%coordinateResolution, cs%target_density, &
1968 cs%interp_CS, rho_scale=us%kg_m3_to_R)
1970 call init_coord_slight(cs%slight_CS, cs%nk, cs%ref_pressure, cs%target_density, &
1971 cs%interp_CS, gv%m_to_H, rho_scale=us%kg_m3_to_R)
1973 call init_coord_adapt(cs%adapt_CS, cs%nk, cs%coordinateResolution, gv%m_to_H)
1980 real,
dimension(:),
intent(in) :: dz
1982 real,
optional,
intent(in) :: scale
1984 if (
size(dz)/=cs%nk)
call mom_error( fatal, &
1985 'setCoordinateResolution: inconsistent number of levels' )
1987 if (
present(scale))
then
1988 cs%coordinateResolution(:) = scale*dz(:)
1990 cs%coordinateResolution(:) = dz(:)
2004 cs%target_density(1) = (gv%Rlay(1) + 0.5*(gv%Rlay(1)-gv%Rlay(2)))
2005 cs%target_density(nz+1) = (gv%Rlay(nz) + 0.5*(gv%Rlay(nz)-gv%Rlay(nz-1)))
2007 cs%target_density(k) = cs%target_density(k-1) + cs%coordinateResolution(k)
2009 cs%target_density_set = .true.
2016 real,
dimension(CS%nk+1),
intent(in) :: rho_int
2018 if (
size(cs%target_density)/=
size(rho_int))
then
2019 call mom_error(fatal,
"set_target_densities inconsistent args!")
2022 cs%target_density(:) = rho_int(:)
2023 cs%target_density_set = .true.
2030 real,
dimension(CS%nk+1),
intent(in) :: max_depths
2031 real,
optional,
intent(in) :: units_to_h
2036 if (.not.
allocated(cs%max_interface_depths))
allocate(cs%max_interface_depths(1:cs%nk+1))
2038 val_to_h = 1.0 ;
if (
present(units_to_h)) val_to_h = units_to_h
2039 if (max_depths(cs%nk+1) < max_depths(1)) val_to_h = -1.0*val_to_h
2042 if (max_depths(cs%nk+1) < max_depths(1))
then
2043 do k=1,cs%nk ;
if (max_depths(k+1) > max_depths(k)) &
2044 call mom_error(fatal,
"Unordered list of maximum depths sent to set_regrid_max_depths!")
2047 do k=1,cs%nk ;
if (max_depths(k+1) < max_depths(k)) &
2048 call mom_error(fatal,
"Unordered list of maximum depths sent to set_regrid_max_depths.")
2053 cs%max_interface_depths(k) = val_to_h * max_depths(k)
2057 select case (cs%regridding_scheme)
2059 call set_hycom_params(cs%hycom_CS, max_interface_depths=cs%max_interface_depths)
2061 call set_slight_params(cs%slight_CS, max_interface_depths=cs%max_interface_depths)
2068 real,
dimension(CS%nk+1),
intent(in) :: max_h
2069 real,
optional,
intent(in) :: units_to_h
2074 if (.not.
allocated(cs%max_layer_thickness))
allocate(cs%max_layer_thickness(1:cs%nk))
2076 val_to_h = 1.0 ;
if (
present( units_to_h)) val_to_h = units_to_h
2079 cs%max_layer_thickness(k) = val_to_h * max_h(k)
2083 select case (cs%regridding_scheme)
2085 call set_hycom_params(cs%hycom_CS, max_layer_thickness=cs%max_layer_thickness)
2087 call set_slight_params(cs%slight_CS, max_layer_thickness=cs%max_layer_thickness)
2096 logical,
optional,
intent(in) :: undo_scaling
2101 unscale = .false. ;
if (
present(undo_scaling)) unscale = undo_scaling
2114 logical,
optional,
intent(in) :: undo_scaling
2120 unscale = .false. ;
if (
present(undo_scaling)) unscale = undo_scaling
2124 if (cs%regridding_scheme == regridding_rho)
then
2125 if (.not. cs%target_density_set) &
2126 call mom_error(fatal,
'MOM_regridding, getCoordinateInterfaces: '//&
2127 'target densities not set!')
2139 cs%coord_scale * cs%coordinateResolution(k)
2145 cs%coordinateResolution(k)
2161 select case ( cs%regridding_scheme )
2164 case ( regridding_sigma_shelf_zstar )
2166 case ( regridding_sigma )
2168 case ( regridding_rho )
2170 case ( regridding_arbitrary )
2173 call mom_error(fatal,
'MOM_regridding, getCoordinateUnits: '//&
2174 'Unknown regridding scheme selected!')
2185 select case ( cs%regridding_scheme )
2186 case ( regridding_zstar )
2190 case ( regridding_sigma_shelf_zstar )
2192 case ( regridding_sigma )
2194 case ( regridding_rho )
2196 case ( regridding_arbitrary )
2205 call mom_error(fatal,
'MOM_regridding, getCoordinateShortName: '//&
2206 'Unknown regridding scheme selected!')
2212 subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_grid_weight, &
2213 interp_scheme, depth_of_time_filter_shallow, depth_of_time_filter_deep, &
2214 compress_fraction, dz_min_surface, nz_fixed_surface, Rho_ML_avg_depth, &
2215 nlay_ML_to_interior, fix_haloclines, halocline_filt_len, &
2216 halocline_strat_tol, integrate_downward_for_e, &
2217 adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, adaptAlpha, adaptDoMin)
2219 logical,
optional,
intent(in) :: boundary_extrapolation
2220 real,
optional,
intent(in) :: min_thickness
2222 real,
optional,
intent(in) :: old_grid_weight
2223 character(len=*),
optional,
intent(in) :: interp_scheme
2224 real,
optional,
intent(in) :: depth_of_time_filter_shallow
2225 real,
optional,
intent(in) :: depth_of_time_filter_deep
2226 real,
optional,
intent(in) :: compress_fraction
2227 real,
optional,
intent(in) :: dz_min_surface
2229 integer,
optional,
intent(in) :: nz_fixed_surface
2230 real,
optional,
intent(in) :: rho_ml_avg_depth
2232 real,
optional,
intent(in) :: nlay_ml_to_interior
2234 logical,
optional,
intent(in) :: fix_haloclines
2235 real,
optional,
intent(in) :: halocline_filt_len
2237 real,
optional,
intent(in) :: halocline_strat_tol
2239 logical,
optional,
intent(in) :: integrate_downward_for_e
2241 real,
optional,
intent(in) :: adapttimeratio
2242 real,
optional,
intent(in) :: adaptzoom
2243 real,
optional,
intent(in) :: adaptzoomcoeff
2244 real,
optional,
intent(in) :: adaptbuoycoeff
2245 real,
optional,
intent(in) :: adaptalpha
2246 logical,
optional,
intent(in) :: adaptdomin
2251 if (
present(boundary_extrapolation))
call set_interp_extrap(cs%interp_CS, boundary_extrapolation)
2253 if (
present(old_grid_weight))
then
2254 if (old_grid_weight<0. .or. old_grid_weight>1.) &
2255 call mom_error(fatal,
'MOM_regridding, set_regrid_params: Weight is out side the range 0..1!')
2256 cs%old_grid_weight = old_grid_weight
2258 if (
present(depth_of_time_filter_shallow)) cs%depth_of_time_filter_shallow = depth_of_time_filter_shallow
2259 if (
present(depth_of_time_filter_deep)) cs%depth_of_time_filter_deep = depth_of_time_filter_deep
2260 if (
present(depth_of_time_filter_shallow) .or.
present(depth_of_time_filter_deep))
then
2261 if (cs%depth_of_time_filter_deep<cs%depth_of_time_filter_shallow)
call mom_error(fatal,
'MOM_regridding, '//&
2262 'set_regrid_params: depth_of_time_filter_deep<depth_of_time_filter_shallow!')
2265 if (
present(min_thickness)) cs%min_thickness = min_thickness
2266 if (
present(compress_fraction)) cs%compressibility_fraction = compress_fraction
2267 if (
present(integrate_downward_for_e)) cs%integrate_downward_for_e = integrate_downward_for_e
2269 select case (cs%regridding_scheme)
2270 case (regridding_zstar)
2271 if (
present(min_thickness))
call set_zlike_params(cs%zlike_CS, min_thickness=min_thickness)
2272 case (regridding_sigma_shelf_zstar)
2273 if (
present(min_thickness))
call set_zlike_params(cs%zlike_CS, min_thickness=min_thickness)
2274 case (regridding_sigma)
2275 if (
present(min_thickness))
call set_sigma_params(cs%sigma_CS, min_thickness=min_thickness)
2276 case (regridding_rho)
2277 if (
present(min_thickness))
call set_rho_params(cs%rho_CS, min_thickness=min_thickness)
2278 if (
present(integrate_downward_for_e)) &
2279 call set_rho_params(cs%rho_CS, integrate_downward_for_e=integrate_downward_for_e)
2280 if (
associated(cs%rho_CS) .and. (
present(interp_scheme) .or.
present(boundary_extrapolation))) &
2281 call set_rho_params(cs%rho_CS, interp_cs=cs%interp_CS)
2283 if (
associated(cs%hycom_CS) .and. (
present(interp_scheme) .or.
present(boundary_extrapolation))) &
2286 if (
present(min_thickness))
call set_slight_params(cs%slight_CS, min_thickness=min_thickness)
2287 if (
present(dz_min_surface))
call set_slight_params(cs%slight_CS, dz_ml_min=dz_min_surface)
2288 if (
present(nz_fixed_surface))
call set_slight_params(cs%slight_CS, nz_fixed_surface=nz_fixed_surface)
2289 if (
present(rho_ml_avg_depth))
call set_slight_params(cs%slight_CS, rho_ml_avg_depth=rho_ml_avg_depth)
2290 if (
present(nlay_ml_to_interior))
call set_slight_params(cs%slight_CS, nlay_ml_offset=nlay_ml_to_interior)
2291 if (
present(fix_haloclines))
call set_slight_params(cs%slight_CS, fix_haloclines=fix_haloclines)
2292 if (
present(halocline_filt_len))
call set_slight_params(cs%slight_CS, halocline_filter_length=halocline_filt_len)
2293 if (
present(halocline_strat_tol))
call set_slight_params(cs%slight_CS, halocline_strat_tol=halocline_strat_tol)
2294 if (
present(compress_fraction))
call set_slight_params(cs%slight_CS, compressibility_fraction=compress_fraction)
2295 if (
associated(cs%slight_CS) .and. (
present(interp_scheme) .or.
present(boundary_extrapolation))) &
2298 if (
present(adapttimeratio))
call set_adapt_params(cs%adapt_CS, adapttimeratio=adapttimeratio)
2299 if (
present(adaptzoom))
call set_adapt_params(cs%adapt_CS, adaptzoom=adaptzoom)
2300 if (
present(adaptzoomcoeff))
call set_adapt_params(cs%adapt_CS, adaptzoomcoeff=adaptzoomcoeff)
2301 if (
present(adaptbuoycoeff))
call set_adapt_params(cs%adapt_CS, adaptbuoycoeff=adaptbuoycoeff)
2302 if (
present(adaptalpha))
call set_adapt_params(cs%adapt_CS, adaptalpha=adaptalpha)
2303 if (
present(adaptdomin))
call set_adapt_params(cs%adapt_CS, adaptdomin=adaptdomin)
2344 real,
intent(in) :: ssh
2345 real,
intent(in) :: depth
2351 select case ( cs%regridding_scheme )
2356 dz = cs%coordinateResolution(k) * ( 1. + ssh/depth )
2358 dz = min(dz, depth - z)
2365 case ( regridding_sigma )
2367 case ( regridding_rho )
2369 case ( regridding_arbitrary )
2372 call mom_error(fatal,
'MOM_regridding, getStaticThickness: '//&
2373 'Unknown regridding scheme selected!')
2380 character(len=*),
intent(in) :: string
2382 real,
dimension(:),
intent(inout) :: dz
2385 real :: dz_min, power, prec, H_total
2389 read( string, *) dz_min, h_total, power, prec
2390 if (prec == -1024.)
call mom_error(fatal,
"dz_function1: "// &
2391 "Problem reading FNC1: string ="//trim(string))
2394 dz(k) = (real(k-1)/real(nk-1))**power
2396 dz(:) = ( h_total - real(nk) * dz_min ) * ( dz(:) / sum(dz) )
2397 dz(:) = anint( dz(:) / prec ) * prec
2398 dz(:) = ( h_total - real(nk) * dz_min ) * ( dz(:) / sum(dz) )
2399 dz(:) = anint( dz(:) / prec ) * prec
2400 dz(nk) = dz(nk) + ( h_total - sum( dz(:) + dz_min ) )
2401 dz(:) = anint( dz(:) / prec ) * prec
2402 dz(:) = dz(:) + dz_min
2409 character(len=*),
intent(in) :: string
2411 real,
dimension(:),
allocatable,
intent(inout) :: rho_target
2413 integer :: nki, k, nk
2414 real :: ddx, dx, rho_1, rho_2, rho_3, drho, rho_4, drho_min
2416 read( string, *) nk, rho_1, rho_2, rho_3, drho, rho_4, drho_min
2417 allocate(rho_target(nk+1))
2419 rho_target(1) = rho_1
2420 rho_target(2) = rho_2
2423 ddx = max( drho_min, real(nki-k)/real(nki*nki) )
2425 rho_target(3+k) = rho_3 + (2. * drho) * dx
2427 rho_target(nki+4) = rho_4