8 use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
25 use mom_io, only : write_field, close_file, single_file, multiple
39 use constants_mod
, only: grav
40 use mpp_mod
, only : mpp_sum, mpp_max, mpp_min, mpp_pe, mpp_npes, mpp_sync
43 use time_interp_external_mod
, only : init_external_field, time_interp_external
44 use time_interp_external_mod
, only : time_interp_external_init
45 use time_manager_mod
, only : print_time, time_type_to_real, real_to_time_type
46 implicit none ;
private 48 #include <MOM_memory.h> 49 #ifdef SYMMETRIC_LAND_ICE 50 # define GRID_SYM_ .true. 51 # define NILIMB_SYM_ NIMEMB_SYM_ 52 # define NJLIMB_SYM_ NJMEMB_SYM_ 53 # define ISUMSTART_INT_ CS%grid%iscB+1 54 # define JSUMSTART_INT_ CS%grid%jscB+1 56 # define GRID_SYM_ .false. 57 # define NILIMB_SYM_ NIMEMB_ 58 # define NJLIMB_SYM_ NJMEMB_ 59 # define ISUMSTART_INT_ CS%grid%iscB 60 # define JSUMSTART_INT_ CS%grid%jscB 74 real :: flux_factor = 1.0
76 character(len=128) :: restart_output_dir =
' ' 77 real,
pointer,
dimension(:,:) :: &
78 mass_shelf => null(), &
80 area_shelf_h => null(), &
84 salt_flux => null(), &
88 exch_vel_t => null(), &
89 exch_vel_s => null(), &
93 tflux_shelf => null(), &
115 u_face_mask => null(), &
116 v_face_mask => null(), &
128 u_face_mask_boundary => null(), v_face_mask_boundary => null(), &
129 u_flux_boundary_values => null(), v_flux_boundary_values => null(), &
131 umask => null(), vmask => null(), &
134 calve_mask => null(), &
142 ice_visc_bilinear => null(), &
143 ice_visc_lower_tri => null(), &
144 ice_visc_upper_tri => null(), &
145 thickness_boundary_values => null(), &
146 u_boundary_values => null(), &
147 v_boundary_values => null(), &
148 h_boundary_values => null(), &
150 t_boundary_values => null(), &
152 taub_beta_eff_bilinear => null(), &
154 taub_beta_eff_lower_tri => null(), &
155 taub_beta_eff_upper_tri => null(), &
157 od_rt => null(), float_frac_rt => null(), &
158 od_av => null(), float_frac => null()
179 real :: kd_molec_salt
180 real :: kd_molec_temp
184 real :: col_thick_melt_threshold
185 logical :: mass_from_file
196 logical :: solo_ice_sheet
198 logical :: gl_regularize
200 integer :: n_sub_regularize
210 real :: a_glen_isothermal
213 real :: c_basal_friction
214 real :: n_basal_friction
215 real :: density_ocean_avg
219 real :: thresh_float_col_depth
220 logical :: moving_shelf_front
222 real :: min_thickness_simple_calve
225 real :: input_thickness
230 real :: velocity_update_time_step
235 integer :: velocity_update_sub_counter
236 integer :: velocity_update_counter
237 integer :: nstep_velocity
239 real :: cg_tolerance, nonlinear_tolerance
240 integer :: cg_max_iterations
241 integer :: nonlin_solve_err_mode
245 logical :: use_reproducing_sums
253 logical :: switch_var
255 type(time_type) :: time
258 logical :: shelf_mass_is_dynamic
259 logical :: override_shelf_movement
267 logical :: const_gamma
268 logical :: find_salt_root
269 logical :: constant_sea_level
272 real :: lambda1, lambda2, lambda3
275 integer :: id_melt = -1, id_exch_vel_s = -1, id_exch_vel_t = -1, &
276 id_tfreeze = -1, id_tfl_shelf = -1, &
277 id_thermal_driving = -1, id_haline_driving = -1, &
278 id_u_ml = -1, id_v_ml = -1, id_sbdry = -1, &
279 id_u_shelf = -1, id_v_shelf = -1, id_h_shelf = -1, id_h_mask = -1, &
280 id_u_mask = -1, id_v_mask = -1, id_t_shelf = -1, id_t_mask = -1, &
281 id_surf_elev = -1, id_bathym = -1, id_float_frac = -1, id_col_thick = -1, &
282 id_area_shelf_h = -1, id_od_av = -1, id_float_frac_rt = -1,&
283 id_ustar_shelf = -1, id_shelf_mass = -1, id_mass_flux = -1
288 integer :: id_read_mass
290 integer :: id_read_area
297 logical :: write_output_to_file
308 real,
intent(in) :: num
309 real,
intent(in) :: denom
310 real :: slope_limiter
313 if (denom .eq. 0)
then 315 elseif (num*denom .le. 0)
then 319 slope_limiter = (r+abs(r))/(1+abs(r))
326 real,
dimension(4),
intent(in) :: X
327 real,
dimension(4),
intent(in) :: Y
328 real :: quad_area, p2, q2, a2, c2, b2, d2
335 p2 = (x(4)-x(1))**2 + (y(4)-y(1))**2 ; q2 = (x(3)-x(2))**2 + (y(3)-y(2))**2
336 a2 = (x(3)-x(4))**2 + (y(3)-y(4))**2 ; c2 = (x(1)-x(2))**2 + (y(1)-y(2))**2
337 b2 = (x(2)-x(4))**2 + (y(2)-y(4))**2 ; d2 = (x(3)-x(1))**2 + (y(3)-y(1))**2
338 quad_area = .25 * sqrt(4*p2*q2-(b2+d2-a2-c2)**2)
348 type(
forcing),
intent(inout) :: fluxes
351 type(time_type),
intent(in) :: Time
352 real,
intent(in) :: time_step
358 real,
dimension(SZI_(CS%grid)) :: &
359 Rhoml, & !< Ocean mixed layer density in kg m-3.
360 dR0_dT, & !< Partial derivative of the mixed layer density
366 real,
dimension(:,:),
allocatable :: mass_flux
367 real,
dimension(:,:),
allocatable :: haline_driving
370 real,
parameter :: VK = 0.40
371 real :: ZETA_N = 0.052
373 real,
parameter :: RC = 0.20
380 real,
dimension(:,:),
allocatable :: Sbdry
383 real :: Sbdry1, Sbdry2, S_a, S_b, S_c
386 real :: hBL_neut_h_molec
392 real :: I_n_star, n_star_term, absf
394 real :: dT_ustar, dS_ustar
397 real :: Gam_mol_t, Gam_mol_s
402 real :: Sb_min, Sb_max
403 real :: dS_min, dS_max
405 real :: wB_flux_new, DwB, dDwB_dwB_in
406 real :: I_Gam_T, I_Gam_S, dG_dwB, iDens
407 real :: u_at_h, v_at_h, Isqrt2
408 logical :: Sb_min_set, Sb_max_set
409 character(4) :: stepnum
410 character(2) :: procnum
413 real,
parameter :: c2_3 = 2.0/3.0
414 integer :: i, j, is, ie, js, je, ied, jed, it1, it3, iters_vel_solve
415 real,
parameter :: rho_fw = 1000.0
416 if (.not.
associated(cs))
call mom_error(fatal,
"shelf_calc_flux: "// &
417 "initialize_ice_shelf must be called before shelf_calc_flux.")
422 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; ied = g%ied ; jed = g%jed
423 i_zeta_n = 1.0 / zeta_n
425 i_rholf = 1.0/(cs%Rho0*lf)
427 sc = cs%kv_molec/cs%kd_molec_salt
428 pr = cs%kv_molec/cs%kd_molec_temp
430 rhocp = cs%Rho0 * cs%Cp
431 isqrt2 = 1.0/sqrt(2.0)
434 gam_mol_t = 12.5 * (pr**c2_3) - 6
435 gam_mol_s = 12.5 * (sc**c2_3) - 6
437 idens = 1.0/cs%density_ocean_avg
443 cs%tflux_shelf(:,:) = 0.0; cs%exch_vel_t(:,:) = 0.0
444 cs%lprec(:,:) = 0.0; cs%exch_vel_s(:,:) = 0.0
445 cs%salt_flux(:,:) = 0.0; cs%t_flux(:,:) = 0.0
446 cs%tfreeze(:,:) = 0.0
448 ALLOCATE ( haline_driving(g%ied,g%jed) ); haline_driving(:,:) = 0.0
449 ALLOCATE ( sbdry(g%ied,g%jed) ); sbdry(:,:) = state%sss(:,:)
454 if (cs%shelf_mass_is_dynamic .and. cs%override_shelf_movement)
then 455 cs%time_step = time_step
461 call hchksum (fluxes%frac_shelf_h,
"frac_shelf_h before apply melting", g%HI, haloshift=0)
462 call hchksum (state%sst,
"sst before apply melting", g%HI, haloshift=0)
463 call hchksum (state%sss,
"sss before apply melting", g%HI, haloshift=0)
464 call hchksum (state%u,
"u_ml before apply melting", g%HI, haloshift=0)
465 call hchksum (state%v,
"v_ml before apply melting", g%HI, haloshift=0)
466 call hchksum (state%ocean_mass,
"ocean_mass before apply melting", g%HI, haloshift=0)
472 do i=is,ie ; p_int(i) = cs%g_Earth * cs%mass_shelf(i,j) ;
enddo 476 rhoml(:), is, ie-is+1, cs%eqn_of_state)
477 call calculate_density_derivs(state%sst(:,j), state%sss(:,j), p_int, &
478 dr0_dt, dr0_ds, is, ie-is+1, cs%eqn_of_state)
483 fluxes%ustar_shelf(i,j)= 0.0
488 if ((idens*state%ocean_mass(i,j) > cs%col_thick_melt_threshold) .and. &
489 (cs%area_shelf_h(i,j) > 0.0) .and. &
490 (cs%isthermo) .and. (state%Hml(i,j) > 0.0) )
then 498 u_at_h = state%u(i,j)
499 v_at_h = state%v(i,j)
501 fluxes%ustar_shelf(i,j)= sqrt(cs%cdrag*((u_at_h**2.0 + v_at_h**2.0) +&
504 ustar_h = max(cs%ustar_bg, fluxes%ustar_shelf(i,j))
506 fluxes%ustar_shelf(i,j) = ustar_h
508 if (
associated(state%taux_shelf) .and.
associated(state%tauy_shelf))
then 509 state%taux_shelf(i,j) = ustar_h*ustar_h*cs%Rho0*isqrt2
510 state%tauy_shelf(i,j) = state%taux_shelf(i,j)
515 absf = 0.25*((abs(g%CoriolisBu(i,j)) + abs(g%CoriolisBu(i-1,j-1))) + &
516 (abs(g%CoriolisBu(i,j-1)) + abs(g%CoriolisBu(i-1,j))))
517 if (absf*state%Hml(i,j) <= vk*ustar_h)
then ; hbl_neut = state%Hml(i,j)
518 else ; hbl_neut = (vk*ustar_h) / absf ;
endif 519 hbl_neut_h_molec = zeta_n * ((hbl_neut * ustar_h) / (5.0 * cs%Kv_molec))
522 db_ds = (cs%g_Earth / rhoml(i)) * dr0_ds(i)
523 db_dt = (cs%g_Earth / rhoml(i)) * dr0_dt(i)
524 ln_neut = 0.0 ;
if (hbl_neut_h_molec > 1.0) ln_neut = log(hbl_neut_h_molec)
526 if (cs%find_salt_root)
then 529 s_a = cs%lambda1 * cs%Gamma_T_3EQ * cs%Cp
533 s_b = cs%Gamma_T_3EQ*cs%Cp*(cs%lambda2+cs%lambda3*p_int(i)- &
534 state%sst(i,j))-lf*cs%Gamma_T_3EQ/35.0
535 s_c = lf*(cs%Gamma_T_3EQ/35.0)*state%sss(i,j)
537 sbdry1 = (-s_b + sqrt(s_b*s_b-4*s_a*s_c))/(2*s_a)
538 sbdry2 = (-s_b - sqrt(s_b*s_b-4*s_a*s_c))/(2*s_a)
539 sbdry(i,j) = max(sbdry1, sbdry2)
541 if (sbdry(i,j) < 0.)
then 542 write(*,*)
'state%sss(i,j)',state%sss(i,j)
543 write(*,*)
'S_a, S_b, S_c',s_a, s_b, s_c
544 write(*,*)
'I,J,Sbdry1,Sbdry2',i,j,sbdry1,sbdry2
546 "shelf_calc_flux: Negative salinity (Sbdry).")
550 sbdry(i,j) = state%sss(i,j) ; sb_max_set = .false.
558 dt_ustar = (state%sst(i,j) - cs%tfreeze(i,j)) * ustar_h
559 ds_ustar = (state%sss(i,j) - sbdry(i,j)) * ustar_h
565 if (cs%const_gamma)
then 567 i_gam_t = cs%Gamma_T_3EQ
568 i_gam_s = cs%Gamma_T_3EQ/35.
570 gam_turb = i_vk * (ln_neut + (0.5 * i_zeta_n - 1.0))
571 i_gam_t = 1.0 / (gam_mol_t + gam_turb)
572 i_gam_s = 1.0 / (gam_mol_s + gam_turb)
575 wt_flux = dt_ustar * i_gam_t
576 wb_flux = db_ds * (ds_ustar * i_gam_s) + db_dt * wt_flux
578 if (wb_flux > 0.0)
then 581 n_star_term = (zeta_n/rc) * (hbl_neut * vk) / ustar_h**3
587 i_n_star = sqrt(1.0 + n_star_term * wb_flux)
588 dins_dwb = 0.5 * n_star_term / i_n_star
589 if (hbl_neut_h_molec > i_n_star**2)
then 590 gam_turb = i_vk * ((ln_neut - 2.0*log(i_n_star)) + &
591 (0.5*i_zeta_n*i_n_star - 1.0))
592 dg_dwb = i_vk * ( -2.0 / i_n_star + (0.5 * i_zeta_n)) * dins_dwb
596 gam_turb = i_vk * (0.5 * i_zeta_n*i_n_star - 1.0)
597 dg_dwb = i_vk * (0.5 * i_zeta_n) * dins_dwb
600 if (cs%const_gamma)
then 602 i_gam_t = cs%Gamma_T_3EQ
603 i_gam_s = cs%Gamma_T_3EQ/35.
605 i_gam_t = 1.0 / (gam_mol_t + gam_turb)
606 i_gam_s = 1.0 / (gam_mol_s + gam_turb)
609 wt_flux = dt_ustar * i_gam_t
610 wb_flux_new = db_ds * (ds_ustar * i_gam_s) + db_dt * wt_flux
613 dwb = wb_flux_new - wb_flux
614 if (abs(wb_flux_new - wb_flux) < &
615 1e-4*(abs(wb_flux_new) + abs(wb_flux)))
exit 617 ddwb_dwb_in = -dg_dwb * (db_ds * (ds_ustar * i_gam_s**2) + &
618 db_dt * (dt_ustar * i_gam_t**2)) - 1.0
621 wb_flux_new = wb_flux - dwb / ddwb_dwb_in
625 cs%t_flux(i,j) = rhocp * wt_flux
626 cs%exch_vel_t(i,j) = ustar_h * i_gam_t
627 cs%exch_vel_s(i,j) = ustar_h * i_gam_s
637 if (cs%t_flux(i,j) <= 0.0)
then 638 cs%lprec(i,j) = i_lf * cs%t_flux(i,j)
639 cs%tflux_shelf(i,j) = 0.0
641 if (cs%insulator)
then 643 cs%tflux_shelf(i,j) = 0.0
644 cs%lprec(i,j) = i_lf * (- cs%tflux_shelf(i,j) + cs%t_flux(i,j))
651 cs%lprec(i,j) = cs%t_flux(i,j) / &
652 (lf + cs%CP_Ice * (cs%Tfreeze(i,j) - cs%Temp_Ice))
654 cs%tflux_shelf(i,j) = cs%t_flux(i,j) - lf*cs%lprec(i,j)
663 if (cs%find_salt_root)
then 667 mass_exch = cs%exch_vel_s(i,j) * cs%Rho0
668 sbdry_it = (state%sss(i,j) * mass_exch + cs%Salin_ice * &
669 cs%lprec(i,j)) / (mass_exch + cs%lprec(i,j))
670 ds_it = sbdry_it - sbdry(i,j)
671 if (abs(ds_it) < 1e-4*(0.5*(state%sss(i,j) + sbdry(i,j) + 1.e-10)))
exit 674 if (ds_it < 0.0)
then 675 if (sb_max_set .and. (sbdry(i,j) > sb_max)) &
676 call mom_error(fatal,
"shelf_calc_flux: Irregular iteration for Sbdry (max).")
677 sb_max = sbdry(i,j) ; ds_max = ds_it ; sb_max_set = .true.
679 if (sb_min_set .and. (sbdry(i,j) < sb_min)) &
681 "shelf_calc_flux: Irregular iteration for Sbdry (min).")
682 sb_min = sbdry(i,j) ; ds_min = ds_it ; sb_min_set = .true.
685 if (sb_min_set .and. sb_max_set)
then 687 sbdry(i,j) = sb_min + (sb_max-sb_min) * &
688 (ds_min / (ds_min - ds_max))
690 sbdry(i,j) = sbdry_it
693 sbdry(i,j) = sbdry_it
704 call calculate_tfreeze(state%sss(i,j), p_int(i), cs%tfreeze(i,j), cs%eqn_of_state)
706 cs%exch_vel_t(i,j) = cs%gamma_t
707 cs%t_flux(i,j) = rhocp * cs%exch_vel_t(i,j) * (state%sst(i,j) - cs%tfreeze(i,j))
708 cs%tflux_shelf(i,j) = 0.0
709 cs%lprec(i,j) = i_lf * cs%t_flux(i,j)
723 if (cs%const_gamma)
then 724 fluxes%iceshelf_melt = cs%lprec * (86400.0*365.0/rho_fw) * cs%flux_factor
726 fluxes%iceshelf_melt = cs%lprec * (86400.0*365.0/cs%density_ice) * cs%flux_factor
731 if ((idens*state%ocean_mass(i,j) > cs%col_thick_melt_threshold) .and. &
732 (cs%area_shelf_h(i,j) > 0.0) .and. &
733 (cs%isthermo) .and. (state%Hml(i,j) > 0.0) )
then 738 if ((cs%g_Earth * cs%mass_shelf(i,j)) < cs%Rho0*cs%cutoff_depth* &
741 fluxes%iceshelf_melt(i,j) = 0.0
744 haline_driving(i,j) = (cs%lprec(i,j) * sbdry(i,j)) / &
745 (cs%Rho0 * cs%exch_vel_s(i,j))
762 if (abs(fluxes%iceshelf_melt(i,j))>0.0)
then 763 if (fluxes%ustar_shelf(i,j) == 0.0)
then 764 write(*,*)
'Something is wrong at i,j',i,j
766 "shelf_calc_flux: |melt| > 0 and star_shelf = 0.")
775 ALLOCATE ( mass_flux(g%ied,g%jed) ); mass_flux(:,:) = 0.0
776 mass_flux = (cs%lprec) * cs%area_shelf_h
778 if (cs%shelf_mass_is_dynamic)
then 780 call pass_var(cs%area_shelf_h, g%domain, complete=.false.)
781 call pass_var(cs%mass_shelf, g%domain)
786 if (cs%shelf_mass_is_dynamic .and. cs%override_shelf_movement)
then 787 if (.not. (cs%mass_from_file))
then 802 if (cs%shelf_mass_is_dynamic .and. .not.cs%override_shelf_movement)
then 811 cs%velocity_update_sub_counter = cs%velocity_update_sub_counter+1
813 if (cs%GL_couple .and. .not. cs%solo_ice_sheet)
then 814 call update_od_ffrac (cs, state%ocean_mass, cs%velocity_update_sub_counter, cs%nstep_velocity, cs%time_step, cs%velocity_update_time_step)
819 if (cs%velocity_update_sub_counter .eq. cs%nstep_velocity)
then 821 if (
is_root_pe())
write(*,*)
"ABOUT TO CALL VELOCITY SOLVER" 825 cs%velocity_update_sub_counter = 0
831 if (cs%id_shelf_mass > 0)
call post_data(cs%id_shelf_mass, cs%mass_shelf, cs%diag)
832 if (cs%id_area_shelf_h > 0)
call post_data(cs%id_area_shelf_h, cs%area_shelf_h, cs%diag)
833 if (cs%id_ustar_shelf > 0)
call post_data(cs%id_ustar_shelf, fluxes%ustar_shelf, cs%diag)
834 if (cs%id_melt > 0)
call post_data(cs%id_melt, fluxes%iceshelf_melt, cs%diag)
835 if (cs%id_thermal_driving > 0)
call post_data(cs%id_thermal_driving, (state%sst-cs%tfreeze), cs%diag)
836 if (cs%id_Sbdry > 0)
call post_data(cs%id_Sbdry, sbdry, cs%diag)
837 if (cs%id_haline_driving > 0)
call post_data(cs%id_haline_driving, haline_driving, cs%diag)
838 if (cs%id_mass_flux > 0)
call post_data(cs%id_mass_flux, mass_flux, cs%diag)
839 if (cs%id_u_ml > 0)
call post_data(cs%id_u_ml,state%u,cs%diag)
840 if (cs%id_v_ml > 0)
call post_data(cs%id_v_ml,state%v,cs%diag)
841 if (cs%id_tfreeze > 0)
call post_data(cs%id_tfreeze, cs%tfreeze, cs%diag)
842 if (cs%id_tfl_shelf > 0)
call post_data(cs%id_tfl_shelf, cs%tflux_shelf, cs%diag)
843 if (cs%id_exch_vel_t > 0)
call post_data(cs%id_exch_vel_t, cs%exch_vel_t, cs%diag)
844 if (cs%id_exch_vel_s > 0)
call post_data(cs%id_exch_vel_s, cs%exch_vel_s, cs%diag)
845 if (cs%id_col_thick > 0)
call post_data(cs%id_col_thick, cs%OD_av, cs%diag)
846 if (cs%id_h_shelf > 0)
call post_data(cs%id_h_shelf,cs%h_shelf,cs%diag)
847 if (cs%id_h_mask > 0)
call post_data(cs%id_h_mask,cs%hmask,cs%diag)
848 if (cs%id_u_shelf > 0)
call post_data(cs%id_u_shelf,cs%u_shelf,cs%diag)
849 if (cs%id_v_shelf > 0)
call post_data(cs%id_v_shelf,cs%v_shelf,cs%diag)
850 if (cs%id_float_frac > 0)
call post_data(cs%id_float_frac,cs%float_frac,cs%diag)
851 if (cs%id_OD_av >0)
call post_data(cs%id_OD_av,cs%OD_av,cs%diag)
852 if (cs%id_float_frac_rt>0)
call post_data(cs%id_float_frac_rt,cs%float_frac_rt,cs%diag)
867 real,
intent(in) :: time_step
868 type(
forcing),
intent(inout) :: fluxes
876 if ((cs%hmask(i,j) .eq. 1) .or. (cs%hmask(i,j) .eq. 2))
then 878 if (
associated(fluxes%lprec)) fluxes%lprec(i,j) = 0.0
879 if (
associated(fluxes%sens)) fluxes%sens(i,j) = 0.0
880 if (
associated(fluxes%frac_shelf_h)) fluxes%frac_shelf_h(i,j) = 0.0
881 if (
associated(fluxes%p_surf)) fluxes%p_surf(i,j) = 0.0
882 if (
associated(fluxes%salt_flux)) fluxes%salt_flux(i,j) = 0.0
884 if (cs%lprec(i,j) / cs%density_ice * time_step .lt. cs%h_shelf (i,j))
then 885 cs%h_shelf (i,j) = cs%h_shelf (i,j) - cs%lprec(i,j) / cs%density_ice * time_step
892 cs%h_shelf(i,j) = 0.0
894 cs%area_shelf_h(i,j) = 0.0
900 call pass_var(cs%area_shelf_h, g%domain)
907 if ((cs%hmask(i,j) .eq. 1) .or. (cs%hmask(i,j) .eq. 2))
then 908 cs%mass_shelf(i,j) = cs%h_shelf(i,j)*cs%density_ice
913 call pass_var(cs%mass_shelf, g%domain)
916 call hchksum (cs%h_shelf,
"h_shelf after change thickness using melt", g%HI, haloshift=0)
917 call hchksum (cs%mass_shelf,
"mass_shelf after change thickness using melt", g%HI, haloshift=0)
926 type(
surface),
intent(inout) :: state
927 type(
forcing),
intent(inout) :: fluxes
934 real :: delta_mass_shelf
939 real :: mean_melt_flux
942 type(time_type) :: Time0
943 real,
dimension(:,:),
allocatable,
target :: last_mass_shelf
945 real,
dimension(:,:),
allocatable,
target :: last_h_shelf
947 real,
dimension(:,:),
allocatable,
target :: last_hmask
949 real,
dimension(:,:),
allocatable,
target :: last_area_shelf_h
952 real,
parameter :: rho_fw = 1000.0
953 integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
954 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
955 isd = g%isd ; jsd = g%jsd ; ied = g%ied ; jed = g%jed
957 irho0 = 1.0 / cs%Rho0
961 if (cs%shelf_mass_is_dynamic)
then 962 do j=jsd,jed ;
do i=isd,ied
963 if (g%areaT(i,j) > 0.0) &
964 fluxes%frac_shelf_h(i,j) = cs%area_shelf_h(i,j) / g%areaT(i,j)
967 do j=jsd,jed ;
do i=isd,ied-1
968 fluxes%frac_shelf_u(i,j) = 0.0
969 if ((g%areaT(i,j) + g%areaT(i+1,j) > 0.0)) &
970 fluxes%frac_shelf_u(i,j) = ((cs%area_shelf_h(i,j) + cs%area_shelf_h(i+1,j)) / &
971 (g%areaT(i,j) + g%areaT(i+1,j)))
972 fluxes%rigidity_ice_u(i,j) = (cs%kv_ice / cs%density_ice) * &
973 min(cs%mass_shelf(i,j), cs%mass_shelf(i+1,j))
975 do j=jsd,jed-1 ;
do i=isd,ied
977 fluxes%frac_shelf_v(i,j) = 0.0
978 if ((g%areaT(i,j) + g%areaT(i,j+1) > 0.0)) &
979 fluxes%frac_shelf_v(i,j) = ((cs%area_shelf_h(i,j) + cs%area_shelf_h(i,j+1)) / &
980 (g%areaT(i,j) + g%areaT(i,j+1)))
981 fluxes%rigidity_ice_v(i,j) = (cs%kv_ice / cs%density_ice) * &
982 max(cs%mass_shelf(i,j), cs%mass_shelf(i,j+1))
984 call pass_vector(fluxes%frac_shelf_u, fluxes%frac_shelf_v, g%domain,
to_all, cgrid_ne)
989 do j=jsd,jed ;
do i=isd,ied-1
990 fluxes%rigidity_ice_u(i,j) = (cs%kv_ice / cs%density_ice) * &
991 min(cs%mass_shelf(i,j), cs%mass_shelf(i+1,j))
994 do j=jsd,jed-1 ;
do i=isd,ied
995 fluxes%rigidity_ice_v(i,j) = (cs%kv_ice / cs%density_ice) * &
996 max(cs%mass_shelf(i,j), cs%mass_shelf(i,j+1))
1001 if (
associated(state%taux_shelf))
then 1002 call uchksum(state%taux_shelf,
"taux_shelf", g%HI, haloshift=0)
1004 if (
associated(state%tauy_shelf))
then 1005 call vchksum(state%tauy_shelf,
"tauy_shelf", g%HI, haloshift=0)
1006 call vchksum(fluxes%rigidity_ice_u,
"rigidity_ice_u", g%HI, haloshift=0)
1007 call vchksum(fluxes%rigidity_ice_v,
"rigidity_ice_v", g%HI, haloshift=0)
1008 call vchksum(fluxes%frac_shelf_u,
"frac_shelf_u", g%HI, haloshift=0)
1009 call vchksum(fluxes%frac_shelf_v,
"frac_shelf_v", g%HI, haloshift=0)
1013 if (
associated(state%taux_shelf) .and.
associated(state%tauy_shelf))
then 1014 call pass_vector(state%taux_shelf, state%tauy_shelf, g%domain,
to_all, cgrid_ne)
1017 if (
associated(fluxes%sw_vis_dir)) fluxes%sw_vis_dir = 0.0
1018 if (
associated(fluxes%sw_vis_dif)) fluxes%sw_vis_dif = 0.0
1019 if (
associated(fluxes%sw_nir_dir)) fluxes%sw_nir_dir = 0.0
1020 if (
associated(fluxes%sw_nir_dif)) fluxes%sw_nir_dif = 0.0
1022 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
1023 frac_area = fluxes%frac_shelf_h(i,j)
1024 if (frac_area > 0.0)
then 1026 taux2 = 0.0 ; tauy2 = 0.0
1027 asu1 = fluxes%frac_shelf_u(i-1,j) * (g%areaT(i-1,j) + g%areaT(i,j))
1028 asu2 = fluxes%frac_shelf_u(i,j) * (g%areaT(i,j) + g%areaT(i+1,j))
1029 asv1 = fluxes%frac_shelf_v(i,j-1) * (g%areaT(i,j-1) + g%areaT(i,j))
1030 asv2 = fluxes%frac_shelf_v(i,j) * (g%areaT(i,j) + g%areaT(i,j+1))
1031 if ((asu1 + asu2 > 0.0) .and.
associated(state%taux_shelf)) &
1032 taux2 = (asu1 * state%taux_shelf(i-1,j)**2 + &
1033 asu2 * state%taux_shelf(i,j)**2 ) / (asu1 + asu2)
1034 if ((asv1 + asv2 > 0.0) .and.
associated(state%tauy_shelf)) &
1035 tauy2 = (asv1 * state%tauy_shelf(i,j-1)**2 + &
1036 asv2 * state%tauy_shelf(i,j)**2 ) / (asv1 + asv2)
1042 if (
associated(fluxes%sw)) fluxes%sw(i,j) = 0.0
1043 if (
associated(fluxes%lw)) fluxes%lw(i,j) = 0.0
1044 if (
associated(fluxes%latent)) fluxes%latent(i,j) = 0.0
1045 if (
associated(fluxes%evap)) fluxes%evap(i,j) = 0.0
1046 if (
associated(fluxes%lprec))
then 1047 if (cs%lprec(i,j) > 0.0 )
then 1048 fluxes%lprec(i,j) = frac_area*cs%lprec(i,j)*cs%flux_factor
1050 fluxes%lprec(i,j) = 0.0
1051 fluxes%evap(i,j) = frac_area*cs%lprec(i,j)*cs%flux_factor
1060 if (
associated(state%frazil))
then 1061 fraz = state%frazil(i,j) / cs%time_step / cs%Lat_fusion
1062 if (
associated(fluxes%evap)) fluxes%evap(i,j) = fluxes%evap(i,j) - fraz
1063 cs%lprec(i,j)=cs%lprec(i,j) - fraz
1064 state%frazil(i,j) = 0.0
1067 if (
associated(fluxes%sens)) fluxes%sens(i,j) = -frac_area*cs%t_flux(i,j)*cs%flux_factor
1068 if (
associated(fluxes%salt_flux)) fluxes%salt_flux(i,j) = frac_area * cs%salt_flux(i,j)*cs%flux_factor
1069 if (
associated(fluxes%p_surf)) fluxes%p_surf(i,j) = frac_area * cs%g_Earth * cs%mass_shelf(i,j)
1071 if (
associated(fluxes%p_surf_full) ) fluxes%p_surf_full(i,j) = &
1072 frac_area * cs%g_Earth * cs%mass_shelf(i,j)
1083 if (cs%constant_sea_level)
then 1085 if (.not.
associated(fluxes%salt_flux))
allocate(fluxes%salt_flux(ie,je))
1086 if (.not.
associated(fluxes%vprec))
allocate(fluxes%vprec(ie,je))
1087 fluxes%salt_flux(:,:) = 0.0; fluxes%vprec(:,:) = 0.0
1089 mean_melt_flux = 0.0; sponge_area = 0.0
1090 do j=js,je ;
do i=is,ie
1091 frac_area = fluxes%frac_shelf_h(i,j)
1092 if (frac_area > 0.0)
then 1093 mean_melt_flux = mean_melt_flux + (cs%lprec(i,j)) * cs%area_shelf_h(i,j)
1096 if (g%geoLonT(i,j) >= 790.0 .AND. g%geoLonT(i,j) <= 800.0)
then 1097 sponge_area = sponge_area + g%areaT(i,j)
1102 if (cs%shelf_mass_is_dynamic .and. cs%override_shelf_movement .and. &
1103 cs%mass_from_file)
then 1104 t0 = time_type_to_real(cs%Time) - cs%time_step
1108 time0 = real_to_time_type(t0)
1109 allocate(last_mass_shelf(isd:ied,jsd:jed))
1110 allocate(last_h_shelf(isd:ied,jsd:jed))
1111 allocate(last_area_shelf_h(isd:ied,jsd:jed))
1112 allocate(last_hmask(isd:ied,jsd:jed))
1113 last_hmask(:,:) = cs%hmask(:,:); last_area_shelf_h(:,:) = cs%area_shelf_h(:,:)
1114 call time_interp_external(cs%id_read_mass, time0, last_mass_shelf)
1115 last_h_shelf = last_mass_shelf/cs%density_ice
1118 if (cs%min_thickness_simple_calve > 0.0)
then 1121 last_mass_shelf = last_h_shelf * cs%density_ice
1124 shelf_mass0 = 0.0; shelf_mass1 = 0.0
1126 do j=js,je ;
do i=is,ie
1128 if (((1.0/cs%density_ocean_avg)*state%ocean_mass(i,j) > 0.1) .and. &
1129 (cs%area_shelf_h(i,j) > 0.0))
then 1131 shelf_mass0 = shelf_mass0 + (last_mass_shelf(i,j) * cs%area_shelf_h(i,j))
1132 shelf_mass1 = shelf_mass1 + (cs%mass_shelf(i,j) * cs%area_shelf_h(i,j))
1136 call mpp_sum(shelf_mass0);
call mpp_sum(shelf_mass1)
1137 delta_mass_shelf = (shelf_mass1 - shelf_mass0)/cs%time_step
1142 delta_mass_shelf = 0.0
1145 delta_mass_shelf = 0.0
1148 call mpp_sum(mean_melt_flux)
1149 call mpp_sum(sponge_area)
1152 mean_melt_flux = (mean_melt_flux+delta_mass_shelf) / sponge_area
1155 do j=js,je ;
do i=is,ie
1157 if (g%geoLonT(i,j) >= 790.0 .AND. g%geoLonT(i,j) <= 800.0)
then 1158 fluxes%vprec(i,j) = -mean_melt_flux * cs%density_ice/1000.
1159 fluxes%sens(i,j) = fluxes%vprec(i,j) * cs%Cp * cs%T0
1160 fluxes%salt_flux(i,j) = fluxes%vprec(i,j) * cs%S0*1.0e-3
1165 if (is_root_pe())
write(*,*)
'Mean melt flux (kg/(m^2 s)),dt',mean_melt_flux,cs%time_step
1166 call mom_forcing_chksum(
"After constant sea level", fluxes, g, haloshift=0)
1174 if (cs%shelf_mass_is_dynamic)
then 1175 do j=g%jsc,g%jec ;
do i=g%isc-1,g%iec
1176 fluxes%rigidity_ice_u(i,j) = (cs%kv_ice / cs%density_ice) * &
1177 max(cs%mass_shelf(i,j), cs%mass_shelf(i+1,j))
1180 do j=g%jsc-1,g%jec ;
do i=g%isc,g%iec
1181 fluxes%rigidity_ice_v(i,j) = (cs%kv_ice / cs%density_ice) * &
1182 max(cs%mass_shelf(i,j), cs%mass_shelf(i,j+1))
1190 subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Time_in, solo_ice_sheet_in)
1191 type(param_file_type),
intent(in) :: param_file
1192 type(ocean_grid_type),
pointer :: ocn_grid
1193 type(time_type),
intent(inout) :: Time
1195 type(diag_ctrl),
target,
intent(in) :: diag
1196 type(forcing),
optional,
intent(inout) :: fluxes
1197 type(time_type),
optional,
intent(in) :: Time_in
1198 logical,
optional,
intent(in) :: solo_ice_sheet_in
1200 type(ocean_grid_type),
pointer :: G, OG
1201 type(directories) :: dirs
1203 type(dyn_horgrid_type),
pointer :: dG => null()
1204 real :: cdrag, drag_bg_vel
1205 logical :: new_sim, save_IC, var_force
1207 #include "version_variable.h" 1208 character(len=200) :: config
1209 character(len=200) :: IC_file,filename,inputdir
1210 character(len=40) :: var_name
1211 character(len=40) :: mdl =
"MOM_ice_shelf" 1212 character(len=2) :: procnum
1213 integer :: i, j, is, ie, js, je, isd, ied, jsd, jed, Isdq, Iedq, Jsdq, Jedq, iters
1214 integer :: wd_halos(2)
1215 logical :: solo_ice_sheet, read_TideAmp
1216 character(len=240) :: Tideamp_file
1218 if (
associated(cs))
then 1219 call mom_error(fatal,
"MOM_ice_shelf.F90, initialize_ice_shelf: "// &
1220 "called with an associated control structure.")
1228 call get_mom_input(dirs=dirs)
1232 call mom_domains_init(cs%grid%domain, param_file, min_halo=wd_halos, symmetric=grid_sym_)
1235 call mom_grid_init(cs%grid, param_file)
1237 call create_dyn_horgrid(dg, cs%grid%HI)
1238 call clone_mom_domain(cs%grid%Domain, dg%Domain)
1240 call set_grid_metrics(dg, param_file)
1244 if (
associated(ocn_grid)) cs%ocn_grid => ocn_grid
1250 if (is_root_pe())
then 1251 write(0,*)
'OG: ', og%isd, og%isc, og%iec, og%ied, og%jsd, og%jsc, og%jsd, og%jed
1252 write(0,*)
'IG: ', g%isd, g%isc, g%iec, g%ied, g%jsd, g%jsc, g%jsd, g%jed
1260 solo_ice_sheet = .false.
1261 if (
present(solo_ice_sheet_in)) solo_ice_sheet = solo_ice_sheet_in
1262 cs%solo_ice_sheet = solo_ice_sheet
1264 if (
present(time_in)) time = time_in
1266 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1267 isd = g%isd ; jsd = g%jsd ; ied = g%ied ; jed = g%jed
1268 isdq = g%IsdB ; iedq = g%IedB ; jsdq = g%JsdB ; jedq = g%JedB
1270 cs%Lat_fusion = 3.34e5
1271 cs%override_shelf_movement = .false.
1273 cs%use_reproducing_sums = .false.
1274 cs%switch_var = .false.
1276 call log_version(param_file, mdl, version,
"")
1277 call get_param(param_file, mdl,
"DEBUG_IS", cs%debug, default=.false.)
1278 call get_param(param_file, mdl,
"DYNAMIC_SHELF_MASS", cs%shelf_mass_is_dynamic, &
1279 "If true, the ice sheet mass can evolve with time.", &
1281 if (cs%shelf_mass_is_dynamic)
then 1282 call get_param(param_file, mdl,
"OVERRIDE_SHELF_MOVEMENT", cs%override_shelf_movement, &
1283 "If true, user provided code specifies the ice-shelf \n"//&
1284 "movement instead of the dynamic ice model.", default=.false.)
1285 call get_param(param_file, mdl,
"GROUNDING_LINE_INTERPOLATE", cs%GL_regularize, &
1286 "THIS PARAMETER NEEDS A DESCRIPTION.", default=.false.)
1287 call get_param(param_file, mdl,
"GROUNDING_LINE_INTERP_SUBGRID_N", cs%n_sub_regularize, &
1288 "THIS PARAMETER NEEDS A DESCRIPTION.", default=0)
1289 call get_param(param_file, mdl,
"GROUNDING_LINE_COUPLE", cs%GL_couple, &
1290 "THIS PARAMETER NEEDS A DESCRIPTION.", default=.false.)
1291 if (cs%GL_regularize) cs%GL_couple = .false.
1292 if (cs%GL_regularize .and. (cs%n_sub_regularize.eq.0))
call mom_error (fatal, &
1293 "GROUNDING_LINE_INTERP_SUBGRID_N must be a positive integer if GL regularization is used")
1295 call get_param(param_file, mdl,
"SHELF_THERMO", cs%isthermo, &
1296 "If true, use a thermodynamically interactive ice shelf.", &
1298 call get_param(param_file, mdl,
"SHELF_THREE_EQN", cs%threeeq, &
1299 "If true, use the three equation expression of \n"//&
1300 "consistency to calculate the fluxes at the ice-ocean \n"//&
1301 "interface.", default=.true.)
1302 call get_param(param_file, mdl,
"SHELF_INSULATOR", cs%insulator, &
1303 "If true, the ice shelf is a perfect insulatior \n"//&
1304 "(no conduction).", default=.false.)
1305 call get_param(param_file, mdl,
"MELTING_CUTOFF_DEPTH", cs%cutoff_depth, &
1306 "Depth above which the melt is set to zero (it must be >= 0) \n"//&
1307 "Default value won't affect the solution.", default=0.0)
1308 if (cs%cutoff_depth < 0.) &
1309 call mom_error(warning,
"Initialize_ice_shelf: MELTING_CUTOFF_DEPTH must be >= 0.")
1311 call get_param(param_file, mdl,
"CONST_SEA_LEVEL", cs%constant_sea_level, &
1312 "If true, apply evaporative, heat and salt fluxes in \n"//&
1313 "the sponge region. This will avoid a large increase \n"//&
1314 "in sea level. This option is needed for some of the \n"//&
1315 "ISOMIP+ experiments (Ocean3 and Ocean4). \n"//&
1316 "IMPORTANT: it is not currently possible to do \n"//&
1317 "prefect restarts using this flag.", default=.false.)
1319 call get_param(param_file, mdl,
"ISOMIP_S_SUR_SPONGE", &
1320 cs%S0,
"Surface salinity in the resoring region.", &
1321 default=33.8, do_not_log=.true.)
1323 call get_param(param_file, mdl,
"ISOMIP_T_SUR_SPONGE", &
1324 cs%T0,
"Surface temperature in the resoring region.", &
1325 default=-1.9, do_not_log=.true.)
1327 call get_param(param_file, mdl,
"SHELF_3EQ_GAMMA", cs%const_gamma, &
1328 "If true, user specifies a constant nondimensional heat-transfer coefficient \n"//&
1329 "(GAMMA_T_3EQ), from which the salt-transfer coefficient is then computed \n"//&
1330 " as GAMMA_T_3EQ/35. This is used with SHELF_THREE_EQN.", default=.false.)
1331 if (cs%const_gamma)
call get_param(param_file, mdl,
"SHELF_3EQ_GAMMA_T", cs%Gamma_T_3EQ, &
1332 "Nondimensional heat-transfer coefficient.",default=2.2e-2, &
1333 units=
"nondim.", fail_if_missing=.true.)
1335 call get_param(param_file, mdl,
"ICE_SHELF_MASS_FROM_FILE", &
1336 cs%mass_from_file,
"Read the mass of the & 1337 ice shelf (every time step) from a file.", default=.false.)
1340 call get_param(param_file, mdl,
"SHELF_S_ROOT", cs%find_salt_root, &
1341 "If SHELF_S_ROOT = True, salinity at the ice/ocean interface (Sbdry) \n "//&
1342 "is computed from a quadratic equation. Otherwise, the previous \n"//&
1343 "interactive method to estimate Sbdry is used.", default=.false.)
1344 if (cs%find_salt_root)
then 1345 call get_param(param_file, mdl,
"TFREEZE_S0_P0",cs%lambda1, &
1346 "this is the freezing potential temperature at \n"//&
1347 "S=0, P=0.", units=
"deg C", default=0.0, do_not_log=.true.)
1348 call get_param(param_file, mdl,
"DTFREEZE_DS",cs%lambda1, &
1349 "this is the derivative of the freezing potential \n"//&
1350 "temperature with salinity.", &
1351 units=
"deg C PSU-1", default=-0.054, do_not_log=.true.)
1352 call get_param(param_file, mdl,
"DTFREEZE_DP",cs%lambda3, &
1353 "this is the derivative of the freezing potential \n"//&
1354 "temperature with pressure.", &
1355 units=
"deg C Pa-1", default=0.0, do_not_log=.true.)
1359 if (.not.cs%threeeq) &
1360 call get_param(param_file, mdl,
"SHELF_2EQ_GAMMA_T", cs%gamma_t, &
1361 "If SHELF_THREE_EQN is false, this the fixed turbulent \n"//&
1362 "exchange velocity at the ice-ocean interface.", &
1363 units=
"m s-1", fail_if_missing=.true.)
1365 call get_param(param_file, mdl,
"G_EARTH", cs%g_Earth, &
1366 "The gravitational acceleration of the Earth.", &
1367 units=
"m s-2", default = 9.80)
1368 call get_param(param_file, mdl,
"C_P", cs%Cp, &
1369 "The heat capacity of sea water.", units=
"J kg-1 K-1", &
1370 fail_if_missing=.true.)
1371 call get_param(param_file, mdl,
"RHO_0", cs%Rho0, &
1372 "The mean ocean density used with BOUSSINESQ true to \n"//&
1373 "calculate accelerations and the mass for conservation \n"//&
1374 "properties, or with BOUSSINSEQ false to convert some \n"//&
1375 "parameters from vertical units of m to kg m-2.", &
1376 units=
"kg m-3", default=1035.0)
1377 call get_param(param_file, mdl,
"C_P_ICE", cs%Cp_ice, &
1378 "The heat capacity of ice.", units=
"J kg-1 K-1", &
1381 call get_param(param_file, mdl,
"ICE_SHELF_FLUX_FACTOR", cs%flux_factor, &
1382 "Non-dimensional factor applied to shelf thermodynamic \n"//&
1383 "fluxes.", units=
"none", default=1.0)
1385 call get_param(param_file, mdl,
"KV_ICE", cs%kv_ice, &
1386 "The viscosity of the ice.", units=
"m2 s-1", default=1.0e10)
1387 call get_param(param_file, mdl,
"KV_MOLECULAR", cs%kv_molec, &
1388 "The molecular kinimatic viscosity of sea water at the \n"//&
1389 "freezing temperature.", units=
"m2 s-1", default=1.95e-6)
1390 call get_param(param_file, mdl,
"ICE_SHELF_SALINITY", cs%Salin_ice, &
1391 "The salinity of the ice inside the ice shelf.", units=
"PSU", &
1393 call get_param(param_file, mdl,
"ICE_SHELF_TEMPERATURE", cs%Temp_ice, &
1394 "The temperature at the center of the ice shelf.", &
1395 units =
"degC", default=-15.0)
1396 call get_param(param_file, mdl,
"KD_SALT_MOLECULAR", cs%kd_molec_salt, &
1397 "The molecular diffusivity of salt in sea water at the \n"//&
1398 "freezing point.", units=
"m2 s-1", default=8.02e-10)
1399 call get_param(param_file, mdl,
"KD_TEMP_MOLECULAR", cs%kd_molec_temp, &
1400 "The molecular diffusivity of heat in sea water at the \n"//&
1401 "freezing point.", units=
"m2 s-1", default=1.41e-7)
1402 call get_param(param_file, mdl,
"RHO_0", cs%density_ocean_avg, &
1403 "avg ocean density used in floatation cond", &
1404 units=
"kg m-3", default=1035.)
1405 call get_param(param_file, mdl,
"DT_FORCING", cs%time_step, &
1406 "The time step for changing forcing, coupling with other \n"//&
1407 "components, or potentially writing certain diagnostics. \n"//&
1408 "The default value is given by DT.", units=
"s", default=0.0)
1409 call get_param(param_file, mdl,
"SHELF_DIAG_TIMESTEP", cs%velocity_update_time_step, &
1410 "A timestep to use for diagnostics of the shelf.", default=0.0)
1412 call get_param(param_file, mdl,
"COL_THICK_MELT_THRESHOLD", cs%col_thick_melt_threshold, &
1413 "The minimum ML thickness where melting is allowed.", units=
"m", &
1416 call get_param(param_file, mdl,
"READ_TIDEAMP", read_tideamp, &
1417 "If true, read a file (given by TIDEAMP_FILE) containing \n"//&
1418 "the tidal amplitude with INT_TIDE_DISSIPATION.", default=.false.)
1420 call safe_alloc_ptr(cs%utide,isd,ied,jsd,jed) ; cs%utide(:,:) = 0.0
1422 if (read_tideamp)
then 1423 call get_param(param_file, mdl,
"TIDEAMP_FILE", tideamp_file, &
1424 "The path to the file containing the spatially varying \n"//&
1425 "tidal amplitudes.", &
1426 default=
"tideamp.nc")
1427 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
1428 inputdir = slasher(inputdir)
1429 tideamp_file = trim(inputdir) // trim(tideamp_file)
1430 call read_data(tideamp_file,
'tideamp',cs%utide,domain=g%domain%mpp_domain,timelevel=1)
1432 call get_param(param_file, mdl,
"UTIDE", utide, &
1433 "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", &
1434 units=
"m s-1", default=0.0)
1438 call eos_init(param_file, cs%eqn_of_state)
1442 if (cs%shelf_mass_is_dynamic .and. .not.cs%override_shelf_movement)
then 1444 call get_param(param_file, mdl,
"A_GLEN_ISOTHERM", cs%A_glen_isothermal, &
1445 "Ice viscosity parameter in Glen's Law", &
1446 units=
"Pa -1/3 a", default=9.461e-18)
1447 call get_param(param_file, mdl,
"GLEN_EXPONENT", cs%n_glen, &
1448 "nonlinearity exponent in Glen's Law", &
1449 units=
"none", default=3.)
1450 call get_param(param_file, mdl,
"MIN_STRAIN_RATE_GLEN", cs%eps_glen_min, &
1451 "min. strain rate to avoid infinite Glen's law viscosity", &
1452 units=
"a-1", default=1.e-12)
1453 call get_param(param_file, mdl,
"BASAL_FRICTION_COEFF", cs%C_basal_friction, &
1454 "ceofficient in sliding law \tau_b = C u^(n_basal_friction)", &
1455 units=
"Pa (m-a)-(n_basal_friction)", fail_if_missing=.true.)
1456 call get_param(param_file, mdl,
"BASAL_FRICTION_EXP", cs%n_basal_friction, &
1457 "exponent in sliding law \tau_b = C u^(m_slide)", &
1458 units=
"none", fail_if_missing=.true.)
1459 call get_param(param_file, mdl,
"DENSITY_ICE", cs%density_ice, &
1460 "A typical density of ice.", units=
"kg m-3", default=917.0)
1462 call get_param(param_file, mdl,
"INPUT_FLUX_ICE_SHELF", cs%input_flux, &
1463 "volume flux at upstream boundary", &
1464 units=
"m2 s-1", default=0.)
1465 call get_param(param_file, mdl,
"INPUT_THICK_ICE_SHELF", cs%input_thickness, &
1466 "flux thickness at upstream boundary", &
1467 units=
"m", default=1000.)
1468 call get_param(param_file, mdl,
"ICE_VELOCITY_TIMESTEP", cs%velocity_update_time_step, &
1469 "seconds between ice velocity calcs", units=
"s", &
1470 fail_if_missing=.true.)
1472 call get_param(param_file, mdl,
"CONJUGATE_GRADIENT_TOLERANCE", cs%cg_tolerance, &
1473 "tolerance in CG solver, relative to initial residual", default=1.e-6)
1474 call get_param(param_file, mdl,
"ICE_NONLINEAR_TOLERANCE", &
1475 cs%nonlinear_tolerance,
"nonlin tolerance in iterative velocity solve",default=1.e-6)
1476 call get_param(param_file, mdl,
"CONJUGATE_GRADIENT_MAXIT", cs%cg_max_iterations, &
1477 "max iteratiions in CG solver", default=2000)
1478 call get_param(param_file, mdl,
"THRESH_FLOAT_COL_DEPTH", cs%thresh_float_col_depth, &
1479 "min ocean thickness to consider ice *floating*; \n"// &
1480 "will only be important with use of tides", &
1481 units=
"m",default=1.e-3)
1483 call get_param(param_file, mdl,
"SHELF_MOVING_FRONT", cs%moving_shelf_front, &
1484 "whether or not to advance shelf front (and calve..)")
1485 call get_param(param_file, mdl,
"CALVE_TO_MASK", cs%calve_to_mask, &
1486 "if true, do not allow an ice shelf where prohibited by a mask")
1487 call get_param(param_file, mdl,
"ICE_SHELF_CFL_FACTOR", cs%CFL_factor, &
1488 "limit timestep as a factor of min (\Delta x / u); \n"// &
1489 "only important for ice-only model", &
1491 call get_param(param_file, mdl,
"NONLIN_SOLVE_ERR_MODE", cs%nonlin_solve_err_mode, &
1492 "choose whether nonlin error in vel solve is based on nonlinear residual (1) \n"// &
1493 "or relative change since last iteration (2)", &
1497 if (cs%debug) cs%use_reproducing_sums = .true.
1499 cs%nstep_velocity = floor(cs%velocity_update_time_step / cs%time_step)
1500 cs%velocity_update_counter = 0
1501 cs%velocity_update_sub_counter = 0
1503 cs%nstep_velocity = 0
1505 call get_param(param_file, mdl,
"DENSITY_ICE", cs%density_ice, &
1506 "A typical density of ice.", units=
"kg m-3", default=900.0)
1509 call get_param(param_file, mdl,
"MIN_THICKNESS_SIMPLE_CALVE", &
1510 cs%min_thickness_simple_calve, &
1511 "min thickness rule for VERY simple calving law",&
1512 units=
"m", default=0.0)
1514 call get_param(param_file, mdl,
"WRITE_OUTPUT_TO_FILE", &
1515 cs%write_output_to_file,
"for debugging purposes",default=.false.)
1517 call get_param(param_file, mdl,
"USTAR_SHELF_BG", cs%ustar_bg, &
1518 "The minimum value of ustar under ice sheves.", units=
"m s-1", &
1520 call get_param(param_file, mdl,
"CDRAG_SHELF", cdrag, &
1521 "CDRAG is the drag coefficient relating the magnitude of \n"//&
1522 "the velocity field to the surface stress.", units=
"nondim", &
1525 if (cs%ustar_bg <= 0.0)
then 1526 call get_param(param_file, mdl,
"DRAG_BG_VEL_SHELF", drag_bg_vel, &
1527 "DRAG_BG_VEL is either the assumed bottom velocity (with \n"//&
1528 "LINEAR_DRAG) or an unresolved velocity that is \n"//&
1529 "combined with the resolved velocity to estimate the \n"//&
1530 "velocity magnitude.", units=
"m s-1", default=0.0)
1531 if (cs%cdrag*drag_bg_vel > 0.0) cs%ustar_bg = sqrt(cs%cdrag)*drag_bg_vel
1536 allocate( cs%mass_shelf(isd:ied,jsd:jed) ) ; cs%mass_shelf(:,:) = 0.0
1537 allocate( cs%area_shelf_h(isd:ied,jsd:jed) ) ; cs%area_shelf_h(:,:) = 0.0
1538 allocate( cs%t_flux(isd:ied,jsd:jed) ) ; cs%t_flux(:,:) = 0.0
1539 allocate( cs%lprec(isd:ied,jsd:jed) ) ; cs%lprec(:,:) = 0.0
1540 allocate( cs%salt_flux(isd:ied,jsd:jed) ) ; cs%salt_flux(:,:) = 0.0
1542 allocate( cs%tflux_shelf(isd:ied,jsd:jed) ) ; cs%tflux_shelf(:,:) = 0.0
1543 allocate( cs%tfreeze(isd:ied,jsd:jed) ) ; cs%tfreeze(:,:) = 0.0
1544 allocate( cs%exch_vel_s(isd:ied,jsd:jed) ) ; cs%exch_vel_s(:,:) = 0.0
1545 allocate( cs%exch_vel_t(isd:ied,jsd:jed) ) ; cs%exch_vel_t(:,:) = 0.0
1547 allocate ( cs%h_shelf(isd:ied,jsd:jed) ) ; cs%h_shelf(:,:) = 0.0
1548 allocate ( cs%hmask(isd:ied,jsd:jed) ) ; cs%hmask(:,:) = -2.0
1552 allocate ( cs%t_shelf(isd:ied,jsd:jed) ) ; cs%t_shelf(:,:) = -10.0
1553 allocate ( cs%t_boundary_values(isd:ied,jsd:jed) ) ; cs%t_boundary_values(:,:) = -15.0
1554 allocate ( cs%tmask(isdq:iedq,jsdq:jedq) ) ; cs%tmask(:,:) = -1.0
1556 if (cs%shelf_mass_is_dynamic .and. .not.cs%override_shelf_movement)
then 1558 allocate ( cs%u_shelf(isdq:iedq,jsdq:jedq) ) ; cs%u_shelf(:,:) = 0.0
1559 allocate ( cs%v_shelf(isdq:iedq,jsdq:jedq) ) ; cs%v_shelf(:,:) = 0.0
1560 allocate ( cs%u_boundary_values(isdq:iedq,jsdq:jedq) ) ; cs%u_boundary_values(:,:) = 0.0
1561 allocate ( cs%v_boundary_values(isdq:iedq,jsdq:jedq) ) ; cs%v_boundary_values(:,:) = 0.0
1562 allocate ( cs%h_boundary_values(isd:ied,jsd:jed) ) ; cs%h_boundary_values(:,:) = 0.0
1563 allocate ( cs%thickness_boundary_values(isd:ied,jsd:jed) ) ; cs%thickness_boundary_values(:,:) = 0.0
1564 allocate ( cs%ice_visc_bilinear(isd:ied,jsd:jed) ) ; cs%ice_visc_bilinear(:,:) = 0.0
1565 allocate ( cs%ice_visc_lower_tri(isd:ied,jsd:jed) ) ; cs%ice_visc_lower_tri = 0.0
1566 allocate ( cs%ice_visc_upper_tri(isd:ied,jsd:jed) ) ; cs%ice_visc_upper_tri = 0.0
1567 allocate ( cs%u_face_mask(isdq:iedq,jsd:jed) ) ; cs%u_face_mask(:,:) = 0.0
1568 allocate ( cs%v_face_mask(isd:ied,jsdq:jedq) ) ; cs%v_face_mask(:,:) = 0.0
1569 allocate ( cs%u_face_mask_boundary(isdq:iedq,jsd:jed) ) ; cs%u_face_mask_boundary(:,:) = -2.0
1570 allocate ( cs%v_face_mask_boundary(isd:ied,jsdq:jedq) ) ; cs%v_face_mask_boundary(:,:) = -2.0
1571 allocate ( cs%u_flux_boundary_values(isdq:iedq,jsd:jed) ) ; cs%u_flux_boundary_values(:,:) = 0.0
1572 allocate ( cs%v_flux_boundary_values(isd:ied,jsdq:jedq) ) ; cs%v_flux_boundary_values(:,:) = 0.0
1573 allocate ( cs%umask(isdq:iedq,jsdq:jedq) ) ; cs%umask(:,:) = -1.0
1574 allocate ( cs%vmask(isdq:iedq,jsdq:jedq) ) ; cs%vmask(:,:) = -1.0
1576 allocate ( cs%taub_beta_eff_bilinear(isd:ied,jsd:jed) ) ; cs%taub_beta_eff_bilinear(:,:) = 0.0
1577 allocate ( cs%taub_beta_eff_upper_tri(isd:ied,jsd:jed) ) ; cs%taub_beta_eff_upper_tri(:,:) = 0.0
1578 allocate ( cs%taub_beta_eff_lower_tri(isd:ied,jsd:jed) ) ; cs%taub_beta_eff_lower_tri(:,:) = 0.0
1579 allocate ( cs%OD_rt(isd:ied,jsd:jed) ) ; cs%OD_rt(:,:) = 0.0
1580 allocate ( cs%OD_av(isd:ied,jsd:jed) ) ; cs%OD_av(:,:) = 0.0
1581 allocate ( cs%float_frac(isd:ied,jsd:jed) ) ; cs%float_frac(:,:) = 0.0
1582 allocate ( cs%float_frac_rt(isd:ied,jsd:jed) ) ; cs%float_frac_rt(:,:) = 0.0
1584 if (cs%calve_to_mask)
then 1585 allocate ( cs%calve_mask (isd:ied,jsd:jed) ) ; cs%calve_mask(:,:) = 0.0
1591 if (.not. solo_ice_sheet)
then 1592 if (is_root_pe()) print *,
"initialize_ice_shelf: allocating fluxes" 1596 call allocate_forcing_type(g, fluxes, ustar=.true., shelf=.true., &
1597 press=.true., water=cs%isthermo, heat=cs%isthermo)
1599 if (is_root_pe()) print *,
"allocating fluxes in solo mode" 1600 call allocate_forcing_type(g, fluxes, ustar=.true., shelf=.true., press=.true.)
1604 call mom_initialize_topography(g%bathyT, g%max_depth, dg, param_file)
1606 call mom_initialize_rotation(g%CoriolisBu, dg, param_file)
1607 call copy_dyngrid_to_mom_grid(dg, cs%grid)
1609 call destroy_dyn_horgrid(dg)
1612 call restart_init(param_file, cs%restart_CSp,
"Shelf.res")
1613 vd = var_desc(
"shelf_mass",
"kg m-2",
"Ice shelf mass",z_grid=
'1')
1614 call register_restart_field(cs%mass_shelf, vd, .true., cs%restart_CSp)
1615 vd = var_desc(
"shelf_area",
"m2",
"Ice shelf area in cell",z_grid=
'1')
1616 call register_restart_field(cs%area_shelf_h, vd, .true., cs%restart_CSp)
1617 vd = var_desc(
"h_shelf",
"m",
"ice sheet/shelf thickness",z_grid=
'1')
1618 call register_restart_field(cs%h_shelf, vd, .true., cs%restart_CSp)
1620 if (cs%shelf_mass_is_dynamic .and. .not.cs%override_shelf_movement)
then 1622 vd = var_desc(
"u_shelf",
"m s-1",
"ice sheet/shelf velocity",
'q',z_grid=
'1')
1623 call register_restart_field(cs%u_shelf, vd, .true., cs%restart_CSp)
1624 vd = var_desc(
"v_shelf",
"m s-1",
"ice sheet/shelf velocity",
'q',z_grid=
'1')
1625 call register_restart_field(cs%v_shelf, vd, .true., cs%restart_CSp)
1629 vd = var_desc(
"h_mask",
"none",
"ice sheet/shelf thickness mask",z_grid=
'1')
1630 call register_restart_field(cs%hmask, vd, .true., cs%restart_CSp)
1633 vd = var_desc(
"t_shelf",
"deg C",
"ice sheet/shelf temperature",z_grid=
'1')
1634 call register_restart_field(cs%t_shelf, vd, .true., cs%restart_CSp)
1640 vd = var_desc(
"OD_av",
"m",
"avg ocean depth in a cell",z_grid=
'1')
1641 call register_restart_field(cs%OD_av, vd, .true., cs%restart_CSp)
1646 vd = var_desc(
"float_frac",
"m",
"degree of grounding",z_grid=
'1')
1647 call register_restart_field(cs%float_frac, vd, .true., cs%restart_CSp)
1652 vd = var_desc(
"viscosity",
"m",
"glens law ice visc",z_grid=
'1')
1653 call register_restart_field(cs%ice_visc_bilinear, vd, .true., cs%restart_CSp)
1654 vd = var_desc(
"tau_b_beta",
"m",
"coefficient of basal traction",z_grid=
'1')
1655 call register_restart_field(cs%taub_beta_eff_bilinear, vd, .true., cs%restart_CSp)
1666 cs%restart_output_dir = dirs%restart_output_dir
1669 if ((dirs%input_filename(1:1) ==
'n') .and. &
1670 (len_trim(dirs%input_filename) == 1)) new_sim = .true.
1672 if (cs%override_shelf_movement .and. cs%mass_from_file)
then 1679 call initialize_ice_thickness (cs%h_shelf, cs%area_shelf_h, cs%hmask, g, param_file)
1684 if ((cs%hmask(i,j) .eq. 1) .or. (cs%hmask(i,j) .eq. 2))
then 1685 cs%mass_shelf(i,j) = cs%h_shelf(i,j)*cs%density_ice
1690 if (cs%min_thickness_simple_calve > 0.0)
then 1703 if (cs%shelf_mass_is_dynamic .and. .not. cs%override_shelf_movement)
then 1713 if (new_sim .and. (.not. (cs%override_shelf_movement .and. cs%mass_from_file)))
then 1716 call initialize_ice_thickness (cs%h_shelf, cs%area_shelf_h, cs%hmask, g, param_file)
1721 if ((cs%hmask(i,j) .eq. 1) .or. (cs%hmask(i,j) .eq. 2))
then 1722 cs%mass_shelf(i,j) = cs%h_shelf(i,j)*cs%density_ice
1728 elseif (.not.new_sim)
then 1731 call restore_state(dirs%input_filename, dirs%restart_input_dir, time, &
1738 if (cs%shelf_mass_is_dynamic .and. .not.cs%override_shelf_movement)
then 1743 if (.not. g%symmetric)
then 1746 if (((i+g%idg_offset) .eq. (g%domain%nihalo+1)).and.(cs%u_face_mask(i-1,j).eq.3))
then 1747 cs%u_shelf (i-1,j-1) = cs%u_boundary_values (i-1,j-1)
1748 cs%u_shelf (i-1,j) = cs%u_boundary_values (i-1,j)
1750 if (((j+g%jdg_offset) .eq. (g%domain%njhalo+1)).and.(cs%v_face_mask(i,j-1).eq.3))
then 1751 cs%u_shelf (i-1,j-1) = cs%u_boundary_values (i-1,j-1)
1752 cs%u_shelf (i,j-1) = cs%u_boundary_values (i,j-1)
1758 call pass_var (cs%OD_av,g%domain)
1759 call pass_var (cs%float_frac,g%domain)
1760 call pass_var (cs%ice_visc_bilinear,g%domain)
1761 call pass_var (cs%taub_beta_eff_bilinear,g%domain)
1762 call pass_vector(cs%u_shelf, cs%v_shelf, g%domain, to_all, bgrid_ne)
1763 call pass_var (cs%area_shelf_h,g%domain)
1764 call pass_var (cs%h_shelf,g%domain)
1765 call pass_var (cs%hmask,g%domain)
1767 if (is_root_pe()) print *,
"RESTORING ICE SHELF FROM FILE!!!!!!!!!!!!!" 1774 call pass_var(cs%area_shelf_h, g%domain)
1775 call pass_var(cs%h_shelf, g%domain)
1776 call pass_var(cs%mass_shelf, g%domain)
1779 if (cs%shelf_mass_is_dynamic .and. .not.cs%override_shelf_movement)
then 1781 call pass_var(g%bathyT, g%domain)
1782 call pass_var(cs%hmask, g%domain)
1787 do j=jsd,jed ;
do i=isd,ied
1788 if (cs%area_shelf_h(i,j) > g%areaT(i,j))
then 1789 call mom_error(warning,
"Initialize_ice_shelf: area_shelf_h exceeds G%areaT.")
1790 cs%area_shelf_h(i,j) = g%areaT(i,j)
1793 if (g%areaT(i,j) > 0.0) fluxes%frac_shelf_h(i,j) = cs%area_shelf_h(i,j) / g%areaT(i,j)
1794 if (
associated(fluxes%p_surf)) &
1795 fluxes%p_surf(i,j) = fluxes%p_surf(i,j) + &
1796 fluxes%frac_shelf_h(i,j) * (cs%g_Earth * cs%mass_shelf(i,j))
1797 if (
associated(fluxes%p_surf_full)) &
1798 fluxes%p_surf_full(i,j) = fluxes%p_surf_full(i,j) + &
1799 fluxes%frac_shelf_h(i,j) * (cs%g_Earth * cs%mass_shelf(i,j))
1804 call hchksum (fluxes%frac_shelf_h,
"IS init: frac_shelf_h", g%HI, haloshift=0)
1807 if (.not. solo_ice_sheet)
then 1808 do j=jsd,jed ;
do i=isd,ied-1
1810 fluxes%frac_shelf_u(i,j) = 0.0
1811 if ((g%areaT(i,j) + g%areaT(i+1,j) > 0.0)) &
1812 fluxes%frac_shelf_u(i,j) = ((cs%area_shelf_h(i,j) + cs%area_shelf_h(i+1,j)) / &
1813 (g%areaT(i,j) + g%areaT(i+1,j)))
1814 fluxes%rigidity_ice_u(i,j) = (cs%kv_ice / cs%density_ice) * &
1815 min(cs%mass_shelf(i,j), cs%mass_shelf(i+1,j))
1819 do j=jsd,jed-1 ;
do i=isd,ied
1821 fluxes%frac_shelf_v(i,j) = 0.0
1822 if ((g%areaT(i,j) + g%areaT(i,j+1) > 0.0)) &
1823 fluxes%frac_shelf_v(i,j) = ((cs%area_shelf_h(i,j) + cs%area_shelf_h(i,j+1)) / &
1824 (g%areaT(i,j) + g%areaT(i,j+1)))
1825 fluxes%rigidity_ice_v(i,j) = (cs%kv_ice / cs%density_ice) * &
1826 min(cs%mass_shelf(i,j), cs%mass_shelf(i,j+1))
1830 if (.not. solo_ice_sheet)
then 1831 call pass_vector(fluxes%frac_shelf_u, fluxes%frac_shelf_v, g%domain, to_all, cgrid_ne)
1841 if (cs%shelf_mass_is_dynamic .and. cs%calve_to_mask .and. &
1842 .not.cs%override_shelf_movement)
then 1844 call mom_mesg(
" MOM_ice_shelf.F90, initialize_ice_shelf: reading calving_mask")
1846 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
1847 inputdir = slasher(inputdir)
1848 call get_param(param_file, mdl,
"CALVING_MASK_FILE", ic_file, &
1849 "The file with a mask for where calving might occur.", &
1850 default=
"ice_shelf_h.nc")
1851 call get_param(param_file, mdl,
"CALVING_MASK_VARNAME", var_name, &
1852 "The variable to use in masking calving.", &
1853 default=
"area_shelf_h")
1855 filename = trim(inputdir)//trim(ic_file)
1856 call log_param(param_file, mdl,
"INPUTDIR/CALVING_MASK_FILE", filename)
1857 if (.not.file_exists(filename, g%Domain))
call mom_error(fatal, &
1858 " calving mask file: Unable to open "//trim(filename))
1860 call read_data(filename,trim(var_name),cs%calve_mask,domain=g%Domain%mpp_domain)
1863 if (cs%calve_mask(i,j) > 0.0) cs%calve_mask(i,j) = 1.0
1867 call pass_var (cs%calve_mask,g%domain)
1870 if (cs%shelf_mass_is_dynamic .and. .not.cs%override_shelf_movement)
then 1873 if (.not. cs%isthermo)
then 1879 if (is_root_pe()) print *,
"NEW SIM: initialize velocity" 1885 if (cs%id_u_shelf > 0)
call post_data(cs%id_u_shelf,cs%u_shelf,cs%diag)
1886 if (cs%id_v_shelf > 0)
call post_data(cs%id_v_shelf,cs%v_shelf,cs%diag)
1890 call get_param(param_file, mdl,
"SAVE_INITIAL_CONDS", save_ic, &
1891 "If true, save the ice shelf initial conditions.", &
1893 if (save_ic)
call get_param(param_file, mdl,
"SHELF_IC_OUTPUT_FILE", ic_file,&
1894 "The name-root of the output file for the ice shelf \n"//&
1895 "initial conditions.", default=
"MOM_Shelf_IC")
1897 if (save_ic .and. .not.((dirs%input_filename(1:1) ==
'r') .and. &
1898 (len_trim(dirs%input_filename) == 1)))
then 1900 call save_restart(dirs%output_directory, cs%Time, g, &
1901 cs%restart_CSp, filename=ic_file)
1905 cs%id_area_shelf_h = register_diag_field(
'ocean_model',
'area_shelf_h', cs%diag%axesT1, cs%Time, &
1906 'Ice Shelf Area in cell',
'meter-2')
1907 cs%id_shelf_mass = register_diag_field(
'ocean_model',
'shelf_mass', cs%diag%axesT1, cs%Time, &
1908 'mass of shelf',
'kg/m^2')
1909 cs%id_mass_flux = register_diag_field(
'ocean_model',
'mass_flux', cs%diag%axesT1,&
1910 cs%Time,
'Total mass flux of freshwater across the ice-ocean interface.',
'kg/s')
1911 cs%id_melt = register_diag_field(
'ocean_model',
'melt', cs%diag%axesT1, cs%Time, &
1912 'Ice Shelf Melt Rate',
'meter year-1')
1913 cs%id_thermal_driving = register_diag_field(
'ocean_model',
'thermal_driving', cs%diag%axesT1, cs%Time, &
1914 'pot. temp. in the boundary layer minus freezing pot. temp. at the ice-ocean interface.',
'Celsius')
1915 cs%id_haline_driving = register_diag_field(
'ocean_model',
'haline_driving', cs%diag%axesT1, cs%Time, &
1916 'salinity in the boundary layer minus salinity at the ice-ocean interface.',
'PPT')
1917 cs%id_Sbdry = register_diag_field(
'ocean_model',
'sbdry', cs%diag%axesT1, cs%Time, &
1918 'salinity at the ice-ocean interface.',
'PPT')
1919 cs%id_u_ml = register_diag_field(
'ocean_model',
'u_ml', cs%diag%axesT1, cs%Time, &
1920 'Eastward vel. in the boundary layer (used to compute ustar)',
'meter second-1')
1921 cs%id_v_ml = register_diag_field(
'ocean_model',
'v_ml', cs%diag%axesT1, cs%Time, &
1922 'Northward vel. in the boundary layer (used to compute ustar)',
'meter second-1')
1923 cs%id_exch_vel_s = register_diag_field(
'ocean_model',
'exch_vel_s', cs%diag%axesT1, cs%Time, &
1924 'Sub-shelf salinity exchange velocity',
'meter second-1')
1925 cs%id_exch_vel_t = register_diag_field(
'ocean_model',
'exch_vel_t', cs%diag%axesT1, cs%Time, &
1926 'Sub-shelf thermal exchange velocity',
'meter second-1')
1927 cs%id_tfreeze = register_diag_field(
'ocean_model',
'tfreeze', cs%diag%axesT1, cs%Time, &
1928 'In Situ Freezing point at ice shelf interface',
'degrees Celsius')
1929 cs%id_tfl_shelf = register_diag_field(
'ocean_model',
'tflux_shelf', cs%diag%axesT1, cs%Time, &
1930 'Heat conduction into ice shelf',
'Watts meter-2')
1931 cs%id_ustar_shelf = register_diag_field(
'ocean_model',
'ustar_shelf', cs%diag%axesT1, cs%Time, &
1932 'Fric vel under shelf',
'm/s')
1934 if (cs%shelf_mass_is_dynamic .and. .not.cs%override_shelf_movement)
then 1935 cs%id_u_shelf = register_diag_field(
'ocean_model',
'u_shelf',cs%diag%axesB1,cs%Time, &
1936 'x-velocity of ice',
'm year')
1937 cs%id_v_shelf = register_diag_field(
'ocean_model',
'v_shelf',cs%diag%axesB1,cs%Time, &
1938 'y-velocity of ice',
'm year')
1939 cs%id_u_mask = register_diag_field(
'ocean_model',
'u_mask',cs%diag%axesB1,cs%Time, &
1940 'mask for u-nodes',
'none')
1941 cs%id_v_mask = register_diag_field(
'ocean_model',
'v_mask',cs%diag%axesB1,cs%Time, &
1942 'mask for v-nodes',
'none')
1943 cs%id_h_mask = register_diag_field(
'ocean_model',
'h_mask',cs%diag%axesT1,cs%Time, &
1944 'ice shelf thickness',
'none')
1945 cs%id_surf_elev = register_diag_field(
'ocean_model',
'ice_surf',cs%diag%axesT1,cs%Time, &
1946 'ice surf elev',
'm')
1947 cs%id_float_frac = register_diag_field(
'ocean_model',
'ice_float_frac',cs%diag%axesT1,cs%Time, &
1948 'fraction of cell that is floating (sort of)',
'none')
1949 cs%id_col_thick = register_diag_field(
'ocean_model',
'col_thick',cs%diag%axesT1,cs%Time, &
1950 'ocean column thickness passed to ice model',
'm')
1951 cs%id_OD_av = register_diag_field(
'ocean_model',
'OD_av',cs%diag%axesT1,cs%Time, &
1952 'intermediate ocean column thickness passed to ice model',
'm')
1953 cs%id_float_frac_rt = register_diag_field(
'ocean_model',
'float_frac_rt',cs%diag%axesT1,cs%Time, &
1954 'timesteps where cell is floating ',
'none')
1963 cs%id_t_shelf = register_diag_field(
'ocean_model',
't_shelf',cs%diag%axesT1,cs%Time, &
1965 cs%id_t_mask = register_diag_field(
'ocean_model',
'tmask',cs%diag%axesT1,cs%Time, &
1966 'mask for T-nodes',
'none')
1969 id_clock_shelf = cpu_clock_id(
'Ice shelf', grain=clock_component)
1970 id_clock_pass = cpu_clock_id(
' Ice shelf halo updates', grain=clock_routine)
1977 type(ocean_grid_type),
intent(in) :: G
1978 type(param_file_type),
intent(in) :: param_file
1980 logical,
optional :: new_sim
1982 integer :: i, j, is, ie, js, je
1983 logical :: read_shelf_area, new_sim_2
1984 character(len=240) :: config, inputdir, shelf_file, filename
1985 character(len=120) :: shelf_mass_var
1986 character(len=120) :: shelf_area_var
1987 character(len=40) :: mdl =
"MOM_ice_shelf" 1988 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1990 if (.not.
present(new_sim))
then 1996 call get_param(param_file, mdl,
"ICE_SHELF_CONFIG", config, &
1997 "A string that specifies how the ice shelf is \n"//&
1998 "initialized. Valid options include:\n"//&
1999 " \tfile\t Read from a file.\n"//&
2000 " \tzero\t Set shelf mass to 0 everywhere.\n"//&
2001 " \tUSER\t Call USER_initialize_shelf_mass.\n", &
2002 fail_if_missing=.true.)
2004 select case ( trim(config) )
2007 call time_interp_external_init()
2009 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
2010 inputdir = slasher(inputdir)
2012 call get_param(param_file, mdl,
"SHELF_FILE", shelf_file, &
2013 "If DYNAMIC_SHELF_MASS = True, OVERRIDE_SHELF_MOVEMENT = True \n"//&
2014 "and ICE_SHELF_MASS_FROM_FILE = True, this is the file from \n"//&
2015 "which to read the shelf mass and area.", &
2016 default=
"shelf_mass.nc")
2017 call get_param(param_file, mdl,
"SHELF_MASS_VAR", shelf_mass_var, &
2018 "The variable in SHELF_FILE with the shelf mass.", &
2019 default=
"shelf_mass")
2020 call get_param(param_file, mdl,
"READ_SHELF_AREA", read_shelf_area, &
2021 "If true, also read the area covered by ice-shelf from SHELF_FILE.", &
2024 filename = trim(slasher(inputdir))//trim(shelf_file)
2025 call log_param(param_file, mdl,
"INPUTDIR/SHELF_FILE", filename)
2028 cs%id_read_mass = init_external_field(filename,shelf_mass_var, &
2029 domain=g%Domain%mpp_domain,verbose=.true.)
2031 cs%id_read_mass = init_external_field(filename,shelf_mass_var, &
2032 domain=g%Domain%mpp_domain)
2036 if (read_shelf_area)
then 2037 call get_param(param_file, mdl,
"SHELF_AREA_VAR", shelf_area_var, &
2038 "The variable in SHELF_FILE with the shelf area.", &
2039 default=
"shelf_area")
2041 cs%id_read_area = init_external_field(filename,shelf_area_var, &
2042 domain=g%Domain%mpp_domain)
2045 if (.not.file_exists(filename, g%Domain))
call mom_error(fatal, &
2046 " initialize_shelf_mass: Unable to open "//trim(filename))
2049 do j=js,je ;
do i=is,ie
2050 cs%mass_shelf(i,j) = 0.0
2051 cs%area_shelf_h(i,j) = 0.0
2055 call user_initialize_shelf_mass(cs%mass_shelf, cs%area_shelf_h, &
2056 cs%h_shelf, cs%hmask, g, cs%user_CS, param_file, new_sim_2)
2058 case default ;
call mom_error(fatal,
"initialize_ice_shelf: "// &
2059 "Unrecognized ice shelf setup "//trim(config))
2066 type(ocean_grid_type),
intent(inout) :: G
2068 type(time_type),
intent(in) :: Time
2069 type(forcing),
intent(inout) :: fluxes
2072 integer :: i, j, is, ie, js, je
2073 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
2076 do j=js,je;
do i=is,ie
2081 call time_interp_external(cs%id_read_mass, time, cs%mass_shelf)
2083 do j=js,je ;
do i=is,ie
2085 if (cs%area_shelf_h(i,j) > 0.0)
then 2086 if (
associated(fluxes%lprec)) fluxes%lprec(i,j) = 0.0
2087 if (
associated(fluxes%sens)) fluxes%sens(i,j) = 0.0
2088 if (
associated(fluxes%frac_shelf_h)) fluxes%frac_shelf_h(i,j) = 0.0
2089 if (
associated(fluxes%p_surf)) fluxes%p_surf(i,j) = 0.0
2090 if (
associated(fluxes%salt_flux)) fluxes%salt_flux(i,j) = 0.0
2092 cs%area_shelf_h(i,j) = 0.0
2094 if (cs%mass_shelf(i,j) > 0.0)
then 2095 cs%area_shelf_h(i,j) = g%areaT(i,j)
2096 cs%h_shelf(i,j) = cs%mass_shelf(i,j)/cs%density_ice
2104 if (cs%min_thickness_simple_calve > 0.0)
then 2108 call pass_var(cs%area_shelf_h, g%domain)
2109 call pass_var(cs%h_shelf, g%domain)
2110 call pass_var(cs%hmask, g%domain)
2111 call pass_var(cs%mass_shelf, g%domain)
2115 do j=js,je ;
do i=is,ie
2116 if (
associated(fluxes%p_surf)) &
2117 fluxes%p_surf(i,j) = (cs%g_Earth * cs%mass_shelf(i,j))
2118 if (
associated(fluxes%p_surf_full)) &
2119 fluxes%p_surf_full(i,j) = (cs%g_Earth * cs%mass_shelf(i,j))
2120 if (g%areaT(i,j) > 0.0) &
2121 fluxes%frac_shelf_h(i,j) = cs%area_shelf_h(i,j) / g%areaT(i,j)
2130 type(time_type),
intent(in) :: Time
2132 type(ocean_grid_type),
pointer :: G
2133 integer :: i, j, iters, isd, ied, jsd, jed
2134 real :: rhoi, rhow, OD
2135 type(time_type) :: dummy_time
2136 real,
dimension(:,:),
pointer :: OD_av, float_frac, h_shelf
2139 rhoi = cs%density_ice
2140 rhow = cs%density_ocean_avg
2141 dummy_time = set_time(0,0)
2143 h_shelf => cs%h_shelf
2144 float_frac => cs%float_frac
2145 isd=g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
2149 od = g%bathyT(i,j) - rhoi/rhow * h_shelf(i,j)
2153 float_frac(i,j) = 0.
2156 float_frac(i,j) = 1.
2168 type(time_type),
intent(in) :: Time
2169 character(len=*),
optional,
intent(in) :: directory
2171 logical,
optional,
intent(in) :: time_stamped
2173 character(len=*),
optional,
intent(in) :: filename_suffix
2176 type(ocean_grid_type),
pointer :: G
2177 character(len=200) :: restart_dir
2178 character(2) :: procnum
2193 if (
present(directory))
then ; restart_dir = directory
2194 else ; restart_dir = cs%restart_output_dir ;
endif 2196 call save_restart(restart_dir, time, cs%grid, cs%restart_CSp, time_stamped)
2203 real,
intent(in) :: time_step
2204 real,
pointer,
dimension(:,:),
intent(in) :: melt_rate
2205 type(time_type) :: Time
2246 type(ocean_grid_type),
pointer :: G
2247 real,
dimension(size(CS%h_shelf,1),size(CS%h_shelf,2)) :: h_after_uflux, h_after_vflux
2248 real,
dimension(size(CS%h_shelf,1),size(CS%h_shelf,2),4) :: flux_enter
2249 integer :: isd, ied, jsd, jed, i, j, isc, iec, jsc, jec
2250 real :: rho, spy, thick_bd
2251 real,
dimension(:,:),
pointer :: hmask
2252 character(len=2) :: procnum
2256 rho = cs%density_ice
2259 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
2260 isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
2261 flux_enter(:,:,:) = 0.0
2263 h_after_uflux(:,:) = 0.0
2264 h_after_vflux(:,:) = 0.0
2269 thick_bd = cs%thickness_boundary_values(i,j)
2270 if (thick_bd .ne. 0.0)
then 2271 cs%h_shelf(i,j) = cs%thickness_boundary_values(i,j)
2292 if (cs%hmask(i,j) .eq. 1)
then 2293 cs%h_shelf (i,j) = h_after_vflux(i,j)
2298 if (cs%moving_shelf_front)
then 2300 if (cs%min_thickness_simple_calve > 0.0)
then 2303 if (cs%calve_to_mask)
then 2304 call calve_to_mask (cs, cs%h_shelf, cs%area_shelf_h, cs%hmask, cs%calve_mask)
2320 real,
dimension(NILIMB_SYM_,NJLIMB_SYM_),
intent(inout) :: u, v
2321 integer,
intent(in) :: FE
2322 integer,
intent(out) :: iters
2323 type(time_type),
intent(in) :: time
2325 real,
dimension(:,:),
pointer :: TAUDX, TAUDY, u_prev_iterate, v_prev_iterate, &
2326 u_bdry_cont, v_bdry_cont, Au, Av, err_u, err_v, &
2327 geolonq, geolatq, u_last, v_last, float_cond, H_node
2328 type(ocean_grid_type),
pointer :: G
2329 integer :: conv_flag, i, j, k,l, iter, isym, &
2330 isdq, iedq, jsdq, jedq, isd, ied, jsd, jed, isumstart, jsumstart, nodefloat, nsub
2331 real :: err_max, err_tempu, err_tempv, err_init, area, max_vel, tempu, tempv, rhoi, rhow
2332 real,
pointer,
dimension(:,:,:,:) :: Phi
2333 real,
pointer,
dimension(:,:,:,:,:,:) :: Phisub
2334 real,
dimension (8,4) :: Phi_temp
2335 real,
dimension (2,2) :: X,Y
2336 character(2) :: iternum
2337 character(2) :: procnum, numproc
2340 nsub = cs%n_sub_regularize
2343 isdq = g%isdB ; iedq = g%iedB ; jsdq = g%jsdB ; jedq = g%jedB
2344 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
2345 rhoi = cs%density_ice
2346 rhow = cs%density_ocean_avg
2347 ALLOCATE (taudx(isdq:iedq,jsdq:jedq) ) ; taudx(:,:)=0
2348 ALLOCATE (taudy(isdq:iedq,jsdq:jedq) ) ; taudy(:,:)=0
2349 ALLOCATE (u_prev_iterate(isdq:iedq,jsdq:jedq) )
2350 ALLOCATE (v_prev_iterate(isdq:iedq,jsdq:jedq) )
2351 ALLOCATE (u_bdry_cont(isdq:iedq,jsdq:jedq) ) ; u_bdry_cont(:,:)=0
2352 ALLOCATE (v_bdry_cont(isdq:iedq,jsdq:jedq) ) ; v_bdry_cont(:,:)=0
2353 ALLOCATE (au(isdq:iedq,jsdq:jedq) ) ; au(:,:)=0
2354 ALLOCATE (av(isdq:iedq,jsdq:jedq) ) ; av(:,:)=0
2355 ALLOCATE (err_u(isdq:iedq,jsdq:jedq) )
2356 ALLOCATE (err_v(isdq:iedq,jsdq:jedq) )
2357 ALLOCATE (u_last(isdq:iedq,jsdq:jedq) )
2358 ALLOCATE (v_last(isdq:iedq,jsdq:jedq) )
2361 ALLOCATE (float_cond(g%isd:g%ied,g%jsd:g%jed)) ; float_cond(:,:)=0
2362 ALLOCATE (h_node(g%isdB:g%iedB,g%jsdB:g%jedB)) ; h_node(:,:)=0
2363 ALLOCATE (phisub(nsub,nsub,2,2,2,2)) ; phisub = 0.0
2365 geolonq => g%geoLonBu ; geolatq => g%geoLatBu
2367 if (g%isc+g%idg_offset==g%isg)
then 2372 isumstart = isumstart_int_
2375 if (g%jsc+g%jdg_offset==g%jsg)
then 2380 jsumstart = jsumstart_int_
2393 if (cs%GL_regularize)
then 2396 call savearray2 (
"H_node",h_node,cs%write_output_to_file)
2403 if ((cs%hmask(i,j) .eq. 1) .and. &
2404 (rhoi/rhow * h_node(i-1+k,j-1+l) - g%bathyT(i,j) .le. 0))
then 2405 nodefloat = nodefloat + 1
2409 if ((nodefloat .gt. 0) .and. (nodefloat .lt. 4))
then 2411 float_cond(i,j) = 1.0
2412 cs%float_frac (i,j) = 1.0
2416 call savearray2 (
"float_cond",float_cond,cs%write_output_to_file)
2418 call pass_var (float_cond, g%Domain)
2422 call savearray2(
"Phisub1111",phisub(:,:,1,1,1,1),cs%write_output_to_file)
2428 u_prev_iterate(:,:) = u(:,:)
2429 v_prev_iterate(:,:) = v(:,:)
2435 allocate (phi(isd:ied,jsd:jed,1:8,1:4)) ; phi(:,:,:,:)=0
2440 if (((i .gt. isd) .and. (j .gt. jsd)) .or. (isym .eq. 1))
then 2441 x(:,:) = geolonq(i-1:i,j-1:j)*1000
2442 y(:,:) = geolatq(i-1:i,j-1:j)*1000
2444 x(2,:) = geolonq(i,j)*1000
2445 x(1,:) = geolonq(i,j)*1000-g%dxT(i,j)
2446 y(:,2) = geolatq(i,j)*1000
2447 y(:,1) = geolatq(i,j)*1000-g%dyT(i,j)
2451 phi(i,j,:,:) = phi_temp
2460 call pass_var (cs%ice_visc_bilinear, g%domain)
2461 call pass_var (cs%taub_beta_eff_bilinear, g%domain)
2465 call pass_var (cs%ice_visc_upper_tri, g%domain)
2466 call pass_var (cs%taub_beta_eff_upper_tri, g%domain)
2467 call pass_var (cs%ice_visc_lower_tri, g%domain)
2468 call pass_var (cs%taub_beta_eff_lower_tri, g%domain)
2476 cs%taub_beta_eff_bilinear (i,j) = cs%taub_beta_eff_bilinear (i,j) * cs%float_frac (i,j)
2478 cs%taub_beta_eff_upper_tri (i,j) = cs%taub_beta_eff_upper_tri (i,j) * cs%float_frac (i,j)
2479 cs%taub_beta_eff_lower_tri (i,j) = cs%taub_beta_eff_lower_tri (i,j) * cs%float_frac (i,j)
2486 rhoi/rhow, u_bdry_cont, v_bdry_cont)
2487 elseif (fe .eq. 2)
then 2491 au(:,:) = 0.0 ; av(:,:) = 0.0
2494 call cg_action_bilinear (au, av, u, v, phi, phisub, cs%umask, cs%vmask, cs%hmask, h_node, &
2495 cs%ice_visc_bilinear, float_cond, g%bathyT, cs%taub_beta_eff_bilinear, g%areaT, &
2496 g%isc-1, g%iec+1, g%jsc-1, g%jec+1, rhoi/rhow)
2497 elseif (fe .eq. 2)
then 2498 call cg_action_triangular (au, av, u, v, cs%umask, cs%vmask, cs%hmask, cs%ice_visc_upper_tri, &
2499 cs%ice_visc_lower_tri, cs%taub_beta_eff_upper_tri, cs%taub_beta_eff_lower_tri, &
2500 g%dxT, g%dyT, g%areaT, g%isc-1, g%iec+1, g%jsc-1, g%jec+1, isym)
2506 err_init = 0 ; err_tempu = 0; err_tempv = 0
2507 do j=jsumstart,g%jecB
2508 do i=isumstart,g%iecB
2509 if (cs%umask(i,j) .eq. 1)
then 2510 err_tempu = abs(au(i,j) + u_bdry_cont(i,j) - taudx(i,j))
2512 if (cs%vmask(i,j) .eq. 1)
then 2513 err_tempv = max(abs(av(i,j) + v_bdry_cont(i,j) - taudy(i,j)), err_tempu)
2515 if (err_tempv .ge. err_init)
then 2516 err_init = err_tempv
2521 call mpp_max (err_init)
2523 if (is_root_pe()) print *,
"INITIAL nonlinear residual: ",err_init
2525 u_last(:,:) = u(:,:) ; v_last(:,:) = v(:,:)
2533 fe, conv_flag, iters, time, phi, phisub)
2537 call qchksum (u,
"u shelf", g%HI, haloshift=2)
2538 call qchksum (v,
"v shelf", g%HI, haloshift=2)
2541 if (is_root_pe()) print *,
"linear solve done",iters,
" iterations" 2545 call pass_var (cs%ice_visc_bilinear, g%domain)
2546 call pass_var (cs%taub_beta_eff_bilinear, g%domain)
2549 call pass_var (cs%ice_visc_upper_tri, g%domain)
2550 call pass_var (cs%taub_beta_eff_upper_tri, g%domain)
2551 call pass_var (cs%ice_visc_lower_tri, g%domain)
2552 call pass_var (cs%taub_beta_eff_lower_tri, g%domain)
2555 if (iter .eq. 1)
then 2564 cs%taub_beta_eff_bilinear (i,j) = cs%taub_beta_eff_bilinear (i,j) * cs%float_frac (i,j)
2566 cs%taub_beta_eff_upper_tri (i,j) = cs%taub_beta_eff_upper_tri (i,j) * cs%float_frac (i,j)
2567 cs%taub_beta_eff_lower_tri (i,j) = cs%taub_beta_eff_lower_tri (i,j) * cs%float_frac (i,j)
2572 u_bdry_cont(:,:) = 0 ; v_bdry_cont(:,:) = 0
2576 rhoi/rhow, u_bdry_cont, v_bdry_cont)
2577 elseif (fe .eq. 2)
then 2581 au(:,:) = 0 ; av(:,:) = 0
2584 call cg_action_bilinear (au, av, u, v, phi, phisub, cs%umask, cs%vmask, cs%hmask, h_node, &
2585 cs%ice_visc_bilinear, float_cond, g%bathyT, cs%taub_beta_eff_bilinear, g%areaT, g%isc-1, &
2586 g%iec+1, g%jsc-1, g%jec+1, rhoi/rhow)
2587 elseif (fe .eq. 2)
then 2588 call cg_action_triangular (au, av, u, v, cs%umask, cs%vmask, cs%hmask, cs%ice_visc_upper_tri, &
2589 cs%ice_visc_lower_tri, cs%taub_beta_eff_upper_tri, cs%taub_beta_eff_lower_tri, &
2590 g%dxT, g%dyT, g%areaT, g%isc-1, g%iec+1, g%jsc-1, g%jec+1, isym)
2595 if (cs%nonlin_solve_err_mode .eq. 1)
then 2597 do j=jsumstart,g%jecB
2598 do i=isumstart,g%iecB
2599 if (cs%umask(i,j) .eq. 1)
then 2600 err_tempu = abs(au(i,j) + u_bdry_cont(i,j) - taudx(i,j))
2602 if (cs%vmask(i,j) .eq. 1)
then 2603 err_tempv = max(abs(av(i,j) + v_bdry_cont(i,j) - taudy(i,j)), err_tempu)
2605 if (err_tempv .ge. err_max)
then 2611 call mpp_max (err_max)
2613 elseif (cs%nonlin_solve_err_mode .eq. 2)
then 2615 max_vel = 0 ; tempu = 0 ; tempv = 0
2617 do j=jsumstart,g%jecB
2618 do i=isumstart,g%iecB
2619 if (cs%umask(i,j) .eq. 1)
then 2620 err_tempu = abs(u_last(i,j)-u(i,j))
2623 if (cs%vmask(i,j) .eq. 1)
then 2624 err_tempv = max(abs(v_last(i,j)- v(i,j)), err_tempu)
2625 tempv = sqrt(v(i,j)**2+tempu**2)
2627 if (err_tempv .ge. err_max)
then 2630 if (tempv .ge. max_vel)
then 2636 u_last(:,:) = u(:,:)
2637 v_last(:,:) = v(:,:)
2639 call mpp_max (max_vel)
2640 call mpp_max (err_max)
2645 if (is_root_pe()) print *,
"nonlinear residual: ",err_max/err_init
2647 if (err_max .le. cs%nonlinear_tolerance * err_init)
then 2649 print *,
"exiting nonlinear solve after ",iter,
" iterations" 2660 DEALLOCATE (u_prev_iterate)
2661 DEALLOCATE (v_prev_iterate)
2662 DEALLOCATE (u_bdry_cont)
2663 DEALLOCATE (v_bdry_cont)
2671 DEALLOCATE (float_cond)
2676 subroutine ice_shelf_solve_inner (CS, u, v, taudx, taudy, H_node, float_cond, FE, conv_flag, iters, time, Phi, Phisub)
2678 real,
dimension(NILIMB_SYM_,NJLIMB_SYM_),
intent(inout) :: u, v
2679 real,
dimension(NILIMB_SYM_,NJLIMB_SYM_),
intent(in) :: taudx, taudy, H_node
2680 real,
dimension(:,:),
intent(in) :: float_cond
2681 integer,
intent(in) :: FE
2682 integer,
intent(out) :: conv_flag, iters
2683 type(time_type) :: time
2684 real,
pointer,
dimension(:,:,:,:) :: Phi
2685 real,
dimension (:,:,:,:,:,:),
pointer :: Phisub
2697 real,
dimension(:,:),
pointer :: hmask, umask, vmask, u_bdry, v_bdry, &
2698 visc, visc_lo, beta, beta_lo, geolonq, geolatq
2699 real,
dimension(LBOUND(u,1):UBOUND(u,1),LBOUND(u,2):UBOUND(u,2)) :: &
2700 Ru, Rv, Zu, Zv, DIAGu, DIAGv, RHSu, RHSv, &
2701 ubd, vbd, Au, Av, Du, Dv, &
2702 Zu_old, Zv_old, Ru_old, Rv_old, &
2704 integer :: iter, i, j, isym, isd, ied, jsd, jed, &
2705 isc, iec, jsc, jec, is, js, ie, je, isumstart, jsumstart, &
2706 isdq, iedq, jsdq, jedq, iscq, iecq, jscq, jecq, nx_halo, ny_halo
2707 real :: tol, beta_k, alpha_k, area, dot_p1, dot_p2, resid0, cg_halo, dot_p1a, dot_p2a
2708 type(ocean_grid_type),
pointer :: G
2709 character(1) :: procnum
2710 character(2) :: gridsize
2712 real,
dimension (8,4) :: Phi_temp
2713 real,
dimension (2,2) :: X,Y
2718 u_bdry => cs%u_boundary_values
2719 v_bdry => cs%v_boundary_values
2722 geolonq => g%geoLonBu
2723 geolatq => g%geoLatBu
2725 isdq = g%isdB ; iedq = g%iedB ; jsdq = g%jsdB ; jedq = g%jedB
2726 iscq = g%iscB ; iecq = g%iecB ; jscq = g%jscB ; jecq = g%jecB
2727 ny_halo = g%domain%njhalo ; nx_halo = g%domain%nihalo
2728 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
2729 isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
2731 zu(:,:) = 0 ; zv(:,:) = 0 ; diagu(:,:) = 0 ; diagv(:,:) = 0
2732 ru(:,:) = 0 ; rv(:,:) = 0 ; au(:,:) = 0 ; av(:,:) = 0
2733 du(:,:) = 0 ; dv(:,:) = 0 ; ubd(:,:) = 0 ; vbd(:,:) = 0
2734 dot_p1 = 0 ; dot_p2 = 0
2744 if (g%isc+g%idg_offset==g%isg)
then 2749 isumstart = isumstart_int_
2752 if (g%jsc+g%jdg_offset==g%jsg)
then 2757 jsumstart = jsumstart_int_
2761 visc => cs%ice_visc_bilinear
2762 beta => cs%taub_beta_eff_bilinear
2763 elseif (fe .eq. 2)
then 2764 visc => cs%ice_visc_upper_tri
2765 visc_lo => cs%ice_visc_lower_tri
2766 beta => cs%taub_beta_eff_upper_tri
2767 beta_lo => cs%taub_beta_eff_lower_tri
2772 cs%density_ice/cs%density_ocean_avg, ubd, vbd)
2773 elseif (fe .eq. 2)
then 2777 rhsu(:,:) = taudx(:,:) - ubd(:,:)
2778 rhsv(:,:) = taudy(:,:) - vbd(:,:)
2781 call pass_vector(rhsu, rhsv, g%domain, to_all, bgrid_ne)
2786 cs%density_ice/cs%density_ocean_avg, phisub, diagu, diagv)
2788 elseif (fe .eq. 2)
then 2790 diagu(:,:) = 1 ; diagv(:,:) = 1
2793 call pass_vector(diagu, diagv, g%domain, to_all, bgrid_ne)
2799 h_node, visc, float_cond, g%bathyT, beta, g%areaT, isc-1, iec+1, jsc-1, &
2800 jec+1, cs%density_ice/cs%density_ocean_avg)
2801 elseif (fe .eq. 2)
then 2803 beta, beta_lo, g%dxT, g%dyT, g%areaT, isc-1, iec+1, jsc-1, jec+1, isym)
2806 call pass_vector(au, av, g%domain, to_all, bgrid_ne)
2808 ru(:,:) = rhsu(:,:) - au(:,:) ; rv(:,:) = rhsv(:,:) - av(:,:)
2810 if (.not. cs%use_reproducing_sums)
then 2814 if (umask(i,j) .eq. 1) dot_p1 = dot_p1 + ru(i,j)**2
2815 if (vmask(i,j) .eq. 1) dot_p1 = dot_p1 + rv(i,j)**2
2819 call mpp_sum (dot_p1)
2825 do j=jsumstart_int_,jecq
2826 do i=isumstart_int_,iecq
2827 if (umask(i,j) .eq. 1) sum_vec(i,j) = ru(i,j)**2
2828 if (vmask(i,j) .eq. 1) sum_vec(i,j) = sum_vec(i,j) + rv(i,j)**2
2832 dot_p1 = reproducing_sum( sum_vec, isumstart_int_, iecq, &
2833 jsumstart_int_, jecq )
2837 resid0 = sqrt(dot_p1)
2841 if (umask(i,j) .eq. 1) zu(i,j) = ru(i,j) / diagu(i,j)
2842 if (vmask(i,j) .eq. 1) zv(i,j) = rv(i,j) / diagv(i,j)
2846 du(:,:) = zu(:,:) ; dv(:,:) = zv(:,:)
2861 do iter = 1,cs%cg_max_iterations
2868 is = isc - cg_halo ; ie = iecq + cg_halo
2869 js = jscq - cg_halo ; je = jecq + cg_halo
2871 au(:,:) = 0 ; av(:,:) = 0
2876 h_node, visc, float_cond, g%bathyT, beta, g%areaT, is, ie, js, &
2877 je, cs%density_ice/cs%density_ocean_avg)
2879 elseif (fe .eq. 2)
then 2882 beta, beta_lo, g%dxT, g%dyT, g%areaT, is, ie, js, je, isym)
2888 if ( .not. cs%use_reproducing_sums)
then 2892 dot_p1 = 0 ; dot_p2 = 0
2895 if (umask(i,j) .eq. 1)
then 2896 dot_p1 = dot_p1 + zu(i,j)*ru(i,j)
2897 dot_p2 = dot_p2 + du(i,j)*au(i,j)
2899 if (vmask(i,j) .eq. 1)
then 2900 dot_p1 = dot_p1 + zv(i,j)*rv(i,j)
2901 dot_p2 = dot_p2 + dv(i,j)*av(i,j)
2905 call mpp_sum (dot_p1) ;
call mpp_sum (dot_p2)
2908 sum_vec(:,:) = 0.0 ; sum_vec_2(:,:) = 0.0
2912 if (umask(i,j) .eq. 1) sum_vec(i,j) = zu(i,j) * ru(i,j)
2913 if (vmask(i,j) .eq. 1) sum_vec(i,j) = sum_vec(i,j) + &
2916 if (umask(i,j) .eq. 1) sum_vec_2(i,j) = du(i,j) * au(i,j)
2917 if (vmask(i,j) .eq. 1) sum_vec_2(i,j) = sum_vec_2(i,j) + &
2922 dot_p1 = reproducing_sum( sum_vec, iscq, iecq, &
2925 dot_p2 = reproducing_sum( sum_vec_2, iscq, iecq, &
2930 alpha_k = dot_p1/dot_p2
2940 if (umask(i,j) .eq. 1) u(i,j) = u(i,j) + alpha_k * du(i,j)
2941 if (vmask(i,j) .eq. 1) v(i,j) = v(i,j) + alpha_k * dv(i,j)
2947 if (umask(i,j) .eq. 1)
then 2948 ru_old(i,j) = ru(i,j) ; zu_old(i,j) = zu(i,j)
2950 if (vmask(i,j) .eq. 1)
then 2951 rv_old(i,j) = rv(i,j) ; zv_old(i,j) = zv(i,j)
2961 if (umask(i,j) .eq. 1) ru(i,j) = ru(i,j) - alpha_k * au(i,j)
2962 if (vmask(i,j) .eq. 1) rv(i,j) = rv(i,j) - alpha_k * av(i,j)
2969 if (umask(i,j) .eq. 1)
then 2970 zu(i,j) = ru(i,j) / diagu(i,j)
2972 if (vmask(i,j) .eq. 1)
then 2973 zv(i,j) = rv(i,j) / diagv(i,j)
2980 if (.not. cs%use_reproducing_sums)
then 2983 dot_p1 = 0 ; dot_p2 = 0
2986 if (umask(i,j) .eq. 1)
then 2987 dot_p1 = dot_p1 + zu(i,j)*ru(i,j)
2988 dot_p2 = dot_p2 + zu_old(i,j)*ru_old(i,j)
2990 if (vmask(i,j) .eq. 1)
then 2991 dot_p1 = dot_p1 + zv(i,j)*rv(i,j)
2992 dot_p2 = dot_p2 + zv_old(i,j)*rv_old(i,j)
2996 call mpp_sum (dot_p1) ;
call mpp_sum (dot_p2)
3001 sum_vec(:,:) = 0.0 ; sum_vec_2(:,:) = 0.0
3003 do j=jsumstart_int_,jecq
3004 do i=isumstart_int_,iecq
3005 if (umask(i,j) .eq. 1) sum_vec(i,j) = zu(i,j) * ru(i,j)
3006 if (vmask(i,j) .eq. 1) sum_vec(i,j) = sum_vec(i,j) + &
3009 if (umask(i,j) .eq. 1) sum_vec_2(i,j) = zu_old(i,j) * ru_old(i,j)
3010 if (vmask(i,j) .eq. 1) sum_vec_2(i,j) = sum_vec_2(i,j) + &
3011 zv_old(i,j) * rv_old(i,j)
3016 dot_p1 = reproducing_sum( sum_vec, isumstart_int_, iecq, &
3017 jsumstart_int_, jecq )
3019 dot_p2 = reproducing_sum( sum_vec_2, isumstart_int_, iecq, &
3020 jsumstart_int_, jecq )
3024 beta_k = dot_p1/dot_p2
3032 if (umask(i,j) .eq. 1) du(i,j) = zu(i,j) + beta_k * du(i,j)
3033 if (vmask(i,j) .eq. 1) dv(i,j) = zv(i,j) + beta_k * dv(i,j)
3041 if (.not. cs%use_reproducing_sums)
then 3045 if (umask(i,j) .eq. 1)
then 3046 dot_p1 = dot_p1 + ru(i,j)**2
3048 if (vmask(i,j) .eq. 1)
then 3049 dot_p1 = dot_p1 + rv(i,j)**2
3053 call mpp_sum (dot_p1)
3059 do j=jsumstart_int_,jecq
3060 do i=isumstart_int_,iecq
3061 if (umask(i,j) .eq. 1) sum_vec(i,j) = ru(i,j)**2
3062 if (vmask(i,j) .eq. 1) sum_vec(i,j) = sum_vec(i,j) + rv(i,j)**2
3066 dot_p1 = reproducing_sum( sum_vec, isumstart_int_, iecq, &
3067 jsumstart_int_, jecq )
3074 dot_p1 = sqrt(dot_p1)
3080 if (dot_p1 .le. cs%cg_tolerance * resid0)
then 3086 cg_halo = cg_halo - 1
3088 if (cg_halo .eq. 0)
then 3090 call pass_vector(du, dv, g%domain, to_all, bgrid_ne)
3091 call pass_vector(u, v, g%domain, to_all, bgrid_ne)
3092 call pass_vector(ru, rv, g%domain, to_all, bgrid_ne)
3100 if (umask(i,j) .eq. 3)
then 3101 u(i,j) = u_bdry(i,j)
3102 elseif (umask(i,j) .eq. 0)
then 3106 if (vmask(i,j) .eq. 3)
then 3107 v(i,j) = v_bdry(i,j)
3108 elseif (vmask(i,j) .eq. 0)
then 3114 call pass_vector (u,v, g%domain, to_all, bgrid_ne)
3116 if (conv_flag .eq. 0)
then 3117 iters = cs%cg_max_iterations
3124 real,
intent(in) :: time_step
3125 real,
dimension(:,:),
intent(in) :: h0
3126 real,
dimension(:,:),
intent(inout) :: h_after_uflux
3127 real,
dimension(:,:,:),
intent(inout) :: flux_enter
3147 integer :: isym, i, j, is, ie, js, je, isd, ied, jsd, jed, gjed, gied
3148 integer :: i_off, j_off
3149 logical :: at_east_bdry, at_west_bdry, one_off_west_bdry, one_off_east_bdry
3150 type(ocean_grid_type),
pointer :: G
3151 real,
dimension(-2:2) :: stencil
3152 real,
dimension(:,:),
pointer :: hmask, u_face_mask, u_flux_boundary_values
3154 flux_diff_cell, phi, dxh, dyh, dxdyh
3156 character (len=1) :: debug_str, procnum
3168 u_face_mask => cs%u_face_mask
3169 u_flux_boundary_values => cs%u_flux_boundary_values
3170 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
3171 i_off = g%idg_offset ; j_off = g%jdg_offset
3174 if (((j+j_off) .le. g%domain%njglobal+g%domain%njhalo) .AND. &
3175 ((j+j_off) .ge. g%domain%njhalo+1))
then 3181 if (((i+i_off) .le. g%domain%niglobal+g%domain%nihalo) .AND. &
3182 ((i+i_off) .ge. g%domain%nihalo+1))
then 3184 if (i+i_off .eq. g%domain%nihalo+1)
then 3187 at_west_bdry=.false.
3190 if (i+i_off .eq. g%domain%niglobal+g%domain%nihalo)
then 3193 at_east_bdry=.false.
3196 if (hmask(i,j) .eq. 1)
then 3198 dxh = g%dxT(i,j) ; dyh = g%dyT(i,j) ; dxdyh = g%areaT(i,j)
3200 h_after_uflux(i,j) = h0(i,j)
3202 stencil(:) = h0(i-2:i+2,j)
3208 if (u_face_mask(i-1,j) .eq. 4.)
then 3210 flux_diff_cell = flux_diff_cell + dyh * time_step * u_flux_boundary_values(i-1,j) / dxdyh
3215 u_face = 0.5 * (cs%u_shelf(i-1,j-1) + cs%u_shelf(i-1,j))
3221 if (u_face .gt. 0)
then 3225 if (at_west_bdry .AND. (hmask(i-1,j).eq.3))
then 3227 stencil(-1) = cs%thickness_boundary_values(i-1,j)
3228 flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step * stencil(-1) / dxdyh
3230 elseif (hmask(i-1,j) * hmask(i-2,j) .eq. 1)
then 3231 phi =
slope_limiter(stencil(-1)-stencil(-2), stencil(0)-stencil(-1))
3232 flux_diff_cell = flux_diff_cell + abs(u_face) * dyh* time_step / dxdyh * &
3233 (stencil(-1) - phi * (stencil(-1)-stencil(0))/2)
3239 flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step / dxdyh * stencil(-1)
3243 elseif (u_face .lt. 0)
then 3244 if (hmask(i-1,j) * hmask(i+1,j) .eq. 1)
then 3245 phi =
slope_limiter(stencil(0)-stencil(1), stencil(-1)-stencil(0))
3246 flux_diff_cell = flux_diff_cell - abs(u_face) * dyh * time_step / dxdyh * &
3247 (stencil(0) - phi * (stencil(0)-stencil(-1))/2)
3250 flux_diff_cell = flux_diff_cell - abs(u_face) * dyh * time_step / dxdyh * stencil(0)
3252 if ((hmask(i-1,j) .eq. 0) .OR. (hmask(i-1,j) .eq. 2))
then 3253 flux_enter(i-1,j,2) = abs(u_face) * dyh * time_step * stencil(0)
3263 if (u_face_mask(i+1,j) .eq. 4.)
then 3265 flux_diff_cell = flux_diff_cell + dyh * time_step * u_flux_boundary_values(i+1,j) / dxdyh
3269 u_face = 0.5 * (cs%u_shelf(i,j-1) + cs%u_shelf(i,j))
3271 if (u_face .lt. 0)
then 3273 if (at_east_bdry .AND. (hmask(i+1,j).eq.3))
then 3276 flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step * stencil(1) / dxdyh
3278 elseif (hmask(i+1,j) * hmask(i+2,j) .eq. 1)
then 3280 phi =
slope_limiter(stencil(1)-stencil(2), stencil(0)-stencil(1))
3281 flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step / dxdyh * &
3282 (stencil(1) - phi * (stencil(1)-stencil(0))/2)
3288 flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step / dxdyh * stencil(1)
3292 elseif (u_face .gt. 0)
then 3294 if (hmask(i-1,j) * hmask(i+1,j) .eq. 1)
then 3296 phi =
slope_limiter(stencil(0)-stencil(-1), stencil(1)-stencil(0))
3297 flux_diff_cell = flux_diff_cell - abs(u_face) * dyh * time_step / dxdyh * &
3298 (stencil(0) - phi * (stencil(0)-stencil(1))/2)
3304 flux_diff_cell = flux_diff_cell - abs(u_face) * dyh * time_step / dxdyh * stencil(0)
3306 if ((hmask(i+1,j) .eq. 0) .OR. (hmask(i+1,j) .eq. 2))
then 3307 flux_enter(i+1,j,1) = abs(u_face) * dyh * time_step * stencil(0)
3314 h_after_uflux(i,j) = h_after_uflux(i,j) + flux_diff_cell
3318 elseif ((hmask(i,j) .eq. 0) .OR. (hmask(i,j) .eq. 2))
then 3320 if (at_west_bdry .AND. (hmask(i-1,j) .EQ. 3))
then 3321 u_face = 0.5 * (cs%u_shelf(i-1,j-1) + cs%u_shelf(i-1,j))
3322 flux_enter(i,j,1) = abs(u_face) * g%dyT(i,j) * time_step * cs%thickness_boundary_values(i-1,j)
3323 elseif (u_face_mask(i-1,j) .eq. 4.)
then 3324 flux_enter(i,j,1) = g%dyT(i,j) * time_step * u_flux_boundary_values(i-1,j)
3327 if (at_east_bdry .AND. (hmask(i+1,j) .EQ. 3))
then 3328 u_face = 0.5 * (cs%u_shelf(i,j-1) + cs%u_shelf(i,j))
3329 flux_enter(i,j,2) = abs(u_face) * g%dyT(i,j) * time_step * cs%thickness_boundary_values(i+1,j)
3330 elseif (u_face_mask(i+1,j) .eq. 4.)
then 3331 flux_enter(i,j,2) = g%dyT(i,j) * time_step * u_flux_boundary_values(i+1,j)
3334 if ((i .eq. is) .AND. (hmask(i,j) .eq. 0) .AND. (hmask(i-1,j) .eq. 1))
then 3339 elseif ((i .eq. ie) .AND. (hmask(i,j) .eq. 0) .AND. (hmask(i+1,j) .eq. 1))
then 3363 real,
intent(in) :: time_step
3364 real,
dimension(:,:),
intent(in) :: h_after_uflux
3365 real,
dimension(:,:),
intent(inout) :: h_after_vflux
3366 real,
dimension(:,:,:),
intent(inout) :: flux_enter
3386 integer :: isym, i, j, is, ie, js, je, isd, ied, jsd, jed, gjed, gied
3387 integer :: i_off, j_off
3388 logical :: at_north_bdry, at_south_bdry, one_off_west_bdry, one_off_east_bdry
3389 type(ocean_grid_type),
pointer :: G
3390 real,
dimension(-2:2) :: stencil
3391 real,
dimension(:,:),
pointer :: hmask, v_face_mask, v_flux_boundary_values
3393 flux_diff_cell, phi, dxh, dyh, dxdyh
3394 character(len=1) :: debug_str, procnum
3406 v_face_mask => cs%v_face_mask
3407 v_flux_boundary_values => cs%v_flux_boundary_values
3408 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
3409 i_off = g%idg_offset ; j_off = g%jdg_offset
3412 if (((i+i_off) .le. g%domain%niglobal+g%domain%nihalo) .AND. &
3413 ((i+i_off) .ge. g%domain%nihalo+1))
then 3419 if (((j+j_off) .le. g%domain%njglobal+g%domain%njhalo) .AND. &
3420 ((j+j_off) .ge. g%domain%njhalo+1))
then 3422 if (j+j_off .eq. g%domain%njhalo+1)
then 3423 at_south_bdry=.true.
3425 at_south_bdry=.false.
3428 if (j+j_off .eq. g%domain%njglobal+g%domain%njhalo)
then 3429 at_north_bdry=.true.
3431 at_north_bdry=.false.
3434 if (hmask(i,j) .eq. 1)
then 3435 dxh = g%dxT(i,j) ; dyh = g%dyT(i,j) ; dxdyh = g%areaT(i,j)
3436 h_after_vflux(i,j) = h_after_uflux(i,j)
3438 stencil(:) = h_after_uflux(i,j-2:j+2)
3443 if (v_face_mask(i,j-1) .eq. 4.)
then 3445 flux_diff_cell = flux_diff_cell + dxh * time_step * v_flux_boundary_values(i,j-1) / dxdyh
3450 v_face = 0.5 * (cs%v_shelf(i-1,j-1) + cs%v_shelf(i,j-1))
3452 if (v_face .gt. 0)
then 3456 if (at_south_bdry .AND. (hmask(i,j-1).eq.3))
then 3458 flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step * stencil(-1) / dxdyh
3460 elseif (hmask(i,j-1) * hmask(i,j-2) .eq. 1)
then 3462 phi =
slope_limiter(stencil(-1)-stencil(-2), stencil(0)-stencil(-1))
3463 flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step / dxdyh * &
3464 (stencil(-1) - phi * (stencil(-1)-stencil(0))/2)
3469 flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step / dxdyh * stencil(-1)
3472 elseif (v_face .lt. 0)
then 3474 if (hmask(i,j-1) * hmask(i,j+1) .eq. 1)
then 3475 phi =
slope_limiter(stencil(0)-stencil(1), stencil(-1)-stencil(0))
3476 flux_diff_cell = flux_diff_cell - abs(v_face) * dxh * time_step / dxdyh * &
3477 (stencil(0) - phi * (stencil(0)-stencil(-1))/2)
3479 flux_diff_cell = flux_diff_cell - abs(v_face) * dxh * time_step / dxdyh * stencil(0)
3481 if ((hmask(i,j-1) .eq. 0) .OR. (hmask(i,j-1) .eq. 2))
then 3482 flux_enter(i,j-1,4) = abs(v_face) * dyh * time_step * stencil(0)
3493 if (v_face_mask(i,j+1) .eq. 4.)
then 3495 flux_diff_cell = flux_diff_cell + dxh * time_step * v_flux_boundary_values(i,j+1) / dxdyh
3500 v_face = 0.5 * (cs%v_shelf(i-1,j) + cs%v_shelf(i,j))
3502 if (v_face .lt. 0)
then 3504 if (at_north_bdry .AND. (hmask(i,j+1).eq.3))
then 3506 flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step * stencil(1) / dxdyh
3507 elseif (hmask(i,j+1) * hmask(i,j+2) .eq. 1)
then 3508 phi =
slope_limiter(stencil(1)-stencil(2), stencil(0)-stencil(1))
3509 flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step / dxdyh * &
3510 (stencil(1) - phi * (stencil(1)-stencil(0))/2)
3514 flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step / dxdyh * stencil(1)
3517 elseif (v_face .gt. 0)
then 3519 if (hmask(i,j-1) * hmask(i,j+1) .eq. 1)
then 3520 phi =
slope_limiter(stencil(0)-stencil(-1), stencil(1)-stencil(0))
3521 flux_diff_cell = flux_diff_cell - abs(v_face) * dxh * time_step / dxdyh * &
3522 (stencil(0) - phi * (stencil(0)-stencil(1))/2)
3526 flux_diff_cell = flux_diff_cell - abs(v_face) * dxh * time_step / dxdyh * stencil(0)
3527 if ((hmask(i,j+1) .eq. 0) .OR. (hmask(i,j+1) .eq. 2))
then 3528 flux_enter(i,j+1,3) = abs(v_face) * dxh * time_step * stencil(0)
3536 h_after_vflux(i,j) = h_after_vflux(i,j) + flux_diff_cell
3538 elseif ((hmask(i,j) .eq. 0) .OR. (hmask(i,j) .eq. 2))
then 3540 if (at_south_bdry .AND. (hmask(i,j-1) .EQ. 3))
then 3541 v_face = 0.5 * (cs%u_shelf(i-1,j-1) + cs%u_shelf(i,j-1))
3542 flux_enter(i,j,3) = abs(v_face) * g%dxT(i,j) * time_step * cs%thickness_boundary_values(i,j-1)
3543 elseif (v_face_mask(i,j-1) .eq. 4.)
then 3544 flux_enter(i,j,3) = g%dxT(i,j) * time_step * v_flux_boundary_values(i,j-1)
3547 if (at_north_bdry .AND. (hmask(i,j+1) .EQ. 3))
then 3548 v_face = 0.5 * (cs%u_shelf(i-1,j) + cs%u_shelf(i,j))
3549 flux_enter(i,j,4) = abs(v_face) * g%dxT(i,j) * time_step * cs%thickness_boundary_values(i,j+1)
3550 elseif (v_face_mask(i,j+1) .eq. 4.)
then 3551 flux_enter(i,j,4) = g%dxT(i,j) * time_step * v_flux_boundary_values(i,j+1)
3554 if ((j .eq. js) .AND. (hmask(i,j) .eq. 0) .AND. (hmask(i,j-1) .eq. 1))
then 3558 elseif ((j .eq. je) .AND. (hmask(i,j) .eq. 0) .AND. (hmask(i,j+1) .eq. 1))
then 3576 real,
dimension(:,:,:),
intent(inout) :: flux_enter
3605 integer :: i, j, isc, iec, jsc, jec, n_flux, k, l, iter_count, isym
3606 integer :: i_off, j_off
3607 integer :: iter_flag
3608 type(ocean_grid_type),
pointer :: G
3609 real,
dimension(:,:),
pointer :: hmask, mass_shelf, area_shelf_h, u_face_mask, v_face_mask, h_shelf
3610 real :: h_reference, dxh, dyh, dxdyh, rho, partial_vol, tot_flux
3611 integer,
dimension(4) :: mapi, mapj, new_partial
3613 real,
dimension (:,:,:),
pointer :: flux_enter_replace => null()
3616 h_shelf => cs%h_shelf
3618 mass_shelf => cs%mass_shelf
3619 area_shelf_h => cs%area_shelf_h
3620 u_face_mask => cs%u_face_mask
3621 v_face_mask => cs%v_face_mask
3622 isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
3623 i_off = g%idg_offset ; j_off = g%jdg_offset
3624 rho = cs%density_ice
3625 iter_count = 0 ; iter_flag = 1
3635 mapi(1) = -1 ; mapi(2) = 1 ; mapi(3:4) = 0
3636 mapj(3) = -1 ; mapj(4) = 1 ; mapj(1:2) = 0
3638 do while (iter_flag .eq. 1)
3642 if (iter_count .gt. 0)
then 3643 flux_enter(:,:,:) = flux_enter_replace(:,:,:)
3644 flux_enter_replace(:,:,:) = 0.0
3647 iter_count = iter_count + 1
3655 if (((j+j_off) .le. g%domain%njglobal+g%domain%njhalo) .AND. &
3656 ((j+j_off) .ge. g%domain%njhalo+1))
then 3660 if (((i+i_off) .le. g%domain%niglobal+g%domain%nihalo) .AND. &
3661 ((i+i_off) .ge. g%domain%nihalo+1))
then 3668 if (flux_enter(i,j,k) .gt. 0)
then 3670 h_reference = h_reference + h_shelf(i+2*k-3,j)
3671 tot_flux = tot_flux + flux_enter(i,j,k)
3672 flux_enter(i,j,k) = 0.0
3677 if (flux_enter(i,j,k+2) .gt. 0)
then 3679 h_reference = h_reference + h_shelf(i,j+2*k-3)
3680 tot_flux = tot_flux + flux_enter(i,j,k+2)
3681 flux_enter(i,j,k+2) = 0.0
3685 if (n_flux .gt. 0)
then 3686 dxdyh = g%areaT(i,j)
3687 h_reference = h_reference /
real(n_flux)
3688 partial_vol = h_shelf(i,j) * area_shelf_h(i,j) + tot_flux
3690 if ((partial_vol / dxdyh) .eq. h_reference)
then 3692 h_shelf(i,j) = h_reference
3693 area_shelf_h(i,j) = dxdyh
3694 elseif ((partial_vol / dxdyh) .lt. h_reference)
then 3697 area_shelf_h(i,j) = partial_vol / h_reference
3698 h_shelf(i,j) = h_reference
3700 if (.not.
associated (flux_enter_replace))
then 3701 allocate ( flux_enter_replace(g%isd:g%ied,g%jsd:g%jed,1:4) )
3702 flux_enter_replace(:,:,:) = 0.0
3706 area_shelf_h(i,j) = dxdyh
3708 partial_vol = partial_vol - h_reference * dxdyh
3712 n_flux = 0 ; new_partial(:) = 0
3715 if (u_face_mask(i-2+k,j) .eq. 2)
then 3717 elseif (hmask(i+2*k-3,j) .eq. 0)
then 3723 if (v_face_mask(i,j-2+k) .eq. 2)
then 3725 elseif (hmask(i,j+2*k-3) .eq. 0)
then 3727 new_partial(k+2) = 1
3731 if (n_flux .eq. 0)
then 3732 h_shelf(i,j) = h_reference + partial_vol / dxdyh
3734 h_shelf(i,j) = h_reference
3737 if (new_partial(k) .eq. 1) &
3738 flux_enter_replace(i+2*k-3,j,3-k) = partial_vol /
real(n_flux)
3741 if (new_partial(k+2) .eq. 1) &
3742 flux_enter_replace(i,j+2*k-3,5-k) = partial_vol /
real(n_flux)
3758 call mpp_max(iter_count)
3760 if(is_root_pe() .and. (iter_count.gt.1)) print *, iter_count,
"MAX ITERATIONS,ADVANCE FRONT" 3762 if (
associated(flux_enter_replace))
DEALLOCATE(flux_enter_replace)
3769 real,
dimension(:,:),
intent(inout) :: h_shelf, area_shelf_h, hmask
3770 type(ocean_grid_type),
pointer :: G
3779 if ((h_shelf(i,j) .lt. cs%min_thickness_simple_calve) .and. (area_shelf_h(i,j).gt. 0.))
then 3781 area_shelf_h(i,j) = 0.0
3789 subroutine calve_to_mask (CS, h_shelf, area_shelf_h, hmask, calve_mask)
3791 real,
dimension(:,:),
intent(inout) :: h_shelf, area_shelf_h, hmask, calve_mask
3793 type(ocean_grid_type),
pointer :: G
3798 if (cs%calve_to_mask)
then 3801 if ((calve_mask(i,j) .eq. 0.0) .and. (hmask(i,j) .ne. 0.0))
then 3803 area_shelf_h(i,j) = 0.0
3814 real,
dimension(:,:),
intent(in) :: OD
3815 real,
dimension(NILIMB_SYM_,NJLIMB_SYM_),
intent(inout) :: TAUD_X, TAUD_Y
3816 integer,
intent(in) :: FE
3831 real,
dimension (:,:),
pointer :: D, &
3833 hmask, u_face_mask, v_face_mask, float_frac
3834 real,
dimension (SIZE(OD,1),SIZE(OD,2)) :: S, &
3836 character(1) :: procnum
3839 real :: rho, rhow, sx, sy, neumann_val, dxh, dyh, dxdyh
3841 type(ocean_grid_type),
pointer :: G
3842 integer :: isym, i, j, iscq, iecq, jscq, jecq, isd, jsd, is, js, iegq, jegq, giec, gjec, gisc, gjsc, cnt, isc, jsc, iec, jec
3843 integer :: i_off, j_off
3848 isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
3849 iscq = g%iscB ; iecq = g%iecB ; jscq = g%jscB ; jecq = g%jecB
3850 isd = g%isd ; jsd = g%jsd
3851 iegq = g%iegB ; jegq = g%jegB
3852 gisc = g%domain%nihalo+1 ; gjsc = g%domain%njhalo+1
3853 giec = g%domain%niglobal+g%domain%nihalo ; gjec = g%domain%njglobal+g%domain%njhalo
3854 is = iscq - (1-isym); js = jscq - (1-isym)
3855 i_off = g%idg_offset ; j_off = g%jdg_offset
3859 float_frac => cs%float_frac
3861 u_face_mask => cs%u_face_mask
3862 v_face_mask => cs%v_face_mask
3863 rho = cs%density_ice
3864 rhow = cs%density_ocean_avg
3866 call savearray2 (
"H",h,cs%write_output_to_file)
3868 call savearray2 (
"u_face_mask", cs%u_face_mask_boundary,cs%write_output_to_file)
3869 call savearray2 (
"umask", cs%umask,cs%write_output_to_file)
3870 call savearray2 (
"v_face_mask", cs%v_face_mask_boundary,cs%write_output_to_file)
3871 call savearray2 (
"vmask", cs%vmask,cs%write_output_to_file)
3884 base(:,:) = -d(:,:) + od(:,:)
3885 s(:,:) = base(:,:) + h(:,:)
3896 dxdyh = g%areaT(i,j)
3899 if (hmask(i,j) .eq. 1)
then 3902 if ((i+i_off) .eq. gisc)
then 3903 if (hmask(i+1,j) .eq. 1)
then 3904 sx = (s(i+1,j)-s(i,j))/dxh
3908 elseif ((i+i_off) .eq. giec)
then 3909 if (hmask(i-1,j) .eq. 1)
then 3910 sx = (s(i,j)-s(i-1,j))/dxh
3915 if (hmask(i+1,j) .eq. 1)
then 3921 if (hmask(i-1,j) .eq. 1)
then 3927 if (cnt .eq. 0)
then 3930 sx = sx / (cnt * dxh)
3937 if ((j+j_off) .eq. gjsc)
then 3938 if (hmask(i,j+1) .eq. 1)
then 3939 sy = (s(i,j+1)-s(i,j))/dyh
3943 elseif ((j+j_off) .eq. gjec)
then 3944 if (hmask(i,j-1) .eq. 1)
then 3945 sy = (s(i,j)-s(i,j-1))/dyh
3950 if (hmask(i,j+1) .eq. 1)
then 3956 if (hmask(i,j-1) .eq. 1)
then 3962 if (cnt .eq. 0)
then 3965 sy = sy / (cnt * dyh)
3973 taud_x(i-1,j-1) = taud_x(i-1,j-1) - .25 * rho * grav * h(i,j) * sx * dxdyh
3974 taud_y(i-1,j-1) = taud_y(i-1,j-1) - .25 * rho * grav * h(i,j) * sy * dxdyh
3977 taud_x(i,j-1) = taud_x(i,j-1) - .25 * rho * grav * h(i,j) * sx * dxdyh
3978 taud_y(i,j-1) = taud_y(i,j-1) - .25 * rho * grav * h(i,j) * sy * dxdyh
3981 taud_x(i-1,j) = taud_x(i-1,j) - .25 * rho * grav * h(i,j) * sx * dxdyh
3982 taud_y(i-1,j) = taud_y(i-1,j) - .25 * rho * grav * h(i,j) * sy * dxdyh
3985 taud_x(i,j) = taud_x(i,j) - .25 * rho * grav * h(i,j) * sx * dxdyh
3986 taud_y(i,j) = taud_y(i,j) - .25 * rho * grav * h(i,j) * sy * dxdyh
3992 taud_x(i-1,j-1) = taud_x(i-1,j-1) - (1./6) * rho * grav * h(i,j) * sx * dxdyh
3993 taud_y(i-1,j-1) = taud_y(i-1,j-1) - (1./6) * rho * grav * h(i,j) * sy * dxdyh
3996 taud_x(i,j-1) = taud_x(i,j-1) - (1./3) * rho * grav * h(i,j) * sx * dxdyh
3997 taud_y(i,j-1) = taud_y(i,j-1) - (1./3) * rho * grav * h(i,j) * sy * dxdyh
4000 taud_x(i-1,j) = taud_x(i-1,j) - (1./3) * rho * grav * h(i,j) * sx * dxdyh
4001 taud_y(i-1,j) = taud_y(i-1,j) - (1./3) * rho * grav * h(i,j) * sy * dxdyh
4004 taud_x(i,j) = taud_x(i,j) - (1./6) * rho * grav * h(i,j) * sx * dxdyh
4005 taud_y(i,j) = taud_y(i,j) - (1./6) * rho * grav * h(i,j) * sy * dxdyh
4009 if (float_frac(i,j) .eq. 1)
then 4010 neumann_val = .5 * grav * (rho * h(i,j) ** 2 - rhow * d(i,j) ** 2)
4012 neumann_val = .5 * grav * (1-rho/rhow) * rho * h(i,j) ** 2
4016 if ((u_face_mask(i-1,j) .eq. 2) .OR. (hmask(i-1,j) .eq. 0) .OR. (hmask(i-1,j) .eq. 2) )
then 4022 taud_x(i-1,j-1) = taud_x(i-1,j-1) - .5 * dyh * neumann_val
4023 taud_x(i-1,j) = taud_x(i-1,j) - .5 * dyh * neumann_val
4026 if ((u_face_mask(i,j) .eq. 2) .OR. (hmask(i+1,j) .eq. 0) .OR. (hmask(i+1,j) .eq. 2) )
then 4027 taud_x(i,j-1) = taud_x(i,j-1) + .5 * dyh * neumann_val
4028 taud_x(i,j) = taud_x(i,j) + .5 * dyh * neumann_val
4031 if ((v_face_mask(i,j-1) .eq. 2) .OR. (hmask(i,j-1) .eq. 0) .OR. (hmask(i,j-1) .eq. 2) )
then 4032 taud_y(i-1,j-1) = taud_y(i-1,j-1) - .5 * dxh * neumann_val
4033 taud_y(i,j-1) = taud_y(i,j-1) - .5 * dxh * neumann_val
4036 if ((v_face_mask(i,j) .eq. 2) .OR. (hmask(i,j+1) .eq. 0) .OR. (hmask(i,j+1) .eq. 2) )
then 4037 taud_y(i-1,j) = taud_y(i-1,j) + .5 * dxh * neumann_val
4038 taud_y(i,j) = taud_y(i,j) + .5 * dxh * neumann_val
4052 type(time_type),
intent(in) :: Time
4054 real,
intent(in) :: input_flux, input_thick
4055 logical,
optional :: new_sim
4065 real,
dimension (:,:) ,
pointer :: thickness_boundary_values, &
4066 u_boundary_values, &
4067 v_boundary_values, &
4068 u_face_mask, v_face_mask, hmask
4069 type(ocean_grid_type),
pointer :: G
4070 integer :: isym, i, j, iscq, iecq, jscq, jecq, isd, jsd, ied, jed, iegq, jegq, giec, gjec, gisc, gjsc, cnt, isc, jsc, iec, jec
4071 integer :: i_off, j_off
4072 real :: A, n, ux, uy, vx, vy, eps_min, domain_width
4084 isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
4086 isd = g%isd ; jsd = g%jsd ; ied = g%ied ; jed = g%jed
4088 i_off = g%idg_offset ; j_off = g%jdg_offset
4090 thickness_boundary_values => cs%thickness_boundary_values
4091 u_boundary_values => cs%u_boundary_values ; v_boundary_values => cs%v_boundary_values
4092 u_face_mask => cs%u_face_mask ; v_face_mask => cs%v_face_mask ; hmask => cs%hmask
4094 domain_width = cs%len_lat
4105 if (hmask(i,j) .eq. 3)
then 4106 thickness_boundary_values(i,j) = input_thick
4109 if ((hmask(i,j) .eq. 0) .or. (hmask(i,j) .eq. 1) .or. (hmask(i,j) .eq. 2))
then 4110 if ((i.le.iec).and.(i.ge.isc))
then 4111 if (u_face_mask(i-1,j) .eq. 3)
then 4112 u_boundary_values(i-1,j-1) = (1 - ((g%geoLatBu(i-1,j-1) - 0.5*cs%len_lat)*2./cs%len_lat)**2) * &
4113 1.5 * input_flux / input_thick
4114 u_boundary_values(i-1,j) = (1 - ((g%geoLatBu(i-1,j) - 0.5*cs%len_lat)*2./cs%len_lat)**2) * &
4115 1.5 * input_flux / input_thick
4120 if (.not.(new_sim))
then 4121 if (.not. g%symmetric)
then 4122 if (((i+i_off) .eq. (g%domain%nihalo+1)).and.(u_face_mask(i-1,j).eq.3))
then 4123 cs%u_shelf (i-1,j-1) = u_boundary_values(i-1,j-1)
4124 cs%u_shelf (i-1,j) = u_boundary_values(i-1,j)
4127 if (((j+j_off) .eq. (g%domain%njhalo+1)).and.(v_face_mask(i,j-1).eq.3))
then 4128 cs%u_shelf (i-1,j-1) = u_boundary_values(i-1,j-1)
4129 cs%u_shelf (i,j-1) = u_boundary_values(i,j-1)
4139 beta_upper, beta_lower, dxh, dyh, dxdyh, is, ie, js, je, isym)
4141 real,
dimension (:,:),
intent (inout) :: uret, vret
4142 real,
dimension (:,:),
intent (in) :: u, v
4143 real,
dimension (:,:),
intent (in) :: umask, vmask
4144 real,
dimension (:,:),
intent (in) :: hmask, nu_upper, nu_lower, beta_upper, beta_lower
4145 real,
dimension (:,:),
intent (in) :: dxh, dyh, dxdyh
4146 integer,
intent(in) :: is, ie, js, je, isym
4156 real :: ux, uy, vx, vy
4162 if (hmask(i,j) .eq. 1)
then 4164 ux = (u(i,j-1)-u(i-1,j-1))/dxh(i,j)
4165 vx = (v(i,j-1)-v(i-1,j-1))/dxh(i,j)
4166 uy = (u(i-1,j)-u(i-1,j-1))/dyh(i,j)
4167 vy = (v(i-1,j)-v(i-1,j-1))/dyh(i,j)
4169 if (umask(i,j-1) .eq. 1)
then 4171 uret(i,j-1) = uret(i,j-1) + &
4172 .5 * dxdyh(i,j) * nu_lower(i,j) * ((4*ux+2*vy) * (1./dxh(i,j)) + (uy+vy) * (0./dyh(i,j)))
4174 vret(i,j-1) = vret(i,j-1) + &
4175 .5 * dxdyh(i,j) * nu_lower(i,j) * ((uy+vx) * (1./dxh(i,j)) + (4*vy+2*ux) * (0./dyh(i,j)))
4177 uret(i,j-1) = uret(i,j-1) + &
4178 beta_lower(i,j) * dxdyh(i,j) * 1./24 * (u(i-1,j-1) + &
4179 u(i-1,j) + u(i,j-1))
4181 vret(i,j-1) = vret(i,j-1) + &
4182 beta_lower(i,j) * dxdyh(i,j) * 1./24 * (v(i-1,j-1) + &
4183 v(i-1,j) + v(i,j-1))
4186 if (umask(i-1,j) .eq. 1)
then 4188 uret(i-1,j) = uret(i-1,j) + &
4189 .5 * dxdyh(i,j) * nu_lower(i,j) * ((4*ux+2*vy) * (0./dxh(i,j)) + (uy+vy) * (1./dyh(i,j)))
4191 vret(i-1,j) = vret(i-1,j) + &
4192 .5 * dxdyh(i,j) * nu_lower(i,j) * ((uy+vx) * (0./dxh(i,j)) + (4*vy+2*ux) * (1./dyh(i,j)))
4194 uret(i,j-1) = uret(i,j-1) + &
4195 beta_lower(i,j) * dxdyh(i,j) * 1./24 * (u(i-1,j-1) + &
4196 u(i-1,j) + u(i,j-1))
4198 vret(i,j-1) = vret(i,j-1) + &
4199 beta_lower(i,j) * dxdyh(i,j) * 1./24 * (v(i-1,j-1) + &
4200 v(i-1,j) + v(i,j-1))
4203 if (umask(i-1,j-1) .eq. 1)
then 4205 uret(i-1,j-1) = uret(i-1,j-1) + &
4206 .5 * dxdyh(i,j) * nu_upper(i,j) * ((4*ux+2*vy) * (-1./dxh(i,j)) + (uy+vy) * (-1./dyh(i,j)))
4208 vret(i-1,j-1) = vret(i-1,j-1) + &
4209 .5 * dxdyh(i,j) * nu_upper(i,j) * ((uy+vx) * (-1./dxh(i,j)) + (4*vy+2*ux) * (-1./dyh(i,j)))
4211 uret(i-1,j-1) = uret(i-1,j-1) + &
4212 beta_lower(i,j) * dxdyh(i,j) * 1./24 * (u(i-1,j-1) + &
4213 u(i-1,j) + u(i,j-1))
4215 vret(i-1,j-1) = vret(i-1,j-1) + &
4216 beta_lower(i,j) * dxdyh(i,j) * 1./24 * (v(i-1,j-1) + &
4217 v(i-1,j) + v(i,j-1))
4221 ux = (u(i,j)-u(i-1,j))/dxh(i,j)
4222 vx = (v(i,j)-v(i-1,j))/dxh(i,j)
4223 uy = (u(i,j)-u(i,j-1))/dyh(i,j)
4224 vy = (v(i,j)-v(i,j-1))/dyh(i,j)
4226 if (umask(i,j-1) .eq. 1)
then 4228 uret(i,j-1) = uret(i,j-1) + &
4229 .5 * dxdyh(i,j) * nu_upper(i,j) * ((4*ux+2*vy) * (0./dxh(i,j)) + (uy+vy) * (-1./dyh(i,j)))
4231 vret(i,j-1) = vret(i,j-1) + &
4232 .5 * dxdyh(i,j) * nu_upper(i,j) * ((uy+vx) * (0./dxh(i,j)) + (4*vy+2*ux) * (-1./dyh(i,j)))
4234 uret(i,j-1) = uret(i,j-1) + &
4235 beta_upper(i,j) * dxdyh(i,j) * 1./24 * (u(i,j) + &
4236 u(i-1,j) + u(i,j-1))
4238 vret(i,j-1) = vret(i,j-1) + &
4239 beta_upper(i,j) * dxdyh(i,j) * 1./24 * (u(i,j) + &
4240 u(i-1,j) + u(i,j-1))
4243 if (umask(i-1,j) .eq. 1)
then 4245 uret(i-1,j) = uret(i-1,j) + &
4246 .5 * dxdyh(i,j) * nu_upper(i,j) * ((4*ux+2*vy) * (-1./dxh(i,j)) + (uy+vy) * (0./dyh(i,j)))
4248 vret(i-1,j) = vret(i-1,j) + &
4249 .5 * dxdyh(i,j) * nu_upper(i,j) * ((uy+vx) * (-1./dxh(i,j)) + (4*vy+2*ux) * (0./dyh(i,j)))
4251 uret(i,j-1) = uret(i,j-1) + &
4252 beta_upper(i,j) * dxdyh(i,j) * 1./24 * (u(i,j) + &
4253 u(i-1,j) + u(i,j-1))
4255 vret(i,j-1) = vret(i,j-1) + &
4256 beta_upper(i,j) * dxdyh(i,j) * 1./24 * (u(i,j) + &
4257 u(i-1,j) + u(i,j-1))
4260 if (umask(i,j) .eq. 1)
then 4262 uret(i,j) = uret(i,j) + &
4263 .5 * dxdyh(i,j) * nu_upper(i,j) * ((4*ux+2*vy) * (1./dxh(i,j)) + (uy+vy) * (1./dyh(i,j)))
4265 vret(i,j) = vret(i,j) + &
4266 .5 * dxdyh(i,j) * nu_upper(i,j) * ((uy+vx) * (1./dxh(i,j)) + (4*vy+2*ux) * (1./dyh(i,j)))
4268 uret(i,j) = uret(i,j) + &
4269 beta_upper(i,j) * dxdyh(i,j) * 1./24 * (u(i,j) + &
4270 u(i-1,j) + u(i,j-1))
4272 vret(i,j) = vret(i,j) + &
4273 beta_upper(i,j) * dxdyh(i,j) * 1./24 * (u(i,j) + &
4274 u(i-1,j) + u(i,j-1))
4284 subroutine cg_action_bilinear (uret, vret, u, v, Phi, Phisub, umask, vmask, hmask, H_node, &
4285 nu, float_cond, D, beta, dxdyh, is, ie, js, je, dens_ratio)
4287 real,
dimension (NILIMB_SYM_,NJLIMB_SYM_),
intent (inout) :: uret, vret
4288 real,
dimension (:,:,:,:),
pointer :: Phi
4289 real,
dimension (:,:,:,:,:,:),
pointer :: Phisub
4290 real,
dimension (NILIMB_SYM_,NJLIMB_SYM_),
intent (in) :: u, v
4291 real,
dimension (NILIMB_SYM_,NJLIMB_SYM_),
intent (in) :: umask, vmask, H_node
4292 real,
dimension (:,:),
intent (in) :: hmask, nu, float_cond, D, beta, dxdyh
4293 real,
intent(in) :: dens_ratio
4294 integer,
intent(in) :: is, ie, js, je
4316 real :: ux, vx, uy, vy, uq, vq, area, basel
4317 integer :: iq, jq, iphi, jphi, i, j, ilq, jlq
4318 real,
dimension(2) :: xquad
4319 real,
dimension(2,2) :: Ucell,Vcell,Hcell,Usubcontr,Vsubcontr,Ucontr
4321 xquad(1) = .5 * (1-sqrt(1./3)) ; xquad(2) = .5 * (1+sqrt(1./3))
4324 do i=is,ie ;
if (hmask(i,j) .eq. 1)
then 4343 do iq=1,2 ;
do jq=1,2
4358 uq = u(i-1,j-1) * xquad(3-iq) * xquad(3-jq) + &
4359 u(i,j-1) * xquad(iq) * xquad(3-jq) + &
4360 u(i-1,j) * xquad(3-iq) * xquad(jq) + &
4361 u(i,j) * xquad(iq) * xquad(jq)
4363 vq = v(i-1,j-1) * xquad(3-iq) * xquad(3-jq) + &
4364 v(i,j-1) * xquad(iq) * xquad(3-jq) + &
4365 v(i-1,j) * xquad(3-iq) * xquad(jq) + &
4366 v(i,j) * xquad(iq) * xquad(jq)
4368 ux = u(i-1,j-1) * phi(i,j,1,2*(jq-1)+iq) + &
4369 u(i,j-1) * phi(i,j,3,2*(jq-1)+iq) + &
4370 u(i-1,j) * phi(i,j,5,2*(jq-1)+iq) + &
4371 u(i,j) * phi(i,j,7,2*(jq-1)+iq)
4373 vx = v(i-1,j-1) * phi(i,j,1,2*(jq-1)+iq) + &
4374 v(i,j-1) * phi(i,j,3,2*(jq-1)+iq) + &
4375 v(i-1,j) * phi(i,j,5,2*(jq-1)+iq) + &
4376 v(i,j) * phi(i,j,7,2*(jq-1)+iq)
4378 uy = u(i-1,j-1) * phi(i,j,2,2*(jq-1)+iq) + &
4379 u(i,j-1) * phi(i,j,4,2*(jq-1)+iq) + &
4380 u(i-1,j) * phi(i,j,6,2*(jq-1)+iq) + &
4381 u(i,j) * phi(i,j,8,2*(jq-1)+iq)
4383 vy = v(i-1,j-1) * phi(i,j,2,2*(jq-1)+iq) + &
4384 v(i,j-1) * phi(i,j,4,2*(jq-1)+iq) + &
4385 v(i-1,j) * phi(i,j,6,2*(jq-1)+iq) + &
4386 v(i,j) * phi(i,j,8,2*(jq-1)+iq)
4388 do iphi=1,2 ;
do jphi=1,2
4389 if (umask(i-2+iphi,j-2+jphi) .eq. 1)
then 4391 uret(i-2+iphi,j-2+jphi) = uret(i-2+iphi,j-2+jphi) + &
4392 .25 * area * nu(i,j) * ((4*ux+2*vy) * phi(i,j,2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + &
4393 (uy+vx) * phi(i,j,2*(2*(jphi-1)+iphi),2*(jq-1)+iq))
4395 if (vmask(i-2+iphi,j-2+jphi) .eq. 1)
then 4397 vret(i-2+iphi,j-2+jphi) = vret(i-2+iphi,j-2+jphi) + &
4398 .25 * area * nu(i,j) * ((uy+vx) * phi(i,j,2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + &
4399 (4*vy+2*ux) * phi(i,j,2*(2*(jphi-1)+iphi),2*(jq-1)+iq))
4402 if (iq .eq. iphi)
then 4408 if (jq .eq. jphi)
then 4414 if (float_cond(i,j) .eq. 0)
then 4416 if (umask(i-2+iphi,j-2+jphi) .eq. 1)
then 4418 uret(i-2+iphi,j-2+jphi) = uret(i-2+iphi,j-2+jphi) + &
4419 .25 * beta(i,j) * area * uq * xquad(ilq) * xquad(jlq)
4423 if (vmask(i-2+iphi,j-2+jphi) .eq. 1)
then 4425 vret(i-2+iphi,j-2+jphi) = vret(i-2+iphi,j-2+jphi) + &
4426 .25 * beta(i,j) * area * vq * xquad(ilq) * xquad(jlq)
4431 ucontr(iphi,jphi) = ucontr(iphi,jphi) + .25 * area * uq * xquad(ilq) * xquad(jlq) * beta(i,j)
4438 if (float_cond(i,j) .eq. 1)
then 4439 usubcontr = 0.0 ; vsubcontr = 0.0 ; basel = d(i,j)
4440 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)
4442 (phisub, hcell, ucell, vcell, area, basel, dens_ratio, usubcontr, vsubcontr, i, j)
4443 do iphi=1,2 ;
do jphi=1,2
4444 if (umask(i-2+iphi,j-2+jphi) .eq. 1)
then 4445 uret(i-2+iphi,j-2+jphi) = uret(i-2+iphi,j-2+jphi) + usubcontr(iphi,jphi) * beta(i,j)
4447 if (vmask(i-2+iphi,j-2+jphi) .eq. 1)
then 4448 vret(i-2+iphi,j-2+jphi) = vret(i-2+iphi,j-2+jphi) + vsubcontr(iphi,jphi) * beta(i,j)
4459 subroutine cg_action_subgrid_basal_bilinear (Phisub, H, U, V, DXDYH, D, dens_ratio, Ucontr, Vcontr, iin, jin)
4460 real,
pointer,
dimension(:,:,:,:,:,:) :: Phisub
4461 real,
dimension(2,2),
intent(in) :: H,U,V
4462 real,
intent(in) :: DXDYH, D, dens_ratio
4463 real,
dimension(2,2),
intent(inout) :: Ucontr, Vcontr
4464 integer,
optional,
intent(in) :: iin, jin
4468 integer :: nsub, i, j, k, l, qx, qy, m, n, i_m, j_m
4469 real :: subarea, hloc, uq, vq
4471 nsub =
size(phisub,1)
4472 subarea = dxdyh / (nsub**2)
4475 if (.not.
present(iin))
then 4481 if (.not.
present(jin))
then 4495 hloc = phisub(i,j,1,1,qx,qy)*h(1,1)+phisub(i,j,1,2,qx,qy)*h(1,2)+&
4496 phisub(i,j,2,1,qx,qy)*h(2,1)+phisub(i,j,2,2,qx,qy)*h(2,2)
4498 if (dens_ratio * hloc - d .gt. 0)
then 4505 uq = uq + phisub(i,j,k,l,qx,qy) * u(k,l) ; vq = vq + phisub(i,j,k,l,qx,qy) * v(k,l)
4509 ucontr(m,n) = ucontr(m,n) + subarea * 0.25 * phisub(i,j,m,n,qx,qy) * uq
4510 vcontr(m,n) = vcontr(m,n) + subarea * 0.25 * phisub(i,j,m,n,qx,qy) * vq
4528 real,
dimension (:,:),
intent(inout) :: u_diagonal, v_diagonal
4532 real,
pointer,
dimension (:,:) :: umask, vmask, &
4533 nu_lower, nu_upper, beta_lower, beta_upper, hmask
4534 type(ocean_grid_type),
pointer :: G
4535 integer :: isym, i, j, is, js, cnt, isc, jsc, iec, jec
4536 real :: A, n, ux, uy, vx, vy, eps_min, domain_width, dxh, dyh, dxdyh
4548 isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
4550 umask => cs%umask ; vmask => cs%vmask ; hmask => cs%hmask
4551 nu_lower => cs%ice_visc_lower_tri ; nu_upper => cs%ice_visc_upper_tri
4552 beta_lower => cs%taub_beta_eff_lower_tri ; beta_upper => cs%taub_beta_eff_upper_tri
4554 do i=isc-1,iec+1 ;
do j=jsc-1,jec+1 ;
if (hmask(i,j) .eq. 1)
then 4557 dxdyh = g%areaT(i,j)
4559 if (umask(i,j-1) .eq. 1)
then 4561 ux = 1./dxh ; uy = 0./dyh
4564 u_diagonal(i,j-1) = u_diagonal(i,j-1) + &
4565 .5 * dxdyh * nu_lower(i,j) * ((4*ux+2*vy) * (1./dxh) + (uy+vy) * (0./dyh))
4567 u_diagonal(i,j-1) = u_diagonal(i,j-1) + &
4568 beta_lower(i,j) * dxdyh * 1./24
4571 vx = 1./dxh ; vy = 0./dyh
4573 v_diagonal(i,j-1) = v_diagonal(i,j-1) + &
4574 .5 * dxdyh * nu_lower(i,j) * ((uy+vx) * (1./dxh) + (4*vy+2*ux) * (0./dyh))
4576 v_diagonal(i,j-1) = v_diagonal(i,j-1) + &
4577 beta_lower(i,j) * dxdyh * 1./24
4579 ux = 0./dxh ; uy = -1./dyh
4582 u_diagonal(i,j-1) = u_diagonal(i,j-1) + &
4583 .5 * dxdyh * nu_upper(i,j) * ((4*ux+2*vy) * (0./dxh) + (uy+vy) * (-1./dyh))
4585 u_diagonal(i,j-1) = u_diagonal(i,j-1) + &
4586 beta_upper(i,j) * dxdyh * 1./24
4588 vx = 0./dxh ; vy = -1./dyh
4591 v_diagonal(i,j-1) = v_diagonal(i,j-1) + &
4592 .5 * dxdyh * nu_upper(i,j) * ((uy+vx) * (0./dxh) + (4*vy+2*ux) * (-1./dyh))
4594 v_diagonal(i,j-1) = v_diagonal(i,j-1) + &
4595 beta_upper(i,j) * dxdyh * 1./24
4599 if (umask(i-1,j) .eq. 1)
then 4601 ux = 0./dxh ; uy = 1./dyh
4604 u_diagonal(i-1,j) = u_diagonal(i-1,j) + &
4605 .5 * dxdyh * nu_lower(i,j) * ((4*ux+2*vy) * (0./dxh) + (uy+vy) * (1./dyh))
4607 u_diagonal(i,j-1) = u_diagonal(i,j-1) + &
4608 beta_lower(i,j) * dxdyh * 1./24
4611 vx = 0./dxh ; vy = 1./dyh
4613 v_diagonal(i-1,j) = v_diagonal(i-1,j) + &
4614 .5 * dxdyh * nu_lower(i,j) * ((uy+vx) * (0./dxh) + (4*vy+2*ux) * (1./dyh))
4616 v_diagonal(i,j-1) = v_diagonal(i,j-1) + &
4617 beta_lower(i,j) * dxdyh * 1./24
4619 ux = -1./dxh ; uy = 0./dyh
4622 u_diagonal(i-1,j) = u_diagonal(i-1,j) + &
4623 .5 * dxdyh * nu_upper(i,j) * ((4*ux+2*vy) * (-1./dxh) + (uy+vy) * (0./dyh))
4625 u_diagonal(i,j-1) = u_diagonal(i,j-1) + &
4626 beta_upper(i,j) * dxdyh * 1./24
4628 vx = -1./dxh ; vy = 0./dyh
4631 v_diagonal(i-1,j) = v_diagonal(i-1,j) + &
4632 .5 * dxdyh * nu_upper(i,j) * ((uy+vx) * (-1./dxh) + (4*vy+2*ux) * (0./dyh))
4634 v_diagonal(i,j-1) = v_diagonal(i,j-1) + &
4635 beta_upper(i,j) * dxdyh * 1./24
4639 if (umask(i-1,j-1) .eq. 1)
then 4641 ux = -1./dxh ; uy = -1./dyh
4644 u_diagonal(i-1,j-1) = u_diagonal(i-1,j-1) + &
4645 .5 * dxdyh * nu_upper(i,j) * ((4*ux+2*vy) * (-1./dxh) + (uy+vy) * (-1./dyh))
4647 u_diagonal(i-1,j-1) = u_diagonal(i-1,j-1) + &
4648 beta_lower(i,j) * dxdyh * 1./24
4650 vx = -1./dxh ; vy = -1./dyh
4653 v_diagonal(i-1,j-1) = v_diagonal(i-1,j-1) + &
4654 .5 * dxdyh * nu_upper(i,j) * ((uy+vx) * (-1./dxh) + (4*vy+2*ux) * (-1./dyh))
4656 v_diagonal(i-1,j-1) = v_diagonal(i-1,j-1) + &
4657 beta_lower(i,j) * dxdyh * 1./24
4660 if (umask(i,j) .eq. 1)
then 4662 ux = 1./ dxh ; uy = 1./dyh
4665 u_diagonal(i,j) = u_diagonal(i,j) + &
4666 .5 * dxdyh * nu_upper(i,j) * ((4*ux+2*vy) * (1./dxh) + (uy+vy) * (1./dyh))
4668 u_diagonal(i,j) = u_diagonal(i,j) + &
4669 beta_upper(i,j) * dxdyh * 1./24
4671 vx = 1./ dxh ; vy = 1./dyh
4674 v_diagonal(i,j) = v_diagonal(i,j) + &
4675 .5 * dxdyh * nu_upper(i,j) * ((uy+vx) * (1./dxh) + (4*vy+2*ux) * (1./dyh))
4677 v_diagonal(i,j) = v_diagonal(i,j) + &
4678 beta_upper(i,j) * dxdyh * 1./24
4681 endif ;
enddo ;
enddo 4688 real,
dimension (NILIMB_SYM_,NJLIMB_SYM_),
intent(in) :: H_node
4690 real,
dimension (:,:),
intent(in) :: float_cond
4691 real,
dimension (:,:,:,:,:,:),
pointer :: Phisub
4692 real,
dimension (NILIMB_SYM_,NJLIMB_SYM_),
intent(inout) :: u_diagonal, v_diagonal
4697 real,
dimension (:,:),
pointer :: umask, vmask, hmask, &
4699 type(ocean_grid_type),
pointer :: G
4700 integer :: isym, i, j, is, js, cnt, isc, jsc, iec, jec, iphi, jphi, iq, jq, ilq, jlq
4701 real :: A, n, ux, uy, vx, vy, eps_min, domain_width, dxh, dyh, dxdyh, area, uq, vq, basel
4702 real,
dimension(8,4) :: Phi
4703 real,
dimension(4) :: X, Y
4704 real,
dimension(2) :: xquad
4705 real,
dimension(2,2) :: Hcell,Usubcontr,Vsubcontr
4717 isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
4719 umask => cs%umask ; vmask => cs%vmask ; hmask => cs%hmask
4720 nu => cs%ice_visc_bilinear
4721 beta => cs%taub_beta_eff_bilinear
4723 xquad(1) = .5 * (1-sqrt(1./3)) ; xquad(2) = .5 * (1+sqrt(1./3))
4732 do j=jsc-1,jec+1 ;
do i=isc-1,iec+1 ;
if (hmask(i,j) .eq. 1)
then 4736 dxdyh = g%areaT(i,j)
4738 x(1:2) = g%geoLonBu (i-1:i,j-1)*1000
4739 x(3:4) = g%geoLonBu (i-1:i,j) *1000
4740 y(1:2) = g%geoLatBu (i-1:i,j-1) *1000
4741 y(3:4) = g%geoLatBu (i-1:i,j)*1000
4752 do iq=1,2 ;
do jq=1,2
4754 do iphi=1,2 ;
do jphi=1,2
4756 if (iq .eq. iphi)
then 4762 if (jq .eq. jphi)
then 4768 if (umask(i-2+iphi,j-2+jphi) .eq. 1)
then 4770 ux = phi(2*(2*(jphi-1)+iphi)-1, 2*(jq-1)+iq)
4771 uy = phi(2*(2*(jphi-1)+iphi), 2*(jq-1)+iq)
4775 u_diagonal(i-2+iphi,j-2+jphi) = u_diagonal(i-2+iphi,j-2+jphi) + &
4776 .25 * dxdyh * nu(i,j) * ((4*ux+2*vy) * phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + &
4777 (uy+vy) * phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq))
4779 uq = xquad(ilq) * xquad(jlq)
4781 if (float_cond(i,j) .eq. 0)
then 4782 u_diagonal(i-2+iphi,j-2+jphi) = u_diagonal(i-2+iphi,j-2+jphi) + &
4783 .25 * beta(i,j) * dxdyh * uq * xquad(ilq) * xquad(jlq)
4788 if (vmask(i-2+iphi,j-2+jphi) .eq. 1)
then 4790 vx = phi(2*(2*(jphi-1)+iphi)-1, 2*(jq-1)+iq)
4791 vy = phi(2*(2*(jphi-1)+iphi), 2*(jq-1)+iq)
4795 v_diagonal(i-2+iphi,j-2+jphi) = v_diagonal(i-2+iphi,j-2+jphi) + &
4796 .25 * dxdyh * nu(i,j) * ((uy+vx) * phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + &
4797 (4*vy+2*ux) * phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq))
4799 vq = xquad(ilq) * xquad(jlq)
4801 if (float_cond(i,j) .eq. 0)
then 4802 v_diagonal(i-2+iphi,j-2+jphi) = v_diagonal(i-2+iphi,j-2+jphi) + &
4803 .25 * beta(i,j) * dxdyh * vq * xquad(ilq) * xquad(jlq)
4809 if (float_cond(i,j) .eq. 1)
then 4810 usubcontr = 0.0 ; vsubcontr = 0.0 ; basel = g%bathyT(i,j)
4811 hcell(:,:) = h_node(i-1:i,j-1:j)
4813 (phisub, hcell, dxdyh, basel, dens_ratio, usubcontr, vsubcontr)
4814 do iphi=1,2 ;
do jphi=1,2
4815 if (umask(i-2+iphi,j-2+jphi) .eq. 1)
then 4816 u_diagonal(i-2+iphi,j-2+jphi) = u_diagonal(i-2+iphi,j-2+jphi) + usubcontr(iphi,jphi) * beta(i,j)
4817 v_diagonal(i-2+iphi,j-2+jphi) = v_diagonal(i-2+iphi,j-2+jphi) + vsubcontr(iphi,jphi) * beta(i,j)
4821 endif ;
enddo ;
enddo 4826 real,
pointer,
dimension(:,:,:,:,:,:) :: Phisub
4827 real,
dimension(2,2),
intent(in) :: H
4828 real,
intent(in) :: DXDYH, D, dens_ratio
4829 real,
dimension(2,2),
intent(inout) :: Ucontr, Vcontr
4833 integer :: nsub, i, j, k, l, qx, qy, m, n
4834 real :: subarea, hloc
4836 nsub =
size(phisub,1)
4837 subarea = dxdyh / (nsub**2)
4846 hloc = phisub(i,j,1,1,qx,qy)*h(1,1)+phisub(i,j,1,2,qx,qy)*h(1,2)+&
4847 phisub(i,j,2,1,qx,qy)*h(2,1)+phisub(i,j,2,2,qx,qy)*h(2,2)
4849 if (dens_ratio * hloc - d .gt. 0)
then 4850 ucontr(m,n) = ucontr(m,n) + subarea * 0.25 * phisub(i,j,m,n,qx,qy)**2
4851 vcontr(m,n) = vcontr(m,n) + subarea * 0.25 * phisub(i,j,m,n,qx,qy)**2
4867 type(time_type),
intent(in) :: Time
4869 real,
dimension (:,:),
intent(inout) :: u_boundary_contr, v_boundary_contr
4874 real,
pointer,
dimension (:,:) :: u_boundary_values, &
4875 v_boundary_values, &
4876 umask, vmask, hmask, &
4877 nu_lower, nu_upper, beta_lower, beta_upper
4878 type(ocean_grid_type),
pointer :: G
4879 integer :: isym, i, j, cnt, isc, jsc, iec, jec
4880 real :: A, n, ux, uy, vx, vy, eps_min, domain_width, dxh, dyh, dxdyh
4892 isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
4894 u_boundary_values => cs%u_boundary_values
4895 v_boundary_values => cs%v_boundary_values
4896 umask => cs%umask ; vmask => cs%vmask ; hmask => cs%hmask
4897 nu_lower => cs%ice_visc_lower_tri ; nu_upper => cs%ice_visc_upper_tri
4898 beta_lower => cs%taub_beta_eff_lower_tri ; beta_upper => cs%taub_beta_eff_upper_tri
4900 domain_width = cs%len_lat
4902 do i=isc-1,iec+1 ;
do j=jsc-1,jec+1 ;
if (hmask(i,j) .eq. 1)
then 4904 if ((umask(i-1,j-1) .eq. 3) .OR. (umask(i,j-1) .eq. 3) .OR. (umask(i-1,j) .eq. 3))
then 4908 dxdyh = g%areaT(i,j)
4910 ux = (u_boundary_values(i,j-1)-u_boundary_values(i-1,j-1))/dxh
4911 vx = (v_boundary_values(i,j-1)-v_boundary_values(i-1,j-1))/dxh
4912 uy = (u_boundary_values(i-1,j)-u_boundary_values(i-1,j-1))/dyh
4913 vy = (v_boundary_values(i-1,j)-v_boundary_values(i-1,j-1))/dyh
4915 if (umask(i,j-1) .eq. 1)
then 4917 u_boundary_contr(i,j-1) = u_boundary_contr(i,j-1) + &
4918 .5 * dxdyh * nu_lower(i,j) * ((4*ux+2*vy) * (1./dxh) + (uy+vy) * (0./dyh))
4920 v_boundary_contr(i,j-1) = v_boundary_contr(i,j-1) + &
4921 .5 * dxdyh * nu_lower(i,j) * ((uy+vx) * (1./dxh) + (4*vy+2*ux) * (0./dyh))
4923 u_boundary_contr(i,j-1) = u_boundary_contr(i,j-1) + &
4924 beta_lower(i,j) * dxdyh * 1./24 * (u_boundary_values(i-1,j-1) + &
4925 u_boundary_values(i-1,j) + u_boundary_values(i,j-1))
4927 v_boundary_contr(i,j-1) = v_boundary_contr(i,j-1) + &
4928 beta_lower(i,j) * dxdyh * 1./24 * (v_boundary_values(i-1,j-1) + &
4929 v_boundary_values(i-1,j) + v_boundary_values(i,j-1))
4932 if (umask(i-1,j) .eq. 1)
then 4934 u_boundary_contr(i-1,j) = u_boundary_contr(i-1,j) + &
4935 .5 * dxdyh * nu_lower(i,j) * ((4*ux+2*vy) * (0./dxh) + (uy+vy) * (1./dyh))
4937 v_boundary_contr(i-1,j) = v_boundary_contr(i-1,j) + &
4938 .5 * dxdyh * nu_lower(i,j) * ((uy+vx) * (0./dxh) + (4*vy+2*ux) * (1./dyh))
4940 u_boundary_contr(i,j-1) = u_boundary_contr(i,j-1) + &
4941 beta_lower(i,j) * dxdyh * 1./24 * (u_boundary_values(i-1,j-1) + &
4942 u_boundary_values(i-1,j) + u_boundary_values(i,j-1))
4944 v_boundary_contr(i,j-1) = v_boundary_contr(i,j-1) + &
4945 beta_lower(i,j) * dxdyh * 1./24 * (v_boundary_values(i-1,j-1) + &
4946 v_boundary_values(i-1,j) + v_boundary_values(i,j-1))
4949 if (umask(i-1,j-1) .eq. 1)
then 4951 u_boundary_contr(i-1,j-1) = u_boundary_contr(i-1,j-1) + &
4952 .5 * dxdyh * nu_upper(i,j) * ((4*ux+2*vy) * (-1./dxh) + (uy+vy) * (-1./dyh))
4954 v_boundary_contr(i-1,j-1) = v_boundary_contr(i-1,j-1) + &
4955 .5 * dxdyh * nu_upper(i,j) * ((uy+vx) * (-1./dxh) + (4*vy+2*ux) * (-1./dyh))
4957 u_boundary_contr(i-1,j-1) = u_boundary_contr(i-1,j-1) + &
4958 beta_lower(i,j) * dxdyh * 1./24 * (u_boundary_values(i-1,j-1) + &
4959 u_boundary_values(i-1,j) + u_boundary_values(i,j-1))
4961 v_boundary_contr(i-1,j-1) = v_boundary_contr(i-1,j-1) + &
4962 beta_lower(i,j) * dxdyh * 1./24 * (v_boundary_values(i-1,j-1) + &
4963 v_boundary_values(i-1,j) + v_boundary_values(i,j-1))
4968 if ((umask(i,j) .eq. 3) .OR. (umask(i,j-1) .eq. 3) .OR. (umask(i-1,j) .eq. 3))
then 4972 dxdyh = g%areaT(i,j)
4974 ux = (u_boundary_values(i,j)-u_boundary_values(i-1,j))/dxh
4975 vx = (v_boundary_values(i,j)-v_boundary_values(i-1,j))/dxh
4976 uy = (u_boundary_values(i,j)-u_boundary_values(i,j-1))/dyh
4977 vy = (v_boundary_values(i,j)-v_boundary_values(i,j-1))/dyh
4979 if (umask(i,j-1) .eq. 1)
then 4981 u_boundary_contr(i,j-1) = u_boundary_contr(i,j-1) + &
4982 .5 * dxdyh * nu_upper(i,j) * ((4*ux+2*vy) * (0./dxh) + (uy+vy) * (-1./dyh))
4984 v_boundary_contr(i,j-1) = v_boundary_contr(i,j-1) + &
4985 .5 * dxdyh * nu_upper(i,j) * ((uy+vx) * (0./dxh) + (4*vy+2*ux) * (-1./dyh))
4987 u_boundary_contr(i,j-1) = u_boundary_contr(i,j-1) + &
4988 beta_upper(i,j) * dxdyh * 1./24 * (u_boundary_values(i,j) + &
4989 u_boundary_values(i-1,j) + &
4990 u_boundary_values(i,j-1))
4992 v_boundary_contr(i,j-1) = v_boundary_contr(i,j-1) + &
4993 beta_upper(i,j) * dxdyh * 1./24 * (u_boundary_values(i,j) + &
4994 u_boundary_values(i-1,j) + &
4995 u_boundary_values(i,j-1))
4998 if (umask(i-1,j) .eq. 1)
then 5000 u_boundary_contr(i-1,j) = u_boundary_contr(i-1,j) + &
5001 .5 * dxdyh * nu_upper(i,j) * ((4*ux+2*vy) * (-1./dxh) + (uy+vy) * (0./dyh))
5003 v_boundary_contr(i-1,j) = v_boundary_contr(i-1,j) + &
5004 .5 * dxdyh * nu_upper(i,j) * ((uy+vx) * (-1./dxh) + (4*vy+2*ux) * (0./dyh))
5006 u_boundary_contr(i,j-1) = u_boundary_contr(i,j-1) + &
5007 beta_upper(i,j) * dxdyh * 1./24 * (u_boundary_values(i,j) + &
5008 u_boundary_values(i-1,j) + &
5009 u_boundary_values(i,j-1))
5011 v_boundary_contr(i,j-1) = v_boundary_contr(i,j-1) + &
5012 beta_upper(i,j) * dxdyh * 1./24 * (u_boundary_values(i,j) + &
5013 u_boundary_values(i-1,j) + &
5014 u_boundary_values(i,j-1))
5017 if (umask(i,j) .eq. 1)
then 5019 u_boundary_contr(i,j) = u_boundary_contr(i,j) + &
5020 .5 * dxdyh * nu_upper(i,j) * ((4*ux+2*vy) * (1./dxh) + (uy+vy) * (1./dyh))
5022 v_boundary_contr(i,j) = v_boundary_contr(i,j) + &
5023 .5 * dxdyh * nu_upper(i,j) * ((uy+vx) * (1./dxh) + (4*vy+2*ux) * (1./dyh))
5025 u_boundary_contr(i,j) = u_boundary_contr(i,j) + &
5026 beta_upper(i,j) * dxdyh * 1./24 * (u_boundary_values(i,j) + &
5027 u_boundary_values(i-1,j) + &
5028 u_boundary_values(i,j-1))
5030 v_boundary_contr(i,j) = v_boundary_contr(i,j) + &
5031 beta_upper(i,j) * dxdyh * 1./24 * (u_boundary_values(i,j) + &
5032 u_boundary_values(i-1,j) + &
5033 u_boundary_values(i,j-1))
5038 endif ;
enddo ;
enddo 5044 type(time_type),
intent(in) :: Time
5045 real,
dimension (:,:,:,:,:,:),
pointer:: Phisub
5047 real,
dimension (NILIMB_SYM_,NJLIMB_SYM_),
intent (in) :: H_node
5048 real,
dimension (:,:),
intent (in) :: float_cond
5050 real,
dimension (NILIMB_SYM_,NJLIMB_SYM_),
intent(inout) :: u_boundary_contr, v_boundary_contr
5055 real,
pointer,
dimension (:,:) :: u_boundary_values, &
5056 v_boundary_values, &
5059 real,
dimension(8,4) :: Phi
5060 real,
dimension(4) :: X, Y
5061 real,
dimension(2) :: xquad
5062 type(ocean_grid_type),
pointer :: G
5063 integer :: isym, i, j, isc, jsc, iec, jec, iq, jq, iphi, jphi, ilq, jlq
5064 real :: A, n, ux, uy, vx, vy, eps_min, domain_width, dxh, dyh, dxdyh, uq, vq, area, basel
5065 real,
dimension(2,2) :: Ucell,Vcell,Hcell,Usubcontr,Vsubcontr
5077 isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
5079 u_boundary_values => cs%u_boundary_values
5080 v_boundary_values => cs%v_boundary_values
5081 umask => cs%umask ; vmask => cs%vmask ; hmask => cs%hmask
5082 nu => cs%ice_visc_bilinear
5083 beta => cs%taub_beta_eff_bilinear
5085 xquad(1) = .5 * (1-sqrt(1./3)) ; xquad(2) = .5 * (1+sqrt(1./3))
5094 do j=jsc-1,jec+1 ;
do i=isc-1,iec+1 ;
if (hmask(i,j) .eq. 1)
then 5099 if ((umask(i-1,j-1) .eq. 3) .OR. (umask(i,j-1) .eq. 3) .OR. &
5100 (umask(i-1,j) .eq. 3) .OR. (umask(i,j) .eq. 3))
then 5105 dxdyh = g%areaT(i,j)
5107 x(1:2) = g%geoLonBu (i-1:i,j-1)*1000
5108 x(3:4) = g%geoLonBu (i-1:i,j)*1000
5109 y(1:2) = g%geoLatBu (i-1:i,j-1)*1000
5110 y(3:4) = g%geoLatBu (i-1:i,j)*1000
5123 do iq=1,2 ;
do jq=1,2
5125 uq = u_boundary_values(i-1,j-1) * xquad(3-iq) * xquad(3-jq) + &
5126 u_boundary_values(i,j-1) * xquad(iq) * xquad(3-jq) + &
5127 u_boundary_values(i-1,j) * xquad(3-iq) * xquad(jq) + &
5128 u_boundary_values(i,j) * xquad(iq) * xquad(jq)
5130 vq = v_boundary_values(i-1,j-1) * xquad(3-iq) * xquad(3-jq) + &
5131 v_boundary_values(i,j-1) * xquad(iq) * xquad(3-jq) + &
5132 v_boundary_values(i-1,j) * xquad(3-iq) * xquad(jq) + &
5133 v_boundary_values(i,j) * xquad(iq) * xquad(jq)
5135 ux = u_boundary_values(i-1,j-1) * phi(1,2*(jq-1)+iq) + &
5136 u_boundary_values(i,j-1) * phi(3,2*(jq-1)+iq) + &
5137 u_boundary_values(i-1,j) * phi(5,2*(jq-1)+iq) + &
5138 u_boundary_values(i,j) * phi(7,2*(jq-1)+iq)
5140 vx = v_boundary_values(i-1,j-1) * phi(1,2*(jq-1)+iq) + &
5141 v_boundary_values(i,j-1) * phi(3,2*(jq-1)+iq) + &
5142 v_boundary_values(i-1,j) * phi(5,2*(jq-1)+iq) + &
5143 v_boundary_values(i,j) * phi(7,2*(jq-1)+iq)
5145 uy = u_boundary_values(i-1,j-1) * phi(2,2*(jq-1)+iq) + &
5146 u_boundary_values(i,j-1) * phi(4,2*(jq-1)+iq) + &
5147 u_boundary_values(i-1,j) * phi(6,2*(jq-1)+iq) + &
5148 u_boundary_values(i,j) * phi(8,2*(jq-1)+iq)
5150 vy = v_boundary_values(i-1,j-1) * phi(2,2*(jq-1)+iq) + &
5151 v_boundary_values(i,j-1) * phi(4,2*(jq-1)+iq) + &
5152 v_boundary_values(i-1,j) * phi(6,2*(jq-1)+iq) + &
5153 v_boundary_values(i,j) * phi(8,2*(jq-1)+iq)
5155 do iphi=1,2 ;
do jphi=1,2
5157 if (iq .eq. iphi)
then 5163 if (jq .eq. jphi)
then 5169 if (umask(i-2+iphi,j-2+jphi) .eq. 1)
then 5172 u_boundary_contr(i-2+iphi,j-2+jphi) = u_boundary_contr(i-2+iphi,j-2+jphi) + &
5173 .25 * dxdyh * nu(i,j) * ( (4*ux+2*vy) * phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + &
5174 (uy+vx) * phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq) )
5176 if (float_cond(i,j) .eq. 0)
then 5177 u_boundary_contr(i-2+iphi,j-2+jphi) = u_boundary_contr(i-2+iphi,j-2+jphi) + &
5178 .25 * beta(i,j) * dxdyh * uq * xquad(ilq) * xquad(jlq)
5183 if (vmask(i-2+iphi,j-2+jphi) .eq. 1)
then 5186 v_boundary_contr(i-2+iphi,j-2+jphi) = v_boundary_contr(i-2+iphi,j-2+jphi) + &
5187 .25 * dxdyh * nu(i,j) * ( (uy+vx) * phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + &
5188 (4*vy+2*ux) * phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq))
5190 if (float_cond(i,j) .eq. 0)
then 5191 v_boundary_contr(i-2+iphi,j-2+jphi) = v_boundary_contr(i-2+iphi,j-2+jphi) + &
5192 .25 * beta(i,j) * dxdyh * vq * xquad(ilq) * xquad(jlq)
5199 if (float_cond(i,j) .eq. 1)
then 5200 usubcontr = 0.0 ; vsubcontr = 0.0 ; basel = g%bathyT(i,j)
5201 ucell(:,:) = u_boundary_values(i-1:i,j-1:j) ; vcell(:,:) = v_boundary_values(i-1:i,j-1:j)
5202 hcell(:,:) = h_node(i-1:i,j-1:j)
5204 (phisub, hcell, ucell, vcell, dxdyh, basel, dens_ratio, usubcontr, vsubcontr)
5205 do iphi=1,2 ;
do jphi = 1,2
5206 if (umask(i-2+iphi,j-2+jphi) .eq. 1)
then 5207 u_boundary_contr(i-2+iphi,j-2+jphi) = u_boundary_contr(i-2+iphi,j-2+jphi) + &
5208 usubcontr(iphi,jphi) * beta(i,j)
5210 if (vmask(i-2+iphi,j-2+jphi) .eq. 1)
then 5211 v_boundary_contr(i-2+iphi,j-2+jphi) = v_boundary_contr(i-2+iphi,j-2+jphi) + &
5212 vsubcontr(iphi,jphi) * beta(i,j)
5217 endif ;
enddo ;
enddo 5223 real,
dimension(:,:),
intent(inout) :: u, v
5232 real,
pointer,
dimension (:,:) :: nu_lower , &
5236 real,
pointer,
dimension (:,:) :: H, &
5239 type(ocean_grid_type),
pointer :: G
5240 integer :: isym, i, j, iscq, iecq, jscq, jecq, isd, jsd, ied, jed, iegq, jegq, giec, gjec, gisc, gjsc, cnt, isc, jsc, iec, jec, is, js
5241 real :: A, n, ux, uy, vx, vy, eps_min, umid, vmid, unorm, C_basal_friction, n_basal_friction, dxh, dyh, dxdyh
5245 if (g%symmetric)
then 5251 isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
5252 iscq = g%iscB ; iecq = g%iecB ; jscq = g%jscB ; jecq = g%jecB
5253 isd = g%isd ; jsd = g%jsd ; ied = g%isd ; jed = g%jsd
5254 iegq = g%iegB ; jegq = g%jegB
5255 gisc = g%domain%nihalo+1 ; gjsc = g%domain%njhalo+1
5256 giec = g%domain%niglobal+gisc ; gjec = g%domain%njglobal+gjsc
5257 is = iscq - (1-isym); js = jscq - (1-isym)
5259 a = cs%A_glen_isothermal ; n = cs%n_glen; eps_min = cs%eps_glen_min
5263 nu_upper => cs%ice_visc_upper_tri
5264 nu_lower => cs%ice_visc_lower_tri
5265 beta_eff_upper => cs%taub_beta_eff_upper_tri
5266 beta_eff_lower => cs%taub_beta_eff_lower_tri
5268 c_basal_friction = cs%C_basal_friction ; n_basal_friction = cs%n_basal_friction
5275 dxdyh = g%areaT(i,j)
5277 if (hmask(i,j) .eq. 1)
then 5278 ux = (u(i,j-1)-u(i-1,j-1)) / dxh
5279 vx = (v(i,j-1)-v(i-1,j-1)) / dxh
5280 uy = (u(i-1,j)-u(i-1,j-1)) / dyh
5281 vy = (v(i-1,j)-v(i-1,j-1)) / dyh
5283 nu_lower(i,j) = a**(-1/n) * (ux**2+vy**2+ux*vy+0.25*(uy+vx)**2+eps_min**2) ** ((1-n)/(2*n)) * h(i,j)
5284 umid = 1./3 * (u(i-1,j-1)+u(i-1,j)+u(i,j-1))
5285 vmid = 1./3 * (v(i-1,j-1)+v(i-1,j)+v(i,j-1))
5286 unorm = sqrt(umid**2+vmid**2+(eps_min*dxh)**2) ; beta_eff_lower(i,j) = c_basal_friction * unorm ** (n_basal_friction-1)
5288 ux = (u(i,j)-u(i-1,j)) / dxh
5289 vx = (v(i,j)-v(i-1,j)) / dxh
5290 uy = (u(i,j)-u(i,j-1)) / dyh
5291 vy = (u(i,j)-u(i,j-1)) / dyh
5293 nu_upper(i,j) = a**(-1/n) * (ux**2+vy**2+ux*vy+0.25*(uy+vx)**2+eps_min**2) ** ((1-n)/(2*n)) * h(i,j)
5294 umid = 1./3 * (u(i,j)+u(i-1,j)+u(i,j-1))
5295 vmid = 1./3 * (v(i,j)+v(i-1,j)+v(i,j-1))
5296 unorm = sqrt(umid**2+vmid**2+(eps_min*dxh)**2) ; beta_eff_upper(i,j) = c_basal_friction * unorm ** (n_basal_friction-1)
5306 real,
dimension(NILIMB_SYM_,NJLIMB_SYM_),
intent(inout) :: u, v
5315 real,
pointer,
dimension (:,:) :: nu, &
5317 real,
pointer,
dimension (:,:) :: H, &
5320 type(ocean_grid_type),
pointer :: G
5321 integer :: isym, i, j, iscq, iecq, jscq, jecq, isd, jsd, ied, jed, iegq, jegq, giec, gjec, gisc, gjsc, cnt, isc, jsc, iec, jec, is, js
5322 real :: A, n, ux, uy, vx, vy, eps_min, umid, vmid, unorm, C_basal_friction, n_basal_friction, dxh, dyh, dxdyh
5327 isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
5328 iscq = g%iscB ; iecq = g%iecB ; jscq = g%jscB ; jecq = g%jecB
5329 isd = g%isd ; jsd = g%jsd ; ied = g%ied ; jed = g%jed
5330 iegq = g%iegB ; jegq = g%jegB
5331 gisc = g%domain%nihalo+1 ; gjsc = g%domain%njhalo+1
5332 giec = g%domain%niglobal+gisc ; gjec = g%domain%njglobal+gjsc
5333 is = iscq - (1-isym); js = jscq - (1-isym)
5335 a = cs%A_glen_isothermal ; n = cs%n_glen; eps_min = cs%eps_glen_min
5336 c_basal_friction = cs%C_basal_friction ; n_basal_friction = cs%n_basal_friction
5340 nu => cs%ice_visc_bilinear
5341 beta => cs%taub_beta_eff_bilinear
5348 dxdyh = g%areaT(i,j)
5350 if (hmask(i,j) .eq. 1)
then 5351 ux = (u(i,j) + u(i,j-1) - u(i-1,j) - u(i-1,j-1)) / (2*dxh)
5352 vx = (v(i,j) + v(i,j-1) - v(i-1,j) - v(i-1,j-1)) / (2*dxh)
5353 uy = (u(i,j) - u(i,j-1) + u(i-1,j) - u(i-1,j-1)) / (2*dyh)
5354 vy = (v(i,j) - v(i,j-1) + v(i-1,j) - v(i-1,j-1)) / (2*dyh)
5356 nu(i,j) = .5 * a**(-1/n) * (ux**2+vy**2+ux*vy+0.25*(uy+vx)**2+eps_min**2) ** ((1-n)/(2*n)) * h(i,j)
5358 umid = (u(i,j) + u(i,j-1) + u(i-1,j) + u(i-1,j-1))/4
5359 vmid = (v(i,j) + v(i,j-1) + v(i-1,j) + v(i-1,j-1))/4
5360 unorm = sqrt(umid**2+vmid**2+(eps_min*dxh)**2) ; beta(i,j) = c_basal_friction * unorm ** (n_basal_friction-1)
5367 subroutine update_od_ffrac (CS, ocean_mass, counter, nstep_velocity, time_step, velocity_update_time_step)
5369 real,
dimension(CS%grid%isd:,CS%grid%jsd:) :: ocean_mass
5370 integer,
intent(in) :: counter
5371 integer,
intent(in) :: nstep_velocity
5372 real,
intent(in) :: time_step
5373 real,
intent(in) :: velocity_update_time_step
5375 type(ocean_grid_type),
pointer :: G
5376 integer :: isc, iec, jsc, jec, i, j
5377 real :: threshold_col_depth, rho_ocean, inv_rho_ocean
5379 threshold_col_depth = cs%thresh_float_col_depth
5383 rho_ocean = cs%density_ocean_avg
5384 inv_rho_ocean = 1./rho_ocean
5386 isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
5390 cs%OD_rt(i,j) = cs%OD_rt(i,j) + ocean_mass(i,j)*inv_rho_ocean
5391 if (ocean_mass(i,j) > threshold_col_depth*rho_ocean)
then 5392 cs%float_frac_rt(i,j) = cs%float_frac_rt(i,j) + 1.0
5397 if (counter .eq. nstep_velocity)
then 5401 cs%float_frac(i,j) = 1.0 - (cs%float_frac_rt(i,j) /
real(nstep_velocity))
5405 cs%OD_av(i,j) = cs%OD_rt(i,j) /
real(nstep_velocity)
5407 cs%OD_rt(i,j) = 0.0 ; cs%float_frac_rt(i,j) = 0.0
5411 call pass_var(cs%float_frac, g%domain)
5412 call pass_var(cs%OD_av, g%domain)
5421 type(ocean_grid_type),
pointer :: G
5422 integer :: i, j, iters, isd, ied, jsd, jed
5423 real :: rhoi, rhow, OD
5424 type(time_type) :: dummy_time
5425 real,
dimension(:,:),
pointer :: OD_av, float_frac, h_shelf
5429 rhoi = cs%density_ice
5430 rhow = cs%density_ocean_avg
5431 dummy_time = set_time(0,0)
5433 h_shelf => cs%h_shelf
5434 float_frac => cs%float_frac
5435 isd=g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
5441 od = g%bathyT(i,j) - rhoi/rhow * h_shelf(i,j)
5445 float_frac(i,j) = 0.
5448 float_frac(i,j) = 1.
5457 real,
dimension(4),
intent(in) :: X, Y
5458 real,
dimension(8,4),
intent (inout) :: Phi
5459 real,
intent (out) :: area
5478 real,
dimension (4) :: xquad, yquad
5479 integer :: node, qpoint, xnode, xq, ynode, yq
5480 real :: a,b,c,d,e,f,xexp,yexp
5482 xquad(1:3:2) = .5 * (1-sqrt(1./3)) ; yquad(1:2) = .5 * (1-sqrt(1./3))
5483 xquad(2:4:2) = .5 * (1+sqrt(1./3)) ; yquad(3:4) = .5 * (1+sqrt(1./3))
5487 a = -x(1)*(1-yquad(qpoint)) + x(2)*(1-yquad(qpoint)) - x(3)*yquad(qpoint) + x(4)*yquad(qpoint)
5488 b = -y(1)*(1-yquad(qpoint)) + y(2)*(1-yquad(qpoint)) - y(3)*yquad(qpoint) + y(4)*yquad(qpoint)
5489 c = -x(1)*(1-xquad(qpoint)) - x(2)*(xquad(qpoint)) + x(3)*(1-xquad(qpoint)) + x(4)*(xquad(qpoint))
5490 d = -y(1)*(1-xquad(qpoint)) - y(2)*(xquad(qpoint)) + y(3)*(1-xquad(qpoint)) + y(4)*(xquad(qpoint))
5494 xnode = 2-mod(node,2) ; ynode = ceiling(
REAL(node)/2)
5496 if (ynode .eq. 1)
then 5497 yexp = 1-yquad(qpoint)
5499 yexp = yquad(qpoint)
5502 if (1 .eq. xnode)
then 5503 xexp = 1-xquad(qpoint)
5505 xexp = xquad(qpoint)
5508 phi(2*node-1,qpoint) = ( d * (2 * xnode - 3) * yexp - b * (2 * ynode - 3) * xexp) / (a*d-b*c)
5509 phi(2*node,qpoint) = ( -c * (2 * xnode - 3) * yexp + a * (2 * ynode - 3) * xexp) / (a*d-b*c)
5520 real,
dimension(nsub,nsub,2,2,2,2),
intent(inout) :: Phisub
5548 integer :: i, j, k, l, qx, qy, indx, indy
5549 real,
dimension(2) :: xquad
5550 real :: x0, y0, x, y, val, fracx
5552 xquad(1) = .5 * (1-sqrt(1./3)) ; xquad(2) = .5 * (1+sqrt(1./3))
5553 fracx = 1.0/
real(nsub)
5557 x0 = (i-1) * fracx ; y0 = (j-1) * fracx
5560 x = x0 + fracx*xquad(qx)
5561 y = y0 + fracx*xquad(qy)
5575 phisub(i,j,k,l,qx,qy) = val
5597 integer :: isym, i, j, iscq, iecq, jscq, jecq, isd, jsd, is, js, iegq, jegq, giec, gjec, gisc, gjsc, isc, jsc, iec, jec, k
5598 integer :: i_off, j_off
5599 type(ocean_grid_type),
pointer :: G
5600 real,
dimension(:,:),
pointer :: umask, vmask, u_face_mask, v_face_mask, hmask, u_face_mask_boundary, v_face_mask_boundary
5603 isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
5604 iscq = g%iscB ; iecq = g%iecB ; jscq = g%jscB ; jecq = g%jecB
5605 i_off = g%idg_offset ; j_off = g%jdg_offset
5606 isd = g%isd ; jsd = g%jsd
5607 iegq = g%iegB ; jegq = g%jegB
5608 gisc = g%Domain%nihalo ; gjsc = g%Domain%njhalo
5609 giec = g%Domain%niglobal+gisc ; gjec = g%Domain%njglobal+gjsc
5613 u_face_mask => cs%u_face_mask
5614 v_face_mask => cs%v_face_mask
5615 u_face_mask_boundary => cs%u_face_mask_boundary
5616 v_face_mask_boundary => cs%v_face_mask_boundary
5628 umask(:,:) = 0 ; vmask(:,:) = 0
5629 u_face_mask(:,:) = 0 ; v_face_mask(:,:) = 0
5631 if (g%symmetric)
then 5634 is = isd+1 ; js = jsd+1
5640 if (hmask(i,j) .eq. 1)
then 5642 umask(i-1:i,j-1:j) = 1.
5643 vmask(i-1:i,j-1:j) = 1.
5647 select case (int(u_face_mask_boundary(i-1+k,j)))
5649 umask(i-1+k,j-1:j)=3.
5650 vmask(i-1+k,j-1:j)=0.
5651 u_face_mask(i-1+k,j)=3.
5653 u_face_mask(i-1+k,j)=2.
5655 umask(i-1+k,j-1:j)=0.
5656 vmask(i-1+k,j-1:j)=0.
5657 u_face_mask(i-1+k,j)=4.
5659 umask(i-1+k,j-1:j)=0.
5660 vmask(i-1+k,j-1:j)=0.
5661 u_face_mask(i-1+k,j)=0.
5663 umask(i-1+k,j-1:j)=0.
5670 select case (int(v_face_mask_boundary(i,j-1+k)))
5672 vmask(i-1:i,j-1+k)=3.
5673 umask(i-1:i,j-1+k)=0.
5674 v_face_mask(i,j-1+k)=3.
5676 v_face_mask(i,j-1+k)=2.
5678 umask(i-1:i,j-1+k)=0.
5679 vmask(i-1:i,j-1+k)=0.
5680 v_face_mask(i,j-1+k)=4.
5682 umask(i-1:i,j-1+k)=0.
5683 vmask(i-1:i,j-1+k)=0.
5684 u_face_mask(i,j-1+k)=0.
5686 vmask(i-1:i,j-1+k)=0.
5707 if (i .lt. g%ied)
then 5708 if ((hmask(i+1,j) .eq. 0) &
5709 .OR. (hmask(i+1,j) .eq. 2))
then 5711 u_face_mask(i,j) = 2.
5715 if (i .gt. g%isd)
then 5716 if ((hmask(i-1,j) .eq. 0) .OR. (hmask(i-1,j) .eq. 2))
then 5718 u_face_mask(i-1,j) = 2.
5722 if (j .gt. g%jsd)
then 5723 if ((hmask(i,j-1) .eq. 0) .OR. (hmask(i,j-1) .eq. 2))
then 5725 v_face_mask(i,j-1) = 2.
5729 if (j .lt. g%jed)
then 5730 if ((hmask(i,j+1) .eq. 0) .OR. (hmask(i,j+1) .eq. 2))
then 5732 v_face_mask(i,j) = 2.
5745 call pass_vector(u_face_mask, v_face_mask, g%domain, to_all, cgrid_ne)
5746 call pass_vector (umask,vmask,g%domain,to_all,bgrid_ne)
5753 real,
dimension (:,:),
intent(in) :: h_shelf, hmask
5754 real,
dimension (NILIMB_SYM_,NJLIMB_SYM_),
intent(inout) :: H_node
5756 type(ocean_grid_type),
pointer :: G
5757 integer :: i, j, isc, iec, jsc, jec, num_h, k, l
5761 isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec
5774 if (hmask(i+k,j+l) .eq. 1.0)
then 5775 summ = summ + h_shelf(i+k,j+l)
5780 if (num_h .gt. 0)
then 5781 h_node(i,j) = summ / num_h
5786 call pass_var(h_node, g%domain)
5794 if (.not.
associated(cs))
return 5796 deallocate(cs%mass_shelf) ;
deallocate(cs%area_shelf_h)
5797 deallocate(cs%t_flux) ;
deallocate(cs%lprec)
5798 deallocate(cs%salt_flux)
5800 deallocate(cs%tflux_shelf) ;
deallocate(cs%tfreeze);
5801 deallocate(cs%exch_vel_t) ;
deallocate(cs%exch_vel_s)
5803 deallocate(cs%h_shelf) ;
deallocate(cs%hmask)
5805 if (cs%shelf_mass_is_dynamic .and. .not.cs%override_shelf_movement)
then 5806 deallocate(cs%u_shelf) ;
deallocate(cs%v_shelf)
5808 deallocate(cs%t_shelf);
deallocate(cs%tmask);
5809 deallocate(cs%t_boundary_values)
5810 deallocate(cs%u_boundary_values) ;
deallocate(cs%v_boundary_values)
5811 deallocate(cs%ice_visc_bilinear)
5812 deallocate(cs%ice_visc_lower_tri) ;
deallocate(cs%ice_visc_upper_tri)
5813 deallocate(cs%u_face_mask) ;
deallocate(cs%v_face_mask)
5814 deallocate(cs%umask) ;
deallocate(cs%vmask)
5816 deallocate(cs%taub_beta_eff_bilinear)
5817 deallocate(cs%taub_beta_eff_upper_tri)
5818 deallocate(cs%taub_beta_eff_lower_tri)
5819 deallocate(cs%OD_rt) ;
deallocate(cs%OD_av)
5820 deallocate(cs%float_frac) ;
deallocate(cs%float_frac_rt)
5833 CHARACTER(*),
intent(in) :: fname
5836 REAL,
DIMENSION(:,:),
intent(in) :: A
5839 INTEGER :: M,N,i,j,iock,lh,FIN
5840 CHARACTER(23000) :: ln
5841 CHARACTER(17) :: sing
5843 CHARACTER(7) :: FMT1
5845 if (.NOT. flag)
then 5849 print *,
"WRITING ARRAY " // fname
5855 OPEN(unit=fin,file=fname,status=
'REPLACE',access=
'SEQUENTIAL',&
5856 action=
'WRITE',iostat=iock)
5858 IF(m .gt. 1300)
THEN 5859 WRITE(fin)
'SECOND DIMENSION TOO LARGE' 5865 WRITE(ln,
'(E17.9)') a(i,1)
5867 WRITE(sing,
'(E17.9)') a(i,j)
5868 ln = trim(ln) //
' ' // trim(sing)
5880 WRITE(fmt1(3:3),
'(I1)') lh
5883 WRITE(fmt1(3:4),
'(I2)') lh
5886 WRITE(fmt1(3:5),
'(I3)') lh
5889 WRITE(fmt1(3:6),
'(I4)') lh
5893 fmt1 = trim(fmt1) //
')' 5897 WRITE(unit=fin,iostat=iock,fmt=trim(fmt1)) trim(ln)
5899 IF(iock .ne. 0)
THEN 5909 subroutine solo_time_step (CS, time_step, n, Time, min_time_step_in)
5911 real,
intent(in) :: time_step
5912 integer,
intent(inout) :: n
5913 type(time_type) :: Time
5914 real,
optional,
intent(in) :: min_time_step_in
5916 type(ocean_grid_type),
pointer :: G
5917 integer :: is, iec, js, jec, i, j, ki, kj, iters
5918 real :: ratio, min_ratio, time_step_remain, local_u_max, &
5919 local_v_max, time_step_int, min_time_step,spy,dumtimeprint
5920 real,
dimension(:,:),
pointer :: u_shelf, v_shelf, hmask, umask, vmask
5922 type(time_type) :: dummy
5923 character(2) :: procnum
5924 character(4) :: stepnum
5926 cs%velocity_update_sub_counter = cs%velocity_update_sub_counter + 1
5929 u_shelf => cs%u_shelf
5930 v_shelf => cs%v_shelf
5934 time_step_remain = time_step
5935 if (.not. (
present (min_time_step_in)))
then 5936 min_time_step = 1000
5938 min_time_step=min_time_step_in
5940 is = g%isc ; iec = g%iec ; js = g%jsc ; jec = g%jec
5944 if (is_root_pe()) print *,
"TIME in ice shelf call, yrs: ", time_type_to_real(time)/spy
5946 do while (time_step_remain .gt. 0.0)
5953 local_u_max = 0 ; local_v_max = 0
5955 if (hmask(i,j) .eq. 1.0)
then 5958 do ki=1,2 ;
do kj = 1,2
5960 local_u_max = max(local_u_max, abs(u_shelf(i-1+ki,j-1+kj)))
5961 local_v_max = max(local_v_max, abs(v_shelf(i-1+ki,j-1+kj)))
5965 ratio = min(g%areaT(i,j) / (local_u_max+1.0e-12), g%areaT(i,j) / (local_v_max+1.0e-12))
5966 min_ratio = min(min_ratio, ratio)
5974 call mpp_min (min_ratio)
5976 time_step_int = min(cs%CFL_factor * min_ratio * (365*86400), time_step)
5978 if (time_step_int .lt. min_time_step)
then 5979 call mom_error (fatal,
"MOM_ice_shelf:solo_time_step: abnormally small timestep")
5981 if (is_root_pe())
then 5982 write(*,*)
"Ice model timestep: ", time_step_int,
" seconds" 5986 if (time_step_int .ge. time_step_remain)
then 5987 time_step_int = time_step_remain
5988 time_step_remain = 0.0
5990 time_step_remain = time_step_remain - time_step_int
5993 write (stepnum,
'(I4)') cs%velocity_update_sub_counter
5997 if (mpp_pe() .eq. 7)
then 5998 call savearray2 (
"hmask",cs%hmask,cs%write_output_to_file)
6005 if (time_step_int .gt. 1000)
then 6018 call enable_averaging(time_step,time,cs%diag)
6019 if (cs%id_area_shelf_h > 0)
call post_data(cs%id_area_shelf_h, cs%area_shelf_h, cs%diag)
6020 if (cs%id_col_thick > 0)
call post_data(cs%id_col_thick, cs%OD_av, cs%diag)
6021 if (cs%id_h_shelf > 0)
call post_data(cs%id_h_shelf,cs%h_shelf,cs%diag)
6022 if (cs%id_h_mask > 0)
call post_data(cs%id_h_mask,cs%hmask,cs%diag)
6023 if (cs%id_u_mask > 0)
call post_data(cs%id_u_mask,cs%umask,cs%diag)
6024 if (cs%id_v_mask > 0)
call post_data(cs%id_v_mask,cs%vmask,cs%diag)
6025 if (cs%id_u_shelf > 0)
call post_data(cs%id_u_shelf,cs%u_shelf,cs%diag)
6026 if (cs%id_v_shelf > 0)
call post_data(cs%id_v_shelf,cs%v_shelf,cs%diag)
6027 if (cs%id_float_frac > 0)
call post_data(cs%id_float_frac,cs%float_frac,cs%diag)
6028 if (cs%id_OD_av >0)
call post_data(cs%id_OD_av,cs%OD_av,cs%diag)
6029 if (cs%id_float_frac_rt>0)
call post_data(cs%id_float_frac_rt,cs%float_frac_rt,cs%diag)
6032 call post_data(cs%id_t_mask,cs%tmask,cs%diag)
6034 call post_data(cs%id_t_shelf,cs%t_shelf,cs%diag)
6036 call disable_averaging(cs%diag)
6045 real,
intent(in) :: time_step
6046 real,
pointer,
dimension(:,:),
intent(in) :: melt_rate
6047 type(time_type) :: Time
6086 type(ocean_grid_type),
pointer :: G
6087 real,
dimension(size(CS%h_shelf,1),size(CS%h_shelf,2)) :: th_after_uflux, th_after_vflux, TH
6088 real,
dimension(size(CS%h_shelf,1),size(CS%h_shelf,2),4) :: flux_enter
6089 integer :: isd, ied, jsd, jed, i, j, isc, iec, jsc, jec
6090 real :: rho, spy, t_bd, Tsurf, adot
6091 real,
dimension(:,:),
pointer :: hmask, Tbot
6092 character(len=2) :: procnum
6096 rho = cs%density_ice
6103 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
6104 isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
6105 flux_enter(:,:,:) = 0.0
6107 th_after_uflux(:,:) = 0.0
6108 th_after_vflux(:,:) = 0.0
6112 t_bd = cs%t_boundary_values(i,j)
6114 if ((cs%hmask(i,j) .eq. 3) .or. (cs%hmask(i,j) .eq. -2))
then 6115 cs%t_shelf(i,j) = cs%t_boundary_values(i,j)
6122 th(i,j) = cs%t_shelf(i,j)*cs%h_shelf (i,j)
6146 if (cs%h_shelf(i,j) .gt. 0.0)
then 6147 cs%t_shelf (i,j) = th_after_vflux(i,j)/cs%h_shelf (i,j)
6149 cs%t_shelf(i,j) = -10.0
6156 t_bd = cs%t_boundary_values(i,j)
6158 if ((cs%hmask(i,j) .eq. 3) .or. (cs%hmask(i,j) .eq. -2))
then 6159 cs%t_shelf(i,j) = t_bd
6167 if ((cs%hmask(i,j) .eq. 1) .or. (cs%hmask(i,j) .eq. 2))
then 6168 if (cs%h_shelf(i,j) .gt. 0.0)
then 6170 cs%t_shelf (i,j) = cs%t_shelf (i,j) + time_step*(adot*tsurf -3/spy*tbot(i,j))/cs%h_shelf (i,j)
6177 cs%t_shelf(i,j) = -10.0
6184 call pass_var(cs%t_shelf, g%domain)
6185 call pass_var(cs%tmask, g%domain)
6188 call hchksum (cs%t_shelf,
"temp after front", g%HI, haloshift=3)
6196 real,
intent(in) :: time_step
6197 real,
dimension(:,:),
intent(in) :: h0
6198 real,
dimension(:,:),
intent(inout) :: h_after_uflux
6199 real,
dimension(:,:,:),
intent(inout) :: flux_enter
6219 integer :: isym, i, j, is, ie, js, je, isd, ied, jsd, jed, gjed, gied
6220 integer :: i_off, j_off
6221 logical :: at_east_bdry, at_west_bdry, one_off_west_bdry, one_off_east_bdry
6222 type(ocean_grid_type),
pointer :: G
6223 real,
dimension(-2:2) :: stencil
6224 real,
dimension(:,:),
pointer :: hmask, u_face_mask, u_flux_boundary_values,u_boundary_values,t_boundary
6226 flux_diff_cell, phi, dxh, dyh, dxdyh
6228 character (len=1) :: debug_str, procnum
6240 u_face_mask => cs%u_face_mask
6241 u_flux_boundary_values => cs%u_flux_boundary_values
6242 u_boundary_values => cs%u_shelf
6244 t_boundary => cs%t_boundary_values
6245 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
6246 i_off = g%idg_offset ; j_off = g%jdg_offset
6249 if (((j+j_off) .le. g%domain%njglobal+g%domain%njhalo) .AND. &
6250 ((j+j_off) .ge. g%domain%njhalo+1))
then 6256 if (((i+i_off) .le. g%domain%niglobal+g%domain%nihalo) .AND. &
6257 ((i+i_off) .ge. g%domain%nihalo+1))
then 6259 if (i+i_off .eq. g%domain%nihalo+1)
then 6262 at_west_bdry=.false.
6265 if (i+i_off .eq. g%domain%niglobal+g%domain%nihalo)
then 6268 at_east_bdry=.false.
6271 if (hmask(i,j) .eq. 1)
then 6273 dxh = g%dxT(i,j) ; dyh = g%dyT(i,j) ; dxdyh = g%areaT(i,j)
6275 h_after_uflux(i,j) = h0(i,j)
6277 stencil(:) = h0(i-2:i+2,j)
6283 if (u_face_mask(i-1,j) .eq. 4.)
then 6285 flux_diff_cell = flux_diff_cell + dyh * time_step * u_flux_boundary_values(i-1,j) * &
6286 t_boundary(i-1,j) / dxdyh
6293 u_face = 0.5 * (cs%u_shelf(i-1,j-1) + cs%u_shelf(i-1,j))
6299 if (u_face .gt. 0)
then 6303 if (at_west_bdry .AND. (hmask(i-1,j).eq.3))
then 6305 stencil(-1) = cs%t_boundary_values(i-1,j)*cs%h_shelf(i-1,j)
6306 flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step * stencil(-1) / dxdyh
6308 elseif (hmask(i-1,j) * hmask(i-2,j) .eq. 1)
then 6309 phi =
slope_limiter(stencil(-1)-stencil(-2), stencil(0)-stencil(-1))
6310 flux_diff_cell = flux_diff_cell + abs(u_face) * dyh* time_step / dxdyh * &
6311 (stencil(-1) - phi * (stencil(-1)-stencil(0))/2)
6317 flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step / dxdyh * stencil(-1)
6321 elseif (u_face .lt. 0)
then 6322 if (hmask(i-1,j) * hmask(i+1,j) .eq. 1)
then 6323 phi =
slope_limiter(stencil(0)-stencil(1), stencil(-1)-stencil(0))
6324 flux_diff_cell = flux_diff_cell - abs(u_face) * dyh * time_step / dxdyh * &
6325 (stencil(0) - phi * (stencil(0)-stencil(-1))/2)
6328 flux_diff_cell = flux_diff_cell - abs(u_face) * dyh * time_step / dxdyh * stencil(0)
6330 if ((hmask(i-1,j) .eq. 0) .OR. (hmask(i-1,j) .eq. 2))
then 6331 flux_enter(i-1,j,2) = abs(u_face) * dyh * time_step * stencil(0)
6341 if (u_face_mask(i+1,j) .eq. 4.)
then 6343 flux_diff_cell = flux_diff_cell + dyh * time_step * u_flux_boundary_values(i+1,j) *&
6344 t_boundary(i+1,j)/ dxdyh
6350 u_face = 0.5 * (cs%u_shelf(i,j-1) + cs%u_shelf(i,j))
6352 if (u_face .lt. 0)
then 6354 if (at_east_bdry .AND. (hmask(i+1,j).eq.3))
then 6357 flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step * stencil(1) / dxdyh
6359 elseif (hmask(i+1,j) * hmask(i+2,j) .eq. 1)
then 6361 phi =
slope_limiter(stencil(1)-stencil(2), stencil(0)-stencil(1))
6362 flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step / dxdyh * &
6363 (stencil(1) - phi * (stencil(1)-stencil(0))/2)
6369 flux_diff_cell = flux_diff_cell + abs(u_face) * dyh * time_step / dxdyh * stencil(1)
6373 elseif (u_face .gt. 0)
then 6375 if (hmask(i-1,j) * hmask(i+1,j) .eq. 1)
then 6377 phi =
slope_limiter(stencil(0)-stencil(-1), stencil(1)-stencil(0))
6378 flux_diff_cell = flux_diff_cell - abs(u_face) * dyh * time_step / dxdyh * &
6379 (stencil(0) - phi * (stencil(0)-stencil(1))/2)
6385 flux_diff_cell = flux_diff_cell - abs(u_face) * dyh * time_step / dxdyh * stencil(0)
6387 if ((hmask(i+1,j) .eq. 0) .OR. (hmask(i+1,j) .eq. 2))
then 6389 flux_enter(i+1,j,1) = abs(u_face) * dyh * time_step * stencil(0)
6396 h_after_uflux(i,j) = h_after_uflux(i,j) + flux_diff_cell
6400 elseif ((hmask(i,j) .eq. 0) .OR. (hmask(i,j) .eq. 2))
then 6402 if (at_west_bdry .AND. (hmask(i-1,j) .EQ. 3))
then 6403 u_face = 0.5 * (cs%u_shelf(i-1,j-1) + cs%u_shelf(i-1,j))
6404 flux_enter(i,j,1) = abs(u_face) * g%dyT(i,j) * time_step * t_boundary(i-1,j)*cs%thickness_boundary_values(i+1,j)
6405 elseif (u_face_mask(i-1,j) .eq. 4.)
then 6406 flux_enter(i,j,1) = g%dyT(i,j) * time_step * u_flux_boundary_values(i-1,j)*t_boundary(i-1,j)
6411 if (at_east_bdry .AND. (hmask(i+1,j) .EQ. 3))
then 6412 u_face = 0.5 * (cs%u_shelf(i,j-1) + cs%u_shelf(i,j))
6413 flux_enter(i,j,2) = abs(u_face) * g%dyT(i,j) * time_step * t_boundary(i+1,j)*cs%thickness_boundary_values(i+1,j)
6414 elseif (u_face_mask(i+1,j) .eq. 4.)
then 6415 flux_enter(i,j,2) = g%dyT(i,j) * time_step * u_flux_boundary_values(i+1,j) * t_boundary(i+1,j)
6449 real,
intent(in) :: time_step
6450 real,
dimension(:,:),
intent(in) :: h_after_uflux
6451 real,
dimension(:,:),
intent(inout) :: h_after_vflux
6452 real,
dimension(:,:,:),
intent(inout) :: flux_enter
6472 integer :: isym, i, j, is, ie, js, je, isd, ied, jsd, jed, gjed, gied
6473 integer :: i_off, j_off
6474 logical :: at_north_bdry, at_south_bdry, one_off_west_bdry, one_off_east_bdry
6475 type(ocean_grid_type),
pointer :: G
6476 real,
dimension(-2:2) :: stencil
6477 real,
dimension(:,:),
pointer :: hmask, v_face_mask, v_flux_boundary_values,t_boundary,v_boundary_values
6479 flux_diff_cell, phi, dxh, dyh, dxdyh
6480 character(len=1) :: debug_str, procnum
6492 v_face_mask => cs%v_face_mask
6493 v_flux_boundary_values => cs%v_flux_boundary_values
6494 t_boundary => cs%t_boundary_values
6495 v_boundary_values => cs%v_shelf
6496 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
6497 i_off = g%idg_offset ; j_off = g%jdg_offset
6500 if (((i+i_off) .le. g%domain%niglobal+g%domain%nihalo) .AND. &
6501 ((i+i_off) .ge. g%domain%nihalo+1))
then 6507 if (((j+j_off) .le. g%domain%njglobal+g%domain%njhalo) .AND. &
6508 ((j+j_off) .ge. g%domain%njhalo+1))
then 6510 if (j+j_off .eq. g%domain%njhalo+1)
then 6511 at_south_bdry=.true.
6513 at_south_bdry=.false.
6515 if (j+j_off .eq. g%domain%njglobal+g%domain%njhalo)
then 6516 at_north_bdry=.true.
6518 at_north_bdry=.false.
6521 if (hmask(i,j) .eq. 1)
then 6522 dxh = g%dxT(i,j) ; dyh = g%dyT(i,j) ; dxdyh = g%areaT(i,j)
6523 h_after_vflux(i,j) = h_after_uflux(i,j)
6525 stencil(:) = h_after_uflux(i,j-2:j+2)
6530 if (v_face_mask(i,j-1) .eq. 4.)
then 6532 flux_diff_cell = flux_diff_cell + dxh * time_step * v_flux_boundary_values(i,j-1) * t_boundary(i,j-1)/ dxdyh
6539 v_face = 0.5 * (cs%v_shelf(i-1,j-1) + cs%v_shelf(i,j-1))
6541 if (v_face .gt. 0)
then 6545 if (at_south_bdry .AND. (hmask(i,j-1).eq.3))
then 6547 flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step * stencil(-1) / dxdyh
6549 elseif (hmask(i,j-1) * hmask(i,j-2) .eq. 1)
then 6551 phi =
slope_limiter(stencil(-1)-stencil(-2), stencil(0)-stencil(-1))
6552 flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step / dxdyh * &
6553 (stencil(-1) - phi * (stencil(-1)-stencil(0))/2)
6558 flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step / dxdyh * stencil(-1)
6561 elseif (v_face .lt. 0)
then 6563 if (hmask(i,j-1) * hmask(i,j+1) .eq. 1)
then 6564 phi =
slope_limiter(stencil(0)-stencil(1), stencil(-1)-stencil(0))
6565 flux_diff_cell = flux_diff_cell - abs(v_face) * dxh * time_step / dxdyh * &
6566 (stencil(0) - phi * (stencil(0)-stencil(-1))/2)
6568 flux_diff_cell = flux_diff_cell - abs(v_face) * dxh * time_step / dxdyh * stencil(0)
6570 if ((hmask(i,j-1) .eq. 0) .OR. (hmask(i,j-1) .eq. 2))
then 6571 flux_enter(i,j-1,4) = abs(v_face) * dyh * time_step * stencil(0)
6582 if (v_face_mask(i,j+1) .eq. 4.)
then 6584 flux_diff_cell = flux_diff_cell + dxh * time_step * v_flux_boundary_values(i,j+1) *&
6585 t_boundary(i,j+1)/ dxdyh
6592 v_face = 0.5 * (cs%v_shelf(i-1,j) + cs%v_shelf(i,j))
6594 if (v_face .lt. 0)
then 6596 if (at_north_bdry .AND. (hmask(i,j+1).eq.3))
then 6598 flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step * stencil(1) / dxdyh
6599 elseif (hmask(i,j+1) * hmask(i,j+2) .eq. 1)
then 6600 phi =
slope_limiter(stencil(1)-stencil(2), stencil(0)-stencil(1))
6601 flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step / dxdyh * &
6602 (stencil(1) - phi * (stencil(1)-stencil(0))/2)
6606 flux_diff_cell = flux_diff_cell + abs(v_face) * dxh * time_step / dxdyh * stencil(1)
6609 elseif (v_face .gt. 0)
then 6611 if (hmask(i,j-1) * hmask(i,j+1) .eq. 1)
then 6612 phi =
slope_limiter(stencil(0)-stencil(-1), stencil(1)-stencil(0))
6613 flux_diff_cell = flux_diff_cell - abs(v_face) * dxh * time_step / dxdyh * &
6614 (stencil(0) - phi * (stencil(0)-stencil(1))/2)
6618 flux_diff_cell = flux_diff_cell - abs(v_face) * dxh * time_step / dxdyh * stencil(0)
6619 if ((hmask(i,j+1) .eq. 0) .OR. (hmask(i,j+1) .eq. 2))
then 6620 flux_enter(i,j+1,3) = abs(v_face) * dxh * time_step * stencil(0)
6628 h_after_vflux(i,j) = h_after_vflux(i,j) + flux_diff_cell
6630 elseif ((hmask(i,j) .eq. 0) .OR. (hmask(i,j) .eq. 2))
then 6632 if (at_south_bdry .AND. (hmask(i,j-1) .EQ. 3))
then 6633 v_face = 0.5 * (cs%v_shelf(i-1,j-1) + cs%v_shelf(i,j-1))
6634 flux_enter(i,j,3) = abs(v_face) * g%dxT(i,j) * time_step * t_boundary(i,j-1)*cs%thickness_boundary_values(i,j-1)
6635 elseif (v_face_mask(i,j-1) .eq. 4.)
then 6636 flux_enter(i,j,3) = g%dxT(i,j) * time_step * v_flux_boundary_values(i,j-1)*t_boundary(i,j-1)
6642 if (at_north_bdry .AND. (hmask(i,j+1) .EQ. 3))
then 6643 v_face = 0.5 * (cs%v_shelf(i-1,j) + cs%v_shelf(i,j))
6644 flux_enter(i,j,4) = abs(v_face) * g%dxT(i,j) * time_step * t_boundary(i,j+1)*cs%thickness_boundary_values(i,j+1)
6645 elseif (v_face_mask(i,j+1) .eq. 4.)
then 6646 flux_enter(i,j,4) = g%dxT(i,j) * time_step * v_flux_boundary_values(i,j+1)*t_boundary(i,j+1)
The following structure contains pointers to various fields which may be used describe the surface st...
subroutine, public user_initialize_shelf_mass(mass_shelf, area_shelf_h, h_shelf, hmask, G, CS, param_file, new_sim)
integer id_clock_pass
Clock for group pass calls.
subroutine, public mom_grid_init(G, param_file, HI, global_indexing, bathymetry_at_vel)
MOM_grid_init initializes the ocean grid array sizes and grid memory.
type(vardesc) function, public var_desc(name, units, longname, hor_grid, z_grid, t_grid, cmor_field_name, cmor_units, conversion, caller)
Returns a vardesc type whose elements have been filled with the provided fields. The argument name is...
This module implements boundary forcing for MOM6.
Initializes fixed aspects of the model, such as horizontal grid metrics, topography and Coriolis...
integer, parameter, public to_all
subroutine calve_to_mask(CS, h_shelf, area_shelf_h, hmask, calve_mask)
Ocean grid type. See mom_grid for details.
subroutine update_shelf_mass(G, CS, Time, fluxes)
Updates the ice shelf mass using data from a file.
Calculates density of sea water from T, S and P.
Calculates the freezing point of sea water from T, S and P.
subroutine ice_shelf_advect_thickness_x(CS, time_step, h0, h_after_uflux, flux_enter)
subroutine, public mom_initialize_topography(D, max_depth, G, PF)
MOM_initialize_topography makes the appropriate call to set up the bathymetry.
Provides the ocean grid type.
subroutine interpolate_h_to_b(CS, h_shelf, hmask, H_node)
subroutine bilinear_shape_functions(X, Y, Phi, area)
subroutine, public mom_forcing_chksum(mesg, fluxes, G, haloshift)
Write out chksums for basic state variables.
subroutine, public allocate_forcing_type(G, fluxes, stress, ustar, water, heat, shelf, press, iceberg)
Conditionally allocate fields within the forcing type.
subroutine shelf_advance_front(CS, flux_enter)
subroutine initialize_shelf_mass(G, param_file, CS, new_sim)
Initializes shelf mass based on three options (file, zero and user)
subroutine, public add_shelf_flux(G, CS, state, fluxes)
Updates suface fluxes that are influenced by sub-ice-shelf melting.
subroutine, public eos_init(param_file, EOS)
Initializes EOS_type by allocating and reading parameters.
This module contains I/O framework code.
subroutine apply_boundary_values_triangle(CS, time, u_boundary_contr, v_boundary_contr)
subroutine ice_shelf_advect_thickness_y(CS, time_step, h_after_uflux, h_after_vflux, flux_enter)
subroutine, public user_update_shelf_mass(mass_shelf, area_shelf_h, h_shelf, hmask, G, CS, Time, new_sim)
subroutine bilinear_shape_functions_subgrid(Phisub, nsub)
By Robert Hallberg, April 1994 - June 2002 *This subroutine initializes the fields for the simulation...
subroutine, public ice_shelf_end(CS)
Deallocates all memory associated with this module.
subroutine ice_shelf_solve_inner(CS, u, v, taudx, taudy, H_node, float_cond, FE, conv_flag, iters, time, Phi, Phisub)
subroutine initialize_diagnostic_fields(CS, FE, Time)
subroutine, public destroy_dyn_horgrid(G)
Release memory used by the dyn_horgrid_type and related structures.
subroutine matrix_diagonal_bilinear(CS, float_cond, H_node, dens_ratio, Phisub, u_diagonal, v_diagonal)
subroutine, public copy_dyngrid_to_mom_grid(dG, oG)
Copies information from a dynamic (shared) horizontal grid type into an ocean_grid_type.
subroutine savearray2(fname, A, flag)
subroutine cg_action_triangular(uret, vret, u, v, umask, vmask, hmask, nu_upper, nu_lower, beta_upper, beta_lower, dxh, dyh, dxdyh, is, ie, js, je, isym)
Implements the thermodynamic aspects of ocean / ice-shelf interactions,.
subroutine, public solo_time_step(CS, time_step, n, Time, min_time_step_in)
subroutine update_od_ffrac_uncoupled(CS)
subroutine, public initialize_ice_thickness(h_shelf, area_shelf_h, hmask, G, PF)
subroutine, public calculate_density_derivs(T, S, pressure, drho_dT, drho_dS, start, npts, EOS)
Calls the appropriate subroutine to calculate density derivatives for 1-D array inputs.
subroutine calc_shelf_visc_bilinear(CS, u, v)
subroutine ice_shelf_advect_temp_y(CS, time_step, h_after_uflux, h_after_vflux, flux_enter)
subroutine, public ice_shelf_save_restart(CS, Time, directory, time_stamped, filename_suffix)
Save the ice shelf restart file.
subroutine cg_action_subgrid_basal_bilinear(Phisub, H, U, V, DXDYH, D, dens_ratio, Ucontr, Vcontr, iin, jin)
subroutine, public mom_domains_init(MOM_dom, param_file, symmetric, static_memory, NIHALO, NJHALO, NIGLOBAL, NJGLOBAL, NIPROC, NJPROC, min_halo, domain_name, include_name, param_suffix)
logical function, public is_root_pe()
subroutine cg_action_bilinear(uret, vret, u, v, Phi, Phisub, umask, vmask, hmask, H_node, nu, float_cond, D, beta, dxdyh, is, ie, js, je, dens_ratio)
subroutine, public set_grid_metrics(G, param_file)
set_grid_metrics is used to set the primary values in the model's horizontal grid. The bathymetry, land-sea mask and any restricted channel widths are not known yet, so these are set later.
subroutine, public initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Time_in, solo_ice_sheet_in)
Initializes shelf model data, parameters and diagnostics.
Structure that contains pointers to the boundary forcing used to drive the liquid ocean simulated by ...
subroutine cg_diagonal_subgrid_basal_bilinear(Phisub, H, DXDYH, D, dens_ratio, Ucontr, Vcontr)
subroutine update_od_ffrac(CS, ocean_mass, counter, nstep_velocity, time_step, velocity_update_time_step)
subroutine, public mom_mesg(message, verb, all_print)
Provides subroutines for quantities specific to the equation of state.
subroutine ice_shelf_advect_temp_x(CS, time_step, h0, h_after_uflux, flux_enter)
subroutine ice_shelf_min_thickness_calve(CS, h_shelf, area_shelf_h, hmask)
Apply a very simple calving law using a minimum thickness rule.
Type for describing a variable, typically a tracer.
subroutine calc_shelf_driving_stress(CS, TAUD_X, TAUD_Y, OD, FE)
subroutine apply_boundary_values_bilinear(CS, time, Phisub, H_node, float_cond, dens_ratio, u_boundary_contr, v_boundary_contr)
Control structure that contains ice shelf parameters and diagnostics handles.
subroutine matrix_diagonal_triangle(CS, u_diagonal, v_diagonal)
subroutine change_thickness_using_melt(CS, G, time_step, fluxes)
Changes the thickness (mass) of the ice shelf based on sub-ice-shelf melting.
subroutine, public save_restart(directory, time, G, CS, time_stamped, filename, GV)
subroutine, public copy_mom_grid_to_dyngrid(oG, dG)
Copies information from an ocean_grid_type into a dynamic (shared) horizontal grid type...
subroutine ice_shelf_solve_outer(CS, u, v, FE, iters, time)
subroutine ice_shelf_temp(CS, time_step, melt_rate, Time)
subroutine, public create_dyn_horgrid(G, HI, bathymetry_at_vel)
Allocate memory used by the dyn_horgrid_type and related structures.
real function slope_limiter(num, denom)
used for flux limiting in advective subroutines Van Leer limiter (source: Wikipedia) ...
real function quad_area(X, Y)
Calculate area of quadrilateral.
subroutine, public restart_init(param_file, CS, restart_root)
subroutine update_velocity_masks(CS)
subroutine, public mom_error(level, message, all_print)
subroutine init_boundary_values(CS, time, input_flux, input_thick, new_sim)
subroutine, public user_initialize_topography(D, G, param_file, max_depth)
Initialize topography.
subroutine calc_shelf_visc_triangular(CS, u, v)
subroutine, public shelf_calc_flux(state, fluxes, Time, time_step, CS)
Calculates fluxes between the ocean and ice-shelf using the three-equations formulation (optional to ...
subroutine, public restore_state(filename, directory, day, G, CS)
A control structure for the equation of state.
subroutine ice_shelf_advect(CS, time_step, melt_rate, Time)