25 use mom_io,
only : create_file, write_field, close_file
27 use mom_regridding,
only : initialize_regridding, regridding_main, end_regridding
30 use mom_regridding,
only : set_target_densities_from_gv, set_target_densities
31 use mom_regridding,
only : regriddingcoordinatemodedoc, default_coordinate_mode
32 use mom_regridding,
only : regriddinginterpschemedoc, regriddingdefaultinterpscheme
36 use mom_regridding,
only : getcoordinateinterfaces, getcoordinateresolution
37 use mom_regridding,
only : getcoordinateunits, getcoordinateshortname
41 use mom_remapping,
only : remappingschemesdoc, remappingdefaultscheme
52 use plm_functions,
only : plm_reconstruction, plm_boundary_extrapolation
53 use ppm_functions,
only : ppm_reconstruction, ppm_boundary_extrapolation
54 use p1m_functions,
only : p1m_interpolation, p1m_boundary_extrapolation
55 use p3m_functions,
only : p3m_interpolation, p3m_boundary_extrapolation
58 implicit none ;
private
59 #include <MOM_memory.h>
64 logical :: remap_uv_using_old_alg
68 real :: regrid_time_scale
76 logical :: remap_after_initialization
78 logical :: show_call_tree
82 integer,
dimension(:),
allocatable :: id_tracer_remap_tendency
83 integer,
dimension(:),
allocatable :: id_htracer_remap_tendency
84 integer,
dimension(:),
allocatable :: id_htracer_remap_tendency_2d
85 logical,
dimension(:),
allocatable :: do_tendency_diag
86 integer :: id_dzregrid = -1
89 integer :: id_u_preale = -1
90 integer :: id_v_preale = -1
91 integer :: id_h_preale = -1
92 integer :: id_t_preale = -1
93 integer :: id_s_preale = -1
94 integer :: id_e_preale = -1
95 integer :: id_vert_remap_h = -1
96 integer :: id_vert_remap_h_tendency = -1
104 public ale_main_offline
105 public ale_offline_inputs
106 public ale_offline_tracer_final
107 public ale_build_grid
108 public ale_regrid_accelerated
109 public ale_remap_scalar
110 public pressure_gradient_plm
111 public pressure_gradient_ppm
112 public adjustgridforintegrity
113 public ale_initregridding
114 public ale_getcoordinate
115 public ale_getcoordinateunits
116 public ale_writecoordinatefile
117 public ale_updateverticalgridtype
118 public ale_initthicknesstocoord
119 public ale_update_regrid_weights
120 public ale_remap_init_conds
121 public ale_register_diags
134 subroutine ale_init( param_file, GV, US, max_depth, CS)
138 real,
intent(in) :: max_depth
139 type(
ale_cs),
pointer :: cs
142 real,
dimension(:),
allocatable :: dz
143 character(len=40) :: mdl =
"MOM_ALE"
144 character(len=80) :: string
145 real :: filter_shallow_depth, filter_deep_depth
146 logical :: check_reconstruction
147 logical :: check_remapping
148 logical :: force_bounds_in_subcell
149 logical :: local_logical
150 logical :: remap_boundary_extrap
152 if (
associated(cs))
then
153 call mom_error(warning,
"ALE_init called with an associated "// &
154 "control structure.")
159 cs%show_call_tree = calltree_showquery()
160 if (cs%show_call_tree)
call calltree_enter(
"ALE_init(), MOM_ALE.F90")
162 call get_param(param_file, mdl,
"REMAP_UV_USING_OLD_ALG", &
163 cs%remap_uv_using_old_alg, &
164 "If true, uses the old remapping-via-a-delta-z method for "//&
165 "remapping u and v. If false, uses the new method that remaps "//&
166 "between grids described by an old and new thickness.", &
170 call ale_initregridding(gv, us, max_depth, param_file, mdl, cs%regridCS)
173 call get_param(param_file, mdl,
"REMAPPING_SCHEME", string, &
174 "This sets the reconstruction scheme used "//&
175 "for vertical remapping for all variables. "//&
176 "It can be one of the following schemes: "//&
177 trim(remappingschemesdoc), default=remappingdefaultscheme)
178 call get_param(param_file, mdl,
"FATAL_CHECK_RECONSTRUCTIONS", check_reconstruction, &
179 "If true, cell-by-cell reconstructions are checked for "//&
180 "consistency and if non-monotonicity or an inconsistency is "//&
181 "detected then a FATAL error is issued.", default=.false.)
182 call get_param(param_file, mdl,
"FATAL_CHECK_REMAPPING", check_remapping, &
183 "If true, the results of remapping are checked for "//&
184 "conservation and new extrema and if an inconsistency is "//&
185 "detected then a FATAL error is issued.", default=.false.)
186 call get_param(param_file, mdl,
"REMAP_BOUND_INTERMEDIATE_VALUES", force_bounds_in_subcell, &
187 "If true, the values on the intermediate grid used for remapping "//&
188 "are forced to be bounded, which might not be the case due to "//&
189 "round off.", default=.false.)
190 call get_param(param_file, mdl,
"REMAP_BOUNDARY_EXTRAP", remap_boundary_extrap, &
191 "If true, values at the interfaces of boundary cells are "//&
192 "extrapolated instead of piecewise constant", default=.false.)
193 call initialize_remapping( cs%remapCS, string, &
194 boundary_extrapolation=remap_boundary_extrap, &
195 check_reconstruction=check_reconstruction, &
196 check_remapping=check_remapping, &
197 force_bounds_in_subcell=force_bounds_in_subcell)
199 call get_param(param_file, mdl,
"REMAP_AFTER_INITIALIZATION", cs%remap_after_initialization, &
200 "If true, applies regridding and remapping immediately after "//&
201 "initialization so that the state is ALE consistent. This is a "//&
202 "legacy step and should not be needed if the initialization is "//&
203 "consistent with the coordinate mode.", default=.true.)
205 call get_param(param_file, mdl,
"REGRID_TIME_SCALE", cs%regrid_time_scale, &
206 "The time-scale used in blending between the current (old) grid "//&
207 "and the target (new) grid. A short time-scale favors the target "//&
208 "grid (0. or anything less than DT_THERM) has no memory of the old "//&
209 "grid. A very long time-scale makes the model more Lagrangian.", &
210 units=
"s", default=0.)
211 call get_param(param_file, mdl,
"REGRID_FILTER_SHALLOW_DEPTH", filter_shallow_depth, &
212 "The depth above which no time-filtering is applied. Above this depth "//&
213 "final grid exactly matches the target (new) grid.", &
214 units=
"m", default=0., scale=gv%m_to_H)
215 call get_param(param_file, mdl,
"REGRID_FILTER_DEEP_DEPTH", filter_deep_depth, &
216 "The depth below which full time-filtering is applied with time-scale "//&
217 "REGRID_TIME_SCALE. Between depths REGRID_FILTER_SHALLOW_DEPTH and "//&
218 "REGRID_FILTER_SHALLOW_DEPTH the filter weights adopt a cubic profile.", &
219 units=
"m", default=0., scale=gv%m_to_H)
220 call set_regrid_params(cs%regridCS, depth_of_time_filter_shallow=filter_shallow_depth, &
221 depth_of_time_filter_deep=filter_deep_depth)
222 call get_param(param_file, mdl,
"REGRID_USE_OLD_DIRECTION", local_logical, &
223 "If true, the regridding ntegrates upwards from the bottom for "//&
224 "interface positions, much as the main model does. If false "//&
225 "regridding integrates downward, consistant with the remapping "//&
226 "code.", default=.true., do_not_log=.true.)
227 call set_regrid_params(cs%regridCS, integrate_downward_for_e=.not.local_logical)
232 if (cs%show_call_tree)
call calltree_leave(
"ALE_init()")
233 end subroutine ale_init
236 subroutine ale_register_diags(Time, G, GV, US, diag, CS)
237 type(time_type),
target,
intent(in) :: time
238 type(ocean_grid_type),
intent(in) :: g
241 type(
diag_ctrl),
target,
intent(in) :: diag
242 type(
ale_cs),
pointer :: cs
248 cs%id_u_preale = register_diag_field(
'ocean_model',
'u_preale', diag%axesCuL, time, &
249 'Zonal velocity before remapping',
'm s-1')
250 cs%id_v_preale = register_diag_field(
'ocean_model',
'v_preale', diag%axesCvL, time, &
251 'Meridional velocity before remapping',
'm s-1')
252 cs%id_h_preale = register_diag_field(
'ocean_model',
'h_preale', diag%axesTL, time, &
253 'Layer Thickness before remapping', get_thickness_units(gv), v_extensive=.true.)
254 cs%id_T_preale = register_diag_field(
'ocean_model',
'T_preale', diag%axesTL, time, &
255 'Temperature before remapping',
'degC')
256 cs%id_S_preale = register_diag_field(
'ocean_model',
'S_preale', diag%axesTL, time, &
257 'Salinity before remapping',
'PSU')
258 cs%id_e_preale = register_diag_field(
'ocean_model',
'e_preale', diag%axesTi, time, &
259 'Interface Heights before remapping',
'm', conversion=us%Z_to_m)
261 cs%id_dzRegrid = register_diag_field(
'ocean_model',
'dzRegrid',diag%axesTi,time, &
262 'Change in interface height due to ALE regridding',
'm')
263 cs%id_vert_remap_h = register_diag_field(
'ocean_model',
'vert_remap_h',diag%axestl,time, &
264 'layer thicknesses after ALE regridding and remapping',
'm', v_extensive = .true.)
265 cs%id_vert_remap_h_tendency = register_diag_field(
'ocean_model',
'vert_remap_h_tendency',diag%axestl,time, &
266 'Layer thicknesses tendency due to ALE regridding and remapping',
'm', v_extensive = .true.)
268 end subroutine ale_register_diags
275 subroutine adjustgridforintegrity( CS, G, GV, h )
277 type(ocean_grid_type),
intent(in) :: g
279 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(inout) :: h
281 call inflate_vanished_layers_old( cs%regridCS, g, gv, h(:,:,:) )
283 end subroutine adjustgridforintegrity
289 subroutine ale_end(CS)
293 call end_remapping( cs%remapCS )
294 call end_regridding( cs%regridCS )
298 end subroutine ale_end
304 subroutine ale_main( G, GV, US, h, u, v, tv, Reg, CS, dt, frac_shelf_h)
305 type(ocean_grid_type),
intent(in) :: g
308 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(inout) :: h
310 real,
dimension(SZIB_(G),SZJ_(G),SZK_(GV)),
intent(inout) :: u
311 real,
dimension(SZI_(G),SZJB_(G),SZK_(GV)),
intent(inout) :: v
314 type(
ale_cs),
pointer :: cs
315 real,
optional,
intent(in) :: dt
316 real,
dimension(:,:),
optional,
pointer :: frac_shelf_h
318 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: dzregrid
319 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: eta_preale
320 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_new
321 integer :: nk, i, j, k, isc, iec, jsc, jec
324 nk = gv%ke; isc = g%isc; iec = g%iec; jsc = g%jsc; jec = g%jec
327 if (
present(frac_shelf_h))
then
328 if (
associated(frac_shelf_h)) ice_shelf = .true.
331 if (cs%show_call_tree)
call calltree_enter(
"ALE_main(), MOM_ALE.F90")
334 if (cs%id_u_preale > 0)
call post_data(cs%id_u_preale, u, cs%diag)
335 if (cs%id_v_preale > 0)
call post_data(cs%id_v_preale, v, cs%diag)
336 if (cs%id_h_preale > 0)
call post_data(cs%id_h_preale, h, cs%diag)
337 if (cs%id_T_preale > 0)
call post_data(cs%id_T_preale, tv%T, cs%diag)
338 if (cs%id_S_preale > 0)
call post_data(cs%id_S_preale, tv%S, cs%diag)
339 if (cs%id_e_preale > 0)
then
340 call find_eta(h, tv, g, gv, us, eta_preale)
341 call post_data(cs%id_e_preale, eta_preale, cs%diag)
344 if (
present(dt))
then
345 call ale_update_regrid_weights( dt, cs )
347 dzregrid(:,:,:) = 0.0
352 call regridding_main( cs%remapCS, cs%regridCS, g, gv, h, tv, h_new, dzregrid, frac_shelf_h)
354 call regridding_main( cs%remapCS, cs%regridCS, g, gv, h, tv, h_new, dzregrid)
357 call check_grid( g, gv, h, 0. )
359 if (cs%show_call_tree)
call calltree_waypoint(
"new grid generated (ALE_main)")
363 if (
present(dt))
then
364 call diag_update_remap_grids(cs%diag)
367 call remap_all_state_vars( cs%remapCS, cs, g, gv, h, h_new, reg, -dzregrid, &
368 u, v, cs%show_call_tree, dt )
370 if (cs%show_call_tree)
call calltree_waypoint(
"state remapped (ALE_main)")
375 do k = 1,nk ;
do j = jsc-1,jec+1 ;
do i = isc-1,iec+1
376 h(i,j,k) = h_new(i,j,k)
377 enddo ;
enddo ;
enddo
379 if (cs%show_call_tree)
call calltree_leave(
"ALE_main()")
381 if (cs%id_dzRegrid>0 .and.
present(dt))
call post_data(cs%id_dzRegrid, dzregrid, cs%diag)
384 end subroutine ale_main
390 subroutine ale_main_offline( G, GV, h, tv, Reg, CS, dt)
391 type(ocean_grid_type),
intent(in) :: g
393 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(inout) :: h
397 type(
ale_cs),
pointer :: cs
398 real,
optional,
intent(in) :: dt
400 real,
dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzregrid
401 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_new
402 integer :: nk, i, j, k, isc, iec, jsc, jec
404 nk = gv%ke; isc = g%isc; iec = g%iec; jsc = g%jsc; jec = g%jec
406 if (cs%show_call_tree)
call calltree_enter(
"ALE_main_offline(), MOM_ALE.F90")
408 if (
present(dt))
then
409 call ale_update_regrid_weights( dt, cs )
411 dzregrid(:,:,:) = 0.0
415 call regridding_main( cs%remapCS, cs%regridCS, g, gv, h, tv, h_new, dzregrid )
417 call check_grid( g, gv, h, 0. )
419 if (cs%show_call_tree)
call calltree_waypoint(
"new grid generated (ALE_main)")
423 call remap_all_state_vars( cs%remapCS, cs, g, gv, h, h_new, reg, debug=cs%show_call_tree, dt=dt )
425 if (cs%show_call_tree)
call calltree_waypoint(
"state remapped (ALE_main)")
430 do k = 1,nk ;
do j = jsc-1,jec+1 ;
do i = isc-1,iec+1
431 h(i,j,k) = h_new(i,j,k)
432 enddo ;
enddo ;
enddo
434 if (cs%show_call_tree)
call calltree_leave(
"ALE_main()")
435 if (cs%id_dzRegrid>0 .and.
present(dt))
call post_data(cs%id_dzRegrid, dzregrid, cs%diag)
437 end subroutine ale_main_offline
442 subroutine ale_offline_inputs(CS, G, GV, h, tv, Reg, uhtr, vhtr, Kd, debug)
444 type(ocean_grid_type),
intent(in ) :: g
446 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(inout) :: h
449 real,
dimension(SZIB_(G),SZJ_(G),SZK_(GV)),
intent(inout) :: uhtr
450 real,
dimension(SZI_(G),SZJB_(G),SZK_(GV)),
intent(inout) :: vhtr
451 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1),
intent(inout) :: kd
452 logical,
intent(in ) :: debug
454 integer :: nk, i, j, k, isc, iec, jsc, jec
455 real,
dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: h_new
456 real,
dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzregrid
457 real,
dimension(SZK_(GV)) :: h_src
458 real,
dimension(SZK_(GV)) :: h_dest, uh_dest
459 real,
dimension(SZK_(GV)) :: temp_vec
461 nk = gv%ke; isc = g%isc; iec = g%iec; jsc = g%jsc; jec = g%jec
462 dzregrid(:,:,:) = 0.0
465 if (debug)
call mom_tracer_chkinv(
"Before ALE_offline_inputs", g, h, reg%Tr, reg%ntr)
470 call regridding_main( cs%remapCS, cs%regridCS, g, gv, h, tv, h_new, dzregrid, conv_adjust = .false. )
471 call check_grid( g, gv, h_new, 0. )
472 if (cs%show_call_tree)
call calltree_waypoint(
"new grid generated (ALE_offline_inputs)")
475 call remap_all_state_vars( cs%remapCS, cs, g, gv, h, h_new, reg, debug=cs%show_call_tree )
476 if (cs%show_call_tree)
call calltree_waypoint(
"state remapped (ALE_inputs)")
479 do j=jsc,jec ;
do i=g%iscB,g%iecB
480 if (g%mask2dCu(i,j)>0.)
then
481 h_src(:) = 0.5 * (h(i,j,:) + h(i+1,j,:))
482 h_dest(:) = 0.5 * (h_new(i,j,:) + h_new(i+1,j,:))
483 call reintegrate_column(nk, h_src, uhtr(i,j,:), nk, h_dest, 0., temp_vec)
484 uhtr(i,j,:) = temp_vec
487 do j=g%jscB,g%jecB ;
do i=isc,iec
488 if (g%mask2dCv(i,j)>0.)
then
489 h_src(:) = 0.5 * (h(i,j,:) + h(i,j+1,:))
490 h_dest(:) = 0.5 * (h_new(i,j,:) + h_new(i,j+1,:))
491 call reintegrate_column(nk, h_src, vhtr(i,j,:), nk, h_dest, 0., temp_vec)
492 vhtr(i,j,:) = temp_vec
496 do j = jsc,jec ;
do i=isc,iec
497 if (g%mask2dT(i,j)>0.)
then
498 if (check_column_integrals(nk, h_src, nk, h_dest))
then
499 call mom_error(fatal,
"ALE_offline_inputs: Kd interpolation columns do not match")
501 call interpolate_column(nk, h(i,j,:), kd(i,j,:), nk, h_new(i,j,:), 0., kd(i,j,:))
505 call ale_remap_scalar(cs%remapCS, g, gv, nk, h, tv%T, h_new, tv%T)
506 call ale_remap_scalar(cs%remapCS, g, gv, nk, h, tv%S, h_new, tv%S)
508 if (debug)
call mom_tracer_chkinv(
"After ALE_offline_inputs", g, h_new, reg%Tr, reg%ntr)
511 do k = 1,nk ;
do j = jsc-1,jec+1 ;
do i = isc-1,iec+1
512 h(i,j,k) = h_new(i,j,k)
513 enddo ;
enddo ;
enddo
515 if (cs%show_call_tree)
call calltree_leave(
"ALE_offline_inputs()")
516 end subroutine ale_offline_inputs
522 subroutine ale_offline_tracer_final( G, GV, h, tv, h_target, Reg, CS)
523 type(ocean_grid_type),
intent(in) :: g
525 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(inout) :: h
528 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(inout) :: h_target
531 type(
ale_cs),
pointer :: cs
534 real,
dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzregrid
535 real,
dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: h_new
536 integer :: nk, i, j, k, isc, iec, jsc, jec
538 nk = gv%ke; isc = g%isc; iec = g%iec; jsc = g%jsc; jec = g%jec
540 if (cs%show_call_tree)
call calltree_enter(
"ALE_offline_tracer_final(), MOM_ALE.F90")
542 call regridding_main( cs%remapCS, cs%regridCS, g, gv, h_target, tv, h_new, dzregrid )
543 call check_grid( g, gv, h_target, 0. )
546 if (cs%show_call_tree)
call calltree_waypoint(
"Source and target grids checked (ALE_offline_tracer_final)")
550 call remap_all_state_vars( cs%remapCS, cs, g, gv, h, h_new, reg, debug=cs%show_call_tree )
552 if (cs%show_call_tree)
call calltree_waypoint(
"state remapped (ALE_offline_tracer_final)")
558 do j = jsc-1,jec+1 ;
do i = isc-1,iec+1
559 h(i,j,k) = h_new(i,j,k)
562 if (cs%show_call_tree)
call calltree_leave(
"ALE_offline_tracer_final()")
563 end subroutine ale_offline_tracer_final
566 subroutine check_grid( G, GV, h, threshold )
567 type(ocean_grid_type),
intent(in) :: G
569 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(in) :: h
571 real,
intent(in) :: threshold
576 do j = g%jsc,g%jec ;
do i = g%isc,g%iec
577 if (g%mask2dT(i,j)>0.)
then
578 if (minval(h(i,j,:)) < threshold)
then
579 write(0,*)
'check_grid: i,j=',i,j,
'h(i,j,:)=',h(i,j,:)
580 if (threshold <= 0.)
then
581 call mom_error(fatal,
"MOM_ALE, check_grid: negative thickness encountered.")
583 call mom_error(fatal,
"MOM_ALE, check_grid: too tiny thickness encountered.")
590 end subroutine check_grid
593 subroutine ale_build_grid( G, GV, regridCS, remapCS, h, tv, debug, frac_shelf_h )
594 type(ocean_grid_type),
intent(in) :: g
599 real,
dimension(SZI_(G),SZJ_(G), SZK_(GV)),
intent(inout) :: h
601 logical,
optional,
intent(in) :: debug
602 real,
dimension(:,:),
optional,
pointer :: frac_shelf_h
604 integer :: nk, i, j, k
605 real,
dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzregrid
606 real,
dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: h_new
607 logical :: show_call_tree, use_ice_shelf
609 show_call_tree = .false.
610 if (
present(debug)) show_call_tree = debug
611 if (show_call_tree)
call calltree_enter(
"ALE_build_grid(), MOM_ALE.F90")
612 use_ice_shelf = .false.
613 if (
present(frac_shelf_h))
then
614 if (
associated(frac_shelf_h)) use_ice_shelf = .true.
619 if (use_ice_shelf)
then
620 call regridding_main( remapcs, regridcs, g, gv, h, tv, h_new, dzregrid, frac_shelf_h )
622 call regridding_main( remapcs, regridcs, g, gv, h, tv, h_new, dzregrid )
628 do j = g%jsc,g%jec ;
do i = g%isc,g%iec
629 if (g%mask2dT(i,j)>0.) h(i,j,:) = h_new(i,j,:)
632 if (show_call_tree)
call calltree_leave(
"ALE_build_grid()")
633 end subroutine ale_build_grid
637 subroutine ale_regrid_accelerated(CS, G, GV, h, tv, n, u, v, Reg, dt, dzRegrid, initial)
639 type(ocean_grid_type),
intent(inout) :: g
641 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
644 integer,
intent(in) :: n
645 real,
dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
647 real,
dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
650 optional,
pointer :: reg
651 real,
optional,
intent(in) :: dt
652 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
653 optional,
intent(inout) :: dzregrid
654 logical,
optional,
intent(in) :: initial
658 integer :: i, j, k, nz
660 type(group_pass_type) :: pass_t_s_h
661 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_loc, h_orig
662 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
target :: t, s
665 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: dzinterface, dzinttotal
670 dzinttotal(:,:,:) = 0.
678 t(:,:,:) = tv%T(:,:,:)
679 s(:,:,:) = tv%S(:,:,:)
684 h_loc(:,:,:) = h(:,:,:)
685 h_orig(:,:,:) = h(:,:,:)
689 call ale_update_regrid_weights(dt, cs)
692 call do_group_pass(pass_t_s_h, g%domain)
695 call regridding_main(cs%remapCS, cs%regridCS, g, gv, h_loc, tv_local, h, dzinterface)
696 dzinttotal(:,:,:) = dzinttotal(:,:,:) + dzinterface(:,:,:)
699 do j = g%jsc-1,g%jec+1 ;
do i = g%isc-1,g%iec+1
700 call remapping_core_h(cs%remapCS, nz, h_orig(i,j,:), tv%S(i,j,:), nz, h(i,j,:), tv_local%S(i,j,:))
701 call remapping_core_h(cs%remapCS, nz, h_orig(i,j,:), tv%T(i,j,:), nz, h(i,j,:), tv_local%T(i,j,:))
705 h_loc(:,:,:) = h(:,:,:)
709 call remap_all_state_vars(cs%remapCS, cs, g, gv, h_orig, h, reg, dzinttotal, u, v, dt=dt)
712 if (
present(dzregrid)) dzregrid(:,:,:) = dzinttotal(:,:,:)
713 end subroutine ale_regrid_accelerated
721 subroutine remap_all_state_vars(CS_remapping, CS_ALE, G, GV, h_old, h_new, Reg, dxInterface, u, v, debug, dt)
723 type(
ale_cs),
intent(in) :: CS_ALE
724 type(ocean_grid_type),
intent(in) :: G
726 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(in) :: h_old
728 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(in) :: h_new
731 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
732 optional,
intent(in) :: dxInterface
734 real,
dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
735 optional,
intent(inout) :: u
736 real,
dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
737 optional,
intent(inout) :: v
738 logical,
optional,
intent(in) :: debug
739 real,
optional,
intent(in) :: dt
741 integer :: i, j, k, m
743 real,
dimension(GV%ke+1) :: dx
744 real,
dimension(GV%ke) :: h1, u_column
745 real,
dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: work_conc
746 real,
dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: work_cont
747 real,
dimension(SZI_(G), SZJ_(G)) :: work_2d
749 real,
dimension(GV%ke) :: h2
750 real :: h_neglect, h_neglect_edge
751 logical :: show_call_tree
754 show_call_tree = .false.
755 if (
present(debug)) show_call_tree = debug
756 if (show_call_tree)
call calltree_enter(
"remap_all_state_vars(), MOM_ALE.F90")
760 if ( .not.
present(dxinterface) .and. (cs_ale%remap_uv_using_old_alg .and. (
present(u) .or.
present(v))) )
then
761 call mom_error(fatal,
"remap_all_state_vars: dxInterface must be present if using old algorithm "// &
762 "and u/v are to be remapped")
766 if (gv%Boussinesq)
then
767 h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
769 h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
775 ntr = 0 ;
if (
associated(reg)) ntr = reg%ntr
777 if (
present(dt))
then
779 work_conc(:,:,:) = 0.0
780 work_cont(:,:,:) = 0.0
785 if (show_call_tree)
call calltree_waypoint(
"remapping tracers (remap_all_state_vars)")
789 do j = g%jsc,g%jec ;
do i = g%isc,g%iec ;
if (g%mask2dT(i,j)>0.)
then
793 call remapping_core_h(cs_remapping, nz, h1, tr%t(i,j,:), nz, h2, &
794 u_column, h_neglect, h_neglect_edge)
797 if (
present(dt))
then
798 if (tr%id_remap_conc>0)
then
800 work_conc(i,j,k) = (u_column(k) - tr%t(i,j,k) ) * idt
803 if (tr%id_remap_cont>0. .or. tr%id_remap_cont_2d>0)
then
805 work_cont(i,j,k) = (u_column(k)*h2(k) - tr%t(i,j,k)*h1(k)) * idt
810 tr%t(i,j,:) = u_column(:)
811 endif ;
enddo ;
enddo
814 if (
present(dt))
then
815 if (tr%id_remap_conc > 0)
then
816 call post_data(tr%id_remap_conc, work_conc, cs_ale%diag)
818 if (tr%id_remap_cont > 0)
then
819 call post_data(tr%id_remap_cont, work_cont, cs_ale%diag)
821 if (tr%id_remap_cont_2d > 0)
then
822 do j = g%jsc,g%jec ;
do i = g%isc,g%iec
825 work_2d(i,j) = work_2d(i,j) + work_cont(i,j,k)
828 call post_data(tr%id_remap_cont_2d, work_2d, cs_ale%diag)
835 if (show_call_tree)
call calltree_waypoint(
"tracers remapped (remap_all_state_vars)")
838 if (
present(u) )
then
840 do j = g%jsc,g%jec ;
do i = g%iscB,g%iecB ;
if (g%mask2dCu(i,j)>0.)
then
842 h1(:) = 0.5 * ( h_old(i,j,:) + h_old(i+1,j,:) )
843 if (cs_ale%remap_uv_using_old_alg)
then
844 dx(:) = 0.5 * ( dxinterface(i,j,:) + dxinterface(i+1,j,:) )
846 h2(k) = max( 0., h1(k) + ( dx(k+1) - dx(k) ) )
849 h2(:) = 0.5 * ( h_new(i,j,:) + h_new(i+1,j,:) )
851 call remapping_core_h(cs_remapping, nz, h1, u(i,j,:), nz, h2, &
852 u_column, h_neglect, h_neglect_edge)
853 u(i,j,:) = u_column(:)
854 endif ;
enddo ;
enddo
857 if (show_call_tree)
call calltree_waypoint(
"u remapped (remap_all_state_vars)")
860 if (
present(v) )
then
862 do j = g%jscB,g%jecB ;
do i = g%isc,g%iec ;
if (g%mask2dCv(i,j)>0.)
then
864 h1(:) = 0.5 * ( h_old(i,j,:) + h_old(i,j+1,:) )
865 if (cs_ale%remap_uv_using_old_alg)
then
866 dx(:) = 0.5 * ( dxinterface(i,j,:) + dxinterface(i,j+1,:) )
868 h2(k) = max( 0., h1(k) + ( dx(k+1) - dx(k) ) )
871 h2(:) = 0.5 * ( h_new(i,j,:) + h_new(i,j+1,:) )
873 call remapping_core_h(cs_remapping, nz, h1, v(i,j,:), nz, h2, &
874 u_column, h_neglect, h_neglect_edge)
875 v(i,j,:) = u_column(:)
876 endif ;
enddo ;
enddo
879 if (cs_ale%id_vert_remap_h > 0)
call post_data(cs_ale%id_vert_remap_h, h_old, cs_ale%diag)
880 if ((cs_ale%id_vert_remap_h_tendency > 0) .and.
present(dt))
then
881 do k = 1, nz ;
do j = g%jsc,g%jec ;
do i = g%isc,g%iec
882 work_cont(i,j,k) = (h_new(i,j,k) - h_old(i,j,k))*idt
883 enddo ;
enddo ;
enddo
884 call post_data(cs_ale%id_vert_remap_h_tendency, work_cont, cs_ale%diag)
886 if (show_call_tree)
call calltree_waypoint(
"v remapped (remap_all_state_vars)")
887 if (show_call_tree)
call calltree_leave(
"remap_all_state_vars()")
889 end subroutine remap_all_state_vars
895 subroutine ale_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_cells, old_remap )
897 type(ocean_grid_type),
intent(in) :: g
899 integer,
intent(in) :: nk_src
900 real,
dimension(SZI_(G),SZJ_(G),nk_src),
intent(in) :: h_src
902 real,
dimension(SZI_(G),SZJ_(G),nk_src),
intent(in) :: s_src
903 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(in) :: h_dst
905 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(inout) :: s_dst
906 logical,
optional,
intent(in) :: all_cells
909 logical,
optional,
intent(in) :: old_remap
912 integer :: i, j, k, n_points
914 real :: h_neglect, h_neglect_edge
915 logical :: ignore_vanished_layers, use_remapping_core_w
917 ignore_vanished_layers = .false.
918 if (
present(all_cells)) ignore_vanished_layers = .not. all_cells
919 use_remapping_core_w = .false.
920 if (
present(old_remap)) use_remapping_core_w = old_remap
924 if (gv%Boussinesq)
then
925 h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
927 h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
931 do j = g%jsc,g%jec ;
do i = g%isc,g%iec
932 if (g%mask2dT(i,j) > 0.)
then
933 if (ignore_vanished_layers)
then
936 if (h_src(i,j,k)>0.) n_points = n_points + 1
940 if (use_remapping_core_w)
then
941 call dzfromh1h2( n_points, h_src(i,j,1:n_points), gv%ke, h_dst(i,j,:), dx )
942 call remapping_core_w(cs, n_points, h_src(i,j,1:n_points), s_src(i,j,1:n_points), &
943 gv%ke, dx, s_dst(i,j,:), h_neglect, h_neglect_edge)
945 call remapping_core_h(cs, n_points, h_src(i,j,1:n_points), s_src(i,j,1:n_points), &
946 gv%ke, h_dst(i,j,:), s_dst(i,j,:), h_neglect, h_neglect_edge)
953 end subroutine ale_remap_scalar
961 subroutine pressure_gradient_plm( CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_extrap )
962 type(ocean_grid_type),
intent(in) :: g
964 type(
ale_cs),
intent(inout) :: cs
965 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
967 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
969 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
971 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
974 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
976 logical,
intent(in) :: bdry_extrap
983 real,
dimension(CS%nk,2) :: ppol_e
984 real,
dimension(CS%nk,2) :: ppol_coefs
988 if (gv%Boussinesq)
then
989 h_neglect = gv%m_to_H*1.0e-30
991 h_neglect = gv%kg_m2_to_H*1.0e-30
996 do j = g%jsc-1,g%jec+1 ;
do i = g%isc-1,g%iec+1
1002 ppol_coefs(:,:) = 0.0
1003 call plm_reconstruction( gv%ke, htmp, tmp, ppol_e, ppol_coefs, h_neglect )
1005 call plm_boundary_extrapolation( gv%ke, htmp, tmp, ppol_e, ppol_coefs, h_neglect )
1008 s_t(i,j,k) = ppol_e(k,1)
1009 s_b(i,j,k) = ppol_e(k,2)
1014 ppol_coefs(:,:) = 0.0
1015 tmp(:) = tv%T(i,j,:)
1016 call plm_reconstruction( gv%ke, htmp, tmp, ppol_e, ppol_coefs, h_neglect )
1018 call plm_boundary_extrapolation( gv%ke, htmp, tmp, ppol_e, ppol_coefs, h_neglect )
1021 t_t(i,j,k) = ppol_e(k,1)
1022 t_b(i,j,k) = ppol_e(k,2)
1027 end subroutine pressure_gradient_plm
1035 subroutine pressure_gradient_ppm( CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_extrap )
1036 type(ocean_grid_type),
intent(in) :: g
1038 type(
ale_cs),
intent(inout) :: cs
1039 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1040 intent(inout) :: s_t
1041 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1042 intent(inout) :: s_b
1043 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1044 intent(inout) :: t_t
1045 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1046 intent(inout) :: t_b
1048 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1050 logical,
intent(in) :: bdry_extrap
1057 real,
dimension(CS%nk,2) :: &
1059 real,
dimension(CS%nk,3) :: &
1061 real :: h_neglect, h_neglect_edge
1064 if (gv%Boussinesq)
then
1065 h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
1067 h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
1072 do j = g%jsc-1,g%jec+1 ;
do i = g%isc-1,g%iec+1
1076 tmp(:) = tv%S(i,j,:)
1080 ppol_coefs(:,:) = 0.0
1082 call edge_values_implicit_h4( gv%ke, htmp, tmp, ppol_e, h_neglect=h_neglect_edge )
1083 call ppm_reconstruction( gv%ke, htmp, tmp, ppol_e, ppol_coefs, h_neglect )
1085 call ppm_boundary_extrapolation( gv%ke, htmp, tmp, ppol_e, ppol_coefs, h_neglect )
1088 s_t(i,j,k) = ppol_e(k,1)
1089 s_b(i,j,k) = ppol_e(k,2)
1094 ppol_coefs(:,:) = 0.0
1095 tmp(:) = tv%T(i,j,:)
1097 call edge_values_implicit_h4( gv%ke, htmp, tmp, ppol_e, h_neglect=1.0e-10*gv%m_to_H )
1098 call ppm_reconstruction( gv%ke, htmp, tmp, ppol_e, ppol_coefs, h_neglect )
1100 call ppm_boundary_extrapolation(gv%ke, htmp, tmp, ppol_e, ppol_coefs, h_neglect )
1103 t_t(i,j,k) = ppol_e(k,1)
1104 t_b(i,j,k) = ppol_e(k,2)
1109 end subroutine pressure_gradient_ppm
1113 subroutine ale_initregridding(GV, US, max_depth, param_file, mdl, regridCS)
1116 real,
intent(in) :: max_depth
1118 character(len=*),
intent(in) :: mdl
1121 character(len=30) :: coord_mode
1123 call get_param(param_file, mdl,
"REGRIDDING_COORDINATE_MODE", coord_mode, &
1124 "Coordinate mode for vertical regridding. "//&
1125 "Choose among the following possibilities: "//&
1126 trim(regriddingcoordinatemodedoc), &
1127 default=default_coordinate_mode, fail_if_missing=.true.)
1129 call initialize_regridding(regridcs, gv, us, max_depth, param_file, mdl, coord_mode,
'',
'')
1131 end subroutine ale_initregridding
1134 function ale_getcoordinate( CS )
1137 real,
dimension(CS%nk+1) :: ale_getcoordinate
1138 ale_getcoordinate(:) = getcoordinateinterfaces( cs%regridCS, undo_scaling=.true. )
1140 end function ale_getcoordinate
1144 function ale_getcoordinateunits( CS )
1147 character(len=20) :: ale_getcoordinateunits
1149 ale_getcoordinateunits = getcoordinateunits( cs%regridCS )
1151 end function ale_getcoordinateunits
1155 logical function ale_remap_init_conds( CS )
1158 ale_remap_init_conds = .false.
1159 if (
associated(cs)) ale_remap_init_conds = cs%remap_after_initialization
1160 end function ale_remap_init_conds
1163 subroutine ale_update_regrid_weights( dt, CS )
1164 real,
intent(in) :: dt
1165 type(
ale_cs),
pointer :: cs
1169 if (
associated(cs))
then
1171 if (cs%regrid_time_scale > 0.0)
then
1172 w = cs%regrid_time_scale / (cs%regrid_time_scale + dt)
1174 call set_regrid_params(cs%regridCS, old_grid_weight=w)
1177 end subroutine ale_update_regrid_weights
1182 subroutine ale_updateverticalgridtype(CS, GV)
1189 gv%sInterface(1:nk+1) = getcoordinateinterfaces( cs%regridCS, undo_scaling=.true. )
1190 gv%sLayer(1:nk) = 0.5*( gv%sInterface(1:nk) + gv%sInterface(2:nk+1) )
1191 gv%zAxisUnits = getcoordinateunits( cs%regridCS )
1192 gv%zAxisLongName = getcoordinateshortname( cs%regridCS )
1196 end subroutine ale_updateverticalgridtype
1202 subroutine ale_writecoordinatefile( CS, GV, directory )
1205 character(len=*),
intent(in) :: directory
1207 character(len=240) :: filepath
1209 type(fieldtype) :: fields(2)
1211 real :: ds(gv%ke), dsi(gv%ke+1)
1213 filepath = trim(directory) // trim(
"Vertical_coordinate")
1214 ds(:) = getcoordinateresolution( cs%regridCS, undo_scaling=.true. )
1216 dsi(2:gv%ke) = 0.5*( ds(1:gv%ke-1) + ds(2:gv%ke) )
1217 dsi(gv%ke+1) = 0.5*ds(gv%ke)
1219 vars(1) = var_desc(
'ds', getcoordinateunits( cs%regridCS ), &
1220 'Layer Coordinate Thickness',
'1',
'L',
'1')
1221 vars(2) = var_desc(
'ds_interface', getcoordinateunits( cs%regridCS ), &
1222 'Layer Center Coordinate Separation',
'1',
'i',
'1')
1224 call create_file(unit, trim(filepath), vars, 2, fields, single_file, gv=gv)
1225 call write_field(unit, fields(1), ds)
1226 call write_field(unit, fields(2), dsi)
1227 call close_file(unit)
1229 end subroutine ale_writecoordinatefile
1233 subroutine ale_initthicknesstocoord( CS, G, GV, h )
1235 type(ocean_grid_type),
intent(in) :: g
1237 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(out) :: h
1242 do j = g%jsd,g%jed ;
do i = g%isd,g%ied
1243 h(i,j,:) = gv%Z_to_H * getstaticthickness( cs%regridCS, 0., g%bathyT(i,j) )
1246 end subroutine ale_initthicknesstocoord