diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index debc63cb46..11b1eb5c16 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -70,7 +70,7 @@ module MOM_dynamics_split_RK2 use MOM_vert_friction, only : updateCFLtruncationValue, vertFPmix use MOM_verticalGrid, only : verticalGrid_type, get_thickness_units use MOM_verticalGrid, only : get_flux_units, get_tr_flux_units -use MOM_wave_interface, only: wave_parameters_CS, Stokes_PGF +use MOM_wave_interface, only : wave_parameters_CS, Stokes_PGF use MOM_CVMix_KPP, only : KPP_get_BLD, KPP_CS use MOM_energetic_PBL, only : energetic_PBL_get_MLD, energetic_PBL_CS use MOM_diabatic_driver, only : diabatic_CS, extract_diabatic_member @@ -284,16 +284,16 @@ module MOM_dynamics_split_RK2 contains !> RK2 splitting for time stepping MOM adiabatic dynamics -subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_surf_begin, p_surf_end, & - uh, vh, uhtr, vhtr, eta_av, G, GV, US, CS, calc_dtbt, VarMix, & - MEKE, thickness_diffuse_CSp, pbv, Waves) +subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, forces, & + p_surf_begin, p_surf_end, uh, vh, uhtr, vhtr, eta_av, G, GV, US, CS, & + calc_dtbt, VarMix, MEKE, thickness_diffuse_CSp, pbv, Waves) type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & - target, intent(inout) :: u !< Zonal velocity [L T-1 ~> m s-1] + target, intent(inout) :: u_inst !< Instantaneous zonal velocity [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & - target, intent(inout) :: v !< Meridional velocity [L T-1 ~> m s-1] + target, intent(inout) :: v_inst !< Instantaneous meridional velocity [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(inout) :: h !< Layer thickness [H ~> m or kg m-2] type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic type @@ -355,9 +355,9 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s ! height or column mass [H T-1 ~> m s-1 or kg m-2 s-1] real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: u_old_rad_OBC ! The starting zonal velocities, which are - ! saved for use in the Flather open boundary condition code [L T-1 ~> m s-1] + ! saved for use in the radiation open boundary condition code [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: v_old_rad_OBC ! The starting meridional velocities, which are - ! saved for use in the Flather open boundary condition code [L T-1 ~> m s-1] + ! saved for use in the radiation open boundary condition code [L T-1 ~> m s-1] ! GMM, TODO: make these allocatable? real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uold ! u-velocity before vert_visc is applied, for fpmix @@ -423,8 +423,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s call updateCFLtruncationValue(Time_local, CS%vertvisc_CSp, US) if (CS%debug) then - call MOM_state_chksum("Start predictor ", u, v, h, uh, vh, G, GV, US, symmetric=sym) - call check_redundant("Start predictor u ", u, v, G, unscale=US%L_T_to_m_s) + call MOM_state_chksum("Start predictor ", u_inst, v_inst, h, uh, vh, G, GV, US, symmetric=sym) + call check_redundant("Start predictor u ", u_inst, v_inst, G, unscale=US%L_T_to_m_s) call check_redundant("Start predictor uh ", uh, vh, G, unscale=GV%H_to_MKS*US%L_to_m**2*US%s_to_T) endif @@ -475,7 +475,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s call create_group_pass(CS%pass_hp_uv, u_av, v_av, G%Domain, halo=max(2,obc_stencil)) call create_group_pass(CS%pass_hp_uv, uh(:,:,:), vh(:,:,:), G%Domain, halo=max(2,obc_stencil)) - call create_group_pass(CS%pass_uv, u, v, G%Domain, halo=max(2,cont_stencil)) + call create_group_pass(CS%pass_uv, u_inst, v_inst, G%Domain, halo=max(2,cont_stencil)) call create_group_pass(CS%pass_h, h, G%Domain, halo=max(2,cont_stencil)) call create_group_pass(CS%pass_av_uvh, u_av, v_av, G%Domain, halo=max(2,obc_stencil)) call create_group_pass(CS%pass_av_uvh, uh(:,:,:), vh(:,:,:), G%Domain, halo=max(2,obc_stencil)) @@ -503,7 +503,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s if (Use_Stokes_PGF) Use_Stokes_PGF = Waves%Stokes_PGF if (Use_Stokes_PGF) then call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1) - call Stokes_PGF(G, GV, US, dz, u, v, CS%PFu_Stokes, CS%PFv_Stokes, Waves) + call Stokes_PGF(G, GV, US, dz, u_inst, v_inst, CS%PFu_Stokes, CS%PFv_Stokes, Waves) ! We are adding Stokes_PGF to hydrostatic PGF here. The diag PFu/PFv ! will therefore report the sum total PGF and we avoid other @@ -576,15 +576,15 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s !$OMP parallel do default(shared) do k=1,nz do j=js,je ; do I=Isq,Ieq - up(I,j,k) = G%mask2dCu(I,j) * (u(I,j,k) + dt * u_bc_accel(I,j,k)) + up(I,j,k) = G%mask2dCu(I,j) * (u_inst(I,j,k) + dt * u_bc_accel(I,j,k)) enddo ; enddo do J=Jsq,Jeq ; do i=is,ie - vp(i,J,k) = G%mask2dCv(i,J) * (v(i,J,k) + dt * v_bc_accel(i,J,k)) + vp(i,J,k) = G%mask2dCv(i,J) * (v_inst(i,J,k) + dt * v_bc_accel(i,J,k)) enddo ; enddo enddo call enable_averages(dt, Time_local, CS%diag) - call set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS%set_visc_CSp) + call set_viscous_ML(u_inst, v_inst, h, tv, forces, visc, dt, G, GV, US, CS%set_visc_CSp) call disable_averaging(CS%diag) if (CS%debug) then @@ -628,7 +628,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s ! u_accel_bt = layer accelerations due to barotropic solver if (associated(CS%BT_cont) .or. CS%BT_use_layer_fluxes) then call cpu_clock_begin(id_clock_continuity) - call continuity(u, v, h, hp, uh_in, vh_in, dt, G, GV, US, CS%continuity_CSp, CS%OBC, pbv, & + call continuity(u_inst, v_inst, h, hp, uh_in, vh_in, dt, G, GV, US, CS%continuity_CSp, CS%OBC, pbv, & visc_rem_u=CS%visc_rem_u, visc_rem_v=CS%visc_rem_v, BT_cont=CS%BT_cont) call cpu_clock_end(id_clock_continuity) if (BT_cont_BT_thick) then @@ -639,7 +639,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s endif if (CS%BT_use_layer_fluxes) then - uh_ptr => uh_in ; vh_ptr => vh_in; u_ptr => u ; v_ptr => v + uh_ptr => uh_in ; vh_ptr => vh_in; u_ptr => u_inst ; v_ptr => v_inst endif call cpu_clock_begin(id_clock_btstep) @@ -647,7 +647,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s if (showCallTree) call callTree_enter("btstep(), MOM_barotropic.F90") ! This is the predictor step call to btstep. ! The CS%ADp argument here stores the weights for certain integrated diagnostics. - call btstep(u, v, eta, dt, u_bc_accel, v_bc_accel, forces, CS%pbce, CS%eta_PF, u_av, v_av, & + call btstep(u_inst, v_inst, eta, dt, u_bc_accel, v_bc_accel, forces, CS%pbce, CS%eta_PF, u_av, v_av, & CS%u_accel_bt, CS%v_accel_bt, eta_pred, CS%uhbt, CS%vhbt, G, GV, US, & CS%barotropic_CSp, CS%visc_rem_u, CS%visc_rem_v, SpV_avg, CS%ADp, CS%OBC, CS%BT_cont, & eta_PF_start, taux_bot, tauy_bot, uh_ptr, vh_ptr, u_ptr, v_ptr) @@ -661,11 +661,11 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s !$OMP parallel do default(shared) do k=1,nz do J=Jsq,Jeq ; do i=is,ie - vp(i,J,k) = G%mask2dCv(i,J) * (v(i,J,k) + dt_pred * & + vp(i,J,k) = G%mask2dCv(i,J) * (v_inst(i,J,k) + dt_pred * & (v_bc_accel(i,J,k) + CS%v_accel_bt(i,J,k))) enddo ; enddo do j=js,je ; do I=Isq,Ieq - up(I,j,k) = G%mask2dCu(I,j) * (u(I,j,k) + dt_pred * & + up(I,j,k) = G%mask2dCu(I,j) * (u_inst(I,j,k) + dt_pred * & (u_bc_accel(I,j,k) + CS%u_accel_bt(I,j,k))) enddo ; enddo enddo @@ -679,7 +679,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s ! call MOM_state_chksum("Predictor 1", up, vp, h, uh, vh, G, GV, US, haloshift=1) call MOM_accel_chksum("Predictor accel", CS%CAu_pred, CS%CAv_pred, CS%PFu, CS%PFv, & CS%diffu, CS%diffv, G, GV, US, CS%pbce, CS%u_accel_bt, CS%v_accel_bt, symmetric=sym) - call MOM_state_chksum("Predictor 1 init", u, v, h, uh, vh, G, GV, US, haloshift=1, & + call MOM_state_chksum("Predictor 1 init", u_inst, v_inst, h, uh, vh, G, GV, US, haloshift=1, & symmetric=sym) call check_redundant("Predictor 1 up", up, vp, G, unscale=US%L_T_to_m_s) call check_redundant("Predictor 1 uh", uh, vh, G, unscale=GV%H_to_MKS*US%L_to_m**2*US%s_to_T) @@ -809,7 +809,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s if (Use_Stokes_PGF) Use_Stokes_PGF = Waves%Stokes_PGF if (Use_Stokes_PGF) then call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1) - call Stokes_PGF(G, GV, US, dz, u, v, CS%PFu_Stokes, CS%PFv_Stokes, Waves) + call Stokes_PGF(G, GV, US, dz, u_inst, v_inst, CS%PFu_Stokes, CS%PFv_Stokes, Waves) if (.not.Waves%Passive_Stokes_PGF) then do k=1,nz do j=js,je ; do I=Isq,Ieq @@ -899,7 +899,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s if (showCallTree) call callTree_enter("btstep(), MOM_barotropic.F90") ! This is the corrector step call to btstep. - call btstep(u, v, eta, dt, u_bc_accel, v_bc_accel, forces, CS%pbce, CS%eta_PF, u_av, v_av, & + call btstep(u_inst, v_inst, eta, dt, u_bc_accel, v_bc_accel, forces, CS%pbce, CS%eta_PF, u_av, v_av, & CS%u_accel_bt, CS%v_accel_bt, eta_pred, CS%uhbt, CS%vhbt, G, GV, US, & CS%barotropic_CSp, CS%visc_rem_u, CS%visc_rem_v, SpV_avg, CS%ADp, CS%OBC, CS%BT_cont, & eta_PF_start, taux_bot, tauy_bot, uh_ptr, vh_ptr, u_ptr, v_ptr, etaav=eta_av) @@ -920,22 +920,22 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s !$OMP parallel do default(shared) do k=1,nz do j=js,je ; do I=Isq,Ieq - u(I,j,k) = G%mask2dCu(I,j) * (u(I,j,k) + dt * & + u_inst(I,j,k) = G%mask2dCu(I,j) * (u_inst(I,j,k) + dt * & (u_bc_accel(I,j,k) + CS%u_accel_bt(I,j,k))) enddo ; enddo do J=Jsq,Jeq ; do i=is,ie - v(i,J,k) = G%mask2dCv(i,J) * (v(i,J,k) + dt * & + v_inst(i,J,k) = G%mask2dCv(i,J) * (v_inst(i,J,k) + dt * & (v_bc_accel(i,J,k) + CS%v_accel_bt(i,J,k))) enddo ; enddo enddo call cpu_clock_end(id_clock_mom_update) if (CS%debug) then - call uvchksum("Corrector 1 [uv]", u, v, G%HI, haloshift=0, symmetric=sym, scale=US%L_T_to_m_s) + call uvchksum("Corrector 1 [uv]", u_inst, v_inst, G%HI, haloshift=0, symmetric=sym, scale=US%L_T_to_m_s) call hchksum(h, "Corrector 1 h", G%HI, haloshift=1, scale=GV%H_to_MKS) call uvchksum("Corrector 1 [uv]h", uh, vh, G%HI, haloshift=2, & symmetric=sym, scale=GV%H_to_MKS*US%L_to_m**2*US%s_to_T) - ! call MOM_state_chksum("Corrector 1", u, v, h, uh, vh, G, GV, US, haloshift=1) + ! call MOM_state_chksum("Corrector 1", u_inst, v_inst, h, uh, vh, G, GV, US, haloshift=1) call MOM_accel_chksum("Corrector accel", CS%CAu, CS%CAv, CS%PFu, CS%PFv, & CS%diffu, CS%diffv, G, GV, US, CS%pbce, CS%u_accel_bt, CS%v_accel_bt, & symmetric=sym) @@ -951,26 +951,26 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s do k = 1, nz do j = js , je do I = Isq, Ieq - uold(I,j,k) = u(I,j,k) + uold(I,j,k) = u_inst(I,j,k) enddo enddo do J = Jsq, Jeq do i = is, ie - vold(i,J,k) = v(i,J,k) + vold(i,J,k) = v_inst(i,J,k) enddo enddo enddo endif call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1) - call vertvisc_coef(u, v, h, dz, forces, visc, tv, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC, VarMix) - call vertvisc(u, v, h, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp, G, GV, US, & + call vertvisc_coef(u_inst, v_inst, h, dz, forces, visc, tv, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC, VarMix) + call vertvisc(u_inst, v_inst, h, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp, G, GV, US, & CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot,waves=waves) if (CS%fpmix) then - call vertFPmix(u, v, uold, vold, hbl, h, forces, dt, & + call vertFPmix(u_inst, v_inst, uold, vold, hbl, h, forces, dt, & G, GV, US, CS%vertvisc_CSp, CS%OBC) - call vertvisc(u, v, h, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp, G, GV, US, & + call vertvisc(u_inst, v_inst, h, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp, G, GV, US, & CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves) endif @@ -1000,7 +1000,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s ! h = h + dt * div . uh ! u_av and v_av adjusted so their mass transports match uhbt and vhbt. call cpu_clock_begin(id_clock_continuity) - call continuity(u, v, h, h, uh, vh, dt, G, GV, US, CS%continuity_CSp, CS%OBC, pbv, & + call continuity(u_inst, v_inst, h, h, uh, vh, dt, G, GV, US, CS%continuity_CSp, CS%OBC, pbv, & CS%uhbt, CS%vhbt, CS%visc_rem_u, CS%visc_rem_v, u_av, v_av) call cpu_clock_end(id_clock_continuity) call do_group_pass(CS%pass_h, G%Domain, clock=id_clock_pass) @@ -1016,7 +1016,9 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s endif if (associated(CS%OBC)) then - call radiation_open_bdry_conds(CS%OBC, u, u_old_rad_OBC, v, v_old_rad_OBC, G, GV, US, dt) + !### I suspect that there is a bug here when u_inst is compared with a previous value of u_av + ! to estimate the dominant outward group velocity, but a fix is not available yet. + call radiation_open_bdry_conds(CS%OBC, u_inst, u_old_rad_OBC, v_inst, v_old_rad_OBC, G, GV, US, dt) endif ! h_av = (h_in + h_out)/2 . Going in to this line, h_av = h_in. @@ -1156,7 +1158,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s if (CS%id_deta_dt > 0) call post_data(CS%id_deta_dt, deta_dt, CS%diag) if (CS%debug) then - call MOM_state_chksum("Corrector ", u, v, h, uh, vh, G, GV, US, symmetric=sym) + call MOM_state_chksum("Corrector ", u_inst, v_inst, h, uh, vh, G, GV, US, symmetric=sym) call uvchksum("Corrector avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, scale=US%L_T_to_m_s) call hchksum(h_av, "Corrector avg h", G%HI, haloshift=1, scale=GV%H_to_MKS) ! call MOM_state_chksum("Corrector avg ", u_av, v_av, h_av, uh, vh, G, GV, US) @@ -1471,8 +1473,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param ! Accel_diag%u_accel_bt => CS%u_accel_bt ; Accel_diag%v_accel_bt => CS%v_accel_bt ! Accel_diag%u_av => CS%u_av ; Accel_diag%v_av => CS%v_av - id_clock_pass_init = cpu_clock_id('(Ocean init message passing)', & - grain=CLOCK_ROUTINE) + id_clock_pass_init = cpu_clock_id('(Ocean init message passing)', grain=CLOCK_ROUTINE) call continuity_init(Time, G, GV, US, param_file, diag, CS%continuity_CSp) cont_stencil = continuity_stencil(CS%continuity_CSp) @@ -1482,17 +1483,14 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param call PressureForce_init(Time, G, GV, US, param_file, diag, CS%PressureForce_CSp, & CS%SAL_CSp, CS%tides_CSp) call hor_visc_init(Time, G, GV, US, param_file, diag, CS%hor_visc, ADp=CS%ADp) - call vertvisc_init(MIS, Time, G, GV, US, param_file, diag, CS%ADp, dirs, & - ntrunc, CS%vertvisc_CSp) + call vertvisc_init(MIS, Time, G, GV, US, param_file, diag, CS%ADp, dirs, ntrunc, CS%vertvisc_CSp) CS%set_visc_CSp => set_visc - call updateCFLtruncationValue(Time, CS%vertvisc_CSp, US, & - activate=is_new_run(restart_CS) ) + call updateCFLtruncationValue(Time, CS%vertvisc_CSp, US, activate=is_new_run(restart_CS) ) if (associated(ALE_CSp)) CS%ALE_CSp => ALE_CSp if (associated(OBC)) then CS%OBC => OBC - if (OBC%ramp) call update_OBC_ramp(Time, CS%OBC, US, & - activate=is_new_run(restart_CS) ) + if (OBC%ramp) call update_OBC_ramp(Time, CS%OBC, US, activate=is_new_run(restart_CS) ) endif if (associated(update_OBC_CSp)) CS%update_OBC_CSp => update_OBC_CSp @@ -1515,8 +1513,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param do j=js,je ; do i=is,ie ; eta(i,j) = CS%eta(i,j) ; enddo ; enddo call barotropic_init(u, v, h, CS%eta, Time, G, GV, US, param_file, diag, & - CS%barotropic_CSp, restart_CS, calc_dtbt, CS%BT_cont, & - CS%SAL_CSp) + CS%barotropic_CSp, restart_CS, calc_dtbt, CS%BT_cont, CS%SAL_CSp) if (.not. query_initialized(CS%diffu, "diffu", restart_CS) .or. & .not. query_initialized(CS%diffv, "diffv", restart_CS)) then diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index 312fa43fe9..ec49081baf 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -84,12 +84,11 @@ module MOM_ice_shelf_dynamics real, pointer, dimension(:,:) :: t_shelf => NULL() !< Vertically integrated temperature in the ice shelf/stream, !! on corner-points (B grid) [C ~> degC] real, pointer, dimension(:,:) :: tmask => NULL() !< A mask on tracer points that is 1 where there is ice. - real, pointer, dimension(:,:) :: ice_visc => NULL() !< Glen's law ice viscosity (Pa s), + real, pointer, dimension(:,:,:) :: ice_visc => NULL() !< Glen's law ice viscosity (Pa s), !! in [R L2 T-1 ~> kg m-1 s-1]. - real, pointer, dimension(:,:,:) :: Ee => NULL() !< Glen's effective strain-rate ** (1-n)/(n) + !! at either 1 (cell-centered) or 4 quadrature points per cell real, pointer, dimension(:,:) :: AGlen_visc => NULL() !< Ice-stiffness parameter in Glen's law ice viscosity, !! often in [Pa-3 s-1] if n_Glen is 3. - real, pointer, dimension(:,:) :: thickness_bdry_val => NULL() !< The ice thickness at an inflowing boundary [Z ~> m]. real, pointer, dimension(:,:) :: u_bdry_val => NULL() !< The zonal ice velocity at inflowing boundaries !! [L yr-1 ~> m yr-1] real, pointer, dimension(:,:) :: v_bdry_val => NULL() !< The meridional ice velocity at inflowing boundaries @@ -114,6 +113,14 @@ module MOM_ice_shelf_dynamics real, pointer, dimension(:,:) :: ground_frac => NULL() !< Fraction of the time a cell is "exposed", i.e. the column !! thickness is below a threshold and interacting with the rock [nondim]. When this !! is 1, the ice-shelf is grounded + real, pointer, dimension(:,:) :: float_cond => NULL() !< If GL_regularize=true, indicates cells containing + !! the grounding line (float_cond=1) or not (float_cond=0) + real, pointer, dimension(:,:,:,:) :: Phi => NULL() !< The gradients of bilinear basis elements at Gaussian + !! 4 quadrature points surrounding the cell vertices [L-1 ~> m-1]. + real, pointer, dimension(:,:,:) :: PhiC => NULL() !< The gradients of bilinear basis elements at 1 cell-centered + !! quadrature point per cell [L-1 ~> m-1]. + real, pointer, dimension(:,:,:,:,:,:) :: Phisub => NULL() !< Quadrature structure weights at subgridscale + !! locations for finite element calculations [nondim] integer :: OD_rt_counter = 0 !< A counter of the number of contributions to OD_rt. real :: velocity_update_time_step !< The time interval over which to update the ice shelf velocity @@ -128,6 +135,13 @@ module MOM_ice_shelf_dynamics real :: density_ice !< A typical density of ice [R ~> kg m-3]. logical :: advect_shelf !< If true (default), advect ice shelf and evolve thickness + logical :: alternate_first_direction_IS !< If true, alternate whether the x- or y-direction + !! updates occur first in directionally split parts of the calculation. + integer :: first_direction_IS !< An integer that indicates which direction is + !! to be updated first in directionally split + !! parts of the ice sheet calculation (e.g. advection). + real :: first_dir_restart_IS = -1.0 !< A real copy of CS%first_direction_IS for use in restart files + integer :: visc_qps !< The number of quadrature points per cell (1 or 4) on which to calculate ice viscosity. character(len=40) :: ice_viscosity_compute !< Specifies whether the ice viscosity is computed internally !! according to Glen's flow law; is constant (for debugging purposes) !! or using observed strain rates and read from a file @@ -150,7 +164,7 @@ module MOM_ice_shelf_dynamics real :: eps_glen_min !< Min. strain rate to avoid infinite Glen's law viscosity, [T-1 ~> s-1]. real :: n_basal_fric !< Exponent in sliding law tau_b = C u^(m_slide) [nondim] logical :: CoulombFriction !< Use Coulomb friction law (Schoof 2005, Gagliardini et al 2007) - real :: CF_MinN !< Minimum Coulomb friction effective pressure [R L2 T-2 ~> Pa] + real :: CF_MinN !< Minimum Coulomb friction effective pressure [Pa] real :: CF_PostPeak !< Coulomb friction post peak exponent [nondim] real :: CF_Max !< Coulomb friction maximum coefficient [nondim] real :: density_ocean_avg !< A typical ocean density [R ~> kg m-3]. This does not affect ocean @@ -202,7 +216,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_bed_elev = -1, & - id_ground_frac = -1, id_col_thick = -1, id_OD_av = -1, & + id_ground_frac = -1, id_col_thick = -1, id_OD_av = -1, id_float_cond = -1, & id_u_mask = -1, id_v_mask = -1, id_ufb_mask =-1, id_vfb_mask = -1, id_t_mask = -1 !>@} ! ids for outputting intermediate thickness in advection subroutine (debugging) @@ -301,11 +315,28 @@ subroutine register_ice_shelf_dyn_restarts(G, US, param_file, CS, restart_CS) call get_param(param_file, mdl, "MISSING_SHELF_TEMPERATURE", T_shelf_missing, & "An ice shelf temperature to use where there is no ice shelf.",& units="degC", default=-10.0, scale=US%degC_to_C, do_not_log=.true.) + + call get_param(param_file, mdl, "NUMBER_OF_ICE_VISCOSITY_QUADRATURE_POINTS", CS%visc_qps, & + "Number of ice viscosity quadrature points. Either 1 (cell-centered) for 4", & + units="none", default=1) + if (CS%visc_qps/=1 .and. CS%visc_qps/=4) call MOM_error (FATAL, & + "NUMBER OF ICE_VISCOSITY_QUADRATURE_POINTS must be 1 or 4") + + call get_param(param_file, mdl, "FIRST_DIRECTION_IS", CS%first_direction_IS, & + "An integer that indicates which direction goes first "//& + "in parts of the code that use directionally split "//& + "updates (e.g. advection), with even numbers (or 0) used for x- first "//& + "and odd numbers used for y-first.", default=0) + call get_param(param_file, mdl, "ALTERNATE_FIRST_DIRECTION_IS", CS%alternate_first_direction_IS, & + "If true, after every advection call, alternate whether the x- or y- "//& + "direction advection updates occur first. "//& + "If this is true, FIRST_DIRECTION applies at the start of a new run or if "//& + "the next first direction can not be found in the restart file.", default=.false.) + allocate(CS%u_shelf(IsdB:IedB,JsdB:JedB), source=0.0) allocate(CS%v_shelf(IsdB:IedB,JsdB:JedB), source=0.0) allocate(CS%t_shelf(isd:ied,jsd:jed), source=T_shelf_missing) ! [C ~> degC] - allocate(CS%ice_visc(isd:ied,jsd:jed), source=0.0) - allocate(CS%Ee(isd:ied,jsd:jed,4), source=0.0) + allocate(CS%ice_visc(isd:ied,jsd:jed,CS%visc_qps), source=0.0) allocate(CS%AGlen_visc(isd:ied,jsd:jed), source=2.261e-25) ! [Pa-3 s-1] allocate(CS%basal_traction(isd:ied,jsd:jed), source=0.0) ! [R L3 T-1 ~> kg s-1] allocate(CS%C_basal_friction(isd:ied,jsd:jed), source=5.0e10) ! [Pa (m-1 s)^n_sliding] @@ -319,6 +350,7 @@ subroutine register_ice_shelf_dyn_restarts(G, US, param_file, CS, restart_CS) allocate(CS%u_face_mask_bdry(IsdB:IedB,JsdB:JedB), source=-2.0) allocate(CS%v_face_mask_bdry(IsdB:iedB,JsdB:JedB), source=-2.0) allocate(CS%h_bdry_val(isd:ied,jsd:jed), source=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", & @@ -349,6 +381,8 @@ subroutine register_ice_shelf_dyn_restarts(G, US, param_file, CS, restart_CS) "ice thickness at the boundary", "m", conversion=US%Z_to_m) call register_restart_field(CS%bed_elev, "bed elevation", .true., restart_CS, & "bed elevation", "m", conversion=US%Z_to_m) + call register_restart_field(CS%first_dir_restart_IS, "first_direction_IS", .false., restart_CS, & + "Indicator of the first direction in split ice shelf calculations.", "nondim") endif end subroutine register_ice_shelf_dyn_restarts @@ -466,7 +500,7 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ units="none", default=.false., fail_if_missing=.false.) call get_param(param_file, mdl, "CF_MinN", CS%CF_MinN, & "Minimum Coulomb friction effective pressure", & - units="Pa", default=1.0, scale=US%Pa_to_RL2_T2, fail_if_missing=.false.) + units="Pa", default=1.0, fail_if_missing=.false.) call get_param(param_file, mdl, "CF_PostPeak", CS%CF_PostPeak, & "Coulomb friction post peak exponent", & units="none", default=1.0, fail_if_missing=.false.) @@ -500,11 +534,14 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ "If true, advect ice shelf and evolve thickness", & default=.true.) call get_param(param_file, mdl, "ICE_VISCOSITY_COMPUTE", CS%ice_viscosity_compute, & - "If MODEL, compute ice viscosity internally at cell centers, if OBS read from a file,"//& - "If MODEL_QUADRATURE, compute at quadrature points (4 per element),"//& + "If MODEL, compute ice viscosity internally using 1 or 4 quadrature points,"//& + "if OBS read from a file,"//& "if CONSTANT a constant value (for debugging).", & default="MODEL") + if ((CS%visc_qps/=1) .and. (trim(CS%ice_viscosity_compute) /= "MODEL")) then + call MOM_error(FATAL, "NUMBER_OF_ICE_VISCOSITY_QUADRATURE_POINTS must be 1 unless ICE_VISCOSITY_COMPUTE==MODEL.") + endif call get_param(param_file, mdl, "INFLOW_SHELF_TEMPERATURE", T_shelf_bdry, & "A default ice shelf temperature to use for ice flowing in through "//& "open boundaries.", units="degC", default=-15.0, scale=US%degC_to_C) @@ -558,7 +595,6 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ if (active_shelf_dynamics) then allocate( CS%t_bdry_val(isd:ied,jsd:jed), source=T_shelf_bdry) ! [C ~> degC] - allocate( CS%thickness_bdry_val(isd:ied,jsd:jed), source=0.0) allocate( CS%u_face_mask(Isdq:Iedq,Jsdq:Jedq), source=0.0) allocate( CS%v_face_mask(Isdq:Iedq,Jsdq:Jedq), source=0.0) allocate( CS%u_flux_bdry_val(Isdq:Iedq,jsd:jed), source=0.0) @@ -566,6 +602,7 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ allocate( CS%umask(Isdq:Iedq,Jsdq:Jedq), source=-1.0) allocate( CS%vmask(Isdq:Iedq,Jsdq:Jedq), source=-1.0) allocate( CS%tmask(Isdq:Iedq,Jsdq:Jedq), source=-1.0) + allocate( CS%float_cond(isd:ied,jsd:jed)) CS%OD_rt_counter = 0 allocate( CS%OD_rt(isd:ied,jsd:jed), source=0.0) @@ -575,6 +612,24 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ allocate( CS%calve_mask(isd:ied,jsd:jed), source=0.0) endif + allocate(CS%Phi(1:8,1:4,isd:ied,jsd:jed), source=0.0) + do j=G%jsd,G%jed ; do i=G%isd,G%ied + call bilinear_shape_fn_grid(G, i, j, CS%Phi(:,:,i,j)) + enddo; enddo + + if (CS%GL_regularize) then + allocate(CS%Phisub(2,2,CS%n_sub_regularize,CS%n_sub_regularize,2,2), source=0.0) + call bilinear_shape_functions_subgrid(CS%Phisub, CS%n_sub_regularize) + endif + + if ((trim(CS%ice_viscosity_compute) == "MODEL") .and. CS%visc_qps==1) then + !for calculating viscosity and 1 cell-centered quadrature point per cell + allocate(CS%PhiC(1:8,G%isc:G%iec,G%jsc:G%jec), source=0.0) + do j=G%jsc,G%jec ; do i=G%isc,G%iec + call bilinear_shape_fn_grid_1qp(G, i, j, CS%PhiC(:,i,j)) + enddo; enddo + endif + CS%elapsed_velocity_time = 0.0 call update_velocity_masks(CS, G, ISS%hmask, CS%umask, CS%vmask, CS%u_face_mask, CS%v_face_mask) @@ -622,15 +677,13 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ endif call pass_var(CS%OD_av,G%domain, complete=.false.) - call pass_var(CS%ground_frac,G%domain, complete=.false.) - call pass_var(CS%ice_visc,G%domain, complete=.false.) + call pass_var(CS%ground_frac, G%domain, complete=.false.) call pass_var(CS%basal_traction, G%domain, complete=.false.) call pass_var(CS%AGlen_visc, G%domain, complete=.false.) call pass_var(CS%bed_elev, G%domain, complete=.false.) call pass_var(CS%C_basal_friction, G%domain, complete=.false.) - call pass_var(CS%h_bdry_val, G%domain, complete=.false.) - call pass_var(CS%thickness_bdry_val, G%domain, complete=.true.) - if (trim(CS%ice_viscosity_compute) == "MODEL_QUADRATURE") call pass_var(CS%Ee,G%domain) + call pass_var(CS%h_bdry_val, G%domain, complete=.true.) + call pass_var(CS%ice_visc, G%domain) call pass_vector(CS%u_shelf, CS%v_shelf, G%domain, TO_ALL, BGRID_NE, complete=.false.) call pass_vector(CS%u_bdry_val, CS%v_bdry_val, G%domain, TO_ALL, BGRID_NE, complete=.false.) @@ -639,6 +692,12 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ endif if (active_shelf_dynamics) then + if (CS%first_dir_restart_IS > -1.0) then + CS%first_direction_IS = modulo(NINT(CS%first_dir_restart_IS), 2) + else + CS%first_dir_restart_IS = real(modulo(CS%first_direction_IS, 2)) + endif + ! If we are calving to a mask, i.e. if a mask exists where a shelf cannot, read the mask from a file. if (CS%calve_to_mask) then call MOM_mesg(" MOM_ice_shelf.F90, initialize_ice_shelf: reading calving_mask") @@ -670,16 +729,15 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ call pass_var(CS%C_basal_friction, G%domain, complete=.false.) ! initialize ice-stiffness AGlen - call initialize_ice_AGlen(CS%AGlen_visc, G, US, param_file) + call initialize_ice_AGlen(CS%AGlen_visc, CS%ice_viscosity_compute, G, US, param_file) call pass_var(CS%AGlen_visc, G%domain, complete=.false.) !initialize boundary conditions call initialize_ice_shelf_boundary_from_file(CS%u_face_mask_bdry, CS%v_face_mask_bdry, & CS%u_bdry_val, CS%v_bdry_val, CS%umask, CS%vmask, CS%h_bdry_val, & - CS%thickness_bdry_val, ISS%hmask, ISS%h_shelf, G, US, param_file ) + ISS%hmask, ISS%h_shelf, G, US, param_file ) call pass_var(ISS%hmask, G%domain, complete=.false.) - call pass_var(CS%h_bdry_val, G%domain, complete=.false.) - call pass_var(CS%thickness_bdry_val, G%domain, complete=.true.) + call pass_var(CS%h_bdry_val, G%domain, complete=.true.) call pass_vector(CS%u_bdry_val, CS%v_bdry_val, G%domain, TO_ALL, BGRID_NE, complete=.false.) call pass_vector(CS%u_face_mask_bdry, CS%v_face_mask_bdry, G%domain, TO_ALL, BGRID_NE, complete=.false.) @@ -696,17 +754,18 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ '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('ice_shelf_model','v_shelf',CS%diag%axesB1, Time, & 'y-velocity of ice', 'm yr-1', conversion=365.0*86400.0*US%L_T_to_m_s) - ! I think that the conversion factors for the next two diagnostics are wrong. - RWH CS%id_taudx_shelf = register_diag_field('ice_shelf_model','taudx_shelf',CS%diag%axesB1, Time, & - 'x-driving stress of ice', 'kPa', conversion=1.e-3*US%RL2_T2_to_Pa) + 'x-driving stress of ice', 'kPa', conversion=1.e-3*US%RLZ_T2_to_Pa) CS%id_taudy_shelf = register_diag_field('ice_shelf_model','taudy_shelf',CS%diag%axesB1, Time, & - 'y-driving stress of ice', 'kPa', conversion=1.e-3*US%RL2_T2_to_Pa) + 'y-driving stress of ice', 'kPa', conversion=1.e-3*US%RLZ_T2_to_Pa) CS%id_u_mask = register_diag_field('ice_shelf_model','u_mask',CS%diag%axesB1, Time, & 'mask for u-nodes', 'none') CS%id_v_mask = register_diag_field('ice_shelf_model','v_mask',CS%diag%axesB1, Time, & 'mask for v-nodes', '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_float_cond = register_diag_field('ice_shelf_model','float_cond',CS%diag%axesT1, Time, & + 'sub-cell grounding cells', 'none') 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, & @@ -716,7 +775,7 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ 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) endif - call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: initialize ice velocity.") + if (new_sim) then call update_OD_ffrac_uncoupled(CS, G, ISS%h_shelf(:,:)) endif @@ -821,6 +880,10 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled ! if (CS%advect_shelf) then call ice_shelf_advect(CS, ISS, G, time_step, Time) + if (CS%alternate_first_direction_IS) then + CS%first_direction_IS = modulo(CS%first_direction_IS+1,2) + CS%first_dir_restart_IS = real(CS%first_direction_IS) + endif endif CS%elapsed_velocity_time = CS%elapsed_velocity_time + time_step if (CS%elapsed_velocity_time >= CS%velocity_update_time_step) update_ice_vel = .true. @@ -832,7 +895,6 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled CS%GL_couple=.false. endif - if (update_ice_vel) then call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) endif @@ -854,12 +916,14 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled call post_data(CS%id_taudy_shelf, taud_y, CS%diag) endif if (CS%id_ground_frac > 0) call post_data(CS%id_ground_frac, CS%ground_frac, CS%diag) + if (CS%id_float_cond > 0) call post_data(CS%id_float_cond, CS%float_cond, 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) then - ice_visc(:,:) = CS%ice_visc(:,:)*G%IareaT(:,:) - if (trim(CS%ice_viscosity_compute) == "MODEL_QUADRATURE") then - ice_visc(:,:) = ice_visc(:,:) * & - 0.25 * (CS%Ee(:,:,1) + CS%Ee(:,:,2) + CS%Ee(:,:,3) + CS%Ee(:,:,4)) + if (CS%visc_qps==4) then + ice_visc(:,:) = (0.25 * G%IareaT(:,:)) * & + ((CS%ice_visc(:,:,1) + CS%ice_visc(:,:,4)) + (CS%ice_visc(:,:,2) + CS%ice_visc(:,:,3))) + else + ice_visc(:,:) = CS%ice_visc(:,:,1)*G%IareaT(:,:) endif call post_data(CS%id_visc_shelf, ice_visc, CS%diag) endif @@ -945,11 +1009,11 @@ subroutine write_ice_shelf_energy(CS, G, US, mass, day, time_step) !calculate KE using cell-centered ice shelf velocity tmp1(:,:)=0.0 - KE_scale_factor = US%L_to_m**2 * US%RZ_to_kg_m2 * US%L_T_to_m_s**2 + KE_scale_factor = US%L_to_m**2 * (US%RZ_to_kg_m2 * US%L_T_to_m_s**2) do j=js,je ; do i=is,ie - tmp1(i,j) = KE_scale_factor * 0.03125 * G%areaT(i,j) * mass(i,j) * & - ((CS%u_shelf(I-1,J-1)+CS%u_shelf(I,J-1)+CS%u_shelf(I,J)+CS%u_shelf(I,J-1))**2 + & - (CS%v_shelf(I-1,J-1)+CS%v_shelf(I,J-1)+CS%v_shelf(I,J)+CS%v_shelf(I,J-1))**2) + tmp1(i,j) = (KE_scale_factor * 0.03125) * (G%areaT(i,j) * mass(i,j)) * & + (((CS%u_shelf(I-1,J-1)+CS%u_shelf(I,J))+(CS%u_shelf(I,J-1)+CS%u_shelf(I-1,J)))**2 + & + ((CS%v_shelf(I-1,J-1)+CS%v_shelf(I,J))+(CS%v_shelf(I,J-1)+CS%v_shelf(I-1,J)))**2) enddo; enddo KE_tot = reproducing_sum(tmp1, isr, ier, jsr, jer) @@ -958,7 +1022,7 @@ subroutine write_ice_shelf_energy(CS, G, US, mass, day, time_step) tmp1(:,:)=0.0 mass_scale_factor = US%L_to_m**2 * US%RZ_to_kg_m2 do j=js,je ; do i=is,ie - tmp1(i,j) = mass_scale_factor * mass(i,j) * G%areaT(i,j) + tmp1(i,j) = mass_scale_factor * (mass(i,j) * G%areaT(i,j)) enddo; enddo mass_tot = reproducing_sum(tmp1, isr, ier, jsr, jer) @@ -1034,7 +1098,7 @@ subroutine ice_shelf_advect(CS, ISS, G, time_step, Time) ! The flux overflows are included here. That is because they will be used to advect 3D scalars ! into partial cells - real, dimension(SZDI_(G),SZDJ_(G)) :: h_after_uflux, h_after_vflux ! Ice thicknesses [Z ~> m]. + real, dimension(SZDI_(G),SZDJ_(G)) :: h_after_flux1, h_after_flux2 ! Ice thicknesses [Z ~> m]. real, dimension(SZDIB_(G),SZDJ_(G)) :: uh_ice ! The accumulated zonal ice volume flux [Z L2 ~> m3] real, dimension(SZDI_(G),SZDJB_(G)) :: vh_ice ! The accumulated meridional ice volume flux [Z L2 ~> m3] type(loop_bounds_type) :: LB @@ -1046,37 +1110,39 @@ subroutine ice_shelf_advect(CS, ISS, G, time_step, Time) uh_ice(:,:) = 0.0 vh_ice(:,:) = 0.0 - h_after_uflux(:,:) = 0.0 - h_after_vflux(:,:) = 0.0 + h_after_flux1(:,:) = 0.0 + h_after_flux2(:,:) = 0.0 ! call MOM_mesg("MOM_ice_shelf.F90: ice_shelf_advect called") - do j=jsd,jed ; do i=isd,ied ; if (CS%thickness_bdry_val(i,j) /= 0.0) then - ISS%h_shelf(i,j) = CS%thickness_bdry_val(i,j) + do j=jsd,jed ; do i=isd,ied ; if (CS%h_bdry_val(i,j) /= 0.0) then + ISS%h_shelf(i,j) = CS%h_bdry_val(i,j) endif ; enddo ; enddo stencil = 2 - LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc-stencil ; LB%jeh = G%jec+stencil - if (LB%jsh < jsd) call MOM_error(FATAL, & - "ice_shelf_advect: Halo is too small for the ice thickness advection stencil.") - - 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) -! 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 - 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) -! 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) + if (modulo(CS%first_direction_IS,2)==0) then + !x first + LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc-stencil ; LB%jeh = G%jec+stencil + if (LB%jsh < jsd) call MOM_error(FATAL, & + "ice_shelf_advect: Halo is too small for the ice thickness advection stencil.") + call ice_shelf_advect_thickness_x(CS, G, LB, time_step, ISS%hmask, ISS%h_shelf, h_after_flux1, uh_ice) + call pass_var(h_after_flux1, G%domain) + LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc ; LB%jeh = G%jec + call ice_shelf_advect_thickness_y(CS, G, LB, time_step, ISS%hmask, h_after_flux1, h_after_flux2, vh_ice) + else + ! y first + LB%ish = G%isc-stencil ; LB%ieh = G%iec+stencil ; LB%jsh = G%jsc ; LB%jeh = G%jec + if (LB%ish < isd) call MOM_error(FATAL, & + "ice_shelf_advect: Halo is too small for the ice thickness advection stencil.") + call ice_shelf_advect_thickness_y(CS, G, LB, time_step, ISS%hmask, ISS%h_shelf, h_after_flux1, vh_ice) + call pass_var(h_after_flux1, G%domain) + LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc ; LB%jeh = G%jec + call ice_shelf_advect_thickness_x(CS, G, LB, time_step, ISS%hmask, h_after_flux1, h_after_flux2, uh_ice) + endif + call pass_var(h_after_flux2, G%domain) do j=jsd,jed do i=isd,ied - if (ISS%hmask(i,j) == 1) ISS%h_shelf(i,j) = h_after_vflux(i,j) + if (ISS%hmask(i,j) == 1) ISS%h_shelf(i,j) = h_after_flux2(i,j) enddo enddo @@ -1092,7 +1158,7 @@ subroutine ice_shelf_advect(CS, ISS, G, time_step, Time) endif do j=jsc,jec; do i=isc,iec - ISS%mass_shelf(i,j) = ISS%h_shelf(i,j) * ISS%area_shelf_h(i,j) * G%IareaT(i,j) * CS%density_ice + ISS%mass_shelf(i,j) = (ISS%h_shelf(i,j) * CS%density_ice) * (ISS%area_shelf_h(i,j) * G%IareaT(i,j)) enddo; enddo call pass_var(ISS%mass_shelf, G%domain, complete=.false.) @@ -1100,12 +1166,6 @@ subroutine ice_shelf_advect(CS, ISS, G, time_step, Time) call pass_var(ISS%area_shelf_h, G%domain, complete=.false.) call pass_var(ISS%hmask, G%domain, complete=.true.) - !call enable_averages(time_step, Time, CS%diag) - !if (CS%id_h_after_adv > 0) call post_data(CS%id_h_after_adv, ISS%h_shelf, CS%diag) - !call disable_averaging(CS%diag) - - !call change_thickness_using_melt(ISS, G, time_step, fluxes, CS%density_ice) - call update_velocity_masks(CS, G, ISS%hmask, CS%umask, CS%vmask, CS%u_face_mask, CS%v_face_mask) end subroutine ice_shelf_advect @@ -1134,40 +1194,30 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i real, dimension(SZDIB_(G),SZDJB_(G)) :: Au, Av ! The retarding lateral stress contributions [R L3 Z T-2 ~> kg m s-2] real, dimension(SZDIB_(G),SZDJB_(G)) :: u_last, v_last ! Previous velocities [L T-1 ~> m s-1] real, dimension(SZDIB_(G),SZDJB_(G)) :: H_node ! Ice shelf thickness at corners [Z ~> m]. - real, dimension(SZDI_(G),SZDJ_(G)) :: float_cond ! If GL_regularize=true, an array indicating where the ice - ! shelf is floating: 0 if floating, 1 if not + real, dimension(SZDI_(G),SZDJ_(G)) :: float_cond ! If GL_regularize=true, indicates cells containing + ! the grounding line (float_cond=1) or not (float_cond=0) real, dimension(SZDIB_(G),SZDJB_(G)) :: Normvec ! Used for convergence character(len=160) :: mesg ! The text of an error message integer :: conv_flag, i, j, k,l, iter - integer :: isdq, iedq, jsdq, jedq, isd, ied, jsd, jed, nodefloat, nsub + integer :: isdq, iedq, jsdq, jedq, isd, ied, jsd, jed, nodefloat real :: err_max, err_tempu, err_tempv, err_init, max_vel, tempu, tempv, Norm, PrevNorm real :: rhoi_rhow ! The density of ice divided by a typical water density [nondim] - real, pointer, dimension(:,:,:,:) :: Phi => NULL() ! The gradients of bilinear basis elements at Gaussian - ! quadrature points surrounding the cell vertices [L-1 ~> m-1]. - real, pointer, dimension(:,:,:,:,:,:) :: Phisub => NULL() ! Quadrature structure weights at subgridscale - ! locations for finite element calculations [nondim] integer :: Is_sum, Js_sum, Ie_sum, Je_sum ! Loop bounds for global sums or arrays starting at 1. - ! for GL interpolation - nsub = CS%n_sub_regularize - isdq = G%isdB ; iedq = G%iedB ; jsdq = G%jsdB ; jedq = G%jedB isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed rhoi_rhow = CS%density_ice / CS%density_ocean_avg taudx(:,:) = 0.0 ; taudy(:,:) = 0.0 - !u_bdry_cont(:,:) = 0.0 ; v_bdry_cont(:,:) = 0.0 Au(:,:) = 0.0 ; Av(:,:) = 0.0 ! need to make these conditional on GL interpolation - float_cond(:,:) = 0.0 ; H_node(:,:) = 0.0 + CS%float_cond(:,:) = 0.0 ; H_node(:,:) = 0.0 !CS%ground_frac(:,:) = 0.0 - allocate(Phisub(nsub,nsub,2,2,2,2), source=0.0) if (.not. CS%GL_couple) then do j=G%jsc,G%jec ; do i=G%isc,G%iec if (rhoi_rhow * ISS%h_shelf(i,j) - CS%bed_elev(i,j) > 0) then - if (CS%GL_regularize) float_cond(i,j) = 1.0 CS%ground_frac(i,j) = 1.0 CS%OD_av(i,j) =0.0 endif @@ -1197,34 +1247,24 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i endif enddo ; enddo if ((nodefloat > 0) .and. (nodefloat < 4)) then - float_cond(i,j) = 1.0 + CS%float_cond(i,j) = 1.0 CS%ground_frac(i,j) = 1.0 endif enddo ; enddo - call pass_var(float_cond, G%Domain, complete=.false.) - - call bilinear_shape_functions_subgrid(Phisub, nsub) + call pass_var(CS%float_cond, G%Domain, complete=.false.) endif - ! must prepare Phi - allocate(Phi(1:8,1:4,isd:ied,jsd:jed), source=0.0) - - do j=jsd,jed ; do i=isd,ied - 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 pass_var(CS%ice_visc, G%domain, complete=.false.) call calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) call pass_var(CS%basal_traction, G%domain, complete=.true.) - if (trim(CS%ice_viscosity_compute) == "MODEL_QUADRATURE") call pass_var(CS%Ee,G%domain) + call calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) + call pass_var(CS%ice_visc, G%domain) ! This makes sure basal stress is only applied when it is supposed to be if (CS%GL_regularize) then do j=G%jsd,G%jed ; do i=G%isd,G%ied - CS%basal_traction(i,j) = CS%basal_traction(i,j) * float_cond(i,j) + if (CS%ground_frac(i,j)/=1.0) CS%basal_traction(i,j) = 0.0 enddo ; enddo else do j=G%jsd,G%jed ; do i=G%isd,G%ied @@ -1233,25 +1273,21 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i endif if (CS%nonlin_solve_err_mode == 1) then - ! call apply_boundary_values(CS, ISS, G, US, time, Phisub, H_node, CS%ice_visc, & - ! CS%basal_traction, float_cond, rhoi_rhow, u_bdry_cont, v_bdry_cont) Au(:,:) = 0.0 ; Av(:,:) = 0.0 - call CG_action(CS, Au, Av, u_shlf, v_shlf, Phi, Phisub, CS%umask, CS%vmask, ISS%hmask, H_node, & - CS%ice_visc, float_cond, CS%bed_elev, CS%basal_traction, & + call CG_action(CS, Au, Av, u_shlf, v_shlf, CS%Phi, CS%Phisub, CS%umask, CS%vmask, ISS%hmask, H_node, & + CS%ice_visc, CS%float_cond, CS%bed_elev, 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, TO_ALL, BGRID_NE) err_init = 0 ; err_tempu = 0 ; err_tempv = 0 do J=G%IscB,G%JecB ; do I=G%IscB,G%IecB if (CS%umask(I,J) == 1) then - !err_tempu = ABS(Au(I,J) + u_bdry_cont(I,J) - taudx(I,J)) err_tempu = ABS(Au(I,J) - taudx(I,J)) if (err_tempu >= err_init) err_init = err_tempu endif if (CS%vmask(I,J) == 1) then - !err_tempv = ABS(Av(I,J) + v_bdry_cont(I,J) - taudy(I,J)) err_tempv = ABS(Av(I,J) - taudy(I,J)) if (err_tempv >= err_init) err_init = err_tempv endif @@ -1271,8 +1307,8 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i ! Include the edge if tile is at the southern bdry; Should add a test to avoid this if reentrant. if (G%jsc+G%jdg_offset==G%jsg) Js_sum = G%JscB + (1-G%JsdB) do J=G%jscB,G%jecB ; do I=G%iscB,G%iecB - if (CS%umask(I,J) == 1) Normvec(I,J) = Normvec(I,J) + u_shlf(I,J)*u_shlf(I,J) - if (CS%vmask(I,J) == 1) Normvec(I,J) = Normvec(I,J) + v_shlf(I,J)*v_shlf(I,J) + if (CS%umask(I,J) == 1) Normvec(I,J) = Normvec(I,J) + (u_shlf(I,J)**2 * US%L_T_to_m_s**2) + if (CS%vmask(I,J) == 1) Normvec(I,J) = Normvec(I,J) + (v_shlf(I,J)**2 * US%L_T_to_m_s**2) enddo; enddo Norm = reproducing_sum( Normvec, Is_sum, Ie_sum, Js_sum, Je_sum ) Norm = sqrt(Norm) @@ -1284,8 +1320,8 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i do iter=1,50 - call ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H_node, float_cond, & - ISS%hmask, conv_flag, iters, time, Phi, Phisub) + call ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H_node, CS%float_cond, & + ISS%hmask, conv_flag, iters, time, CS%Phi, CS%Phisub) if (CS%debug) then call qchksum(u_shlf, "u shelf", G%HI, haloshift=2, scale=US%L_T_to_m_s) @@ -1295,16 +1331,15 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i 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 pass_var(CS%ice_visc, G%domain, complete=.false.) call calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) call pass_var(CS%basal_traction, G%domain, complete=.true.) - if (trim(CS%ice_viscosity_compute) == "MODEL_QUADRATURE") call pass_var(CS%Ee,G%domain) + call calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) + call pass_var(CS%ice_visc, G%domain) ! makes sure basal stress is only applied when it is supposed to be if (CS%GL_regularize) then do j=G%jsd,G%jed ; do i=G%isd,G%ied - CS%basal_traction(i,j) = CS%basal_traction(i,j) * float_cond(i,j) + if (CS%ground_frac(i,j)/=1.0) CS%basal_traction(i,j) = 0.0 enddo ; enddo else do j=G%jsd,G%jed ; do i=G%isd,G%ied @@ -1313,15 +1348,11 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i endif if (CS%nonlin_solve_err_mode == 1) then - !u_bdry_cont(:,:) = 0 ; v_bdry_cont(:,:) = 0 - - ! call apply_boundary_values(CS, ISS, G, US, time, Phisub, H_node, CS%ice_visc, & - ! CS%basal_traction, float_cond, rhoi_rhow, u_bdry_cont, v_bdry_cont) Au(:,:) = 0 ; Av(:,:) = 0 - call CG_action(CS, Au, Av, u_shlf, v_shlf, Phi, Phisub, CS%umask, CS%vmask, ISS%hmask, H_node, & - CS%ice_visc, float_cond, CS%bed_elev, CS%basal_traction, & + call CG_action(CS, Au, Av, u_shlf, v_shlf, CS%Phi, CS%Phisub, CS%umask, CS%vmask, ISS%hmask, H_node, & + CS%ice_visc, CS%float_cond, CS%bed_elev, 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, TO_ALL, BGRID_NE) @@ -1330,12 +1361,10 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i do J=G%jscB,G%jecB ; do I=G%jscB,G%iecB if (CS%umask(I,J) == 1) then - !err_tempu = ABS(Au(I,J) + u_bdry_cont(I,J) - taudx(I,J)) err_tempu = ABS(Au(I,J) - taudx(I,J)) if (err_tempu >= err_max) err_max = err_tempu endif if (CS%vmask(I,J) == 1) then - !err_tempv = ABS(Av(I,J) + v_bdry_cont(I,J) - taudy(I,J)) err_tempv = ABS(Av(I,J) - taudy(I,J)) if (err_tempv >= err_max) err_max = err_tempv endif @@ -1372,8 +1401,8 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i elseif (CS%nonlin_solve_err_mode == 3) then PrevNorm=Norm; Norm=0.0; Normvec=0.0 do J=G%jscB,G%jecB ; do I=G%iscB,G%iecB - if (CS%umask(I,J) == 1) Normvec(I,J) = Normvec(I,J) + u_shlf(I,J)*u_shlf(I,J) - if (CS%vmask(I,J) == 1) Normvec(I,J) = Normvec(I,J) + v_shlf(I,J)*v_shlf(I,J) + if (CS%umask(I,J) == 1) Normvec(I,J) = Normvec(I,J) + (u_shlf(I,J)**2 * US%L_T_to_m_s**2) + if (CS%vmask(I,J) == 1) Normvec(I,J) = Normvec(I,J) + (v_shlf(I,J)**2 * US%L_T_to_m_s**2) enddo; enddo Norm = reproducing_sum( Normvec, Is_sum, Ie_sum, Js_sum, Je_sum ) Norm = sqrt(Norm) @@ -1393,8 +1422,6 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i call MOM_mesg(mesg) write(mesg,*) "ice_shelf_solve_outer: exiting nonlinear solve after ",iter," iterations" call MOM_mesg(mesg) - deallocate(Phi) - deallocate(Phisub) end subroutine ice_shelf_solve_outer @@ -1417,8 +1444,8 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H intent(in) :: H_node !< The ice shelf thickness at nodal (corner) !! points [Z ~> m]. real, dimension(SZDI_(G),SZDJ_(G)), & - intent(in) :: float_cond !< If GL_regularize=true, an array indicating where the ice - !! shelf is floating: 0 if floating, 1 if not + intent(in) :: float_cond !< If GL_regularize=true, indicates cells containing + !! the grounding line (float_cond=1) or not (float_cond=0) real, dimension(SZDI_(G),SZDJ_(G)), & intent(in) :: hmask !< A mask indicating which tracer points are !! partly or fully covered by an ice-shelf @@ -1473,7 +1500,7 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H Zu(:,:) = 0 ; Zv(:,:) = 0 ; DIAGu(:,:) = 0 ; DIAGv(:,:) = 0 Ru(:,:) = 0 ; Rv(:,:) = 0 ; Au(:,:) = 0 ; Av(:,:) = 0 ; RHSu(:,:) = 0 ; RHSv(:,:) = 0 - Du(:,:) = 0 ; Dv(:,:) = 0 !; ubd(:,:) = 0 ; vbd(:,:) = 0 + Du(:,:) = 0 ; Dv(:,:) = 0 dot_p1 = 0 ! Determine the loop limits for sums, bearing in mind that the arrays will be starting at 1. @@ -1487,16 +1514,13 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H ! Include the edge if tile is at the southern bdry; Should add a test to avoid this if reentrant. if (G%jsc+G%jdg_offset==G%jsg) Js_sum = G%JscB + (1-G%JsdB) - !call apply_boundary_values(CS, ISS, G, US, time, Phisub, H_node, CS%ice_visc, & - ! CS%basal_traction, float_cond, rhoi_rhow, ubd, vbd) - - RHSu(:,:) = taudx(:,:) !- ubd(:,:) - RHSv(:,:) = taudy(:,:) !- vbd(:,:) + RHSu(:,:) = taudx(:,:) + RHSv(:,:) = taudy(:,:) call pass_vector(RHSu, RHSv, G%domain, TO_ALL, BGRID_NE, complete=.false.) call matrix_diagonal(CS, G, US, float_cond, H_node, CS%ice_visc, CS%basal_traction, & - hmask, rhoi_rhow, Phisub, DIAGu, DIAGv) + hmask, rhoi_rhow, Phi, Phisub, DIAGu, DIAGv) call pass_vector(DIAGu, DIAGv, G%domain, TO_ALL, BGRID_NE, complete=.false.) @@ -1508,9 +1532,8 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H Ru(:,:) = (RHSu(:,:) - Au(:,:)) Rv(:,:) = (RHSv(:,:) - Av(:,:)) - - resid_scale = US%L_to_m**2*US%s_to_T*US%RZ_to_kg_m2*US%L_T_to_m_s**2 - resid2_scale = (US%RZ_to_kg_m2*US%L_to_m*US%L_T_to_m_s**2)**2 + resid_scale = (US%L_to_m**2*US%s_to_T)*(US%RZ_to_kg_m2*US%L_T_to_m_s**2) + resid2_scale = ((US%RZ_to_kg_m2*US%L_to_m)*US%L_T_to_m_s**2)**2 sum_vec(:,:) = 0.0 do j=jscq,jecq ; do i=iscq,iecq @@ -1518,7 +1541,7 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H if (CS%vmask(I,J) == 1) sum_vec(I,J) = sum_vec(I,J) + resid2_scale*Rv(I,J)**2 enddo ; enddo - dot_p1 = reproducing_sum( sum_vec, Js_sum, Ie_sum, Js_sum, Je_sum ) + dot_p1 = reproducing_sum( sum_vec, Is_sum, Ie_sum, Js_sum, Je_sum ) resid0 = sqrt(dot_p1) @@ -1654,7 +1677,8 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H dot_p1 = reproducing_sum( sum_vec, Is_sum, Ie_sum, Js_sum, Je_sum ) dot_p1 = sqrt(dot_p1) - if (dot_p1 <= CS%cg_tolerance * resid0) then + + if (dot_p1 <= (CS%cg_tolerance * resid0)) then iters = iter conv_flag = 1 exit @@ -1732,37 +1756,39 @@ subroutine ice_shelf_advect_thickness_x(CS, G, LB, time_step, hmask, h0, h_after do j=jsh,jeh ; do I=ish-1,ieh if (CS%u_face_mask(I,j) == 4.) then ! The flux itself is a specified boundary condition. - uh_ice(I,j) = time_step * G%dyCu(I,j) * CS%u_flux_bdry_val(I,j) + uh_ice(I,j) = (time_step * G%dyCu(I,j)) * CS%u_flux_bdry_val(I,j) elseif ((hmask(i,j) == 1 .or. hmask(i,j) == 3) .or. (hmask(i+1,j) == 1 .or. hmask(i+1,j) == 3)) then u_face = 0.5 * (CS%u_shelf(I,J-1) + CS%u_shelf(I,J)) h_face = 0.0 ! This will apply when the source cell is iceless or not fully ice covered. if (u_face > 0) then if (hmask(i,j) == 3) then ! This is a open boundary inflow from the west - h_face = CS%thickness_bdry_val(i,j) + h_face = CS%h_bdry_val(i,j) elseif (hmask(i,j) == 1) then ! There can be eastward flow through this face. - if ((hmask(i-1,j) == 1) .and. (hmask(i+1,j) == 1)) then + if ((hmask(i-1,j) == 1 .or. hmask(i-1,j) == 3) .and. & + (hmask(i+1,j) == 1 .or. hmask(i+1,j) == 3)) then slope_lim = slope_limiter(h0(i,j)-h0(i-1,j), h0(i+1,j)-h0(i,j)) ! This is a 2nd-order centered scheme with a slope limiter. We could try PPM here. - h_face = h0(i,j) - slope_lim * 0.5 * (h0(i,j)-h0(i+1,j)) + h_face = h0(i,j) - slope_lim * (0.5 * (h0(i,j)-h0(i+1,j))) else h_face = h0(i,j) endif endif else if (hmask(i+1,j) == 3) then ! This is a open boundary inflow from the east - h_face = CS%thickness_bdry_val(i+1,j) + h_face = CS%h_bdry_val(i+1,j) elseif (hmask(i+1,j) == 1) then - if ((hmask(i,j) == 1) .and. (hmask(i+2,j) == 1)) then + if ((hmask(i,j) == 1 .or. hmask(i,j) == 3) .and. & + (hmask(i+2,j) == 1 .or. hmask(i+2,j) == 3)) then slope_lim = slope_limiter(h0(i+1,j)-h0(i,j), h0(i+2,j)-h0(i+1,j)) - h_face = h0(i+1,j) - slope_lim * 0.5 * (h0(i+1,j)-h0(i,j)) + h_face = h0(i+1,j) - slope_lim * (0.5 * (h0(i+2,j)-h0(i+1,j))) else h_face = h0(i+1,j) endif endif endif - uh_ice(I,j) = time_step * G%dyCu(I,j) * u_face * h_face + uh_ice(I,j) = (time_step * G%dyCu(I,j)) * (u_face * h_face) else uh_ice(I,j) = 0.0 endif @@ -1811,37 +1837,39 @@ subroutine ice_shelf_advect_thickness_y(CS, G, LB, time_step, hmask, h0, h_after do J=jsh-1,jeh ; do i=ish,ieh if (CS%v_face_mask(i,J) == 4.) then ! The flux itself is a specified boundary condition. - vh_ice(i,J) = time_step * G%dxCv(i,J) * CS%v_flux_bdry_val(i,J) + vh_ice(i,J) = (time_step * G%dxCv(i,J)) * CS%v_flux_bdry_val(i,J) elseif ((hmask(i,j) == 1 .or. hmask(i,j) == 3) .or. (hmask(i,j+1) == 1 .or. hmask(i,j+1) == 3)) then v_face = 0.5 * (CS%v_shelf(I-1,J) + CS%v_shelf(I,J)) h_face = 0.0 ! This will apply when the source cell is iceless or not fully ice covered. if (v_face > 0) then if (hmask(i,j) == 3) then ! This is a open boundary inflow from the south - h_face = CS%thickness_bdry_val(i,j) - elseif (hmask(i,j) == 1) then ! There can be northtward flow through this face. - if ((hmask(i,j-1) == 1) .and. (hmask(i,j+1) == 1)) then + h_face = CS%h_bdry_val(i,j) + elseif (hmask(i,j) == 1) then ! There can be northward flow through this face. + if ((hmask(i,j-1) == 1 .or. hmask(i,j-1) == 3) .and. & + (hmask(i,j+1) == 1 .or. hmask(i,j+1) == 3)) then slope_lim = slope_limiter(h0(i,j)-h0(i,j-1), h0(i,j+1)-h0(i,j)) ! This is a 2nd-order centered scheme with a slope limiter. We could try PPM here. - h_face = h0(i,j) - slope_lim * 0.5 * (h0(i,j)-h0(i,j+1)) + h_face = h0(i,j) - slope_lim * (0.5 * (h0(i,j)-h0(i,j+1))) else h_face = h0(i,j) endif endif else if (hmask(i,j+1) == 3) then ! This is a open boundary inflow from the north - h_face = CS%thickness_bdry_val(i,j+1) + h_face = CS%h_bdry_val(i,j+1) elseif (hmask(i,j+1) == 1) then - if ((hmask(i,j) == 1) .and. (hmask(i,j+2) == 1)) then + if ((hmask(i,j) == 1 .or. hmask(i,j) == 3) .and. & + (hmask(i,j+2) == 1 .or. hmask(i,j+2) == 3)) then slope_lim = slope_limiter(h0(i,j+1)-h0(i,j), h0(i,j+2)-h0(i,j+1)) - h_face = h0(i,j+1) - slope_lim * 0.5 * (h0(i,j+1)-h0(i,j)) + h_face = h0(i,j+1) - slope_lim * (0.5 * (h0(i,j+2)-h0(i,j+1))) else h_face = h0(i,j+1) endif endif endif - vh_ice(i,J) = time_step * G%dxCv(i,J) * v_face * h_face + vh_ice(i,J) = (time_step * G%dxCv(i,J)) * (v_face * h_face) else vh_ice(i,J) = 0.0 endif @@ -1902,7 +1930,11 @@ subroutine shelf_advance_front(CS, ISS, G, hmask, uh_ice, vh_ice) integer :: iter_flag real :: h_reference ! A reference thicknesss based on neighboring cells [Z ~> m] + real :: h_reference_ew !contribution to reference thickness from east + west cells [Z ~> m] + real :: h_reference_ns !contribution to reference thickness from north + south cells [Z ~> m] real :: tot_flux ! The total ice mass flux [Z L2 ~> m3] + real :: tot_flux_ew ! The contribution to total ice mass flux from east + west cells [Z L2 ~> m3] + real :: tot_flux_ns ! The contribution to total ice mass flux from north + south cells [Z L2 ~> m3] real :: partial_vol ! The volume covered by ice shelf [Z L2 ~> m3] real :: dxdyh ! Cell area [L2 ~> m2] character(len=160) :: mesg ! The text of an error message @@ -1953,15 +1985,17 @@ subroutine shelf_advance_front(CS, ISS, G, hmask, uh_ice, vh_ice) ((i+i_off) >= 1)) then ! first get reference thickness by averaging over cells that are fluxing into this cell n_flux = 0 - h_reference = 0.0 - tot_flux = 0.0 + h_reference_ew = 0.0 + h_reference_ns = 0.0 + tot_flux_ew = 0.0 + tot_flux_ns = 0.0 do k=1,2 if (flux_enter(i,j,k) > 0) then n_flux = n_flux + 1 - h_reference = h_reference + flux_enter(i,j,k) * ISS%h_shelf(i+2*k-3,j) + h_reference_ew = h_reference_ew + flux_enter(i,j,k) * ISS%h_shelf(i+2*k-3,j) !h_reference = h_reference + ISS%h_shelf(i+2*k-3,j) - tot_flux = tot_flux + flux_enter(i,j,k) + tot_flux_ew = tot_flux_ew + flux_enter(i,j,k) flux_enter(i,j,k) = 0.0 endif enddo @@ -1969,13 +2003,16 @@ subroutine shelf_advance_front(CS, ISS, G, hmask, uh_ice, vh_ice) do k=1,2 if (flux_enter(i,j,k+2) > 0) then n_flux = n_flux + 1 - h_reference = h_reference + flux_enter(i,j,k+2) * ISS%h_shelf(i,j+2*k-3) + h_reference_ns = h_reference_ns + flux_enter(i,j,k+2) * ISS%h_shelf(i,j+2*k-3) !h_reference = h_reference + ISS%h_shelf(i,j+2*k-3) - tot_flux = tot_flux + flux_enter(i,j,k+2) + tot_flux_ns = tot_flux_ns + flux_enter(i,j,k+2) flux_enter(i,j,k+2) = 0.0 endif enddo + h_reference = h_reference_ew + h_reference_ns + tot_flux = tot_flux_ew + tot_flux_ns + if (n_flux > 0) then dxdyh = G%areaT(i,j) h_reference = h_reference / tot_flux @@ -1983,7 +2020,7 @@ subroutine shelf_advance_front(CS, ISS, G, hmask, uh_ice, vh_ice) partial_vol = ISS%h_shelf(i,j) * ISS%area_shelf_h(i,j) + tot_flux if ((partial_vol / G%areaT(i,j)) == h_reference) then ! cell is exactly covered, no overflow - if (ISS%hmask(i,j).ne.3) ISS%hmask(i,j) = 1 + if (ISS%hmask(i,j)/=3) ISS%hmask(i,j) = 1 ISS%h_shelf(i,j) = h_reference ISS%area_shelf_h(i,j) = G%areaT(i,j) elseif ((partial_vol / G%areaT(i,j)) < h_reference) then @@ -1993,7 +2030,7 @@ subroutine shelf_advance_front(CS, ISS, G, hmask, uh_ice, vh_ice) ISS%h_shelf(i,j) = h_reference else - if (ISS%hmask(i,j).ne.3) ISS%hmask(i,j) = 1 + if (ISS%hmask(i,j)/=3) ISS%hmask(i,j) = 1 ISS%area_shelf_h(i,j) = G%areaT(i,j) !h_temp(i,j) = h_reference partial_vol = partial_vol - h_reference * G%areaT(i,j) @@ -2131,7 +2168,7 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) ! (it is assumed that base_ice = bed + OD) real, dimension(SIZE(OD,1),SIZE(OD,2)) :: S ! surface elevation [Z ~> m]. - + real, dimension(SZDI_(G),SZDJ_(G)) :: sx_e, sy_e !element contributions to driving stress real :: rho, rhow, rhoi_rhow ! Ice and ocean densities [R ~> kg m-3] real :: sx, sy ! Ice shelf top slopes [Z L-1 ~> nondim] real :: neumann_val ! [R Z L2 T-2 ~> kg s-2] @@ -2143,15 +2180,14 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) integer :: i_off, j_off 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 +! 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 +! iegq = G%iegB ; jegq = G%jegB ! gisc = G%domain%nihalo+1 ; gjsc = 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 ; gjec = G%domain%njglobal - is = iscq - 1; js = jscq - 1 +! is = iscq - 1; js = jscq - 1 i_off = G%idg_offset ; j_off = G%jdg_offset @@ -2182,6 +2218,8 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) call pass_var(S, G%domain) + sx_e(:,:)=0.0; sy_e(:,:)=0.0 + do j=jsc-1,jec+1 do i=isc-1,iec+1 cnt = 0 @@ -2210,14 +2248,14 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) else ! interior if (ISS%hmask(i+1,j) == 1 .or. ISS%hmask(i+1,j) == 3) then cnt = cnt+1 - Dx =dxh+ G%dxT(i+1,j) + Dx = dxh + G%dxT(i+1,j) sx = S(i+1,j) else sx = S(i,j) endif if (ISS%hmask(i-1,j) == 1 .or. ISS%hmask(i-1,j) == 3) then cnt = cnt+1 - Dx =dxh+ G%dxT(i-1,j) + Dx = dxh + G%dxT(i-1,j) sx = sx - S(i-1,j) else sx = sx - S(i,j) @@ -2225,7 +2263,7 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) if (cnt == 0) then sx = 0 else - sx = sx / ( Dx) + sx = sx / Dx endif endif @@ -2247,54 +2285,36 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) else ! interior if (ISS%hmask(i,j+1) == 1 .or. ISS%hmask(i,j+1) == 3) then cnt = cnt+1 - Dy =dyh+ G%dyT(i,j+1) + Dy = dyh + G%dyT(i,j+1) sy = S(i,j+1) else sy = S(i,j) endif if (ISS%hmask(i,j-1) == 1 .or. ISS%hmask(i,j-1) == 3) then cnt = cnt+1 + Dy = dyh + G%dyT(i,j-1) sy = sy - S(i,j-1) - Dy =dyh+ G%dyT(i,j-1) else sy = sy - S(i,j) endif if (cnt == 0) then sy = 0 else - sy = sy / (Dy) + sy = sy / Dy endif endif - ! SW vertex - !if (ISS%hmask(I-1,J-1) == 1) 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 - ! SE vertex - !if (ISS%hmask(I,J-1) == 1) 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 - ! NW vertex - !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 - ! 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 + sx_e(i,j) = (-.25 * G%areaT(i,j)) * ((rho * grav) * (ISS%h_shelf(i,j) * sx)) + sy_e(i,j) = (-.25 * G%areaT(i,j)) * ((rho * grav) * (ISS%h_shelf(i,j) * sy)) !Stress (Neumann) boundary conditions if (CS%ground_frac(i,j) == 1) then - neumann_val = (.5 * grav * (rho * ISS%h_shelf(i,j)**2 - rhow * CS%bed_elev(i,j)**2)) + neumann_val = ((.5 * grav) * (rho * ISS%h_shelf(i,j)**2 - rhow * CS%bed_elev(i,j)**2)) else - neumann_val = (.5 * grav * (1-rho/rhow) * rho * ISS%h_shelf(i,j)**2) + neumann_val = (.5 * grav) * ((1-rho/rhow) * (rho * ISS%h_shelf(i,j)**2)) endif if ((CS%u_face_mask_bdry(I-1,j) == 2) .OR. & - ((ISS%hmask(i-1,j) == 0 .OR. ISS%hmask(i-1,j) == 2) .AND. (i+i_off .ne. gisc))) then + ((ISS%hmask(i-1,j) == 0 .OR. ISS%hmask(i-1,j) == 2) .AND. (i+i_off /= gisc))) then ! left face of the cell is at a stress boundary ! the depth-integrated longitudinal stress is equal to the difference of depth-integrated ! pressure on either side of the face @@ -2309,30 +2329,33 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) endif if ((CS%u_face_mask_bdry(I,j) == 2) .OR. & - ((ISS%hmask(i+1,j) == 0 .OR. ISS%hmask(i+1,j) == 2) .and. (i+i_off .ne. giec))) then + ((ISS%hmask(i+1,j) == 0 .OR. ISS%hmask(i+1,j) == 2) .and. (i+i_off /= giec))) then ! east face of the cell is at a stress boundary taudx(I,J-1) = taudx(I,J-1) + .5 * dyh * neumann_val taudx(I,J) = taudx(I,J) + .5 * dyh * neumann_val endif if ((CS%v_face_mask_bdry(i,J-1) == 2) .OR. & - ((ISS%hmask(i,j-1) == 0 .OR. ISS%hmask(i,j-1) == 2) .and. (j+j_off .ne. gjsc))) then + ((ISS%hmask(i,j-1) == 0 .OR. ISS%hmask(i,j-1) == 2) .and. (j+j_off /= gjsc))) then ! south face of the cell is at a stress boundary taudy(I-1,J-1) = taudy(I-1,J-1) - .5 * dxh * neumann_val taudy(I,J-1) = taudy(I,J-1) - .5 * dxh * neumann_val endif if ((CS%v_face_mask_bdry(i,J) == 2) .OR. & - ((ISS%hmask(i,j+1) == 0 .OR. ISS%hmask(i,j+1) == 2) .and. (j+j_off .ne. gjec))) then + ((ISS%hmask(i,j+1) == 0 .OR. ISS%hmask(i,j+1) == 2) .and. (j+j_off /= gjec))) then ! north face of the cell is at a stress boundary taudy(I-1,J) = taudy(I-1,J) + .5 * dxh * neumann_val taudy(I,J) = taudy(I,J) + .5 * dxh * neumann_val endif - endif enddo enddo + do J=jsc-2,jec+1; do I=isc-2,iec+1 + taudx(I,J) = taudx(I,J) + ((sx_e(i,j)+sx_e(i+1,j+1)) + (sx_e(i+1,j)+sx_e(i,j+1))) + taudy(I,J) = taudy(I,J) + ((sy_e(i,j)+sy_e(i+1,j+1)) + (sy_e(i+1,j)+sy_e(i,j+1))) + enddo; enddo end subroutine calc_shelf_driving_stress ! Not used? Seems to be only set up to work for a specific test case with u_face_mask==3 @@ -2369,7 +2392,7 @@ subroutine init_boundary_values(CS, G, time, hmask, input_flux, input_thick, new do i=isd,ied if (hmask(i,j) == 3) then - CS%thickness_bdry_val(i,j) = input_thick + CS%h_bdry_val(i,j) = input_thick endif if ((hmask(i,j) == 0) .or. (hmask(i,j) == 1) .or. (hmask(i,j) == 2)) then @@ -2436,7 +2459,7 @@ subroutine CG_action(CS, uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask, 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)), & + real, dimension(SZDI_(G),SZDJ_(G),CS%visc_qps), & intent(in) :: ice_visc !< A field related to the ice viscosity from Glen's !! flow law [R L4 Z T-1 ~> kg m2 s-1]. The exact form !! and units depend on the basal law exponent. @@ -2479,88 +2502,124 @@ subroutine CG_action(CS, uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask, real :: ux, uy, vx, vy ! Components of velocity shears or divergence [T-1 ~> s-1] real :: uq, vq ! Interpolated velocities [L T-1 ~> m s-1] - integer :: iq, jq, iphi, jphi, i, j, ilq, jlq, Itgt, Jtgt + integer :: iq, jq, iphi, jphi, i, j, ilq, jlq, Itgt, Jtgt, qp, qpv + logical :: visc_qp4 real, dimension(2) :: xquad real, dimension(2,2) :: Ucell, Vcell, Hcell, Usub, Vsub - real :: Ee + real, dimension(2,2,4) :: uret_qp, vret_qp + real, dimension(SZDIB_(G),SZDJB_(G),4) :: uret_b, vret_b xquad(1) = .5 * (1-sqrt(1./3)) ; xquad(2) = .5 * (1+sqrt(1./3)) - Ee=1.0 + if (CS%visc_qps == 4) then + visc_qp4=.true. + else + visc_qp4=.false. + qpv = 1 + endif + + uret(:,:) = 0.0; vret(:,:)=0.0 + uret_b(:,:,:)=0.0 ; vret_b(:,:,:)=0.0 do j=js,je ; do i=is,ie ; if (hmask(i,j) == 1 .or. hmask(i,j)==3) then + uret_qp(:,:,:)=0.0; vret_qp(:,:,:)=0.0 + do iq=1,2 ; do jq=1,2 - uq = u_shlf(I-1,J-1) * (xquad(3-iq) * xquad(3-jq)) + & - u_shlf(I,J-1) * (xquad(iq) * xquad(3-jq)) + & - u_shlf(I-1,J) * (xquad(3-iq) * xquad(jq)) + & - u_shlf(I,J) * (xquad(iq) * xquad(jq)) - - vq = v_shlf(I-1,J-1) * (xquad(3-iq) * xquad(3-jq)) + & - v_shlf(I,J-1) * (xquad(iq) * xquad(3-jq)) + & - v_shlf(I-1,J) * (xquad(3-iq) * xquad(jq)) + & - v_shlf(I,J) * (xquad(iq) * xquad(jq)) - - ux = u_shlf(I-1,J-1) * Phi(1,2*(jq-1)+iq,i,j) + & - u_shlf(I,J-1) * Phi(3,2*(jq-1)+iq,i,j) + & - u_shlf(I-1,J) * Phi(5,2*(jq-1)+iq,i,j) + & - u_shlf(I,J) * Phi(7,2*(jq-1)+iq,i,j) - - vx = v_shlf(I-1,J-1) * Phi(1,2*(jq-1)+iq,i,j) + & - v_shlf(I,J-1) * Phi(3,2*(jq-1)+iq,i,j) + & - v_shlf(I-1,J) * Phi(5,2*(jq-1)+iq,i,j) + & - v_shlf(I,J) * Phi(7,2*(jq-1)+iq,i,j) - - uy = u_shlf(I-1,J-1) * Phi(2,2*(jq-1)+iq,i,j) + & - u_shlf(I,J-1) * Phi(4,2*(jq-1)+iq,i,j) + & - u_shlf(I-1,J) * Phi(6,2*(jq-1)+iq,i,j) + & - u_shlf(I,J) * Phi(8,2*(jq-1)+iq,i,j) - - vy = v_shlf(I-1,J-1) * Phi(2,2*(jq-1)+iq,i,j) + & - v_shlf(I,J-1) * Phi(4,2*(jq-1)+iq,i,j) + & - 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) - - if (trim(CS%ice_viscosity_compute) == "MODEL_QUADRATURE") Ee = CS%Ee(i,j,2*(jq-1)+iq) - - 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 * Ee * 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)) - if (vmask(Itgt,Jtgt) == 1) vret(Itgt,Jtgt) = vret(Itgt,Jtgt) + 0.25 * Ee * ice_visc(i,j) * & - ((uy+vx) * Phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq,i,j) + & - (4*vy+2*ux) * Phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq,i,j)) + qp = 2*(jq-1)+iq !current quad point + + uq = (u_shlf(I-1,J-1) * (xquad(3-iq) * xquad(3-jq)) + & + u_shlf(I,J) * (xquad(iq) * xquad(jq))) + & + (u_shlf(I,J-1) * (xquad(iq) * xquad(3-jq)) + & + u_shlf(I-1,J) * (xquad(3-iq) * xquad(jq))) + + vq = (v_shlf(I-1,J-1) * (xquad(3-iq) * xquad(3-jq)) + & + v_shlf(I,J) * (xquad(iq) * xquad(jq))) + & + (v_shlf(I,J-1) * (xquad(iq) * xquad(3-jq)) + & + v_shlf(I-1,J) * (xquad(3-iq) * xquad(jq))) + + ux = (u_shlf(I-1,J-1) * Phi(1,qp,i,j) + & + u_shlf(I,J) * Phi(7,qp,i,j)) + & + (u_shlf(I,J-1) * Phi(3,qp,i,j) + & + u_shlf(I-1,J) * Phi(5,qp,i,j)) + + vx = (v_shlf(I-1,J-1) * Phi(1,qp,i,j) + & + v_shlf(I,J) * Phi(7,qp,i,j)) + & + (v_shlf(I,J-1) * Phi(3,qp,i,j) + & + v_shlf(I-1,J) * Phi(5,qp,i,j)) + + uy = (u_shlf(I-1,J-1) * Phi(2,qp,i,j) + & + u_shlf(I,J) * Phi(8,qp,i,j)) + & + (u_shlf(I,J-1) * Phi(4,qp,i,j) + & + u_shlf(I-1,J) * Phi(6,qp,i,j)) + + vy = (v_shlf(I-1,J-1) * Phi(2,qp,i,j) + & + v_shlf(I,J) * Phi(8,qp,i,j)) + & + (v_shlf(I,J-1) * Phi(4,qp,i,j) + & + v_shlf(I-1,J) * Phi(6,qp,i,j)) + + if (visc_qp4) qpv = qp !current quad point for viscosity + + do jphi=1,2 ; Jtgt = J-2+jphi ; do iphi=1,2 ; Itgt = I-2+iphi + if (umask(Itgt,Jtgt) == 1) uret_qp(iphi,jphi,qp) = ice_visc(i,j,qpv) * & + ((4*ux+2*vy) * Phi(2*(2*(jphi-1)+iphi)-1,qp,i,j) + & + (uy+vx) * Phi(2*(2*(jphi-1)+iphi),qp,i,j)) + if (vmask(Itgt,Jtgt) == 1) vret_qp(iphi,jphi,qp) = ice_visc(i,j,qpv) * & + ((uy+vx) * Phi(2*(2*(jphi-1)+iphi)-1,qp,i,j) + & + (4*vy+2*ux) * Phi(2*(2*(jphi-1)+iphi),qp,i,j)) if (float_cond(i,j) == 0) then ilq = 1 ; if (iq == iphi) ilq = 2 jlq = 1 ; if (jq == jphi) jlq = 2 - if (umask(Itgt,Jtgt) == 1) uret(Itgt,Jtgt) = uret(Itgt,Jtgt) + & - 0.25 * basal_trac(i,j) * uq * (xquad(ilq) * xquad(jlq)) - if (vmask(Itgt,Jtgt) == 1) vret(Itgt,Jtgt) = vret(Itgt,Jtgt) + & - 0.25 * basal_trac(i,j) * vq * (xquad(ilq) * xquad(jlq)) + if (umask(Itgt,Jtgt) == 1) uret_qp(iphi,jphi,qp) = uret_qp(iphi,jphi,qp) + & + (basal_trac(i,j) * uq) * (xquad(ilq) * xquad(jlq)) + if (vmask(Itgt,Jtgt) == 1) vret_qp(iphi,jphi,qp) = vret_qp(iphi,jphi,qp) + & + (basal_trac(i,j) * vq) * (xquad(ilq) * xquad(jlq)) endif enddo ; enddo enddo ; enddo + !element contribution to SW node (node 1, which sees the current element as element 4) + uret_b(I-1,J-1,4) = 0.25*((uret_qp(1,1,1)+uret_qp(1,1,4))+(uret_qp(1,1,2)+uret_qp(1,1,3))) + vret_b(I-1,J-1,4) = 0.25*((vret_qp(1,1,1)+vret_qp(1,1,4))+(vret_qp(1,1,2)+vret_qp(1,1,3))) + + !element contribution to NW node (node 3, which sees the current element as element 2) + uret_b(I-1,J ,2) = 0.25*((uret_qp(1,2,1)+uret_qp(1,2,4))+(uret_qp(1,2,2)+uret_qp(1,2,3))) + vret_b(I-1,J ,2) = 0.25*((vret_qp(1,2,1)+vret_qp(1,2,4))+(vret_qp(1,2,2)+vret_qp(1,2,3))) + + !element contribution to SE node (node 2, which sees the current element as element 3) + uret_b(I ,J-1,3) = 0.25*((uret_qp(2,1,1)+uret_qp(2,1,4))+(uret_qp(2,1,2)+uret_qp(2,1,3))) + vret_b(I ,J-1,3) = 0.25*((vret_qp(2,1,1)+vret_qp(2,1,4))+(vret_qp(2,1,2)+vret_qp(2,1,3))) + + !element contribution to NE node (node 4, which sees the current element as element 1) + uret_b(I ,J ,1) = 0.25*((uret_qp(2,2,1)+uret_qp(2,2,4))+(uret_qp(2,2,2)+uret_qp(2,2,3))) + vret_b(I ,J ,1) = 0.25*((vret_qp(2,2,1)+vret_qp(2,2,4))+(vret_qp(2,2,2)+vret_qp(2,2,3))) + if (float_cond(i,j) == 1) then Ucell(:,:) = u_shlf(I-1:I,J-1:J) ; Vcell(:,:) = v_shlf(I-1:I,J-1:J) - Hcell(:,:) = H_node(i-1:i,j-1:j) - call CG_action_subgrid_basal(Phisub, Hcell, Ucell, Vcell, bathyT(i,j), dens_ratio, Usub, Vsub) - - if (umask(I-1,J-1) == 1) uret(I-1,J-1) = uret(I-1,J-1) + Usub(1,1) * basal_trac(i,j) - if (umask(I-1,J) == 1) uret(I-1,J) = uret(I-1,J) + Usub(1,2) * basal_trac(i,j) - if (umask(I,J-1) == 1) uret(I,J-1) = uret(I,J-1) + Usub(2,1) * basal_trac(i,j) - if (umask(I,J) == 1) uret(I,J) = uret(I,J) + Usub(2,2) * basal_trac(i,j) - - if (vmask(I-1,J-1) == 1) vret(I-1,J-1) = vret(I-1,J-1) + Vsub(1,1) * basal_trac(i,j) - if (vmask(I-1,J) == 1) vret(I-1,J) = vret(I-1,J) + Vsub(1,2) * basal_trac(i,j) - if (vmask(I,J-1) == 1) vret(I,J-1) = vret(I,J-1) + Vsub(2,1) * basal_trac(i,j) - if (vmask(I,J) == 1) vret(I,J) = vret(I,J) + Vsub(2,2) * basal_trac(i,j) - endif + Hcell(:,:) = H_node(I-1:I,J-1:J) + + call CG_action_subgrid_basal(Phisub, Hcell, Ucell, Vcell, & + bathyT(i,j), dens_ratio, Usub, Vsub) + if (umask(I-1,J-1) == 1) uret_b(I-1,J-1,4) = uret_b(I-1,J-1,4) + (Usub(1,1) * basal_trac(i,j)) + if (umask(I-1,J ) == 1) uret_b(I-1,J ,2) = uret_b(I-1,J ,2) + (Usub(1,2) * basal_trac(i,j)) + if (umask(I ,J-1) == 1) uret_b(I ,J-1,3) = uret_b(I ,J-1,3) + (Usub(2,1) * basal_trac(i,j)) + if (umask(I ,J ) == 1) uret_b(I ,J ,1) = uret_b(I ,J ,1) + (Usub(2,2) * basal_trac(i,j)) + + if (vmask(I-1,J-1) == 1) vret_b(I-1,J-1,4) = vret_b(I-1,J-1,4) + (Vsub(1,1) * basal_trac(i,j)) + if (vmask(I-1,J ) == 1) vret_b(I-1,J ,2) = vret_b(I-1,J ,2) + (Vsub(1,2) * basal_trac(i,j)) + if (vmask(I ,J-1) == 1) vret_b(I ,J-1,3) = vret_b(I ,J-1,3) + (Vsub(2,1) * basal_trac(i,j)) + if (vmask(I ,J ) == 1) vret_b(I ,J ,1) = vret_b(I ,J ,1) + (Vsub(2,2) * basal_trac(i,j)) + endif endif ; enddo ; enddo + do J=js-1,je ; do I=is-1,ie + uret(I,J) = (uret_b(I,J,1)+uret_b(I,J,4)) + (uret_b(I,J,2)+uret_b(I,J,3)) + vret(I,J) = (vret_b(I,J,1)+vret_b(I,J,4)) + (vret_b(I,J,2)+vret_b(I,J,3)) + enddo; enddo + end subroutine CG_action subroutine CG_action_subgrid_basal(Phisub, H, U, V, bathyT, dens_ratio, Ucontr, Vcontr) @@ -2579,45 +2638,105 @@ subroutine CG_action_subgrid_basal(Phisub, H, U, V, bathyT, dens_ratio, Ucontr, real, dimension(2,2), intent(out) :: Vcontr !< The areal average of v-velocities where the ice shelf !! is grounded, or 0 where it is floating [L T-1 ~> m s-1]. + real, dimension(SIZE(Phisub,3),SIZE(Phisub,3),2,2) :: Ucontr_sub, Vcontr_sub ! The contributions to Ucontr and Vcontr + !! at each sub-cell + real, dimension(2,2,SIZE(Phisub,3),SIZE(Phisub,3)) :: uloc_arr !The local sub-cell u-velocity [L T-1 ~> m s-1] + real, dimension(2,2,SIZE(Phisub,3),SIZE(Phisub,3)) :: vloc_arr !The local sub-cell v-velocity [L T-1 ~> m s-1] + real, dimension(2,2) :: Ucontr_q, Vcontr_q !Contributions to a node from each quadrature point in a sub-grid cell real :: subarea ! The fractional sub-cell area [nondim] real :: hloc ! The local sub-cell ice thickness [Z ~> m] integer :: nsub, i, j, qx, qy, m, n - nsub = size(Phisub,1) + nsub = size(Phisub,3) subarea = 1.0 / (nsub**2) + uloc_arr(:,:,:,:) = 0.0; vloc_arr(:,:,:,:)=0.0 + + do j=1,nsub ; do i=1,nsub; do qy=1,2 ; do qx=1,2 + hloc = (Phisub(qx,qy,i,j,1,1)*H(1,1) + Phisub(qx,qy,i,j,2,2)*H(2,2)) + & + (Phisub(qx,qy,i,j,1,2)*H(1,2) + Phisub(qx,qy,i,j,2,1)*H(2,1)) + if (dens_ratio * hloc - bathyT > 0) then + uloc_arr(qx,qy,i,j) = ((Phisub(qx,qy,i,j,1,1) * U(1,1) + Phisub(qx,qy,i,j,2,2) * U(2,2)) + & + (Phisub(qx,qy,i,j,1,2) * U(1,2) + Phisub(qx,qy,i,j,2,1) * U(2,1))) + vloc_arr(qx,qy,i,j) = ((Phisub(qx,qy,i,j,1,1) * V(1,1) + Phisub(qx,qy,i,j,2,2) * V(2,2)) + & + (Phisub(qx,qy,i,j,1,2) * V(1,2) + Phisub(qx,qy,i,j,2,1) * V(2,1))) + endif + enddo; enddo ; enddo ; enddo + + do n=1,2 ; do m=1,2 ; do j=1,nsub ; do i=1,nsub + do qy=1,2 ; do qx=1,2 + !calculate quadrature point contributions for the sub-cell, to each node + Ucontr_q(qx,qy) = Phisub(qx,qy,i,j,m,n) * uloc_arr(qx,qy,i,j) + Vcontr_q(qx,qy) = Phisub(qx,qy,i,j,m,n) * vloc_arr(qx,qy,i,j) + enddo; enddo + + !calculate sub-cell contribution to each node by summing up quadrature point contributions from the sub-cell + Ucontr_sub(i,j,m,n) = (subarea * 0.25) * ((Ucontr_q(1,1) + Ucontr_q(2,2)) + (Ucontr_q(1,2)+Ucontr_q(2,1))) + Vcontr_sub(i,j,m,n) = (subarea * 0.25) * ((Vcontr_q(1,1) + Vcontr_q(2,2)) + (Vcontr_q(1,2)+Vcontr_q(2,1))) + enddo; enddo ; enddo ; enddo + + !sum up the sub-cell contributions to each node do n=1,2 ; do m=1,2 - Ucontr(m,n) = 0.0 ; Vcontr(m,n) = 0.0 - do qy=1,2 ; do qx=1,2 ; do j=1,nsub ; do i=1,nsub - hloc = (Phisub(i,j,1,1,qx,qy)*H(1,1) + Phisub(i,j,2,2,qx,qy)*H(2,2)) + & - (Phisub(i,j,1,2,qx,qy)*H(1,2) + Phisub(i,j,2,1,qx,qy)*H(2,1)) - if (dens_ratio * hloc - bathyT > 0) then - Ucontr(m,n) = Ucontr(m,n) + subarea * 0.25 * Phisub(i,j,m,n,qx,qy) * & - ((Phisub(i,j,1,1,qx,qy) * U(1,1) + Phisub(i,j,2,2,qx,qy) * U(2,2)) + & - (Phisub(i,j,1,2,qx,qy) * U(1,2) + Phisub(i,j,2,1,qx,qy) * U(2,1))) - Vcontr(m,n) = Vcontr(m,n) + subarea * 0.25 * Phisub(i,j,m,n,qx,qy) * & - ((Phisub(i,j,1,1,qx,qy) * V(1,1) + Phisub(i,j,2,2,qx,qy) * V(2,2)) + & - (Phisub(i,j,1,2,qx,qy) * V(1,2) + Phisub(i,j,2,1,qx,qy) * V(2,1))) - endif - enddo ; enddo ; enddo ; enddo + call sum_square_matrix(Ucontr(m,n),Ucontr_sub(:,:,m,n),nsub) + call sum_square_matrix(Vcontr(m,n),Vcontr_sub(:,:,m,n),nsub) enddo ; enddo end subroutine CG_action_subgrid_basal + +!! Returns the sum of the elements in a square matrix. This sum is bitwise identical even if the matrices are rotated. +subroutine sum_square_matrix(sum_out, mat_in, n) + integer, intent(in) :: n !< The length and width of each matrix in mat_in + real, dimension(n,n), intent(in) :: mat_in !< The n x n matrix whose elements will be summed + real, intent(out) :: sum_out !< The sum of the elements of matrix mat_in + integer :: s0,e0,s1,e1 + + sum_out=0.0 + + s0=1; e0=n + + !start by summing elements on outer edges of matrix + do while (s0 returns the diagonal entries of the matrix for a Jacobi preconditioning subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, hmask, dens_ratio, & - Phisub, u_diagonal, v_diagonal) + Phi, Phisub, u_diagonal, v_diagonal) type(ice_shelf_dyn_CS), intent(in) :: CS !< A pointer to the ice shelf control structure 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(SZDI_(G),SZDJ_(G)), & - intent(in) :: float_cond !< If GL_regularize=true, an array indicating where the ice - !! shelf is floating: 0 if floating, 1 if not + intent(in) :: float_cond !< If GL_regularize=true, indicates cells containing + !! the grounding line (float_cond=1) or not (float_cond=0) real, dimension(SZDIB_(G),SZDJB_(G)), & intent(in) :: H_node !< The ice shelf thickness at nodal !! (corner) points [Z ~> m]. - real, dimension(SZDI_(G),SZDJ_(G)), & + real, dimension(SZDI_(G),SZDJ_(G),CS%visc_qps), & intent(in) :: ice_visc !< A field related to the ice viscosity from Glen's !! flow law [R L4 Z T-1 ~> kg m2 s-1]. The exact form !! and units depend on the basal law exponent. @@ -2630,6 +2749,9 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, !! partly or fully covered by an ice-shelf real, intent(in) :: dens_ratio !< The density of ice divided by the density !! of seawater [nondim] + 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(:,:,:,:,:,:), intent(in) :: Phisub !< Quadrature structure weights at subgridscale !! locations for finite element calculations [nondim] real, dimension(SZDIB_(G),SZDJB_(G)), & @@ -2644,89 +2766,122 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, real :: ux, uy, vx, vy ! Interpolated weight gradients [L-1 ~> m-1] real :: uq, vq - real, dimension(8,4) :: Phi ! Weight gradients [L-1 ~> m-1] real, dimension(2) :: xquad real, dimension(2,2) :: Hcell, sub_ground - integer :: i, j, isc, jsc, iec, jec, iphi, jphi, iq, jq, ilq, jlq, Itgt, Jtgt - real :: Ee + real, dimension(2,2,4) :: u_diag_qp, v_diag_qp + real, dimension(SZDIB_(G),SZDJB_(G),4) :: u_diag_b, v_diag_b + logical :: visc_qp4 + integer :: i, j, isc, jsc, iec, jec, iphi, jphi, iq, jq, ilq, jlq, Itgt, Jtgt, qp, qpv isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec xquad(1) = .5 * (1-sqrt(1./3)) ; xquad(2) = .5 * (1+sqrt(1./3)) - Ee=1.0 + if (CS%visc_qps == 4) then + visc_qp4=.true. + else + visc_qp4=.false. + qpv = 1 + endif - do j=jsc-1,jec+1 ; do i=isc-1,iec+1 ; if (hmask(i,j) == 1 .or. hmask(i,j)==3) then + u_diag_b(:,:,:)=0.0 + v_diag_b(:,:,:)=0.0 - call bilinear_shape_fn_grid(G, i, j, Phi) + do j=jsc-1,jec+1 ; do i=isc-1,iec+1 ; if (hmask(i,j) == 1 .or. hmask(i,j)==3) then ! 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 + u_diag_qp(:,:,:)=0.0; v_diag_qp(:,:,:)=0.0 + do iq=1,2 ; do jq=1,2 - if (trim(CS%ice_viscosity_compute) == "MODEL_QUADRATURE") Ee = CS%Ee(i,j,2*(jq-1)+iq) - do iphi=1,2 ; do jphi=1,2 + qp = 2*(jq-1)+iq !current quad point + if (visc_qp4) qpv = qp !current quad point for viscosity + + do jphi=1,2 ; Jtgt = J-2+jphi ; do iphi=1,2 ; Itgt = I-2+iphi - Itgt = I-2+iphi ; Jtgt = J-2+jphi ilq = 1 ; if (iq == iphi) ilq = 2 jlq = 1 ; if (jq == jphi) jlq = 2 if (CS%umask(Itgt,Jtgt) == 1) then - ux = Phi(2*(2*(jphi-1)+iphi)-1, 2*(jq-1)+iq) - uy = Phi(2*(2*(jphi-1)+iphi), 2*(jq-1)+iq) + ux = Phi(2*(2*(jphi-1)+iphi)-1,qp,i,j) + uy = Phi(2*(2*(jphi-1)+iphi),qp,i,j) vx = 0. vy = 0. - u_diagonal(Itgt,Jtgt) = u_diagonal(Itgt,Jtgt) + & - 0.25 * Ee * ice_visc(i,j) * ((4*ux+2*vy) * Phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + & - (uy+vx) * Phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq)) + u_diag_qp(iphi,jphi,qp) = & + ice_visc(i,j,qpv) * ((4*ux+2*vy) * Phi(2*(2*(jphi-1)+iphi)-1,qp,i,j) + & + (uy+vx) * Phi(2*(2*(jphi-1)+iphi),qp,i,j)) if (float_cond(i,j) == 0) then uq = xquad(ilq) * xquad(jlq) - u_diagonal(Itgt,Jtgt) = u_diagonal(Itgt,Jtgt) + & - 0.25 * basal_trac(i,j) * uq * (xquad(ilq) * xquad(jlq)) + u_diag_qp(iphi,jphi,qp) = u_diag_qp(iphi,jphi,qp) + & + (basal_trac(i,j) * uq) * (xquad(ilq) * xquad(jlq)) endif endif if (CS%vmask(Itgt,Jtgt) == 1) then - vx = Phi(2*(2*(jphi-1)+iphi)-1, 2*(jq-1)+iq) - vy = Phi(2*(2*(jphi-1)+iphi), 2*(jq-1)+iq) + vx = Phi(2*(2*(jphi-1)+iphi)-1,qp,i,j) + vy = Phi(2*(2*(jphi-1)+iphi),qp,i,j) ux = 0. uy = 0. - v_diagonal(Itgt,Jtgt) = v_diagonal(Itgt,Jtgt) + & - 0.25 * Ee * ice_visc(i,j) * ((uy+vx) * Phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + & - (4*vy+2*ux) * Phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq)) + v_diag_qp(iphi,jphi,qp) = & + ice_visc(i,j,qpv) * ((uy+vx) * Phi(2*(2*(jphi-1)+iphi)-1,qp,i,j) + & + (4*vy+2*ux) * Phi(2*(2*(jphi-1)+iphi),qp,i,j)) if (float_cond(i,j) == 0) then vq = xquad(ilq) * xquad(jlq) - v_diagonal(Itgt,Jtgt) = v_diagonal(Itgt,Jtgt) + & - 0.25 * basal_trac(i,j) * vq * (xquad(ilq) * xquad(jlq)) + v_diag_qp(iphi,jphi,qp) = v_diag_qp(iphi,jphi,qp) + & + (basal_trac(i,j) * vq) * (xquad(ilq) * xquad(jlq)) endif endif enddo ; enddo enddo ; enddo + !element contribution to SW node (node 1, which sees the current element as element 4) + u_diag_b(I-1,J-1,4) = 0.25*((u_diag_qp(1,1,1)+u_diag_qp(1,1,4))+(u_diag_qp(1,1,2)+u_diag_qp(1,1,3))) + v_diag_b(I-1,J-1,4) = 0.25*((v_diag_qp(1,1,1)+v_diag_qp(1,1,4))+(v_diag_qp(1,1,2)+v_diag_qp(1,1,3))) + + !element contribution to NW node (node 3, which sees the current element as element 2) + u_diag_b(I-1,J ,2) = 0.25*((u_diag_qp(1,2,1)+u_diag_qp(1,2,4))+(u_diag_qp(1,2,2)+u_diag_qp(1,2,3))) + v_diag_b(I-1,J ,2) = 0.25*((v_diag_qp(1,2,1)+v_diag_qp(1,2,4))+(v_diag_qp(1,2,2)+v_diag_qp(1,2,3))) + + !element contribution to SE node (node 2, which sees the current element as element 3) + u_diag_b(I ,J-1,3) = 0.25*((u_diag_qp(2,1,1)+u_diag_qp(2,1,4))+(u_diag_qp(2,1,2)+u_diag_qp(2,1,3))) + v_diag_b(I ,J-1,3) = 0.25*((v_diag_qp(2,1,1)+v_diag_qp(2,1,4))+(v_diag_qp(2,1,2)+v_diag_qp(2,1,3))) + + !element contribution to NE node (node 4, which sees the current element as element 1) + u_diag_b(I ,J ,1) = 0.25*((u_diag_qp(2,2,1)+u_diag_qp(2,2,4))+(u_diag_qp(2,2,2)+u_diag_qp(2,2,3))) + v_diag_b(I ,J ,1) = 0.25*((v_diag_qp(2,2,1)+v_diag_qp(2,2,4))+(v_diag_qp(2,2,2)+v_diag_qp(2,2,3))) + if (float_cond(i,j) == 1) then Hcell(:,:) = H_node(i-1:i,j-1:j) call CG_diagonal_subgrid_basal(Phisub, Hcell, CS%bed_elev(i,j), dens_ratio, sub_ground) - 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) - endif - if (CS%vmask(Itgt,Jtgt) == 1) then - v_diagonal(Itgt,Jtgt) = v_diagonal(Itgt,Jtgt) + sub_ground(iphi,jphi) * basal_trac(i,j) - endif - enddo ; enddo + + if (CS%umask(I-1,J-1) == 1) u_diag_b(I-1,J-1,4) = u_diag_b(I-1,J-1,4) + sub_ground(1,1) * basal_trac(i,j) + if (CS%umask(I-1,J ) == 1) u_diag_b(I-1,J ,2) = u_diag_b(I-1,J ,2) + sub_ground(1,2) * basal_trac(i,j) + if (CS%umask(I ,J-1) == 1) u_diag_b(I ,J-1,3) = u_diag_b(I ,J-1,3) + sub_ground(2,1) * basal_trac(i,j) + if (CS%umask(I ,J ) == 1) u_diag_b(I ,J ,1) = u_diag_b(I ,J ,1) + sub_ground(2,2) * basal_trac(i,j) + + if (CS%vmask(I-1,J-1) == 1) v_diag_b(I-1,J-1,4) = v_diag_b(I-1,J-1,4) + sub_ground(1,1) * basal_trac(i,j) + if (CS%vmask(I-1,J ) == 1) v_diag_b(I-1,J ,2) = v_diag_b(I-1,J ,2) + sub_ground(1,2) * basal_trac(i,j) + if (CS%vmask(I ,J-1) == 1) v_diag_b(I ,J-1,3) = v_diag_b(I ,J-1,3) + sub_ground(2,1) * basal_trac(i,j) + if (CS%vmask(I ,J ) == 1) v_diag_b(I ,J ,1) = v_diag_b(I ,J ,1) + sub_ground(2,2) * basal_trac(i,j) endif endif ; enddo ; enddo + do J=jsc-2,jec+1 ; do I=isc-2,iec+1 + u_diagonal(I,J) = (u_diag_b(I,J,1)+u_diag_b(I,J,4)) + (u_diag_b(I,J,2)+u_diag_b(I,J,3)) + v_diagonal(I,J) = (v_diag_b(I,J,1)+v_diag_b(I,J,4)) + (v_diag_b(I,J,2)+v_diag_b(I,J,3)) + enddo ; enddo + end subroutine matrix_diagonal -subroutine CG_diagonal_subgrid_basal (Phisub, H_node, bathyT, dens_ratio, sub_grnd) +subroutine CG_diagonal_subgrid_basal (Phisub, H_node, bathyT, dens_ratio, f_grnd) real, dimension(:,:,:,:,:,:), & intent(in) :: Phisub !< Quadrature structure weights at subgridscale !! locations for finite element calculations [nondim] @@ -2735,184 +2890,44 @@ subroutine CG_diagonal_subgrid_basal (Phisub, H_node, bathyT, dens_ratio, sub_gr real, intent(in) :: bathyT !< The depth of ocean bathymetry at tracer points [Z ~> m]. real, intent(in) :: dens_ratio !< The density of ice divided by the density !! of seawater [nondim] - real, dimension(2,2), intent(out) :: sub_grnd !< The weighted fraction of the sub-cell where the ice shelf - !! is grounded [nondim] - - ! bathyT = cellwise-constant bed elevation + real, dimension(2,2), intent(out) :: f_grnd !< The weighted fraction of the sub-cell where the ice shelf + !! is grounded [nondim] + real, dimension(SIZE(Phisub,3),SIZE(Phisub,3),2,2) :: f_grnd_sub ! The contributions to nodal f_grnd + !! from each sub-cell + integer, dimension(2,2,SIZE(Phisub,3),SIZE(Phisub,3)) :: grnd_stat !0 at floating quad points, 1 at grounded + real, dimension(2,2) :: f_grnd_q !Contributions to a node from each quadrature point in a sub-grid cell real :: subarea ! The fractional sub-cell area [nondim] real :: hloc ! The local sub-region thickness [Z ~> m] integer :: nsub, i, j, qx, qy, m, n - nsub = size(Phisub,1) + nsub = size(Phisub,3) subarea = 1.0 / (nsub**2) - sub_grnd(:,:) = 0.0 - do m=1,2 ; do n=1,2 ; do j=1,nsub ; do i=1,nsub ; do qx=1,2 ; do qy = 1,2 + grnd_stat(:,:,:,:)=0 - hloc = (Phisub(i,j,1,1,qx,qy)*H_node(1,1) + Phisub(i,j,2,2,qx,qy)*H_node(2,2)) + & - (Phisub(i,j,1,2,qx,qy)*H_node(1,2) + Phisub(i,j,2,1,qx,qy)*H_node(2,1)) + do j=1,nsub ; do i=1,nsub; do qy=1,2 ; do qx=1,2 + hloc = (Phisub(qx,qy,i,j,1,1)*H_node(1,1) + Phisub(qx,qy,i,j,2,2)*H_node(2,2)) + & + (Phisub(qx,qy,i,j,1,2)*H_node(1,2) + Phisub(qx,qy,i,j,2,1)*H_node(2,1)) + if (dens_ratio * hloc - bathyT > 0) grnd_stat(qx,qy,i,j) = 1 + enddo; enddo ; enddo ; enddo - if (dens_ratio * hloc - bathyT > 0) then - sub_grnd(m,n) = sub_grnd(m,n) + subarea * 0.25 * Phisub(i,j,m,n,qx,qy)**2 - endif + do n=1,2 ; do m=1,2 ; do j=1,nsub ; do i=1,nsub + do qy=1,2 ; do qx = 1,2 + f_grnd_q(qx,qy) = grnd_stat(qx,qy,i,j) * Phisub(qx,qy,i,j,m,n)**2 + enddo ; enddo + !calculate sub-cell contribution to each node by summing up quadrature point contributions from the sub-cell + f_grnd_sub(i,j,m,n) = (subarea * 0.25) * ((f_grnd_q(1,1) + f_grnd_q(2,2)) + (f_grnd_q(1,2)+f_grnd_q(2,1))) + enddo ; enddo ; enddo ; enddo - enddo ; enddo ; enddo ; enddo ; enddo ; enddo + !sum up the sub-cell contributions to each node + do n=1,2 ; do m=1,2 + call sum_square_matrix(f_grnd(m,n),f_grnd_sub(:,:,m,n),nsub) + enddo ; enddo end subroutine CG_diagonal_subgrid_basal -subroutine apply_boundary_values(CS, ISS, G, US, time, Phisub, H_node, ice_visc, basal_trac, float_cond, & - dens_ratio, u_bdry_contr, v_bdry_contr) - - type(ice_shelf_dyn_CS), intent(in) :: 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 - type(time_type), intent(in) :: Time !< The current model time - real, dimension(:,:,:,:,:,:), & - intent(in) :: Phisub !< Quadrature structure weights at subgridscale - !! locations for finite element calculations [nondim] - real, dimension(SZDIB_(G),SZDJB_(G)), & - intent(in) :: H_node !< The ice shelf thickness at nodal - !! (corner) points [Z ~> m]. - real, dimension(SZDI_(G),SZDJ_(G)), & - intent(in) :: ice_visc !< A field related to the ice viscosity from Glen's - !! flow law. The exact form and units depend on the - !! basal law exponent. [R L4 Z T-1 ~> kg m2 s-1]. - real, dimension(SZDI_(G),SZDJ_(G)), & - intent(in) :: basal_trac !< Area-integrated taub_beta field related to the nonlinear - !! part of the "linearized" basal stress [R L3 T-1 ~> kg s-1]. - - real, dimension(SZDI_(G),SZDJ_(G)), & - intent(in) :: float_cond !< An array indicating where the ice - !! shelf is floating: 0 if floating, 1 if not. - real, intent(in) :: dens_ratio !< The density of ice divided by the density - !! of seawater, nondimensional - real, dimension(SZDIB_(G),SZDJB_(G)), & - intent(inout) :: u_bdry_contr !< Zonal force contributions due to the - !! open boundaries [R L3 Z T-2 ~> kg m s-2] - real, dimension(SZDIB_(G),SZDJB_(G)), & - intent(inout) :: v_bdry_contr !< Meridional force contributions due to the - !! open boundaries [R L3 Z T-2 ~> kg m s-2] - -! this will be a per-setup function. the boundary values of thickness and velocity -! (and possibly other variables) will be updated in this function - - real, dimension(8,4) :: Phi - real, dimension(2) :: xquad - real :: ux, uy, vx, vy ! Components of velocity shears or divergence [T-1 ~> s-1] - real :: uq, vq ! Interpolated velocities [L T-1 ~> m s-1] - real, dimension(2,2) :: Ucell,Vcell,Hcell,Usubcontr,Vsubcontr - integer :: i, j, isc, jsc, iec, jec, iq, jq, iphi, jphi, ilq, jlq, Itgt, Jtgt - real :: Ee - - isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec - - xquad(1) = .5 * (1-sqrt(1./3)) ; xquad(2) = .5 * (1+sqrt(1./3)) - - Ee=1.0 - - do j=jsc-1,jec+1 ; do i=isc-1,iec+1 ; if (ISS%hmask(i,j) == 1 .or. ISS%hmask(i,j) == 3) then - - ! process this cell if any corners have umask set to non-dirichlet bdry. - - if ((CS%umask(I-1,J-1) == 3) .OR. (CS%umask(I,J-1) == 3) .OR. & - (CS%umask(I-1,J) == 3) .OR. (CS%umask(I,J) == 3) .OR. & - (CS%vmask(I-1,J-1) == 3) .OR. (CS%vmask(I,J-1) == 3) .OR. & - (CS%vmask(I-1,J) == 3) .OR. (CS%vmask(I,J) == 3)) then - - call bilinear_shape_fn_grid(G, i, j, Phi) - - ! 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 - - uq = CS%u_bdry_val(I-1,J-1) * (xquad(3-iq) * xquad(3-jq)) + & - CS%u_bdry_val(I,J-1) * (xquad(iq) * xquad(3-jq)) + & - CS%u_bdry_val(I-1,J) * (xquad(3-iq) * xquad(jq)) + & - CS%u_bdry_val(I,J) * (xquad(iq) * xquad(jq)) - - vq = CS%v_bdry_val(I-1,J-1) * (xquad(3-iq) * xquad(3-jq)) + & - CS%v_bdry_val(I,J-1) * (xquad(iq) * xquad(3-jq)) + & - CS%v_bdry_val(I-1,J) * (xquad(3-iq) * xquad(jq)) + & - CS%v_bdry_val(I,J) * (xquad(iq) * xquad(jq)) - - ux = CS%u_bdry_val(I-1,J-1) * Phi(1,2*(jq-1)+iq) + & - CS%u_bdry_val(I,J-1) * Phi(3,2*(jq-1)+iq) + & - CS%u_bdry_val(I-1,J) * Phi(5,2*(jq-1)+iq) + & - CS%u_bdry_val(I,J) * Phi(7,2*(jq-1)+iq) - - vx = CS%v_bdry_val(I-1,J-1) * Phi(1,2*(jq-1)+iq) + & - CS%v_bdry_val(I,J-1) * Phi(3,2*(jq-1)+iq) + & - CS%v_bdry_val(I-1,J) * Phi(5,2*(jq-1)+iq) + & - CS%v_bdry_val(I,J) * Phi(7,2*(jq-1)+iq) - - uy = CS%u_bdry_val(I-1,J-1) * Phi(2,2*(jq-1)+iq) + & - CS%u_bdry_val(I,J-1) * Phi(4,2*(jq-1)+iq) + & - CS%u_bdry_val(I-1,J) * Phi(6,2*(jq-1)+iq) + & - CS%u_bdry_val(I,J) * Phi(8,2*(jq-1)+iq) - - vy = CS%v_bdry_val(I-1,J-1) * Phi(2,2*(jq-1)+iq) + & - CS%v_bdry_val(I,J-1) * Phi(4,2*(jq-1)+iq) + & - 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) - - if (trim(CS%ice_viscosity_compute) == "MODEL_QUADRATURE") Ee = CS%Ee(i,j,2*(jq-1)+iq) - - 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 - - if (CS%umask(Itgt,Jtgt) == 1) then - u_bdry_contr(Itgt,Jtgt) = u_bdry_contr(Itgt,Jtgt) + & - 0.25 * Ee * ice_visc(i,j) * ( (4*ux+2*vy) * Phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + & - (uy+vx) * Phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq) ) - - if (float_cond(i,j) == 0) then - u_bdry_contr(Itgt,Jtgt) = u_bdry_contr(Itgt,Jtgt) + & - 0.25 * basal_trac(i,j) * uq * (xquad(ilq) * xquad(jlq)) - endif - endif - - if (CS%vmask(Itgt,Jtgt) == 1) then - v_bdry_contr(Itgt,Jtgt) = v_bdry_contr(Itgt,Jtgt) + & - 0.25 * Ee * ice_visc(i,j) * ( (uy+vx) * Phi(2*(2*(jphi-1)+iphi)-1,2*(jq-1)+iq) + & - (4*vy+2*ux) * Phi(2*(2*(jphi-1)+iphi),2*(jq-1)+iq) ) - - if (float_cond(i,j) == 0) then - v_bdry_contr(Itgt,Jtgt) = v_bdry_contr(Itgt,Jtgt) + & - 0.25 * basal_trac(i,j) * vq * (xquad(ilq) * xquad(jlq)) - endif - endif - enddo ; enddo - enddo ; enddo - - if (float_cond(i,j) == 1) then - Ucell(:,:) = CS%u_bdry_val(i-1:i,j-1:j) ; Vcell(:,:) = CS%v_bdry_val(i-1:i,j-1:j) - Hcell(:,:) = H_node(i-1:i,j-1:j) - call CG_action_subgrid_basal(Phisub, Hcell, Ucell, Vcell, CS%bed_elev(i,j), & - dens_ratio, Usubcontr, Vsubcontr) - - if (CS%umask(I-1,J-1) == 1) u_bdry_contr(I-1,J-1) = u_bdry_contr(I-1,J-1) + Usubcontr(1,1) * basal_trac(i,j) - if (CS%umask(I-1,J) == 1) u_bdry_contr(I-1,J) = u_bdry_contr(I-1,J) + Usubcontr(1,2) * basal_trac(i,j) - if (CS%umask(I,J-1) == 1) u_bdry_contr(I,J-1) = u_bdry_contr(I,J-1) + Usubcontr(2,1) * basal_trac(i,j) - if (CS%umask(I,J) == 1) u_bdry_contr(I,J) = u_bdry_contr(I,J) + Usubcontr(2,2) * basal_trac(i,j) - - if (CS%vmask(I-1,J-1) == 1) v_bdry_contr(I-1,J-1) = v_bdry_contr(I-1,J-1) + Vsubcontr(1,1) * basal_trac(i,j) - if (CS%vmask(I-1,J) == 1) v_bdry_contr(I-1,J) = v_bdry_contr(I-1,J) + Vsubcontr(1,2) * basal_trac(i,j) - if (CS%vmask(I,J-1) == 1) v_bdry_contr(I,J-1) = v_bdry_contr(I,J-1) + Vsubcontr(2,1) * basal_trac(i,j) - if (CS%vmask(I,J) == 1) v_bdry_contr(I,J) = v_bdry_contr(I,J) + Vsubcontr(2,2) * basal_trac(i,j) - endif - endif - endif ; enddo ; enddo - - call pass_vector(u_bdry_contr, v_bdry_contr, G%domain, TO_ALL, BGRID_NE) -end subroutine apply_boundary_values - - !> Update depth integrated viscosity, based on horizontal strain rates subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) type(ice_shelf_dyn_CS), intent(inout) :: CS !< A pointer to the ice shelf control structure @@ -2924,9 +2939,6 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) 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]. - real, pointer, dimension(:,:,:,:) :: Phi => NULL() ! The gradients of bilinear basis elements at Gaussian - ! quadrature points surrounding the cell vertices [L-1 ~> m-1]. - real, pointer, dimension(:,:,:) :: PhiC => NULL() ! Same as Phi, but 1 quadrature point per cell (rather than 4) ! update DEPTH_INTEGRATED viscosity, based on horizontal strain rates - this is for bilinear FEM solve @@ -2938,8 +2950,7 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) real :: Visc_coef, n_g real :: ux, uy, vx, vy real :: eps_min ! Velocity shears [T-1 ~> s-1] -! real, dimension(8,4) :: Phi -! real :: umid, vmid, unorm ! Velocities [L T-1 ~> m s-1] + logical :: model_qp1, model_qp4 isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec iscq = G%iscB ; iecq = G%iecB ; jscq = G%jscB ; jecq = G%jecB @@ -2948,99 +2959,94 @@ 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 + i_off = G%idg_offset ; j_off = G%jdg_offset if (trim(CS%ice_viscosity_compute) == "MODEL") then - allocate(PhiC(1:8,isc:iec,jsc:jec), source=0.0) - do j=jsc,jec ; do i=isc,iec - call bilinear_shape_fn_grid_1qp(G, i, j, PhiC(:,i,j)) - enddo; enddo - elseif (trim(CS%ice_viscosity_compute) == "MODEL_QUADRATURE") then - allocate(Phi(1:8,1:4,isc:iec,jsc:jec), source=0.0) - do j=jsc,jec ; do i=isc,iec - call bilinear_shape_fn_grid(G, i, j, Phi(:,:,i,j)) - enddo; enddo + if (CS%visc_qps==1) then + model_qp1=.true. + model_qp4=.false. + else + model_qp1=.false. + model_qp4=.true. + endif endif n_g = CS%n_glen; eps_min = CS%eps_glen_min - CS%ice_visc(:,:) = 1.0e22 -! 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 if (trim(CS%ice_viscosity_compute) == "CONSTANT") then - CS%ice_visc(i,j) = 1e15 * US%kg_m3_to_R*US%m_to_L*US%m_s_to_L_T * (G%areaT(i,j) * ISS%h_shelf(i,j)) + CS%ice_visc(i,j,1) = 1e15 * (US%kg_m3_to_R*US%m_to_L*US%m_s_to_L_T) * (G%areaT(i,j) * ISS%h_shelf(i,j)) ! constant viscocity for debugging elseif (trim(CS%ice_viscosity_compute) == "OBS") then - if (CS%AGlen_visc(i,j) >0) CS%ice_visc(i,j) = CS%AGlen_visc(i,j)*(G%areaT(i,j) * ISS%h_shelf(i,j)) - ! Here CS%Aglen_visc(i,j) is the ice viscocity [Pa s-1] computed from obs and read from a file - elseif (trim(CS%ice_viscosity_compute) == "MODEL") then - - Visc_coef = ( (US%RL2_T2_to_Pa)**(-CS%n_glen)*US%T_to_s )**(-1./CS%n_glen) * & - (CS%AGlen_visc(i,j))**(-1./CS%n_glen) - ! Units of Aglen_visc [Pa-3 s-1] - - ux = u_shlf(I-1,J-1) * PhiC(1,i,j) + & - u_shlf(I,J) * PhiC(7,i,j) + & - u_shlf(I-1,J) * PhiC(5,i,j) + & - u_shlf(I,J-1) * PhiC(3,i,j) - - vx = v_shlf(I-1,J-1) * PhiC(1,i,j) + & - v_shlf(I,J) * PhiC(7,i,j) + & - v_shlf(I-1,J) * PhiC(5,i,j) + & - v_shlf(I,J-1) * PhiC(3,i,j) - - uy = u_shlf(I-1,J-1) * PhiC(2,i,j) + & - u_shlf(I,J) * PhiC(8,i,j) + & - u_shlf(I-1,J) * PhiC(6,i,j) + & - u_shlf(I,J-1) * PhiC(4,i,j) - - vy = v_shlf(I-1,J-1) * PhiC(2,i,j) + & - v_shlf(I,J) * PhiC(8,i,j) + & - v_shlf(I-1,J) * PhiC(6,i,j) + & - v_shlf(I,J-1) * PhiC(4,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)) - elseif (trim(CS%ice_viscosity_compute) == "MODEL_QUADRATURE") then - !in this case, we will compute viscosity at quadrature points within subroutines CG_action - !and apply_boundary_values. CS%ice_visc(i,j) will include everything except the effective strain rate term: - Visc_coef = ( (US%RL2_T2_to_Pa)**(-CS%n_glen)*US%T_to_s )**(-1./CS%n_glen) * & - (CS%AGlen_visc(i,j))**(-1./CS%n_glen) - CS%ice_visc(i,j) = 0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) + if (CS%AGlen_visc(i,j) >0) CS%ice_visc(i,j,1) = CS%AGlen_visc(i,j) * (G%areaT(i,j) * ISS%h_shelf(i,j)) + ! Here CS%Aglen_visc(i,j) is the ice viscosity [Pa s ~> R L2 T-1] computed from obs and read from a file + elseif (model_qp1) then + !calculate viscosity at 1 cell-centered quadrature point per cell + + Visc_coef = (CS%AGlen_visc(i,j))**(-1./n_g) + ! Units of Aglen_visc [Pa-(n_g) s-1] + + ux = (u_shlf(I-1,J-1) * CS%PhiC(1,i,j) + & + u_shlf(I,J) * CS%PhiC(7,i,j)) + & + (u_shlf(I-1,J) * CS%PhiC(5,i,j) + & + u_shlf(I,J-1) * CS%PhiC(3,i,j)) + + vx = (v_shlf(I-1,J-1) * CS%PhiC(1,i,j) + & + v_shlf(I,J) * CS%PhiC(7,i,j)) + & + (v_shlf(I-1,J) * CS%PhiC(5,i,j) + & + v_shlf(I,J-1) * CS%PhiC(3,i,j)) + + uy = (u_shlf(I-1,J-1) * CS%PhiC(2,i,j) + & + u_shlf(I,J) * CS%PhiC(8,i,j)) + & + (u_shlf(I-1,J) * CS%PhiC(6,i,j) + & + u_shlf(I,J-1) * CS%PhiC(4,i,j)) + + vy = (v_shlf(I-1,J-1) * CS%PhiC(2,i,j) + & + v_shlf(I,J) * CS%PhiC(8,i,j)) + & + (v_shlf(I-1,J) * CS%PhiC(6,i,j) + & + v_shlf(I,J-1) * CS%PhiC(4,i,j)) + + CS%ice_visc(i,j,1) = 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)) * & + (US%Pa_to_RL2_T2*US%s_to_T) + elseif (model_qp4) then + !calculate viscosity at 4 quadrature points per cell + + Visc_coef = (CS%AGlen_visc(i,j))**(-1./n_g) do iq=1,2 ; do jq=1,2 - ux = u_shlf(I-1,J-1) * Phi(1,2*(jq-1)+iq,i,j) + & - u_shlf(I,J-1) * Phi(3,2*(jq-1)+iq,i,j) + & - u_shlf(I-1,J) * Phi(5,2*(jq-1)+iq,i,j) + & - u_shlf(I,J) * Phi(7,2*(jq-1)+iq,i,j) - - vx = v_shlf(I-1,J-1) * Phi(1,2*(jq-1)+iq,i,j) + & - v_shlf(I,J-1) * Phi(3,2*(jq-1)+iq,i,j) + & - v_shlf(I-1,J) * Phi(5,2*(jq-1)+iq,i,j) + & - v_shlf(I,J) * Phi(7,2*(jq-1)+iq,i,j) - - uy = u_shlf(I-1,J-1) * Phi(2,2*(jq-1)+iq,i,j) + & - u_shlf(I,J-1) * Phi(4,2*(jq-1)+iq,i,j) + & - u_shlf(I-1,J) * Phi(6,2*(jq-1)+iq,i,j) + & - u_shlf(I,J) * Phi(8,2*(jq-1)+iq,i,j) - - vy = v_shlf(I-1,J-1) * Phi(2,2*(jq-1)+iq,i,j) + & - v_shlf(I,J-1) * Phi(4,2*(jq-1)+iq,i,j) + & - 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) - - CS%Ee(i,j,2*(jq-1)+iq) = & - (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-1,J-1) * CS%Phi(1,2*(jq-1)+iq,i,j) + & + u_shlf(I,J) * CS%Phi(7,2*(jq-1)+iq,i,j)) + & + (u_shlf(I,J-1) * CS%Phi(3,2*(jq-1)+iq,i,j) + & + u_shlf(I-1,J) * CS%Phi(5,2*(jq-1)+iq,i,j)) + + vx = (v_shlf(I-1,J-1) * CS%Phi(1,2*(jq-1)+iq,i,j) + & + v_shlf(I,J) * CS%Phi(7,2*(jq-1)+iq,i,j)) + & + (v_shlf(I,J-1) * CS%Phi(3,2*(jq-1)+iq,i,j) + & + v_shlf(I-1,J) * CS%Phi(5,2*(jq-1)+iq,i,j)) + + uy = (u_shlf(I-1,J-1) * CS%Phi(2,2*(jq-1)+iq,i,j) + & + u_shlf(I,J) * CS%Phi(8,2*(jq-1)+iq,i,j)) + & + (u_shlf(I,J-1) * CS%Phi(4,2*(jq-1)+iq,i,j) + & + u_shlf(I-1,J) * CS%Phi(6,2*(jq-1)+iq,i,j)) + + vy = (v_shlf(I-1,J-1) * CS%Phi(2,2*(jq-1)+iq,i,j) + & + v_shlf(I,J) * CS%Phi(8,2*(jq-1)+iq,i,j)) + & + (v_shlf(I,J-1) * CS%Phi(4,2*(jq-1)+iq,i,j) + & + v_shlf(I-1,J) * CS%Phi(6,2*(jq-1)+iq,i,j)) + + CS%ice_visc(i,j,2*(jq-1)+iq) = 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)) * & + (US%Pa_to_RL2_T2*US%s_to_T) enddo; enddo endif endif enddo ; enddo - if (trim(CS%ice_viscosity_compute) == "MODEL") deallocate(PhiC) - if (trim(CS%ice_viscosity_compute) == "MODEL_QUADRATURE") deallocate(Phi) end subroutine calc_shelf_visc @@ -3065,8 +3071,9 @@ subroutine calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) real :: umid, vmid, unorm, eps_min ! Velocities [L T-1 ~> m s-1] real :: alpha !Coulomb coefficient [nondim] real :: Hf !"floatation thickness" for Coulomb friction [Z ~> m] - real :: fN !Effective pressure (ice pressure - ocean pressure) for Coulomb friction [R L2 T-2 ~> Pa] - real :: fB !for Coulomb Friction [(L T-1)^CS%CF_PostPeak ~> (m s-1)^CS%CF_PostPeak] + real :: fN !Effective pressure (ice pressure - ocean pressure) for Coulomb friction [Pa] + real :: fB !for Coulomb Friction [(T L-1)^CS%CF_PostPeak ~> (s m-1)^CS%CF_PostPeak] + real :: fN_scale !To convert effective pressure to mks units during Coulomb friction [Pa T2 R-1 L-2 ~> 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 @@ -3079,11 +3086,12 @@ subroutine calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) eps_min = CS%eps_glen_min if (CS%CoulombFriction) then - if (CS%CF_PostPeak.ne.1.0) THEN + if (CS%CF_PostPeak/=1.0) THEN alpha = (CS%CF_PostPeak-1.0)**(CS%CF_PostPeak-1.0) / CS%CF_PostPeak**CS%CF_PostPeak ![nondim] else alpha = 1.0 endif + fN_scale = US%R_to_kg_m3 * US%L_T_to_m_s**2 endif do j=jsd+1,jed @@ -3091,20 +3099,22 @@ subroutine calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) 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 = US%L_T_to_m_s*sqrt(umid**2 + vmid**2 + eps_min**2*(G%dxT(i,j)**2 + G%dyT(i,j)**2)) + unorm = US%L_T_to_m_s * sqrt( (umid**2 + vmid**2) + (eps_min**2 * (G%dxT(i,j)**2 + G%dyT(i,j)**2)) ) !Coulomb friction (Schoof 2005, Gagliardini et al 2007) if (CS%CoulombFriction) then !Effective pressure - Hf = max(CS%density_ocean_avg * CS%bed_elev(i,j)/CS%density_ice, 0.0) - fN = max(CS%density_ice * CS%g_Earth * (ISS%h_shelf(i,j) - Hf),CS%CF_MinN) - + Hf = max((CS%density_ocean_avg/CS%density_ice) * CS%bed_elev(i,j), 0.0) + fN = max(fN_scale*((CS%density_ice * CS%g_Earth) * (ISS%h_shelf(i,j) - Hf)),CS%CF_MinN) fB = alpha * (CS%C_basal_friction(i,j) / (CS%CF_Max * fN))**(CS%CF_PostPeak/CS%n_basal_fric) - CS%basal_traction(i,j) = G%areaT(i,j) * CS%C_basal_friction(i,j) * & - unorm**(CS%n_basal_fric-1.0) / (1.0 + fB * unorm**CS%CF_PostPeak)**(CS%n_basal_fric) + + CS%basal_traction(i,j) = ((G%areaT(i,j) * CS%C_basal_friction(i,j)) * & + (unorm**(CS%n_basal_fric-1.0) / (1.0 + fB * unorm**CS%CF_PostPeak)**(CS%n_basal_fric))) * & + (US%Pa_to_RLZ_T2*US%L_T_to_m_s) else - !linear (CS%n_basal_fric=1) or "Weertman"/power-law (CS%n_basal_fric .ne. 1) - CS%basal_traction(i,j) = G%areaT(i,j) * CS%C_basal_friction(i,j) * unorm**(CS%n_basal_fric-1) + !linear (CS%n_basal_fric=1) or "Weertman"/power-law (CS%n_basal_fric /= 1) + CS%basal_traction(i,j) = ((G%areaT(i,j) * CS%C_basal_friction(i,j)) * (unorm**(CS%n_basal_fric-1))) * & + (US%Pa_to_RLZ_T2*US%L_T_to_m_s) endif endif enddo @@ -3353,8 +3363,8 @@ subroutine bilinear_shape_fn_grid(G, i, j, Phi) xexp = xquad(qpoint) endif - Phi(2*node-1,qpoint) = ( d * (2 * xnode - 3) * yexp ) / (a*d) - Phi(2*node,qpoint) = ( a * (2 * ynode - 3) * xexp ) / (a*d) + Phi(2*node-1,qpoint) = ( (d * (2 * xnode - 3)) * yexp ) / (a*d) + Phi(2*node,qpoint) = ( (a * (2 * ynode - 3)) * xexp ) / (a*d) enddo enddo @@ -3400,15 +3410,15 @@ subroutine bilinear_shape_fn_grid_1qp(G, i, j, Phi) do node=1,4 xnode = 2-mod(node,2) ; ynode = ceiling(REAL(node)/2) - Phi(2*node-1) = ( d * (2 * xnode - 3) * yexp ) / (a*d) - Phi(2*node) = ( a * (2 * ynode - 3) * xexp ) / (a*d) + Phi(2*node-1) = ( (d * (2 * xnode - 3)) * yexp ) / (a*d) + Phi(2*node) = ( (a * (2 * ynode - 3)) * xexp ) / (a*d) enddo end subroutine bilinear_shape_fn_grid_1qp subroutine bilinear_shape_functions_subgrid(Phisub, nsub) integer, intent(in) :: nsub !< The number of subgridscale quadrature locations in each direction - real, dimension(nsub,nsub,2,2,2,2), & + real, dimension(2,2,nsub,nsub,2,2), & intent(inout) :: Phisub !< Quadrature structure weights at subgridscale !! locations for finite element calculations [nondim] @@ -3420,13 +3430,13 @@ subroutine bilinear_shape_functions_subgrid(Phisub, nsub) ! i think this general approach may not work for nonrectangular elements... ! - ! Phisub(i,j,k,l,q1,q2) + ! Phisub(q1,q2,i,j,k,l) + ! q1: quad point x-index + ! q2: quad point y-index ! i: subgrid index in x-direction ! j: subgrid index in y-direction ! k: basis function x-index ! l: basis function y-index - ! q1: quad point x-index - ! q2: quad point y-index ! e.g. k=1,l=1 => node 1 ! q1=2,q2=1 => quad point 2 @@ -3447,10 +3457,10 @@ subroutine bilinear_shape_functions_subgrid(Phisub, nsub) do qy=1,2 ; do qx=1,2 x = x0 + fracx*xquad(qx) y = y0 + fracx*xquad(qy) - Phisub(i,j,1,1,qx,qy) = (1.0-x) * (1.0-y) - Phisub(i,j,1,2,qx,qy) = (1.0-x) * y - Phisub(i,j,2,1,qx,qy) = x * (1.0-y) - Phisub(i,j,2,2,qx,qy) = x * y + Phisub(qx,qy,i,j,1,1) = (1.0-x) * (1.0-y) + Phisub(qx,qy,i,j,1,2) = (1.0-x) * y + Phisub(qx,qy,i,j,2,1) = x * (1.0-y) + Phisub(qx,qy,i,j,2,2) = x * y enddo ; enddo enddo ; enddo @@ -3623,8 +3633,8 @@ subroutine interpolate_H_to_B(G, h_shelf, hmask, H_node) intent(inout) :: H_node !< The ice shelf thickness at nodal (corner) !! points [Z ~> m]. - integer :: i, j, isc, iec, jsc, jec, num_h, k, l - real :: summ + integer :: i, j, isc, iec, jsc, jec, num_h, k, l, ic, jc + real :: h_arr(2,2) isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec @@ -3635,19 +3645,18 @@ subroutine interpolate_H_to_B(G, h_shelf, hmask, H_node) do j=jsc-1,jec do i=isc-1,iec - summ = 0.0 num_h = 0 - do k=0,1 - do l=0,1 - if (hmask(i+k,j+l) == 1.0 .or. hmask(i+k,j+l) == 3.0) then - summ = summ + h_shelf(i+k,j+l) - num_h = num_h + 1 - endif - enddo - enddo - if (num_h > 0) then - H_node(i,j) = summ / num_h - endif + do l=1,2; jc=j-1+l; do k=1,2; ic=i-1+k + if (hmask(ic,jc) == 1.0 .or. hmask(ic,jc) == 3.0) then + h_arr(k,l)=h_shelf(ic,jc) + num_h = num_h + 1 + else + h_arr(k,l)=0.0 + endif + if (num_h > 0) then + H_node(i,j) = ((h_arr(1,1)+h_arr(2,2))+(h_arr(1,2)+h_arr(2,1))) / num_h + endif + enddo; enddo enddo enddo @@ -3667,9 +3676,11 @@ subroutine ice_shelf_dyn_end(CS) deallocate(CS%u_bdry_val, CS%v_bdry_val) deallocate(CS%u_face_mask, CS%v_face_mask) deallocate(CS%umask, CS%vmask) + deallocate(CS%u_face_mask_bdry, CS%v_face_mask_bdry) + deallocate(CS%h_bdry_val) + deallocate(CS%float_cond) deallocate(CS%ice_visc, CS%AGlen_visc) - deallocate(CS%Ee) deallocate(CS%basal_traction,CS%C_basal_friction) deallocate(CS%OD_rt, CS%OD_av) deallocate(CS%t_bdry_val, CS%bed_elev) diff --git a/src/ice_shelf/MOM_ice_shelf_initialize.F90 b/src/ice_shelf/MOM_ice_shelf_initialize.F90 index 1e2076f889..f976187c2b 100644 --- a/src/ice_shelf/MOM_ice_shelf_initialize.F90 +++ b/src/ice_shelf/MOM_ice_shelf_initialize.F90 @@ -274,8 +274,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,& - US, PF ) + hmask, h_shelf, G, US, PF ) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure real, dimension(SZIB_(G),SZJB_(G)), & @@ -299,10 +298,6 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b 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]. - 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)), & @@ -350,8 +345,13 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b if (G%geoLonBu(i,j) == westlon) then hmask(i+1,j) = 3.0 - h_bdry_val(i+1,j) = h_shelf(i+1,j) - thickness_bdry_val(i+1,j) = h_bdry_val(i+0*1,j) + !--- + !OLD: thickness_bdry_val was used for ice dynamics, and h_bdry_val was not used anywhere except here: + !h_bdry_val(i+1,j) = h_shelf(i+1,j) ; thickness_bdry_val(i+1,j) = h_bdry_val(i+0*1,j) + !--- + !NEW: h_bdry_val is used for ice dynamics instead of thickness_bdry_val, which was removed + h_bdry_val(i+1,j) = h_shelf(i+0*1,j) !why 0*1 + !--- u_face_mask_bdry(i+1,j) = 5.0 u_bdry_val(i+1,j) = input_vel*(1-16.0*((G%geoLatBu(i-1,j)/lenlat-0.5))**4) !velocity distribution endif @@ -393,8 +393,6 @@ end subroutine initialize_ice_shelf_boundary_channel !> Initialize ice shelf flow from file subroutine initialize_ice_flow_from_file(bed_elev,u_shelf, v_shelf,float_cond,& G, US, PF) -!subroutine initialize_ice_flow_from_file(bed_elev,u_shelf, v_shelf,float_cond,ice_visc,& -! G, US, PF) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure real, dimension(SZDI_(G),SZDJ_(G)), & intent(inout) :: bed_elev !< The bed elevation [Z ~> m]. @@ -412,9 +410,8 @@ subroutine initialize_ice_flow_from_file(bed_elev,u_shelf, v_shelf,float_cond,& ! h_shelf [Z ~> m] and area_shelf_h [L2 ~> m2] (and dimensionless) and updates hmask character(len=200) :: filename,vel_file,inputdir,bed_topo_file ! Strings for file/path character(len=200) :: ushelf_varname, vshelf_varname, & - ice_visc_varname, floatfr_varname, bed_varname ! Variable name in file + floatfr_varname, bed_varname ! Variable name in file character(len=40) :: mdl = "initialize_ice_velocity_from_file" ! This subroutine's name. - real :: len_sidestress call MOM_mesg(" MOM_ice_shelf_init_profile.F90, initialize_velocity_from_file: reading velocity") @@ -423,9 +420,6 @@ subroutine initialize_ice_flow_from_file(bed_elev,u_shelf, v_shelf,float_cond,& 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) @@ -435,9 +429,6 @@ subroutine initialize_ice_flow_from_file(bed_elev,u_shelf, v_shelf,float_cond,& call get_param(PF, mdl, "ICE_V_VEL_VARNAME", vshelf_varname, & "The name of the v velocity variable in ICE_VELOCITY_FILE.", & default="v_shelf") - call get_param(PF, mdl, "ICE_VISC_VARNAME", ice_visc_varname, & - "The name of the ice viscosity variable in ICE_VELOCITY_FILE.", & - default="viscosity") call get_param(PF, mdl, "ICE_FLOAT_FRAC_VARNAME", floatfr_varname, & "The name of the ice float fraction (grounding fraction) variable in ICE_VELOCITY_FILE.", & default="float_frac") @@ -462,7 +453,7 @@ end subroutine initialize_ice_flow_from_file !> Initialize ice shelf b.c.s from file subroutine initialize_ice_shelf_boundary_from_file(u_face_mask_bdry, v_face_mask_bdry, & - u_bdry_val, v_bdry_val, umask, vmask, h_bdry_val, thickness_bdry_val, & + u_bdry_val, v_bdry_val, umask, vmask, h_bdry_val, & hmask, h_shelf, G, US, PF ) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure @@ -480,8 +471,6 @@ subroutine initialize_ice_shelf_boundary_from_file(u_face_mask_bdry, v_face_mask intent(inout) :: umask !< A mask for ice shelf velocity [nondim] real, dimension(SZDIB_(G),SZDJB_(G)), & intent(inout) :: vmask !< A mask for ice shelf velocity [nondim] - real, dimension(SZDI_(G),SZDJ_(G)), & - intent(inout) :: thickness_bdry_val !< The ice shelf thickness at open boundaries [Z ~> m] 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)), & @@ -501,7 +490,6 @@ subroutine initialize_ice_shelf_boundary_from_file(u_face_mask_bdry, v_face_mask integer :: i, j, isc, jsc, iec, jec h_bdry_val(:,:) = 0. - thickness_bdry_val(:,:) = 0. call MOM_mesg(" MOM_ice_shelf_init_profile.F90, initialize_b_c_s_from_file: reading b.c.s") @@ -542,9 +530,9 @@ subroutine initialize_ice_shelf_boundary_from_file(u_face_mask_bdry, v_face_mask call MOM_read_data(filename, trim(ufcmskbdry_varname), u_face_mask_bdry, G%Domain, position=CORNER, & - scale=US%m_s_to_L_T) + scale=1.) call MOM_read_data(filename, trim(vfcmskbdry_varname), v_face_mask_bdry, G%Domain, position=CORNER, & - scale=US%m_s_to_L_T) + scale=1.) call MOM_read_data(filename, trim(ubdryv_varname), u_bdry_val, G%Domain, position=CORNER, scale=US%m_s_to_L_T) call MOM_read_data(filename, trim(vbdryv_varname), v_bdry_val, G%Domain, position=CORNER, scale=US%m_s_to_L_T) call MOM_read_data(filename, trim(umask_varname), umask, G%Domain, position=CORNER, scale=1.) @@ -557,7 +545,6 @@ subroutine initialize_ice_shelf_boundary_from_file(u_face_mask_bdry, v_face_mask do j=jsc,jec do i=isc,iec if (hmask(i,j) == 3.) then - thickness_bdry_val(i,j) = h_shelf(i,j) h_bdry_val(i,j) = h_shelf(i,j) endif enddo @@ -587,7 +574,7 @@ subroutine initialize_ice_C_basal_friction(C_basal_friction, G, US, PF) if (trim(config)=="CONSTANT") then call get_param(PF, mdl, "BASAL_FRICTION_COEFF", C_friction, & - "Coefficient in sliding law.", units="Pa (m s-1)^(n_basal_fric)", default=5.e10) + "Coefficient in sliding law.", units="Pa (s m-1)^(n_basal_fric)", default=5.e10) C_basal_friction(:,:) = C_friction elseif (trim(config)=="FILE") then @@ -615,10 +602,13 @@ subroutine initialize_ice_C_basal_friction(C_basal_friction, G, US, PF) !> Initialize ice-stiffness parameter -subroutine initialize_ice_AGlen(AGlen, G, US, PF) +subroutine initialize_ice_AGlen(AGlen, ice_viscosity_compute, G, US, PF) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure real, dimension(SZDI_(G),SZDJ_(G)), & intent(inout) :: AGlen !< The ice-stiffness parameter A_Glen, often in [Pa-3 s-1] + character(len=40) :: ice_viscosity_compute !< Specifies whether the ice viscosity is computed internally + !! according to Glen's flow law; is constant (for debugging purposes) + !! or using observed strain rates and read from a file 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 @@ -635,7 +625,7 @@ subroutine initialize_ice_AGlen(AGlen, G, US, PF) if (trim(config)=="CONSTANT") then call get_param(PF, mdl, "A_GLEN", A_Glen, & - "Ice-stiffness parameter.", units="Pa-3 s-1", default=2.261e-25) + "Ice-stiffness parameter.", units="Pa-n_g s-1", default=2.261e-25) AGlen(:,:) = A_Glen @@ -655,8 +645,14 @@ subroutine initialize_ice_AGlen(AGlen, G, US, PF) if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, & " initialize_ice_stiffness_from_file: Unable to open "//trim(filename)) - call MOM_read_data(filename,trim(varname), AGlen, G%Domain) + if (trim(ice_viscosity_compute) == "OBS") then + !AGlen is the ice viscosity [Pa s ~> R L2 T-1] computed from obs and read from a file + call MOM_read_data(filename, trim(varname), AGlen, G%Domain, scale=US%Pa_to_RL2_T2*US%s_to_T) + else + !AGlen is the ice stiffness parameter [Pa-n_g s-1] + call MOM_read_data(filename, trim(varname), AGlen, G%Domain) + endif endif end subroutine initialize_ice_AGlen @@ -681,7 +677,7 @@ subroutine initialize_ice_SMB(SMB, G, US, PF) if (trim(config)=="CONSTANT") then call get_param(PF, mdl, "SMB", SMB_val, & - "Surface mass balance.", units="kg m-2 s-1", default=0.0) + "Surface mass balance.", units="kg m-2 s-1", default=0.0, scale=US%kg_m2s_to_RZ_T) SMB(:,:) = SMB_val @@ -701,7 +697,7 @@ subroutine initialize_ice_SMB(SMB, G, US, PF) if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, & " initialize_ice_SMV_from_file: Unable to open "//trim(filename)) - call MOM_read_data(filename,trim(varname), SMB, G%Domain) + call MOM_read_data(filename,trim(varname), SMB, G%Domain, scale=US%kg_m2s_to_RZ_T) endif end subroutine initialize_ice_SMB