60 implicit none ;
private
61 #include <MOM_memory.h>
66 logical :: remap_uv_using_old_alg
70 real :: regrid_time_scale
78 logical :: remap_after_initialization
80 logical :: answers_2018
84 logical :: show_call_tree
88 integer,
dimension(:),
allocatable :: id_tracer_remap_tendency
89 integer,
dimension(:),
allocatable :: id_htracer_remap_tendency
90 integer,
dimension(:),
allocatable :: id_htracer_remap_tendency_2d
91 logical,
dimension(:),
allocatable :: do_tendency_diag
92 integer :: id_dzregrid = -1
95 integer :: id_u_preale = -1
96 integer :: id_v_preale = -1
97 integer :: id_h_preale = -1
98 integer :: id_t_preale = -1
99 integer :: id_s_preale = -1
100 integer :: id_e_preale = -1
101 integer :: id_vert_remap_h = -1
102 integer :: id_vert_remap_h_tendency = -1
140 subroutine ale_init( param_file, GV, US, max_depth, CS)
144 real,
intent(in) :: max_depth
145 type(
ale_cs),
pointer :: cs
148 real,
dimension(:),
allocatable :: dz
149 character(len=40) :: mdl =
"MOM_ALE"
150 character(len=80) :: string
151 real :: filter_shallow_depth, filter_deep_depth
152 logical :: default_2018_answers
153 logical :: check_reconstruction
154 logical :: check_remapping
155 logical :: force_bounds_in_subcell
156 logical :: local_logical
157 logical :: remap_boundary_extrap
159 if (
associated(cs))
then
160 call mom_error(warning,
"ALE_init called with an associated "// &
161 "control structure.")
166 cs%show_call_tree = calltree_showquery()
167 if (cs%show_call_tree)
call calltree_enter(
"ALE_init(), MOM_ALE.F90")
169 call get_param(param_file, mdl,
"REMAP_UV_USING_OLD_ALG", &
170 cs%remap_uv_using_old_alg, &
171 "If true, uses the old remapping-via-a-delta-z method for "//&
172 "remapping u and v. If false, uses the new method that remaps "//&
173 "between grids described by an old and new thickness.", &
180 call get_param(param_file, mdl,
"REMAPPING_SCHEME", string, &
181 "This sets the reconstruction scheme used "//&
182 "for vertical remapping for all variables. "//&
183 "It can be one of the following schemes: "//&
184 trim(remappingschemesdoc), default=remappingdefaultscheme)
185 call get_param(param_file, mdl,
"FATAL_CHECK_RECONSTRUCTIONS", check_reconstruction, &
186 "If true, cell-by-cell reconstructions are checked for "//&
187 "consistency and if non-monotonicity or an inconsistency is "//&
188 "detected then a FATAL error is issued.", default=.false.)
189 call get_param(param_file, mdl,
"FATAL_CHECK_REMAPPING", check_remapping, &
190 "If true, the results of remapping are checked for "//&
191 "conservation and new extrema and if an inconsistency is "//&
192 "detected then a FATAL error is issued.", default=.false.)
193 call get_param(param_file, mdl,
"REMAP_BOUND_INTERMEDIATE_VALUES", force_bounds_in_subcell, &
194 "If true, the values on the intermediate grid used for remapping "//&
195 "are forced to be bounded, which might not be the case due to "//&
196 "round off.", default=.false.)
197 call get_param(param_file, mdl,
"REMAP_BOUNDARY_EXTRAP", remap_boundary_extrap, &
198 "If true, values at the interfaces of boundary cells are "//&
199 "extrapolated instead of piecewise constant", default=.false.)
200 call get_param(param_file, mdl,
"DEFAULT_2018_ANSWERS", default_2018_answers, &
201 "This sets the default value for the various _2018_ANSWERS parameters.", &
203 call get_param(param_file, mdl,
"REMAPPING_2018_ANSWERS", cs%answers_2018, &
204 "If true, use the order of arithmetic and expressions that recover the "//&
205 "answers from the end of 2018. Otherwise, use updated and more robust "//&
206 "forms of the same expressions.", default=default_2018_answers)
207 call initialize_remapping( cs%remapCS, string, &
208 boundary_extrapolation=remap_boundary_extrap, &
209 check_reconstruction=check_reconstruction, &
210 check_remapping=check_remapping, &
211 force_bounds_in_subcell=force_bounds_in_subcell, &
212 answers_2018=cs%answers_2018)
214 call get_param(param_file, mdl,
"REMAP_AFTER_INITIALIZATION", cs%remap_after_initialization, &
215 "If true, applies regridding and remapping immediately after "//&
216 "initialization so that the state is ALE consistent. This is a "//&
217 "legacy step and should not be needed if the initialization is "//&
218 "consistent with the coordinate mode.", default=.true.)
220 call get_param(param_file, mdl,
"REGRID_TIME_SCALE", cs%regrid_time_scale, &
221 "The time-scale used in blending between the current (old) grid "//&
222 "and the target (new) grid. A short time-scale favors the target "//&
223 "grid (0. or anything less than DT_THERM) has no memory of the old "//&
224 "grid. A very long time-scale makes the model more Lagrangian.", &
225 units=
"s", default=0., scale=us%s_to_T)
226 call get_param(param_file, mdl,
"REGRID_FILTER_SHALLOW_DEPTH", filter_shallow_depth, &
227 "The depth above which no time-filtering is applied. Above this depth "//&
228 "final grid exactly matches the target (new) grid.", &
229 units=
"m", default=0., scale=gv%m_to_H)
230 call get_param(param_file, mdl,
"REGRID_FILTER_DEEP_DEPTH", filter_deep_depth, &
231 "The depth below which full time-filtering is applied with time-scale "//&
232 "REGRID_TIME_SCALE. Between depths REGRID_FILTER_SHALLOW_DEPTH and "//&
233 "REGRID_FILTER_SHALLOW_DEPTH the filter weights adopt a cubic profile.", &
234 units=
"m", default=0., scale=gv%m_to_H)
235 call set_regrid_params(cs%regridCS, depth_of_time_filter_shallow=filter_shallow_depth, &
236 depth_of_time_filter_deep=filter_deep_depth)
237 call get_param(param_file, mdl,
"REGRID_USE_OLD_DIRECTION", local_logical, &
238 "If true, the regridding ntegrates upwards from the bottom for "//&
239 "interface positions, much as the main model does. If false "//&
240 "regridding integrates downward, consistant with the remapping "//&
241 "code.", default=.true., do_not_log=.true.)
242 call set_regrid_params(cs%regridCS, integrate_downward_for_e=.not.local_logical)
252 type(time_type),
target,
intent(in) :: time
253 type(ocean_grid_type),
intent(in) :: g
256 type(
diag_ctrl),
target,
intent(in) :: diag
257 type(
ale_cs),
pointer :: cs
263 cs%id_u_preale = register_diag_field(
'ocean_model',
'u_preale', diag%axesCuL, time, &
264 'Zonal velocity before remapping',
'm s-1', conversion=us%L_T_to_m_s)
265 cs%id_v_preale = register_diag_field(
'ocean_model',
'v_preale', diag%axesCvL, time, &
266 'Meridional velocity before remapping',
'm s-1', conversion=us%L_T_to_m_s)
267 cs%id_h_preale = register_diag_field(
'ocean_model',
'h_preale', diag%axesTL, time, &
269 conversion=gv%H_to_MKS, v_extensive=.true.)
270 cs%id_T_preale = register_diag_field(
'ocean_model',
'T_preale', diag%axesTL, time, &
271 'Temperature before remapping',
'degC')
272 cs%id_S_preale = register_diag_field(
'ocean_model',
'S_preale', diag%axesTL, time, &
273 'Salinity before remapping',
'PSU')
274 cs%id_e_preale = register_diag_field(
'ocean_model',
'e_preale', diag%axesTi, time, &
275 'Interface Heights before remapping',
'm', conversion=us%Z_to_m)
277 cs%id_dzRegrid = register_diag_field(
'ocean_model',
'dzRegrid',diag%axesTi,time, &
278 'Change in interface height due to ALE regridding',
'm', &
279 conversion=gv%H_to_m)
280 cs%id_vert_remap_h = register_diag_field(
'ocean_model',
'vert_remap_h', &
281 diag%axestl, time,
'layer thicknesses after ALE regridding and remapping',
'm', &
282 conversion=gv%H_to_m, v_extensive=.true.)
283 cs%id_vert_remap_h_tendency = register_diag_field(
'ocean_model',
'vert_remap_h_tendency',diag%axestl,time, &
284 'Layer thicknesses tendency due to ALE regridding and remapping',
'm', &
285 conversion=gv%H_to_m*us%s_to_T, v_extensive = .true.)
296 type(ocean_grid_type),
intent(in) :: g
298 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(inout) :: h
300 call inflate_vanished_layers_old( cs%regridCS, g, gv, h(:,:,:) )
312 call end_remapping( cs%remapCS )
313 call end_regridding( cs%regridCS )
323 subroutine ale_main( G, GV, US, h, u, v, tv, Reg, CS, OBC, dt, frac_shelf_h)
324 type(ocean_grid_type),
intent(in) :: g
327 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(inout) :: h
329 real,
dimension(SZIB_(G),SZJ_(G),SZK_(GV)),
intent(inout) :: u
330 real,
dimension(SZI_(G),SZJB_(G),SZK_(GV)),
intent(inout) :: v
333 type(
ale_cs),
pointer :: cs
335 real,
optional,
intent(in) :: dt
336 real,
dimension(:,:),
optional,
pointer :: frac_shelf_h
338 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: dzregrid
339 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: eta_preale
340 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_new
341 integer :: nk, i, j, k, isc, iec, jsc, jec
344 nk = gv%ke; isc = g%isc; iec = g%iec; jsc = g%jsc; jec = g%jec
347 if (
present(frac_shelf_h))
then
348 if (
associated(frac_shelf_h)) ice_shelf = .true.
351 if (cs%show_call_tree)
call calltree_enter(
"ALE_main(), MOM_ALE.F90")
354 if (cs%id_u_preale > 0)
call post_data(cs%id_u_preale, u, cs%diag)
355 if (cs%id_v_preale > 0)
call post_data(cs%id_v_preale, v, cs%diag)
356 if (cs%id_h_preale > 0)
call post_data(cs%id_h_preale, h, cs%diag)
357 if (cs%id_T_preale > 0)
call post_data(cs%id_T_preale, tv%T, cs%diag)
358 if (cs%id_S_preale > 0)
call post_data(cs%id_S_preale, tv%S, cs%diag)
359 if (cs%id_e_preale > 0)
then
360 call find_eta(h, tv, g, gv, us, eta_preale)
361 call post_data(cs%id_e_preale, eta_preale, cs%diag)
364 if (
present(dt))
then
367 dzregrid(:,:,:) = 0.0
372 call regridding_main( cs%remapCS, cs%regridCS, g, gv, h, tv, h_new, dzregrid, frac_shelf_h)
374 call regridding_main( cs%remapCS, cs%regridCS, g, gv, h, tv, h_new, dzregrid)
383 if (
present(dt))
then
388 u, v, cs%show_call_tree, dt )
395 do k = 1,nk ;
do j = jsc-1,jec+1 ;
do i = isc-1,iec+1
396 h(i,j,k) = h_new(i,j,k)
397 enddo ;
enddo ;
enddo
401 if (cs%id_dzRegrid>0 .and.
present(dt))
call post_data(cs%id_dzRegrid, dzregrid, cs%diag)
411 type(ocean_grid_type),
intent(in) :: g
413 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(inout) :: h
417 type(
ale_cs),
pointer :: cs
419 real,
optional,
intent(in) :: dt
421 real,
dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzregrid
422 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_new
423 integer :: nk, i, j, k, isc, iec, jsc, jec
425 nk = gv%ke; isc = g%isc; iec = g%iec; jsc = g%jsc; jec = g%jec
427 if (cs%show_call_tree)
call calltree_enter(
"ALE_main_offline(), MOM_ALE.F90")
429 if (
present(dt))
then
432 dzregrid(:,:,:) = 0.0
436 call regridding_main( cs%remapCS, cs%regridCS, g, gv, h, tv, h_new, dzregrid )
445 debug=cs%show_call_tree, dt=dt )
452 do k = 1,nk ;
do j = jsc-1,jec+1 ;
do i = isc-1,iec+1
453 h(i,j,k) = h_new(i,j,k)
454 enddo ;
enddo ;
enddo
457 if (cs%id_dzRegrid>0 .and.
present(dt))
call post_data(cs%id_dzRegrid, dzregrid, cs%diag)
464 subroutine ale_offline_inputs(CS, G, GV, h, tv, Reg, uhtr, vhtr, Kd, debug, OBC)
466 type(ocean_grid_type),
intent(in ) :: g
468 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(inout) :: h
471 real,
dimension(SZIB_(G),SZJ_(G),SZK_(GV)),
intent(inout) :: uhtr
472 real,
dimension(SZI_(G),SZJB_(G),SZK_(GV)),
intent(inout) :: vhtr
473 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1),
intent(inout) :: kd
474 logical,
intent(in ) :: debug
477 integer :: nk, i, j, k, isc, iec, jsc, jec
478 real,
dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: h_new
479 real,
dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzregrid
480 real,
dimension(SZK_(GV)) :: h_src
481 real,
dimension(SZK_(GV)) :: h_dest, uh_dest
482 real,
dimension(SZK_(GV)) :: temp_vec
484 nk = gv%ke; isc = g%isc; iec = g%iec; jsc = g%jsc; jec = g%jec
485 dzregrid(:,:,:) = 0.0
488 if (debug)
call mom_tracer_chkinv(
"Before ALE_offline_inputs", g, h, reg%Tr, reg%ntr)
493 call regridding_main( cs%remapCS, cs%regridCS, g, gv, h, tv, h_new, dzregrid, conv_adjust = .false. )
495 if (cs%show_call_tree)
call calltree_waypoint(
"new grid generated (ALE_offline_inputs)")
502 do j=jsc,jec ;
do i=g%iscB,g%iecB
503 if (g%mask2dCu(i,j)>0.)
then
504 h_src(:) = 0.5 * (h(i,j,:) + h(i+1,j,:))
505 h_dest(:) = 0.5 * (h_new(i,j,:) + h_new(i+1,j,:))
507 uhtr(i,j,:) = temp_vec
510 do j=g%jscB,g%jecB ;
do i=isc,iec
511 if (g%mask2dCv(i,j)>0.)
then
512 h_src(:) = 0.5 * (h(i,j,:) + h(i,j+1,:))
513 h_dest(:) = 0.5 * (h_new(i,j,:) + h_new(i,j+1,:))
515 vhtr(i,j,:) = temp_vec
519 do j = jsc,jec ;
do i=isc,iec
520 if (g%mask2dT(i,j)>0.)
then
522 call mom_error(fatal,
"ALE_offline_inputs: Kd interpolation columns do not match")
524 call interpolate_column(nk, h(i,j,:), kd(i,j,:), nk, h_new(i,j,:), 0., kd(i,j,:))
531 if (debug)
call mom_tracer_chkinv(
"After ALE_offline_inputs", g, h_new, reg%Tr, reg%ntr)
534 do k = 1,nk ;
do j = jsc-1,jec+1 ;
do i = isc-1,iec+1
535 h(i,j,k) = h_new(i,j,k)
536 enddo ;
enddo ;
enddo
538 if (cs%show_call_tree)
call calltree_leave(
"ALE_offline_inputs()")
546 type(ocean_grid_type),
intent(in) :: g
548 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(inout) :: h
551 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(inout) :: h_target
554 type(
ale_cs),
pointer :: cs
558 real,
dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzregrid
559 real,
dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: h_new
560 integer :: nk, i, j, k, isc, iec, jsc, jec
562 nk = gv%ke; isc = g%isc; iec = g%iec; jsc = g%jsc; jec = g%jec
564 if (cs%show_call_tree)
call calltree_enter(
"ALE_offline_tracer_final(), MOM_ALE.F90")
566 call regridding_main( cs%remapCS, cs%regridCS, g, gv, h_target, tv, h_new, dzregrid )
570 if (cs%show_call_tree)
call calltree_waypoint(
"Source and target grids checked (ALE_offline_tracer_final)")
576 if (cs%show_call_tree)
call calltree_waypoint(
"state remapped (ALE_offline_tracer_final)")
582 do j = jsc-1,jec+1 ;
do i = isc-1,iec+1
583 h(i,j,k) = h_new(i,j,k)
586 if (cs%show_call_tree)
call calltree_leave(
"ALE_offline_tracer_final()")
591 type(ocean_grid_type),
intent(in) :: G
593 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(in) :: h
595 real,
intent(in) :: threshold
600 do j = g%jsc,g%jec ;
do i = g%isc,g%iec
601 if (g%mask2dT(i,j)>0.)
then
602 if (minval(h(i,j,:)) < threshold)
then
603 write(0,*)
'check_grid: i,j=',i,j,
'h(i,j,:)=',h(i,j,:)
604 if (threshold <= 0.)
then
605 call mom_error(fatal,
"MOM_ALE, check_grid: negative thickness encountered.")
607 call mom_error(fatal,
"MOM_ALE, check_grid: too tiny thickness encountered.")
617 subroutine ale_build_grid( G, GV, regridCS, remapCS, h, tv, debug, frac_shelf_h )
618 type(ocean_grid_type),
intent(in) :: g
623 real,
dimension(SZI_(G),SZJ_(G), SZK_(GV)),
intent(inout) :: h
625 logical,
optional,
intent(in) :: debug
626 real,
dimension(:,:),
optional,
pointer :: frac_shelf_h
628 integer :: nk, i, j, k
629 real,
dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzregrid
630 real,
dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: h_new
631 logical :: show_call_tree, use_ice_shelf
633 show_call_tree = .false.
634 if (
present(debug)) show_call_tree = debug
635 if (show_call_tree)
call calltree_enter(
"ALE_build_grid(), MOM_ALE.F90")
636 use_ice_shelf = .false.
637 if (
present(frac_shelf_h))
then
638 if (
associated(frac_shelf_h)) use_ice_shelf = .true.
643 if (use_ice_shelf)
then
644 call regridding_main( remapcs, regridcs, g, gv, h, tv, h_new, dzregrid, frac_shelf_h )
646 call regridding_main( remapcs, regridcs, g, gv, h, tv, h_new, dzregrid )
652 do j = g%jsc,g%jec ;
do i = g%isc,g%iec
653 if (g%mask2dT(i,j)>0.) h(i,j,:) = h_new(i,j,:)
661 subroutine ale_regrid_accelerated(CS, G, GV, h, tv, n, u, v, OBC, Reg, dt, dzRegrid, initial)
663 type(ocean_grid_type),
intent(inout) :: g
665 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
668 integer,
intent(in) :: n
669 real,
dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
671 real,
dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
675 optional,
pointer :: reg
676 real,
optional,
intent(in) :: dt
677 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
678 optional,
intent(inout) :: dzregrid
679 logical,
optional,
intent(in) :: initial
683 integer :: i, j, k, nz
685 type(group_pass_type) :: pass_t_s_h
686 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_loc, h_orig
687 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
target :: t, s
690 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: dzinterface, dzinttotal
695 dzinttotal(:,:,:) = 0.
703 t(:,:,:) = tv%T(:,:,:)
704 s(:,:,:) = tv%S(:,:,:)
709 h_loc(:,:,:) = h(:,:,:)
710 h_orig(:,:,:) = h(:,:,:)
720 call regridding_main(cs%remapCS, cs%regridCS, g, gv, h_loc, tv_local, h, dzinterface)
721 dzinttotal(:,:,:) = dzinttotal(:,:,:) + dzinterface(:,:,:)
724 do j = g%jsc-1,g%jec+1 ;
do i = g%isc-1,g%iec+1
725 call remapping_core_h(cs%remapCS, nz, h_orig(i,j,:), tv%S(i,j,:), nz, h(i,j,:), tv_local%S(i,j,:))
726 call remapping_core_h(cs%remapCS, nz, h_orig(i,j,:), tv%T(i,j,:), nz, h(i,j,:), tv_local%T(i,j,:))
730 h_loc(:,:,:) = h(:,:,:)
734 call remap_all_state_vars(cs%remapCS, cs, g, gv, h_orig, h, reg, obc, dzinttotal, u, v)
737 if (
present(dzregrid)) dzregrid(:,:,:) = dzinttotal(:,:,:)
747 dxInterface, u, v, debug, dt)
749 type(
ale_cs),
intent(in) :: CS_ALE
750 type(ocean_grid_type),
intent(in) :: G
752 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(in) :: h_old
754 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(in) :: h_new
758 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
759 optional,
intent(in) :: dxInterface
761 real,
dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
762 optional,
intent(inout) :: u
763 real,
dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
764 optional,
intent(inout) :: v
765 logical,
optional,
intent(in) :: debug
766 real,
optional,
intent(in) :: dt
768 integer :: i, j, k, m
770 real,
dimension(GV%ke+1) :: dx
771 real,
dimension(GV%ke) :: h1, u_column
772 real,
dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: work_conc
773 real,
dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: work_cont
774 real,
dimension(SZI_(G), SZJ_(G)) :: work_2d
777 real,
dimension(GV%ke) :: h2
778 real :: h_neglect, h_neglect_edge
779 logical :: show_call_tree
782 show_call_tree = .false.
783 if (
present(debug)) show_call_tree = debug
784 if (show_call_tree)
call calltree_enter(
"remap_all_state_vars(), MOM_ALE.F90")
788 if ( .not.
present(dxinterface) .and. (cs_ale%remap_uv_using_old_alg .and. (
present(u) .or.
present(v))) )
then
789 call mom_error(fatal,
"remap_all_state_vars: dxInterface must be present if using old algorithm "// &
790 "and u/v are to be remapped")
794 if (gv%Boussinesq)
then
795 h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
797 h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
803 ntr = 0 ;
if (
associated(reg)) ntr = reg%ntr
805 if (
present(dt))
then
807 work_conc(:,:,:) = 0.0
808 work_cont(:,:,:) = 0.0
813 if (show_call_tree)
call calltree_waypoint(
"remapping tracers (remap_all_state_vars)")
817 do j = g%jsc,g%jec ;
do i = g%isc,g%iec ;
if (g%mask2dT(i,j)>0.)
then
821 call remapping_core_h(cs_remapping, nz, h1, tr%t(i,j,:), nz, h2, &
822 u_column, h_neglect, h_neglect_edge)
825 if (
present(dt))
then
826 if (tr%id_remap_conc > 0)
then
828 work_conc(i,j,k) = (u_column(k) - tr%t(i,j,k)) * idt
831 if (tr%id_remap_cont > 0 .or. tr%id_remap_cont_2d > 0)
then
833 work_cont(i,j,k) = (u_column(k)*h2(k) - tr%t(i,j,k)*h1(k)) * idt
838 tr%t(i,j,:) = u_column(:)
839 endif ;
enddo ;
enddo
842 if (
present(dt))
then
843 if (tr%id_remap_conc > 0)
then
844 call post_data(tr%id_remap_conc, work_conc, cs_ale%diag)
846 if (tr%id_remap_cont > 0)
then
847 call post_data(tr%id_remap_cont, work_cont, cs_ale%diag)
849 if (tr%id_remap_cont_2d > 0)
then
850 do j = g%jsc,g%jec ;
do i = g%isc,g%iec
853 work_2d(i,j) = work_2d(i,j) + work_cont(i,j,k)
856 call post_data(tr%id_remap_cont_2d, work_2d, cs_ale%diag)
863 if (show_call_tree)
call calltree_waypoint(
"tracers remapped (remap_all_state_vars)")
866 if (
present(u) )
then
868 do j = g%jsc,g%jec ;
do i = g%iscB,g%iecB ;
if (g%mask2dCu(i,j)>0.)
then
870 h1(:) = 0.5 * ( h_old(i,j,:) + h_old(i+1,j,:) )
871 if (cs_ale%remap_uv_using_old_alg)
then
872 dx(:) = 0.5 * ( dxinterface(i,j,:) + dxinterface(i+1,j,:) )
874 h2(k) = max( 0., h1(k) + ( dx(k+1) - dx(k) ) )
877 h2(:) = 0.5 * ( h_new(i,j,:) + h_new(i+1,j,:) )
879 if (
associated(obc))
then
880 if (obc%segnum_u(i,j) .ne. 0)
then
881 if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e)
then
885 h1(:) = h_old(i+1,j,:)
886 h2(:) = h_new(i+1,j,:)
890 call remapping_core_h(cs_remapping, nz, h1, u(i,j,:), nz, h2, &
891 u_column, h_neglect, h_neglect_edge)
892 u(i,j,:) = u_column(:)
893 endif ;
enddo ;
enddo
899 if (
present(v) )
then
901 do j = g%jscB,g%jecB ;
do i = g%isc,g%iec ;
if (g%mask2dCv(i,j)>0.)
then
903 h1(:) = 0.5 * ( h_old(i,j,:) + h_old(i,j+1,:) )
904 if (cs_ale%remap_uv_using_old_alg)
then
905 dx(:) = 0.5 * ( dxinterface(i,j,:) + dxinterface(i,j+1,:) )
907 h2(k) = max( 0., h1(k) + ( dx(k+1) - dx(k) ) )
910 h2(:) = 0.5 * ( h_new(i,j,:) + h_new(i,j+1,:) )
912 if (
associated(obc))
then
913 if (obc%segnum_v(i,j) .ne. 0)
then
918 h1(:) = h_old(i,j+1,:)
919 h2(:) = h_new(i,j+1,:)
923 call remapping_core_h(cs_remapping, nz, h1, v(i,j,:), nz, h2, &
924 u_column, h_neglect, h_neglect_edge)
925 v(i,j,:) = u_column(:)
926 endif ;
enddo ;
enddo
929 if (cs_ale%id_vert_remap_h > 0)
call post_data(cs_ale%id_vert_remap_h, h_old, cs_ale%diag)
930 if ((cs_ale%id_vert_remap_h_tendency > 0) .and.
present(dt))
then
931 do k = 1, nz ;
do j = g%jsc,g%jec ;
do i = g%isc,g%iec
932 work_cont(i,j,k) = (h_new(i,j,k) - h_old(i,j,k))*idt
933 enddo ;
enddo ;
enddo
934 call post_data(cs_ale%id_vert_remap_h_tendency, work_cont, cs_ale%diag)
945 subroutine ale_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_cells, old_remap )
947 type(ocean_grid_type),
intent(in) :: g
949 integer,
intent(in) :: nk_src
950 real,
dimension(SZI_(G),SZJ_(G),nk_src),
intent(in) :: h_src
952 real,
dimension(SZI_(G),SZJ_(G),nk_src),
intent(in) :: s_src
953 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(in) :: h_dst
955 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(inout) :: s_dst
956 logical,
optional,
intent(in) :: all_cells
959 logical,
optional,
intent(in) :: old_remap
962 integer :: i, j, k, n_points
964 real :: h_neglect, h_neglect_edge
965 logical :: ignore_vanished_layers, use_remapping_core_w
967 ignore_vanished_layers = .false.
968 if (
present(all_cells)) ignore_vanished_layers = .not. all_cells
969 use_remapping_core_w = .false.
970 if (
present(old_remap)) use_remapping_core_w = old_remap
974 if (gv%Boussinesq)
then
975 h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
977 h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
981 do j = g%jsc,g%jec ;
do i = g%isc,g%iec
982 if (g%mask2dT(i,j) > 0.)
then
983 if (ignore_vanished_layers)
then
986 if (h_src(i,j,k)>0.) n_points = n_points + 1
990 if (use_remapping_core_w)
then
991 call dzfromh1h2( n_points, h_src(i,j,1:n_points), gv%ke, h_dst(i,j,:), dx )
992 call remapping_core_w(cs, n_points, h_src(i,j,1:n_points), s_src(i,j,1:n_points), &
993 gv%ke, dx, s_dst(i,j,:), h_neglect, h_neglect_edge)
995 call remapping_core_h(cs, n_points, h_src(i,j,1:n_points), s_src(i,j,1:n_points), &
996 gv%ke, h_dst(i,j,:), s_dst(i,j,:), h_neglect, h_neglect_edge)
1012 type(ocean_grid_type),
intent(in) :: g
1014 type(
ale_cs),
intent(inout) :: cs
1015 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1016 intent(inout) :: s_t
1017 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1018 intent(inout) :: s_b
1019 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1020 intent(inout) :: t_t
1021 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1022 intent(inout) :: t_b
1024 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1026 logical,
intent(in) :: bdry_extrap
1033 real,
dimension(CS%nk,2) :: ppol_e
1034 real,
dimension(CS%nk,2) :: ppol_coefs
1038 if (gv%Boussinesq)
then
1039 h_neglect = gv%m_to_H*1.0e-30
1041 h_neglect = gv%kg_m2_to_H*1.0e-30
1046 do j = g%jsc-1,g%jec+1 ;
do i = g%isc-1,g%iec+1
1049 tmp(:) = tv%S(i,j,:)
1052 ppol_coefs(:,:) = 0.0
1058 s_t(i,j,k) = ppol_e(k,1)
1059 s_b(i,j,k) = ppol_e(k,2)
1064 ppol_coefs(:,:) = 0.0
1065 tmp(:) = tv%T(i,j,:)
1071 t_t(i,j,k) = ppol_e(k,1)
1072 t_b(i,j,k) = ppol_e(k,2)
1086 type(ocean_grid_type),
intent(in) :: g
1088 type(
ale_cs),
intent(inout) :: cs
1089 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1090 intent(inout) :: s_t
1091 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1092 intent(inout) :: s_b
1093 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1094 intent(inout) :: t_t
1095 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1096 intent(inout) :: t_b
1098 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1100 logical,
intent(in) :: bdry_extrap
1107 real,
dimension(CS%nk,2) :: &
1109 real,
dimension(CS%nk,3) :: &
1111 real :: h_neglect, h_neglect_edge
1114 if (gv%Boussinesq)
then
1115 h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
1117 h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
1122 do j = g%jsc-1,g%jec+1 ;
do i = g%isc-1,g%iec+1
1126 tmp(:) = tv%S(i,j,:)
1130 ppol_coefs(:,:) = 0.0
1133 answers_2018=cs%answers_2018 )
1139 s_t(i,j,k) = ppol_e(k,1)
1140 s_b(i,j,k) = ppol_e(k,2)
1145 ppol_coefs(:,:) = 0.0
1146 tmp(:) = tv%T(i,j,:)
1149 answers_2018=cs%answers_2018 )
1155 t_t(i,j,k) = ppol_e(k,1)
1156 t_b(i,j,k) = ppol_e(k,2)
1168 real,
intent(in) :: max_depth
1170 character(len=*),
intent(in) :: mdl
1173 character(len=30) :: coord_mode
1175 call get_param(param_file, mdl,
"REGRIDDING_COORDINATE_MODE", coord_mode, &
1176 "Coordinate mode for vertical regridding. "//&
1177 "Choose among the following possibilities: "//&
1178 trim(regriddingcoordinatemodedoc), &
1179 default=default_coordinate_mode, fail_if_missing=.true.)
1181 call initialize_regridding(regridcs, gv, us, max_depth, param_file, mdl, coord_mode,
'',
'')
1190 ale_getcoordinate(:) = getcoordinateinterfaces( cs%regridCS, undo_scaling=.true. )
1216 real,
intent(in) :: dt
1217 type(
ale_cs),
pointer :: cs
1221 if (
associated(cs))
then
1223 if (cs%regrid_time_scale > 0.0)
then
1224 w = cs%regrid_time_scale / (cs%regrid_time_scale + dt)
1226 call set_regrid_params(cs%regridCS, old_grid_weight=w)
1241 gv%sInterface(1:nk+1) = getcoordinateinterfaces( cs%regridCS, undo_scaling=.true. )
1242 gv%sLayer(1:nk) = 0.5*( gv%sInterface(1:nk) + gv%sInterface(2:nk+1) )
1243 gv%zAxisUnits = getcoordinateunits( cs%regridCS )
1244 gv%zAxisLongName = getcoordinateshortname( cs%regridCS )
1257 character(len=*),
intent(in) :: directory
1259 character(len=240) :: filepath
1261 type(fieldtype) :: fields(2)
1263 real :: ds(gv%ke), dsi(gv%ke+1)
1265 filepath = trim(directory) // trim(
"Vertical_coordinate")
1266 ds(:) = getcoordinateresolution( cs%regridCS, undo_scaling=.true. )
1268 dsi(2:gv%ke) = 0.5*( ds(1:gv%ke-1) + ds(2:gv%ke) )
1269 dsi(gv%ke+1) = 0.5*ds(gv%ke)
1271 vars(1) = var_desc(
'ds', getcoordinateunits( cs%regridCS ), &
1272 'Layer Coordinate Thickness',
'1',
'L',
'1')
1273 vars(2) = var_desc(
'ds_interface', getcoordinateunits( cs%regridCS ), &
1274 'Layer Center Coordinate Separation',
'1',
'i',
'1')
1276 call create_file(unit, trim(filepath), vars, 2, fields, single_file, gv=gv)
1277 call write_field(unit, fields(1), ds)
1278 call write_field(unit, fields(2), dsi)
1279 call close_file(unit)
1287 type(ocean_grid_type),
intent(in) :: g
1289 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(out) :: h
1294 do j = g%jsd,g%jed ;
do i = g%isd,g%ied