7 use mom_cpu_clock,
only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
27 implicit none ;
private
29 #include <MOM_memory.h>
42 real,
pointer,
dimension(:,:) :: u_shelf => null()
44 real,
pointer,
dimension(:,:) :: v_shelf => null()
47 real,
pointer,
dimension(:,:) :: u_face_mask => null()
55 real,
pointer,
dimension(:,:) :: v_face_mask => null()
57 real,
pointer,
dimension(:,:) :: u_face_mask_bdry => null()
58 real,
pointer,
dimension(:,:) :: v_face_mask_bdry => null()
59 real,
pointer,
dimension(:,:) :: u_flux_bdry_val => null()
61 real,
pointer,
dimension(:,:) :: v_flux_bdry_val => null()
64 real,
pointer,
dimension(:,:) :: umask => null()
67 real,
pointer,
dimension(:,:) :: vmask => null()
70 real,
pointer,
dimension(:,:) :: calve_mask => null()
72 real,
pointer,
dimension(:,:) :: t_shelf => null()
74 real,
pointer,
dimension(:,:) :: tmask => null()
75 real,
pointer,
dimension(:,:) :: ice_visc => null()
76 real,
pointer,
dimension(:,:) :: thickness_bdry_val => null()
77 real,
pointer,
dimension(:,:) :: u_bdry_val => null()
78 real,
pointer,
dimension(:,:) :: v_bdry_val => null()
79 real,
pointer,
dimension(:,:) :: h_bdry_val => null()
80 real,
pointer,
dimension(:,:) :: t_bdry_val => null()
82 real,
pointer,
dimension(:,:) :: taub_beta_eff => null()
85 real,
pointer,
dimension(:,:) :: od_rt => null()
86 real,
pointer,
dimension(:,:) :: float_frac_rt => null()
87 real,
pointer,
dimension(:,:) :: od_av => null()
88 real,
pointer,
dimension(:,:) :: float_frac => null()
91 integer :: od_rt_counter = 0
93 real :: velocity_update_time_step
99 real :: elapsed_velocity_time
104 logical :: gl_regularize
106 integer :: n_sub_regularize
119 real :: a_glen_isothermal
122 real :: c_basal_friction
124 real :: n_basal_friction
125 real :: density_ocean_avg
129 real :: thresh_float_col_depth
130 logical :: moving_shelf_front
132 real :: min_thickness_simple_calve
136 real :: nonlinear_tolerance
138 integer :: cg_max_iterations
139 integer :: nonlin_solve_err_mode
141 logical :: use_reproducing_sums
148 logical :: module_is_initialized = .false.
151 integer :: id_u_shelf = -1, id_v_shelf = -1, id_t_shelf = -1, &
152 id_float_frac = -1, id_col_thick = -1, id_od_av = -1, &
153 id_u_mask = -1, id_v_mask = -1, id_t_mask = -1
166 real,
intent(in) :: num
167 real,
intent(in) :: denom
173 elseif (num*denom <= 0)
then
184 real,
dimension(4),
intent(in) :: x
185 real,
dimension(4),
intent(in) :: y
186 real ::
quad_area, p2, q2, a2, c2, b2, d2
193 p2 = (x(4)-x(1))**2 + (y(4)-y(1))**2 ; q2 = (x(3)-x(2))**2 + (y(3)-y(2))**2
194 a2 = (x(3)-x(4))**2 + (y(3)-y(4))**2 ; c2 = (x(1)-x(2))**2 + (y(1)-y(2))**2
195 b2 = (x(2)-x(4))**2 + (y(2)-y(4))**2 ; d2 = (x(3)-x(1))**2 + (y(3)-y(1))**2
196 quad_area = .25 * sqrt(4*p2*q2-(b2+d2-a2-c2)**2)
208 logical :: shelf_mass_is_dynamic, override_shelf_movement, active_shelf_dynamics
209 character(len=40) :: mdl =
"MOM_ice_shelf_dyn"
210 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
212 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
213 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
215 if (
associated(cs))
then
216 call mom_error(fatal,
"MOM_ice_shelf_dyn.F90, register_ice_shelf_dyn_restarts: "// &
217 "called with an associated control structure.")
222 override_shelf_movement = .false. ; active_shelf_dynamics = .false.
223 call get_param(param_file, mdl,
"DYNAMIC_SHELF_MASS", shelf_mass_is_dynamic, &
224 "If true, the ice sheet mass can evolve with time.", &
225 default=.false., do_not_log=.true.)
226 if (shelf_mass_is_dynamic)
then
227 call get_param(param_file, mdl,
"OVERRIDE_SHELF_MOVEMENT", override_shelf_movement, &
228 "If true, user provided code specifies the ice-shelf "//&
229 "movement instead of the dynamic ice model.", default=.false., do_not_log=.true.)
230 active_shelf_dynamics = .not.override_shelf_movement
233 if (active_shelf_dynamics)
then
234 allocate( cs%u_shelf(isdb:iedb,jsdb:jedb) ) ; cs%u_shelf(:,:) = 0.0
235 allocate( cs%v_shelf(isdb:iedb,jsdb:jedb) ) ; cs%v_shelf(:,:) = 0.0
236 allocate( cs%t_shelf(isd:ied,jsd:jed) ) ; cs%t_shelf(:,:) = -10.0
237 allocate( cs%ice_visc(isd:ied,jsd:jed) ) ; cs%ice_visc(:,:) = 0.0
238 allocate( cs%taub_beta_eff(isd:ied,jsd:jed) ) ; cs%taub_beta_eff(:,:) = 0.0
239 allocate( cs%OD_av(isd:ied,jsd:jed) ) ; cs%OD_av(:,:) = 0.0
240 allocate( cs%float_frac(isd:ied,jsd:jed) ) ; cs%float_frac(:,:) = 0.0
244 "ice sheet/shelf u-velocity",
"m s-1", hor_grid=
'Bu')
246 "ice sheet/shelf v-velocity",
"m s-1", hor_grid=
'Bu')
248 "ice sheet/shelf vertically averaged temperature",
"deg C")
250 "Average open ocean depth in a cell",
"m")
252 "fractional degree of grounding",
"nondim")
254 "Glens law ice viscosity",
"m (seems wrong)")
256 "Coefficient of basal traction",
"m (seems wrong)")
265 type(time_type),
intent(inout) :: time
271 type(
diag_ctrl),
target,
intent(in) :: diag
272 logical,
intent(in) :: new_sim
274 logical,
optional,
intent(in) :: solo_ice_sheet_in
281 # include "version_variable.h"
282 character(len=200) :: config
283 character(len=200) :: ic_file,filename,inputdir
284 character(len=40) :: var_name
285 character(len=40) :: mdl =
"MOM_ice_shelf_dyn"
286 logical :: shelf_mass_is_dynamic, override_shelf_movement, active_shelf_dynamics
288 integer :: i, j, is, ie, js, je, isd, ied, jsd, jed, isdq, iedq, jsdq, jedq, iters
290 isdq = g%isdB ; iedq = g%iedB ; jsdq = g%jsdB ; jedq = g%jedB
291 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
293 if (.not.
associated(cs))
then
294 call mom_error(fatal,
"MOM_ice_shelf_dyn.F90, initialize_ice_shelf_dyn: "// &
295 "called with an associated control structure.")
298 if (cs%module_is_initialized)
then
299 call mom_error(warning,
"MOM_ice_shelf_dyn.F90, initialize_ice_shelf_dyn was "//&
300 "called with a control structure that has already been initialized.")
302 cs%module_is_initialized = .true.
308 call get_param(param_file, mdl,
"DEBUG", debug, default=.false.)
309 call get_param(param_file, mdl,
"DEBUG_IS", cs%debug, &
310 "If true, write verbose debugging messages for the ice shelf.", &
312 call get_param(param_file, mdl,
"DYNAMIC_SHELF_MASS", shelf_mass_is_dynamic, &
313 "If true, the ice sheet mass can evolve with time.", &
315 override_shelf_movement = .false. ; active_shelf_dynamics = .false.
316 if (shelf_mass_is_dynamic)
then
317 call get_param(param_file, mdl,
"OVERRIDE_SHELF_MOVEMENT", override_shelf_movement, &
318 "If true, user provided code specifies the ice-shelf "//&
319 "movement instead of the dynamic ice model.", default=.false., do_not_log=.true.)
320 active_shelf_dynamics = .not.override_shelf_movement
322 call get_param(param_file, mdl,
"GROUNDING_LINE_INTERPOLATE", cs%GL_regularize, &
323 "If true, regularize the floatation condition at the "//&
324 "grounding line as in Goldberg Holland Schoof 2009.", default=.false.)
325 call get_param(param_file, mdl,
"GROUNDING_LINE_INTERP_SUBGRID_N", cs%n_sub_regularize, &
326 "The number of sub-partitions of each cell over which to "//&
327 "integrate for the interpolated grounding line. Each cell "//&
328 "is divided into NxN equally-sized rectangles, over which the "//&
329 "basal contribution is integrated by iterative quadrature.", &
331 call get_param(param_file, mdl,
"GROUNDING_LINE_COUPLE", cs%GL_couple, &
332 "If true, let the floatation condition be determined by "//&
333 "ocean column thickness. This means that update_OD_ffrac "//&
334 "will be called. GL_REGULARIZE and GL_COUPLE are exclusive.", &
335 default=.false., do_not_log=cs%GL_regularize)
336 if (cs%GL_regularize) cs%GL_couple = .false.
337 if (cs%GL_regularize .and. (cs%n_sub_regularize == 0))
call mom_error (fatal, &
338 "GROUNDING_LINE_INTERP_SUBGRID_N must be a positive integer if GL regularization is used")
339 call get_param(param_file, mdl,
"ICE_SHELF_CFL_FACTOR", cs%CFL_factor, &
340 "A factor used to limit timestep as CFL_FACTOR * min (\Delta x / u). "//&
341 "This is only used with an ice-only model.", default=0.25)
343 call get_param(param_file, mdl,
"RHO_0", cs%density_ocean_avg, &
344 "avg ocean density used in floatation cond", &
345 units=
"kg m-3", default=1035.)
346 if (active_shelf_dynamics)
then
347 call get_param(param_file, mdl,
"ICE_VELOCITY_TIMESTEP", cs%velocity_update_time_step, &
348 "seconds between ice velocity calcs", units=
"s", &
349 fail_if_missing=.true.)
350 call get_param(param_file, mdl,
"G_EARTH", cs%g_Earth, &
351 "The gravitational acceleration of the Earth.", &
352 units=
"m s-2", default = 9.80)
354 call get_param(param_file, mdl,
"A_GLEN_ISOTHERM", cs%A_glen_isothermal, &
355 "Ice viscosity parameter in Glen's Law", &
356 units=
"Pa -1/3 a", default=9.461e-18)
357 call get_param(param_file, mdl,
"GLEN_EXPONENT", cs%n_glen, &
358 "nonlinearity exponent in Glen's Law", &
359 units=
"none", default=3.)
360 call get_param(param_file, mdl,
"MIN_STRAIN_RATE_GLEN", cs%eps_glen_min, &
361 "min. strain rate to avoid infinite Glen's law viscosity", &
362 units=
"a-1", default=1.e-12)
363 call get_param(param_file, mdl,
"BASAL_FRICTION_COEFF", cs%C_basal_friction, &
364 "ceofficient in sliding law \tau_b = C u^(n_basal_friction)", &
365 units=
"Pa (m-a)-(n_basal_friction)", fail_if_missing=.true.)
366 call get_param(param_file, mdl,
"BASAL_FRICTION_EXP", cs%n_basal_friction, &
367 "exponent in sliding law \tau_b = C u^(m_slide)", &
368 units=
"none", fail_if_missing=.true.)
369 call get_param(param_file, mdl,
"DENSITY_ICE", cs%density_ice, &
370 "A typical density of ice.", units=
"kg m-3", default=917.0)
371 call get_param(param_file, mdl,
"CONJUGATE_GRADIENT_TOLERANCE", cs%cg_tolerance, &
372 "tolerance in CG solver, relative to initial residual", default=1.e-6)
373 call get_param(param_file, mdl,
"ICE_NONLINEAR_TOLERANCE", cs%nonlinear_tolerance, &
374 "nonlin tolerance in iterative velocity solve",default=1.e-6)
375 call get_param(param_file, mdl,
"CONJUGATE_GRADIENT_MAXIT", cs%cg_max_iterations, &
376 "max iteratiions in CG solver", default=2000)
377 call get_param(param_file, mdl,
"THRESH_FLOAT_COL_DEPTH", cs%thresh_float_col_depth, &
378 "min ocean thickness to consider ice *floating*; "//&
379 "will only be important with use of tides", &
380 units=
"m", default=1.e-3, scale=us%m_to_Z)
381 call get_param(param_file, mdl,
"NONLIN_SOLVE_ERR_MODE", cs%nonlin_solve_err_mode, &
382 "Choose whether nonlin error in vel solve is based on nonlinear "//&
383 "residual (1) or relative change since last iteration (2)", default=1)
384 call get_param(param_file, mdl,
"SHELF_DYN_REPRODUCING_SUMS", cs%use_reproducing_sums, &
385 "If true, use the reproducing extended-fixed-point sums in "//&
386 "the ice shelf dynamics solvers.", default=.true.)
388 call get_param(param_file, mdl,
"SHELF_MOVING_FRONT", cs%moving_shelf_front, &
389 "Specify whether to advance shelf front (and calve).", &
391 call get_param(param_file, mdl,
"CALVE_TO_MASK", cs%calve_to_mask, &
392 "If true, do not allow an ice shelf where prohibited by a mask.", &
395 call get_param(param_file, mdl,
"MIN_THICKNESS_SIMPLE_CALVE", &
396 cs%min_thickness_simple_calve, &
397 "Min thickness rule for the VERY simple calving law",&
398 units=
"m", default=0.0, scale=us%m_to_Z)
404 if (active_shelf_dynamics)
then
406 allocate( cs%u_bdry_val(isdq:iedq,jsdq:jedq) ) ; cs%u_bdry_val(:,:) = 0.0
407 allocate( cs%v_bdry_val(isdq:iedq,jsdq:jedq) ) ; cs%v_bdry_val(:,:) = 0.0
408 allocate( cs%t_bdry_val(isd:ied,jsd:jed) ) ; cs%t_bdry_val(:,:) = -15.0
409 allocate( cs%h_bdry_val(isd:ied,jsd:jed) ) ; cs%h_bdry_val(:,:) = 0.0
410 allocate( cs%thickness_bdry_val(isd:ied,jsd:jed) ) ; cs%thickness_bdry_val(:,:) = 0.0
411 allocate( cs%u_face_mask(isdq:iedq,jsd:jed) ) ; cs%u_face_mask(:,:) = 0.0
412 allocate( cs%v_face_mask(isd:ied,jsdq:jedq) ) ; cs%v_face_mask(:,:) = 0.0
413 allocate( cs%u_face_mask_bdry(isdq:iedq,jsd:jed) ) ; cs%u_face_mask_bdry(:,:) = -2.0
414 allocate( cs%v_face_mask_bdry(isd:ied,jsdq:jedq) ) ; cs%v_face_mask_bdry(:,:) = -2.0
415 allocate( cs%u_flux_bdry_val(isdq:iedq,jsd:jed) ) ; cs%u_flux_bdry_val(:,:) = 0.0
416 allocate( cs%v_flux_bdry_val(isd:ied,jsdq:jedq) ) ; cs%v_flux_bdry_val(:,:) = 0.0
417 allocate( cs%umask(isdq:iedq,jsdq:jedq) ) ; cs%umask(:,:) = -1.0
418 allocate( cs%vmask(isdq:iedq,jsdq:jedq) ) ; cs%vmask(:,:) = -1.0
419 allocate( cs%tmask(isdq:iedq,jsdq:jedq) ) ; cs%tmask(:,:) = -1.0
422 allocate( cs%OD_rt(isd:ied,jsd:jed) ) ; cs%OD_rt(:,:) = 0.0
423 allocate( cs%float_frac_rt(isd:ied,jsd:jed) ) ; cs%float_frac_rt(:,:) = 0.0
425 if (cs%calve_to_mask)
then
426 allocate( cs%calve_mask(isd:ied,jsd:jed) ) ; cs%calve_mask(:,:) = 0.0
429 cs%elapsed_velocity_time = 0.0
435 if (active_shelf_dynamics .and. .not.new_sim)
then
436 if ((us%m_to_Z_restart /= 0.0) .and. (us%m_to_Z_restart /= us%m_to_Z))
then
437 z_rescale = us%m_to_Z / us%m_to_Z_restart
438 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
439 cs%OD_av(i,j) = z_rescale * cs%OD_av(i,j)
448 if (.not. g%symmetric)
then
449 do j=g%jsd,g%jed ;
do i=g%isd,g%ied
450 if (((i+g%idg_offset) == (g%domain%nihalo+1)).and.(cs%u_face_mask(i-1,j) == 3))
then
451 cs%u_shelf(i-1,j-1) = cs%u_bdry_val(i-1,j-1)
452 cs%u_shelf(i-1,j) = cs%u_bdry_val(i-1,j)
454 if (((j+g%jdg_offset) == (g%domain%njhalo+1)).and.(cs%v_face_mask(i,j-1) == 3))
then
455 cs%u_shelf(i-1,j-1) = cs%u_bdry_val(i-1,j-1)
456 cs%u_shelf(i,j-1) = cs%u_bdry_val(i,j-1)
462 call pass_var(cs%float_frac,g%domain)
464 call pass_var(cs%taub_beta_eff,g%domain)
468 if (active_shelf_dynamics)
then
470 if (cs%calve_to_mask)
then
471 call mom_mesg(
" MOM_ice_shelf.F90, initialize_ice_shelf: reading calving_mask")
473 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
474 inputdir = slasher(inputdir)
475 call get_param(param_file, mdl,
"CALVING_MASK_FILE", ic_file, &
476 "The file with a mask for where calving might occur.", &
477 default=
"ice_shelf_h.nc")
478 call get_param(param_file, mdl,
"CALVING_MASK_VARNAME", var_name, &
479 "The variable to use in masking calving.", &
480 default=
"area_shelf_h")
482 filename = trim(inputdir)//trim(ic_file)
483 call log_param(param_file, mdl,
"INPUTDIR/CALVING_MASK_FILE", filename)
485 " calving mask file: Unable to open "//trim(filename))
487 call mom_read_data(filename,trim(var_name),cs%calve_mask,g%Domain)
488 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
489 if (cs%calve_mask(i,j) > 0.0) cs%calve_mask(i,j) = 1.0
491 call pass_var(cs%calve_mask,g%domain)
497 call mom_mesg(
"MOM_ice_shelf.F90, initialize_ice_shelf: initialize ice velocity.")
501 if (cs%id_u_shelf > 0)
call post_data(cs%id_u_shelf,cs%u_shelf,cs%diag)
502 if (cs%id_v_shelf > 0)
call post_data(cs%id_v_shelf,cs%v_shelf,cs%diag)
506 cs%id_u_shelf = register_diag_field(
'ocean_model',
'u_shelf',cs%diag%axesCu1, time, &
507 'x-velocity of ice',
'm yr-1')
508 cs%id_v_shelf = register_diag_field(
'ocean_model',
'v_shelf',cs%diag%axesCv1, time, &
509 'y-velocity of ice',
'm yr-1')
510 cs%id_u_mask = register_diag_field(
'ocean_model',
'u_mask',cs%diag%axesCu1, time, &
511 'mask for u-nodes',
'none')
512 cs%id_v_mask = register_diag_field(
'ocean_model',
'v_mask',cs%diag%axesCv1, time, &
513 'mask for v-nodes',
'none')
516 cs%id_float_frac = register_diag_field(
'ocean_model',
'ice_float_frac',cs%diag%axesT1, time, &
517 'fraction of cell that is floating (sort of)',
'none')
518 cs%id_col_thick = register_diag_field(
'ocean_model',
'col_thick',cs%diag%axesT1, time, &
519 'ocean column thickness passed to ice model',
'm', conversion=us%Z_to_m)
520 cs%id_OD_av = register_diag_field(
'ocean_model',
'OD_av',cs%diag%axesT1, time, &
521 'intermediate ocean column thickness passed to ice model',
'm', conversion=us%Z_to_m)
530 cs%id_t_shelf = register_diag_field(
'ocean_model',
't_shelf',cs%diag%axesT1, time, &
532 cs%id_t_mask = register_diag_field(
'ocean_model',
'tmask',cs%diag%axesT1, time, &
533 'mask for T-nodes',
'none')
545 type(time_type),
intent(in) :: Time
547 integer :: i, j, iters, isd, ied, jsd, jed
548 real :: rhoi_rhow, OD
549 type(time_type) :: dummy_time
551 rhoi_rhow = cs%density_ice / cs%density_ocean_avg
552 dummy_time = set_time(0,0)
553 isd=g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
557 od = g%bathyT(i,j) - rhoi_rhow * iss%h_shelf(i,j)
561 cs%float_frac(i,j) = 0.
564 cs%float_frac(i,j) = 1.
582 real :: ratio, min_ratio
583 real :: local_u_max, local_v_max
587 do j=g%jsc,g%jec ;
do i=g%isc,g%iec ;
if (iss%hmask(i,j) == 1.0)
then
588 local_u_max = max(abs(cs%u_shelf(i,j)), abs(cs%u_shelf(i+1,j+1)), &
589 abs(cs%u_shelf(i+1,j)), abs(cs%u_shelf(i,j+1)))
590 local_v_max = max(abs(cs%v_shelf(i,j)), abs(cs%v_shelf(i+1,j+1)), &
591 abs(cs%v_shelf(i+1,j)), abs(cs%v_shelf(i,j+1)))
594 ratio = g%US%L_to_m**2*min(g%areaT(i,j) / (local_u_max + 1.0e-12), &
595 g%areaT(i,j) / (local_v_max + 1.0e-12))
596 min_ratio = min(min_ratio, ratio)
597 endif ;
enddo ;
enddo
599 call min_across_pes(min_ratio)
608 subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled_grounding, must_update_vel)
614 real,
intent(in) :: time_step
615 type(time_type),
intent(in) :: time
616 real,
dimension(SZDI_(G),SZDJ_(G)), &
617 optional,
intent(in) :: ocean_mass
619 logical,
optional,
intent(in) :: coupled_grounding
621 logical,
optional,
intent(in) :: must_update_vel
624 logical :: update_ice_vel, coupled_gl
626 update_ice_vel = .false.
627 if (
present(must_update_vel)) update_ice_vel = must_update_vel
630 if (
present(ocean_mass) .and.
present(coupled_grounding)) coupled_gl = coupled_grounding
633 cs%elapsed_velocity_time = cs%elapsed_velocity_time + time_step
634 if (cs%elapsed_velocity_time >= cs%velocity_update_time_step) update_ice_vel = .true.
638 elseif (update_ice_vel)
then
642 if (update_ice_vel)
then
646 call ice_shelf_temp(cs, iss, g, us, time_step, iss%water_flux, time)
648 if (update_ice_vel)
then
650 if (cs%id_col_thick > 0)
call post_data(cs%id_col_thick, cs%OD_av, cs%diag)
651 if (cs%id_u_shelf > 0)
call post_data(cs%id_u_shelf,cs%u_shelf,cs%diag)
652 if (cs%id_v_shelf > 0)
call post_data(cs%id_v_shelf,cs%v_shelf,cs%diag)
653 if (cs%id_t_shelf > 0)
call post_data(cs%id_t_shelf,cs%t_shelf,cs%diag)
654 if (cs%id_float_frac > 0)
call post_data(cs%id_float_frac,cs%float_frac,cs%diag)
655 if (cs%id_OD_av >0)
call post_data(cs%id_OD_av, cs%OD_av,cs%diag)
657 if (cs%id_u_mask > 0)
call post_data(cs%id_u_mask,cs%umask,cs%diag)
658 if (cs%id_v_mask > 0)
call post_data(cs%id_v_mask,cs%vmask,cs%diag)
659 if (cs%id_t_mask > 0)
call post_data(cs%id_t_mask,cs%tmask,cs%diag)
663 cs%elapsed_velocity_time = 0.0
676 real,
intent(in) :: time_step
677 type(time_type),
intent(in) :: Time
711 real,
dimension(SZDI_(G),SZDJ_(G)) :: h_after_uflux, h_after_vflux
712 real,
dimension(SZDI_(G),SZDJ_(G),4) :: flux_enter
713 integer :: isd, ied, jsd, jed, i, j, isc, iec, jsc, jec
719 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
720 isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
721 flux_enter(:,:,:) = 0.0
723 h_after_uflux(:,:) = 0.0
724 h_after_vflux(:,:) = 0.0
727 do j=jsd,jed ;
do i=isd,ied ;
if (cs%thickness_bdry_val(i,j) /= 0.0)
then
728 iss%h_shelf(i,j) = cs%thickness_bdry_val(i,j)
729 endif ;
enddo ;
enddo
747 if (iss%hmask(i,j) == 1) iss%h_shelf(i,j) = h_after_vflux(i,j)
751 if (cs%moving_shelf_front)
then
753 if (cs%min_thickness_simple_calve > 0.0)
then
755 cs%min_thickness_simple_calve)
757 if (cs%calve_to_mask)
then
758 call calve_to_mask(g, iss%h_shelf, iss%area_shelf_h, iss%hmask, cs%calve_mask)
778 real,
dimension(SZDIB_(G),SZDJB_(G)), &
780 real,
dimension(SZDIB_(G),SZDJB_(G)), &
782 integer,
intent(out) :: iters
783 type(time_type),
intent(in) :: Time
785 real,
dimension(SZDIB_(G),SZDJB_(G)) :: TAUDX, TAUDY, u_prev_iterate, v_prev_iterate, &
786 u_bdry_cont, v_bdry_cont, Au, Av, err_u, err_v, &
788 real,
dimension(SZDIB_(G),SZDJB_(G)) :: H_node
789 real,
dimension(SZDI_(G),SZDJ_(G)) :: float_cond
791 character(len=160) :: mesg
792 integer :: conv_flag, i, j, k,l, iter
793 integer :: isdq, iedq, jsdq, jedq, isd, ied, jsd, jed, isumstart, jsumstart, nodefloat, nsub
794 real :: err_max, err_tempu, err_tempv, err_init, area, max_vel, tempu, tempv, rhoi_rhow
795 real,
pointer,
dimension(:,:,:,:) :: Phi => null()
796 real,
pointer,
dimension(:,:,:,:,:,:) :: Phisub => null()
797 real,
dimension(8,4) :: Phi_temp
798 real,
dimension(2,2) :: X,Y
799 character(2) :: iternum
800 character(2) :: numproc
803 nsub = cs%n_sub_regularize
805 isdq = g%isdB ; iedq = g%iedB ; jsdq = g%jsdB ; jedq = g%jedB
806 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
807 rhoi_rhow = cs%density_ice / cs%density_ocean_avg
809 taudx(:,:) = 0.0 ; taudy(:,:) = 0.0
810 u_bdry_cont(:,:) = 0.0 ; v_bdry_cont(:,:) = 0.0
811 au(:,:) = 0.0 ; av(:,:) = 0.0
814 float_cond(:,:) = 0.0 ; h_node(:,:)=0
815 allocate(phisub(nsub,nsub,2,2,2,2)) ; phisub = 0.0
819 if (g%isc+g%idg_offset==g%isg) isumstart = g%iscB
823 if (g%jsc+g%jdg_offset==g%jsg) jsumstart = g%jscB
835 if (cs%GL_regularize)
then
844 if ((iss%hmask(i,j) == 1) .and. &
845 (rhoi_rhow * h_node(i-1+k,j-1+l) - g%bathyT(i,j) <= 0))
then
846 nodefloat = nodefloat + 1
850 if ((nodefloat > 0) .and. (nodefloat < 4))
then
851 float_cond(i,j) = 1.0
852 cs%float_frac(i,j) = 1.0
865 u_prev_iterate(:,:) = u(:,:)
866 v_prev_iterate(:,:) = v(:,:)
869 allocate(phi(isd:ied,jsd:jed,1:8,1:4)) ; phi(:,:,:,:) = 0.0
871 do j=jsd,jed ;
do i=isd,ied
872 if (((i > isd) .and. (j > jsd)))
then
873 x(:,:) = g%geoLonBu(i-1:i,j-1:j)*1000
874 y(:,:) = g%geoLatBu(i-1:i,j-1:j)*1000
876 x(2,:) = g%geoLonBu(i,j)*1000
877 x(1,:) = g%geoLonBu(i,j)*1000 - us%L_to_m*g%dxT(i,j)
878 y(:,2) = g%geoLatBu(i,j)*1000
879 y(:,1) = g%geoLatBu(i,j)*1000 - us%L_to_m*g%dyT(i,j)
883 phi(i,j,:,:) = phi_temp
888 call pass_var(cs%ice_visc, g%domain)
889 call pass_var(cs%taub_beta_eff, g%domain)
893 do j=g%jsd,g%jed ;
do i=g%isd,g%ied
894 cs%taub_beta_eff(i,j) = cs%taub_beta_eff(i,j) * cs%float_frac(i,j)
898 cs%taub_beta_eff, float_cond, &
899 rhoi_rhow, u_bdry_cont, v_bdry_cont)
901 au(:,:) = 0.0 ; av(:,:) = 0.0
903 call cg_action(au, av, u, v, phi, phisub, cs%umask, cs%vmask, iss%hmask, h_node, &
904 cs%ice_visc, float_cond, g%bathyT(:,:), cs%taub_beta_eff, g%US%L_to_m**2*g%areaT, &
905 g, g%isc-1, g%iec+1, g%jsc-1, g%jec+1, rhoi_rhow)
907 err_init = 0 ; err_tempu = 0; err_tempv = 0
908 do j=jsumstart,g%jecB
909 do i=isumstart,g%iecB
910 if (cs%umask(i,j) == 1)
then
911 err_tempu = abs(au(i,j) + u_bdry_cont(i,j) - taudx(i,j))
913 if (cs%vmask(i,j) == 1)
then
914 err_tempv = max(abs(av(i,j) + v_bdry_cont(i,j) - taudy(i,j)), err_tempu)
916 if (err_tempv >= err_init)
then
922 call max_across_pes(err_init)
924 write(mesg,*)
"ice_shelf_solve_outer: INITIAL nonlinear residual = ",err_init
927 u_last(:,:) = u(:,:) ; v_last(:,:) = v(:,:)
934 iss%hmask, conv_flag, iters, time, phi, phisub)
937 call qchksum(u,
"u shelf", g%HI, haloshift=2)
938 call qchksum(v,
"v shelf", g%HI, haloshift=2)
941 write(mesg,*)
"ice_shelf_solve_outer: linear solve done in ",iters,
" iterations"
945 call pass_var(cs%ice_visc, g%domain)
946 call pass_var(cs%taub_beta_eff, g%domain)
950 do j=g%jsd,g%jed ;
do i=g%isd,g%ied
951 cs%taub_beta_eff(i,j) = cs%taub_beta_eff(i,j) * cs%float_frac(i,j)
954 u_bdry_cont(:,:) = 0 ; v_bdry_cont(:,:) = 0
957 cs%taub_beta_eff, float_cond, &
958 rhoi_rhow, u_bdry_cont, v_bdry_cont)
960 au(:,:) = 0 ; av(:,:) = 0
962 call cg_action(au, av, u, v, phi, phisub, cs%umask, cs%vmask, iss%hmask, h_node, &
963 cs%ice_visc, float_cond, g%bathyT(:,:), cs%taub_beta_eff, g%US%L_to_m**2*g%areaT, &
964 g, g%isc-1, g%iec+1, g%jsc-1, g%jec+1, rhoi_rhow)
968 if (cs%nonlin_solve_err_mode == 1)
then
970 do j=jsumstart,g%jecB
971 do i=isumstart,g%iecB
972 if (cs%umask(i,j) == 1)
then
973 err_tempu = abs(au(i,j) + u_bdry_cont(i,j) - taudx(i,j))
975 if (cs%vmask(i,j) == 1)
then
976 err_tempv = max(abs(av(i,j) + v_bdry_cont(i,j) - taudy(i,j)), err_tempu)
978 if (err_tempv >= err_max)
then
984 call max_across_pes(err_max)
986 elseif (cs%nonlin_solve_err_mode == 2)
then
988 max_vel = 0 ; tempu = 0 ; tempv = 0
990 do j=jsumstart,g%jecB
991 do i=isumstart,g%iecB
992 if (cs%umask(i,j) == 1)
then
993 err_tempu = abs(u_last(i,j)-u(i,j))
996 if (cs%vmask(i,j) == 1)
then
997 err_tempv = max(abs(v_last(i,j)- v(i,j)), err_tempu)
998 tempv = sqrt(v(i,j)**2+tempu**2)
1000 if (err_tempv >= err_max)
then
1003 if (tempv >= max_vel)
then
1009 u_last(:,:) = u(:,:)
1010 v_last(:,:) = v(:,:)
1012 call max_across_pes(max_vel)
1013 call max_across_pes(err_max)
1018 write(mesg,*)
"ice_shelf_solve_outer: nonlinear residual = ",err_max/err_init
1021 if (err_max <= cs%nonlinear_tolerance * err_init)
then
1022 write(mesg,*)
"ice_shelf_solve_outer: exiting nonlinear solve after ",iter,
" iterations"
1035 hmask, conv_flag, iters, time, Phi, Phisub)
1040 real,
dimension(SZDIB_(G),SZDJB_(G)), &
1042 real,
dimension(SZDIB_(G),SZDJB_(G)), &
1044 real,
dimension(SZDIB_(G),SZDJB_(G)), &
1046 real,
dimension(SZDIB_(G),SZDJB_(G)), &
1048 real,
dimension(SZDIB_(G),SZDJB_(G)), &
1049 intent(in) :: H_node
1051 real,
dimension(SZDI_(G),SZDJ_(G)), &
1052 intent(in) :: float_cond
1054 real,
dimension(SZDI_(G),SZDJ_(G)), &
1057 integer,
intent(out) :: conv_flag
1059 integer,
intent(out) :: iters
1060 type(time_type),
intent(in) :: Time
1061 real,
dimension(SZDI_(G),SZDJ_(G),8,4), &
1064 real,
dimension(:,:,:,:,:,:), &
1065 intent(in) :: Phisub
1076 real,
dimension(SZDIB_(G),SZDJB_(G)) :: &
1077 Ru, Rv, Zu, Zv, DIAGu, DIAGv, RHSu, RHSv, &
1078 ubd, vbd, Au, Av, Du, Dv, &
1079 Zu_old, Zv_old, Ru_old, Rv_old, &
1081 integer :: iter, i, j, isd, ied, jsd, jed, &
1082 isc, iec, jsc, jec, is, js, ie, je, isumstart, jsumstart, &
1083 isdq, iedq, jsdq, jedq, iscq, iecq, jscq, jecq, nx_halo, ny_halo
1084 real :: tol, beta_k, alpha_k, area, dot_p1, dot_p2, resid0, cg_halo, dot_p1a, dot_p2a
1085 character(2) :: gridsize
1087 real,
dimension(8,4) :: Phi_temp
1088 real,
dimension(2,2) :: X,Y
1090 isdq = g%isdB ; iedq = g%iedB ; jsdq = g%jsdB ; jedq = g%jedB
1091 iscq = g%iscB ; iecq = g%iecB ; jscq = g%jscB ; jecq = g%jecB
1092 ny_halo = g%domain%njhalo ; nx_halo = g%domain%nihalo
1093 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1094 isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
1096 zu(:,:) = 0 ; zv(:,:) = 0 ; diagu(:,:) = 0 ; diagv(:,:) = 0
1097 ru(:,:) = 0 ; rv(:,:) = 0 ; au(:,:) = 0 ; av(:,:) = 0
1098 du(:,:) = 0 ; dv(:,:) = 0 ; ubd(:,:) = 0 ; vbd(:,:) = 0
1099 dot_p1 = 0 ; dot_p2 = 0
1103 if (g%isc+g%idg_offset==g%isg) isumstart = g%iscB
1107 if (g%jsc+g%jdg_offset==g%jsg) jsumstart = g%jscB
1110 cs%taub_beta_eff, float_cond, &
1111 cs%density_ice/cs%density_ocean_avg, ubd, vbd)
1113 rhsu(:,:) = taudx(:,:) - ubd(:,:)
1114 rhsv(:,:) = taudy(:,:) - vbd(:,:)
1120 cs%taub_beta_eff, hmask, &
1121 cs%density_ice/cs%density_ocean_avg, phisub, diagu, diagv)
1126 call cg_action(au, av, u, v, phi, phisub, cs%umask, cs%vmask, hmask, &
1127 h_node, cs%ice_visc, float_cond, g%bathyT(:,:), cs%taub_beta_eff, &
1128 g%US%L_to_m**2*g%areaT, g, isc-1, iec+1, jsc-1, jec+1, cs%density_ice/cs%density_ocean_avg)
1132 ru(:,:) = rhsu(:,:) - au(:,:) ; rv(:,:) = rhsv(:,:) - av(:,:)
1134 if (.not. cs%use_reproducing_sums)
then
1138 if (cs%umask(i,j) == 1) dot_p1 = dot_p1 + ru(i,j)**2
1139 if (cs%vmask(i,j) == 1) dot_p1 = dot_p1 + rv(i,j)**2
1143 call sum_across_pes(dot_p1)
1151 if (cs%umask(i,j) == 1) sum_vec(i,j) = ru(i,j)**2
1152 if (cs%vmask(i,j) == 1) sum_vec(i,j) = sum_vec(i,j) + rv(i,j)**2
1156 dot_p1 =
reproducing_sum( sum_vec, isumstart+(1-g%IsdB), iecq+(1-g%IsdB), &
1157 jsumstart+(1-g%JsdB), jecq+(1-g%JsdB) )
1161 resid0 = sqrt(dot_p1)
1165 if (cs%umask(i,j) == 1) zu(i,j) = ru(i,j) / diagu(i,j)
1166 if (cs%vmask(i,j) == 1) zv(i,j) = rv(i,j) / diagv(i,j)
1170 du(:,:) = zu(:,:) ; dv(:,:) = zv(:,:)
1185 do iter = 1,cs%cg_max_iterations
1192 is = isc - cg_halo ; ie = iecq + cg_halo
1193 js = jscq - cg_halo ; je = jecq + cg_halo
1195 au(:,:) = 0 ; av(:,:) = 0
1197 call cg_action(au, av, du, dv, phi, phisub, cs%umask, cs%vmask, hmask, &
1198 h_node, cs%ice_visc, float_cond, g%bathyT(:,:), cs%taub_beta_eff, &
1199 g%US%L_to_m**2*g%areaT, g, is, ie, js, je, cs%density_ice/cs%density_ocean_avg)
1203 if ( .not. cs%use_reproducing_sums)
then
1207 dot_p1 = 0 ; dot_p2 = 0
1210 if (cs%umask(i,j) == 1)
then
1211 dot_p1 = dot_p1 + zu(i,j)*ru(i,j)
1212 dot_p2 = dot_p2 + du(i,j)*au(i,j)
1214 if (cs%vmask(i,j) == 1)
then
1215 dot_p1 = dot_p1 + zv(i,j)*rv(i,j)
1216 dot_p2 = dot_p2 + dv(i,j)*av(i,j)
1220 call sum_across_pes(dot_p1) ;
call sum_across_pes(dot_p2)
1223 sum_vec(:,:) = 0.0 ; sum_vec_2(:,:) = 0.0
1227 if (cs%umask(i,j) == 1) sum_vec(i,j) = zu(i,j) * ru(i,j)
1228 if (cs%vmask(i,j) == 1) sum_vec(i,j) = sum_vec(i,j) + zv(i,j) * rv(i,j)
1230 if (cs%umask(i,j) == 1) sum_vec_2(i,j) = du(i,j) * au(i,j)
1231 if (cs%vmask(i,j) == 1) sum_vec_2(i,j) = sum_vec_2(i,j) + dv(i,j) * av(i,j)
1235 dot_p1 =
reproducing_sum( sum_vec, isumstart+(1-g%IsdB), iecq+(1-g%IsdB), &
1236 jsumstart+(1-g%JsdB), jecq+(1-g%JsdB) )
1238 dot_p2 =
reproducing_sum( sum_vec_2, isumstart+(1-g%IsdB), iecq+(1-g%IsdB), &
1239 jsumstart+(1-g%JsdB), jecq+(1-g%JsdB) )
1242 alpha_k = dot_p1/dot_p2
1246 if (cs%umask(i,j) == 1) u(i,j) = u(i,j) + alpha_k * du(i,j)
1247 if (cs%vmask(i,j) == 1) v(i,j) = v(i,j) + alpha_k * dv(i,j)
1253 if (cs%umask(i,j) == 1)
then
1254 ru_old(i,j) = ru(i,j) ; zu_old(i,j) = zu(i,j)
1256 if (cs%vmask(i,j) == 1)
then
1257 rv_old(i,j) = rv(i,j) ; zv_old(i,j) = zv(i,j)
1267 if (cs%umask(i,j) == 1) ru(i,j) = ru(i,j) - alpha_k * au(i,j)
1268 if (cs%vmask(i,j) == 1) rv(i,j) = rv(i,j) - alpha_k * av(i,j)
1275 if (cs%umask(i,j) == 1)
then
1276 zu(i,j) = ru(i,j) / diagu(i,j)
1278 if (cs%vmask(i,j) == 1)
then
1279 zv(i,j) = rv(i,j) / diagv(i,j)
1286 if (.not. cs%use_reproducing_sums)
then
1289 dot_p1 = 0 ; dot_p2 = 0
1292 if (cs%umask(i,j) == 1)
then
1293 dot_p1 = dot_p1 + zu(i,j)*ru(i,j)
1294 dot_p2 = dot_p2 + zu_old(i,j)*ru_old(i,j)
1296 if (cs%vmask(i,j) == 1)
then
1297 dot_p1 = dot_p1 + zv(i,j)*rv(i,j)
1298 dot_p2 = dot_p2 + zv_old(i,j)*rv_old(i,j)
1302 call sum_across_pes(dot_p1) ;
call sum_across_pes(dot_p2)
1307 sum_vec(:,:) = 0.0 ; sum_vec_2(:,:) = 0.0
1311 if (cs%umask(i,j) == 1) sum_vec(i,j) = zu(i,j) * ru(i,j)
1312 if (cs%vmask(i,j) == 1) sum_vec(i,j) = sum_vec(i,j) + &
1315 if (cs%umask(i,j) == 1) sum_vec_2(i,j) = zu_old(i,j) * ru_old(i,j)
1316 if (cs%vmask(i,j) == 1) sum_vec_2(i,j) = sum_vec_2(i,j) + &
1317 zv_old(i,j) * rv_old(i,j)
1322 dot_p1 =
reproducing_sum(sum_vec, isumstart+(1-g%IsdB), iecq+(1-g%IsdB), &
1323 jsumstart+(1-g%JsdB), jecq+(1-g%JsdB) )
1325 dot_p2 =
reproducing_sum(sum_vec_2, isumstart+(1-g%IsdB), iecq+(1-g%IsdB), &
1326 jsumstart+(1-g%JsdB), jecq+(1-g%JsdB) )
1330 beta_k = dot_p1/dot_p2
1338 if (cs%umask(i,j) == 1) du(i,j) = zu(i,j) + beta_k * du(i,j)
1339 if (cs%vmask(i,j) == 1) dv(i,j) = zv(i,j) + beta_k * dv(i,j)
1347 if (.not. cs%use_reproducing_sums)
then
1351 if (cs%umask(i,j) == 1)
then
1352 dot_p1 = dot_p1 + ru(i,j)**2
1354 if (cs%vmask(i,j) == 1)
then
1355 dot_p1 = dot_p1 + rv(i,j)**2
1359 call sum_across_pes(dot_p1)
1367 if (cs%umask(i,j) == 1) sum_vec(i,j) = ru(i,j)**2
1368 if (cs%vmask(i,j) == 1) sum_vec(i,j) = sum_vec(i,j) + rv(i,j)**2
1372 dot_p1 =
reproducing_sum( sum_vec, isumstart+(1-g%IsdB), iecq+(1-g%IsdB), &
1373 jsumstart+(1-g%JsdB), jecq+(1-g%JsdB) )
1376 dot_p1 = sqrt(dot_p1)
1378 if (dot_p1 <= cs%cg_tolerance * resid0)
then
1384 cg_halo = cg_halo - 1
1386 if (cg_halo == 0)
then
1398 if (cs%umask(i,j) == 3)
then
1399 u(i,j) = cs%u_bdry_val(i,j)
1400 elseif (cs%umask(i,j) == 0)
then
1404 if (cs%vmask(i,j) == 3)
then
1405 v(i,j) = cs%v_bdry_val(i,j)
1406 elseif (cs%vmask(i,j) == 0)
then
1414 if (conv_flag == 0)
then
1415 iters = cs%cg_max_iterations
1423 real,
intent(in) :: time_step
1424 real,
dimension(SZDI_(G),SZDJ_(G)), &
1425 intent(inout) :: hmask
1427 real,
dimension(SZDI_(G),SZDJ_(G)), &
1429 real,
dimension(SZDI_(G),SZDJ_(G)), &
1430 intent(inout) :: h_after_uflux
1432 real,
dimension(SZDI_(G),SZDJ_(G),4), &
1433 intent(inout) :: flux_enter
1454 integer :: i, j, is, ie, js, je, isd, ied, jsd, jed, gjed, gied
1455 integer :: i_off, j_off
1456 logical :: at_east_bdry, at_west_bdry, one_off_west_bdry, one_off_east_bdry
1457 real,
dimension(-2:2) :: stencil
1459 flux_diff_cell, phi, dxh, dyh, dxdyh
1460 character (len=1) :: debug_str
1462 is = g%isc-2 ; ie = g%iec+2 ; js = g%jsc ; je = g%jec
1463 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1464 i_off = g%idg_offset ; j_off = g%jdg_offset
1467 if (((j+j_off) <= g%domain%njglobal+g%domain%njhalo) .AND. &
1468 ((j+j_off) >= g%domain%njhalo+1))
then
1474 if (((i+i_off) <= g%domain%niglobal+g%domain%nihalo) .AND. &
1475 ((i+i_off) >= g%domain%nihalo+1))
then
1477 if (i+i_off == g%domain%nihalo+1)
then
1480 at_west_bdry=.false.
1483 if (i+i_off == g%domain%niglobal+g%domain%nihalo)
then
1486 at_east_bdry=.false.
1489 if (hmask(i,j) == 1)
then
1491 dxh = g%US%L_to_m*g%dxT(i,j) ; dyh = g%US%L_to_m*g%dyT(i,j) ; dxdyh = g%US%L_to_m**2*g%areaT(i,j)
1493 h_after_uflux(i,j) = h0(i,j)
1495 stencil(:) = h0(i-2:i+2,j)
1501 if (cs%u_face_mask(i-1,j) == 4.)
then
1503 flux_diff_cell = flux_diff_cell + dyh * time_step * cs%u_flux_bdry_val(i-1,j) / dxdyh
1508 u_face = 0.5 * (cs%u_shelf(i-1,j-1) + cs%u_shelf(i-1,j))
1510 if (u_face > 0)
then
1514 if (at_west_bdry .AND. (hmask(i-1,j) == 3))
then
1516 stencil(-1) = cs%thickness_bdry_val(i-1,j)
1517 flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step * stencil(-1) / dxdyh
1519 elseif (hmask(i-1,j) * hmask(i-2,j) == 1)
then
1520 phi =
slope_limiter(stencil(-1)-stencil(-2), stencil(0)-stencil(-1))
1521 flux_diff_cell = flux_diff_cell + abs(u_face) * dyh* time_step / dxdyh * &
1522 (stencil(-1) - phi * (stencil(-1)-stencil(0))/2)
1528 flux_diff_cell = flux_diff_cell + abs(u_face) * (dyh * time_step / dxdyh) * stencil(-1)
1532 elseif (u_face < 0)
then
1533 if (hmask(i-1,j) * hmask(i+1,j) == 1)
then
1534 phi =
slope_limiter(stencil(0)-stencil(1), stencil(-1)-stencil(0))
1535 flux_diff_cell = flux_diff_cell - abs(u_face) * dyh * time_step / dxdyh * &
1536 (stencil(0) - phi * (stencil(0)-stencil(-1))/2)
1539 flux_diff_cell = flux_diff_cell - abs(u_face) * (dyh * time_step / dxdyh) * stencil(0)
1541 if ((hmask(i-1,j) == 0) .OR. (hmask(i-1,j) == 2))
then
1542 flux_enter(i-1,j,2) = abs(u_face) * dyh * time_step * stencil(0)
1552 if (cs%u_face_mask(i+1,j) == 4.)
then
1554 flux_diff_cell = flux_diff_cell + dyh * time_step * cs%u_flux_bdry_val(i+1,j) / dxdyh
1558 u_face = 0.5 * (cs%u_shelf(i,j-1) + cs%u_shelf(i,j))
1560 if (u_face < 0)
then
1562 if (at_east_bdry .AND. (hmask(i+1,j) == 3))
then
1565 flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step * stencil(1) / dxdyh
1567 elseif (hmask(i+1,j) * hmask(i+2,j) == 1)
then
1569 phi =
slope_limiter(stencil(1)-stencil(2), stencil(0)-stencil(1))
1570 flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step / dxdyh * &
1571 (stencil(1) - phi * (stencil(1)-stencil(0))/2)
1577 flux_diff_cell = flux_diff_cell + abs(u_face) * (dyh * time_step / dxdyh) * stencil(1)
1581 elseif (u_face > 0)
then
1583 if (hmask(i-1,j) * hmask(i+1,j) == 1)
then
1585 phi =
slope_limiter(stencil(0)-stencil(-1), stencil(1)-stencil(0))
1586 flux_diff_cell = flux_diff_cell - abs(u_face) * dyh * time_step / dxdyh * &
1587 (stencil(0) - phi * (stencil(0)-stencil(1))/2)
1593 flux_diff_cell = flux_diff_cell - abs(u_face) * (dyh * time_step / dxdyh) * stencil(0)
1595 if ((hmask(i+1,j) == 0) .OR. (hmask(i+1,j) == 2))
then
1596 flux_enter(i+1,j,1) = abs(u_face) * dyh * time_step * stencil(0)
1603 h_after_uflux(i,j) = h_after_uflux(i,j) + flux_diff_cell
1607 elseif ((hmask(i,j) == 0) .OR. (hmask(i,j) == 2))
then
1609 if (at_west_bdry .AND. (hmask(i-1,j) == 3))
then
1610 u_face = 0.5 * (cs%u_shelf(i-1,j-1) + cs%u_shelf(i-1,j))
1611 flux_enter(i,j,1) = abs(u_face) * g%US%L_to_m*g%dyT(i,j) * time_step * cs%thickness_bdry_val(i-1,j)
1612 elseif (cs%u_face_mask(i-1,j) == 4.)
then
1613 flux_enter(i,j,1) = g%US%L_to_m*g%dyT(i,j) * time_step * cs%u_flux_bdry_val(i-1,j)
1616 if (at_east_bdry .AND. (hmask(i+1,j) == 3))
then
1617 u_face = 0.5 * (cs%u_shelf(i,j-1) + cs%u_shelf(i,j))
1618 flux_enter(i,j,2) = abs(u_face) * g%US%L_to_m*g%dyT(i,j) * time_step * cs%thickness_bdry_val(i+1,j)
1619 elseif (cs%u_face_mask(i+1,j) == 4.)
then
1620 flux_enter(i,j,2) = g%US%L_to_m*g%dyT(i,j) * time_step * cs%u_flux_bdry_val(i+1,j)
1623 if ((i == is) .AND. (hmask(i,j) == 0) .AND. (hmask(i-1,j) == 1))
then
1629 elseif ((i == ie) .AND. (hmask(i,j) == 0) .AND. (hmask(i+1,j) == 1))
then
1653 real,
intent(in) :: time_step
1654 real,
dimension(SZDI_(G),SZDJ_(G)), &
1655 intent(inout) :: hmask
1657 real,
dimension(SZDI_(G),SZDJ_(G)), &
1658 intent(in) :: h_after_uflux
1660 real,
dimension(SZDI_(G),SZDJ_(G)), &
1661 intent(inout) :: h_after_vflux
1663 real,
dimension(SZDI_(G),SZDJ_(G),4), &
1664 intent(inout) :: flux_enter
1685 integer :: i, j, is, ie, js, je, isd, ied, jsd, jed, gjed, gied
1686 integer :: i_off, j_off
1687 logical :: at_north_bdry, at_south_bdry, one_off_west_bdry, one_off_east_bdry
1688 real,
dimension(-2:2) :: stencil
1690 flux_diff_cell, phi, dxh, dyh, dxdyh
1691 character(len=1) :: debug_str
1693 is = g%isc ; ie = g%iec ; js = g%jsc-1 ; je = g%jec+1 ; isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1694 i_off = g%idg_offset ; j_off = g%jdg_offset
1697 if (((i+i_off) <= g%domain%niglobal+g%domain%nihalo) .AND. &
1698 ((i+i_off) >= g%domain%nihalo+1))
then
1704 if (((j+j_off) <= g%domain%njglobal+g%domain%njhalo) .AND. &
1705 ((j+j_off) >= g%domain%njhalo+1))
then
1707 if (j+j_off == g%domain%njhalo+1)
then
1708 at_south_bdry=.true.
1710 at_south_bdry=.false.
1713 if (j+j_off == g%domain%njglobal+g%domain%njhalo)
then
1714 at_north_bdry=.true.
1716 at_north_bdry=.false.
1719 if (hmask(i,j) == 1)
then
1720 dxh = g%US%L_to_m*g%dxT(i,j) ; dyh = g%US%L_to_m*g%dyT(i,j) ; dxdyh = g%US%L_to_m**2*g%areaT(i,j)
1721 h_after_vflux(i,j) = h_after_uflux(i,j)
1723 stencil(:) = h_after_uflux(i,j-2:j+2)
1728 if (cs%v_face_mask(i,j-1) == 4.)
then
1730 flux_diff_cell = flux_diff_cell + dxh * time_step * cs%v_flux_bdry_val(i,j-1) / dxdyh
1735 v_face = 0.5 * (cs%v_shelf(i-1,j-1) + cs%v_shelf(i,j-1))
1737 if (v_face > 0)
then
1741 if (at_south_bdry .AND. (hmask(i,j-1) == 3))
then
1743 flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step * stencil(-1) / dxdyh
1745 elseif (hmask(i,j-1) * hmask(i,j-2) == 1)
then
1747 phi =
slope_limiter(stencil(-1)-stencil(-2), stencil(0)-stencil(-1))
1748 flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step / dxdyh * &
1749 (stencil(-1) - phi * (stencil(-1)-stencil(0))/2)
1754 flux_diff_cell = flux_diff_cell + abs(v_face) * (dxh * time_step / dxdyh) * stencil(-1)
1757 elseif (v_face < 0)
then
1759 if (hmask(i,j-1) * hmask(i,j+1) == 1)
then
1760 phi =
slope_limiter(stencil(0)-stencil(1), stencil(-1)-stencil(0))
1761 flux_diff_cell = flux_diff_cell - abs(v_face) * dxh * time_step / dxdyh * &
1762 (stencil(0) - phi * (stencil(0)-stencil(-1))/2)
1764 flux_diff_cell = flux_diff_cell - abs(v_face) * (dxh * time_step / dxdyh) * stencil(0)
1766 if ((hmask(i,j-1) == 0) .OR. (hmask(i,j-1) == 2))
then
1767 flux_enter(i,j-1,4) = abs(v_face) * dyh * time_step * stencil(0)
1778 if (cs%v_face_mask(i,j+1) == 4.)
then
1780 flux_diff_cell = flux_diff_cell + dxh * time_step * cs%v_flux_bdry_val(i,j+1) / dxdyh
1785 v_face = 0.5 * (cs%v_shelf(i-1,j) + cs%v_shelf(i,j))
1787 if (v_face < 0)
then
1789 if (at_north_bdry .AND. (hmask(i,j+1) == 3))
then
1791 flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step * stencil(1) / dxdyh
1792 elseif (hmask(i,j+1) * hmask(i,j+2) == 1)
then
1793 phi =
slope_limiter(stencil(1)-stencil(2), stencil(0)-stencil(1))
1794 flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step / dxdyh * &
1795 (stencil(1) - phi * (stencil(1)-stencil(0))/2)
1799 flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step / dxdyh * stencil(1)
1802 elseif (v_face > 0)
then
1804 if (hmask(i,j-1) * hmask(i,j+1) == 1)
then
1805 phi =
slope_limiter(stencil(0)-stencil(-1), stencil(1)-stencil(0))
1806 flux_diff_cell = flux_diff_cell - abs(v_face) * dxh * time_step / dxdyh * &
1807 (stencil(0) - phi * (stencil(0)-stencil(1))/2)
1811 flux_diff_cell = flux_diff_cell - abs(v_face) * dxh * time_step / dxdyh * stencil(0)
1812 if ((hmask(i,j+1) == 0) .OR. (hmask(i,j+1) == 2))
then
1813 flux_enter(i,j+1,3) = abs(v_face) * dxh * time_step * stencil(0)
1821 h_after_vflux(i,j) = h_after_vflux(i,j) + flux_diff_cell
1823 elseif ((hmask(i,j) == 0) .OR. (hmask(i,j) == 2))
then
1825 if (at_south_bdry .AND. (hmask(i,j-1) == 3))
then
1826 v_face = 0.5 * (cs%u_shelf(i-1,j-1) + cs%u_shelf(i,j-1))
1827 flux_enter(i,j,3) = abs(v_face) * g%US%L_to_m*g%dxT(i,j) * time_step * cs%thickness_bdry_val(i,j-1)
1828 elseif (cs%v_face_mask(i,j-1) == 4.)
then
1829 flux_enter(i,j,3) = g%US%L_to_m*g%dxT(i,j) * time_step * cs%v_flux_bdry_val(i,j-1)
1832 if (at_north_bdry .AND. (hmask(i,j+1) == 3))
then
1833 v_face = 0.5 * (cs%u_shelf(i-1,j) + cs%u_shelf(i,j))
1834 flux_enter(i,j,4) = abs(v_face) * g%US%L_to_m*g%dxT(i,j) * time_step * cs%thickness_bdry_val(i,j+1)
1835 elseif (cs%v_face_mask(i,j+1) == 4.)
then
1836 flux_enter(i,j,4) = g%US%L_to_m*g%dxT(i,j) * time_step * cs%v_flux_bdry_val(i,j+1)
1839 if ((j == js) .AND. (hmask(i,j) == 0) .AND. (hmask(i,j-1) == 1))
then
1844 elseif ((j == je) .AND. (hmask(i,j) == 0) .AND. (hmask(i,j+1) == 1))
then
1864 real,
dimension(SZDI_(G),SZDJ_(G),4), &
1865 intent(inout) :: flux_enter
1895 integer :: i, j, isc, iec, jsc, jec, n_flux, k, l, iter_count
1896 integer :: i_off, j_off
1897 integer :: iter_flag
1899 real :: h_reference, dxh, dyh, dxdyh, rho, partial_vol, tot_flux
1900 character(len=160) :: mesg
1901 integer,
dimension(4) :: mapi, mapj, new_partial
1903 real,
dimension(SZDI_(G),SZDJ_(G),4) :: flux_enter_replace
1905 isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
1906 i_off = g%idg_offset ; j_off = g%jdg_offset
1907 rho = cs%density_ice
1908 iter_count = 0 ; iter_flag = 1
1911 mapi(1) = -1 ; mapi(2) = 1 ; mapi(3:4) = 0
1912 mapj(3) = -1 ; mapj(4) = 1 ; mapj(1:2) = 0
1914 do while (iter_flag == 1)
1918 if (iter_count > 0)
then
1919 flux_enter(:,:,:) = flux_enter_replace(:,:,:)
1921 flux_enter_replace(:,:,:) = 0.0
1923 iter_count = iter_count + 1
1929 if (((j+j_off) <= g%domain%njglobal+g%domain%njhalo) .AND. &
1930 ((j+j_off) >= g%domain%njhalo+1))
then
1934 if (((i+i_off) <= g%domain%niglobal+g%domain%nihalo) .AND. &
1935 ((i+i_off) >= g%domain%nihalo+1))
then
1942 if (flux_enter(i,j,k) > 0)
then
1944 h_reference = h_reference + iss%h_shelf(i+2*k-3,j)
1945 tot_flux = tot_flux + flux_enter(i,j,k)
1946 flux_enter(i,j,k) = 0.0
1951 if (flux_enter(i,j,k+2) > 0)
then
1953 h_reference = h_reference + iss%h_shelf(i,j+2*k-3)
1954 tot_flux = tot_flux + flux_enter(i,j,k+2)
1955 flux_enter(i,j,k+2) = 0.0
1959 if (n_flux > 0)
then
1960 dxdyh = g%US%L_to_m**2*g%areaT(i,j)
1961 h_reference = h_reference / real(n_flux)
1962 partial_vol = iss%h_shelf(i,j) * iss%area_shelf_h(i,j) + tot_flux
1964 if ((partial_vol / dxdyh) == h_reference)
then
1966 iss%h_shelf(i,j) = h_reference
1967 iss%area_shelf_h(i,j) = dxdyh
1968 elseif ((partial_vol / dxdyh) < h_reference)
then
1971 iss%area_shelf_h(i,j) = partial_vol / h_reference
1972 iss%h_shelf(i,j) = h_reference
1976 iss%area_shelf_h(i,j) = dxdyh
1978 partial_vol = partial_vol - h_reference * dxdyh
1982 n_flux = 0 ; new_partial(:) = 0
1985 if (cs%u_face_mask(i-2+k,j) == 2)
then
1987 elseif (iss%hmask(i+2*k-3,j) == 0)
then
1993 if (cs%v_face_mask(i,j-2+k) == 2)
then
1995 elseif (iss%hmask(i,j+2*k-3) == 0)
then
1997 new_partial(k+2) = 1
2001 if (n_flux == 0)
then
2002 iss%h_shelf(i,j) = h_reference + partial_vol / dxdyh
2004 iss%h_shelf(i,j) = h_reference
2007 if (new_partial(k) == 1) &
2008 flux_enter_replace(i+2*k-3,j,3-k) = partial_vol / real(n_flux)
2011 if (new_partial(k+2) == 1) &
2012 flux_enter_replace(i,j+2*k-3,5-k) = partial_vol / real(n_flux)
2028 call max_across_pes(iter_count)
2030 if (
is_root_pe() .and. (iter_count > 1))
then
2031 write(mesg,*)
"shelf_advance_front: ", iter_count,
" max iterations"
2040 real,
dimension(SZDI_(G),SZDJ_(G)), &
2041 intent(inout) :: h_shelf
2042 real,
dimension(SZDI_(G),SZDJ_(G)), &
2043 intent(inout) :: area_shelf_h
2044 real,
dimension(SZDI_(G),SZDJ_(G)), &
2045 intent(inout) :: hmask
2047 real,
intent(in) :: thickness_calve
2055 if ((h_shelf(i,j) < thickness_calve) .and. (area_shelf_h(i,j) > 0.))
then
2057 area_shelf_h(i,j) = 0.0
2065 subroutine calve_to_mask(G, h_shelf, area_shelf_h, hmask, calve_mask)
2067 real,
dimension(SZDI_(G),SZDJ_(G)), &
2068 intent(inout) :: h_shelf
2069 real,
dimension(SZDI_(G),SZDJ_(G)), &
2070 intent(inout) :: area_shelf_h
2071 real,
dimension(SZDI_(G),SZDJ_(G)), &
2072 intent(inout) :: hmask
2074 real,
dimension(SZDI_(G),SZDJ_(G)), &
2075 intent(in) :: calve_mask
2080 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
2081 if ((calve_mask(i,j) == 0.0) .and. (hmask(i,j) /= 0.0))
then
2083 area_shelf_h(i,j) = 0.0
2096 real,
dimension(SZDI_(G),SZDJ_(G)), &
2098 real,
dimension(SZDIB_(G),SZDJB_(G)), &
2099 intent(inout) :: TAUD_X
2100 real,
dimension(SZDIB_(G),SZDJB_(G)), &
2101 intent(inout) :: TAUD_Y
2114 real,
dimension(SIZE(OD,1),SIZE(OD,2)) :: S, &
2118 real :: rho, rhow, sx, sy, neumann_val, dxh, dyh, dxdyh, grav
2120 integer :: i, j, iscq, iecq, jscq, jecq, isd, jsd, is, js, iegq, jegq
2121 integer :: giec, gjec, gisc, gjsc, cnt, isc, jsc, iec, jec
2122 integer :: i_off, j_off
2124 isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
2125 iscq = g%iscB ; iecq = g%iecB ; jscq = g%jscB ; jecq = g%jecB
2126 isd = g%isd ; jsd = g%jsd
2127 iegq = g%iegB ; jegq = g%jegB
2128 gisc = g%domain%nihalo+1 ; gjsc = g%domain%njhalo+1
2129 giec = g%domain%niglobal+g%domain%nihalo ; gjec = g%domain%njglobal+g%domain%njhalo
2130 is = iscq - 1; js = jscq - 1
2131 i_off = g%idg_offset ; j_off = g%jdg_offset
2133 rho = cs%density_ice
2134 rhow = cs%density_ocean_avg
2135 grav = us%Z_to_m**2 * cs%g_Earth
2140 base(:,:) = -g%bathyT(:,:) + od(:,:)
2141 s(:,:) = base(:,:) + iss%h_shelf(:,:)
2148 dxh = us%L_to_m*g%dxT(i,j)
2149 dyh = us%L_to_m*g%dyT(i,j)
2150 dxdyh = us%L_to_m**2*g%areaT(i,j)
2152 if (iss%hmask(i,j) == 1)
then
2155 if ((i+i_off) == gisc)
then
2156 if (iss%hmask(i+1,j) == 1)
then
2157 sx = (s(i+1,j)-s(i,j))/dxh
2161 elseif ((i+i_off) == giec)
then
2162 if (iss%hmask(i-1,j) == 1)
then
2163 sx = (s(i,j)-s(i-1,j))/dxh
2168 if (iss%hmask(i+1,j) == 1)
then
2174 if (iss%hmask(i-1,j) == 1)
then
2183 sx = sx / (cnt * dxh)
2190 if ((j+j_off) == gjsc)
then
2191 if (iss%hmask(i,j+1) == 1)
then
2192 sy = (s(i,j+1)-s(i,j))/dyh
2196 elseif ((j+j_off) == gjec)
then
2197 if (iss%hmask(i,j-1) == 1)
then
2198 sy = (s(i,j)-s(i,j-1))/dyh
2203 if (iss%hmask(i,j+1) == 1)
then
2209 if (iss%hmask(i,j-1) == 1)
then
2218 sy = sy / (cnt * dyh)
2223 taud_x(i-1,j-1) = taud_x(i-1,j-1) - .25 * rho * grav * iss%h_shelf(i,j) * sx * dxdyh
2224 taud_y(i-1,j-1) = taud_y(i-1,j-1) - .25 * rho * grav * iss%h_shelf(i,j) * sy * dxdyh
2227 taud_x(i,j-1) = taud_x(i,j-1) - .25 * rho * grav * iss%h_shelf(i,j) * sx * dxdyh
2228 taud_y(i,j-1) = taud_y(i,j-1) - .25 * rho * grav * iss%h_shelf(i,j) * sy * dxdyh
2231 taud_x(i-1,j) = taud_x(i-1,j) - .25 * rho * grav * iss%h_shelf(i,j) * sx * dxdyh
2232 taud_y(i-1,j) = taud_y(i-1,j) - .25 * rho * grav * iss%h_shelf(i,j) * sy * dxdyh
2235 taud_x(i,j) = taud_x(i,j) - .25 * rho * grav * iss%h_shelf(i,j) * sx * dxdyh
2236 taud_y(i,j) = taud_y(i,j) - .25 * rho * grav * iss%h_shelf(i,j) * sy * dxdyh
2238 if (cs%float_frac(i,j) == 1)
then
2239 neumann_val = .5 * grav * (rho * iss%h_shelf(i,j)**2 - rhow * g%bathyT(i,j)**2)
2241 neumann_val = .5 * grav * (1-rho/rhow) * rho * iss%h_shelf(i,j)**2
2245 if ((cs%u_face_mask(i-1,j) == 2) .OR. (iss%hmask(i-1,j) == 0) .OR. (iss%hmask(i-1,j) == 2) )
then
2255 taud_x(i-1,j-1) = taud_x(i-1,j-1) - .5 * dyh * neumann_val
2256 taud_x(i-1,j) = taud_x(i-1,j) - .5 * dyh * neumann_val
2259 if ((cs%u_face_mask(i,j) == 2) .OR. (iss%hmask(i+1,j) == 0) .OR. (iss%hmask(i+1,j) == 2) )
then
2261 taud_x(i,j-1) = taud_x(i,j-1) + .5 * dyh * neumann_val
2262 taud_x(i,j) = taud_x(i,j) + .5 * dyh * neumann_val
2265 if ((cs%v_face_mask(i,j-1) == 2) .OR. (iss%hmask(i,j-1) == 0) .OR. (iss%hmask(i,j-1) == 2) )
then
2267 taud_y(i-1,j-1) = taud_y(i-1,j-1) - .5 * dxh * neumann_val
2268 taud_y(i,j-1) = taud_y(i,j-1) - .5 * dxh * neumann_val
2271 if ((cs%v_face_mask(i,j) == 2) .OR. (iss%hmask(i,j+1) == 0) .OR. (iss%hmask(i,j+1) == 2) )
then
2273 taud_y(i-1,j) = taud_y(i-1,j) + .5 * dxh * neumann_val
2274 taud_y(i,j) = taud_y(i,j) + .5 * dxh * neumann_val
2286 type(time_type),
intent(in) :: Time
2287 real,
dimension(SZDI_(G),SZDJ_(G)), &
2290 real,
intent(in) :: input_flux
2291 real,
intent(in) :: input_thick
2292 logical,
optional,
intent(in) :: new_sim
2301 integer :: i, j, iscq, iecq, jscq, jecq, isd, jsd, ied, jed, iegq, jegq
2302 integer :: gjec, gisc, gjsc, cnt, isc, jsc, iec, jec
2303 integer :: i_off, j_off
2304 real :: A, n, ux, uy, vx, vy, eps_min, domain_width
2307 isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
2309 isd = g%isd ; jsd = g%jsd ; ied = g%ied ; jed = g%jed
2311 i_off = g%idg_offset ; j_off = g%jdg_offset
2313 domain_width = g%len_lat
2320 if (hmask(i,j) == 3)
then
2321 cs%thickness_bdry_val(i,j) = input_thick
2324 if ((hmask(i,j) == 0) .or. (hmask(i,j) == 1) .or. (hmask(i,j) == 2))
then
2325 if ((i <= iec).and.(i >= isc))
then
2326 if (cs%u_face_mask(i-1,j) == 3)
then
2327 cs%u_bdry_val(i-1,j-1) = (1 - ((g%geoLatBu(i-1,j-1) - 0.5*g%len_lat)*2./g%len_lat)**2) * &
2328 1.5 * input_flux / input_thick
2329 cs%u_bdry_val(i-1,j) = (1 - ((g%geoLatBu(i-1,j) - 0.5*g%len_lat)*2./g%len_lat)**2) * &
2330 1.5 * input_flux / input_thick
2335 if (.not.(new_sim))
then
2336 if (.not. g%symmetric)
then
2337 if (((i+i_off) == (g%domain%nihalo+1)).and.(cs%u_face_mask(i-1,j) == 3))
then
2338 cs%u_shelf(i-1,j-1) = cs%u_bdry_val(i-1,j-1)
2339 cs%u_shelf(i-1,j) = cs%u_bdry_val(i-1,j)
2341 if (((j+j_off) == (g%domain%njhalo+1)).and.(cs%v_face_mask(i,j-1) == 3))
then
2342 cs%u_shelf(i-1,j-1) = cs%u_bdry_val(i-1,j-1)
2343 cs%u_shelf(i,j-1) = cs%u_bdry_val(i,j-1)
2353 subroutine cg_action(uret, vret, u, v, Phi, Phisub, umask, vmask, hmask, H_node, &
2354 nu, float_cond, bathyT, beta, dxdyh, G, is, ie, js, je, dens_ratio)
2357 real,
dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), &
2358 intent(inout) :: uret
2359 real,
dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), &
2360 intent(inout) :: vret
2361 real,
dimension(SZDI_(G),SZDJ_(G),8,4), &
2364 real,
dimension(:,:,:,:,:,:), &
2365 intent(in) :: Phisub
2367 real,
dimension(SZDIB_(G),SZDJB_(G)), &
2369 real,
dimension(SZDIB_(G),SZDJB_(G)), &
2371 real,
dimension(SZDIB_(G),SZDJB_(G)), &
2374 real,
dimension(SZDIB_(G),SZDJB_(G)), &
2377 real,
dimension(SZDIB_(G),SZDJB_(G)), &
2378 intent(in) :: H_node
2380 real,
dimension(SZDI_(G),SZDJ_(G)), &
2383 real,
dimension(SZDIB_(G),SZDJB_(G)), &
2387 real,
dimension(SZDI_(G),SZDJ_(G)), &
2388 intent(in) :: float_cond
2390 real,
dimension(SZDI_(G),SZDJ_(G)), &
2391 intent(in) :: bathyT
2392 real,
dimension(SZDIB_(G),SZDJB_(G)), &
2397 real,
dimension(SZDI_(G),SZDJ_(G)), &
2399 real,
intent(in) :: dens_ratio
2401 integer,
intent(in) :: is
2402 integer,
intent(in) :: ie
2403 integer,
intent(in) :: js
2404 integer,
intent(in) :: je
2425 real :: ux, vx, uy, vy, uq, vq, area, basel
2426 integer :: iq, jq, iphi, jphi, i, j, ilq, jlq
2427 real,
dimension(2) :: xquad
2428 real,
dimension(2,2) :: Ucell,Vcell,Hcell,Usubcontr,Vsubcontr,Ucontr
2430 xquad(1) = .5 * (1-sqrt(1./3)) ; xquad(2) = .5 * (1+sqrt(1./3))
2433 do i=is,ie ;
if (hmask(i,j) == 1)
then
2452 do iq=1,2 ;
do jq=1,2
2467 uq = u(i-1,j-1) * xquad(3-iq) * xquad(3-jq) + &
2468 u(i,j-1) * xquad(iq) * xquad(3-jq) + &
2469 u(i-1,j) * xquad(3-iq) * xquad(jq) + &
2470 u(i,j) * xquad(iq) * xquad(jq)
2472 vq = v(i-1,j-1) * xquad(3-iq) * xquad(3-jq) + &
2473 v(i,j-1) * xquad(iq) * xquad(3-jq) + &
2474 v(i-1,j) * xquad(3-iq) * xquad(jq) + &
2475 v(i,j) * xquad(iq) * xquad(jq)
2477 ux = u(i-1,j-1) * phi(i,j,1,2*(jq-1)+iq) + &
2478 u(i,j-1) * phi(i,j,3,2*(jq-1)+iq) + &
2479 u(i-1,j) * phi(i,j,5,2*(jq-1)+iq) + &
2480 u(i,j) * phi(i,j,7,2*(jq-1)+iq)
2482 vx = v(i-1,j-1) * phi(i,j,1,2*(jq-1)+iq) + &
2483 v(i,j-1) * phi(i,j,3,2*(jq-1)+iq) + &
2484 v(i-1,j) * phi(i,j,5,2*(jq-1)+iq) + &
2485 v(i,j) * phi(i,j,7,2*(jq-1)+iq)
2487 uy = u(i-1,j-1) * phi(i,j,2,2*(jq-1)+iq) + &
2488 u(i,j-1) * phi(i,j,4,2*(jq-1)+iq) + &
2489 u(i-1,j) * phi(i,j,6,2*(jq-1)+iq) + &
2490 u(i,j) * phi(i,j,8,2*(jq-1)+iq)
2492 vy = v(i-1,j-1) * phi(i,j,2,2*(jq-1)+iq) + &
2493 v(i,j-1) * phi(i,j,4,2*(jq-1)+iq) + &
2494 v(i-1,j) * phi(i,j,6,2*(jq-1)+iq) + &
2495 v(i,j) * phi(i,j,8,2*(jq-1)+iq)
2497 do iphi=1,2 ;
do jphi=1,2
2498 if (umask(i-2+iphi,j-2+jphi) == 1)
then
2500 uret(i-2+iphi,j-2+jphi) = uret(i-2+iphi,j-2+jphi) + &
2501 .25 * area * nu(i,j) * ((4*ux+2*vy) * phi(i,j,2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + &
2502 (uy+vx) * phi(i,j,2*(2*(jphi-1)+iphi),2*(jq-1)+iq))
2504 if (vmask(i-2+iphi,j-2+jphi) == 1)
then
2506 vret(i-2+iphi,j-2+jphi) = vret(i-2+iphi,j-2+jphi) + &
2507 .25 * area * nu(i,j) * ((uy+vx) * phi(i,j,2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + &
2508 (4*vy+2*ux) * phi(i,j,2*(2*(jphi-1)+iphi),2*(jq-1)+iq))
2511 if (iq == iphi)
then
2517 if (jq == jphi)
then
2523 if (float_cond(i,j) == 0)
then
2525 if (umask(i-2+iphi,j-2+jphi) == 1)
then
2527 uret(i-2+iphi,j-2+jphi) = uret(i-2+iphi,j-2+jphi) + &
2528 .25 * beta(i,j) * area * uq * xquad(ilq) * xquad(jlq)
2532 if (vmask(i-2+iphi,j-2+jphi) == 1)
then
2534 vret(i-2+iphi,j-2+jphi) = vret(i-2+iphi,j-2+jphi) + &
2535 .25 * beta(i,j) * area * vq * xquad(ilq) * xquad(jlq)
2540 ucontr(iphi,jphi) = ucontr(iphi,jphi) + .25 * area * uq * xquad(ilq) * xquad(jlq) * beta(i,j)
2544 if (float_cond(i,j) == 1)
then
2545 usubcontr = 0.0 ; vsubcontr = 0.0 ; basel = bathyt(i,j)
2546 ucell(:,:) = u(i-1:i,j-1:j) ; vcell(:,:) = v(i-1:i,j-1:j) ; hcell(:,:) = h_node(i-1:i,j-1:j)
2548 dens_ratio, usubcontr, vsubcontr)
2549 do iphi=1,2 ;
do jphi=1,2
2550 if (umask(i-2+iphi,j-2+jphi) == 1)
then
2551 uret(i-2+iphi,j-2+jphi) = uret(i-2+iphi,j-2+jphi) + usubcontr(iphi,jphi) * beta(i,j)
2553 if (vmask(i-2+iphi,j-2+jphi) == 1)
then
2554 vret(i-2+iphi,j-2+jphi) = vret(i-2+iphi,j-2+jphi) + vsubcontr(iphi,jphi) * beta(i,j)
2565 real,
dimension(:,:,:,:,:,:), &
2566 intent(in) :: Phisub
2568 real,
dimension(2,2),
intent(in) :: H
2569 real,
dimension(2,2),
intent(in) :: U
2570 real,
dimension(2,2),
intent(in) :: V
2571 real,
intent(in) :: DXDYH
2572 real,
intent(in) :: bathyT
2573 real,
intent(in) :: dens_ratio
2575 real,
dimension(2,2),
intent(inout) :: Ucontr
2577 real,
dimension(2,2),
intent(inout) :: Vcontr
2580 integer :: nsub, i, j, k, l, qx, qy, m, n
2581 real :: subarea, hloc, uq, vq
2583 nsub =
size(phisub,1)
2584 subarea = dxdyh / (nsub**2)
2593 hloc = phisub(i,j,1,1,qx,qy)*h(1,1) + phisub(i,j,1,2,qx,qy)*h(1,2) + &
2594 phisub(i,j,2,1,qx,qy)*h(2,1) + phisub(i,j,2,2,qx,qy)*h(2,2)
2596 if (dens_ratio * hloc - bathyt > 0)
then
2603 uq = uq + phisub(i,j,k,l,qx,qy) * u(k,l) ; vq = vq + phisub(i,j,k,l,qx,qy) * v(k,l)
2607 ucontr(m,n) = ucontr(m,n) + subarea * 0.25 * phisub(i,j,m,n,qx,qy) * uq
2608 vcontr(m,n) = vcontr(m,n) + subarea * 0.25 * phisub(i,j,m,n,qx,qy) * vq
2622 subroutine matrix_diagonal(CS, G, float_cond, H_node, nu, beta, hmask, dens_ratio, &
2623 Phisub, u_diagonal, v_diagonal)
2627 real,
dimension(SZDI_(G),SZDJ_(G)), &
2628 intent(in) :: float_cond
2630 real,
dimension(SZDIB_(G),SZDJB_(G)), &
2631 intent(in) :: H_node
2633 real,
dimension(SZDIB_(G),SZDJB_(G)), &
2637 real,
dimension(SZDIB_(G),SZDJB_(G)), &
2641 real,
dimension(SZDI_(G),SZDJ_(G)), &
2644 real,
intent(in) :: dens_ratio
2646 real,
dimension(:,:,:,:,:,:),
intent(in) :: Phisub
2648 real,
dimension(SZDIB_(G),SZDJB_(G)), &
2649 intent(inout) :: u_diagonal
2651 real,
dimension(SZDIB_(G),SZDJB_(G)), &
2652 intent(inout) :: v_diagonal
2658 integer :: i, j, is, js, cnt, isc, jsc, iec, jec, iphi, jphi, iq, jq, ilq, jlq
2659 real :: A, n, ux, uy, vx, vy, eps_min, domain_width, dxh, dyh, dxdyh, area, uq, vq, basel
2660 real,
dimension(8,4) :: Phi
2661 real,
dimension(4) :: X, Y
2662 real,
dimension(2) :: xquad
2663 real,
dimension(2,2) :: Hcell,Usubcontr,Vsubcontr
2666 isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
2668 xquad(1) = .5 * (1-sqrt(1./3)) ; xquad(2) = .5 * (1+sqrt(1./3))
2677 do j=jsc-1,jec+1 ;
do i=isc-1,iec+1 ;
if (hmask(i,j) == 1)
then
2679 dxh = g%US%L_to_m*g%dxT(i,j)
2680 dyh = g%US%L_to_m*g%dyT(i,j)
2681 dxdyh = g%US%L_to_m**2*g%areaT(i,j)
2683 x(1:2) = g%geoLonBu(i-1:i,j-1)*1000
2684 x(3:4) = g%geoLonBu(i-1:i,j) *1000
2685 y(1:2) = g%geoLatBu(i-1:i,j-1) *1000
2686 y(3:4) = g%geoLatBu(i-1:i,j)*1000
2697 do iq=1,2 ;
do jq=1,2
2699 do iphi=1,2 ;
do jphi=1,2
2701 if (iq == iphi)
then
2707 if (jq == jphi)
then
2713 if (cs%umask(i-2+iphi,j-2+jphi) == 1)
then
2715 ux = phi(2*(2*(jphi-1)+iphi)-1, 2*(jq-1)+iq)
2716 uy = phi(2*(2*(jphi-1)+iphi), 2*(jq-1)+iq)
2720 u_diagonal(i-2+iphi,j-2+jphi) = u_diagonal(i-2+iphi,j-2+jphi) + &
2721 .25 * dxdyh * nu(i,j) * ((4*ux+2*vy) * phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + &
2722 (uy+vy) * phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq))
2724 uq = xquad(ilq) * xquad(jlq)
2726 if (float_cond(i,j) == 0)
then
2727 u_diagonal(i-2+iphi,j-2+jphi) = u_diagonal(i-2+iphi,j-2+jphi) + &
2728 .25 * beta(i,j) * dxdyh * uq * xquad(ilq) * xquad(jlq)
2733 if (cs%vmask(i-2+iphi,j-2+jphi) == 1)
then
2735 vx = phi(2*(2*(jphi-1)+iphi)-1, 2*(jq-1)+iq)
2736 vy = phi(2*(2*(jphi-1)+iphi), 2*(jq-1)+iq)
2740 v_diagonal(i-2+iphi,j-2+jphi) = v_diagonal(i-2+iphi,j-2+jphi) + &
2741 .25 * dxdyh * nu(i,j) * ((uy+vx) * phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + &
2742 (4*vy+2*ux) * phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq))
2744 vq = xquad(ilq) * xquad(jlq)
2746 if (float_cond(i,j) == 0)
then
2747 v_diagonal(i-2+iphi,j-2+jphi) = v_diagonal(i-2+iphi,j-2+jphi) + &
2748 .25 * beta(i,j) * dxdyh * vq * xquad(ilq) * xquad(jlq)
2754 if (float_cond(i,j) == 1)
then
2755 usubcontr = 0.0 ; vsubcontr = 0.0 ; basel = g%bathyT(i,j)
2756 hcell(:,:) = h_node(i-1:i,j-1:j)
2758 do iphi=1,2 ;
do jphi=1,2
2759 if (cs%umask(i-2+iphi,j-2+jphi) == 1)
then
2760 u_diagonal(i-2+iphi,j-2+jphi) = u_diagonal(i-2+iphi,j-2+jphi) + usubcontr(iphi,jphi) * beta(i,j)
2761 v_diagonal(i-2+iphi,j-2+jphi) = v_diagonal(i-2+iphi,j-2+jphi) + vsubcontr(iphi,jphi) * beta(i,j)
2765 endif ;
enddo ;
enddo
2770 real,
dimension(:,:,:,:,:,:), &
2771 intent(in) :: Phisub
2773 real,
dimension(2,2),
intent(in) :: H_node
2775 real,
intent(in) :: DXDYH
2776 real,
intent(in) :: bathyT
2777 real,
intent(in) :: dens_ratio
2779 real,
dimension(2,2),
intent(inout) :: Ucontr
2781 real,
dimension(2,2),
intent(inout) :: Vcontr
2786 integer :: nsub, i, j, k, l, qx, qy, m, n
2787 real :: subarea, hloc
2789 nsub =
size(phisub,1)
2790 subarea = dxdyh / (nsub**2)
2792 do m=1,2 ;
do n=1,2 ;
do j=1,nsub ;
do i=1,nsub ;
do qx=1,2 ;
do qy = 1,2
2794 hloc = phisub(i,j,1,1,qx,qy)*h_node(1,1) + phisub(i,j,1,2,qx,qy)*h_node(1,2) + &
2795 phisub(i,j,2,1,qx,qy)*h_node(2,1) + phisub(i,j,2,2,qx,qy)*h_node(2,2)
2797 if (dens_ratio * hloc - bathyt > 0)
then
2798 ucontr(m,n) = ucontr(m,n) + subarea * 0.25 * phisub(i,j,m,n,qx,qy)**2
2799 vcontr(m,n) = vcontr(m,n) + subarea * 0.25 * phisub(i,j,m,n,qx,qy)**2
2802 enddo ;
enddo ;
enddo ;
enddo ;
enddo ;
enddo
2808 dens_ratio, u_bdry_contr, v_bdry_contr)
2814 type(time_type),
intent(in) :: Time
2815 real,
dimension(:,:,:,:,:,:), &
2816 intent(in) :: Phisub
2818 real,
dimension(SZDIB_(G),SZDJB_(G)), &
2819 intent(in) :: H_node
2821 real,
dimension(SZDIB_(G),SZDJB_(G)), &
2825 real,
dimension(SZDIB_(G),SZDJB_(G)), &
2829 real,
dimension(SZDI_(G),SZDJ_(G)), &
2830 intent(in) :: float_cond
2832 real,
intent(in) :: dens_ratio
2834 real,
dimension(SZDIB_(G),SZDJB_(G)), &
2835 intent(inout) :: u_bdry_contr
2837 real,
dimension(SZDIB_(G),SZDJB_(G)), &
2838 intent(inout) :: v_bdry_contr
2844 real,
dimension(8,4) :: Phi
2845 real,
dimension(4) :: X, Y
2846 real,
dimension(2) :: xquad
2847 integer :: i, j, isc, jsc, iec, jec, iq, jq, iphi, jphi, ilq, jlq
2848 real :: A, n, ux, uy, vx, vy, eps_min, domain_width, dxh, dyh, dxdyh, uq, vq, area, basel
2849 real,
dimension(2,2) :: Ucell,Vcell,Hcell,Usubcontr,Vsubcontr
2852 isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
2854 xquad(1) = .5 * (1-sqrt(1./3)) ; xquad(2) = .5 * (1+sqrt(1./3))
2863 do j=jsc-1,jec+1 ;
do i=isc-1,iec+1 ;
if (iss%hmask(i,j) == 1)
then
2868 if ((cs%umask(i-1,j-1) == 3) .OR. (cs%umask(i,j-1) == 3) .OR. &
2869 (cs%umask(i-1,j) == 3) .OR. (cs%umask(i,j) == 3))
then
2871 dxh = g%US%L_to_m*g%dxT(i,j)
2872 dyh = g%US%L_to_m*g%dyT(i,j)
2873 dxdyh = g%US%L_to_m**2*g%areaT(i,j)
2875 x(1:2) = g%geoLonBu(i-1:i,j-1)*1000
2876 x(3:4) = g%geoLonBu(i-1:i,j)*1000
2877 y(1:2) = g%geoLatBu(i-1:i,j-1)*1000
2878 y(3:4) = g%geoLatBu(i-1:i,j)*1000
2889 do iq=1,2 ;
do jq=1,2
2891 uq = cs%u_bdry_val(i-1,j-1) * xquad(3-iq) * xquad(3-jq) + &
2892 cs%u_bdry_val(i,j-1) * xquad(iq) * xquad(3-jq) + &
2893 cs%u_bdry_val(i-1,j) * xquad(3-iq) * xquad(jq) + &
2894 cs%u_bdry_val(i,j) * xquad(iq) * xquad(jq)
2896 vq = cs%v_bdry_val(i-1,j-1) * xquad(3-iq) * xquad(3-jq) + &
2897 cs%v_bdry_val(i,j-1) * xquad(iq) * xquad(3-jq) + &
2898 cs%v_bdry_val(i-1,j) * xquad(3-iq) * xquad(jq) + &
2899 cs%v_bdry_val(i,j) * xquad(iq) * xquad(jq)
2901 ux = cs%u_bdry_val(i-1,j-1) * phi(1,2*(jq-1)+iq) + &
2902 cs%u_bdry_val(i,j-1) * phi(3,2*(jq-1)+iq) + &
2903 cs%u_bdry_val(i-1,j) * phi(5,2*(jq-1)+iq) + &
2904 cs%u_bdry_val(i,j) * phi(7,2*(jq-1)+iq)
2906 vx = cs%v_bdry_val(i-1,j-1) * phi(1,2*(jq-1)+iq) + &
2907 cs%v_bdry_val(i,j-1) * phi(3,2*(jq-1)+iq) + &
2908 cs%v_bdry_val(i-1,j) * phi(5,2*(jq-1)+iq) + &
2909 cs%v_bdry_val(i,j) * phi(7,2*(jq-1)+iq)
2911 uy = cs%u_bdry_val(i-1,j-1) * phi(2,2*(jq-1)+iq) + &
2912 cs%u_bdry_val(i,j-1) * phi(4,2*(jq-1)+iq) + &
2913 cs%u_bdry_val(i-1,j) * phi(6,2*(jq-1)+iq) + &
2914 cs%u_bdry_val(i,j) * phi(8,2*(jq-1)+iq)
2916 vy = cs%v_bdry_val(i-1,j-1) * phi(2,2*(jq-1)+iq) + &
2917 cs%v_bdry_val(i,j-1) * phi(4,2*(jq-1)+iq) + &
2918 cs%v_bdry_val(i-1,j) * phi(6,2*(jq-1)+iq) + &
2919 cs%v_bdry_val(i,j) * phi(8,2*(jq-1)+iq)
2921 do iphi=1,2 ;
do jphi=1,2
2923 if (iq == iphi)
then
2929 if (jq == jphi)
then
2935 if (cs%umask(i-2+iphi,j-2+jphi) == 1)
then
2938 u_bdry_contr(i-2+iphi,j-2+jphi) = u_bdry_contr(i-2+iphi,j-2+jphi) + &
2939 .25 * dxdyh * nu(i,j) * ( (4*ux+2*vy) * phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + &
2940 (uy+vx) * phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq) )
2942 if (float_cond(i,j) == 0)
then
2943 u_bdry_contr(i-2+iphi,j-2+jphi) = u_bdry_contr(i-2+iphi,j-2+jphi) + &
2944 .25 * beta(i,j) * dxdyh * uq * xquad(ilq) * xquad(jlq)
2949 if (cs%vmask(i-2+iphi,j-2+jphi) == 1)
then
2952 v_bdry_contr(i-2+iphi,j-2+jphi) = v_bdry_contr(i-2+iphi,j-2+jphi) + &
2953 .25 * dxdyh * nu(i,j) * ( (uy+vx) * phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + &
2954 (4*vy+2*ux) * phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq))
2956 if (float_cond(i,j) == 0)
then
2957 v_bdry_contr(i-2+iphi,j-2+jphi) = v_bdry_contr(i-2+iphi,j-2+jphi) + &
2958 .25 * beta(i,j) * dxdyh * vq * xquad(ilq) * xquad(jlq)
2965 if (float_cond(i,j) == 1)
then
2966 usubcontr = 0.0 ; vsubcontr = 0.0 ; basel = g%bathyT(i,j)
2967 ucell(:,:) = cs%u_bdry_val(i-1:i,j-1:j) ; vcell(:,:) = cs%v_bdry_val(i-1:i,j-1:j)
2968 hcell(:,:) = h_node(i-1:i,j-1:j)
2970 dens_ratio, usubcontr, vsubcontr)
2971 do iphi=1,2 ;
do jphi = 1,2
2972 if (cs%umask(i-2+iphi,j-2+jphi) == 1)
then
2973 u_bdry_contr(i-2+iphi,j-2+jphi) = u_bdry_contr(i-2+iphi,j-2+jphi) + &
2974 usubcontr(iphi,jphi) * beta(i,j)
2976 if (cs%vmask(i-2+iphi,j-2+jphi) == 1)
then
2977 v_bdry_contr(i-2+iphi,j-2+jphi) = v_bdry_contr(i-2+iphi,j-2+jphi) + &
2978 vsubcontr(iphi,jphi) * beta(i,j)
2983 endif ;
enddo ;
enddo
2995 real,
dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), &
2997 real,
dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), &
3007 integer :: i, j, iscq, iecq, jscq, jecq, isd, jsd, ied, jed, iegq, jegq
3008 integer :: giec, gjec, gisc, gjsc, cnt, isc, jsc, iec, jec, is, js
3009 real :: A, n, ux, uy, vx, vy, eps_min, umid, vmid, unorm, C_basal_friction, n_basal_friction, dxh, dyh, dxdyh
3011 isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
3012 iscq = g%iscB ; iecq = g%iecB ; jscq = g%jscB ; jecq = g%jecB
3013 isd = g%isd ; jsd = g%jsd ; ied = g%ied ; jed = g%jed
3014 iegq = g%iegB ; jegq = g%jegB
3015 gisc = g%domain%nihalo+1 ; gjsc = g%domain%njhalo+1
3016 giec = g%domain%niglobal+gisc ; gjec = g%domain%njglobal+gjsc
3017 is = iscq - 1; js = jscq - 1
3019 a = cs%A_glen_isothermal ; n = cs%n_glen; eps_min = cs%eps_glen_min
3020 c_basal_friction = cs%C_basal_friction ; n_basal_friction = cs%n_basal_friction
3025 dxh = us%L_to_m*g%dxT(i,j)
3026 dyh = us%L_to_m*g%dyT(i,j)
3027 dxdyh = us%L_to_m**2*g%areaT(i,j)
3029 if (iss%hmask(i,j) == 1)
then
3030 ux = (u(i,j) + u(i,j-1) - u(i-1,j) - u(i-1,j-1)) / (2*dxh)
3031 vx = (v(i,j) + v(i,j-1) - v(i-1,j) - v(i-1,j-1)) / (2*dxh)
3032 uy = (u(i,j) - u(i,j-1) + u(i-1,j) - u(i-1,j-1)) / (2*dyh)
3033 vy = (v(i,j) - v(i,j-1) + v(i-1,j) - v(i-1,j-1)) / (2*dyh)
3035 cs%ice_visc(i,j) = .5 * a**(-1/n) * &
3036 (ux**2+vy**2+ux*vy+0.25*(uy+vx)**2+eps_min**2) ** ((1-n)/(2*n)) * &
3037 us%Z_to_m*iss%h_shelf(i,j)
3039 umid = (u(i,j) + u(i,j-1) + u(i-1,j) + u(i-1,j-1))/4
3040 vmid = (v(i,j) + v(i,j-1) + v(i-1,j) + v(i-1,j-1))/4
3041 unorm = sqrt(umid**2+vmid**2+(eps_min*dxh)**2)
3042 cs%taub_beta_eff(i,j) = c_basal_friction * unorm ** (n_basal_friction-1)
3053 real,
dimension(SZDI_(G),SZDJ_(G)), &
3054 intent(in) :: ocean_mass
3055 logical,
intent(in) :: find_avg
3058 integer :: isc, iec, jsc, jec, i, j
3062 i_rho_ocean = 1.0 / (us%Z_to_m*cs%density_ocean_avg)
3064 isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
3066 do j=jsc,jec ;
do i=isc,iec
3067 cs%OD_rt(i,j) = cs%OD_rt(i,j) + ocean_mass(i,j)*i_rho_ocean
3068 if (ocean_mass(i,j)*i_rho_ocean > cs%thresh_float_col_depth)
then
3069 cs%float_frac_rt(i,j) = cs%float_frac_rt(i,j) + 1.0
3072 cs%OD_rt_counter = cs%OD_rt_counter + 1
3075 i_counter = 1.0 / real(cs%OD_rt_counter)
3076 do j=jsc,jec ;
do i=isc,iec
3077 cs%float_frac(i,j) = 1.0 - (cs%float_frac_rt(i,j) * i_counter)
3078 cs%OD_av(i,j) = cs%OD_rt(i,j) * i_counter
3080 cs%OD_rt(i,j) = 0.0 ; cs%float_frac_rt(i,j) = 0.0
3083 call pass_var(cs%float_frac, g%domain)
3092 real,
dimension(SZDI_(G),SZDJ_(G)), &
3093 intent(in) :: h_shelf
3095 integer :: i, j, iters, isd, ied, jsd, jed
3096 real :: rhoi_rhow, OD
3098 rhoi_rhow = cs%density_ice / cs%density_ocean_avg
3099 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
3103 od = g%bathyT(i,j) - rhoi_rhow * h_shelf(i,j)
3107 cs%float_frac(i,j) = 0.
3110 cs%float_frac(i,j) = 1.
3121 real,
dimension(4),
intent(in) :: X
3122 real,
dimension(4),
intent(in) :: Y
3123 real,
dimension(8,4),
intent(inout) :: Phi
3125 real,
intent(out) :: area
3144 real,
dimension(4) :: xquad, yquad
3145 integer :: node, qpoint, xnode, xq, ynode, yq
3146 real :: a,b,c,d,e,f,xexp,yexp
3148 xquad(1:3:2) = .5 * (1-sqrt(1./3)) ; yquad(1:2) = .5 * (1-sqrt(1./3))
3149 xquad(2:4:2) = .5 * (1+sqrt(1./3)) ; yquad(3:4) = .5 * (1+sqrt(1./3))
3153 a = -x(1)*(1-yquad(qpoint)) + x(2)*(1-yquad(qpoint)) - x(3)*yquad(qpoint) + x(4)*yquad(qpoint)
3154 b = -y(1)*(1-yquad(qpoint)) + y(2)*(1-yquad(qpoint)) - y(3)*yquad(qpoint) + y(4)*yquad(qpoint)
3155 c = -x(1)*(1-xquad(qpoint)) - x(2)*(xquad(qpoint)) + x(3)*(1-xquad(qpoint)) + x(4)*(xquad(qpoint))
3156 d = -y(1)*(1-xquad(qpoint)) - y(2)*(xquad(qpoint)) + y(3)*(1-xquad(qpoint)) + y(4)*(xquad(qpoint))
3160 xnode = 2-mod(node,2) ; ynode = ceiling(real(node)/2)
3162 if (ynode == 1)
then
3163 yexp = 1-yquad(qpoint)
3165 yexp = yquad(qpoint)
3168 if (1 == xnode)
then
3169 xexp = 1-xquad(qpoint)
3171 xexp = xquad(qpoint)
3174 phi(2*node-1,qpoint) = ( d * (2 * xnode - 3) * yexp - b * (2 * ynode - 3) * xexp) / (a*d-b*c)
3175 phi(2*node,qpoint) = ( -c * (2 * xnode - 3) * yexp + a * (2 * ynode - 3) * xexp) / (a*d-b*c)
3186 real,
dimension(nsub,nsub,2,2,2,2), &
3187 intent(inout) :: Phisub
3189 integer,
intent(in) :: nsub
3214 integer :: i, j, k, l, qx, qy, indx, indy
3215 real,
dimension(2) :: xquad
3216 real :: x0, y0, x, y, val, fracx
3218 xquad(1) = .5 * (1-sqrt(1./3)) ; xquad(2) = .5 * (1+sqrt(1./3))
3219 fracx = 1.0/real(nsub)
3223 x0 = (i-1) * fracx ; y0 = (j-1) * fracx
3226 x = x0 + fracx*xquad(qx)
3227 y = y0 + fracx*xquad(qy)
3241 phisub(i,j,k,l,qx,qy) = val
3255 real,
dimension(SZDI_(G),SZDJ_(G)), &
3258 real,
dimension(SZDIB_(G),SZDJB_(G)), &
3259 intent(out) :: umask
3261 real,
dimension(SZDIB_(G),SZDJB_(G)), &
3262 intent(out) :: vmask
3264 real,
dimension(SZDIB_(G),SZDJ_(G)), &
3265 intent(out) :: u_face_mask
3266 real,
dimension(SZDI_(G),SZDJB_(G)), &
3267 intent(out) :: v_face_mask
3273 integer :: i, j, k, iscq, iecq, jscq, jecq, isd, jsd, is, js, iegq, jegq
3274 integer :: giec, gjec, gisc, gjsc, isc, jsc, iec, jec
3275 integer :: i_off, j_off
3277 isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
3278 iscq = g%iscB ; iecq = g%iecB ; jscq = g%jscB ; jecq = g%jecB
3279 i_off = g%idg_offset ; j_off = g%jdg_offset
3280 isd = g%isd ; jsd = g%jsd
3281 iegq = g%iegB ; jegq = g%jegB
3282 gisc = g%Domain%nihalo ; gjsc = g%Domain%njhalo
3283 giec = g%Domain%niglobal+gisc ; gjec = g%Domain%njglobal+gjsc
3285 umask(:,:) = 0 ; vmask(:,:) = 0
3286 u_face_mask(:,:) = 0 ; v_face_mask(:,:) = 0
3288 if (g%symmetric)
then
3291 is = isd+1 ; js = jsd+1
3297 if (hmask(i,j) == 1)
then
3299 umask(i-1:i,j-1:j) = 1.
3300 vmask(i-1:i,j-1:j) = 1.
3304 select case (int(cs%u_face_mask_bdry(i-1+k,j)))
3306 umask(i-1+k,j-1:j)=3.
3307 vmask(i-1+k,j-1:j)=0.
3308 u_face_mask(i-1+k,j)=3.
3310 u_face_mask(i-1+k,j)=2.
3312 umask(i-1+k,j-1:j)=0.
3313 vmask(i-1+k,j-1:j)=0.
3314 u_face_mask(i-1+k,j)=4.
3316 umask(i-1+k,j-1:j)=0.
3317 vmask(i-1+k,j-1:j)=0.
3318 u_face_mask(i-1+k,j)=0.
3320 umask(i-1+k,j-1:j)=0.
3327 select case (int(cs%v_face_mask_bdry(i,j-1+k)))
3329 vmask(i-1:i,j-1+k)=3.
3330 umask(i-1:i,j-1+k)=0.
3331 v_face_mask(i,j-1+k)=3.
3333 v_face_mask(i,j-1+k)=2.
3335 umask(i-1:i,j-1+k)=0.
3336 vmask(i-1:i,j-1+k)=0.
3337 v_face_mask(i,j-1+k)=4.
3339 umask(i-1:i,j-1+k)=0.
3340 vmask(i-1:i,j-1+k)=0.
3341 u_face_mask(i,j-1+k)=0.
3343 vmask(i-1:i,j-1+k)=0.
3365 if ((hmask(i+1,j) == 0) &
3366 .OR. (hmask(i+1,j) == 2))
then
3368 u_face_mask(i,j) = 2.
3373 if ((hmask(i-1,j) == 0) .OR. (hmask(i-1,j) == 2))
then
3375 u_face_mask(i-1,j) = 2.
3380 if ((hmask(i,j-1) == 0) .OR. (hmask(i,j-1) == 2))
then
3382 v_face_mask(i,j-1) = 2.
3387 if ((hmask(i,j+1) == 0) .OR. (hmask(i,j+1) == 2))
then
3389 v_face_mask(i,j) = 2.
3411 real,
dimension(SZDI_(G),SZDJ_(G)), &
3412 intent(in) :: h_shelf
3413 real,
dimension(SZDI_(G),SZDJ_(G)), &
3416 real,
dimension(SZDIB_(G),SZDJB_(G)), &
3417 intent(inout) :: H_node
3420 integer :: i, j, isc, iec, jsc, jec, num_h, k, l
3423 isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
3436 if (hmask(i+k,j+l) == 1.0)
then
3437 summ = summ + h_shelf(i+k,j+l)
3443 h_node(i,j) = summ / num_h
3448 call pass_var(h_node, g%domain, position=corner)
3456 if (.not.
associated(cs))
return
3458 deallocate(cs%u_shelf, cs%v_shelf)
3459 deallocate(cs%t_shelf, cs%tmask)
3460 deallocate(cs%u_bdry_val, cs%v_bdry_val, cs%t_bdry_val)
3461 deallocate(cs%u_face_mask, cs%v_face_mask)
3462 deallocate(cs%umask, cs%vmask)
3464 deallocate(cs%ice_visc, cs%taub_beta_eff)
3465 deallocate(cs%OD_rt, cs%OD_av)
3466 deallocate(cs%float_frac, cs%float_frac_rt)
3474 subroutine ice_shelf_temp(CS, ISS, G, US, time_step, melt_rate, Time)
3480 real,
intent(in) :: time_step
3481 real,
dimension(SZDI_(G),SZDJ_(G)), &
3482 intent(in) :: melt_rate
3483 type(time_type),
intent(in) :: Time
3514 real,
dimension(SZDI_(G),SZDJ_(G)) :: th_after_uflux, th_after_vflux, TH
3515 real,
dimension(SZDI_(G),SZDJ_(G),4) :: flux_enter
3516 integer :: isd, ied, jsd, jed, i, j, isc, iec, jsc, jec
3517 real :: rho, spy, t_bd, Tsurf, adot
3519 rho = cs%density_ice
3522 adot = 0.1*us%m_to_Z/spy
3525 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
3526 isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
3527 flux_enter(:,:,:) = 0.0
3529 th_after_uflux(:,:) = 0.0
3530 th_after_vflux(:,:) = 0.0
3534 t_bd = cs%t_bdry_val(i,j)
3536 if ((iss%hmask(i,j) == 3) .or. (iss%hmask(i,j) == -2))
then
3537 cs%t_shelf(i,j) = cs%t_bdry_val(i,j)
3544 th(i,j) = cs%t_shelf(i,j)*iss%h_shelf(i,j)
3562 if (iss%h_shelf(i,j) > 0.0)
then
3563 cs%t_shelf(i,j) = th_after_vflux(i,j)/(iss%h_shelf(i,j))
3565 cs%t_shelf(i,j) = -10.0
3572 t_bd = cs%t_bdry_val(i,j)
3574 if ((iss%hmask(i,j) == 3) .or. (iss%hmask(i,j) == -2))
then
3575 cs%t_shelf(i,j) = t_bd
3583 if ((iss%hmask(i,j) == 1) .or. (iss%hmask(i,j) == 2))
then
3584 if (iss%h_shelf(i,j) > 0.0)
then
3587 cs%t_shelf(i,j) = cs%t_shelf(i,j) + &
3588 time_step*(adot*tsurf - (3.0*us%m_to_Z/spy)*iss%tfreeze(i,j)) / iss%h_shelf(i,j)
3595 cs%t_shelf(i,j) = -10.0
3602 call pass_var(cs%t_shelf, g%domain)
3606 call hchksum(cs%t_shelf,
"temp after front", g%HI, haloshift=3)
3615 real,
intent(in) :: time_step
3616 real,
dimension(SZDI_(G),SZDJ_(G)), &
3619 real,
dimension(SZDI_(G),SZDJ_(G)), &
3621 real,
dimension(SZDI_(G),SZDJ_(G)), &
3622 intent(inout) :: h_after_uflux
3624 real,
dimension(SZDI_(G),SZDJ_(G),4), &
3625 intent(inout) :: flux_enter
3646 integer :: i, j, is, ie, js, je, isd, ied, jsd, jed, gjed, gied
3647 integer :: i_off, j_off
3648 logical :: at_east_bdry, at_west_bdry, one_off_west_bdry, one_off_east_bdry
3649 real,
dimension(-2:2) :: stencil
3651 flux_diff_cell, phi, dxh, dyh, dxdyh
3653 character (len=1) :: debug_str
3656 is = g%isc-2 ; ie = g%iec+2 ; js = g%jsc ; je = g%jec ; isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
3657 i_off = g%idg_offset ; j_off = g%jdg_offset
3660 if (((j+j_off) <= g%domain%njglobal+g%domain%njhalo) .AND. &
3661 ((j+j_off) >= g%domain%njhalo+1))
then
3667 if (((i+i_off) <= g%domain%niglobal+g%domain%nihalo) .AND. &
3668 ((i+i_off) >= g%domain%nihalo+1))
then
3670 if (i+i_off == g%domain%nihalo+1)
then
3673 at_west_bdry=.false.
3676 if (i+i_off == g%domain%niglobal+g%domain%nihalo)
then
3679 at_east_bdry=.false.
3682 if (hmask(i,j) == 1)
then
3684 dxh = g%US%L_to_m*g%dxT(i,j) ; dyh = g%US%L_to_m*g%dyT(i,j) ; dxdyh = g%US%L_to_m**2*g%areaT(i,j)
3686 h_after_uflux(i,j) = h0(i,j)
3688 stencil(:) = h0(i-2:i+2,j)
3694 if (cs%u_face_mask(i-1,j) == 4.)
then
3696 flux_diff_cell = flux_diff_cell + dyh * time_step * cs%u_flux_bdry_val(i-1,j) * &
3697 cs%t_bdry_val(i-1,j) / dxdyh
3701 u_face = 0.5 * (cs%u_shelf(i-1,j-1) + cs%u_shelf(i-1,j))
3703 if (u_face > 0)
then
3707 if (at_west_bdry .AND. (hmask(i-1,j) == 3))
then
3709 flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step * stencil(-1) / dxdyh
3711 elseif (hmask(i-1,j) * hmask(i-2,j) == 1)
then
3712 phi =
slope_limiter(stencil(-1)-stencil(-2), stencil(0)-stencil(-1))
3713 flux_diff_cell = flux_diff_cell + abs(u_face) * dyh* time_step / dxdyh * &
3714 (stencil(-1) - phi * (stencil(-1)-stencil(0))/2)
3720 flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step / dxdyh * stencil(-1)
3724 elseif (u_face < 0)
then
3725 if (hmask(i-1,j) * hmask(i+1,j) == 1)
then
3726 phi =
slope_limiter(stencil(0)-stencil(1), stencil(-1)-stencil(0))
3727 flux_diff_cell = flux_diff_cell - abs(u_face) * dyh * time_step / dxdyh * &
3728 (stencil(0) - phi * (stencil(0)-stencil(-1))/2)
3731 flux_diff_cell = flux_diff_cell - abs(u_face) * dyh * time_step / dxdyh * stencil(0)
3733 if ((hmask(i-1,j) == 0) .OR. (hmask(i-1,j) == 2))
then
3734 flux_enter(i-1,j,2) = abs(u_face) * dyh * time_step * stencil(0)
3744 if (cs%u_face_mask(i+1,j) == 4.)
then
3746 flux_diff_cell = flux_diff_cell + dyh * time_step * cs%u_flux_bdry_val(i+1,j) *&
3747 cs%t_bdry_val(i+1,j)/ dxdyh
3750 u_face = 0.5 * (cs%u_shelf(i,j-1) + cs%u_shelf(i,j))
3752 if (u_face < 0)
then
3754 if (at_east_bdry .AND. (hmask(i+1,j) == 3))
then
3757 flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step * stencil(1) / dxdyh
3759 elseif (hmask(i+1,j) * hmask(i+2,j) == 1)
then
3761 phi =
slope_limiter(stencil(1)-stencil(2), stencil(0)-stencil(1))
3762 flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step / dxdyh * &
3763 (stencil(1) - phi * (stencil(1)-stencil(0))/2)
3769 flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step / dxdyh * stencil(1)
3773 elseif (u_face > 0)
then
3775 if (hmask(i-1,j) * hmask(i+1,j) == 1)
then
3777 phi =
slope_limiter(stencil(0)-stencil(-1), stencil(1)-stencil(0))
3778 flux_diff_cell = flux_diff_cell - abs(u_face) * dyh * time_step / dxdyh * &
3779 (stencil(0) - phi * (stencil(0)-stencil(1))/2)
3785 flux_diff_cell = flux_diff_cell - abs(u_face) * dyh * time_step / dxdyh * stencil(0)
3787 if ((hmask(i+1,j) == 0) .OR. (hmask(i+1,j) == 2))
then
3789 flux_enter(i+1,j,1) = abs(u_face) * dyh * time_step * stencil(0)
3796 h_after_uflux(i,j) = h_after_uflux(i,j) + flux_diff_cell
3800 elseif ((hmask(i,j) == 0) .OR. (hmask(i,j) == 2))
then
3802 if (at_west_bdry .AND. (hmask(i-1,j) == 3))
then
3803 u_face = 0.5 * (cs%u_shelf(i-1,j-1) + cs%u_shelf(i-1,j))
3804 flux_enter(i,j,1) = abs(u_face) * g%US%L_to_m*g%dyT(i,j) * time_step * cs%t_bdry_val(i-1,j)* &
3805 cs%thickness_bdry_val(i+1,j)
3806 elseif (cs%u_face_mask(i-1,j) == 4.)
then
3807 flux_enter(i,j,1) = g%US%L_to_m*g%dyT(i,j) * time_step * cs%u_flux_bdry_val(i-1,j)*cs%t_bdry_val(i-1,j)
3810 if (at_east_bdry .AND. (hmask(i+1,j) == 3))
then
3811 u_face = 0.5 * (cs%u_shelf(i,j-1) + cs%u_shelf(i,j))
3812 flux_enter(i,j,2) = abs(u_face) * g%US%L_to_m*g%dyT(i,j) * time_step * cs%t_bdry_val(i+1,j)* &
3813 cs%thickness_bdry_val(i+1,j)
3814 elseif (cs%u_face_mask(i+1,j) == 4.)
then
3815 flux_enter(i,j,2) = g%US%L_to_m*g%dyT(i,j) * time_step * cs%u_flux_bdry_val(i+1,j) * cs%t_bdry_val(i+1,j)
3846 real,
intent(in) :: time_step
3847 real,
dimension(SZDI_(G),SZDJ_(G)), &
3850 real,
dimension(SZDI_(G),SZDJ_(G)), &
3851 intent(in) :: h_after_uflux
3853 real,
dimension(SZDI_(G),SZDJ_(G)), &
3854 intent(inout) :: h_after_vflux
3856 real,
dimension(SZDI_(G),SZDJ_(G),4), &
3857 intent(inout) :: flux_enter
3878 integer :: i, j, is, ie, js, je, isd, ied, jsd, jed, gjed, gied
3879 integer :: i_off, j_off
3880 logical :: at_north_bdry, at_south_bdry, one_off_west_bdry, one_off_east_bdry
3881 real,
dimension(-2:2) :: stencil
3883 flux_diff_cell, phi, dxh, dyh, dxdyh
3884 character(len=1) :: debug_str
3886 is = g%isc ; ie = g%iec ; js = g%jsc-1 ; je = g%jec+1 ; isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
3887 i_off = g%idg_offset ; j_off = g%jdg_offset
3890 if (((i+i_off) <= g%domain%niglobal+g%domain%nihalo) .AND. &
3891 ((i+i_off) >= g%domain%nihalo+1))
then
3897 if (((j+j_off) <= g%domain%njglobal+g%domain%njhalo) .AND. &
3898 ((j+j_off) >= g%domain%njhalo+1))
then
3900 if (j+j_off == g%domain%njhalo+1)
then
3901 at_south_bdry=.true.
3903 at_south_bdry=.false.
3905 if (j+j_off == g%domain%njglobal+g%domain%njhalo)
then
3906 at_north_bdry=.true.
3908 at_north_bdry=.false.
3911 if (hmask(i,j) == 1)
then
3912 dxh = g%US%L_to_m*g%dxT(i,j) ; dyh = g%US%L_to_m*g%dyT(i,j) ; dxdyh = g%US%L_to_m**2*g%areaT(i,j)
3913 h_after_vflux(i,j) = h_after_uflux(i,j)
3915 stencil(:) = h_after_uflux(i,j-2:j+2)
3920 if (cs%v_face_mask(i,j-1) == 4.)
then
3922 flux_diff_cell = flux_diff_cell + dxh * time_step * cs%v_flux_bdry_val(i,j-1) * &
3923 cs%t_bdry_val(i,j-1)/ dxdyh
3927 v_face = 0.5 * (cs%v_shelf(i-1,j-1) + cs%v_shelf(i,j-1))
3929 if (v_face > 0)
then
3933 if (at_south_bdry .AND. (hmask(i,j-1) == 3))
then
3935 flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step * stencil(-1) / dxdyh
3937 elseif (hmask(i,j-1) * hmask(i,j-2) == 1)
then
3939 phi =
slope_limiter(stencil(-1)-stencil(-2), stencil(0)-stencil(-1))
3940 flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step / dxdyh * &
3941 (stencil(-1) - phi * (stencil(-1)-stencil(0))/2)
3946 flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step / dxdyh * stencil(-1)
3949 elseif (v_face < 0)
then
3951 if (hmask(i,j-1) * hmask(i,j+1) == 1)
then
3952 phi =
slope_limiter(stencil(0)-stencil(1), stencil(-1)-stencil(0))
3953 flux_diff_cell = flux_diff_cell - abs(v_face) * dxh * time_step / dxdyh * &
3954 (stencil(0) - phi * (stencil(0)-stencil(-1))/2)
3956 flux_diff_cell = flux_diff_cell - abs(v_face) * dxh * time_step / dxdyh * stencil(0)
3958 if ((hmask(i,j-1) == 0) .OR. (hmask(i,j-1) == 2))
then
3959 flux_enter(i,j-1,4) = abs(v_face) * dyh * time_step * stencil(0)
3970 if (cs%v_face_mask(i,j+1) == 4.)
then
3972 flux_diff_cell = flux_diff_cell + dxh * time_step * cs%v_flux_bdry_val(i,j+1) *&
3973 cs%t_bdry_val(i,j+1)/ dxdyh
3977 v_face = 0.5 * (cs%v_shelf(i-1,j) + cs%v_shelf(i,j))
3979 if (v_face < 0)
then
3981 if (at_north_bdry .AND. (hmask(i,j+1) == 3))
then
3983 flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step * stencil(1) / dxdyh
3984 elseif (hmask(i,j+1) * hmask(i,j+2) == 1)
then
3985 phi =
slope_limiter(stencil(1)-stencil(2), stencil(0)-stencil(1))
3986 flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step / dxdyh * &
3987 (stencil(1) - phi * (stencil(1)-stencil(0))/2)
3991 flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step / dxdyh * stencil(1)
3994 elseif (v_face > 0)
then
3996 if (hmask(i,j-1) * hmask(i,j+1) == 1)
then
3997 phi =
slope_limiter(stencil(0)-stencil(-1), stencil(1)-stencil(0))
3998 flux_diff_cell = flux_diff_cell - abs(v_face) * dxh * time_step / dxdyh * &
3999 (stencil(0) - phi * (stencil(0)-stencil(1))/2)
4003 flux_diff_cell = flux_diff_cell - abs(v_face) * dxh * time_step / dxdyh * stencil(0)
4004 if ((hmask(i,j+1) == 0) .OR. (hmask(i,j+1) == 2))
then
4005 flux_enter(i,j+1,3) = abs(v_face) * dxh * time_step * stencil(0)
4013 h_after_vflux(i,j) = h_after_vflux(i,j) + flux_diff_cell
4015 elseif ((hmask(i,j) == 0) .OR. (hmask(i,j) == 2))
then
4017 if (at_south_bdry .AND. (hmask(i,j-1) == 3))
then
4018 v_face = 0.5 * (cs%v_shelf(i-1,j-1) + cs%v_shelf(i,j-1))
4019 flux_enter(i,j,3) = abs(v_face) * g%US%L_to_m*g%dxT(i,j) * time_step * cs%t_bdry_val(i,j-1)* &
4020 cs%thickness_bdry_val(i,j-1)
4021 elseif (cs%v_face_mask(i,j-1) == 4.)
then
4022 flux_enter(i,j,3) = g%US%L_to_m*g%dxT(i,j) * time_step * cs%v_flux_bdry_val(i,j-1)*cs%t_bdry_val(i,j-1)
4025 if (at_north_bdry .AND. (hmask(i,j+1) == 3))
then
4026 v_face = 0.5 * (cs%v_shelf(i-1,j) + cs%v_shelf(i,j))
4027 flux_enter(i,j,4) = abs(v_face) * g%US%L_to_m*g%dxT(i,j) * time_step * cs%t_bdry_val(i,j+1)* &
4028 cs%thickness_bdry_val(i,j+1)
4029 elseif (cs%v_face_mask(i,j+1) == 4.)
then
4030 flux_enter(i,j,4) = g%US%L_to_m*g%dxT(i,j) * time_step * cs%v_flux_bdry_val(i,j+1)*cs%t_bdry_val(i,j+1)