From f1041d1f90c2380c5543b75a65ea53c4b57a5060 Mon Sep 17 00:00:00 2001 From: Olga Sergienko Date: Tue, 15 Dec 2020 20:32:06 -0500 Subject: [PATCH 01/17] Fixed bugs in CG_action, matrix_diagonal and calc_shelf_visc in MOM_ice_shelf_dynamics.F90 modified initialize_ice_shelf_boundary_channel in MOM_ice_shelf_initialze.F90 --- src/ice_shelf/MOM_ice_shelf_dynamics.F90 | 212 ++++++++++++++++----- src/ice_shelf/MOM_ice_shelf_initialize.F90 | 125 ++++++++++++ 2 files changed, 286 insertions(+), 51 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index f038190753..8480906de8 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -23,6 +23,7 @@ module MOM_ice_shelf_dynamics use MOM_ice_shelf_state, only : ice_shelf_state use MOM_coms, only : reproducing_sum, sum_across_PEs, max_across_PEs, min_across_PEs use MOM_checksums, only : hchksum, qchksum +use MOM_ice_shelf_initialize, only : initialize_ice_shelf_boundary_channel !OVS intializing b.c.s implicit none ; private @@ -366,20 +367,23 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ call get_param(param_file, mdl, "A_GLEN_ISOTHERM", CS%A_glen_isothermal, & "Ice viscosity parameter in Glen's Law", & - units="Pa-3 yr-1", default=9.461e-18, scale=1.0/(365.0*86400.0)) + units="Pa-3 s-1", default=2.2261e-25, scale=1.0) !OVS change units to Pa-3 s-1 +! units="Pa-3 yr-1", default=9.461e-18, scale=1.0/(365.0*86400.0)) ! This default is equivalent to 3.0001e-25 Pa-3 s-1, appropriate at about -10 C. call get_param(param_file, mdl, "GLEN_EXPONENT", CS%n_glen, & "nonlinearity exponent in Glen's Law", & units="none", default=3.) call get_param(param_file, mdl, "MIN_STRAIN_RATE_GLEN", CS%eps_glen_min, & "min. strain rate to avoid infinite Glen's law viscosity", & - units="a-1", default=1.e-12, scale=US%T_to_s/(365.0*86400.0)) + units="s-1", default=1.e-19, scale=US%T_to_s) !OVS change units to s-1 + !units="a-1", default=1.e-12, scale=US%T_to_s/(365.0*86400.0)) call get_param(param_file, mdl, "BASAL_FRICTION_EXP", CS%n_basal_fric, & "Exponent in sliding law \tau_b = C u^(n_basal_fric)", & units="none", fail_if_missing=.true.) call get_param(param_file, mdl, "BASAL_FRICTION_COEFF", CS%C_basal_friction, & "Coefficient in sliding law \tau_b = C u^(n_basal_fric)", & - units="Pa (m yr-1)-(n_basal_fric)", scale=US%kg_m2s_to_RZ_T*((365.0*86400.0)**CS%n_basal_fric), & + units="Pa (m s-1)^(n_basal_fric)", scale=US%kg_m2s_to_RZ_T**CS%n_basal_fric, & ! OVS change units to s-1 + !units="Pa (m yr-1)-(n_basal_fric)", scale=US%kg_m2s_to_RZ_T*((365.0*86400.0)**CS%n_basal_fric), & fail_if_missing=.true.) call get_param(param_file, mdl, "DENSITY_ICE", CS%density_ice, & "A typical density of ice.", units="kg m-3", default=917.0, scale=US%kg_m3_to_R) @@ -399,10 +403,11 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ call get_param(param_file, mdl, "SHELF_MOVING_FRONT", CS%moving_shelf_front, & "Specify whether to advance shelf front (and calve).", & - default=.true.) + default=.false.) call get_param(param_file, mdl, "CALVE_TO_MASK", CS%calve_to_mask, & "If true, do not allow an ice shelf where prohibited by a mask.", & default=.false.) + endif call get_param(param_file, mdl, "MIN_THICKNESS_SIMPLE_CALVE", CS%min_thickness_simple_calve, & "Min thickness rule for the VERY simple calving law",& @@ -515,8 +520,13 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ call pass_var(CS%calve_mask,G%domain) endif + call initialize_ice_shelf_boundary_channel(CS%u_face_mask_bdry, CS%v_face_mask_bdry, & + CS%u_flux_bdry_val, CS%v_flux_bdry_val, CS%u_bdry_val, CS%v_bdry_val, CS%h_bdry_val, & + CS%thickness_bdry_val, ISS%hmask, ISS%h_shelf, G, & +! CS%flux_bdry, & + US, param_file ) !OVS initialize b.c.s ! call init_boundary_values(CS, G, time, ISS%hmask, CS%input_flux, CS%input_thickness, new_sim) - + call update_velocity_masks(CS, G, ISS%hmask, CS%umask, CS%vmask, CS%u_face_mask, CS%v_face_mask) if (new_sim) then call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: initialize ice velocity.") call update_OD_ffrac_uncoupled(CS, G, ISS%h_shelf(:,:)) @@ -823,8 +833,18 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, iters, time) ! need to make these conditional on GL interpolation float_cond(:,:) = 0.0 ; H_node(:,:) = 0.0 + CS%ground_frac(:,:) = 0.0 allocate(Phisub(nsub,nsub,2,2,2,2)) ; Phisub(:,:,:,:,:,:) = 0.0 + do j=G%jsc,G%jec + do i=G%isc,G%iec + if (rhoi_rhow * ISS%h_shelf(i,j) - G%bathyT(i,j) > 0) then + float_cond(i,j) = 1.0 + CS%ground_frac(i,j) = 1.0 + endif + enddo + enddo + call calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, CS%OD_av) ! This is to determine which cells contain the grounding line, the criterion being that the cell @@ -867,8 +887,9 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, iters, time) enddo ; enddo call calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) - call pass_var(CS%ice_visc, G%domain) + + call calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) call pass_var(CS%basal_traction, G%domain) ! This makes sure basal stress is only applied when it is supposed to be @@ -884,7 +905,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, iters, time) call CG_action(Au, Av, u_shlf, v_shlf, Phi, Phisub, CS%umask, CS%vmask, ISS%hmask, H_node, & CS%ice_visc, float_cond, G%bathyT, CS%basal_traction, & G, US, G%isc-1, G%iec+1, G%jsc-1, G%jec+1, rhoi_rhow) - + call pass_vector(Au,Av,G%domain) !OVS pass Au and Av if (CS%nonlin_solve_err_mode == 1) then err_init = 0 ; err_tempu = 0 ; err_tempv = 0 do J=G%IscB,G%JecB ; do I=G%IscB,G%IecB @@ -920,6 +941,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, iters, time) call calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) call pass_var(CS%ice_visc, G%domain) + call calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) call pass_var(CS%basal_traction, G%domain) ! makes sure basal stress is only applied when it is supposed to be @@ -986,8 +1008,11 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, iters, time) call MOM_mesg(mesg, 5) if (err_max <= CS%nonlinear_tolerance * err_init) then + write(mesg,*) "ice_shelf_solve_outer: nonlinear fractional residual = ", err_max/err_init + call MOM_mesg(mesg) write(mesg,*) "ice_shelf_solve_outer: exiting nonlinear solve after ",iter," iterations" - call MOM_mesg(mesg, 5) +! call MOM_mesg(mesg, 5) + call MOM_mesg(mesg) exit endif @@ -1073,7 +1098,7 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H rhoi_rhow = CS%density_ice / CS%density_ocean_avg Zu(:,:) = 0 ; Zv(:,:) = 0 ; DIAGu(:,:) = 0 ; DIAGv(:,:) = 0 - Ru(:,:) = 0 ; Rv(:,:) = 0 ; Au(:,:) = 0 ; Av(:,:) = 0 + Ru(:,:) = 0 ; Rv(:,:) = 0 ; Au(:,:) = 0 ; Av(:,:) = 0 ; RHSu(:,:) = 0 ; RHSv(:,:) = 0 Du(:,:) = 0 ; Dv(:,:) = 0 ; ubd(:,:) = 0 ; vbd(:,:) = 0 dot_p1 = 0 @@ -1125,8 +1150,8 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H do j=jsdq,jedq do i=isdq,iedq - if (CS%umask(I,J) == 1) Zu(I,J) = Ru(I,J) / DIAGu(I,J) - if (CS%vmask(I,J) == 1) Zv(I,J) = Rv(I,J) / DIAGv(I,J) + if (CS%umask(I,J) == 1 .AND.(DIAGu(I,J)/=0)) Zu(I,J) = Ru(I,J) / DIAGu(I,J) + if (CS%vmask(I,J) == 1 .AND.(DIAGv(I,J)/=0)) Zv(I,J) = Rv(I,J) / DIAGv(I,J) enddo enddo @@ -1161,7 +1186,7 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H ! Au, Av valid region moves in by 1 - + call pass_vector(Au,Av,G%domain, TO_ALL, BGRID_NE) sum_vec(:,:) = 0.0 ; sum_vec_2(:,:) = 0.0 do j=jscq,jecq ; do i=iscq,iecq @@ -1205,10 +1230,10 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H do j=jsdq,jedq do i=isdq,iedq - if (CS%umask(I,J) == 1) then + if (CS%umask(I,J) == 1 .AND.(DIAGu(I,J)/=0)) then Zu(I,J) = Ru(I,J) / DIAGu(I,J) endif - if (CS%vmask(I,J) == 1) then + if (CS%vmask(I,J) == 1 .AND.(DIAGv(I,J)/=0)) then Zv(I,J) = Rv(I,J) / DIAGv(I,J) endif enddo @@ -1732,7 +1757,7 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) BASE ! basal elevation of shelf/stream [Z ~> m]. - real :: rho, rhow ! Ice and ocean densities [R ~> kg m-3] + real :: rho, rhow, rhoi_rhow ! Ice and ocean densities [R ~> kg m-3] real :: sx, sy ! Ice shelf top slopes [Z L-1 ~> m s-1] real :: neumann_val ! [R Z L2 T-2 ~> kg s-2] real :: dxh, dyh ! Local grid spacing [L ~> m] @@ -1754,13 +1779,26 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) rho = CS%density_ice rhow = CS%density_ocean_avg grav = CS%g_Earth - + rhoi_rhow = rho/rhow ! prelim - go through and calculate S ! or is this faster? BASE(:,:) = -G%bathyT(:,:) + OD(:,:) S(:,:) = BASE(:,:) + ISS%h_shelf(:,:) + ! check whether the ice is floating or grounded + do j=jsc-1,jec+1 + do i=isc-1,iec+1 +! do i=isc-G%domain%nihalo,iec+G%domain%nihalo + +! if (ISS%h_shelf(i,j) < rhow/rho * G%bathyT(i,j)) then + if (rhoi_rhow * ISS%h_shelf(i,j) - G%bathyT(i,j) <= 0) then + S(i,j)=(1 - rhoi_rhow)*ISS%h_shelf(i,j) + endif + + + enddo + enddo do j=jsc-1,jec+1 do i=isc-1,iec+1 cnt = 0 @@ -1840,23 +1878,34 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) endif ! SW vertex - taudx(I-1,J-1) = taudx(I-1,J-1) - .25 * rho * grav * ISS%h_shelf(i,j) * sx * G%areaT(i,j) - taudy(I-1,J-1) = taudy(I-1,J-1) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j) - + if (ISS%hmask(I-1,J-1) == 1) then + if (CS%u_face_mask(I-1,J-1) /= 3) then + taudx(I-1,J-1) = taudx(I-1,J-1) - .25 * rho * grav * ISS%h_shelf(i,j) * sx * G%areaT(i,j) + taudy(I-1,J-1) = taudy(I-1,J-1) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j) + endif + endif ! SE vertex - taudx(I,J-1) = taudx(I,J-1) - .25 * rho * grav * ISS%h_shelf(i,j) * sx * G%areaT(i,j) - taudy(I,J-1) = taudy(I,J-1) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j) - + if (ISS%hmask(I,J-1) == 1) then + if (CS%u_face_mask(I,J-1) /= 3) then + taudx(I,J-1) = taudx(I,J-1) - .25 * rho * grav * ISS%h_shelf(i,j) * sx * G%areaT(i,j) + taudy(I,J-1) = taudy(I,J-1) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j) + endif + endif ! NW vertex - taudx(I-1,J) = taudx(I-1,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sx * G%areaT(i,j) - taudy(I-1,J) = taudy(I-1,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j) - + if (CS%u_face_mask(I-1,J) /= 3) then + taudx(I-1,J) = taudx(I-1,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sx * G%areaT(i,j) + taudy(I-1,J) = taudy(I-1,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j) + endif ! NE vertex - taudx(I,J) = taudx(I,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sx * G%areaT(i,j) - taudy(I,J) = taudy(I,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j) - - if (CS%ground_frac(i,j) == 1) then - neumann_val = .5 * grav * (rho * ISS%h_shelf(i,j)**2 - rhow * G%bathyT(i,j)**2) + if (ISS%hmask(I,J) == 1) then + if (CS%u_face_mask(I,J) /= 3) then + taudx(I,J) = taudx(I,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sx * G%areaT(i,j) + taudy(I,J) = taudy(I,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j) + endif + endif + if (CS%ground_frac(i,j) == 1) then +! neumann_val = .5 * grav * (rho * ISS%h_shelf(i,j)**2 - rhow * G%bathyT(i,j)**2) + neumann_val = .5 * grav * (1-rho/rhow) * rho * ISS%h_shelf(i,j)**2 else neumann_val = .5 * grav * (1-rho/rhow) * rho * ISS%h_shelf(i,j)**2 endif @@ -1976,7 +2025,7 @@ subroutine CG_action(uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask, hmas intent(inout) :: uret !< The retarding stresses working at u-points [R L3 Z T-2 ~> kg m s-2]. real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), & intent(inout) :: vret !< The retarding stresses working at v-points [R L3 Z T-2 ~> kg m s-2]. - real, dimension(SZDI_(G),SZDJ_(G),8,4), & + real, dimension(8,4,SZDI_(G),SZDJ_(G)), & intent(in) :: Phi !< The gradients of bilinear basis elements at Gaussian !! quadrature points surrounding the cell vertices [L-1 ~> m-1]. real, dimension(:,:,:,:,:,:), & @@ -2080,7 +2129,7 @@ subroutine CG_action(uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask, hmas v_shlf(I-1,J) * Phi(6,2*(jq-1)+iq,i,j) + & v_shlf(I,J) * Phi(8,2*(jq-1)+iq,i,j) - do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; Jtgt = J-2-jphi + do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; ;Jtgt = J-2+jphi !Jtgt = J-2-jphi !OVS fix index if (umask(Itgt,Jtgt) == 1) uret(Itgt,Jtgt) = uret(Itgt,Jtgt) + 0.25 * ice_visc(i,j) * & ((4*ux+2*vy) * Phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq,i,j) + & (uy+vx) * Phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq,i,j)) @@ -2214,7 +2263,7 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, ! Phi(2*i-1,j) gives d(Phi_i)/dx at quadrature point j ! Phi(2*i,j) gives d(Phi_i)/dy at quadrature point j - do iq=1,2 ; do jq=1,2 ; do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; Jtgt = J-2-jphi + do iq=1,2 ; do jq=1,2 ; do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; Jtgt = J-2+jphi !Jtgt = J-2-jphi !OVS fix index ilq = 1 ; if (iq == iphi) ilq = 2 jlq = 1 ; if (jq == jphi) jlq = 2 @@ -2258,7 +2307,7 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, if (float_cond(i,j) == 1) then Hcell(:,:) = H_node(i-1:i,j-1:j) call CG_diagonal_subgrid_basal(Phisub, Hcell, G%bathyT(i,j), dens_ratio, sub_ground) - do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; Jtgt = J-2-jphi + do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; Jtgt = J-2+jphi !Jtgt = J-2-jphi !OVS fix index if (CS%umask(Itgt,Jtgt) == 1) then u_diagonal(Itgt,Jtgt) = u_diagonal(Itgt,Jtgt) + sub_ground(iphi,jphi) * basal_trac(i,j) v_diagonal(Itgt,Jtgt) = v_diagonal(Itgt,Jtgt) + sub_ground(iphi,jphi) * basal_trac(i,j) @@ -2399,7 +2448,7 @@ subroutine apply_boundary_values(CS, ISS, G, US, time, Phisub, H_node, ice_visc, CS%v_bdry_val(I-1,J) * Phi(6,2*(jq-1)+iq) + & CS%v_bdry_val(I,J) * Phi(8,2*(jq-1)+iq) - do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; Jtgt = J-2-jphi + do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; Jtgt = J-2+jphi !Jtgt = J-2-jphi !OVS fix index ilq = 1 ; if (iq == iphi) ilq = 2 jlq = 1 ; if (jq == jphi) jlq = 2 @@ -2472,7 +2521,7 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) integer :: giec, gjec, gisc, gjsc, cnt, isc, jsc, iec, jec, is, js real :: Visc_coef, n_g real :: ux, uy, vx, vy, eps_min ! Velocity shears [T-1 ~> s-1] - real :: umid, vmid, unorm ! Velocities [L T-1 ~> m s-1] +! real :: umid, vmid, unorm ! Velocities [L T-1 ~> m s-1] isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec iscq = G%iscB ; iecq = G%iecB ; jscq = G%jscB ; jecq = G%jecB @@ -2484,7 +2533,7 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) n_g = CS%n_glen; eps_min = CS%eps_glen_min - Visc_coef = US%kg_m2s_to_RZ_T*US%m_to_L*US%Z_to_L*(CS%A_glen_isothermal)**(1./CS%n_glen) + Visc_coef = US%kg_m2s_to_RZ_T*US%m_to_L*US%Z_to_L*(CS%A_glen_isothermal)**(-1./CS%n_glen) !OVS '-' in the exponent do j=jsd+1,jed-1 do i=isd+1,ied-1 @@ -2497,6 +2546,50 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) CS%ice_visc(i,j) = 0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * & (US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2))**((1.-n_g)/(2.*n_g)) +! umid = ((u_shlf(I,J) + u_shlf(I-1,J-1)) + (u_shlf(I,J-1) + u_shlf(I-1,J))) * 0.25 +! vmid = ((v_shlf(I,J) + v_shlf(I-1,J-1)) + (v_shlf(I,J-1) + v_shlf(I-1,J))) * 0.25 +! unorm = sqrt(umid**2 + vmid**2 + eps_min**2*(G%dxT(i,j)**2 + G%dyT(i,j)**2)) +! CS%basal_traction(i,j) = G%areaT(i,j) * CS%C_basal_friction * (US%L_T_to_m_s*unorm)**(CS%n_basal_fric-1) + endif + enddo + enddo + +end subroutine calc_shelf_visc + +subroutine calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) + type(ice_shelf_dyn_CS), intent(inout) :: CS !< A pointer to the ice shelf control structure + type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe + !! the ice-shelf state + type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf. + type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors + real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), & + intent(inout) :: u_shlf !< The zonal ice shelf velocity [L T-1 ~> m s-1]. + real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), & + intent(inout) :: v_shlf !< The meridional ice shelf velocity [L T-1 ~> m s-1]. + +! also this subroutine updates the nonlinear part of the basal traction + +! this may be subject to change later... to make it "hybrid" + + integer :: i, j, iscq, iecq, jscq, jecq, isd, jsd, ied, jed, iegq, jegq + integer :: giec, gjec, gisc, gjsc, cnt, isc, jsc, iec, jec, is, js + real :: umid, vmid, unorm, eps_min ! Velocities [L T-1 ~> m s-1] + + isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec + iscq = G%iscB ; iecq = G%iecB ; jscq = G%jscB ; jecq = G%jecB + isd = G%isd ; jsd = G%jsd ; ied = G%ied ; jed = G%jed + iegq = G%iegB ; jegq = G%jegB + gisc = G%domain%nihalo+1 ; gjsc = G%domain%njhalo+1 + giec = G%domain%niglobal+gisc ; gjec = G%domain%njglobal+gjsc + is = iscq - 1; js = jscq - 1 + + eps_min = CS%eps_glen_min + + + do j=jsd+1,jed-1 + do i=isd+1,ied-1 + + if (ISS%hmask(i,j) == 1) then umid = ((u_shlf(I,J) + u_shlf(I-1,J-1)) + (u_shlf(I,J-1) + u_shlf(I-1,J))) * 0.25 vmid = ((v_shlf(I,J) + v_shlf(I-1,J-1)) + (v_shlf(I,J-1) + v_shlf(I-1,J))) * 0.25 unorm = sqrt(umid**2 + vmid**2 + eps_min**2*(G%dxT(i,j)**2 + G%dyT(i,j)**2)) @@ -2505,7 +2598,7 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) enddo enddo -end subroutine calc_shelf_visc +end subroutine calc_shelf_taub subroutine update_OD_ffrac(CS, G, US, ocean_mass, find_avg) type(ice_shelf_dyn_CS), intent(inout) :: CS !< A pointer to the ice shelf control structure @@ -2673,8 +2766,18 @@ subroutine bilinear_shape_fn_grid(G, i, j, Phi) xquad(2:4:2) = .5 * (1+sqrt(1./3)) ; yquad(3:4) = .5 * (1+sqrt(1./3)) do qpoint=1,4 - a = G%dxCv(i,J-1) * (1-yquad(qpoint)) + G%dxCv(i,J) * yquad(qpoint) ! d(x)/d(x*) - d = G%dyCu(I-1,j) * (1-xquad(qpoint)) + G%dyCu(I,j) * xquad(qpoint) ! d(y)/d(y*) + if (J>1) then + a = G%dxCv(i,J-1) * (1-yquad(qpoint)) + G%dxCv(i,J) * yquad(qpoint) ! d(x)/d(x*) + else + a= G%dxCv(i,J) !* yquad(qpoint) ! d(x)/d(x*) + endif + if (I>1) then + d = G%dyCu(I-1,j) * (1-xquad(qpoint)) + G%dyCu(I,j) * xquad(qpoint) ! d(y)/d(y*) + else + d = G%dyCu(I,j) !* xquad(qpoint) + endif +! a = G%dxCv(i,J-1) * (1-yquad(qpoint)) + G%dxCv(i,J) * yquad(qpoint) ! d(x)/d(x*) +! d = G%dyCu(I-1,j) * (1-xquad(qpoint)) + G%dyCu(I,j) * xquad(qpoint) ! d(y)/d(y*) do node=1,4 xnode = 2-mod(node,2) ; ynode = ceiling(REAL(node)/2) @@ -2793,21 +2896,28 @@ subroutine update_velocity_masks(CS, G, hmask, umask, vmask, u_face_mask, v_face is = isd+1 ; js = jsd+1 endif - do j=js,G%jed +! do j=js,G%jed + do j=js-1,G%jed !OVS change index do i=is,G%ied if (hmask(i,j) == 1) then - umask(I-1:I,j-1:j) = 1. - vmask(I-1:I,j-1:j) = 1. + umask(I,j) = 1. + vmask(I,j) = 1. do k=0,1 select case (int(CS%u_face_mask_bdry(I-1+k,j))) case (3) - umask(I-1+k,J-1:J)=3. - vmask(I-1+k,J-1:J)=0. + ! vmask(I-1+k,J-1)=0. u_face_mask(I-1+k,j)=3. + umask(I-1+k,J)=3. + !vmask(I-1+k,J)=0. + vmask(I-1+k,J)=3. + !u_face_mask(I-1+k,j-1)=3. +! umask(I-1+k,J-1:J)=3. +! vmask(I-1+k,J-1:J)=0. +! u_face_mask(I-1+k,j)=3. case (2) u_face_mask(I-1+k,j)=2. case (4) @@ -2815,9 +2925,9 @@ subroutine update_velocity_masks(CS, G, hmask, umask, vmask, u_face_mask, v_face vmask(I-1+k,J-1:J)=0. u_face_mask(I-1+k,j)=4. case (0) - umask(I-1+k,J-1:J)=0. - vmask(I-1+k,J-1:J)=0. - u_face_mask(I-1+k,j)=0. +! umask(I-1+k,J-1:J)=0. +! vmask(I-1+k,J-1:J)=0. +! u_face_mask(I-1+k,j)=0. case (1) ! stress free x-boundary umask(I-1+k,J-1:J)=0. case default @@ -2838,9 +2948,9 @@ subroutine update_velocity_masks(CS, G, hmask, umask, vmask, u_face_mask, v_face vmask(I-1:I,J-1+k)=0. v_face_mask(i,J-1+k)=4. case (0) - umask(I-1:I,J-1+k)=0. - vmask(I-1:I,J-1+k)=0. - v_face_mask(i,J-1+k)=0. +! umask(I-1:I,J-1+k)=0. +! vmask(I-1:I,J-1+k)=0. +! v_face_mask(i,J-1+k)=0. case (1) ! stress free y-boundary vmask(I-1:I,J-1+k)=0. case default diff --git a/src/ice_shelf/MOM_ice_shelf_initialize.F90 b/src/ice_shelf/MOM_ice_shelf_initialize.F90 index 90c98fa487..367f8d7dce 100644 --- a/src/ice_shelf/MOM_ice_shelf_initialize.F90 +++ b/src/ice_shelf/MOM_ice_shelf_initialize.F90 @@ -18,6 +18,7 @@ module MOM_ice_shelf_initialize !MJHpublic initialize_ice_shelf_boundary, initialize_ice_thickness public initialize_ice_thickness +public initialize_ice_shelf_boundary_channel ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with @@ -263,6 +264,130 @@ subroutine initialize_ice_thickness_channel(h_shelf, area_shelf_h, hmask, G, US, endif ; enddo end subroutine initialize_ice_thickness_channel +subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_bdry, & + u_flux_bdry_val, v_flux_bdry_val, u_bdry_val, v_bdry_val, h_bdry_val, & + thickness_bdry_val, hmask, h_shelf, G,& ! OVS h_shelf 11/10/20 +! flux_bdry, & + US, PF ) + + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + real, dimension(SZIB_(G),SZJ_(G)), & + intent(inout) :: u_face_mask_bdry !< A boundary-type mask at C-grid u faces + real, dimension(SZIB_(G),SZJ_(G)), & + intent(inout) :: u_flux_bdry_val !< The boundary thickness flux through + !! C-grid u faces [L Z T-1 ~> m2 s-1]. + real, dimension(SZI_(G),SZJB_(G)), & + intent(inout) :: v_face_mask_bdry !< A boundary-type mask at C-grid v faces + real, dimension(SZI_(G),SZJB_(G)), & + intent(inout) :: v_flux_bdry_val !< The boundary thickness flux through + !! C-grid v faces [L Z T-1 ~> m2 s-1]. + real, dimension(SZIB_(G),SZJB_(G)), & + intent(inout) :: u_bdry_val !< The zonal ice shelf velocity at open + !! boundary vertices [L T-1 ~> m s-1]. + real, dimension(SZIB_(G),SZJB_(G)), & + intent(inout) :: v_bdry_val !< The meridional ice shelf velocity at open + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(inout) :: thickness_bdry_val !< The ice shelf thickness at open boundaries [Z ~> m] + !! boundary vertices [L T-1 ~> m s-1]. + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(inout) :: h_bdry_val !< The ice shelf thickness at open boundaries [Z ~> m] + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(inout) :: hmask !< A mask indicating which tracer points are + !! partly or fully covered by an ice-shelf + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(inout) :: h_shelf !< Ice-shelf thickness OVS 11/10/20 +! logical, intent(in) :: flux_bdry !< If true, use mass fluxes as the boundary value. + type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors + type(param_file_type), intent(in) :: PF !< A structure to parse for run-time parameters + + character(len=40) :: mdl = "initialize_ice_shelf_boundary_channel" ! This subroutine's name. + integer :: i, j, isd, jsd, is, js, iegq, jegq, giec, gjec, gisc, gjsc,gisd,gjsd, isc, jsc, iec, jec, ied, jed + real :: input_thick ! The input ice shelf thickness [Z ~> m] +! real :: input_flux ! The input ice flux per unit length [L Z T-1 ~> m2 s-1] + real :: input_vel ! The input ice velocity per [L Z T-1 ~> m s-1] + real :: lenlat, len_stress, westlon, lenlon, southlat + + + call get_param(PF, mdl, "LENLAT", lenlat, fail_if_missing=.true.) + + call get_param(PF, mdl, "LENLON", lenlon, fail_if_missing=.true.) + + call get_param(PF, mdl, "WESTLON", westlon, fail_if_missing=.true.) + + call get_param(PF, mdl, "SOUTHLAT", southlat, fail_if_missing=.true.) + + call get_param(PF, mdl, "INPUT_VEL_ICE_SHELF", input_vel, & + "inflow ice velocity at upstream boundary", & + units="m s-1", default=0., scale=US%m_s_to_L_T*US%m_to_Z) + call get_param(PF, mdl, "INPUT_THICK_ICE_SHELF", input_thick, & + "flux thickness at upstream boundary", & + units="m", default=1000., scale=US%m_to_Z) + call get_param(PF, mdl, "LEN_SIDE_STRESS", len_stress, & + "maximum position of no-flow condition in along-flow direction", & + units="km", default=0.) + + call MOM_mesg(mdl//": setting boundary") + + isd = G%isd ; ied = G%ied + jsd = G%jsd ; jed = G%jed + isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec + gjsd = G%Domain%njglobal ; gisd = G%Domain%niglobal + gisc = G%Domain%nihalo ; gjsc = G%Domain%njhalo + giec = G%Domain%niglobal+gisc ; gjec = G%Domain%njglobal+gjsc + +!-----------b.c.s based on geopositions ----------------- + do j=jsc-1,jec+1 + do i=isc-1,iec+1 + ! upstream boundary - set either dirichlet or flux condition + + if (G%geoLonBu(i,j) == westlon) then + ! if (flux_bdry) then + ! u_face_mask_bdry(i-1,j) = 4.0 + ! u_flux_bdry_val(i-1,j) = input_flux + ! else + hmask(i+1,j) = 3.0 + h_bdry_val(i+1,j) = h_shelf(i+1,j) !OVS 11/10/20 !input_thick + thickness_bdry_val(i+1,j) = h_bdry_val(i+1,j) + u_face_mask_bdry(i+1,j) = 3.0 + u_bdry_val(i+1,j) = input_vel*(1-16.0*((G%geoLatBu(i-1,j)/lenlat-0.5))**4) !OVS 11/09/20 U b.c. + ! u_bdry_val(i+1,j) = (1 - ((G%geoLatBu(i,j) - 0.5*lenlat)*2./lenlat)**2) * & + ! 1.5 * input_flux / input_thick + ! endif + endif + + + ! side boundaries: no flow + if (G%geoLatBu(i,j-1) == southlat) then !bot boundary + if (len_stress == 0. .OR. G%geoLonCv(i,j) <= len_stress) then + v_face_mask_bdry(i,j+1) = 0. + u_face_mask_bdry(i,j-1) = 3. !OVS 11/25/20 + else + v_face_mask_bdry(i,j+1) = 1. + u_face_mask_bdry(i,j-1) = 3. !OVS 11/25/20 + u_bdry_val(i,j) = 0. + !hmask(i,j) = 0.0 !OVS 11/25/20 + endif + elseif (G%geoLatBu(i,j-1) == southlat+lenlat) then !top boundary + if (len_stress == 0. .OR. G%geoLonCv(i,j) <= len_stress) then + v_face_mask_bdry(i,j-1) = 0. + u_face_mask_bdry(i,j-1) = 3. !OVS 11/25/20 + else + v_face_mask_bdry(i,j-1) = 1. + u_face_mask_bdry(i,j-1) = 3. !OVS 11/25/20 + !u_bdry_val(i,j) = 0. !OVS 11/25/20 + !hmask(i,j) = 0.0 !OVS 11/25/20 + endif + endif + + ! downstream boundary - CFBC + if (G%geoLonBu(i,j) == westlon+lenlon) then + u_face_mask_bdry(i-1,j) = 2.0 + endif + + enddo + enddo + +end subroutine initialize_ice_shelf_boundary_channel !BEGIN MJH ! subroutine initialize_ice_shelf_boundary(u_face_mask_bdry, v_face_mask_bdry, & From 39dd3e3221bb19ef90d191f3c69f77cd50de03cb Mon Sep 17 00:00:00 2001 From: Olga Sergienko Date: Tue, 15 Dec 2020 21:05:04 -0500 Subject: [PATCH 02/17] Modified MOM_ice_shelf_dynamics.F90 --- src/ice_shelf/MOM_ice_shelf_dynamics.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index 8480906de8..174891582a 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -672,6 +672,7 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled call update_OD_ffrac_uncoupled(CS, G, ISS%h_shelf(:,:)) endif + if (update_ice_vel) then call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf, iters, Time) endif From ebac0adb8524b450882d679f5c74af71e50d50b7 Mon Sep 17 00:00:00 2001 From: Olga Sergienko Date: Tue, 29 Dec 2020 10:38:55 -0500 Subject: [PATCH 03/17] Modifications to register_diag_field in MOM_ice_shelf_dynamics to make ice-shelf_fields consistent with diag_table Modifications to MOM_ice_shelf.F90 to apply melting to the case of a dynamic ice shelf. --- src/ice_shelf/MOM_ice_shelf.F90 | 11 ++++++++++ src/ice_shelf/MOM_ice_shelf_dynamics.F90 | 28 ++++++++++++++++++------ 2 files changed, 32 insertions(+), 7 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 56461dbc3d..5663b326b7 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -710,6 +710,17 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step, CS) endif endif + ! Melting has been computed, now is time to update thickness and mass with dynamic ice shelf + if (CS%active_shelf_dynamics) then !OVS 12/10/20 + call change_thickness_using_melt(ISS, G, US, US%s_to_T*time_step, fluxes, CS%density_ice, CS%debug) !OVS 12/10/20 + + if (CS%debug) then + call hchksum(ISS%h_shelf, "h_shelf after change thickness using melt", G%HI, haloshift=0, scale=US%Z_to_m) + call hchksum(ISS%mass_shelf, "mass_shelf after change thickness using melt", G%HI, haloshift=0, & + scale=US%RZ_to_kg_m2) + endif + endif !OVS 12/10/20 + if (CS%debug) call MOM_forcing_chksum("Before add shelf flux", fluxes, G, CS%US, haloshift=0) call add_shelf_flux(G, US, CS, sfc_state, fluxes) diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index 174891582a..ead882dd75 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -537,21 +537,35 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ endif ! Register diagnostics. - CS%id_u_shelf = register_diag_field('ocean_model','u_shelf',CS%diag%axesCu1, Time, & +! CS%id_u_shelf = register_diag_field('ocean_model','u_shelf',CS%diag%axesCu1, Time, & +! 'x-velocity of ice', 'm yr-1', conversion=365.0*86400.0*US%L_T_to_m_s) + CS%id_u_shelf = register_diag_field('ice_shelf_model','u_shelf',CS%diag%axesCu1, Time, & 'x-velocity of ice', 'm yr-1', conversion=365.0*86400.0*US%L_T_to_m_s) - CS%id_v_shelf = register_diag_field('ocean_model','v_shelf',CS%diag%axesCv1, Time, & +! CS%id_v_shelf = register_diag_field('ocean_model','v_shelf',CS%diag%axesCv1, Time, & +! 'y-velocity of ice', 'm yr-1', conversion=365.0*86400.0*US%L_T_to_m_s) + CS%id_v_shelf = register_diag_field('ice_shelf_model','v_shelf',CS%diag%axesCv1, Time, & 'y-velocity of ice', 'm yr-1', conversion=365.0*86400.0*US%L_T_to_m_s) - CS%id_u_mask = register_diag_field('ocean_model','u_mask',CS%diag%axesCu1, Time, & +! CS%id_u_mask = register_diag_field('ocean_model','u_mask',CS%diag%axesCu1, Time, & +! 'mask for u-nodes', 'none') + CS%id_u_mask = register_diag_field('ice_shelf_model','u_mask',CS%diag%axesCu1, Time, & 'mask for u-nodes', 'none') - CS%id_v_mask = register_diag_field('ocean_model','v_mask',CS%diag%axesCv1, Time, & +! CS%id_v_mask = register_diag_field('ocean_model','v_mask',CS%diag%axesCv1, Time, & +! 'mask for v-nodes', 'none') + CS%id_v_mask = register_diag_field('ice_shelf_model','v_mask',CS%diag%axesCv1, Time, & 'mask for v-nodes', 'none') ! CS%id_surf_elev = register_diag_field('ocean_model','ice_surf',CS%diag%axesT1, Time, & ! 'ice surf elev', 'm') - CS%id_ground_frac = register_diag_field('ocean_model','ice_ground_frac',CS%diag%axesT1, Time, & +! CS%id_ground_frac = register_diag_field('ocean_model','ice_ground_frac',CS%diag%axesT1, Time, & +! 'fraction of cell that is grounded', 'none') + CS%id_ground_frac = register_diag_field('ice_shelf_model','ice_ground_frac',CS%diag%axesT1, Time, & 'fraction of cell that is grounded', 'none') - CS%id_col_thick = register_diag_field('ocean_model','col_thick',CS%diag%axesT1, Time, & +! CS%id_col_thick = register_diag_field('ocean_model','col_thick',CS%diag%axesT1, Time, & +! 'ocean column thickness passed to ice model', 'm', conversion=US%Z_to_m) + CS%id_col_thick = register_diag_field('ice_shelf_model','col_thick',CS%diag%axesT1, Time, & 'ocean column thickness passed to ice model', 'm', conversion=US%Z_to_m) - CS%id_OD_av = register_diag_field('ocean_model','OD_av',CS%diag%axesT1, Time, & +! CS%id_OD_av = register_diag_field('ocean_model','OD_av',CS%diag%axesT1, Time, & +! 'intermediate ocean column thickness passed to ice model', 'm', conversion=US%Z_to_m) + CS%id_OD_av = register_diag_field('ice_shelf_model','OD_av',CS%diag%axesT1, Time, & 'intermediate ocean column thickness passed to ice model', 'm', conversion=US%Z_to_m) !CS%id_h_after_uflux = register_diag_field('ocean_model','h_after_uflux',CS%diag%axesh1, Time, & ! 'thickness after u flux ', 'none') From f30f636b2d65853180b125bd1f935c6d956c816b Mon Sep 17 00:00:00 2001 From: Olga Sergienko Date: Mon, 8 Feb 2021 19:35:38 -0500 Subject: [PATCH 04/17] corrected indecises in computation of driving stresses --- src/ice_shelf/MOM_ice_shelf_dynamics.F90 | 114 ++++++++++++++------- src/ice_shelf/MOM_ice_shelf_initialize.F90 | 38 ++++++- 2 files changed, 110 insertions(+), 42 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index ead882dd75..ef884dc434 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -44,7 +44,10 @@ module MOM_ice_shelf_dynamics !! on q-points (B grid) [L T-1 ~> m s-1] real, pointer, dimension(:,:) :: v_shelf => NULL() !< the meridional velocity of the ice shelf/sheet !! on q-points (B grid) [L T-1 ~> m s-1] - + real, pointer, dimension(:,:) :: taudx_shelf => NULL() !< the driving stress of the ice shelf/sheet + !! on q-points (C grid) [Pa ~> Pa] + real, pointer, dimension(:,:) :: taudy_shelf => NULL() !< the meridional stress of the ice shelf/sheet + !! on q-points (C grid) [Pa ~> Pa] real, pointer, dimension(:,:) :: u_face_mask => NULL() !< mask for velocity boundary conditions on the C-grid !! u-face - this is because the FEM cares about FACES THAT GET INTEGRATED OVER, !! not vertices. Will represent boundary conditions on computational boundary @@ -152,6 +155,7 @@ module MOM_ice_shelf_dynamics !>@{ Diagnostic handles integer :: id_u_shelf = -1, id_v_shelf = -1, id_t_shelf = -1, & + id_taudx_shelf = -1, id_taudy_shelf = -1, & id_ground_frac = -1, id_col_thick = -1, id_OD_av = -1, & id_u_mask = -1, id_v_mask = -1, id_t_mask = -1 !>@} @@ -250,7 +254,8 @@ subroutine register_ice_shelf_dyn_restarts(G, param_file, CS, restart_CS) allocate( CS%basal_traction(isd:ied,jsd:jed) ) ; CS%basal_traction(:,:) = 0.0 allocate( CS%OD_av(isd:ied,jsd:jed) ) ; CS%OD_av(:,:) = 0.0 allocate( CS%ground_frac(isd:ied,jsd:jed) ) ; CS%ground_frac(:,:) = 0.0 - + allocate( CS%taudx_shelf(Isd:Ied,Jsd:Jed) ) ; CS%taudx_shelf(:,:) = 0.0 + allocate( CS%taudy_shelf(Isd:Ied,Jsd:Jed) ) ; CS%taudy_shelf(:,:) = 0.0 ! additional restarts for ice shelf state call register_restart_field(CS%u_shelf, "u_shelf", .false., restart_CS, & "ice sheet/shelf u-velocity", "m s-1", hor_grid='Bu') @@ -258,6 +263,10 @@ subroutine register_ice_shelf_dyn_restarts(G, param_file, CS, restart_CS) "ice sheet/shelf v-velocity", "m s-1", hor_grid='Bu') call register_restart_field(CS%t_shelf, "t_shelf", .true., restart_CS, & "ice sheet/shelf vertically averaged temperature", "deg C") + call register_restart_field(CS%taudx_shelf, "taudx_shelf", .true., restart_CS, & !OVS 02/8/21 + "ice sheet/shelf taudx-driving stress", "kPa") + call register_restart_field(CS%taudy_shelf, "taudy_shelf", .true., restart_CS, & !OVS 02/08/21 + "ice sheet/shelf taudy-driving stress", "kPa") call register_restart_field(CS%OD_av, "OD_av", .true., restart_CS, & "Average open ocean depth in a cell","m") call register_restart_field(CS%ground_frac, "ground_frac", .true., restart_CS, & @@ -521,7 +530,8 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ endif call initialize_ice_shelf_boundary_channel(CS%u_face_mask_bdry, CS%v_face_mask_bdry, & - CS%u_flux_bdry_val, CS%v_flux_bdry_val, CS%u_bdry_val, CS%v_bdry_val, CS%h_bdry_val, & + CS%u_flux_bdry_val, CS%v_flux_bdry_val, CS%u_bdry_val, CS%v_bdry_val, CS%u_shelf, CS%v_shelf,& + CS%h_bdry_val, & CS%thickness_bdry_val, ISS%hmask, ISS%h_shelf, G, & ! CS%flux_bdry, & US, param_file ) !OVS initialize b.c.s @@ -530,10 +540,12 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ if (new_sim) then call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: initialize ice velocity.") call update_OD_ffrac_uncoupled(CS, G, ISS%h_shelf(:,:)) - call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf, iters, Time) - +! call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf, iters, Time) + call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) !OVS 02/08/21 if (CS%id_u_shelf > 0) call post_data(CS%id_u_shelf, CS%u_shelf, CS%diag) if (CS%id_v_shelf > 0) call post_data(CS%id_v_shelf, CS%v_shelf,CS%diag) + if (CS%id_taudx_shelf > 0) call post_data(CS%id_taudx_shelf, CS%taudx_shelf,CS%diag) + if (CS%id_taudy_shelf > 0) call post_data(CS%id_taudy_shelf, CS%taudy_shelf,CS%diag) endif ! Register diagnostics. @@ -545,6 +557,10 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ ! 'y-velocity of ice', 'm yr-1', conversion=365.0*86400.0*US%L_T_to_m_s) CS%id_v_shelf = register_diag_field('ice_shelf_model','v_shelf',CS%diag%axesCv1, Time, & 'y-velocity of ice', 'm yr-1', conversion=365.0*86400.0*US%L_T_to_m_s) + CS%id_taudx_shelf = register_diag_field('ice_shelf_model','taudx_shelf',CS%diag%axesT1, Time, & + 'x-driving stress of ice', 'kPa', conversion=1.e-3*US%L_T_to_m_s) + CS%id_taudy_shelf = register_diag_field('ice_shelf_model','taudy_shelf',CS%diag%axesT1, Time, & + 'x-driving stress of ice', 'kPa', conversion=1.e-3*US%L_T_to_m_s) ! CS%id_u_mask = register_diag_field('ocean_model','u_mask',CS%diag%axesCu1, Time, & ! 'mask for u-nodes', 'none') CS%id_u_mask = register_diag_field('ice_shelf_model','u_mask',CS%diag%axesCu1, Time, & @@ -559,6 +575,7 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ ! 'fraction of cell that is grounded', 'none') CS%id_ground_frac = register_diag_field('ice_shelf_model','ice_ground_frac',CS%diag%axesT1, Time, & 'fraction of cell that is grounded', 'none') + ! CS%id_col_thick = register_diag_field('ocean_model','col_thick',CS%diag%axesT1, Time, & ! 'ocean column thickness passed to ice model', 'm', conversion=US%Z_to_m) CS%id_col_thick = register_diag_field('ice_shelf_model','col_thick',CS%diag%axesT1, Time, & @@ -575,10 +592,10 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ ! 'thickness after front adv ', 'none') !!! OVS vertically integrated temperature - CS%id_t_shelf = register_diag_field('ocean_model','t_shelf',CS%diag%axesT1, Time, & - 'T of ice', 'oC') - CS%id_t_mask = register_diag_field('ocean_model','tmask',CS%diag%axesT1, Time, & - 'mask for T-nodes', 'none') +! CS%id_t_shelf = register_diag_field('ocean_model','t_shelf',CS%diag%axesT1, Time, & +! 'T of ice', 'oC') +! CS%id_t_mask = register_diag_field('ocean_model','tmask',CS%diag%axesT1, Time, & +! 'mask for T-nodes', 'none') endif end subroutine initialize_ice_shelf_dyn @@ -615,8 +632,8 @@ subroutine initialize_diagnostic_fields(CS, ISS, G, US, Time) enddo enddo - call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf, iters, dummy_time) - +! call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf, iters, dummy_time) + call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) !OVS 02/08/21 end subroutine initialize_diagnostic_fields !> This function returns the global maximum advective timestep that can be taken based on the current @@ -676,7 +693,7 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled coupled_GL = .false. if (present(ocean_mass) .and. present(coupled_grounding)) coupled_GL = coupled_grounding - call ice_shelf_advect(CS, ISS, G, time_step, Time) +! call ice_shelf_advect(CS, ISS, G, time_step, Time) !OVS 02/08/21 CS%elapsed_velocity_time = CS%elapsed_velocity_time + time_step if (CS%elapsed_velocity_time >= CS%velocity_update_time_step) update_ice_vel = .true. @@ -688,7 +705,8 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled if (update_ice_vel) then - call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf, iters, Time) +! call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf, iters, Time) + call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) !OVS 02/08/21 endif call ice_shelf_temp(CS, ISS, G, US, time_step, ISS%water_flux, Time) @@ -699,6 +717,8 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled if (CS%id_u_shelf > 0) call post_data(CS%id_u_shelf, CS%u_shelf, CS%diag) if (CS%id_v_shelf > 0) call post_data(CS%id_v_shelf, CS%v_shelf, CS%diag) if (CS%id_t_shelf > 0) call post_data(CS%id_t_shelf,CS%t_shelf,CS%diag) + if (CS%id_taudx_shelf > 0) call post_data(CS%id_taudx_shelf, CS%taudx_shelf, CS%diag) + if (CS%id_taudy_shelf > 0) call post_data(CS%id_taudy_shelf, CS%taudy_shelf, CS%diag) if (CS%id_ground_frac > 0) call post_data(CS%id_ground_frac, CS%ground_frac,CS%diag) if (CS%id_OD_av >0) call post_data(CS%id_OD_av, CS%OD_av,CS%diag) @@ -801,7 +821,8 @@ subroutine ice_shelf_advect(CS, ISS, G, time_step, Time) end subroutine ice_shelf_advect -subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, iters, time) +!subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, iters, time) + subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, iters, Time) !OVS 02/08/21 type(ice_shelf_dyn_CS), intent(inout) :: CS !< The ice shelf dynamics control structure type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe !! the ice-shelf state @@ -861,7 +882,9 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, iters, time) enddo call calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, CS%OD_av) - + call pass_vector(taudx, taudy, G%domain, TO_ALL, BGRID_NE) !OVS 02/01/21 +! call pass_var(taudx, G%Domain) !OVS 01/21/21 +! call pass_var(taudy, G%Domain) !OVS 01/21/21 ! This is to determine which cells contain the grounding line, the criterion being that the cell ! is ice-covered, with some nodes floating and some grounded flotation condition is estimated by ! assuming topography is cellwise constant and H is bilinear in a cell; floating where @@ -1303,7 +1326,7 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H cg_halo = cg_halo - 1 if (cg_halo == 0) then - ! pass vectors + ! pass vectors call pass_vector(Du, Dv, G%domain, TO_ALL, BGRID_NE) call pass_vector(u_shlf, v_shlf, G%domain, TO_ALL, BGRID_NE) call pass_vector(Ru, Rv, G%domain, TO_ALL, BGRID_NE) @@ -1786,8 +1809,10 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) iscq = G%iscB ; iecq = G%iecB ; jscq = G%jscB ; jecq = G%jecB isd = G%isd ; jsd = G%jsd iegq = G%iegB ; jegq = G%jegB - gisc = G%domain%nihalo+1 ; gjsc = G%domain%njhalo+1 - giec = G%domain%niglobal+G%domain%nihalo ; gjec = G%domain%njglobal+G%domain%njhalo +! gisc = G%domain%nihalo+1 ; gjsc = G%domain%njhalo+1 + gisc = 0*G%domain%nihalo+1 ; gjsc = 0*G%domain%njhalo+1 +! giec = G%domain%niglobal+G%domain%nihalo ; gjec = G%domain%njglobal+G%domain%njhalo + giec = G%domain%niglobal+0*G%domain%nihalo ; gjec = G%domain%njglobal+0*G%domain%njhalo is = iscq - 1; js = jscq - 1 i_off = G%idg_offset ; j_off = G%jdg_offset @@ -1802,9 +1827,10 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) S(:,:) = BASE(:,:) + ISS%h_shelf(:,:) ! check whether the ice is floating or grounded - do j=jsc-1,jec+1 - do i=isc-1,iec+1 -! do i=isc-G%domain%nihalo,iec+G%domain%nihalo +! do j=jsc-1,jec+1 !OVS 02/02/21 +! do i=isc-1,iec+1 !OVS 02/02/21 + do j=jsc-G%domain%njhalo,jec+G%domain%njhalo !OVS 02/02/21 + do i=isc-G%domain%nihalo,iec+G%domain%nihalo !OVS 02/02/21 ! if (ISS%h_shelf(i,j) < rhow/rho * G%bathyT(i,j)) then if (rhoi_rhow * ISS%h_shelf(i,j) - G%bathyT(i,j) <= 0) then @@ -1816,6 +1842,8 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) enddo do j=jsc-1,jec+1 do i=isc-1,iec+1 +! do j=jsc-G%domain%njhalo+1,jec+G%domain%njhalo-1 !OVS 02/02/21 +! do i=isc-G%domain%nihalo+1,iec+G%domain%nihalo-1 !OVS 02/02/21 cnt = 0 sx = 0 sy = 0 @@ -1826,12 +1854,14 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) ! calculate sx if ((i+i_off) == gisc) then ! at left computational bdry - if (ISS%hmask(i+1,j) == 1) then +! if ((i-i_off) == gisc) then ! at left computational bdry !OVS 02/02/21 + if (ISS%hmask(i+1,j) == 1) then sx = (S(i+1,j)-S(i,j))/dxh else sx = 0 endif elseif ((i+i_off) == giec) then ! at east computational bdry +! elseif ((i-i_off) == giec) then ! at east computational bdry !OVS 02/02/21 if (ISS%hmask(i-1,j) == 1) then sx = (S(i,j)-S(i-1,j))/dxh else @@ -1861,12 +1891,14 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) ! calculate sy, similarly if ((j+j_off) == gjsc) then ! at south computational bdry +! if ((j-j_off) == gjsc) then ! at south computational bdry !OVS 02/02/21 if (ISS%hmask(i,j+1) == 1) then sy = (S(i,j+1)-S(i,j))/dyh else sy = 0 endif elseif ((j+j_off) == gjec) then ! at nprth computational bdry +! elseif ((j-j_off) == gjec) then ! at nprth computational bdry !OVS 02/02/21 if (ISS%hmask(i,j-1) == 1) then sy = (S(i,j)-S(i,j-1))/dyh else @@ -1894,29 +1926,31 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) ! SW vertex if (ISS%hmask(I-1,J-1) == 1) then - if (CS%u_face_mask(I-1,J-1) /= 3) then +! if (CS%u_face_mask(I-1,J-1) /= 3) then taudx(I-1,J-1) = taudx(I-1,J-1) - .25 * rho * grav * ISS%h_shelf(i,j) * sx * G%areaT(i,j) taudy(I-1,J-1) = taudy(I-1,J-1) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j) - endif +! endif endif ! SE vertex if (ISS%hmask(I,J-1) == 1) then - if (CS%u_face_mask(I,J-1) /= 3) then +! if (CS%u_face_mask(I,J-1) /= 3) then taudx(I,J-1) = taudx(I,J-1) - .25 * rho * grav * ISS%h_shelf(i,j) * sx * G%areaT(i,j) taudy(I,J-1) = taudy(I,J-1) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j) - endif +! endif endif ! NW vertex - if (CS%u_face_mask(I-1,J) /= 3) then + if (ISS%hmask(I-1,J) == 1) then +! if (CS%u_face_mask(I-1,J) /= 3) then taudx(I-1,J) = taudx(I-1,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sx * G%areaT(i,j) taudy(I-1,J) = taudy(I-1,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j) - endif +! endif + endif ! NE vertex if (ISS%hmask(I,J) == 1) then - if (CS%u_face_mask(I,J) /= 3) then +! if (CS%u_face_mask(I,J) /= 3) then taudx(I,J) = taudx(I,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sx * G%areaT(i,j) taudy(I,J) = taudy(I,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j) - endif +! endif endif if (CS%ground_frac(i,j) == 1) then ! neumann_val = .5 * grav * (rho * ISS%h_shelf(i,j)**2 - rhow * G%bathyT(i,j)**2) @@ -2550,8 +2584,8 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) Visc_coef = US%kg_m2s_to_RZ_T*US%m_to_L*US%Z_to_L*(CS%A_glen_isothermal)**(-1./CS%n_glen) !OVS '-' in the exponent - do j=jsd+1,jed-1 - do i=isd+1,ied-1 + do j=jsd+1,jed!-1 OVS 02/01/21 + do i=isd+1,ied!-1 OVS 02/01/21 if (ISS%hmask(i,j) == 1) then ux = ((u_shlf(I,J) + u_shlf(I,J-1)) - (u_shlf(I-1,J) + u_shlf(I-1,J-1))) / (2*G%dxT(i,j)) @@ -2601,8 +2635,8 @@ subroutine calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) eps_min = CS%eps_glen_min - do j=jsd+1,jed-1 - do i=isd+1,ied-1 + do j=jsd+1,jed!-1 OVS 02/01/21 + do i=isd+1,ied!-1 OVS 02/01/21 if (ISS%hmask(i,j) == 1) then umid = ((u_shlf(I,J) + u_shlf(I-1,J-1)) + (u_shlf(I,J-1) + u_shlf(I-1,J))) * 0.25 @@ -2911,8 +2945,8 @@ subroutine update_velocity_masks(CS, G, hmask, umask, vmask, u_face_mask, v_face is = isd+1 ; js = jsd+1 endif -! do j=js,G%jed - do j=js-1,G%jed !OVS change index + do j=js,G%jed +! do j=js-1,G%jed !OVS change index do i=is,G%ied if (hmask(i,j) == 1) then @@ -2953,8 +2987,12 @@ subroutine update_velocity_masks(CS, G, hmask, umask, vmask, u_face_mask, v_face select case (int(CS%v_face_mask_bdry(i,J-1+k))) case (3) - vmask(I-1:I,J-1+k)=3. - umask(I-1:I,J-1+k)=0. +! vmask(I-1:I,J-1+k)=3. +! umask(I-1:I,J-1+k)=0. + vmask(I-1,J-1+k)=3. + umask(I-1,J-1+k)=0. + vmask(I,J-1+k)=3. + umask(I,J-1+k)=0. v_face_mask(i,J-1+k)=3. case (2) v_face_mask(i,J-1+k)=2. diff --git a/src/ice_shelf/MOM_ice_shelf_initialize.F90 b/src/ice_shelf/MOM_ice_shelf_initialize.F90 index 367f8d7dce..7025e53981 100644 --- a/src/ice_shelf/MOM_ice_shelf_initialize.F90 +++ b/src/ice_shelf/MOM_ice_shelf_initialize.F90 @@ -266,7 +266,7 @@ subroutine initialize_ice_thickness_channel(h_shelf, area_shelf_h, hmask, G, US, end subroutine initialize_ice_thickness_channel subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_bdry, & u_flux_bdry_val, v_flux_bdry_val, u_bdry_val, v_bdry_val, h_bdry_val, & - thickness_bdry_val, hmask, h_shelf, G,& ! OVS h_shelf 11/10/20 + thickness_bdry_val, hmask, h_shelf, u_shelf, v_shelf, G,& ! OVS h_shelf 11/10/20 ! flux_bdry, & US, PF ) @@ -286,6 +286,10 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b !! boundary vertices [L T-1 ~> m s-1]. real, dimension(SZIB_(G),SZJB_(G)), & intent(inout) :: v_bdry_val !< The meridional ice shelf velocity at open + real, dimension(SZIB_(G),SZJB_(G)), & + intent(inout) :: u_shelf !< The zonal ice shelf velocity [L T-1 ~> m s-1]. + real, dimension(SZIB_(G),SZJB_(G)), & + intent(inout) :: v_shelf !< The meridional ice shelf velocity [L T-1 ~> m s-1]. real, dimension(SZDI_(G),SZDJ_(G)), & intent(inout) :: thickness_bdry_val !< The ice shelf thickness at open boundaries [Z ~> m] !! boundary vertices [L T-1 ~> m s-1]. @@ -362,9 +366,11 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b v_face_mask_bdry(i,j+1) = 0. u_face_mask_bdry(i,j-1) = 3. !OVS 11/25/20 else - v_face_mask_bdry(i,j+1) = 1. +! v_face_mask_bdry(i,j+1) = 1. + v_face_mask_bdry(i,j-1) = 3. !OVS 01/20/21 u_face_mask_bdry(i,j-1) = 3. !OVS 11/25/20 - u_bdry_val(i,j) = 0. +! u_bdry_val(i,j) = 0. +! v_bdry_val(i,j) = 0. !OVS 01/20/21 !hmask(i,j) = 0.0 !OVS 11/25/20 endif elseif (G%geoLatBu(i,j-1) == southlat+lenlat) then !top boundary @@ -372,7 +378,8 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b v_face_mask_bdry(i,j-1) = 0. u_face_mask_bdry(i,j-1) = 3. !OVS 11/25/20 else - v_face_mask_bdry(i,j-1) = 1. +! v_face_mask_bdry(i,j-1) = 1. + v_face_mask_bdry(i,j-1) = 3. !OVS 01/20/21 u_face_mask_bdry(i,j-1) = 3. !OVS 11/25/20 !u_bdry_val(i,j) = 0. !OVS 11/25/20 !hmask(i,j) = 0.0 !OVS 11/25/20 @@ -387,6 +394,29 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b enddo enddo + +! if (.not. G%symmetric) then +!! do j=G%jsd,G%jed +!! do i=G%isd,G%ied +! do j=jsc-1,jec+1 +! do i=isc-1,iec+1 +!! if (((i+G%idg_offset) == (G%domain%nihalo+1)).and.(u_face_mask_bdry(I-1,j) == 3)) then +! if (u_face_mask_bdry(I-1,j) == 3) then +! u_shelf(I-1,J-1) = u_bdry_val(I-1,J-1) +! u_shelf(I-1,J) = u_bdry_val(I-1,J) +! v_shelf(I-1,J-1) = v_bdry_val(I-1,J-1) +! v_shelf(I-1,J) = v_bdry_val(I-1,J) +! endif +!! if (((j+G%jdg_offset) == (G%domain%njhalo+1)).and.(v_face_mask_bdry(i,J-1) == 3)) then +! if (v_face_mask_bdry(I,j-1) == 3) then +! u_shelf(I-1,J-1) = u_bdry_val(I-1,J-1) +! u_shelf(I,J-1) = u_bdry_val(I,J-1) +! v_shelf(I-1,J-1) = v_bdry_val(I-1,J-1) +! v_shelf(I,J-1) = v_bdry_val(I,J-1) +! endif +! enddo +! enddo +! endif end subroutine initialize_ice_shelf_boundary_channel !BEGIN MJH From 775205208d4df98879e42ea66cd0e3b520f9a646 Mon Sep 17 00:00:00 2001 From: Olga Sergienko Date: Wed, 10 Feb 2021 12:39:08 -0500 Subject: [PATCH 05/17] fixed ice-shelf advection --- src/ice_shelf/MOM_ice_shelf_dynamics.F90 | 16 +++++++++------- src/ice_shelf/MOM_ice_shelf_initialize.F90 | 12 +++++++----- 2 files changed, 16 insertions(+), 12 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index ef884dc434..89c91172e1 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -558,9 +558,9 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ CS%id_v_shelf = register_diag_field('ice_shelf_model','v_shelf',CS%diag%axesCv1, Time, & 'y-velocity of ice', 'm yr-1', conversion=365.0*86400.0*US%L_T_to_m_s) CS%id_taudx_shelf = register_diag_field('ice_shelf_model','taudx_shelf',CS%diag%axesT1, Time, & - 'x-driving stress of ice', 'kPa', conversion=1.e-3*US%L_T_to_m_s) + 'x-driving stress of ice', 'kPa', conversion=1.e-9*US%L_T_to_m_s) CS%id_taudy_shelf = register_diag_field('ice_shelf_model','taudy_shelf',CS%diag%axesT1, Time, & - 'x-driving stress of ice', 'kPa', conversion=1.e-3*US%L_T_to_m_s) + 'x-driving stress of ice', 'kPa', conversion=1.e-9*US%L_T_to_m_s) ! CS%id_u_mask = register_diag_field('ocean_model','u_mask',CS%diag%axesCu1, Time, & ! 'mask for u-nodes', 'none') CS%id_u_mask = register_diag_field('ice_shelf_model','u_mask',CS%diag%axesCu1, Time, & @@ -693,7 +693,7 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled coupled_GL = .false. if (present(ocean_mass) .and. present(coupled_grounding)) coupled_GL = coupled_grounding -! call ice_shelf_advect(CS, ISS, G, time_step, Time) !OVS 02/08/21 + call ice_shelf_advect(CS, ISS, G, time_step, Time) !OVS 02/08/21 CS%elapsed_velocity_time = CS%elapsed_velocity_time + time_step if (CS%elapsed_velocity_time >= CS%velocity_update_time_step) update_ice_vel = .true. @@ -782,7 +782,7 @@ subroutine ice_shelf_advect(CS, ISS, G, time_step, Time) call ice_shelf_advect_thickness_x(CS, G, LB, time_step, ISS%hmask, ISS%h_shelf, h_after_uflux, uh_ice) ! call enable_averages(time_step, Time, CS%diag) -! call pass_var(h_after_uflux, G%domain) + call pass_var(h_after_uflux, G%domain) ! if (CS%id_h_after_uflux > 0) call post_data(CS%id_h_after_uflux, h_after_uflux, CS%diag) ! call disable_averaging(CS%diag) @@ -790,7 +790,7 @@ subroutine ice_shelf_advect(CS, ISS, G, time_step, Time) call ice_shelf_advect_thickness_y(CS, G, LB, time_step, ISS%hmask, h_after_uflux, h_after_vflux, vh_ice) ! call enable_averages(time_step, Time, CS%diag) -! call pass_var(h_after_vflux, G%domain) + call pass_var(h_after_vflux, G%domain) ! if (CS%id_h_after_vflux > 0) call post_data(CS%id_h_after_vflux, h_after_vflux, CS%diag) ! call disable_averaging(CS%diag) @@ -882,7 +882,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, ite enddo call calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, CS%OD_av) - call pass_vector(taudx, taudy, G%domain, TO_ALL, BGRID_NE) !OVS 02/01/21 +! call pass_vector(taudx, taudy, G%domain, TO_ALL, BGRID_NE) !OVS 02/01/21 ! call pass_var(taudx, G%Domain) !OVS 01/21/21 ! call pass_var(taudy, G%Domain) !OVS 01/21/21 ! This is to determine which cells contain the grounding line, the criterion being that the cell @@ -1842,6 +1842,8 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) enddo do j=jsc-1,jec+1 do i=isc-1,iec+1 +! do j=G%jsd+1,G%jed-1 !OVS 02/08/21 +! do i=G%isd+1,G%ied-1 !OVS 02/08/21 ! do j=jsc-G%domain%njhalo+1,jec+G%domain%njhalo-1 !OVS 02/02/21 ! do i=isc-G%domain%nihalo+1,iec+G%domain%nihalo-1 !OVS 02/02/21 cnt = 0 @@ -2594,7 +2596,7 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) vy = ((v_shlf(I,J) + v_shlf(I-1,J)) - (v_shlf(I,J-1) + v_shlf(I-1,J-1))) / (2*G%dyT(i,j)) CS%ice_visc(i,j) = 0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * & (US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2))**((1.-n_g)/(2.*n_g)) - +! CS%ice_visc(i,j) =1e15*(G%areaT(i,j) * ISS%h_shelf(i,j)) !OVS 02/09/21 constvisc ! umid = ((u_shlf(I,J) + u_shlf(I-1,J-1)) + (u_shlf(I,J-1) + u_shlf(I-1,J))) * 0.25 ! vmid = ((v_shlf(I,J) + v_shlf(I-1,J-1)) + (v_shlf(I,J-1) + v_shlf(I-1,J))) * 0.25 ! unorm = sqrt(umid**2 + vmid**2 + eps_min**2*(G%dxT(i,j)**2 + G%dyT(i,j)**2)) diff --git a/src/ice_shelf/MOM_ice_shelf_initialize.F90 b/src/ice_shelf/MOM_ice_shelf_initialize.F90 index 7025e53981..2bfe64677c 100644 --- a/src/ice_shelf/MOM_ice_shelf_initialize.F90 +++ b/src/ice_shelf/MOM_ice_shelf_initialize.F90 @@ -349,11 +349,13 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b ! u_face_mask_bdry(i-1,j) = 4.0 ! u_flux_bdry_val(i-1,j) = input_flux ! else - hmask(i+1,j) = 3.0 - h_bdry_val(i+1,j) = h_shelf(i+1,j) !OVS 11/10/20 !input_thick - thickness_bdry_val(i+1,j) = h_bdry_val(i+1,j) - u_face_mask_bdry(i+1,j) = 3.0 - u_bdry_val(i+1,j) = input_vel*(1-16.0*((G%geoLatBu(i-1,j)/lenlat-0.5))**4) !OVS 11/09/20 U b.c. +! hmask(i+1,j) = 3.0 + hmask(i,j) = 3.0 +! h_bdry_val(i+1,j) = h_shelf(i+1,j) !OVS 11/10/20 !input_thick + h_bdry_val(i,j) = h_shelf(i,j) + thickness_bdry_val(i+0*1,j) = h_bdry_val(i+0*1,j) + u_face_mask_bdry(i+0*1,j) = 3.0 + u_bdry_val(i+0*1,j) = input_vel*(1-16.0*((G%geoLatBu(i-1,j)/lenlat-0.5))**4) !OVS 11/09/20 U b.c. ! u_bdry_val(i+1,j) = (1 - ((G%geoLatBu(i,j) - 0.5*lenlat)*2./lenlat)**2) * & ! 1.5 * input_flux / input_thick ! endif From f0ae41c0caf262607b0c28b9e5e54bb313c7809a Mon Sep 17 00:00:00 2001 From: Olga Sergienko Date: Thu, 11 Feb 2021 17:33:02 -0500 Subject: [PATCH 06/17] modified viscosity computations --- src/ice_shelf/MOM_ice_shelf_dynamics.F90 | 204 +++++++++++++++++++---- 1 file changed, 168 insertions(+), 36 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index 89c91172e1..db3a49cfe9 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -925,10 +925,11 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, ite enddo ; enddo call calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) - call pass_var(CS%ice_visc, G%domain) - +! call pass_var(CS%ice_visc, G%domain) +! call pass_vector(CS%ice_visc, G%domain, TO_ALL, BGRID_NE) !OVS 02/11/21 call calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) call pass_var(CS%basal_traction, G%domain) +! call pass_vector(CS%ice_visc,CS%basal_traction, G%domain, TO_ALL, BGRID_NE) !OVS 02/11/21 ! This makes sure basal stress is only applied when it is supposed to be do j=G%jsd,G%jed ; do i=G%isd,G%ied @@ -1842,10 +1843,6 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) enddo do j=jsc-1,jec+1 do i=isc-1,iec+1 -! do j=G%jsd+1,G%jed-1 !OVS 02/08/21 -! do i=G%isd+1,G%ied-1 !OVS 02/08/21 -! do j=jsc-G%domain%njhalo+1,jec+G%domain%njhalo-1 !OVS 02/02/21 -! do i=isc-G%domain%nihalo+1,iec+G%domain%nihalo-1 !OVS 02/02/21 cnt = 0 sx = 0 sy = 0 @@ -1856,14 +1853,12 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) ! calculate sx if ((i+i_off) == gisc) then ! at left computational bdry -! if ((i-i_off) == gisc) then ! at left computational bdry !OVS 02/02/21 if (ISS%hmask(i+1,j) == 1) then sx = (S(i+1,j)-S(i,j))/dxh else sx = 0 endif elseif ((i+i_off) == giec) then ! at east computational bdry -! elseif ((i-i_off) == giec) then ! at east computational bdry !OVS 02/02/21 if (ISS%hmask(i-1,j) == 1) then sx = (S(i,j)-S(i-1,j))/dxh else @@ -1893,14 +1888,12 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) ! calculate sy, similarly if ((j+j_off) == gjsc) then ! at south computational bdry -! if ((j-j_off) == gjsc) then ! at south computational bdry !OVS 02/02/21 if (ISS%hmask(i,j+1) == 1) then sy = (S(i,j+1)-S(i,j))/dyh else sy = 0 endif elseif ((j+j_off) == gjec) then ! at nprth computational bdry -! elseif ((j-j_off) == gjec) then ! at nprth computational bdry !OVS 02/02/21 if (ISS%hmask(i,j-1) == 1) then sy = (S(i,j)-S(i,j-1))/dyh else @@ -1928,31 +1921,23 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) ! SW vertex if (ISS%hmask(I-1,J-1) == 1) then -! if (CS%u_face_mask(I-1,J-1) /= 3) then taudx(I-1,J-1) = taudx(I-1,J-1) - .25 * rho * grav * ISS%h_shelf(i,j) * sx * G%areaT(i,j) taudy(I-1,J-1) = taudy(I-1,J-1) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j) -! endif endif ! SE vertex if (ISS%hmask(I,J-1) == 1) then -! if (CS%u_face_mask(I,J-1) /= 3) then taudx(I,J-1) = taudx(I,J-1) - .25 * rho * grav * ISS%h_shelf(i,j) * sx * G%areaT(i,j) taudy(I,J-1) = taudy(I,J-1) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j) -! endif endif ! NW vertex if (ISS%hmask(I-1,J) == 1) then -! if (CS%u_face_mask(I-1,J) /= 3) then taudx(I-1,J) = taudx(I-1,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sx * G%areaT(i,j) taudy(I-1,J) = taudy(I-1,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j) -! endif endif ! NE vertex if (ISS%hmask(I,J) == 1) then -! if (CS%u_face_mask(I,J) /= 3) then taudx(I,J) = taudx(I,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sx * G%areaT(i,j) taudy(I,J) = taudy(I,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j) -! endif endif if (CS%ground_frac(i,j) == 1) then ! neumann_val = .5 * grav * (rho * ISS%h_shelf(i,j)**2 - rhow * G%bathyT(i,j)**2) @@ -2567,11 +2552,11 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) ! also this subroutine updates the nonlinear part of the basal traction ! this may be subject to change later... to make it "hybrid" - + real, dimension(SZDIB_(G),SZDJB_(G)) :: eII integer :: i, j, iscq, iecq, jscq, jecq, isd, jsd, ied, jed, iegq, jegq - integer :: giec, gjec, gisc, gjsc, cnt, isc, jsc, iec, jec, is, js + integer :: giec, gjec, gisc, gjsc, cnt, isc, jsc, iec, jec, is, js, i_off, j_off real :: Visc_coef, n_g - real :: ux, uy, vx, vy, eps_min ! Velocity shears [T-1 ~> s-1] + real :: ux, uy, vx, vy, eps_min, dxh, dyh ! Velocity shears [T-1 ~> s-1] ! real :: umid, vmid, unorm ! Velocities [L T-1 ~> m s-1] isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec @@ -2581,30 +2566,177 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) gisc = G%domain%nihalo+1 ; gjsc = G%domain%njhalo+1 giec = G%domain%niglobal+gisc ; gjec = G%domain%njglobal+gjsc is = iscq - 1; js = jscq - 1 + i_off = G%idg_offset ; j_off = G%jdg_offset n_g = CS%n_glen; eps_min = CS%eps_glen_min + CS%ice_visc(:,:) = 0.0 + eII(:,:) = (US%s_to_T**2 * (eps_min**2)) Visc_coef = US%kg_m2s_to_RZ_T*US%m_to_L*US%Z_to_L*(CS%A_glen_isothermal)**(-1./CS%n_glen) !OVS '-' in the exponent - - do j=jsd+1,jed!-1 OVS 02/01/21 - do i=isd+1,ied!-1 OVS 02/01/21 - - if (ISS%hmask(i,j) == 1) then - ux = ((u_shlf(I,J) + u_shlf(I,J-1)) - (u_shlf(I-1,J) + u_shlf(I-1,J-1))) / (2*G%dxT(i,j)) - vx = ((v_shlf(I,J) + v_shlf(I,J-1)) - (v_shlf(I-1,J) + v_shlf(I-1,J-1))) / (2*G%dxT(i,j)) - uy = ((u_shlf(I,J) + u_shlf(I-1,J)) - (u_shlf(I,J-1) + u_shlf(I-1,J-1))) / (2*G%dyT(i,j)) - vy = ((v_shlf(I,J) + v_shlf(I-1,J)) - (v_shlf(I,J-1) + v_shlf(I-1,J-1))) / (2*G%dyT(i,j)) - CS%ice_visc(i,j) = 0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * & - (US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2))**((1.-n_g)/(2.*n_g)) -! CS%ice_visc(i,j) =1e15*(G%areaT(i,j) * ISS%h_shelf(i,j)) !OVS 02/09/21 constvisc + call pass_vector(u_shlf, v_shlf, G%domain, TO_ALL, BGRID_NE) +! do j=jsc-1,jec+1 +! do i=isc-1,iec+1 +!! do j=jsd+1,jed!-1 OVS 02/01/21 +!! do i=isd+1,ied!-1 OVS 02/01/21 + +! if (ISS%hmask(i,j) == 1) then +! ux = ((u_shlf(I,J) + 0*u_shlf(I,J-1)) - (u_shlf(I-1,J) + 0*u_shlf(I-1,J-1))) / (G%dxT(i,j)) +! vx = ((v_shlf(I,J) + v_shlf(I,J-1)) - (v_shlf(I-1,J) + v_shlf(I-1,J-1))) / (2*G%dxT(i,j)) +! uy = ((u_shlf(I,J) + u_shlf(I-1,J)) - (u_shlf(I,J-1) + u_shlf(I-1,J-1))) / (2*G%dyT(i,j)) +! vy = ((v_shlf(I,J) + v_shlf(I-1,J)) - (v_shlf(I,J-1) + v_shlf(I-1,J-1))) / (2*G%dyT(i,j)) +!! ux = ((u_shlf(I,J) + u_shlf(I,J-1)) - (u_shlf(I-1,J) + u_shlf(I-1,J-1))) / (2*G%dxT(i,j)) +!! vx = ((v_shlf(I,J) + v_shlf(I,J-1)) - (v_shlf(I-1,J) + v_shlf(I-1,J-1))) / (2*G%dxT(i,j)) +!! uy = ((u_shlf(I,J) + u_shlf(I-1,J)) - (u_shlf(I,J-1) + u_shlf(I-1,J-1))) / (2*G%dyT(i,j)) +!! vy = ((v_shlf(I,J) + v_shlf(I-1,J)) - (v_shlf(I,J-1) + v_shlf(I-1,J-1))) / (2*G%dyT(i,j)) +! CS%ice_visc(i,j) = 0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * & +! (US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2))**((1.-n_g)/(2.*n_g)) +! CS%ice_visc(i,j) =1e15*(G%areaT(i,j) * ISS%h_shelf(i,j)) ! constant viscocity for debugging ! umid = ((u_shlf(I,J) + u_shlf(I-1,J-1)) + (u_shlf(I,J-1) + u_shlf(I-1,J))) * 0.25 ! vmid = ((v_shlf(I,J) + v_shlf(I-1,J-1)) + (v_shlf(I,J-1) + v_shlf(I-1,J))) * 0.25 ! unorm = sqrt(umid**2 + vmid**2 + eps_min**2*(G%dxT(i,j)**2 + G%dyT(i,j)**2)) ! CS%basal_traction(i,j) = G%areaT(i,j) * CS%C_basal_friction * (US%L_T_to_m_s*unorm)**(CS%n_basal_fric-1) - endif - enddo - enddo +! endif +! enddo +! enddo + + + do j=jsc-1,jec+1 + do i=isc-1,iec+1 + cnt = 0 + ux = 0 + uy = 0 + vx = 0 + vy = 0 + dxh = G%dxT(i,j) + dyh = G%dyT(i,j) + + if (ISS%hmask(i,j) == 1) then ! we are inside the global computational bdry, at an ice-filled cell + ! calculate sx +! if ((i+i_off) == gisc) then ! at left computational bdry +! if (ISS%hmask(i+1,j) == 1) then +! ux = (u_shlf(i+1,j)-u_shlf(i,j))/dxh +! vx = (v_shlf(i+1,j)-v_shlf(i,j))/dxh +! else +! ux = 0 +! vx = 0 +! endif +! elseif ((i+i_off) == giec) then ! at east computational bdry +! if (ISS%hmask(i-1,j) == 1) then +! ux = (u_shlf(i,j)-u_shlf(i-1,j))/dxh +! vx = (v_shlf(i,j)-v_shlf(i-1,j))/dxh +! else +! ux = 0 +! vx = 0 +! endif +! else ! interior + if (ISS%hmask(i+1,j) == 1) then + cnt = cnt+1 + ux = u_shlf(i+1,j) + vx = v_shlf(i+1,j) + else + ux = u_shlf(i,j) + vx = v_shlf(i,j) + endif + if (ISS%hmask(i-1,j) == 1) then + cnt = cnt+1 + ux = ux - u_shlf(i-1,j) + vx = vx - v_shlf(i-1,j) + else + ux = ux - u_shlf(i,j) + vx = vx - v_shlf(i,j) + endif + if (cnt == 0) then + ux = 0 + vx = 0 + else + ux = ux / (cnt * dxh) + vx = vx / (cnt * dxh) + endif +! endif + cnt = 0 + + ! calculate sy, similarly +! if ((j+j_off) == gjsc) then ! at south computational bdry +! if (ISS%hmask(i,j+1) == 1) then +! uy = (u_shlf(i,j+1)-u_shlf(i,j))/dyh +! vy = (v_shlf(i,j+1)-v_shlf(i,j))/dyh +! else +! vy = 0 +! endif +! elseif ((j+j_off) == gjec) then ! at nprth computational bdry +! if (ISS%hmask(i,j-1) == 1) then +! uy = (u_shlf(i,j)-u_shlf(i,j-1))/dyh +! vy = (v_shlf(i,j)-v_shlf(i,j-1))/dyh +! else +! uy = 0 +! vy = 0 +! endif +! else ! interior + if (ISS%hmask(i,j+1) == 1) then + cnt = cnt+1 + uy = u_shlf(i,j+1) + vy = v_shlf(i,j+1) + else + uy = u_shlf(i,j) + vy = v_shlf(i,j) + endif + if (ISS%hmask(i,j-1) == 1) then + cnt = cnt+1 + uy = uy - u_shlf(i,j-1) + vy = vy - v_shlf(i,j-1) + else + uy = uy - u_shlf(i,j) + vy = vy - v_shlf(i,j) + endif + if (cnt == 0) then + uy = 0 + vy = 0 + else + uy = uy / (cnt * dyh) + vy = vy / (cnt * dyh) + endif +! endif + + ! SW vertex + if (ISS%hmask(I-1,J-1) == 1) then + eII(i-1,j-1) = eII(i-1,j-1)+.125*(US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2)) + endif + ! SE vertex + if (ISS%hmask(I,J-1) == 1) then + eII(i,j-1) = eII(i,j-1)+.125*(US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2)) + +! CS%ice_visc(i,j-1) = CS%ice_visc(i,j-1)+.25*0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * & +! (US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2))**((1.-n_g)/(2.*n_g)) + endif + ! NW vertex + if (ISS%hmask(I-1,J) == 1) then + eII(i-1,j) = eII(i-1,j)+.125*(US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2)) + +! CS%ice_visc(i-1,j) = CS%ice_visc(i-1,j)+.25*0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * & +! (US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2))**((1.-n_g)/(2.*n_g)) + endif + ! NE vertex + if (ISS%hmask(I,J) == 1) then + eII(i,j) = eII(i,j)+.25*(US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2)) + +! CS%ice_visc(i,j) = CS%ice_visc(i,j)+.25*0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * & +! (US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2))**((1.-n_g)/(2.*n_g)) + endif + if (ISS%hmask(I+1,J+1) == 1) then + eII(i+1,j+1) = eII(i+1,j+1)+.125*(US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2)) + endif + if (ISS%hmask(I,J+1) == 1) then + eII(i,j+1) = eII(i,j+1)+.125*(US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2)) + endif + if (ISS%hmask(I+1,J) == 1) then + eII(i+1,j) = eII(i+1,j)+.125*(US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2)) + endif + endif + CS%ice_visc(i,j) =0.5 * Visc_coef*(G%areaT(i,j) * ISS%h_shelf(i,j))*eII(i,j)**((1.-n_g)/(2.*n_g)) +! CS%ice_visc(i,j) =1e15*(G%areaT(i,j) * ISS%h_shelf(i,j)) !constant viscosity for debugging + enddo + enddo end subroutine calc_shelf_visc subroutine calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) From 89f4386eb153c63e55b942c4ec0fc8187d1ae8d5 Mon Sep 17 00:00:00 2001 From: Olga Sergienko Date: Thu, 18 Feb 2021 17:16:00 -0500 Subject: [PATCH 07/17] corrected initialize_boundary_channel call --- src/ice_shelf/MOM_ice_shelf_initialize.F90 | 28 ++++++++++++---------- 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf_initialize.F90 b/src/ice_shelf/MOM_ice_shelf_initialize.F90 index 2bfe64677c..532729c58c 100644 --- a/src/ice_shelf/MOM_ice_shelf_initialize.F90 +++ b/src/ice_shelf/MOM_ice_shelf_initialize.F90 @@ -265,9 +265,8 @@ subroutine initialize_ice_thickness_channel(h_shelf, area_shelf_h, hmask, G, US, end subroutine initialize_ice_thickness_channel subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_bdry, & - u_flux_bdry_val, v_flux_bdry_val, u_bdry_val, v_bdry_val, h_bdry_val, & - thickness_bdry_val, hmask, h_shelf, u_shelf, v_shelf, G,& ! OVS h_shelf 11/10/20 -! flux_bdry, & + u_flux_bdry_val, v_flux_bdry_val, u_bdry_val, v_bdry_val, u_shelf, v_shelf, h_bdry_val, & + thickness_bdry_val, hmask, h_shelf, G,& US, PF ) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure @@ -340,7 +339,8 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b giec = G%Domain%niglobal+gisc ; gjec = G%Domain%njglobal+gjsc !-----------b.c.s based on geopositions ----------------- - do j=jsc-1,jec+1 +! do j=jsc-1,jec+1 + do j=jsc-0*1,jec+1 do i=isc-1,iec+1 ! upstream boundary - set either dirichlet or flux condition @@ -349,13 +349,13 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b ! u_face_mask_bdry(i-1,j) = 4.0 ! u_flux_bdry_val(i-1,j) = input_flux ! else -! hmask(i+1,j) = 3.0 - hmask(i,j) = 3.0 -! h_bdry_val(i+1,j) = h_shelf(i+1,j) !OVS 11/10/20 !input_thick - h_bdry_val(i,j) = h_shelf(i,j) - thickness_bdry_val(i+0*1,j) = h_bdry_val(i+0*1,j) - u_face_mask_bdry(i+0*1,j) = 3.0 - u_bdry_val(i+0*1,j) = input_vel*(1-16.0*((G%geoLatBu(i-1,j)/lenlat-0.5))**4) !OVS 11/09/20 U b.c. + hmask(i+1,j) = 3.0 +! hmask(i,j) = 3.0 + h_bdry_val(i+1,j) = h_shelf(i+1,j) !OVS 11/10/20 !input_thick +! h_bdry_val(i,j) = h_shelf(i,j) + thickness_bdry_val(i+1,j) = h_bdry_val(i+0*1,j) + u_face_mask_bdry(i+1,j) = 3.0 + u_bdry_val(i+1,j) = input_vel*(1-16.0*((G%geoLatBu(i-1,j)/lenlat-0.5))**4) !OVS 11/09/20 U b.c. ! u_bdry_val(i+1,j) = (1 - ((G%geoLatBu(i,j) - 0.5*lenlat)*2./lenlat)**2) * & ! 1.5 * input_flux / input_thick ! endif @@ -396,7 +396,11 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b enddo enddo - +! call pass_var(hmask, G%domain) +! call pass_var(h_bdry_val, G%domain) +! call pass_var(thickness_bdry_val, G%domain) + + ! if (.not. G%symmetric) then !! do j=G%jsd,G%jed !! do i=G%isd,G%ied From 271bfce9402728937542e296e986c6aa4c172337 Mon Sep 17 00:00:00 2001 From: Olga Sergienko Date: Mon, 22 Feb 2021 15:10:32 -0500 Subject: [PATCH 08/17] corrected boundary mask in init_boundary_channel and updated u_ and v_bdry_val through halo --- src/ice_shelf/MOM_ice_shelf_dynamics.F90 | 342 +++++++++++++-------- src/ice_shelf/MOM_ice_shelf_initialize.F90 | 15 +- 2 files changed, 219 insertions(+), 138 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index db3a49cfe9..84605e3092 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -160,7 +160,7 @@ module MOM_ice_shelf_dynamics id_u_mask = -1, id_v_mask = -1, id_t_mask = -1 !>@} ! ids for outputting intermediate thickness in advection subroutine (debugging) - !integer :: id_h_after_uflux = -1, id_h_after_vflux = -1, id_h_after_adv = -1 + integer :: id_h_after_uflux = -1, id_h_after_vflux = -1, id_h_after_adv = -1, id_visc_shelf = -1 type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to control diagnostic output. @@ -535,18 +535,27 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ CS%thickness_bdry_val, ISS%hmask, ISS%h_shelf, G, & ! CS%flux_bdry, & US, param_file ) !OVS initialize b.c.s + + call pass_var(ISS%hmask, G%domain) + call pass_var(CS%h_bdry_val, G%domain) + call pass_var(CS%thickness_bdry_val, G%domain) + call pass_var(CS%u_bdry_val, G%domain) + call pass_var(CS%v_bdry_val, G%domain) + call pass_var(CS%u_face_mask_bdry, G%domain) + call pass_var(CS%v_face_mask_bdry, G%domain) ! call init_boundary_values(CS, G, time, ISS%hmask, CS%input_flux, CS%input_thickness, new_sim) call update_velocity_masks(CS, G, ISS%hmask, CS%umask, CS%vmask, CS%u_face_mask, CS%v_face_mask) - if (new_sim) then - call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: initialize ice velocity.") - call update_OD_ffrac_uncoupled(CS, G, ISS%h_shelf(:,:)) -! call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf, iters, Time) - call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) !OVS 02/08/21 - if (CS%id_u_shelf > 0) call post_data(CS%id_u_shelf, CS%u_shelf, CS%diag) - if (CS%id_v_shelf > 0) call post_data(CS%id_v_shelf, CS%v_shelf,CS%diag) - if (CS%id_taudx_shelf > 0) call post_data(CS%id_taudx_shelf, CS%taudx_shelf,CS%diag) - if (CS%id_taudy_shelf > 0) call post_data(CS%id_taudy_shelf, CS%taudy_shelf,CS%diag) - endif +! if (new_sim) then +! call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: initialize ice velocity.") +! call update_OD_ffrac_uncoupled(CS, G, ISS%h_shelf(:,:)) +!! call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf, iters, Time) +! call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) !OVS 02/08/21 +! if (CS%id_u_shelf > 0) call post_data(CS%id_u_shelf, CS%u_shelf, CS%diag) +! if (CS%id_v_shelf > 0) call post_data(CS%id_v_shelf, CS%v_shelf,CS%diag) +! if (CS%id_taudx_shelf > 0) call post_data(CS%id_taudx_shelf, CS%taudx_shelf,CS%diag) +! if (CS%id_taudy_shelf > 0) call post_data(CS%id_taudy_shelf, CS%taudy_shelf,CS%diag) +! if (CS%id_visc_shelf > 0) call post_data(CS%id_visc_shelf, CS%ice_visc,CS%diag) +! endif ! Register diagnostics. ! CS%id_u_shelf = register_diag_field('ocean_model','u_shelf',CS%diag%axesCu1, Time, & @@ -580,17 +589,29 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ ! 'ocean column thickness passed to ice model', 'm', conversion=US%Z_to_m) CS%id_col_thick = register_diag_field('ice_shelf_model','col_thick',CS%diag%axesT1, Time, & 'ocean column thickness passed to ice model', 'm', conversion=US%Z_to_m) + CS%id_visc_shelf = register_diag_field('ice_shelf_model','ice_visc',CS%diag%axesT1, Time, & + 'viscosity', 'm', conversion=1e-6*US%Z_to_m) ! CS%id_OD_av = register_diag_field('ocean_model','OD_av',CS%diag%axesT1, Time, & ! 'intermediate ocean column thickness passed to ice model', 'm', conversion=US%Z_to_m) CS%id_OD_av = register_diag_field('ice_shelf_model','OD_av',CS%diag%axesT1, Time, & 'intermediate ocean column thickness passed to ice model', 'm', conversion=US%Z_to_m) - !CS%id_h_after_uflux = register_diag_field('ocean_model','h_after_uflux',CS%diag%axesh1, Time, & - ! 'thickness after u flux ', 'none') - !CS%id_h_after_vflux = register_diag_field('ocean_model','h_after_vflux',CS%diag%axesh1, Time, & - ! 'thickness after v flux ', 'none') - !CS%id_h_after_adv = register_diag_field('ocean_model','h_after_adv',CS%diag%axesh1, Time, & - ! 'thickness after front adv ', 'none') - + CS%id_h_after_uflux = register_diag_field('ice_shelf_model','h_after_uflux',CS%diag%axesT1, Time, & + 'thickness after u flux ', 'none') + CS%id_h_after_vflux = register_diag_field('ice_shelf_model','h_after_vflux',CS%diag%axesT1, Time, & + 'thickness after v flux ', 'none') + CS%id_h_after_adv = register_diag_field('ice_shelf_model','h_after_adv',CS%diag%axesT1, Time, & + 'thickness after front adv ', 'none') + if (new_sim) then + call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: initialize ice velocity.") + call update_OD_ffrac_uncoupled(CS, G, ISS%h_shelf(:,:)) +! call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf, iters, Time) + call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) !OVS 02/08/21 + if (CS%id_u_shelf > 0) call post_data(CS%id_u_shelf, CS%u_shelf, CS%diag) + if (CS%id_v_shelf > 0) call post_data(CS%id_v_shelf, CS%v_shelf,CS%diag) + if (CS%id_taudx_shelf > 0) call post_data(CS%id_taudx_shelf, CS%taudx_shelf,CS%diag) + if (CS%id_taudy_shelf > 0) call post_data(CS%id_taudy_shelf, CS%taudy_shelf,CS%diag) + if (CS%id_visc_shelf > 0) call post_data(CS%id_visc_shelf, CS%ice_visc,CS%diag) + endif !!! OVS vertically integrated temperature ! CS%id_t_shelf = register_diag_field('ocean_model','t_shelf',CS%diag%axesT1, Time, & ! 'T of ice', 'oC') @@ -693,7 +714,7 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled coupled_GL = .false. if (present(ocean_mass) .and. present(coupled_grounding)) coupled_GL = coupled_grounding - call ice_shelf_advect(CS, ISS, G, time_step, Time) !OVS 02/08/21 +! call ice_shelf_advect(CS, ISS, G, time_step, Time) !OVS 02/08/21 CS%elapsed_velocity_time = CS%elapsed_velocity_time + time_step if (CS%elapsed_velocity_time >= CS%velocity_update_time_step) update_ice_vel = .true. @@ -721,6 +742,7 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled if (CS%id_taudy_shelf > 0) call post_data(CS%id_taudy_shelf, CS%taudy_shelf, CS%diag) if (CS%id_ground_frac > 0) call post_data(CS%id_ground_frac, CS%ground_frac,CS%diag) if (CS%id_OD_av >0) call post_data(CS%id_OD_av, CS%OD_av,CS%diag) + if (CS%id_visc_shelf > 0) call post_data(CS%id_visc_shelf, CS%ice_visc,CS%diag) if (CS%id_u_mask > 0) call post_data(CS%id_u_mask,CS%umask,CS%diag) if (CS%id_v_mask > 0) call post_data(CS%id_v_mask,CS%vmask,CS%diag) @@ -783,7 +805,7 @@ subroutine ice_shelf_advect(CS, ISS, G, time_step, Time) ! call enable_averages(time_step, Time, CS%diag) call pass_var(h_after_uflux, G%domain) -! if (CS%id_h_after_uflux > 0) call post_data(CS%id_h_after_uflux, h_after_uflux, CS%diag) + if (CS%id_h_after_uflux > 0) call post_data(CS%id_h_after_uflux, h_after_uflux, CS%diag) ! call disable_averaging(CS%diag) LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc ; LB%jeh = G%jec @@ -791,7 +813,7 @@ subroutine ice_shelf_advect(CS, ISS, G, time_step, Time) ! call enable_averages(time_step, Time, CS%diag) call pass_var(h_after_vflux, G%domain) -! if (CS%id_h_after_vflux > 0) call post_data(CS%id_h_after_vflux, h_after_vflux, CS%diag) + if (CS%id_h_after_vflux > 0) call post_data(CS%id_h_after_vflux, h_after_vflux, CS%diag) ! call disable_averaging(CS%diag) do j=jsd,jed @@ -882,7 +904,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, ite enddo call calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, CS%OD_av) -! call pass_vector(taudx, taudy, G%domain, TO_ALL, BGRID_NE) !OVS 02/01/21 + call pass_vector(taudx, taudy, G%domain, TO_ALL, BGRID_NE) !OVS 02/01/21 ! call pass_var(taudx, G%Domain) !OVS 01/21/21 ! call pass_var(taudy, G%Domain) !OVS 01/21/21 ! This is to determine which cells contain the grounding line, the criterion being that the cell @@ -925,7 +947,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, ite enddo ; enddo call calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) -! call pass_var(CS%ice_visc, G%domain) + call pass_var(CS%ice_visc, G%domain) ! call pass_vector(CS%ice_visc, G%domain, TO_ALL, BGRID_NE) !OVS 02/11/21 call calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) call pass_var(CS%basal_traction, G%domain) @@ -1329,7 +1351,9 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H if (cg_halo == 0) then ! pass vectors call pass_vector(Du, Dv, G%domain, TO_ALL, BGRID_NE) - call pass_vector(u_shlf, v_shlf, G%domain, TO_ALL, BGRID_NE) +! call pass_vector(u_shlf, v_shlf, G%domain, TO_ALL, BGRID_NE) + call pass_var(u_shlf, G%domain) + call pass_var(v_shlf, G%domain) call pass_vector(Ru, Rv, G%domain, TO_ALL, BGRID_NE) cg_halo = 3 endif @@ -2531,6 +2555,8 @@ subroutine apply_boundary_values(CS, ISS, G, US, time, Phisub, H_node, ice_visc, endif endif ; enddo ; enddo + call pass_vector(u_bdry_contr, v_bdry_contr, G%domain, TO_ALL, BGRID_NE) !OVS 02/19/21 + end subroutine apply_boundary_values !> Update depth integrated viscosity, based on horizontal strain rates, and also update the @@ -2552,11 +2578,14 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) ! also this subroutine updates the nonlinear part of the basal traction ! this may be subject to change later... to make it "hybrid" - real, dimension(SZDIB_(G),SZDJB_(G)) :: eII - integer :: i, j, iscq, iecq, jscq, jecq, isd, jsd, ied, jed, iegq, jegq +! real, dimension(SZDIB_(G),SZDJB_(G)) :: eII, ux, uy, vx, vy + integer :: i, j, iscq, iecq, jscq, jecq, isd, jsd, ied, jed, iegq, jegq, iq, jq integer :: giec, gjec, gisc, gjsc, cnt, isc, jsc, iec, jec, is, js, i_off, j_off real :: Visc_coef, n_g - real :: ux, uy, vx, vy, eps_min, dxh, dyh ! Velocity shears [T-1 ~> s-1] + real :: ux, uy, vx, vy + real :: eps_min, dxh, dyh ! Velocity shears [T-1 ~> s-1] + real, dimension(8,4) :: Phi + real, dimension(2) :: xquad ! real :: umid, vmid, unorm ! Velocities [L T-1 ~> m s-1] isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec @@ -2570,48 +2599,95 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) n_g = CS%n_glen; eps_min = CS%eps_glen_min - CS%ice_visc(:,:) = 0.0 - eII(:,:) = (US%s_to_T**2 * (eps_min**2)) +! CS%ice_visc(:,:) = 0.0 +! ux(:,:) = 0.0; uy(:,:) = 0.0; vx(:,:) =0.0; vy(:,:) =0.0 +! eII(:,:) = (US%s_to_T**2 * (eps_min**2)) Visc_coef = US%kg_m2s_to_RZ_T*US%m_to_L*US%Z_to_L*(CS%A_glen_isothermal)**(-1./CS%n_glen) !OVS '-' in the exponent - call pass_vector(u_shlf, v_shlf, G%domain, TO_ALL, BGRID_NE) +! call pass_vector(u_shlf, v_shlf, G%domain, TO_ALL, BGRID_NE) ! do j=jsc-1,jec+1 ! do i=isc-1,iec+1 -!! do j=jsd+1,jed!-1 OVS 02/01/21 -!! do i=isd+1,ied!-1 OVS 02/01/21 - -! if (ISS%hmask(i,j) == 1) then + do j=jsd+1,jed-1 !OVS 02/01/21 + do i=isd+1,ied-1 !OVS 02/01/21 + + if ((ISS%hmask(i,j) == 1) .OR. (ISS%hmask(i,j) == 3)) then +! ux(i,j) = ((u_shlf(I,J) + u_shlf(I,J-1)) - (u_shlf(I-1,J) + u_shlf(I-1,J-1))) / (2*G%dxT(i,j)) +! vx(i,j) = ((v_shlf(I,J) + v_shlf(I,J-1)) - (v_shlf(I-1,J) + v_shlf(I-1,J-1))) / (2*G%dxT(i,j)) +! uy(i,j) = ((u_shlf(I,J) + u_shlf(I-1,J)) - (u_shlf(I,J-1) + u_shlf(I-1,J-1))) / (2*G%dyT(i,j)) +! vy(i,j) = ((v_shlf(I,J) + v_shlf(I-1,J)) - (v_shlf(I,J-1) + v_shlf(I-1,J-1))) / (2*G%dyT(i,j)) +! endif +! enddo +! enddo +! call pass_vector(ux, uy, G%domain, TO_ALL, BGRID_NE) +! call pass_vector(vx, vy, G%domain, TO_ALL, BGRID_NE) ! ux = ((u_shlf(I,J) + 0*u_shlf(I,J-1)) - (u_shlf(I-1,J) + 0*u_shlf(I-1,J-1))) / (G%dxT(i,j)) ! vx = ((v_shlf(I,J) + v_shlf(I,J-1)) - (v_shlf(I-1,J) + v_shlf(I-1,J-1))) / (2*G%dxT(i,j)) ! uy = ((u_shlf(I,J) + u_shlf(I-1,J)) - (u_shlf(I,J-1) + u_shlf(I-1,J-1))) / (2*G%dyT(i,j)) ! vy = ((v_shlf(I,J) + v_shlf(I-1,J)) - (v_shlf(I,J-1) + v_shlf(I-1,J-1))) / (2*G%dyT(i,j)) -!! ux = ((u_shlf(I,J) + u_shlf(I,J-1)) - (u_shlf(I-1,J) + u_shlf(I-1,J-1))) / (2*G%dxT(i,j)) -!! vx = ((v_shlf(I,J) + v_shlf(I,J-1)) - (v_shlf(I-1,J) + v_shlf(I-1,J-1))) / (2*G%dxT(i,j)) -!! uy = ((u_shlf(I,J) + u_shlf(I-1,J)) - (u_shlf(I,J-1) + u_shlf(I-1,J-1))) / (2*G%dyT(i,j)) -!! vy = ((v_shlf(I,J) + v_shlf(I-1,J)) - (v_shlf(I,J-1) + v_shlf(I-1,J-1))) / (2*G%dyT(i,j)) -! CS%ice_visc(i,j) = 0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * & -! (US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2))**((1.-n_g)/(2.*n_g)) + ux = ((u_shlf(I,J) + u_shlf(I,J-1)) - (u_shlf(I-1,J) + u_shlf(I-1,J-1))) / (2*G%dxT(i,j)) + vx = ((v_shlf(I,J) + v_shlf(I,J-1)) - (v_shlf(I-1,J) + v_shlf(I-1,J-1))) / (2*G%dxT(i,j)) + uy = ((u_shlf(I,J) + u_shlf(I-1,J)) - (u_shlf(I,J-1) + u_shlf(I-1,J-1))) / (2*G%dyT(i,j)) + vy = ((v_shlf(I,J) + v_shlf(I-1,J)) - (v_shlf(I,J-1) + v_shlf(I-1,J-1))) / (2*G%dyT(i,j)) + CS%ice_visc(i,j) = 0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * & + (US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2))**((1.-n_g)/(2.*n_g)) ! CS%ice_visc(i,j) =1e15*(G%areaT(i,j) * ISS%h_shelf(i,j)) ! constant viscocity for debugging ! umid = ((u_shlf(I,J) + u_shlf(I-1,J-1)) + (u_shlf(I,J-1) + u_shlf(I-1,J))) * 0.25 ! vmid = ((v_shlf(I,J) + v_shlf(I-1,J-1)) + (v_shlf(I,J-1) + v_shlf(I-1,J))) * 0.25 ! unorm = sqrt(umid**2 + vmid**2 + eps_min**2*(G%dxT(i,j)**2 + G%dyT(i,j)**2)) ! CS%basal_traction(i,j) = G%areaT(i,j) * CS%C_basal_friction * (US%L_T_to_m_s*unorm)**(CS%n_basal_fric-1) -! endif -! enddo -! enddo + endif + enddo + enddo - - do j=jsc-1,jec+1 - do i=isc-1,iec+1 - cnt = 0 - ux = 0 - uy = 0 - vx = 0 - vy = 0 - dxh = G%dxT(i,j) - dyh = G%dyT(i,j) +! do j=jsc-1,jec+1 +! do i=isc-1,iec+1 +!! do j=jsd+1,jed!-1 OVS 02/01/21 +!! do i=isd+1,ied!-1 OVS 02/01/21 - if (ISS%hmask(i,j) == 1) then ! we are inside the global computational bdry, at an ice-filled cell +! if (ISS%hmask(i,j) == 1) then +! CS%ice_visc(i,j) = 0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * & +! (US%s_to_T**2 * (ux(i,j)**2 + vy(i,j)**2 + ux(i,j)*vy(i,j) + 0.25*(uy(i,j)+vx(i,j))**2 + eps_min**2))**((1.-n_g)/(2.*n_g)) +! endif +! enddo +! enddo +! xquad(1) = .5 * (1-sqrt(1./3)) ; xquad(2) = .5 * (1+sqrt(1./3)) +! do j=jsc-1,jec+1 +! do i=isc-1,iec+1 +! cnt = 0 +! ux = 0 +! uy = 0 +! vx = 0 +! vy = 0 +! dxh = G%dxT(i,j) +! dyh = G%dyT(i,j) + +! if (ISS%hmask(i,j) == 1) then ! we are inside the global computational bdry, at an ice-filled cell + +! call bilinear_shape_fn_grid(G, i, j, Phi) +! do jq = 1,2 +! do iq = 1,2 + +! ux = u_shlf(I-1,J-1) * Phi(1,2*(jq-1)+iq) + & +! u_shlf(I,J-1) * Phi(3,2*(jq-1)+iq) + & +! u_shlf(I-1,J) * Phi(5,2*(jq-1)+iq) + & +! u_shlf(I,J) * Phi(7,2*(jq-1)+iq) + +! vx = v_shlf(I-1,J-1) * Phi(1,2*(jq-1)+iq) + & +! v_shlf(I,J-1) * Phi(3,2*(jq-1)+iq) + & +! v_shlf(I-1,J) * Phi(5,2*(jq-1)+iq) + & +! v_shlf(I,J) * Phi(7,2*(jq-1)+iq) + +! uy = u_shlf(I-1,J-1) * Phi(2,2*(jq-1)+iq) + & +! u_shlf(I,J-1) * Phi(4,2*(jq-1)+iq) + & +! u_shlf(I-1,J) * Phi(6,2*(jq-1)+iq) + & +! u_shlf(I,J) * Phi(8,2*(jq-1)+iq) + +! vy = v_shlf(I-1,J-1) * Phi(2,2*(jq-1)+iq) + & +! v_shlf(I,J-1) * Phi(4,2*(jq-1)+iq) + & +! v_shlf(I-1,J) * Phi(6,2*(jq-1)+iq) + & +! v_shlf(I,J) * Phi(8,2*(jq-1)+iq) +! enddo +! enddo ! calculate sx ! if ((i+i_off) == gisc) then ! at left computational bdry ! if (ISS%hmask(i+1,j) == 1) then @@ -2630,31 +2706,31 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) ! vx = 0 ! endif ! else ! interior - if (ISS%hmask(i+1,j) == 1) then - cnt = cnt+1 - ux = u_shlf(i+1,j) - vx = v_shlf(i+1,j) - else - ux = u_shlf(i,j) - vx = v_shlf(i,j) - endif - if (ISS%hmask(i-1,j) == 1) then - cnt = cnt+1 - ux = ux - u_shlf(i-1,j) - vx = vx - v_shlf(i-1,j) - else - ux = ux - u_shlf(i,j) - vx = vx - v_shlf(i,j) - endif - if (cnt == 0) then - ux = 0 - vx = 0 - else - ux = ux / (cnt * dxh) - vx = vx / (cnt * dxh) - endif -! endif - cnt = 0 +! if (ISS%hmask(i+1,j) == 1) then +! cnt = cnt+1 +! ux = u_shlf(i+1,j) +! vx = v_shlf(i+1,j) +! else +! ux = u_shlf(i,j) +! vx = v_shlf(i,j) +! endif +! if (ISS%hmask(i-1,j) == 1) then +! cnt = cnt+1 +! ux = ux - u_shlf(i-1,j) +! vx = vx - v_shlf(i-1,j) +! else +! ux = ux - u_shlf(i,j) +! vx = vx - v_shlf(i,j) +! endif +! if (cnt == 0) then +! ux = 0 +! vx = 0 +! else +! ux = ux / (cnt * dxh) +! vx = vx / (cnt * dxh) +! endif +!! endif +! cnt = 0 ! calculate sy, similarly ! if ((j+j_off) == gjsc) then ! at south computational bdry @@ -2673,70 +2749,72 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) ! vy = 0 ! endif ! else ! interior - if (ISS%hmask(i,j+1) == 1) then - cnt = cnt+1 - uy = u_shlf(i,j+1) - vy = v_shlf(i,j+1) - else - uy = u_shlf(i,j) - vy = v_shlf(i,j) - endif - if (ISS%hmask(i,j-1) == 1) then - cnt = cnt+1 - uy = uy - u_shlf(i,j-1) - vy = vy - v_shlf(i,j-1) - else - uy = uy - u_shlf(i,j) - vy = vy - v_shlf(i,j) - endif - if (cnt == 0) then - uy = 0 - vy = 0 - else - uy = uy / (cnt * dyh) - vy = vy / (cnt * dyh) - endif -! endif +! if (ISS%hmask(i,j+1) == 1) then +! cnt = cnt+1 +! uy = u_shlf(i,j+1) +! vy = v_shlf(i,j+1) +! else +! uy = u_shlf(i,j) +! vy = v_shlf(i,j) +! endif +! if (ISS%hmask(i,j-1) == 1) then +! cnt = cnt+1 +! uy = uy - u_shlf(i,j-1) +! vy = vy - v_shlf(i,j-1) +! else +! uy = uy - u_shlf(i,j) +! vy = vy - v_shlf(i,j) +! endif +! if (cnt == 0) then +! uy = 0 +! vy = 0 +! else +! uy = uy / (cnt * dyh) +! vy = vy / (cnt * dyh) +! endif +!! endif - ! SW vertex - if (ISS%hmask(I-1,J-1) == 1) then - eII(i-1,j-1) = eII(i-1,j-1)+.125*(US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2)) - endif +! ! SW vertex +! if (ISS%hmask(I-1,J-1) == 1) then +! eII(i-1,j-1) = eII(i-1,j-1)+.125*(US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2)) +! endif ! SE vertex - if (ISS%hmask(I,J-1) == 1) then - eII(i,j-1) = eII(i,j-1)+.125*(US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2)) +! if (ISS%hmask(I,J-1) == 1) then +! eII(i,j-1) = eII(i,j-1)+.125*(US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2)) ! CS%ice_visc(i,j-1) = CS%ice_visc(i,j-1)+.25*0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * & ! (US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2))**((1.-n_g)/(2.*n_g)) - endif +! endif ! NW vertex - if (ISS%hmask(I-1,J) == 1) then - eII(i-1,j) = eII(i-1,j)+.125*(US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2)) +! if (ISS%hmask(I-1,J) == 1) then +! eII(i-1,j) = eII(i-1,j)+.125*(US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2)) ! CS%ice_visc(i-1,j) = CS%ice_visc(i-1,j)+.25*0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * & ! (US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2))**((1.-n_g)/(2.*n_g)) - endif +! endif ! NE vertex - if (ISS%hmask(I,J) == 1) then - eII(i,j) = eII(i,j)+.25*(US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2)) - +! if (ISS%hmask(I,J) == 1) then +! eII(i,j) = eII(i,j)+.25*(US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2)) +! eII(i,j) = (US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2)) + ! CS%ice_visc(i,j) = CS%ice_visc(i,j)+.25*0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * & ! (US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2))**((1.-n_g)/(2.*n_g)) - endif - if (ISS%hmask(I+1,J+1) == 1) then - eII(i+1,j+1) = eII(i+1,j+1)+.125*(US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2)) - endif - if (ISS%hmask(I,J+1) == 1) then - eII(i,j+1) = eII(i,j+1)+.125*(US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2)) - endif - if (ISS%hmask(I+1,J) == 1) then - eII(i+1,j) = eII(i+1,j)+.125*(US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2)) - endif - endif - CS%ice_visc(i,j) =0.5 * Visc_coef*(G%areaT(i,j) * ISS%h_shelf(i,j))*eII(i,j)**((1.-n_g)/(2.*n_g)) -! CS%ice_visc(i,j) =1e15*(G%areaT(i,j) * ISS%h_shelf(i,j)) !constant viscosity for debugging - enddo - enddo +! endif +! if (ISS%hmask(I+1,J+1) == 1) then +! eII(i+1,j+1) = eII(i+1,j+1)+.125*(US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2)) +! endif +! if (ISS%hmask(I,J+1) == 1) then +! eII(i,j+1) = eII(i,j+1)+.125*(US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2)) +! endif +! if (ISS%hmask(I+1,J) == 1) then +! eII(i+1,j) = eII(i+1,j)+.125*(US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2)) +! endif +! CS%ice_visc(i,j) =0.5 * Visc_coef*(G%areaT(i,j) * ISS%h_shelf(i,j))*eII(i,j)**((1.-n_g)/(2.*n_g)) +! endif +! CS%ice_visc(i,j) =0.5 * Visc_coef*(G%areaT(i,j) * ISS%h_shelf(i,j))*eII(i,j)**((1.-n_g)/(2.*n_g)) + ! CS%ice_visc(i,j) =1e15*(G%areaT(i,j) * ISS%h_shelf(i,j)) !constant viscosity for debugging +! enddo +! enddo end subroutine calc_shelf_visc subroutine calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) @@ -2772,7 +2850,7 @@ subroutine calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) do j=jsd+1,jed!-1 OVS 02/01/21 do i=isd+1,ied!-1 OVS 02/01/21 - if (ISS%hmask(i,j) == 1) then + if ((ISS%hmask(i,j) == 1) .OR. (ISS%hmask(i,j) == 3)) then umid = ((u_shlf(I,J) + u_shlf(I-1,J-1)) + (u_shlf(I,J-1) + u_shlf(I-1,J))) * 0.25 vmid = ((v_shlf(I,J) + v_shlf(I-1,J-1)) + (v_shlf(I,J-1) + v_shlf(I-1,J))) * 0.25 unorm = sqrt(umid**2 + vmid**2 + eps_min**2*(G%dxT(i,j)**2 + G%dyT(i,j)**2)) diff --git a/src/ice_shelf/MOM_ice_shelf_initialize.F90 b/src/ice_shelf/MOM_ice_shelf_initialize.F90 index 532729c58c..3b6926e58f 100644 --- a/src/ice_shelf/MOM_ice_shelf_initialize.F90 +++ b/src/ice_shelf/MOM_ice_shelf_initialize.F90 @@ -366,13 +366,16 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b if (G%geoLatBu(i,j-1) == southlat) then !bot boundary if (len_stress == 0. .OR. G%geoLonCv(i,j) <= len_stress) then v_face_mask_bdry(i,j+1) = 0. - u_face_mask_bdry(i,j-1) = 3. !OVS 11/25/20 +! u_face_mask_bdry(i,j-1) = 3. !OVS 11/25/20 + u_face_mask_bdry(i,j) = 3. !OVS 11/25/20 + u_bdry_val(i,j) = 0. + v_bdry_val(i,j) = 0. else -! v_face_mask_bdry(i,j+1) = 1. - v_face_mask_bdry(i,j-1) = 3. !OVS 01/20/21 - u_face_mask_bdry(i,j-1) = 3. !OVS 11/25/20 -! u_bdry_val(i,j) = 0. -! v_bdry_val(i,j) = 0. !OVS 01/20/21 + v_face_mask_bdry(i,j+1) = 1. +! v_face_mask_bdry(i,j) = 3. !OVS 01/20/21 + u_face_mask_bdry(i,j) = 3. !OVS 11/25/20 + u_bdry_val(i,j) = 0. + v_bdry_val(i,j) = 0. !OVS 01/20/21 !hmask(i,j) = 0.0 !OVS 11/25/20 endif elseif (G%geoLatBu(i,j-1) == southlat+lenlat) then !top boundary From fdd83e6583415d6355d4bc9dfb3de421ddb66f9e Mon Sep 17 00:00:00 2001 From: Olga Sergienko Date: Mon, 22 Feb 2021 17:15:13 -0500 Subject: [PATCH 09/17] dynamic ice shelf with non-linear viscosity and evolving ice thickness due to sub-ice-shelf melting --- src/ice_shelf/MOM_ice_shelf_dynamics.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index 84605e3092..c63a42beaf 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -714,7 +714,7 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled coupled_GL = .false. if (present(ocean_mass) .and. present(coupled_grounding)) coupled_GL = coupled_grounding -! call ice_shelf_advect(CS, ISS, G, time_step, Time) !OVS 02/08/21 + call ice_shelf_advect(CS, ISS, G, time_step, Time) !OVS 02/08/21 CS%elapsed_velocity_time = CS%elapsed_velocity_time + time_step if (CS%elapsed_velocity_time >= CS%velocity_update_time_step) update_ice_vel = .true. From 32cfe35a44a40fd2947c089946b525d795b670b5 Mon Sep 17 00:00:00 2001 From: Olga Sergienko Date: Wed, 24 Feb 2021 15:41:54 -0500 Subject: [PATCH 10/17] modified MOM_ice_shelf_initialize for testing with viscosity from a file --- src/ice_shelf/MOM_ice_shelf_dynamics.F90 | 38 ++++++----- src/ice_shelf/MOM_ice_shelf_initialize.F90 | 77 ++++++++++++++++++++++ 2 files changed, 99 insertions(+), 16 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index c63a42beaf..5e6ba60a1a 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -23,7 +23,7 @@ module MOM_ice_shelf_dynamics use MOM_ice_shelf_state, only : ice_shelf_state use MOM_coms, only : reproducing_sum, sum_across_PEs, max_across_PEs, min_across_PEs use MOM_checksums, only : hchksum, qchksum -use MOM_ice_shelf_initialize, only : initialize_ice_shelf_boundary_channel !OVS intializing b.c.s +use MOM_ice_shelf_initialize, only : initialize_ice_shelf_boundary_channel,initialize_ice_flow_from_file !OVS intializing b.c.s implicit none ; private @@ -535,7 +535,6 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ CS%thickness_bdry_val, ISS%hmask, ISS%h_shelf, G, & ! CS%flux_bdry, & US, param_file ) !OVS initialize b.c.s - call pass_var(ISS%hmask, G%domain) call pass_var(CS%h_bdry_val, G%domain) call pass_var(CS%thickness_bdry_val, G%domain) @@ -545,6 +544,9 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ call pass_var(CS%v_face_mask_bdry, G%domain) ! call init_boundary_values(CS, G, time, ISS%hmask, CS%input_flux, CS%input_thickness, new_sim) call update_velocity_masks(CS, G, ISS%hmask, CS%umask, CS%vmask, CS%u_face_mask, CS%v_face_mask) +! call initialize_ice_flow_from_file(CS%u_shelf, CS%v_shelf,CS%ice_visc,CS%ground_frac, ISS%hmask,ISS%h_shelf, & +! G, US, param_file) !spacially variable viscosity from a file for debugging +! call pass_var(CS%ice_visc, G%domain) ! if (new_sim) then ! call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: initialize ice velocity.") ! call update_OD_ffrac_uncoupled(CS, G, ISS%h_shelf(:,:)) @@ -713,7 +715,7 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled coupled_GL = .false. if (present(ocean_mass) .and. present(coupled_grounding)) coupled_GL = coupled_grounding - +! call ice_shelf_advect(CS, ISS, G, time_step, Time) !OVS 02/08/21 CS%elapsed_velocity_time = CS%elapsed_velocity_time + time_step if (CS%elapsed_velocity_time >= CS%velocity_update_time_step) update_ice_vel = .true. @@ -946,7 +948,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, ite call bilinear_shape_fn_grid(G, i, j, Phi(:,:,i,j)) enddo ; enddo - call calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) + call calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) !OVS 02/24/21 call pass_var(CS%ice_visc, G%domain) ! call pass_vector(CS%ice_visc, G%domain, TO_ALL, BGRID_NE) !OVS 02/11/21 call calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) @@ -1000,7 +1002,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, ite write(mesg,*) "ice_shelf_solve_outer: linear solve done in ",iters," iterations" call MOM_mesg(mesg, 5) - call calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) + call calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) !OVS 02/24/21 call pass_var(CS%ice_visc, G%domain) call calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) call pass_var(CS%basal_traction, G%domain) @@ -1351,7 +1353,7 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H if (cg_halo == 0) then ! pass vectors call pass_vector(Du, Dv, G%domain, TO_ALL, BGRID_NE) -! call pass_vector(u_shlf, v_shlf, G%domain, TO_ALL, BGRID_NE) + call pass_vector(u_shlf, v_shlf, G%domain, TO_ALL, BGRID_NE) call pass_var(u_shlf, G%domain) call pass_var(v_shlf, G%domain) call pass_vector(Ru, Rv, G%domain, TO_ALL, BGRID_NE) @@ -2604,10 +2606,10 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) ! eII(:,:) = (US%s_to_T**2 * (eps_min**2)) Visc_coef = US%kg_m2s_to_RZ_T*US%m_to_L*US%Z_to_L*(CS%A_glen_isothermal)**(-1./CS%n_glen) !OVS '-' in the exponent ! call pass_vector(u_shlf, v_shlf, G%domain, TO_ALL, BGRID_NE) -! do j=jsc-1,jec+1 -! do i=isc-1,iec+1 - do j=jsd+1,jed-1 !OVS 02/01/21 - do i=isd+1,ied-1 !OVS 02/01/21 + do j=jsc-0*1,jec+0*1 + do i=isc-0*1,iec+0*1 +! do j=jsd+1,jed-1 !OVS 02/01/21 +! do i=isd+1,ied-1 !OVS 02/01/21 if ((ISS%hmask(i,j) == 1) .OR. (ISS%hmask(i,j) == 3)) then ! ux(i,j) = ((u_shlf(I,J) + u_shlf(I,J-1)) - (u_shlf(I-1,J) + u_shlf(I-1,J-1))) / (2*G%dxT(i,j)) @@ -2619,14 +2621,18 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) ! enddo ! call pass_vector(ux, uy, G%domain, TO_ALL, BGRID_NE) ! call pass_vector(vx, vy, G%domain, TO_ALL, BGRID_NE) -! ux = ((u_shlf(I,J) + 0*u_shlf(I,J-1)) - (u_shlf(I-1,J) + 0*u_shlf(I-1,J-1))) / (G%dxT(i,j)) + ux = ((u_shlf(I,J) + u_shlf(I,J-1) + u_shlf(I,J+1)) - & + (u_shlf(I-1,J) + u_shlf(I-1,J-1) + u_shlf(I-1,J+1))) / (3*G%dxT(i,j)) + vx = ((v_shlf(I,J) + v_shlf(I,J-1) + v_shlf(I,J+1)) - & + (v_shlf(I-1,J) + v_shlf(I-1,J-1) + v_shlf(I-1,J+1))) / (3*G%dxT(i,j)) + uy = ((u_shlf(I,J) + u_shlf(I-1,J) + u_shlf(I+1,J)) - & + (u_shlf(I,J-1) + u_shlf(I-1,J-1) + u_shlf(I+1,J-1))) / (3*G%dyT(i,j)) + vy = ((v_shlf(I,J) + v_shlf(I-1,J)+ v_shlf(I+1,J)) - & + (v_shlf(I,J-1) + v_shlf(I-1,J-1)+ v_shlf(I+1,J-1))) / (3*G%dyT(i,j)) +! ux = ((u_shlf(I,J) + u_shlf(I,J-1)) - (u_shlf(I-1,J) + u_shlf(I-1,J-1))) / (2*G%dxT(i,j)) ! vx = ((v_shlf(I,J) + v_shlf(I,J-1)) - (v_shlf(I-1,J) + v_shlf(I-1,J-1))) / (2*G%dxT(i,j)) ! uy = ((u_shlf(I,J) + u_shlf(I-1,J)) - (u_shlf(I,J-1) + u_shlf(I-1,J-1))) / (2*G%dyT(i,j)) -! vy = ((v_shlf(I,J) + v_shlf(I-1,J)) - (v_shlf(I,J-1) + v_shlf(I-1,J-1))) / (2*G%dyT(i,j)) - ux = ((u_shlf(I,J) + u_shlf(I,J-1)) - (u_shlf(I-1,J) + u_shlf(I-1,J-1))) / (2*G%dxT(i,j)) - vx = ((v_shlf(I,J) + v_shlf(I,J-1)) - (v_shlf(I-1,J) + v_shlf(I-1,J-1))) / (2*G%dxT(i,j)) - uy = ((u_shlf(I,J) + u_shlf(I-1,J)) - (u_shlf(I,J-1) + u_shlf(I-1,J-1))) / (2*G%dyT(i,j)) - vy = ((v_shlf(I,J) + v_shlf(I-1,J)) - (v_shlf(I,J-1) + v_shlf(I-1,J-1))) / (2*G%dyT(i,j)) +! vy = ((v_shlf(I,J) + v_shlf(I-1,J)) - (v_shlf(I,J-1) + v_shlf(I-1,J-1))) / (2*G%dyT(i,j)) CS%ice_visc(i,j) = 0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * & (US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2))**((1.-n_g)/(2.*n_g)) ! CS%ice_visc(i,j) =1e15*(G%areaT(i,j) * ISS%h_shelf(i,j)) ! constant viscocity for debugging diff --git a/src/ice_shelf/MOM_ice_shelf_initialize.F90 b/src/ice_shelf/MOM_ice_shelf_initialize.F90 index 3b6926e58f..f2e01c461b 100644 --- a/src/ice_shelf/MOM_ice_shelf_initialize.F90 +++ b/src/ice_shelf/MOM_ice_shelf_initialize.F90 @@ -19,6 +19,7 @@ module MOM_ice_shelf_initialize !MJHpublic initialize_ice_shelf_boundary, initialize_ice_thickness public initialize_ice_thickness public initialize_ice_shelf_boundary_channel +public initialize_ice_flow_from_file ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with @@ -589,4 +590,80 @@ end subroutine initialize_ice_shelf_boundary_channel !END MJH end subroutine initialize_ice_shelf_boundary_channel +subroutine initialize_ice_flow_from_file(u_shelf, v_shelf,ice_visc,float_cond, hmask,h_shelf, G, US, PF) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(inout) :: u_shelf !< The ice shelf u velocity [Z ~> m T ~>s]. + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(inout) :: v_shelf !< The ice shelf v velocity [Z ~> m T ~> s]. + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(inout) :: ice_visc !< The ice shelf viscosity [Pa ~> m T ~> s]. + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(inout) :: float_cond !< An array indicating where the ice + !! shelf is floating: 0 if floating, 1 if not. + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(in) :: hmask !< A mask indicating which tracer points are + !! partly or fully covered by an ice-shelf + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(in) :: h_shelf !< A mask indicating which tracer points are + !! partly or fully covered by an ice-shelf + type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors + type(param_file_type), intent(in) :: PF !< A structure to parse for run-time parameters + + ! This subroutine reads ice thickness and area from a file and puts it into + ! h_shelf [Z ~> m] and area_shelf_h [L2 ~> m2] (and dimensionless) and updates hmask + character(len=200) :: filename,vel_file,inputdir ! Strings for file/path + character(len=200) :: ushelf_varname, vshelf_varname, ice_visc_varname, floatfr_varname ! Variable name in file + character(len=40) :: mdl = "initialize_ice_velocity_from_file" ! This subroutine's name. + integer :: i, j, isc, jsc, iec, jec + real :: len_sidestress, mask, udh + + call MOM_mesg(" MOM_ice_shelf_init_profile.F90, initialize_velocity_from_file: reading velocity") + + call get_param(PF, mdl, "INPUTDIR", inputdir, default=".") + inputdir = slasher(inputdir) + call get_param(PF, mdl, "ICE_VELOCITY_FILE", vel_file, & + "The file from which the velocity is read.", & + default="ice_shelf_vel.nc") + call get_param(PF, mdl, "LEN_SIDE_STRESS", len_sidestress, & + "position past which shelf sides are stress free.", & + default=0.0, units="axis_units") + + filename = trim(inputdir)//trim(vel_file) + call log_param(PF, mdl, "INPUTDIR/THICKNESS_FILE", filename) + call get_param(PF, mdl, "ICE_U_VEL_VARNAME", ushelf_varname, & + "The name of the thickness variable in ICE_VELOCITY_FILE.", & + default="u_shelf") + call get_param(PF, mdl, "ICE_V_VEL_VARNAME", vshelf_varname, & + "The name of the thickness variable in ICE_VELOCITY_FILE.", & + default="v_shelf") + call get_param(PF, mdl, "ICE_VISC_VARNAME", ice_visc_varname, & + "The name of the thickness variable in ICE_VELOCITY_FILE.", & + default="viscosity") + + if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, & + " initialize_ice_shelf_velocity_from_file: Unable to open "//trim(filename)) + + !hmask_varname = "hmask" + floatfr_varname = "float_frac" + +! call MOM_read_data(filename, trim(ushelf_varname), u_shelf, G%Domain, scale=1.0) !/(365.0*86400.0)) +! call MOM_read_data(filename,trim(vshelf_varname), v_shelf, G%Domain, scale=1.0) !/(365.0*86400.0)) + call MOM_read_data(filename,trim(ice_visc_varname), ice_visc, G%Domain, scale=1.0) !*(365.0*86400.0)) + call MOM_read_data(filename,trim(floatfr_varname), float_cond, G%Domain, scale=1.) +! call get_param(PF, mdl, "ICE_BOUNDARY_CONFIG", config, & +! "This specifies how the ice domain boundary is specified", & +! fail_if_missing=.true.) + + isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec + + do j=jsc,jec + do i=isc,iec + if (hmask(i,j) == 1.) then + ice_visc(i,j) = ice_visc(i,j) * (G%areaT(i,j) * h_shelf(i,j)) + endif + enddo + enddo + +end subroutine initialize_ice_flow_from_file end module MOM_ice_shelf_initialize From 5483bfed1243fd765b34059f96a68f4a0dc2b5ed Mon Sep 17 00:00:00 2001 From: Olga Sergienko Date: Wed, 3 Mar 2021 10:46:48 -0500 Subject: [PATCH 11/17] Cleaned initialize_ice_shelf_boundary_channel --- src/ice_shelf/MOM_ice_shelf_initialize.F90 | 33 ++++++++++------------ 1 file changed, 15 insertions(+), 18 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf_initialize.F90 b/src/ice_shelf/MOM_ice_shelf_initialize.F90 index f2e01c461b..7ba1ab7076 100644 --- a/src/ice_shelf/MOM_ice_shelf_initialize.F90 +++ b/src/ice_shelf/MOM_ice_shelf_initialize.F90 @@ -265,6 +265,8 @@ subroutine initialize_ice_thickness_channel(h_shelf, area_shelf_h, hmask, G, US, endif ; enddo end subroutine initialize_ice_thickness_channel + +!> Initialize ice shelf boundary conditions for a channel configuration subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_bdry, & u_flux_bdry_val, v_flux_bdry_val, u_bdry_val, v_bdry_val, u_shelf, v_shelf, h_bdry_val, & thickness_bdry_val, hmask, h_shelf, G,& @@ -309,7 +311,7 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b real :: input_thick ! The input ice shelf thickness [Z ~> m] ! real :: input_flux ! The input ice flux per unit length [L Z T-1 ~> m2 s-1] real :: input_vel ! The input ice velocity per [L Z T-1 ~> m s-1] - real :: lenlat, len_stress, westlon, lenlon, southlat + real :: lenlat, len_stress, westlon, lenlon, southlat ! The input positions of the channel boundarises call get_param(PF, mdl, "LENLAT", lenlat, fail_if_missing=.true.) @@ -352,11 +354,11 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b ! else hmask(i+1,j) = 3.0 ! hmask(i,j) = 3.0 - h_bdry_val(i+1,j) = h_shelf(i+1,j) !OVS 11/10/20 !input_thick + h_bdry_val(i+1,j) = h_shelf(i+1,j) ! h_bdry_val(i,j) = h_shelf(i,j) thickness_bdry_val(i+1,j) = h_bdry_val(i+0*1,j) u_face_mask_bdry(i+1,j) = 3.0 - u_bdry_val(i+1,j) = input_vel*(1-16.0*((G%geoLatBu(i-1,j)/lenlat-0.5))**4) !OVS 11/09/20 U b.c. + u_bdry_val(i+1,j) = input_vel*(1-16.0*((G%geoLatBu(i-1,j)/lenlat-0.5))**4) !velocity distribution ! u_bdry_val(i+1,j) = (1 - ((G%geoLatBu(i,j) - 0.5*lenlat)*2./lenlat)**2) * & ! 1.5 * input_flux / input_thick ! endif @@ -367,28 +369,26 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b if (G%geoLatBu(i,j-1) == southlat) then !bot boundary if (len_stress == 0. .OR. G%geoLonCv(i,j) <= len_stress) then v_face_mask_bdry(i,j+1) = 0. -! u_face_mask_bdry(i,j-1) = 3. !OVS 11/25/20 - u_face_mask_bdry(i,j) = 3. !OVS 11/25/20 +! u_face_mask_bdry(i,j-1) = 3. + u_face_mask_bdry(i,j) = 3. u_bdry_val(i,j) = 0. v_bdry_val(i,j) = 0. else v_face_mask_bdry(i,j+1) = 1. -! v_face_mask_bdry(i,j) = 3. !OVS 01/20/21 - u_face_mask_bdry(i,j) = 3. !OVS 11/25/20 + u_face_mask_bdry(i,j) = 3. u_bdry_val(i,j) = 0. - v_bdry_val(i,j) = 0. !OVS 01/20/21 - !hmask(i,j) = 0.0 !OVS 11/25/20 + v_bdry_val(i,j) = 0. endif elseif (G%geoLatBu(i,j-1) == southlat+lenlat) then !top boundary if (len_stress == 0. .OR. G%geoLonCv(i,j) <= len_stress) then v_face_mask_bdry(i,j-1) = 0. - u_face_mask_bdry(i,j-1) = 3. !OVS 11/25/20 + u_face_mask_bdry(i,j-1) = 3. else ! v_face_mask_bdry(i,j-1) = 1. - v_face_mask_bdry(i,j-1) = 3. !OVS 01/20/21 - u_face_mask_bdry(i,j-1) = 3. !OVS 11/25/20 - !u_bdry_val(i,j) = 0. !OVS 11/25/20 - !hmask(i,j) = 0.0 !OVS 11/25/20 + v_face_mask_bdry(i,j-1) = 3. + u_face_mask_bdry(i,j-1) = 3. + !u_bdry_val(i,j) = 0. + !hmask(i,j) = 0.0 endif endif @@ -400,9 +400,6 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b enddo enddo -! call pass_var(hmask, G%domain) -! call pass_var(h_bdry_val, G%domain) -! call pass_var(thickness_bdry_val, G%domain) ! if (.not. G%symmetric) then @@ -590,6 +587,7 @@ end subroutine initialize_ice_shelf_boundary_channel !END MJH end subroutine initialize_ice_shelf_boundary_channel +!> Initialize ice shelf flow from file subroutine initialize_ice_flow_from_file(u_shelf, v_shelf,ice_visc,float_cond, hmask,h_shelf, G, US, PF) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure real, dimension(SZDI_(G),SZDJ_(G)), & @@ -644,7 +642,6 @@ subroutine initialize_ice_flow_from_file(u_shelf, v_shelf,ice_visc,float_cond, h if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, & " initialize_ice_shelf_velocity_from_file: Unable to open "//trim(filename)) - !hmask_varname = "hmask" floatfr_varname = "float_frac" ! call MOM_read_data(filename, trim(ushelf_varname), u_shelf, G%Domain, scale=1.0) !/(365.0*86400.0)) From 9aa75c8691c8c4e321a5894875d04a519bc4f0a1 Mon Sep 17 00:00:00 2001 From: Olga Sergienko Date: Wed, 3 Mar 2021 11:16:53 -0500 Subject: [PATCH 12/17] Modified MOM_ice_shelf_initialize.F90 --- src/ice_shelf/MOM_ice_shelf_initialize.F90 | 43 ++++++++++------------ 1 file changed, 19 insertions(+), 24 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf_initialize.F90 b/src/ice_shelf/MOM_ice_shelf_initialize.F90 index 7ba1ab7076..f9f31a373e 100644 --- a/src/ice_shelf/MOM_ice_shelf_initialize.F90 +++ b/src/ice_shelf/MOM_ice_shelf_initialize.F90 @@ -302,14 +302,12 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b !! partly or fully covered by an ice-shelf real, dimension(SZDI_(G),SZDJ_(G)), & intent(inout) :: h_shelf !< Ice-shelf thickness OVS 11/10/20 -! logical, intent(in) :: flux_bdry !< If true, use mass fluxes as the boundary value. type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors type(param_file_type), intent(in) :: PF !< A structure to parse for run-time parameters character(len=40) :: mdl = "initialize_ice_shelf_boundary_channel" ! This subroutine's name. integer :: i, j, isd, jsd, is, js, iegq, jegq, giec, gjec, gisc, gjsc,gisd,gjsd, isc, jsc, iec, jec, ied, jed - real :: input_thick ! The input ice shelf thickness [Z ~> m] -! real :: input_flux ! The input ice flux per unit length [L Z T-1 ~> m2 s-1] + real :: input_thick ! The input ice shelf thickness [Z ~> m] real :: input_vel ! The input ice velocity per [L Z T-1 ~> m s-1] real :: lenlat, len_stress, westlon, lenlon, southlat ! The input positions of the channel boundarises @@ -341,27 +339,27 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b gisc = G%Domain%nihalo ; gjsc = G%Domain%njhalo giec = G%Domain%niglobal+gisc ; gjec = G%Domain%njglobal+gjsc -!-----------b.c.s based on geopositions ----------------- -! do j=jsc-1,jec+1 + !---------b.c.s based on geopositions ----------------- + ! do j=jsc-1,jec+1 do j=jsc-0*1,jec+1 do i=isc-1,iec+1 ! upstream boundary - set either dirichlet or flux condition if (G%geoLonBu(i,j) == westlon) then - ! if (flux_bdry) then - ! u_face_mask_bdry(i-1,j) = 4.0 - ! u_flux_bdry_val(i-1,j) = input_flux - ! else + ! if (flux_bdry) then + ! u_face_mask_bdry(i-1,j) = 4.0 + ! u_flux_bdry_val(i-1,j) = input_flux + ! else hmask(i+1,j) = 3.0 -! hmask(i,j) = 3.0 + ! hmask(i,j) = 3.0 h_bdry_val(i+1,j) = h_shelf(i+1,j) -! h_bdry_val(i,j) = h_shelf(i,j) + ! h_bdry_val(i,j) = h_shelf(i,j) thickness_bdry_val(i+1,j) = h_bdry_val(i+0*1,j) u_face_mask_bdry(i+1,j) = 3.0 u_bdry_val(i+1,j) = input_vel*(1-16.0*((G%geoLatBu(i-1,j)/lenlat-0.5))**4) !velocity distribution - ! u_bdry_val(i+1,j) = (1 - ((G%geoLatBu(i,j) - 0.5*lenlat)*2./lenlat)**2) * & - ! 1.5 * input_flux / input_thick - ! endif + ! u_bdry_val(i+1,j) = (1 - ((G%geoLatBu(i,j) - 0.5*lenlat)*2./lenlat)**2) * & + ! 1.5 * input_flux / input_thick + ! endif endif @@ -369,7 +367,7 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b if (G%geoLatBu(i,j-1) == southlat) then !bot boundary if (len_stress == 0. .OR. G%geoLonCv(i,j) <= len_stress) then v_face_mask_bdry(i,j+1) = 0. -! u_face_mask_bdry(i,j-1) = 3. + ! u_face_mask_bdry(i,j-1) = 3. u_face_mask_bdry(i,j) = 3. u_bdry_val(i,j) = 0. v_bdry_val(i,j) = 0. @@ -384,7 +382,7 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b v_face_mask_bdry(i,j-1) = 0. u_face_mask_bdry(i,j-1) = 3. else -! v_face_mask_bdry(i,j-1) = 1. + !v_face_mask_bdry(i,j-1) = 1. v_face_mask_bdry(i,j-1) = 3. u_face_mask_bdry(i,j-1) = 3. !u_bdry_val(i,j) = 0. @@ -398,10 +396,7 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b endif enddo - enddo - - - + enddo ! if (.not. G%symmetric) then !! do j=G%jsd,G%jed !! do i=G%isd,G%ied @@ -623,7 +618,7 @@ subroutine initialize_ice_flow_from_file(u_shelf, v_shelf,ice_visc,float_cond, h call get_param(PF, mdl, "ICE_VELOCITY_FILE", vel_file, & "The file from which the velocity is read.", & default="ice_shelf_vel.nc") - call get_param(PF, mdl, "LEN_SIDE_STRESS", len_sidestress, & + call get_param(PF, mdl, "LEN_SIDE_STRESS", len_sidestress, & "position past which shelf sides are stress free.", & default=0.0, units="axis_units") @@ -642,12 +637,12 @@ subroutine initialize_ice_flow_from_file(u_shelf, v_shelf,ice_visc,float_cond, h if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, & " initialize_ice_shelf_velocity_from_file: Unable to open "//trim(filename)) - floatfr_varname = "float_frac" + floatfr_varname = "float_frac" ! call MOM_read_data(filename, trim(ushelf_varname), u_shelf, G%Domain, scale=1.0) !/(365.0*86400.0)) ! call MOM_read_data(filename,trim(vshelf_varname), v_shelf, G%Domain, scale=1.0) !/(365.0*86400.0)) - call MOM_read_data(filename,trim(ice_visc_varname), ice_visc, G%Domain, scale=1.0) !*(365.0*86400.0)) - call MOM_read_data(filename,trim(floatfr_varname), float_cond, G%Domain, scale=1.) + call MOM_read_data(filename,trim(ice_visc_varname), ice_visc, G%Domain, scale=1.0) !*(365.0*86400.0)) + call MOM_read_data(filename,trim(floatfr_varname), float_cond, G%Domain, scale=1.) ! call get_param(PF, mdl, "ICE_BOUNDARY_CONFIG", config, & ! "This specifies how the ice domain boundary is specified", & ! fail_if_missing=.true.) From 2232fa2882806f8ebabb324d93cc580c107bb17a Mon Sep 17 00:00:00 2001 From: Olga Sergienko Date: Wed, 3 Mar 2021 15:51:11 -0500 Subject: [PATCH 13/17] corrected style errors in MOM_ice_shelf.F90; MOM_ice_shelf_dynamics.F90; MOM_ice_shelf_initialize.F90 --- src/ice_shelf/MOM_ice_shelf.F90 | 2 +- src/ice_shelf/MOM_ice_shelf_dynamics.F90 | 295 +++------------------ src/ice_shelf/MOM_ice_shelf_initialize.F90 | 41 +-- 3 files changed, 63 insertions(+), 275 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 5663b326b7..5d2bc88b4c 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -719,7 +719,7 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step, CS) call hchksum(ISS%mass_shelf, "mass_shelf after change thickness using melt", G%HI, haloshift=0, & scale=US%RZ_to_kg_m2) endif - endif !OVS 12/10/20 + endif !OVS 12/10/20 if (CS%debug) call MOM_forcing_chksum("Before add shelf flux", fluxes, G, CS%US, haloshift=0) diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index 5e6ba60a1a..8640902989 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -23,7 +23,7 @@ module MOM_ice_shelf_dynamics use MOM_ice_shelf_state, only : ice_shelf_state use MOM_coms, only : reproducing_sum, sum_across_PEs, max_across_PEs, min_across_PEs use MOM_checksums, only : hchksum, qchksum -use MOM_ice_shelf_initialize, only : initialize_ice_shelf_boundary_channel,initialize_ice_flow_from_file !OVS intializing b.c.s +use MOM_ice_shelf_initialize, only : initialize_ice_shelf_boundary_channel,initialize_ice_flow_from_file implicit none ; private @@ -47,7 +47,7 @@ module MOM_ice_shelf_dynamics real, pointer, dimension(:,:) :: taudx_shelf => NULL() !< the driving stress of the ice shelf/sheet !! on q-points (C grid) [Pa ~> Pa] real, pointer, dimension(:,:) :: taudy_shelf => NULL() !< the meridional stress of the ice shelf/sheet - !! on q-points (C grid) [Pa ~> Pa] + !! on q-points (C grid) [Pa ~> Pa] real, pointer, dimension(:,:) :: u_face_mask => NULL() !< mask for velocity boundary conditions on the C-grid !! u-face - this is because the FEM cares about FACES THAT GET INTEGRATED OVER, !! not vertices. Will represent boundary conditions on computational boundary @@ -263,10 +263,10 @@ subroutine register_ice_shelf_dyn_restarts(G, param_file, CS, restart_CS) "ice sheet/shelf v-velocity", "m s-1", hor_grid='Bu') call register_restart_field(CS%t_shelf, "t_shelf", .true., restart_CS, & "ice sheet/shelf vertically averaged temperature", "deg C") - call register_restart_field(CS%taudx_shelf, "taudx_shelf", .true., restart_CS, & !OVS 02/8/21 + call register_restart_field(CS%taudx_shelf, "taudx_shelf", .true., restart_CS, & "ice sheet/shelf taudx-driving stress", "kPa") - call register_restart_field(CS%taudy_shelf, "taudy_shelf", .true., restart_CS, & !OVS 02/08/21 - "ice sheet/shelf taudy-driving stress", "kPa") + call register_restart_field(CS%taudy_shelf, "taudy_shelf", .true., restart_CS, & + "ice sheet/shelf taudy-driving stress", "kPa") call register_restart_field(CS%OD_av, "OD_av", .true., restart_CS, & "Average open ocean depth in a cell","m") call register_restart_field(CS%ground_frac, "ground_frac", .true., restart_CS, & @@ -376,23 +376,20 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ call get_param(param_file, mdl, "A_GLEN_ISOTHERM", CS%A_glen_isothermal, & "Ice viscosity parameter in Glen's Law", & - units="Pa-3 s-1", default=2.2261e-25, scale=1.0) !OVS change units to Pa-3 s-1 -! units="Pa-3 yr-1", default=9.461e-18, scale=1.0/(365.0*86400.0)) + units="Pa-3 s-1", default=2.2261e-25, scale=1.0) ! This default is equivalent to 3.0001e-25 Pa-3 s-1, appropriate at about -10 C. call get_param(param_file, mdl, "GLEN_EXPONENT", CS%n_glen, & "nonlinearity exponent in Glen's Law", & units="none", default=3.) call get_param(param_file, mdl, "MIN_STRAIN_RATE_GLEN", CS%eps_glen_min, & "min. strain rate to avoid infinite Glen's law viscosity", & - units="s-1", default=1.e-19, scale=US%T_to_s) !OVS change units to s-1 - !units="a-1", default=1.e-12, scale=US%T_to_s/(365.0*86400.0)) + units="s-1", default=1.e-19, scale=US%T_to_s) call get_param(param_file, mdl, "BASAL_FRICTION_EXP", CS%n_basal_fric, & "Exponent in sliding law \tau_b = C u^(n_basal_fric)", & units="none", fail_if_missing=.true.) call get_param(param_file, mdl, "BASAL_FRICTION_COEFF", CS%C_basal_friction, & "Coefficient in sliding law \tau_b = C u^(n_basal_fric)", & - units="Pa (m s-1)^(n_basal_fric)", scale=US%kg_m2s_to_RZ_T**CS%n_basal_fric, & ! OVS change units to s-1 - !units="Pa (m yr-1)-(n_basal_fric)", scale=US%kg_m2s_to_RZ_T*((365.0*86400.0)**CS%n_basal_fric), & + units="Pa (m s-1)^(n_basal_fric)", scale=US%kg_m2s_to_RZ_T**CS%n_basal_fric, & fail_if_missing=.true.) call get_param(param_file, mdl, "DENSITY_ICE", CS%density_ice, & "A typical density of ice.", units="kg m-3", default=917.0, scale=US%kg_m3_to_R) @@ -416,7 +413,7 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ call get_param(param_file, mdl, "CALVE_TO_MASK", CS%calve_to_mask, & "If true, do not allow an ice shelf where prohibited by a mask.", & default=.false.) - + endif call get_param(param_file, mdl, "MIN_THICKNESS_SIMPLE_CALVE", CS%min_thickness_simple_calve, & "Min thickness rule for the VERY simple calving law",& @@ -424,7 +421,6 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ ! Allocate memory in the ice shelf dynamics control structure that was not ! previously allocated for registration for restarts. - ! OVS vertically integrated Temperature if (active_shelf_dynamics) then ! DNG @@ -533,68 +529,37 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ CS%u_flux_bdry_val, CS%v_flux_bdry_val, CS%u_bdry_val, CS%v_bdry_val, CS%u_shelf, CS%v_shelf,& CS%h_bdry_val, & CS%thickness_bdry_val, ISS%hmask, ISS%h_shelf, G, & -! CS%flux_bdry, & - US, param_file ) !OVS initialize b.c.s + US, param_file ) call pass_var(ISS%hmask, G%domain) call pass_var(CS%h_bdry_val, G%domain) - call pass_var(CS%thickness_bdry_val, G%domain) - call pass_var(CS%u_bdry_val, G%domain) - call pass_var(CS%v_bdry_val, G%domain) - call pass_var(CS%u_face_mask_bdry, G%domain) - call pass_var(CS%v_face_mask_bdry, G%domain) -! call init_boundary_values(CS, G, time, ISS%hmask, CS%input_flux, CS%input_thickness, new_sim) + call pass_var(CS%thickness_bdry_val, G%domain) + call pass_var(CS%u_bdry_val, G%domain) + call pass_var(CS%v_bdry_val, G%domain) + call pass_var(CS%u_face_mask_bdry, G%domain) + call pass_var(CS%v_face_mask_bdry, G%domain) + !call init_boundary_values(CS, G, time, ISS%hmask, CS%input_flux, CS%input_thickness, new_sim) call update_velocity_masks(CS, G, ISS%hmask, CS%umask, CS%vmask, CS%u_face_mask, CS%v_face_mask) -! call initialize_ice_flow_from_file(CS%u_shelf, CS%v_shelf,CS%ice_visc,CS%ground_frac, ISS%hmask,ISS%h_shelf, & -! G, US, param_file) !spacially variable viscosity from a file for debugging -! call pass_var(CS%ice_visc, G%domain) -! if (new_sim) then -! call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: initialize ice velocity.") -! call update_OD_ffrac_uncoupled(CS, G, ISS%h_shelf(:,:)) -!! call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf, iters, Time) -! call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) !OVS 02/08/21 -! if (CS%id_u_shelf > 0) call post_data(CS%id_u_shelf, CS%u_shelf, CS%diag) -! if (CS%id_v_shelf > 0) call post_data(CS%id_v_shelf, CS%v_shelf,CS%diag) -! if (CS%id_taudx_shelf > 0) call post_data(CS%id_taudx_shelf, CS%taudx_shelf,CS%diag) -! if (CS%id_taudy_shelf > 0) call post_data(CS%id_taudy_shelf, CS%taudy_shelf,CS%diag) -! if (CS%id_visc_shelf > 0) call post_data(CS%id_visc_shelf, CS%ice_visc,CS%diag) -! endif ! Register diagnostics. -! CS%id_u_shelf = register_diag_field('ocean_model','u_shelf',CS%diag%axesCu1, Time, & -! 'x-velocity of ice', 'm yr-1', conversion=365.0*86400.0*US%L_T_to_m_s) CS%id_u_shelf = register_diag_field('ice_shelf_model','u_shelf',CS%diag%axesCu1, Time, & 'x-velocity of ice', 'm yr-1', conversion=365.0*86400.0*US%L_T_to_m_s) -! CS%id_v_shelf = register_diag_field('ocean_model','v_shelf',CS%diag%axesCv1, Time, & -! 'y-velocity of ice', 'm yr-1', conversion=365.0*86400.0*US%L_T_to_m_s) CS%id_v_shelf = register_diag_field('ice_shelf_model','v_shelf',CS%diag%axesCv1, Time, & 'y-velocity of ice', 'm yr-1', conversion=365.0*86400.0*US%L_T_to_m_s) CS%id_taudx_shelf = register_diag_field('ice_shelf_model','taudx_shelf',CS%diag%axesT1, Time, & 'x-driving stress of ice', 'kPa', conversion=1.e-9*US%L_T_to_m_s) CS%id_taudy_shelf = register_diag_field('ice_shelf_model','taudy_shelf',CS%diag%axesT1, Time, & 'x-driving stress of ice', 'kPa', conversion=1.e-9*US%L_T_to_m_s) -! CS%id_u_mask = register_diag_field('ocean_model','u_mask',CS%diag%axesCu1, Time, & -! 'mask for u-nodes', 'none') CS%id_u_mask = register_diag_field('ice_shelf_model','u_mask',CS%diag%axesCu1, Time, & 'mask for u-nodes', 'none') -! CS%id_v_mask = register_diag_field('ocean_model','v_mask',CS%diag%axesCv1, Time, & -! 'mask for v-nodes', 'none') CS%id_v_mask = register_diag_field('ice_shelf_model','v_mask',CS%diag%axesCv1, Time, & 'mask for v-nodes', 'none') -! CS%id_surf_elev = register_diag_field('ocean_model','ice_surf',CS%diag%axesT1, Time, & -! 'ice surf elev', 'm') -! CS%id_ground_frac = register_diag_field('ocean_model','ice_ground_frac',CS%diag%axesT1, Time, & -! 'fraction of cell that is grounded', 'none') CS%id_ground_frac = register_diag_field('ice_shelf_model','ice_ground_frac',CS%diag%axesT1, Time, & 'fraction of cell that is grounded', 'none') -! CS%id_col_thick = register_diag_field('ocean_model','col_thick',CS%diag%axesT1, Time, & -! 'ocean column thickness passed to ice model', 'm', conversion=US%Z_to_m) CS%id_col_thick = register_diag_field('ice_shelf_model','col_thick',CS%diag%axesT1, Time, & 'ocean column thickness passed to ice model', 'm', conversion=US%Z_to_m) CS%id_visc_shelf = register_diag_field('ice_shelf_model','ice_visc',CS%diag%axesT1, Time, & 'viscosity', 'm', conversion=1e-6*US%Z_to_m) -! CS%id_OD_av = register_diag_field('ocean_model','OD_av',CS%diag%axesT1, Time, & -! 'intermediate ocean column thickness passed to ice model', 'm', conversion=US%Z_to_m) CS%id_OD_av = register_diag_field('ice_shelf_model','OD_av',CS%diag%axesT1, Time, & 'intermediate ocean column thickness passed to ice model', 'm', conversion=US%Z_to_m) CS%id_h_after_uflux = register_diag_field('ice_shelf_model','h_after_uflux',CS%diag%axesT1, Time, & @@ -606,8 +571,8 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ if (new_sim) then call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: initialize ice velocity.") call update_OD_ffrac_uncoupled(CS, G, ISS%h_shelf(:,:)) -! call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf, iters, Time) - call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) !OVS 02/08/21 + !call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf, iters, Time) + call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) if (CS%id_u_shelf > 0) call post_data(CS%id_u_shelf, CS%u_shelf, CS%diag) if (CS%id_v_shelf > 0) call post_data(CS%id_v_shelf, CS%v_shelf,CS%diag) if (CS%id_taudx_shelf > 0) call post_data(CS%id_taudx_shelf, CS%taudx_shelf,CS%diag) @@ -655,8 +620,7 @@ subroutine initialize_diagnostic_fields(CS, ISS, G, US, Time) enddo enddo -! call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf, iters, dummy_time) - call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) !OVS 02/08/21 + call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) end subroutine initialize_diagnostic_fields !> This function returns the global maximum advective timestep that can be taken based on the current @@ -716,7 +680,7 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled coupled_GL = .false. if (present(ocean_mass) .and. present(coupled_grounding)) coupled_GL = coupled_grounding ! - call ice_shelf_advect(CS, ISS, G, time_step, Time) !OVS 02/08/21 + call ice_shelf_advect(CS, ISS, G, time_step, Time) CS%elapsed_velocity_time = CS%elapsed_velocity_time + time_step if (CS%elapsed_velocity_time >= CS%velocity_update_time_step) update_ice_vel = .true. @@ -729,7 +693,7 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled if (update_ice_vel) then ! call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf, iters, Time) - call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) !OVS 02/08/21 + call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) endif call ice_shelf_temp(CS, ISS, G, US, time_step, ISS%water_flux, Time) @@ -741,10 +705,10 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled if (CS%id_v_shelf > 0) call post_data(CS%id_v_shelf, CS%v_shelf, CS%diag) if (CS%id_t_shelf > 0) call post_data(CS%id_t_shelf,CS%t_shelf,CS%diag) if (CS%id_taudx_shelf > 0) call post_data(CS%id_taudx_shelf, CS%taudx_shelf, CS%diag) - if (CS%id_taudy_shelf > 0) call post_data(CS%id_taudy_shelf, CS%taudy_shelf, CS%diag) + if (CS%id_taudy_shelf > 0) call post_data(CS%id_taudy_shelf, CS%taudy_shelf, CS%diag) if (CS%id_ground_frac > 0) call post_data(CS%id_ground_frac, CS%ground_frac,CS%diag) if (CS%id_OD_av >0) call post_data(CS%id_OD_av, CS%OD_av,CS%diag) - if (CS%id_visc_shelf > 0) call post_data(CS%id_visc_shelf, CS%ice_visc,CS%diag) + if (CS%id_visc_shelf > 0) call post_data(CS%id_visc_shelf, CS%ice_visc,CS%diag) if (CS%id_u_mask > 0) call post_data(CS%id_u_mask,CS%umask,CS%diag) if (CS%id_v_mask > 0) call post_data(CS%id_v_mask,CS%vmask,CS%diag) @@ -908,7 +872,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, ite call calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, CS%OD_av) call pass_vector(taudx, taudy, G%domain, TO_ALL, BGRID_NE) !OVS 02/01/21 ! call pass_var(taudx, G%Domain) !OVS 01/21/21 -! call pass_var(taudy, G%Domain) !OVS 01/21/21 +! call pass_var(taudy, G%Domain) !OVS 01/21/21 ! This is to determine which cells contain the grounding line, the criterion being that the cell ! is ice-covered, with some nodes floating and some grounded flotation condition is estimated by ! assuming topography is cellwise constant and H is bilinear in a cell; floating where @@ -1072,7 +1036,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, ite if (err_max <= CS%nonlinear_tolerance * err_init) then write(mesg,*) "ice_shelf_solve_outer: nonlinear fractional residual = ", err_max/err_init - call MOM_mesg(mesg) + call MOM_mesg(mesg) write(mesg,*) "ice_shelf_solve_outer: exiting nonlinear solve after ",iter," iterations" ! call MOM_mesg(mesg, 5) call MOM_mesg(mesg) @@ -1354,8 +1318,8 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H ! pass vectors call pass_vector(Du, Dv, G%domain, TO_ALL, BGRID_NE) call pass_vector(u_shlf, v_shlf, G%domain, TO_ALL, BGRID_NE) - call pass_var(u_shlf, G%domain) - call pass_var(v_shlf, G%domain) + call pass_var(u_shlf, G%domain) + call pass_var(v_shlf, G%domain) call pass_vector(Ru, Rv, G%domain, TO_ALL, BGRID_NE) cg_halo = 3 endif @@ -1856,7 +1820,7 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) ! check whether the ice is floating or grounded ! do j=jsc-1,jec+1 !OVS 02/02/21 ! do i=isc-1,iec+1 !OVS 02/02/21 - do j=jsc-G%domain%njhalo,jec+G%domain%njhalo !OVS 02/02/21 + do j=jsc-G%domain%njhalo,jec+G%domain%njhalo !OVS 02/02/21 do i=isc-G%domain%nihalo,iec+G%domain%nihalo !OVS 02/02/21 ! if (ISS%h_shelf(i,j) < rhow/rho * G%bathyT(i,j)) then @@ -1866,7 +1830,7 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) enddo - enddo + enddo do j=jsc-1,jec+1 do i=isc-1,iec+1 cnt = 0 @@ -1956,16 +1920,16 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) taudy(I,J-1) = taudy(I,J-1) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j) endif ! NW vertex - if (ISS%hmask(I-1,J) == 1) then + if (ISS%hmask(I-1,J) == 1) then taudx(I-1,J) = taudx(I-1,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sx * G%areaT(i,j) taudy(I-1,J) = taudy(I-1,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j) - endif + endif ! NE vertex if (ISS%hmask(I,J) == 1) then taudx(I,J) = taudx(I,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sx * G%areaT(i,j) taudy(I,J) = taudy(I,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j) endif - if (CS%ground_frac(i,j) == 1) then + if (CS%ground_frac(i,j) == 1) then ! neumann_val = .5 * grav * (rho * ISS%h_shelf(i,j)**2 - rhow * G%bathyT(i,j)**2) neumann_val = .5 * grav * (1-rho/rhow) * rho * ISS%h_shelf(i,j)**2 else @@ -2325,7 +2289,7 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, ! Phi(2*i-1,j) gives d(Phi_i)/dx at quadrature point j ! Phi(2*i,j) gives d(Phi_i)/dy at quadrature point j - do iq=1,2 ; do jq=1,2 ; do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; Jtgt = J-2+jphi !Jtgt = J-2-jphi !OVS fix index + do iq=1,2 ; do jq=1,2 ; do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; Jtgt = J-2+jphi ilq = 1 ; if (iq == iphi) ilq = 2 jlq = 1 ; if (jq == jphi) jlq = 2 @@ -2584,7 +2548,7 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) integer :: i, j, iscq, iecq, jscq, jecq, isd, jsd, ied, jed, iegq, jegq, iq, jq integer :: giec, gjec, gisc, gjsc, cnt, isc, jsc, iec, jec, is, js, i_off, j_off real :: Visc_coef, n_g - real :: ux, uy, vx, vy + real :: ux, uy, vx, vy real :: eps_min, dxh, dyh ! Velocity shears [T-1 ~> s-1] real, dimension(8,4) :: Phi real, dimension(2) :: xquad @@ -2602,7 +2566,7 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) n_g = CS%n_glen; eps_min = CS%eps_glen_min ! CS%ice_visc(:,:) = 0.0 -! ux(:,:) = 0.0; uy(:,:) = 0.0; vx(:,:) =0.0; vy(:,:) =0.0 +! ux(:,:) = 0.0; uy(:,:) = 0.0; vx(:,:) =0.0; vy(:,:) =0.0 ! eII(:,:) = (US%s_to_T**2 * (eps_min**2)) Visc_coef = US%kg_m2s_to_RZ_T*US%m_to_L*US%Z_to_L*(CS%A_glen_isothermal)**(-1./CS%n_glen) !OVS '-' in the exponent ! call pass_vector(u_shlf, v_shlf, G%domain, TO_ALL, BGRID_NE) @@ -2616,11 +2580,11 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) ! vx(i,j) = ((v_shlf(I,J) + v_shlf(I,J-1)) - (v_shlf(I-1,J) + v_shlf(I-1,J-1))) / (2*G%dxT(i,j)) ! uy(i,j) = ((u_shlf(I,J) + u_shlf(I-1,J)) - (u_shlf(I,J-1) + u_shlf(I-1,J-1))) / (2*G%dyT(i,j)) ! vy(i,j) = ((v_shlf(I,J) + v_shlf(I-1,J)) - (v_shlf(I,J-1) + v_shlf(I-1,J-1))) / (2*G%dyT(i,j)) -! endif +! endif +! enddo ! enddo -! enddo ! call pass_vector(ux, uy, G%domain, TO_ALL, BGRID_NE) -! call pass_vector(vx, vy, G%domain, TO_ALL, BGRID_NE) +! call pass_vector(vx, vy, G%domain, TO_ALL, BGRID_NE) ux = ((u_shlf(I,J) + u_shlf(I,J-1) + u_shlf(I,J+1)) - & (u_shlf(I-1,J) + u_shlf(I-1,J-1) + u_shlf(I-1,J+1))) / (3*G%dxT(i,j)) vx = ((v_shlf(I,J) + v_shlf(I,J-1) + v_shlf(I,J+1)) - & @@ -2628,14 +2592,14 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) uy = ((u_shlf(I,J) + u_shlf(I-1,J) + u_shlf(I+1,J)) - & (u_shlf(I,J-1) + u_shlf(I-1,J-1) + u_shlf(I+1,J-1))) / (3*G%dyT(i,j)) vy = ((v_shlf(I,J) + v_shlf(I-1,J)+ v_shlf(I+1,J)) - & - (v_shlf(I,J-1) + v_shlf(I-1,J-1)+ v_shlf(I+1,J-1))) / (3*G%dyT(i,j)) + (v_shlf(I,J-1) + v_shlf(I-1,J-1)+ v_shlf(I+1,J-1))) / (3*G%dyT(i,j)) ! ux = ((u_shlf(I,J) + u_shlf(I,J-1)) - (u_shlf(I-1,J) + u_shlf(I-1,J-1))) / (2*G%dxT(i,j)) ! vx = ((v_shlf(I,J) + v_shlf(I,J-1)) - (v_shlf(I-1,J) + v_shlf(I-1,J-1))) / (2*G%dxT(i,j)) ! uy = ((u_shlf(I,J) + u_shlf(I-1,J)) - (u_shlf(I,J-1) + u_shlf(I-1,J-1))) / (2*G%dyT(i,j)) ! vy = ((v_shlf(I,J) + v_shlf(I-1,J)) - (v_shlf(I,J-1) + v_shlf(I-1,J-1))) / (2*G%dyT(i,j)) CS%ice_visc(i,j) = 0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * & (US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2))**((1.-n_g)/(2.*n_g)) -! CS%ice_visc(i,j) =1e15*(G%areaT(i,j) * ISS%h_shelf(i,j)) ! constant viscocity for debugging + CS%ice_visc(i,j) =1e14*(G%areaT(i,j) * ISS%h_shelf(i,j)) ! constant viscocity for debugging ! umid = ((u_shlf(I,J) + u_shlf(I-1,J-1)) + (u_shlf(I,J-1) + u_shlf(I-1,J))) * 0.25 ! vmid = ((v_shlf(I,J) + v_shlf(I-1,J-1)) + (v_shlf(I,J-1) + v_shlf(I-1,J))) * 0.25 ! unorm = sqrt(umid**2 + vmid**2 + eps_min**2*(G%dxT(i,j)**2 + G%dyT(i,j)**2)) @@ -2644,183 +2608,6 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) enddo enddo -! do j=jsc-1,jec+1 -! do i=isc-1,iec+1 -!! do j=jsd+1,jed!-1 OVS 02/01/21 -!! do i=isd+1,ied!-1 OVS 02/01/21 - -! if (ISS%hmask(i,j) == 1) then -! CS%ice_visc(i,j) = 0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * & -! (US%s_to_T**2 * (ux(i,j)**2 + vy(i,j)**2 + ux(i,j)*vy(i,j) + 0.25*(uy(i,j)+vx(i,j))**2 + eps_min**2))**((1.-n_g)/(2.*n_g)) -! endif -! enddo -! enddo - -! xquad(1) = .5 * (1-sqrt(1./3)) ; xquad(2) = .5 * (1+sqrt(1./3)) -! do j=jsc-1,jec+1 -! do i=isc-1,iec+1 -! cnt = 0 -! ux = 0 -! uy = 0 -! vx = 0 -! vy = 0 -! dxh = G%dxT(i,j) -! dyh = G%dyT(i,j) - -! if (ISS%hmask(i,j) == 1) then ! we are inside the global computational bdry, at an ice-filled cell - -! call bilinear_shape_fn_grid(G, i, j, Phi) -! do jq = 1,2 -! do iq = 1,2 - -! ux = u_shlf(I-1,J-1) * Phi(1,2*(jq-1)+iq) + & -! u_shlf(I,J-1) * Phi(3,2*(jq-1)+iq) + & -! u_shlf(I-1,J) * Phi(5,2*(jq-1)+iq) + & -! u_shlf(I,J) * Phi(7,2*(jq-1)+iq) - -! vx = v_shlf(I-1,J-1) * Phi(1,2*(jq-1)+iq) + & -! v_shlf(I,J-1) * Phi(3,2*(jq-1)+iq) + & -! v_shlf(I-1,J) * Phi(5,2*(jq-1)+iq) + & -! v_shlf(I,J) * Phi(7,2*(jq-1)+iq) - -! uy = u_shlf(I-1,J-1) * Phi(2,2*(jq-1)+iq) + & -! u_shlf(I,J-1) * Phi(4,2*(jq-1)+iq) + & -! u_shlf(I-1,J) * Phi(6,2*(jq-1)+iq) + & -! u_shlf(I,J) * Phi(8,2*(jq-1)+iq) - -! vy = v_shlf(I-1,J-1) * Phi(2,2*(jq-1)+iq) + & -! v_shlf(I,J-1) * Phi(4,2*(jq-1)+iq) + & -! v_shlf(I-1,J) * Phi(6,2*(jq-1)+iq) + & -! v_shlf(I,J) * Phi(8,2*(jq-1)+iq) -! enddo -! enddo - ! calculate sx -! if ((i+i_off) == gisc) then ! at left computational bdry -! if (ISS%hmask(i+1,j) == 1) then -! ux = (u_shlf(i+1,j)-u_shlf(i,j))/dxh -! vx = (v_shlf(i+1,j)-v_shlf(i,j))/dxh -! else -! ux = 0 -! vx = 0 -! endif -! elseif ((i+i_off) == giec) then ! at east computational bdry -! if (ISS%hmask(i-1,j) == 1) then -! ux = (u_shlf(i,j)-u_shlf(i-1,j))/dxh -! vx = (v_shlf(i,j)-v_shlf(i-1,j))/dxh -! else -! ux = 0 -! vx = 0 -! endif -! else ! interior -! if (ISS%hmask(i+1,j) == 1) then -! cnt = cnt+1 -! ux = u_shlf(i+1,j) -! vx = v_shlf(i+1,j) -! else -! ux = u_shlf(i,j) -! vx = v_shlf(i,j) -! endif -! if (ISS%hmask(i-1,j) == 1) then -! cnt = cnt+1 -! ux = ux - u_shlf(i-1,j) -! vx = vx - v_shlf(i-1,j) -! else -! ux = ux - u_shlf(i,j) -! vx = vx - v_shlf(i,j) -! endif -! if (cnt == 0) then -! ux = 0 -! vx = 0 -! else -! ux = ux / (cnt * dxh) -! vx = vx / (cnt * dxh) -! endif -!! endif -! cnt = 0 - - ! calculate sy, similarly -! if ((j+j_off) == gjsc) then ! at south computational bdry -! if (ISS%hmask(i,j+1) == 1) then -! uy = (u_shlf(i,j+1)-u_shlf(i,j))/dyh -! vy = (v_shlf(i,j+1)-v_shlf(i,j))/dyh -! else -! vy = 0 -! endif -! elseif ((j+j_off) == gjec) then ! at nprth computational bdry -! if (ISS%hmask(i,j-1) == 1) then -! uy = (u_shlf(i,j)-u_shlf(i,j-1))/dyh -! vy = (v_shlf(i,j)-v_shlf(i,j-1))/dyh -! else -! uy = 0 -! vy = 0 -! endif -! else ! interior -! if (ISS%hmask(i,j+1) == 1) then -! cnt = cnt+1 -! uy = u_shlf(i,j+1) -! vy = v_shlf(i,j+1) -! else -! uy = u_shlf(i,j) -! vy = v_shlf(i,j) -! endif -! if (ISS%hmask(i,j-1) == 1) then -! cnt = cnt+1 -! uy = uy - u_shlf(i,j-1) -! vy = vy - v_shlf(i,j-1) -! else -! uy = uy - u_shlf(i,j) -! vy = vy - v_shlf(i,j) -! endif -! if (cnt == 0) then -! uy = 0 -! vy = 0 -! else -! uy = uy / (cnt * dyh) -! vy = vy / (cnt * dyh) -! endif -!! endif - -! ! SW vertex -! if (ISS%hmask(I-1,J-1) == 1) then -! eII(i-1,j-1) = eII(i-1,j-1)+.125*(US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2)) -! endif - ! SE vertex -! if (ISS%hmask(I,J-1) == 1) then -! eII(i,j-1) = eII(i,j-1)+.125*(US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2)) - -! CS%ice_visc(i,j-1) = CS%ice_visc(i,j-1)+.25*0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * & -! (US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2))**((1.-n_g)/(2.*n_g)) -! endif - ! NW vertex -! if (ISS%hmask(I-1,J) == 1) then -! eII(i-1,j) = eII(i-1,j)+.125*(US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2)) - -! CS%ice_visc(i-1,j) = CS%ice_visc(i-1,j)+.25*0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * & -! (US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2))**((1.-n_g)/(2.*n_g)) -! endif - ! NE vertex -! if (ISS%hmask(I,J) == 1) then -! eII(i,j) = eII(i,j)+.25*(US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2)) -! eII(i,j) = (US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2)) - -! CS%ice_visc(i,j) = CS%ice_visc(i,j)+.25*0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * & -! (US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2))**((1.-n_g)/(2.*n_g)) -! endif -! if (ISS%hmask(I+1,J+1) == 1) then -! eII(i+1,j+1) = eII(i+1,j+1)+.125*(US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2)) -! endif -! if (ISS%hmask(I,J+1) == 1) then -! eII(i,j+1) = eII(i,j+1)+.125*(US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2)) -! endif -! if (ISS%hmask(I+1,J) == 1) then -! eII(i+1,j) = eII(i+1,j)+.125*(US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2)) -! endif -! CS%ice_visc(i,j) =0.5 * Visc_coef*(G%areaT(i,j) * ISS%h_shelf(i,j))*eII(i,j)**((1.-n_g)/(2.*n_g)) -! endif -! CS%ice_visc(i,j) =0.5 * Visc_coef*(G%areaT(i,j) * ISS%h_shelf(i,j))*eII(i,j)**((1.-n_g)/(2.*n_g)) - ! CS%ice_visc(i,j) =1e15*(G%areaT(i,j) * ISS%h_shelf(i,j)) !constant viscosity for debugging -! enddo -! enddo end subroutine calc_shelf_visc subroutine calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) @@ -3181,7 +2968,7 @@ subroutine update_velocity_masks(CS, G, hmask, umask, vmask, u_face_mask, v_face umask(I-1+k,J)=3. !vmask(I-1+k,J)=0. vmask(I-1+k,J)=3. - !u_face_mask(I-1+k,j-1)=3. + !u_face_mask(I-1+k,j-1)=3. ! umask(I-1+k,J-1:J)=3. ! vmask(I-1+k,J-1:J)=0. ! u_face_mask(I-1+k,j)=3. diff --git a/src/ice_shelf/MOM_ice_shelf_initialize.F90 b/src/ice_shelf/MOM_ice_shelf_initialize.F90 index f9f31a373e..d77efa358b 100644 --- a/src/ice_shelf/MOM_ice_shelf_initialize.F90 +++ b/src/ice_shelf/MOM_ice_shelf_initialize.F90 @@ -269,7 +269,7 @@ end subroutine initialize_ice_thickness_channel !> Initialize ice shelf boundary conditions for a channel configuration subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_bdry, & u_flux_bdry_val, v_flux_bdry_val, u_bdry_val, v_bdry_val, u_shelf, v_shelf, h_bdry_val, & - thickness_bdry_val, hmask, h_shelf, G,& + thickness_bdry_val, hmask, h_shelf, G,& US, PF ) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure @@ -289,9 +289,9 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b real, dimension(SZIB_(G),SZJB_(G)), & intent(inout) :: v_bdry_val !< The meridional ice shelf velocity at open real, dimension(SZIB_(G),SZJB_(G)), & - intent(inout) :: u_shelf !< The zonal ice shelf velocity [L T-1 ~> m s-1]. + intent(inout) :: u_shelf !< The zonal ice shelf velocity [L T-1 ~> m s-1]. real, dimension(SZIB_(G),SZJB_(G)), & - intent(inout) :: v_shelf !< The meridional ice shelf velocity [L T-1 ~> m s-1]. + intent(inout) :: v_shelf !< The meridional ice shelf velocity [L T-1 ~> m s-1]. real, dimension(SZDI_(G),SZDJ_(G)), & intent(inout) :: thickness_bdry_val !< The ice shelf thickness at open boundaries [Z ~> m] !! boundary vertices [L T-1 ~> m s-1]. @@ -301,7 +301,7 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b intent(inout) :: hmask !< A mask indicating which tracer points are !! partly or fully covered by an ice-shelf real, dimension(SZDI_(G),SZDJ_(G)), & - intent(inout) :: h_shelf !< Ice-shelf thickness OVS 11/10/20 + intent(inout) :: h_shelf !< Ice-shelf thickness OVS 11/10/20 type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors type(param_file_type), intent(in) :: PF !< A structure to parse for run-time parameters @@ -367,26 +367,26 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b if (G%geoLatBu(i,j-1) == southlat) then !bot boundary if (len_stress == 0. .OR. G%geoLonCv(i,j) <= len_stress) then v_face_mask_bdry(i,j+1) = 0. - ! u_face_mask_bdry(i,j-1) = 3. - u_face_mask_bdry(i,j) = 3. + ! u_face_mask_bdry(i,j-1) = 3. + u_face_mask_bdry(i,j) = 3. u_bdry_val(i,j) = 0. - v_bdry_val(i,j) = 0. + v_bdry_val(i,j) = 0. else v_face_mask_bdry(i,j+1) = 1. - u_face_mask_bdry(i,j) = 3. + u_face_mask_bdry(i,j) = 3. u_bdry_val(i,j) = 0. - v_bdry_val(i,j) = 0. + v_bdry_val(i,j) = 0. endif elseif (G%geoLatBu(i,j-1) == southlat+lenlat) then !top boundary if (len_stress == 0. .OR. G%geoLonCv(i,j) <= len_stress) then v_face_mask_bdry(i,j-1) = 0. - u_face_mask_bdry(i,j-1) = 3. + u_face_mask_bdry(i,j-1) = 3. else !v_face_mask_bdry(i,j-1) = 1. - v_face_mask_bdry(i,j-1) = 3. - u_face_mask_bdry(i,j-1) = 3. - !u_bdry_val(i,j) = 0. - !hmask(i,j) = 0.0 + v_face_mask_bdry(i,j-1) = 3. + u_face_mask_bdry(i,j-1) = 3. + !u_bdry_val(i,j) = 0. + !hmask(i,j) = 0.0 endif endif @@ -396,9 +396,9 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b endif enddo - enddo + enddo ! if (.not. G%symmetric) then -!! do j=G%jsd,G%jed +!! do j=G%jsd,G%jed !! do i=G%isd,G%ied ! do j=jsc-1,jec+1 ! do i=isc-1,iec+1 @@ -416,7 +416,7 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b ! v_shelf(I-1,J-1) = v_bdry_val(I-1,J-1) ! v_shelf(I,J-1) = v_bdry_val(I,J-1) ! endif -! enddo +! enddo ! enddo ! endif end subroutine initialize_ice_shelf_boundary_channel @@ -583,7 +583,8 @@ end subroutine initialize_ice_shelf_boundary_channel !END MJH end subroutine initialize_ice_shelf_boundary_channel !> Initialize ice shelf flow from file -subroutine initialize_ice_flow_from_file(u_shelf, v_shelf,ice_visc,float_cond, hmask,h_shelf, G, US, PF) +subroutine initialize_ice_flow_from_file(u_shelf, v_shelf,ice_visc,float_cond,& + hmask,h_shelf, G, US, PF) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure real, dimension(SZDI_(G),SZDJ_(G)), & intent(inout) :: u_shelf !< The ice shelf u velocity [Z ~> m T ~>s]. @@ -593,13 +594,13 @@ subroutine initialize_ice_flow_from_file(u_shelf, v_shelf,ice_visc,float_cond, h intent(inout) :: ice_visc !< The ice shelf viscosity [Pa ~> m T ~> s]. real, dimension(SZDI_(G),SZDJ_(G)), & intent(inout) :: float_cond !< An array indicating where the ice - !! shelf is floating: 0 if floating, 1 if not. + !! shelf is floating: 0 if floating, 1 if not. real, dimension(SZDI_(G),SZDJ_(G)), & intent(in) :: hmask !< A mask indicating which tracer points are !! partly or fully covered by an ice-shelf real, dimension(SZDI_(G),SZDJ_(G)), & intent(in) :: h_shelf !< A mask indicating which tracer points are - !! partly or fully covered by an ice-shelf + !! partly or fully covered by an ice-shelf type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors type(param_file_type), intent(in) :: PF !< A structure to parse for run-time parameters From aed4f0ee71e4c6e179ebe196e6b4d48935c76ee1 Mon Sep 17 00:00:00 2001 From: Olga Sergienko Date: Wed, 3 Mar 2021 16:35:41 -0500 Subject: [PATCH 14/17] More style errors --- src/ice_shelf/MOM_ice_shelf_dynamics.F90 | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index 8640902989..4abcda0aa0 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -160,8 +160,9 @@ module MOM_ice_shelf_dynamics id_u_mask = -1, id_v_mask = -1, id_t_mask = -1 !>@} ! ids for outputting intermediate thickness in advection subroutine (debugging) + !>@{ Diagnostic handles for debugging integer :: id_h_after_uflux = -1, id_h_after_vflux = -1, id_h_after_adv = -1, id_visc_shelf = -1 - + !>@} type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to control diagnostic output. end type ice_shelf_dyn_CS @@ -809,8 +810,9 @@ subroutine ice_shelf_advect(CS, ISS, G, time_step, Time) end subroutine ice_shelf_advect +!>This subroutine computes u- and v-velocities of the ice shelf iterating on non-linear ice viscosity !subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, iters, time) - subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, iters, Time) !OVS 02/08/21 + subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, iters, Time) type(ice_shelf_dyn_CS), intent(inout) :: CS !< The ice shelf dynamics control structure type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe !! the ice-shelf state @@ -823,7 +825,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, ite integer, intent(out) :: iters !< The number of iterations used in the solver. type(time_type), intent(in) :: Time !< The current model time - real, dimension(SZDIB_(G),SZDJB_(G)) :: taudx, taudy ! Driving stresses at q-points [R L3 Z T-2 ~> kg m s-2] + real, dimension(SZDIB_(G),SZDJB_(G)) :: taudx, taudy !< Driving stresses at q-points [R L3 Z T-2 ~> kg m s-2] real, dimension(SZDIB_(G),SZDJB_(G)) :: u_bdry_cont ! Boundary u-stress contribution [R L3 Z T-2 ~> kg m s-2] real, dimension(SZDIB_(G),SZDJB_(G)) :: v_bdry_cont ! Boundary v-stress contribution [R L3 Z T-2 ~> kg m s-2] real, dimension(SZDIB_(G),SZDJB_(G)) :: Au, Av ! The retarding lateral stress contributions [R L3 Z T-2 ~> kg m s-2] From 43dadc16a357c891df640400c4cb7902e3064e5d Mon Sep 17 00:00:00 2001 From: Olga Sergienko Date: Wed, 3 Mar 2021 17:11:13 -0500 Subject: [PATCH 15/17] Defined variables in ice_shelf_solve_outer --- src/ice_shelf/MOM_ice_shelf_dynamics.F90 | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index 4abcda0aa0..299eda4f33 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -825,7 +825,10 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, ite integer, intent(out) :: iters !< The number of iterations used in the solver. type(time_type), intent(in) :: Time !< The current model time - real, dimension(SZDIB_(G),SZDJB_(G)) :: taudx, taudy !< Driving stresses at q-points [R L3 Z T-2 ~> kg m s-2] + real, dimension(SZDIB_(G),SZDJB_(G)), & + intent(out) :: taudx !< Driving x-stress at q-points [R L3 Z T-2 ~> kg m s-2] + real, dimension(SZDIB_(G),SZDJB_(G)), & + intent(out) :: taudy !< Driving y-stress at q-points [R L3 Z T-2 ~> kg m s-2] real, dimension(SZDIB_(G),SZDJB_(G)) :: u_bdry_cont ! Boundary u-stress contribution [R L3 Z T-2 ~> kg m s-2] real, dimension(SZDIB_(G),SZDJB_(G)) :: v_bdry_cont ! Boundary v-stress contribution [R L3 Z T-2 ~> kg m s-2] real, dimension(SZDIB_(G),SZDJB_(G)) :: Au, Av ! The retarding lateral stress contributions [R L3 Z T-2 ~> kg m s-2] @@ -2601,7 +2604,7 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) ! vy = ((v_shlf(I,J) + v_shlf(I-1,J)) - (v_shlf(I,J-1) + v_shlf(I-1,J-1))) / (2*G%dyT(i,j)) CS%ice_visc(i,j) = 0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * & (US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2))**((1.-n_g)/(2.*n_g)) - CS%ice_visc(i,j) =1e14*(G%areaT(i,j) * ISS%h_shelf(i,j)) ! constant viscocity for debugging +! CS%ice_visc(i,j) =1e14*(G%areaT(i,j) * ISS%h_shelf(i,j)) ! constant viscocity for debugging ! umid = ((u_shlf(I,J) + u_shlf(I-1,J-1)) + (u_shlf(I,J-1) + u_shlf(I-1,J))) * 0.25 ! vmid = ((v_shlf(I,J) + v_shlf(I-1,J-1)) + (v_shlf(I,J-1) + v_shlf(I-1,J))) * 0.25 ! unorm = sqrt(umid**2 + vmid**2 + eps_min**2*(G%dxT(i,j)**2 + G%dyT(i,j)**2)) From b47e493d3ba28ac5544c22231a8cfb4e88feeed8 Mon Sep 17 00:00:00 2001 From: Olga Sergienko Date: Tue, 9 Mar 2021 13:15:10 -0500 Subject: [PATCH 16/17] Removed blocks of commented code. Added parentheses in calc_shelf_visc --- src/ice_shelf/MOM_ice_shelf_dynamics.F90 | 56 ++----- src/ice_shelf/MOM_ice_shelf_initialize.F90 | 184 +-------------------- 2 files changed, 18 insertions(+), 222 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index 299eda4f33..6b30b2d83d 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -579,12 +579,11 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ if (CS%id_taudx_shelf > 0) call post_data(CS%id_taudx_shelf, CS%taudx_shelf,CS%diag) if (CS%id_taudy_shelf > 0) call post_data(CS%id_taudy_shelf, CS%taudy_shelf,CS%diag) if (CS%id_visc_shelf > 0) call post_data(CS%id_visc_shelf, CS%ice_visc,CS%diag) - endif -!!! OVS vertically integrated temperature ! CS%id_t_shelf = register_diag_field('ocean_model','t_shelf',CS%diag%axesT1, Time, & ! 'T of ice', 'oC') ! CS%id_t_mask = register_diag_field('ocean_model','tmask',CS%diag%axesT1, Time, & ! 'mask for T-nodes', 'none') + endif endif end subroutine initialize_ice_shelf_dyn @@ -875,9 +874,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, ite enddo call calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, CS%OD_av) - call pass_vector(taudx, taudy, G%domain, TO_ALL, BGRID_NE) !OVS 02/01/21 -! call pass_var(taudx, G%Domain) !OVS 01/21/21 -! call pass_var(taudy, G%Domain) !OVS 01/21/21 + call pass_vector(taudx, taudy, G%domain, TO_ALL, BGRID_NE) ! This is to determine which cells contain the grounding line, the criterion being that the cell ! is ice-covered, with some nodes floating and some grounded flotation condition is estimated by ! assuming topography is cellwise constant and H is bilinear in a cell; floating where @@ -917,12 +914,10 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, ite call bilinear_shape_fn_grid(G, i, j, Phi(:,:,i,j)) enddo ; enddo - call calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) !OVS 02/24/21 + call calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) call pass_var(CS%ice_visc, G%domain) -! call pass_vector(CS%ice_visc, G%domain, TO_ALL, BGRID_NE) !OVS 02/11/21 call calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) call pass_var(CS%basal_traction, G%domain) -! call pass_vector(CS%ice_visc,CS%basal_traction, G%domain, TO_ALL, BGRID_NE) !OVS 02/11/21 ! This makes sure basal stress is only applied when it is supposed to be do j=G%jsd,G%jed ; do i=G%isd,G%ied @@ -937,7 +932,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, ite call CG_action(Au, Av, u_shlf, v_shlf, Phi, Phisub, CS%umask, CS%vmask, ISS%hmask, H_node, & CS%ice_visc, float_cond, G%bathyT, CS%basal_traction, & G, US, G%isc-1, G%iec+1, G%jsc-1, G%jec+1, rhoi_rhow) - call pass_vector(Au,Av,G%domain) !OVS pass Au and Av + call pass_vector(Au,Av,G%domain) if (CS%nonlin_solve_err_mode == 1) then err_init = 0 ; err_tempu = 0 ; err_tempv = 0 do J=G%IscB,G%JecB ; do I=G%IscB,G%IecB @@ -971,7 +966,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, ite write(mesg,*) "ice_shelf_solve_outer: linear solve done in ",iters," iterations" call MOM_mesg(mesg, 5) - call calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) !OVS 02/24/21 + call calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) call pass_var(CS%ice_visc, G%domain) call calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) call pass_var(CS%basal_traction, G%domain) @@ -1823,10 +1818,8 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) S(:,:) = BASE(:,:) + ISS%h_shelf(:,:) ! check whether the ice is floating or grounded -! do j=jsc-1,jec+1 !OVS 02/02/21 -! do i=isc-1,iec+1 !OVS 02/02/21 - do j=jsc-G%domain%njhalo,jec+G%domain%njhalo !OVS 02/02/21 - do i=isc-G%domain%nihalo,iec+G%domain%nihalo !OVS 02/02/21 + do j=jsc-G%domain%njhalo,jec+G%domain%njhalo + do i=isc-G%domain%nihalo,iec+G%domain%nihalo ! if (ISS%h_shelf(i,j) < rhow/rho * G%bathyT(i,j)) then if (rhoi_rhow * ISS%h_shelf(i,j) - G%bathyT(i,j) <= 0) then @@ -2160,7 +2153,7 @@ subroutine CG_action(uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask, hmas v_shlf(I-1,J) * Phi(6,2*(jq-1)+iq,i,j) + & v_shlf(I,J) * Phi(8,2*(jq-1)+iq,i,j) - do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; ;Jtgt = J-2+jphi !Jtgt = J-2-jphi !OVS fix index + do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; ;Jtgt = J-2+jphi if (umask(Itgt,Jtgt) == 1) uret(Itgt,Jtgt) = uret(Itgt,Jtgt) + 0.25 * ice_visc(i,j) * & ((4*ux+2*vy) * Phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq,i,j) + & (uy+vx) * Phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq,i,j)) @@ -2338,7 +2331,7 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, if (float_cond(i,j) == 1) then Hcell(:,:) = H_node(i-1:i,j-1:j) call CG_diagonal_subgrid_basal(Phisub, Hcell, G%bathyT(i,j), dens_ratio, sub_ground) - do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; Jtgt = J-2+jphi !Jtgt = J-2-jphi !OVS fix index + do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; Jtgt = J-2+jphi if (CS%umask(Itgt,Jtgt) == 1) then u_diagonal(Itgt,Jtgt) = u_diagonal(Itgt,Jtgt) + sub_ground(iphi,jphi) * basal_trac(i,j) v_diagonal(Itgt,Jtgt) = v_diagonal(Itgt,Jtgt) + sub_ground(iphi,jphi) * basal_trac(i,j) @@ -2479,7 +2472,7 @@ subroutine apply_boundary_values(CS, ISS, G, US, time, Phisub, H_node, ice_visc, CS%v_bdry_val(I-1,J) * Phi(6,2*(jq-1)+iq) + & CS%v_bdry_val(I,J) * Phi(8,2*(jq-1)+iq) - do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; Jtgt = J-2+jphi !Jtgt = J-2-jphi !OVS fix index + do iphi=1,2 ; do jphi=1,2 ; Itgt = I-2+iphi ; Jtgt = J-2+jphi ilq = 1 ; if (iq == iphi) ilq = 2 jlq = 1 ; if (jq == jphi) jlq = 2 @@ -2590,14 +2583,14 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) ! enddo ! call pass_vector(ux, uy, G%domain, TO_ALL, BGRID_NE) ! call pass_vector(vx, vy, G%domain, TO_ALL, BGRID_NE) - ux = ((u_shlf(I,J) + u_shlf(I,J-1) + u_shlf(I,J+1)) - & - (u_shlf(I-1,J) + u_shlf(I-1,J-1) + u_shlf(I-1,J+1))) / (3*G%dxT(i,j)) + ux = ((u_shlf(I,J) + (u_shlf(I,J-1) + u_shlf(I,J+1))) - & + (u_shlf(I-1,J) + (u_shlf(I-1,J-1) + u_shlf(I-1,J+1)))) / (3*G%dxT(i,j)) vx = ((v_shlf(I,J) + v_shlf(I,J-1) + v_shlf(I,J+1)) - & - (v_shlf(I-1,J) + v_shlf(I-1,J-1) + v_shlf(I-1,J+1))) / (3*G%dxT(i,j)) - uy = ((u_shlf(I,J) + u_shlf(I-1,J) + u_shlf(I+1,J)) - & - (u_shlf(I,J-1) + u_shlf(I-1,J-1) + u_shlf(I+1,J-1))) / (3*G%dyT(i,j)) - vy = ((v_shlf(I,J) + v_shlf(I-1,J)+ v_shlf(I+1,J)) - & - (v_shlf(I,J-1) + v_shlf(I-1,J-1)+ v_shlf(I+1,J-1))) / (3*G%dyT(i,j)) + (v_shlf(I-1,J) + (v_shlf(I-1,J-1) + v_shlf(I-1,J+1)))) / (3*G%dxT(i,j)) + uy = ((u_shlf(I,J) + (u_shlf(I-1,J) + u_shlf(I+1,J))) - & + (u_shlf(I,J-1) + (u_shlf(I-1,J-1) + u_shlf(I+1,J-1)))) / (3*G%dyT(i,j)) + vy = ((v_shlf(I,J) + (v_shlf(I-1,J)+ v_shlf(I+1,J))) - & + (v_shlf(I,J-1) + (v_shlf(I-1,J-1)+ v_shlf(I+1,J-1)))) / (3*G%dyT(i,j)) ! ux = ((u_shlf(I,J) + u_shlf(I,J-1)) - (u_shlf(I-1,J) + u_shlf(I-1,J-1))) / (2*G%dxT(i,j)) ! vx = ((v_shlf(I,J) + v_shlf(I,J-1)) - (v_shlf(I-1,J) + v_shlf(I-1,J-1))) / (2*G%dxT(i,j)) ! uy = ((u_shlf(I,J) + u_shlf(I-1,J)) - (u_shlf(I,J-1) + u_shlf(I-1,J-1))) / (2*G%dyT(i,j)) @@ -3020,21 +3013,6 @@ subroutine update_velocity_masks(CS, G, hmask, umask, vmask, u_face_mask, v_face end select enddo - !if (CS%u_face_mask_bdry(I-1,j) >= 0) then ! Western boundary - ! u_face_mask(I-1,j) = CS%u_face_mask_bdry(I-1,j) - ! umask(I-1,J-1:J) = 3. - ! vmask(I-1,J-1:J) = 0. - !endif - - !if (j_off+j == gjsc+1) then ! SoutherN boundary - ! v_face_mask(i,J-1) = 0. - ! umask(I-1:I,J-1) = 0. - ! vmask(I-1:I,J-1) = 0. - !elseif (j_off+j == gjec) then ! Northern boundary - ! v_face_mask(i,J) = 0. - ! umask(I-1:I,J) = 0. - ! vmask(I-1:I,J) = 0. - !endif if (i < G%ied) then if ((hmask(i+1,j) == 0) .OR. (hmask(i+1,j) == 2)) then diff --git a/src/ice_shelf/MOM_ice_shelf_initialize.F90 b/src/ice_shelf/MOM_ice_shelf_initialize.F90 index d77efa358b..ff05ed7c6a 100644 --- a/src/ice_shelf/MOM_ice_shelf_initialize.F90 +++ b/src/ice_shelf/MOM_ice_shelf_initialize.F90 @@ -397,190 +397,8 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b enddo enddo -! if (.not. G%symmetric) then -!! do j=G%jsd,G%jed -!! do i=G%isd,G%ied -! do j=jsc-1,jec+1 -! do i=isc-1,iec+1 -!! if (((i+G%idg_offset) == (G%domain%nihalo+1)).and.(u_face_mask_bdry(I-1,j) == 3)) then -! if (u_face_mask_bdry(I-1,j) == 3) then -! u_shelf(I-1,J-1) = u_bdry_val(I-1,J-1) -! u_shelf(I-1,J) = u_bdry_val(I-1,J) -! v_shelf(I-1,J-1) = v_bdry_val(I-1,J-1) -! v_shelf(I-1,J) = v_bdry_val(I-1,J) -! endif -!! if (((j+G%jdg_offset) == (G%domain%njhalo+1)).and.(v_face_mask_bdry(i,J-1) == 3)) then -! if (v_face_mask_bdry(I,j-1) == 3) then -! u_shelf(I-1,J-1) = u_bdry_val(I-1,J-1) -! u_shelf(I,J-1) = u_bdry_val(I,J-1) -! v_shelf(I-1,J-1) = v_bdry_val(I-1,J-1) -! v_shelf(I,J-1) = v_bdry_val(I,J-1) -! endif -! enddo -! enddo -! endif end subroutine initialize_ice_shelf_boundary_channel -!BEGIN MJH -! subroutine initialize_ice_shelf_boundary(u_face_mask_bdry, v_face_mask_bdry, & -! u_flux_bdry_val, v_flux_bdry_val, u_bdry_val, v_bdry_val, h_bdry_val, & -! hmask, G, US, PF ) - -! type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure -! real, dimension(SZIB_(G),SZJ_(G)), & -! intent(inout) :: u_face_mask_bdry !< A boundary-type mask at C-grid u faces -! real, dimension(SZIB_(G),SZJ_(G)), & -! intent(inout) :: u_flux_bdry_val !< The boundary thickness flux through -! !! C-grid u faces [L Z T-1 ~> m2 s-1]. -! real, dimension(SZI_(G),SZJB_(G)), & -! intent(inout) :: v_face_mask_bdry !< A boundary-type mask at C-grid v faces -! real, dimension(SZI_(G),SZJB_(G)), & -! intent(inout) :: v_flux_bdry_val !< The boundary thickness flux through -! !! C-grid v faces [L Z T-1 ~> m2 s-1]. -! real, dimension(SZIB_(G),SZJB_(G)), & -! intent(inout) :: u_bdry_val !< The zonal ice shelf velocity at open -! !! boundary vertices [L T-1 ~> m s-1]. -! real, dimension(SZIB_(G),SZJB_(G)), & -! intent(inout) :: v_bdry_val !< The meridional ice shelf velocity at open -! !! boundary vertices [L T-1 ~> m s-1]. -! real, dimension(SZDI_(G),SZDJ_(G)), & -! intent(inout) :: h_bdry_val !< The ice shelf thickness at open boundaries [Z ~> m] -! real, dimension(SZDI_(G),SZDJ_(G)), & -! intent(inout) :: hmask !< A mask indicating which tracer points are -! !! partly or fully covered by an ice-shelf -! type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors -! type(param_file_type), intent(in) :: PF !< A structure to parse for run-time parameters - -! character(len=40) :: mdl = "initialize_ice_shelf_boundary" ! This subroutine's name. -! character(len=200) :: config -! logical flux_bdry - -! call get_param(PF, mdl, "ICE_BOUNDARY_CONFIG", config, & -! "This specifies how the ice domain boundary is specified. "//& -! "valid values include CHANNEL, FILE and USER.", & -! fail_if_missing=.true.) -! call get_param(PF, mdl, "ICE_BOUNDARY_FLUX_CONDITION", flux_bdry, & -! "This specifies whether mass input is a dirichlet or "//& -! "flux condition", default=.true.) - -! select case ( trim(config) ) -! case ("CHANNEL") -! call initialize_ice_shelf_boundary_channel(u_face_mask_bdry, & -! v_face_mask_bdry, u_flux_bdry_val, v_flux_bdry_val, & -! u_bdry_val, v_bdry_val, h_bdry_val, hmask, G, & -! flux_bdry, PF) -! case ("FILE"); call MOM_error(FATAL,"MOM_initialize: "// & -! "Unrecognized topography setup "//trim(config)) -! case ("USER"); call MOM_error(FATAL,"MOM_initialize: "// & -! "Unrecognized topography setup "//trim(config)) -! case default ; call MOM_error(FATAL,"MOM_initialize: "// & -! "Unrecognized topography setup "//trim(config)) -! end select - -! end subroutine initialize_ice_shelf_boundary - -! subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_bdry, & -! u_flux_bdry_val, v_flux_bdry_val, u_bdry_val, v_bdry_val, h_bdry_val, & -! hmask, G, flux_bdry, US, PF ) - -! type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure -! real, dimension(SZIB_(G),SZJ_(G)), & -! intent(inout) :: u_face_mask_bdry !< A boundary-type mask at C-grid u faces -! real, dimension(SZIB_(G),SZJ_(G)), & -! intent(inout) :: u_flux_bdry_val !< The boundary thickness flux through -! !! C-grid u faces [L Z T-1 ~> m2 s-1]. -! real, dimension(SZI_(G),SZJB_(G)), & -! intent(inout) :: v_face_mask_bdry !< A boundary-type mask at C-grid v faces -! real, dimension(SZI_(G),SZJB_(G)), & -! intent(inout) :: v_flux_bdry_val !< The boundary thickness flux through -! !! C-grid v faces [L Z T-1 ~> m2 s-1]. -! real, dimension(SZIB_(G),SZJB_(G)), & -! intent(inout) :: u_bdry_val !< The zonal ice shelf velocity at open - !! boundary vertices [L T-1 ~> m s-1]. -! real, dimension(SZIB_(G),SZJB_(G)), & -! intent(inout) :: v_bdry_val !< The meridional ice shelf velocity at open - !! boundary vertices [L T-1 ~> m s-1]. -! real, dimension(SZDI_(G),SZDJ_(G)), & -! intent(inout) :: h_bdry_val !< The ice shelf thickness at open boundaries [Z ~> m] -! real, dimension(SZDI_(G),SZDJ_(G)), & -! intent(inout) :: hmask !< A mask indicating which tracer points are -! !! partly or fully covered by an ice-shelf -! logical, intent(in) :: flux_bdry !< If true, use mass fluxes as the boundary value. -! type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors -! type(param_file_type), intent(in) :: PF !< A structure to parse for run-time parameters - -! character(len=40) :: mdl = "initialize_ice_shelf_boundary_channel" ! This subroutine's name. -! integer :: i, j, isd, jsd, is, js, iegq, jegq, giec, gjec, gisc, gjsc, isc, jsc, iec, jec, ied, jed -! real :: input_thick ! The input ice shelf thickness [Z ~> m] -! real :: input_flux ! The input ice flux per unit length [L Z T-1 ~> m2 s-1] -! real :: lenlat, len_stress - -! call get_param(PF, mdl, "LENLAT", lenlat, fail_if_missing=.true.) - -! call get_param(PF, mdl, "INPUT_FLUX_ICE_SHELF", input_flux, & -! "volume flux at upstream boundary", & -! units="m2 s-1", default=0., scale=US%m_s_to_L_T*US%m_to_Z) -! call get_param(PF, mdl, "INPUT_THICK_ICE_SHELF", input_thick, & -! "flux thickness at upstream boundary", & -! units="m", default=1000., scale=US%m_to_Z) -! call get_param(PF, mdl, "LEN_SIDE_STRESS", len_stress, & -! "maximum position of no-flow condition in along-flow direction", & -! units="km", default=0.) - -! call MOM_mesg(mdl//": setting boundary") - -! isd = G%isd ; ied = G%ied -! jsd = G%jsd ; jed = G%jed -! isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec -! gisc = G%Domain%nihalo ; gjsc = G%Domain%njhalo -! giec = G%Domain%niglobal+gisc ; gjec = G%Domain%njglobal+gjsc - -! do j=jsd,jed -! do i=isd,ied - -! ! upstream boundary - set either dirichlet or flux condition - -! if ((i+G%idg_offset) == G%domain%nihalo+1) then -! if (flux_bdry) then -! u_face_mask_bdry(i-1,j) = 4.0 -! u_flux_bdry_val(i-1,j) = input_flux -! else -! hmask(i-1,j) = 3.0 -! h_bdry_val(i-1,j) = input_thick -! u_face_mask_bdry(i-1,j) = 3.0 -! u_bdry_val(i-1,j-1) = (1 - ((G%geoLatBu(i-1,j-1) - 0.5*lenlat)*2./lenlat)**2) * & -! 1.5 * input_flux / input_thick -! u_bdry_val(i-1,j) = (1 - ((G%geoLatBu(i-1,j) - 0.5*lenlat)*2./lenlat)**2) * & -! 1.5 * input_flux / input_thick -! endif -! endif - -! ! side boundaries: no flow - -! if (G%jdg_offset+j == gjsc+1) then !bot boundary -! if (len_stress == 0. .OR. G%geoLonCv(i,j-1) <= len_stress) then -! v_face_mask_bdry(i,j-1) = 0. -! else -! v_face_mask_bdry(i,j-1) = 1. -! endif -! elseif (G%jdg_offset+j == gjec) then !top boundary -! if (len_stress == 0. .OR. G%geoLonCv(i,j-1) <= len_stress) then -! v_face_mask_bdry(i,j) = 0. -! else -! v_face_mask_bdry(i,j) = 1. -! endif -! endif - -! ! downstream boundary - CFBC - -! if (i+G%idg_offset == giec) then -! u_face_mask_bdry(i,j) = 2.0 -! endif - -! enddo -! enddo - -!END MJH end subroutine initialize_ice_shelf_boundary_channel !> Initialize ice shelf flow from file subroutine initialize_ice_flow_from_file(u_shelf, v_shelf,ice_visc,float_cond,& @@ -642,7 +460,7 @@ subroutine initialize_ice_flow_from_file(u_shelf, v_shelf,ice_visc,float_cond,& ! call MOM_read_data(filename, trim(ushelf_varname), u_shelf, G%Domain, scale=1.0) !/(365.0*86400.0)) ! call MOM_read_data(filename,trim(vshelf_varname), v_shelf, G%Domain, scale=1.0) !/(365.0*86400.0)) - call MOM_read_data(filename,trim(ice_visc_varname), ice_visc, G%Domain, scale=1.0) !*(365.0*86400.0)) + call MOM_read_data(filename,trim(ice_visc_varname), ice_visc, G%Domain, scale=1.0) call MOM_read_data(filename,trim(floatfr_varname), float_cond, G%Domain, scale=1.) ! call get_param(PF, mdl, "ICE_BOUNDARY_CONFIG", config, & ! "This specifies how the ice domain boundary is specified", & From abc8fe46324f4f7711100da37ad4b6bec0162cb0 Mon Sep 17 00:00:00 2001 From: Olga Sergienko Date: Thu, 11 Mar 2021 12:30:24 -0500 Subject: [PATCH 17/17] Removed blocks of commented text and multiplications by 0 --- src/ice_shelf/MOM_ice_shelf.F90 | 6 +- src/ice_shelf/MOM_ice_shelf_dynamics.F90 | 69 +++++----------------- src/ice_shelf/MOM_ice_shelf_initialize.F90 | 39 ++---------- 3 files changed, 23 insertions(+), 91 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 5d2bc88b4c..b2cb9f9c29 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -711,15 +711,15 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step, CS) endif ! Melting has been computed, now is time to update thickness and mass with dynamic ice shelf - if (CS%active_shelf_dynamics) then !OVS 12/10/20 - call change_thickness_using_melt(ISS, G, US, US%s_to_T*time_step, fluxes, CS%density_ice, CS%debug) !OVS 12/10/20 + if (CS%active_shelf_dynamics) then + call change_thickness_using_melt(ISS, G, US, US%s_to_T*time_step, fluxes, CS%density_ice, CS%debug) if (CS%debug) then call hchksum(ISS%h_shelf, "h_shelf after change thickness using melt", G%HI, haloshift=0, scale=US%Z_to_m) call hchksum(ISS%mass_shelf, "mass_shelf after change thickness using melt", G%HI, haloshift=0, & scale=US%RZ_to_kg_m2) endif - endif !OVS 12/10/20 + endif if (CS%debug) call MOM_forcing_chksum("Before add shelf flux", fluxes, G, CS%US, haloshift=0) diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index 6b30b2d83d..8360530f21 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -424,7 +424,6 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ ! previously allocated for registration for restarts. if (active_shelf_dynamics) then - ! DNG allocate( CS%u_bdry_val(Isdq:Iedq,Jsdq:Jedq) ) ; CS%u_bdry_val(:,:) = 0.0 allocate( CS%v_bdry_val(Isdq:Iedq,Jsdq:Jedq) ) ; CS%v_bdry_val(:,:) = 0.0 allocate( CS%t_bdry_val(isd:ied,jsd:jed) ) ; CS%t_bdry_val(:,:) = -15.0 @@ -572,7 +571,6 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ if (new_sim) then call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: initialize ice velocity.") call update_OD_ffrac_uncoupled(CS, G, ISS%h_shelf(:,:)) - !call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf, iters, Time) call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) if (CS%id_u_shelf > 0) call post_data(CS%id_u_shelf, CS%u_shelf, CS%diag) if (CS%id_v_shelf > 0) call post_data(CS%id_v_shelf, CS%v_shelf,CS%diag) @@ -692,7 +690,6 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled if (update_ice_vel) then -! call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf, iters, Time) call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) endif @@ -1801,9 +1798,9 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) isd = G%isd ; jsd = G%jsd iegq = G%iegB ; jegq = G%jegB ! gisc = G%domain%nihalo+1 ; gjsc = G%domain%njhalo+1 - gisc = 0*G%domain%nihalo+1 ; gjsc = 0*G%domain%njhalo+1 + gisc = 1 ; gjsc = 1 ! giec = G%domain%niglobal+G%domain%nihalo ; gjec = G%domain%njglobal+G%domain%njhalo - giec = G%domain%niglobal+0*G%domain%nihalo ; gjec = G%domain%njglobal+0*G%domain%njhalo + giec = G%domain%niglobal ; gjec = G%domain%njglobal is = iscq - 1; js = jscq - 1 i_off = G%idg_offset ; j_off = G%jdg_offset @@ -2519,7 +2516,7 @@ subroutine apply_boundary_values(CS, ISS, G, US, time, Phisub, H_node, ice_visc, endif endif ; enddo ; enddo - call pass_vector(u_bdry_contr, v_bdry_contr, G%domain, TO_ALL, BGRID_NE) !OVS 02/19/21 + call pass_vector(u_bdry_contr, v_bdry_contr, G%domain, TO_ALL, BGRID_NE) end subroutine apply_boundary_values @@ -2563,26 +2560,11 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) n_g = CS%n_glen; eps_min = CS%eps_glen_min -! CS%ice_visc(:,:) = 0.0 -! ux(:,:) = 0.0; uy(:,:) = 0.0; vx(:,:) =0.0; vy(:,:) =0.0 -! eII(:,:) = (US%s_to_T**2 * (eps_min**2)) - Visc_coef = US%kg_m2s_to_RZ_T*US%m_to_L*US%Z_to_L*(CS%A_glen_isothermal)**(-1./CS%n_glen) !OVS '-' in the exponent -! call pass_vector(u_shlf, v_shlf, G%domain, TO_ALL, BGRID_NE) - do j=jsc-0*1,jec+0*1 - do i=isc-0*1,iec+0*1 -! do j=jsd+1,jed-1 !OVS 02/01/21 -! do i=isd+1,ied-1 !OVS 02/01/21 + Visc_coef = US%kg_m2s_to_RZ_T*US%m_to_L*US%Z_to_L*(CS%A_glen_isothermal)**(-1./CS%n_glen) + do j=jsc,jec + do i=isc,iec if ((ISS%hmask(i,j) == 1) .OR. (ISS%hmask(i,j) == 3)) then -! ux(i,j) = ((u_shlf(I,J) + u_shlf(I,J-1)) - (u_shlf(I-1,J) + u_shlf(I-1,J-1))) / (2*G%dxT(i,j)) -! vx(i,j) = ((v_shlf(I,J) + v_shlf(I,J-1)) - (v_shlf(I-1,J) + v_shlf(I-1,J-1))) / (2*G%dxT(i,j)) -! uy(i,j) = ((u_shlf(I,J) + u_shlf(I-1,J)) - (u_shlf(I,J-1) + u_shlf(I-1,J-1))) / (2*G%dyT(i,j)) -! vy(i,j) = ((v_shlf(I,J) + v_shlf(I-1,J)) - (v_shlf(I,J-1) + v_shlf(I-1,J-1))) / (2*G%dyT(i,j)) -! endif -! enddo -! enddo -! call pass_vector(ux, uy, G%domain, TO_ALL, BGRID_NE) -! call pass_vector(vx, vy, G%domain, TO_ALL, BGRID_NE) ux = ((u_shlf(I,J) + (u_shlf(I,J-1) + u_shlf(I,J+1))) - & (u_shlf(I-1,J) + (u_shlf(I-1,J-1) + u_shlf(I-1,J+1)))) / (3*G%dxT(i,j)) vx = ((v_shlf(I,J) + v_shlf(I,J-1) + v_shlf(I,J+1)) - & @@ -2591,17 +2573,8 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) (u_shlf(I,J-1) + (u_shlf(I-1,J-1) + u_shlf(I+1,J-1)))) / (3*G%dyT(i,j)) vy = ((v_shlf(I,J) + (v_shlf(I-1,J)+ v_shlf(I+1,J))) - & (v_shlf(I,J-1) + (v_shlf(I-1,J-1)+ v_shlf(I+1,J-1)))) / (3*G%dyT(i,j)) -! ux = ((u_shlf(I,J) + u_shlf(I,J-1)) - (u_shlf(I-1,J) + u_shlf(I-1,J-1))) / (2*G%dxT(i,j)) -! vx = ((v_shlf(I,J) + v_shlf(I,J-1)) - (v_shlf(I-1,J) + v_shlf(I-1,J-1))) / (2*G%dxT(i,j)) -! uy = ((u_shlf(I,J) + u_shlf(I-1,J)) - (u_shlf(I,J-1) + u_shlf(I-1,J-1))) / (2*G%dyT(i,j)) -! vy = ((v_shlf(I,J) + v_shlf(I-1,J)) - (v_shlf(I,J-1) + v_shlf(I-1,J-1))) / (2*G%dyT(i,j)) CS%ice_visc(i,j) = 0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * & (US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2))**((1.-n_g)/(2.*n_g)) -! CS%ice_visc(i,j) =1e14*(G%areaT(i,j) * ISS%h_shelf(i,j)) ! constant viscocity for debugging -! umid = ((u_shlf(I,J) + u_shlf(I-1,J-1)) + (u_shlf(I,J-1) + u_shlf(I-1,J))) * 0.25 -! vmid = ((v_shlf(I,J) + v_shlf(I-1,J-1)) + (v_shlf(I,J-1) + v_shlf(I-1,J))) * 0.25 -! unorm = sqrt(umid**2 + vmid**2 + eps_min**2*(G%dxT(i,j)**2 + G%dyT(i,j)**2)) -! CS%basal_traction(i,j) = G%areaT(i,j) * CS%C_basal_friction * (US%L_T_to_m_s*unorm)**(CS%n_basal_fric-1) endif enddo enddo @@ -2638,8 +2611,8 @@ subroutine calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) eps_min = CS%eps_glen_min - do j=jsd+1,jed!-1 OVS 02/01/21 - do i=isd+1,ied!-1 OVS 02/01/21 + do j=jsd+1,jed + do i=isd+1,ied if ((ISS%hmask(i,j) == 1) .OR. (ISS%hmask(i,j) == 3)) then umid = ((u_shlf(I,J) + u_shlf(I-1,J-1)) + (u_shlf(I,J-1) + u_shlf(I-1,J))) * 0.25 @@ -2949,7 +2922,6 @@ subroutine update_velocity_masks(CS, G, hmask, umask, vmask, u_face_mask, v_face endif do j=js,G%jed -! do j=js-1,G%jed !OVS change index do i=is,G%ied if (hmask(i,j) == 1) then @@ -2966,10 +2938,6 @@ subroutine update_velocity_masks(CS, G, hmask, umask, vmask, u_face_mask, v_face umask(I-1+k,J)=3. !vmask(I-1+k,J)=0. vmask(I-1+k,J)=3. - !u_face_mask(I-1+k,j-1)=3. -! umask(I-1+k,J-1:J)=3. -! vmask(I-1+k,J-1:J)=0. -! u_face_mask(I-1+k,j)=3. case (2) u_face_mask(I-1+k,j)=2. case (4) @@ -2977,9 +2945,9 @@ subroutine update_velocity_masks(CS, G, hmask, umask, vmask, u_face_mask, v_face vmask(I-1+k,J-1:J)=0. u_face_mask(I-1+k,j)=4. case (0) -! umask(I-1+k,J-1:J)=0. -! vmask(I-1+k,J-1:J)=0. -! u_face_mask(I-1+k,j)=0. + umask(I-1+k,J-1:J)=0. + vmask(I-1+k,J-1:J)=0. + u_face_mask(I-1+k,j)=0. case (1) ! stress free x-boundary umask(I-1+k,J-1:J)=0. case default @@ -2990,8 +2958,6 @@ subroutine update_velocity_masks(CS, G, hmask, umask, vmask, u_face_mask, v_face select case (int(CS%v_face_mask_bdry(i,J-1+k))) case (3) -! vmask(I-1:I,J-1+k)=3. -! umask(I-1:I,J-1+k)=0. vmask(I-1,J-1+k)=3. umask(I-1,J-1+k)=0. vmask(I,J-1+k)=3. @@ -3004,9 +2970,9 @@ subroutine update_velocity_masks(CS, G, hmask, umask, vmask, u_face_mask, v_face vmask(I-1:I,J-1+k)=0. v_face_mask(i,J-1+k)=4. case (0) -! umask(I-1:I,J-1+k)=0. -! vmask(I-1:I,J-1+k)=0. -! v_face_mask(i,J-1+k)=0. + umask(I-1:I,J-1+k)=0. + vmask(I-1:I,J-1+k)=0. + v_face_mask(i,J-1+k)=0. case (1) ! stress free y-boundary vmask(I-1:I,J-1+k)=0. case default @@ -3134,7 +3100,6 @@ subroutine ice_shelf_temp(CS, ISS, G, US, time_step, melt_rate, Time) intent(in) :: melt_rate !< basal melt rate [R Z T-1 ~> kg m-2 s-1] type(time_type), intent(in) :: Time !< The current model time -! 5/23/12 OVS ! This subroutine takes the velocity (on the Bgrid) and timesteps ! (HT)_t = - div (uHT) + (adot Tsurf -bdot Tbot) once and then calculates T=HT/H ! @@ -3170,12 +3135,6 @@ subroutine ice_shelf_temp(CS, ISS, G, US, time_step, melt_rate, Time) TH(i,j) = CS%t_shelf(i,j)*ISS%h_shelf(i,j) enddo ; enddo -! call enable_averages(time_step, Time, CS%diag) -! call pass_var(h_after_uflux, G%domain) -! call pass_var(h_after_vflux, G%domain) -! if (CS%id_h_after_uflux > 0) call post_data(CS%id_h_after_uflux, h_after_uflux, CS%diag) -! if (CS%id_h_after_vflux > 0) call post_data(CS%id_h_after_vflux, h_after_vflux, CS%diag) -! call disable_averaging(CS%diag) call ice_shelf_advect_temp_x(CS, G, time_step, ISS%hmask, TH, th_after_uflux) call ice_shelf_advect_temp_y(CS, G, time_step, ISS%hmask, th_after_uflux, th_after_vflux) diff --git a/src/ice_shelf/MOM_ice_shelf_initialize.F90 b/src/ice_shelf/MOM_ice_shelf_initialize.F90 index ff05ed7c6a..1f5ddf909b 100644 --- a/src/ice_shelf/MOM_ice_shelf_initialize.F90 +++ b/src/ice_shelf/MOM_ice_shelf_initialize.F90 @@ -16,7 +16,6 @@ module MOM_ice_shelf_initialize #include -!MJHpublic initialize_ice_shelf_boundary, initialize_ice_thickness public initialize_ice_thickness public initialize_ice_shelf_boundary_channel public initialize_ice_flow_from_file @@ -132,10 +131,6 @@ subroutine initialize_ice_thickness_from_file(h_shelf, area_shelf_h, hmask, G, U call MOM_read_data(filename, trim(thickness_varname), h_shelf, G%Domain, scale=US%m_to_Z) call MOM_read_data(filename,trim(area_varname), area_shelf_h, G%Domain, scale=US%m_to_L**2) -! call get_param(PF, mdl, "ICE_BOUNDARY_CONFIG", config, & -! "This specifies how the ice domain boundary is specified", & -! fail_if_missing=.true.) - isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec do j=jsc,jec @@ -228,7 +223,6 @@ subroutine initialize_ice_thickness_channel(h_shelf, area_shelf_h, hmask, G, US, if (G%geoLonCu(i-1,j) >= edge_pos) then ! Everything past the edge is open ocean. -! mass_shelf(i,j) = 0.0 area_shelf_h(i,j) = 0.0 hmask (i,j) = 0.0 h_shelf (i,j) = 0.0 @@ -244,11 +238,7 @@ subroutine initialize_ice_thickness_channel(h_shelf, area_shelf_h, hmask, G, US, if (G%geoLonT(i,j) > slope_pos) then h_shelf(i,j) = min_draft -! mass_shelf(i,j) = Rho_ocean * min_draft else -! mass_shelf(i,j) = Rho_ocean * (min_draft + & -! (CS%max_draft - CS%min_draft) * & -! min(1.0, (c1*(slope_pos - G%geoLonT(i,j)))**2) ) h_shelf(i,j) = (min_draft + & (max_draft - min_draft) * & min(1.0, (c1*(slope_pos - G%geoLonT(i,j)))**2) ) @@ -301,7 +291,7 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b intent(inout) :: hmask !< A mask indicating which tracer points are !! partly or fully covered by an ice-shelf real, dimension(SZDI_(G),SZDJ_(G)), & - intent(inout) :: h_shelf !< Ice-shelf thickness OVS 11/10/20 + intent(inout) :: h_shelf !< Ice-shelf thickness type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors type(param_file_type), intent(in) :: PF !< A structure to parse for run-time parameters @@ -340,26 +330,16 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b giec = G%Domain%niglobal+gisc ; gjec = G%Domain%njglobal+gjsc !---------b.c.s based on geopositions ----------------- - ! do j=jsc-1,jec+1 - do j=jsc-0*1,jec+1 + do j=jsc,jec+1 do i=isc-1,iec+1 ! upstream boundary - set either dirichlet or flux condition if (G%geoLonBu(i,j) == westlon) then - ! if (flux_bdry) then - ! u_face_mask_bdry(i-1,j) = 4.0 - ! u_flux_bdry_val(i-1,j) = input_flux - ! else hmask(i+1,j) = 3.0 - ! hmask(i,j) = 3.0 h_bdry_val(i+1,j) = h_shelf(i+1,j) - ! h_bdry_val(i,j) = h_shelf(i,j) thickness_bdry_val(i+1,j) = h_bdry_val(i+0*1,j) u_face_mask_bdry(i+1,j) = 3.0 u_bdry_val(i+1,j) = input_vel*(1-16.0*((G%geoLatBu(i-1,j)/lenlat-0.5))**4) !velocity distribution - ! u_bdry_val(i+1,j) = (1 - ((G%geoLatBu(i,j) - 0.5*lenlat)*2./lenlat)**2) * & - ! 1.5 * input_flux / input_thick - ! endif endif @@ -367,7 +347,6 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b if (G%geoLatBu(i,j-1) == southlat) then !bot boundary if (len_stress == 0. .OR. G%geoLonCv(i,j) <= len_stress) then v_face_mask_bdry(i,j+1) = 0. - ! u_face_mask_bdry(i,j-1) = 3. u_face_mask_bdry(i,j) = 3. u_bdry_val(i,j) = 0. v_bdry_val(i,j) = 0. @@ -382,11 +361,8 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b v_face_mask_bdry(i,j-1) = 0. u_face_mask_bdry(i,j-1) = 3. else - !v_face_mask_bdry(i,j-1) = 1. v_face_mask_bdry(i,j-1) = 3. u_face_mask_bdry(i,j-1) = 3. - !u_bdry_val(i,j) = 0. - !hmask(i,j) = 0.0 endif endif @@ -458,13 +434,10 @@ subroutine initialize_ice_flow_from_file(u_shelf, v_shelf,ice_visc,float_cond,& floatfr_varname = "float_frac" -! call MOM_read_data(filename, trim(ushelf_varname), u_shelf, G%Domain, scale=1.0) !/(365.0*86400.0)) -! call MOM_read_data(filename,trim(vshelf_varname), v_shelf, G%Domain, scale=1.0) !/(365.0*86400.0)) - call MOM_read_data(filename,trim(ice_visc_varname), ice_visc, G%Domain, scale=1.0) - call MOM_read_data(filename,trim(floatfr_varname), float_cond, G%Domain, scale=1.) -! call get_param(PF, mdl, "ICE_BOUNDARY_CONFIG", config, & -! "This specifies how the ice domain boundary is specified", & -! fail_if_missing=.true.) + call MOM_read_data(filename, trim(ushelf_varname), u_shelf, G%Domain, scale=1.0) + call MOM_read_data(filename,trim(vshelf_varname), v_shelf, G%Domain, scale=1.0) + call MOM_read_data(filename,trim(ice_visc_varname), ice_visc, G%Domain, scale=1.0) + call MOM_read_data(filename,trim(floatfr_varname), float_cond, G%Domain, scale=1.) isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec