From a6c5348fd6393af900f4aae79f51a38981961a8b Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Tue, 14 Dec 2021 00:11:47 -0500 Subject: [PATCH 1/4] Hydrostatic initialization in ice cavities - Iteratively solve for the initial ice shelf displacement in cavities by calculating the pressure at the current displacement depth using the unperturbed profile. - This change should obsolete TRIM_IC_FOR_PSURF and DEPRESS_INITIAL_SURFACE for ice shelf applications and should work for arbitrary equations of state. - Existing implementations (e.g. ISOMIP) should turn off the above options in order to exercise this feature. - This code change should not impact non ice-shelf configurations or those with either of the above two options. --- src/core/MOM.F90 | 28 ++- src/ice_shelf/MOM_ice_shelf.F90 | 14 +- src/ice_shelf/MOM_ice_shelf_dynamics.F90 | 118 +++++++++---- src/ice_shelf/MOM_ice_shelf_initialize.F90 | 45 ++--- .../MOM_state_initialization.F90 | 159 +++++++++++++++--- 5 files changed, 261 insertions(+), 103 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 0e036d9d8f..dfef140497 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -278,7 +278,8 @@ module MOM type(time_type) :: dtbt_reset_interval !< A time_time representation of dtbt_reset_period. type(time_type) :: dtbt_reset_time !< The next time DTBT should be calculated. real, dimension(:,:), pointer :: frac_shelf_h => NULL() !< fraction of total area occupied - !! by ice shelf [nondim] + !! by ice shelf [nondim] + real, dimension(:,:), pointer :: mass_shelf => NULL() !< Mass of ice shelf [R Z ~> kg m-2] real, dimension(:,:,:), pointer :: & h_pre_dyn => NULL(), & !< The thickness before the transports [H ~> m or kg m-2]. T_pre_dyn => NULL(), & !< Temperature before the transports [degC]. @@ -1724,9 +1725,12 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & real, allocatable :: v_in(:,:,:) ! Initial meridional velocities [L T-1 ~> m s-1] real, allocatable :: h_in(:,:,:) ! Initial layer thicknesses [H ~> m or kg m-2] real, allocatable, target :: frac_shelf_in(:,:) ! Initial fraction of the total cell area occupied - ! by an ice shelf [nondim] + ! by an ice shelf [nondim] + real, allocatable, target :: mass_shelf_in(:,:) ! Initial mass of ice shelf contained within a grid cell + ! [RZ ~> kg m-2] real, allocatable, target :: T_in(:,:,:) ! Initial temperatures [degC] real, allocatable, target :: S_in(:,:,:) ! Initial salinities [ppt] + type(ocean_OBC_type), pointer :: OBC_in => NULL() type(sponge_CS), pointer :: sponge_in_CSp => NULL() type(ALE_sponge_CS), pointer :: ALE_sponge_in_CSp => NULL() @@ -2491,14 +2495,18 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & ! for legacy reasons. The actual ice shelf diag CS is internal to the ice shelf call initialize_ice_shelf(param_file, G_in, Time, ice_shelf_CSp, diag_ptr) allocate(frac_shelf_in(G_in%isd:G_in%ied, G_in%jsd:G_in%jed), source=0.0) + allocate(mass_shelf_in(G_in%isd:G_in%ied, G_in%jsd:G_in%jed), source=0.0) allocate(CS%frac_shelf_h(isd:ied, jsd:jed), source=0.0) - call ice_shelf_query(ice_shelf_CSp,G,CS%frac_shelf_h) + allocate(CS%mass_shelf(isd:ied, jsd:jed), source=0.0) + call ice_shelf_query(ice_shelf_CSp,G,CS%frac_shelf_h, CS%mass_shelf) + CS%mass_shelf = CS%mass_shelf*US%kg_m3_to_R*US%m_to_Z ! MOM_initialize_state is using the unrotated metric call rotate_array(CS%frac_shelf_h, -turns, frac_shelf_in) + call rotate_array(CS%mass_shelf, -turns, mass_shelf_in) call MOM_initialize_state(u_in, v_in, h_in, CS%tv, Time, G_in, GV, US, & param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, & sponge_in_CSp, ALE_sponge_in_CSp, oda_incupd_in_CSp, OBC_in, Time_in, & - frac_shelf_h=frac_shelf_in) + frac_shelf_h=frac_shelf_in, mass_shelf = mass_shelf_in) else call MOM_initialize_state(u_in, v_in, h_in, CS%tv, Time, G_in, GV, US, & param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, & @@ -2542,16 +2550,18 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & deallocate(S_in) endif if (use_ice_shelf) & - deallocate(frac_shelf_in) + deallocate(frac_shelf_in,mass_shelf_in) else if (use_ice_shelf) then call initialize_ice_shelf(param_file, G, Time, ice_shelf_CSp, diag_ptr) allocate(CS%frac_shelf_h(isd:ied, jsd:jed), source=0.0) - call ice_shelf_query(ice_shelf_CSp,G,CS%frac_shelf_h) + allocate(CS%mass_shelf(isd:ied, jsd:jed), source=0.0) + call ice_shelf_query(ice_shelf_CSp,G,CS%frac_shelf_h, CS%mass_shelf) + CS%mass_shelf = CS%mass_shelf*US%kg_m3_to_R*US%m_to_Z call MOM_initialize_state(CS%u, CS%v, CS%h, CS%tv, Time, G, GV, US, & param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, & CS%sponge_CSp, CS%ALE_sponge_CSp,CS%oda_incupd_CSp, CS%OBC, Time_in, & - frac_shelf_h=CS%frac_shelf_h) + frac_shelf_h=CS%frac_shelf_h, mass_shelf=CS%mass_shelf) else call MOM_initialize_state(CS%u, CS%v, CS%h, CS%tv, Time, G, GV, US, & param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, & @@ -2566,8 +2576,10 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & endif endif - if (use_ice_shelf .and. CS%debug) & + if (use_ice_shelf .and. CS%debug) then call hchksum(CS%frac_shelf_h, "MOM:frac_shelf_h", G%HI, haloshift=0) + call hchksum(CS%mass_shelf, "MOM:mass_shelf", G%HI, haloshift=0) + endif call cpu_clock_end(id_clock_MOM_init) call callTree_waypoint("returned from MOM_initialize_state() (initialize_MOM)") diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 9dd3791211..88f51571cd 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -2032,11 +2032,12 @@ subroutine update_shelf_mass(G, US, CS, ISS, Time) end subroutine update_shelf_mass !> Save the ice shelf restart file -subroutine ice_shelf_query(CS, G, frac_shelf_h) +subroutine ice_shelf_query(CS, G, frac_shelf_h, mass_shelf) type(ice_shelf_CS), pointer :: CS !< ice shelf control structure type(ocean_grid_type), intent(in) :: G !< A pointer to an ocean grid control structure. - real, optional, dimension(SZI_(G),SZJ_(G)) :: frac_shelf_h !< - !< Ice shelf area fraction [nodim]. + real, optional, dimension(SZI_(G),SZJ_(G)) :: frac_shelf_h !< Ice shelf area fraction [nodim]. + real, optional, dimension(SZI_(G),SZJ_(G)) :: mass_shelf ! kg m-2] + integer :: i, j @@ -2047,6 +2048,13 @@ subroutine ice_shelf_query(CS, G, frac_shelf_h) enddo ; enddo endif + if (present(mass_shelf)) then + do j=G%jsd,G%jed ; do i=G%isd,G%ied + mass_shelf(i,j) = 0.0 + if (G%areaT(i,j)>0.) mass_shelf(i,j) = CS%ISS%mass_shelf(i,j) + enddo ; enddo + endif + end subroutine ice_shelf_query !> Save the ice shelf restart file diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index df2e801613..76da7ecff3 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -260,9 +260,9 @@ subroutine register_ice_shelf_dyn_restarts(G, param_file, CS, restart_CS) allocate( CS%v_shelf(IsdB:IedB,JsdB:JedB), source=0.0 ) allocate( CS%t_shelf(isd:ied,jsd:jed), source=-10.0 ) ! [degC] allocate( CS%ice_visc(isd:ied,jsd:jed), source=0.0 ) - allocate( CS%AGlen_visc(isd:ied,jsd:jed), source=2.261e-25 ) ! [Units?] + allocate( CS%AGlen_visc(isd:ied,jsd:jed), source=2.261e-25 ) ! [Pa-3s-1] allocate( CS%basal_traction(isd:ied,jsd:jed), source=0.0 ) - allocate( CS%C_basal_friction(isd:ied,jsd:jed), source=5.0e10 ) ! [Units?] + allocate( CS%C_basal_friction(isd:ied,jsd:jed), source=5.0e10 ) ! [Pa (m-1 s)^n_sliding] allocate( CS%OD_av(isd:ied,jsd:jed), source=0.0 ) allocate( CS%ground_frac(isd:ied,jsd:jed), source=0.0 ) allocate( CS%taudx_shelf(IsdB:IedB,JsdB:JedB), source=0.0 ) @@ -553,8 +553,8 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ call pass_vector(CS%u_bdry_val, CS%v_bdry_val, G%domain, TO_ALL, BGRID_NE) call pass_vector(CS%u_face_mask_bdry, CS%v_face_mask_bdry, G%domain, TO_ALL, BGRID_NE) - !initialize ice flow velocities from file - call initialize_ice_flow_from_file(CS%bed_elev,CS%u_shelf, CS%v_shelf,CS%ground_frac, ISS%hmask,ISS%h_shelf, & + !initialize ice flow characteristic (velocities, bed elevation under the grounded part, etc) from file + call initialize_ice_flow_from_file(CS%bed_elev,CS%u_shelf, CS%v_shelf,CS%ground_frac, & G, US, param_file) call pass_vector(CS%u_shelf, CS%v_shelf, G%domain, TO_ALL, BGRID_NE) call pass_var(CS%bed_elev, G%domain,CENTER) @@ -567,9 +567,9 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ 'y-velocity of ice', 'm yr-1', conversion=365.0*86400.0*US%L_T_to_m_s) ! 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-9*US%RL2_T2_to_Pa) + 'x-driving stress of ice', 'kPa', conversion=1.e-3*US%RL2_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-9*US%RL2_T2_to_Pa) + 'y-driving stress of ice', 'kPa', conversion=1.e-3*US%RL2_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, & @@ -579,9 +579,9 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ CS%id_col_thick = register_diag_field('ice_shelf_model','col_thick',CS%diag%axesT1, Time, & 'ocean column thickness passed to ice model', 'm', conversion=US%Z_to_m) CS%id_visc_shelf = register_diag_field('ice_shelf_model','ice_visc',CS%diag%axesT1, Time, & - 'viscosity', 'm', conversion=1e-6*US%Z_to_m) + 'vi-viscosity', 'Pa s-1 m', conversion=US%Z_to_m) !vertically integrated viscosity CS%id_taub = register_diag_field('ice_shelf_model','taub_beta',CS%diag%axesT1, Time, & - 'taub', 'Pa yr m-1', conversion=1e-6*US%Z_to_m) + 'taub', 'MPa', conversion=1e-6*US%RL2_T2_to_Pa) 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 @@ -673,7 +673,10 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled logical, optional, intent(in) :: coupled_grounding !< If true, the grounding line is !! determined by coupled ice-ocean dynamics logical, optional, intent(in) :: must_update_vel !< Always update the ice velocities if true. - + real, dimension(SZDIB_(G),SZDJB_(G)) ::taud_x,taud_y ! Pa] + real, dimension(SZDI_(G),SZDJ_(G)) :: ice_visc ! Pa s-1 m] + real, dimension(SZDI_(G),SZDJ_(G)) :: basal_tr ! Pa] integer :: iters logical :: update_ice_vel, coupled_GL @@ -706,12 +709,24 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled if (CS%id_u_shelf > 0) call post_data(CS%id_u_shelf, CS%u_shelf, CS%diag) if (CS%id_v_shelf > 0) call post_data(CS%id_v_shelf, CS%v_shelf, CS%diag) ! if (CS%id_t_shelf > 0) call post_data(CS%id_t_shelf,CS%t_shelf,CS%diag) - if (CS%id_taudx_shelf > 0) call post_data(CS%id_taudx_shelf, CS%taudx_shelf, CS%diag) - if (CS%id_taudy_shelf > 0) call post_data(CS%id_taudy_shelf, CS%taudy_shelf, CS%diag) + if (CS%id_taudx_shelf > 0) then + taud_x(:,:) = CS%taudx_shelf(:,:)/G%areaT(:,:) + call post_data(CS%id_taudx_shelf,taud_x , CS%diag) + endif + if (CS%id_taudy_shelf > 0) then + taud_y(:,:) = CS%taudy_shelf(:,:)/G%areaT(:,:) + 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_OD_av >0) call post_data(CS%id_OD_av, CS%OD_av,CS%diag) - if (CS%id_visc_shelf > 0) call post_data(CS%id_visc_shelf, CS%ice_visc,CS%diag) - if (CS%id_taub > 0) call post_data(CS%id_taub, CS%basal_traction,CS%diag) + if (CS%id_visc_shelf > 0) then + ice_visc(:,:)=CS%ice_visc(:,:)/G%areaT(:,:) + call post_data(CS%id_visc_shelf, ice_visc,CS%diag) + endif + if (CS%id_taub > 0) then + basal_tr(:,:) = CS%basal_traction(:,:)/G%areaT(:,:) + call post_data(CS%id_taub, basal_tr,CS%diag) + endif !! if (CS%id_u_mask > 0) call post_data(CS%id_u_mask,CS%umask,CS%diag) if (CS%id_v_mask > 0) call post_data(CS%id_v_mask,CS%vmask,CS%diag) @@ -874,6 +889,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i if (rhoi_rhow * ISS%h_shelf(i,j) - CS%bed_elev(i,j) > 0) then float_cond(i,j) = 1.0 CS%ground_frac(i,j) = 1.0 + CS%OD_av(i,j) =0.0 endif enddo enddo @@ -960,7 +976,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i !! begin loop - do iter=1,100 + 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) @@ -1775,7 +1791,7 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) intent(inout) :: taudx !< X-direction driving stress at q-points [kg L s-2 ~> kg m s-2] real, dimension(SZDIB_(G),SZDJB_(G)), & intent(inout) :: taudy !< Y-direction driving stress at q-points [kg L s-2 ~> kg m s-2] - ! This will become [R L3 Z T-2 ~> kg m s-2] + ! This will become [R L3 Z T-2 ~> kg m s-2] ! driving stress! @@ -1790,12 +1806,14 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) real, dimension(SIZE(OD,1),SIZE(OD,2)) :: S, & ! surface elevation [Z ~> m]. BASE ! basal elevation of shelf/stream [Z ~> m]. + real, pointer, dimension(:,:,:,:) :: Phi => NULL() ! The gradients of bilinear basis elements at Gaussian + ! quadrature points surrounding the cell vertices [m-1]. 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] - real :: dxh, dyh ! Local grid spacing [L ~> m] + real :: dxh, dyh,Dx,Dy ! Local grid spacing [L ~> m] real :: grav ! The gravitational acceleration [L2 Z-1 T-2 ~> m s-2] integer :: i, j, iscq, iecq, jscq, jecq, isd, jsd, is, js, iegq, jegq @@ -1813,6 +1831,7 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) is = iscq - 1; js = jscq - 1 i_off = G%idg_offset ; j_off = G%jdg_offset + rho = CS%density_ice rhow = CS%density_ocean_avg grav = CS%g_Earth @@ -1821,13 +1840,14 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) ! or is this faster? BASE(:,:) = -CS%bed_elev(:,:) + OD(:,:) - S(:,:) = BASE(:,:) + ISS%h_shelf(:,:) - + S(:,:) = -CS%bed_elev(:,:) + ISS%h_shelf(:,:) ! check whether the ice is floating or grounded do j=jsc-G%domain%njhalo,jec+G%domain%njhalo do i=isc-G%domain%nihalo,iec+G%domain%nihalo if (rhoi_rhow * ISS%h_shelf(i,j) - CS%bed_elev(i,j) <= 0) then S(i,j) = (1 - rhoi_rhow)*ISS%h_shelf(i,j) + else + S(i,j)=ISS%h_shelf(i,j)-CS%bed_elev(i,j) endif enddo enddo @@ -1838,7 +1858,8 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) sy = 0 dxh = G%dxT(i,j) dyh = G%dyT(i,j) - + Dx=dxh + Dy=dyh if (ISS%hmask(i,j) == 1) then ! we are inside the global computational bdry, at an ice-filled cell ! calculate sx @@ -1857,12 +1878,14 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) else ! interior if (ISS%hmask(i+1,j) == 1) then cnt = cnt+1 + 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) then cnt = cnt+1 + Dx =dxh+ G%dxT(i-1,j) sx = sx - S(i-1,j) else sx = sx - S(i,j) @@ -1870,7 +1893,7 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) if (cnt == 0) then sx = 0 else - sx = sx / (cnt * dxh) + sx = sx / ( Dx) endif endif @@ -1892,6 +1915,7 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) else ! interior if (ISS%hmask(i,j+1) == 1) then cnt = cnt+1 + Dy =dyh+ G%dyT(i,j+1) sy = S(i,j+1) else sy = S(i,j) @@ -1899,13 +1923,14 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) if (ISS%hmask(i,j-1) == 1) then cnt = cnt+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 / (cnt * dyh) + sy = sy / (Dy) endif endif @@ -1930,10 +1955,10 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) taudy(I,J) = taudy(I,J) - .25 * rho * grav * ISS%h_shelf(i,j) * sy * G%areaT(i,j) endif if (CS%ground_frac(i,j) == 1) then -! 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)) neumann_val = .5 * grav * (1-rho/rhow) * rho * ISS%h_shelf(i,j)**2 else - neumann_val = .5 * grav * (1-rho/rhow) * rho * ISS%h_shelf(i,j)**2 + neumann_val = (.5 * grav * (1-rho/rhow) * rho * ISS%h_shelf(i,j)**2) endif if ((CS%u_face_mask(I-1,j) == 2) .OR. (ISS%hmask(i-1,j) == 0) .OR. (ISS%hmask(i-1,j) == 2) ) then @@ -1971,7 +1996,6 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) endif enddo enddo - end subroutine calc_shelf_driving_stress subroutine init_boundary_values(CS, G, time, hmask, input_flux, input_thick, new_sim) @@ -2528,8 +2552,8 @@ subroutine apply_boundary_values(CS, ISS, G, US, time, Phisub, H_node, ice_visc, 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, and also update the -!! nonlinear part of the basal traction. 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 type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe @@ -2540,9 +2564,9 @@ 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 [m-1]. ! update DEPTH_INTEGRATED viscosity, based on horizontal strain rates - this is for bilinear FEM solve -! so there is an "upper" and "lower" bilinear viscosity ! also this subroutine updates the nonlinear part of the basal traction @@ -2553,7 +2577,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, dxh, dyh ! Velocity shears [T-1 ~> s-1] - real, dimension(8,4) :: Phi +! real, dimension(8,4) :: Phi real, dimension(2) :: xquad ! real :: umid, vmid, unorm ! Velocities [L T-1 ~> m s-1] @@ -2566,6 +2590,12 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) is = iscq - 1; js = jscq - 1 i_off = G%idg_offset ; j_off = G%jdg_offset + allocate(Phi(1:8,1:4,isd:ied,jsd:jed), source=0.0) + + do j=jsc,jec ; do i=isc,iec + call bilinear_shape_fn_grid(G, i, j, Phi(:,:,i,j)) + enddo ; enddo + n_g = CS%n_glen; eps_min = CS%eps_glen_min CS%ice_visc(:,:)=1e22 ! Visc_coef = US%kg_m2s_to_RZ_T*US%m_to_L*US%Z_to_L*(CS%A_glen_isothermal)**(-1./CS%n_glen) @@ -2575,21 +2605,35 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) if ((ISS%hmask(i,j) == 1) .OR. (ISS%hmask(i,j) == 3)) then Visc_coef = US%kg_m2s_to_RZ_T*US%m_to_L*US%Z_to_L*(CS%AGlen_visc(i,j))**(-1./CS%n_glen) - ux = ((u_shlf(I,J) + (u_shlf(I,J-1) + u_shlf(I,J+1))) - & - (u_shlf(I-1,J) + (u_shlf(I-1,J-1) + u_shlf(I-1,J+1)))) / (3*G%dxT(i,j)) - vx = ((v_shlf(I,J) + v_shlf(I,J-1) + v_shlf(I,J+1)) - & - (v_shlf(I-1,J) + (v_shlf(I-1,J-1) + v_shlf(I-1,J+1)))) / (3*G%dxT(i,j)) - uy = ((u_shlf(I,J) + (u_shlf(I-1,J) + u_shlf(I+1,J))) - & - (u_shlf(I,J-1) + (u_shlf(I-1,J-1) + u_shlf(I+1,J-1)))) / (3*G%dyT(i,j)) - vy = ((v_shlf(I,J) + (v_shlf(I-1,J)+ v_shlf(I+1,J))) - & - (v_shlf(I,J-1) + (v_shlf(I-1,J-1)+ v_shlf(I+1,J-1)))) / (3*G%dyT(i,j)) + 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) + enddo ; enddo ! CS%ice_visc(i,j) =1e15*(G%areaT(i,j) * ISS%h_shelf(i,j)) ! constant viscocity for debugging 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)) endif enddo enddo - + deallocate(Phi) end subroutine calc_shelf_visc subroutine calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) diff --git a/src/ice_shelf/MOM_ice_shelf_initialize.F90 b/src/ice_shelf/MOM_ice_shelf_initialize.F90 index 469cba39ce..7cc3c020a3 100644 --- a/src/ice_shelf/MOM_ice_shelf_initialize.F90 +++ b/src/ice_shelf/MOM_ice_shelf_initialize.F90 @@ -149,14 +149,13 @@ subroutine initialize_ice_thickness_from_file(h_shelf, area_shelf_h, hmask, G, U enddo ; enddo endif - if (len_sidestress > 0.) then do j=jsc,jec do i=isc,iec ! taper ice shelf in area where there is no sidestress - ! but do not interfere with hmask - if (G%geoLonCv(i,j) > len_sidestress) then + if ((len_sidestress > 0.) .and. (G%geoLonCv(i,j) > len_sidestress)) then udh = exp(-(G%geoLonCv(i,j)-len_sidestress)/5.0) * h_shelf(i,j) if (udh <= 25.0) then h_shelf(i,j) = 0.0 @@ -180,7 +179,6 @@ subroutine initialize_ice_thickness_from_file(h_shelf, area_shelf_h, hmask, G, U endif enddo enddo - endif end subroutine initialize_ice_thickness_from_file !> Initialize ice shelf thickness for a channel configuration @@ -397,13 +395,13 @@ 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,& - hmask,h_shelf, G, US, PF) -!subroutine initialize_ice_flow_from_file(bed_elev,u_shelf, v_shelf,ice_visc,float_cond,& +!subroutine initialize_ice_flow_from_file(bed_elev,u_shelf, v_shelf,float_cond,& ! hmask,h_shelf, G, US, PF) +subroutine initialize_ice_flow_from_file(bed_elev,u_shelf, v_shelf,float_cond,& + 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 ice shelf u velocity [Z ~> m]. + intent(inout) :: bed_elev !< The bed elevation [Z ~> m]. real, dimension(SZIB_(G),SZJB_(G)), & intent(inout) :: u_shelf !< The zonal ice shelf velocity [L T-1 ~> m s-1]. real, dimension(SZIB_(G),SZJB_(G)), & @@ -412,12 +410,12 @@ subroutine initialize_ice_flow_from_file(bed_elev,u_shelf, v_shelf,float_cond,& real, dimension(SZDI_(G),SZDJ_(G)), & intent(inout) :: float_cond !< An array indicating where the ice !! shelf is floating: 0 if floating, 1 if not. - real, dimension(SZDI_(G),SZDJ_(G)), & - intent(in) :: hmask !< A mask indicating which tracer points are - !! partly or fully covered by an ice-shelf - real, dimension(SZDI_(G),SZDJ_(G)), & - intent(in) :: h_shelf !< A mask indicating which tracer points are +! real, dimension(SZDI_(G),SZDJ_(G)), & +! intent(in) :: hmask !< A mask indicating which tracer points are !! partly or fully covered by an ice-shelf +! real, dimension(SZDI_(G),SZDJ_(G)), & +! intent(in) :: h_shelf !< A mask indicating which tracer points are +! !! partly or fully covered by an ice-shelf type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors type(param_file_type), intent(in) :: PF !< A structure to parse for run-time parameters @@ -453,10 +451,10 @@ subroutine initialize_ice_flow_from_file(bed_elev,u_shelf, v_shelf,float_cond,& "The name of the thickness variable in ICE_VELOCITY_FILE.", & default="viscosity") call get_param(PF, mdl, "BED_TOPO_FILE", bed_topo_file, & - "The file from which the velocity is read.", & + "The file from which the bed elevation is read.", & default="ice_shelf_vel.nc") call get_param(PF, mdl, "BED_TOPO_VARNAME", bed_varname, & - "The name of the thickness variable in ICE_VELOCITY_FILE.", & + "The name of the thickness variable in ICE_INPUT_FILE.", & default="depth") if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, & " initialize_ice_shelf_velocity_from_file: Unable to open "//trim(filename)) @@ -470,15 +468,8 @@ subroutine initialize_ice_flow_from_file(bed_elev,u_shelf, v_shelf,float_cond,& filename = trim(inputdir)//trim(bed_topo_file) call MOM_read_data(filename,trim(bed_varname), bed_elev, G%Domain, scale=1.) - isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec +! isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec -! do j=jsc,jec -! do i=isc,iec -! if (hmask(i,j) == 1.) then -! ice_visc(i,j) = ice_visc(i,j) * (G%areaT(i,j) * h_shelf(i,j)) -! endif -! enddo -! enddo end subroutine initialize_ice_flow_from_file @@ -510,7 +501,7 @@ subroutine initialize_ice_shelf_boundary_from_file(u_face_mask_bdry, v_face_mask intent(inout) :: hmask !< A mask indicating which tracer points are !! partly or fully covered by an ice-shelf real, dimension(SZDI_(G),SZDJ_(G)), & - intent(inout) :: h_shelf !< Ice-shelf thickness + intent(in) :: h_shelf !< Ice-shelf thickness type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors type(param_file_type), intent(in) :: PF !< A structure to parse for run-time parameters @@ -535,9 +526,9 @@ subroutine initialize_ice_shelf_boundary_from_file(u_face_mask_bdry, v_face_mask call get_param(PF, mdl, "ICE_THICKNESS_FILE", icethick_file, & "The file from which the ice-shelf thickness is read.", & default="ice_shelf_thick.nc") - call get_param(PF, mdl, "ICE_THICKNESS_VARNAME", h_varname, & - "The name of the thickness variable in ICE_THICKNESS_FILE.", & - default="h_shelf") +! call get_param(PF, mdl, "ICE_THICKNESS_VARNAME", h_varname, & +! "The name of the thickness variable in ICE_THICKNESS_FILE.", & +! default="h_shelf") call get_param(PF, mdl, "ICE_THICKNESS_MASK_VARNAME", hmsk_varname, & "The name of the icethickness mask variable in ICE_THICKNESS_FILE.", & default="h_mask") @@ -574,7 +565,7 @@ subroutine initialize_ice_shelf_boundary_from_file(u_face_mask_bdry, v_face_mask call MOM_read_data(filename,trim(vmask_varname), vmask, G%Domain, position=CORNER,scale=1.) filename = trim(inputdir)//trim(icethick_file) - call MOM_read_data(filename, trim(h_varname), h_shelf, G%Domain, scale=US%m_to_Z) +! call MOM_read_data(filename, trim(h_varname), h_shelf, G%Domain, scale=US%m_to_Z) call MOM_read_data(filename,trim(hmsk_varname), hmask, G%Domain, scale=1.) isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 2aab378b4a..d7cc3b6c2d 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -117,7 +117,7 @@ module MOM_state_initialization !! conditions or by reading them from a restart (or saves) file. subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & restart_CS, ALE_CSp, tracer_Reg, sponge_CSp, & - ALE_sponge_CSp, oda_incupd_CSp, OBC, Time_in, frac_shelf_h) + ALE_sponge_CSp, oda_incupd_CSp, OBC, Time_in, frac_shelf_h, mass_shelf) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -147,6 +147,9 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & real, dimension(SZI_(G),SZJ_(G)), & optional, intent(in) :: frac_shelf_h !< The fraction of the grid cell covered !! by a floating ice shelf [nondim]. + real, dimension(SZI_(G),SZJ_(G)), & + optional, intent(in) :: mass_shelf !< The mass per unit area of the overlying + !! ice shelf [ RZ ~> kg m-2 ] ! Local variables real :: depth_tot(SZI_(G),SZJ_(G)) ! The nominal total depth of the ocean [Z ~> m] character(len=200) :: filename ! The name of an input file. @@ -158,6 +161,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & real :: vel_rescale ! A rescaling factor for velocities from the representation in ! a restart file to the internal representation in this run. real :: dt ! The baroclinic dynamics timestep for this run [T ~> s]. + logical :: from_Z_file, useALE logical :: new_sim integer :: write_geom @@ -404,6 +408,23 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & if (use_temperature .and. use_OBC) & call fill_temp_salt_segments(G, GV, OBC, tv) + ! Calculate the initial surface displacement under ice shelf + + if (new_sim) then + call get_param(PF, mdl, "DEPRESS_INITIAL_SURFACE", depress_sfc, & + "If true, depress the initial surface to avoid huge "//& + "tsunamis when a large surface pressure is applied.", & + default=.false., do_not_log=just_read) + call get_param(PF, mdl, "TRIM_IC_FOR_P_SURF", trim_ic_for_p_surf, & + "If true, cuts way the top of the column for initial conditions "//& + "at the depth where the hydrostatic pressure matches the imposed "//& + "surface pressure which is read from file.", default=.false., & + do_not_log=just_read) + + if (use_ice_shelf .and. present(mass_shelf) .and. .not. (trim_ic_for_p_surf .or. depress_sfc)) & + call calc_sfc_displacement(PF, G, GV, US, mass_shelf, tv, h) + endif + ! The thicknesses in halo points might be needed to initialize the velocities. if (new_sim) call pass_var(h, G%Domain) @@ -458,15 +479,6 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & call convert_thickness(h, G, GV, US, tv) ! Remove the mass that would be displaced by an ice shelf or inverse barometer. - call get_param(PF, mdl, "DEPRESS_INITIAL_SURFACE", depress_sfc, & - "If true, depress the initial surface to avoid huge "//& - "tsunamis when a large surface pressure is applied.", & - default=.false., do_not_log=just_read) - call get_param(PF, mdl, "TRIM_IC_FOR_P_SURF", trim_ic_for_p_surf, & - "If true, cuts way the top of the column for initial conditions "//& - "at the depth where the hydrostatic pressure matches the imposed "//& - "surface pressure which is read from file.", default=.false., & - do_not_log=just_read) if (depress_sfc .and. trim_ic_for_p_surf) call MOM_error(FATAL, "MOM_initialize_state: "//& "DEPRESS_INITIAL_SURFACE and TRIM_IC_FOR_P_SURF are exclusive and cannot both be True") if (new_sim .and. debug .and. (depress_sfc .or. trim_ic_for_p_surf)) & @@ -1035,7 +1047,7 @@ subroutine convert_thickness(h, G, GV, US, tv) end subroutine convert_thickness !> Depress the sea-surface based on an initial condition file -subroutine depress_surface(h, G, GV, US, param_file, tv, just_read) +subroutine depress_surface(h, G, GV, US, param_file, tv, just_read, z_top_shelf) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -1045,6 +1057,8 @@ subroutine depress_surface(h, G, GV, US, param_file, tv, just_read) type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables logical, intent(in) :: just_read !< If true, this call will only read !! parameters without changing h. + real, dimension(SZI_(G),SZJ_(G)), & + optional, intent(in) :: z_top_shelf !< Top interface position under ice shelf [Z ~> m] ! Local variables real, dimension(SZI_(G),SZJ_(G)) :: & eta_sfc ! The free surface height that the model should use [Z ~> m]. @@ -1057,30 +1071,40 @@ subroutine depress_surface(h, G, GV, US, param_file, tv, just_read) character(len=200) :: inputdir, eta_srf_file ! Strings for file/path character(len=200) :: filename, eta_srf_var ! Strings for file/path integer :: i, j, k, is, ie, js, je, nz + logical :: use_z_shelf is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke - ! Read the surface height (or pressure) from a file. + use_z_shelf = present(z_top_shelf) - call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") - inputdir = slasher(inputdir) - call get_param(param_file, mdl, "SURFACE_HEIGHT_IC_FILE", eta_srf_file,& - "The initial condition file for the surface height.", & - fail_if_missing=.not.just_read, do_not_log=just_read) - call get_param(param_file, mdl, "SURFACE_HEIGHT_IC_VAR", eta_srf_var, & - "The initial condition variable for the surface height.",& - default="SSH", do_not_log=just_read) - filename = trim(inputdir)//trim(eta_srf_file) - if (.not.just_read) & - call log_param(param_file, mdl, "INPUTDIR/SURFACE_HEIGHT_IC_FILE", filename) - call get_param(param_file, mdl, "SURFACE_HEIGHT_IC_SCALE", scale_factor, & - "A scaling factor to convert SURFACE_HEIGHT_IC_VAR into units of m", & - units="variable", default=1.0, scale=US%m_to_Z, do_not_log=just_read) + if (.not. use_z_shelf) then + ! Read the surface height (or pressure) from a file. - if (just_read) return ! All run-time parameters have been read, so return. + call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") + inputdir = slasher(inputdir) + call get_param(param_file, mdl, "SURFACE_HEIGHT_IC_FILE", eta_srf_file,& + "The initial condition file for the surface height.", & + fail_if_missing=.not.just_read, do_not_log=just_read) + call get_param(param_file, mdl, "SURFACE_HEIGHT_IC_VAR", eta_srf_var, & + "The initial condition variable for the surface height.",& + default="SSH", do_not_log=just_read) + filename = trim(inputdir)//trim(eta_srf_file) + if (.not.just_read) & + call log_param(param_file, mdl, "INPUTDIR/SURFACE_HEIGHT_IC_FILE", filename) + + call get_param(param_file, mdl, "SURFACE_HEIGHT_IC_SCALE", scale_factor, & + "A scaling factor to convert SURFACE_HEIGHT_IC_VAR into units of m", & + units="variable", default=1.0, scale=US%m_to_Z, do_not_log=just_read) - call MOM_read_data(filename, eta_srf_var, eta_sfc, G%Domain, scale=scale_factor) + if (just_read) return ! All run-time parameters have been read, so return. + + call MOM_read_data(filename, eta_srf_var, eta_sfc, G%Domain, scale=scale_factor) + else + do j=js,je ; do i=is,ie + eta_sfc(i,j) = z_top_shelf(i,j) + enddo; enddo + endif ! Convert thicknesses to interface heights. call find_eta(h, tv, G, GV, US, eta, dZref=G%Z_ref) @@ -1201,6 +1225,84 @@ subroutine trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read) end subroutine trim_for_ice +!> Calculate the hydrostatic equilibrium position of the surface under an ice shelf +subroutine calc_sfc_displacement(PF, G, GV, US, mass_shelf, tv, h) + type(param_file_type), intent(in) :: PF !< Parameter file structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: mass_shelf !< Ice shelf mass [RZ ~> kg m-2] + type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamics structure + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(inout) :: h !< Layer thickness [H ~> m or kg m-2] + + real :: z_top_shelf(SZI_(G),SZJ_(G)) ! The depth of the top interface under ice shelves [Z ~> m] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: & + eta ! The free surface height that the model should use [Z ~> m]. + ! temporary arrays + real, dimension(SZK_(GV)) :: rho_col ! potential density in the column for use in ice + real, dimension(SZK_(GV)) :: rho_h ! potential density multiplied by thickness [RZ ~> kg m-2 ] + real, dimension(SZK_(GV)) :: h_tmp ! temporary storage for thicknesses [H ~> m] + real, dimension(SZK_(GV)) :: p_ref ! pressure for density [RZ ~> kg m-2] + real, dimension(SZK_(GV)+1) :: ei_tmp, ei_orig ! temporary storage for interface positions [Z ~> m] + real :: z_top, z_col, mass_disp, residual, tol + integer :: is, ie, js, je, k, nz, i, j, max_iter, iter + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + + tol = US%m_to_Z*0.001 ! 1mm tolerance + max_iter = 1e3 + + call MOM_mesg("Started calculating initial interface position under ice shelf ") + ! Convert thicknesses to interface heights. + call find_eta(h, tv, G, GV, US, eta, dZref=G%Z_ref) + do j = js, je ; do i = is, ie + iter=1 + z_top_shelf(i,j) = 0.0 + p_ref(:) = tv%p_ref + if (G%mask2dT(i,j) > 0. .and. mass_shelf(i,j) .gt. 0.) then + call calculate_density(tv%T(i,j,:), tv%S(i,j,:), P_Ref, rho_col, tv%eqn_of_state) + z_top = min(max(-1.0*mass_shelf(i,j)/rho_col(1),-G%bathyT(i,j)),0.) + h_tmp = 0.0 + z_col = 0.0 + ei_tmp(1:nz+1)=eta(i,j,1:nz+1) + ei_orig(1:nz+1)=eta(i,j,1:nz+1) + do k=1,nz+1 + if (ei_tmp(k)tol .and. z_top>-G%bathyT(i,j) .and. iter .lt. max_iter) + z_top=min(max(z_top-(residual*0.5e-3),-G%bathyT(i,j)),0.0) + h_tmp = 0.0 + z_col = 0.0 + ei_tmp(1:nz+1)=ei_orig(1:nz+1) + do k=1,nz+1 + if (ei_tmp(k) Adjust the layer thicknesses by removing the top of the water column above the !! depth where the hydrostatic pressure matches p_surf @@ -2597,6 +2699,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just old_remap=remap_old_alg, answers_2018=answers_2018 ) call ALE_remap_scalar(remapCS, G, GV, nkd, h1, tmpS1dIn, h, tv%S, all_cells=remap_full_column, & old_remap=remap_old_alg, answers_2018=answers_2018 ) + deallocate( h1 ) deallocate( tmpT1dIn ) deallocate( tmpS1dIn ) From 4e349bffdda2f9426b730fb13623dd97235879bd Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Wed, 29 Dec 2021 12:48:46 -0500 Subject: [PATCH 2/4] Addresses a number of issues identified in code review. --- src/core/MOM.F90 | 6 +-- .../MOM_state_initialization.F90 | 46 ++++++++++--------- 2 files changed, 27 insertions(+), 25 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 93093f4eb2..dae65ad6da 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1732,7 +1732,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & real, allocatable, target :: frac_shelf_in(:,:) ! Initial fraction of the total cell area occupied ! by an ice shelf [nondim] real, allocatable, target :: mass_shelf_in(:,:) ! Initial mass of ice shelf contained within a grid cell - ! [RZ ~> kg m-2] + ! [R Z ~> kg m-2] real, allocatable, target :: T_in(:,:,:) ! Initial temperatures [degC] real, allocatable, target :: S_in(:,:,:) ! Initial salinities [ppt] @@ -2504,7 +2504,6 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & allocate(CS%frac_shelf_h(isd:ied, jsd:jed), source=0.0) allocate(CS%mass_shelf(isd:ied, jsd:jed), source=0.0) call ice_shelf_query(ice_shelf_CSp,G,CS%frac_shelf_h, CS%mass_shelf) - CS%mass_shelf = CS%mass_shelf*US%kg_m3_to_R*US%m_to_Z ! MOM_initialize_state is using the unrotated metric call rotate_array(CS%frac_shelf_h, -turns, frac_shelf_in) call rotate_array(CS%mass_shelf, -turns, mass_shelf_in) @@ -2562,7 +2561,6 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & allocate(CS%frac_shelf_h(isd:ied, jsd:jed), source=0.0) allocate(CS%mass_shelf(isd:ied, jsd:jed), source=0.0) call ice_shelf_query(ice_shelf_CSp,G,CS%frac_shelf_h, CS%mass_shelf) - CS%mass_shelf = CS%mass_shelf*US%kg_m3_to_R*US%m_to_Z call MOM_initialize_state(CS%u, CS%v, CS%h, CS%tv, Time, G, GV, US, & param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, & CS%sponge_CSp, CS%ALE_sponge_CSp,CS%oda_incupd_CSp, CS%OBC, Time_in, & @@ -2583,7 +2581,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & if (use_ice_shelf .and. CS%debug) then call hchksum(CS%frac_shelf_h, "MOM:frac_shelf_h", G%HI, haloshift=0) - call hchksum(CS%mass_shelf, "MOM:mass_shelf", G%HI, haloshift=0) + call hchksum(CS%mass_shelf, "MOM:mass_shelf", G%HI, haloshift=0,scale=US%RZ_to_kg_m2) endif call cpu_clock_end(id_clock_MOM_init) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 08e8d43130..b2fabf19d1 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -410,17 +410,17 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & ! Calculate the initial surface displacement under ice shelf - if (new_sim) then - call get_param(PF, mdl, "DEPRESS_INITIAL_SURFACE", depress_sfc, & - "If true, depress the initial surface to avoid huge "//& - "tsunamis when a large surface pressure is applied.", & - default=.false., do_not_log=just_read) - call get_param(PF, mdl, "TRIM_IC_FOR_P_SURF", trim_ic_for_p_surf, & - "If true, cuts way the top of the column for initial conditions "//& - "at the depth where the hydrostatic pressure matches the imposed "//& - "surface pressure which is read from file.", default=.false., & - do_not_log=just_read) + call get_param(PF, mdl, "DEPRESS_INITIAL_SURFACE", depress_sfc, & + "If true, depress the initial surface to avoid huge "//& + "tsunamis when a large surface pressure is applied.", & + default=.false., do_not_log=just_read) + call get_param(PF, mdl, "TRIM_IC_FOR_P_SURF", trim_ic_for_p_surf, & + "If true, cuts way the top of the column for initial conditions "//& + "at the depth where the hydrostatic pressure matches the imposed "//& + "surface pressure which is read from file.", default=.false., & + do_not_log=just_read) + if (new_sim) then if (use_ice_shelf .and. present(mass_shelf) .and. .not. (trim_ic_for_p_surf .or. depress_sfc)) & call calc_sfc_displacement(PF, G, GV, US, mass_shelf, tv, h) endif @@ -1251,17 +1251,21 @@ subroutine calc_sfc_displacement(PF, G, GV, US, mass_shelf, tv, h) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke - tol = US%m_to_Z*0.001 ! 1mm tolerance - max_iter = 1e3 + tol = 0.001 ! The initialization tolerance for ice shelf initialization (m) + call get_param(PF, mdl, "ICE_SHELF_INITIALIZATION_Z_TOLERANCE", tol, & + "A initialization tolerance for the calculation of the static "// & + "ice shelf displacement (m) using initial temperature and salinity profile.",& + default=tol, units="m", scale=US%m_to_Z) + max_iter = 1e3 call MOM_mesg("Started calculating initial interface position under ice shelf ") ! Convert thicknesses to interface heights. call find_eta(h, tv, G, GV, US, eta, dZref=G%Z_ref) do j = js, je ; do i = is, ie - iter=1 + iter = 1 z_top_shelf(i,j) = 0.0 p_ref(:) = tv%p_ref - if (G%mask2dT(i,j) > 0. .and. mass_shelf(i,j) .gt. 0.) then + if (G%mask2dT(i,j) .gt. 0. .and. mass_shelf(i,j) .gt. 0.) then call calculate_density(tv%T(i,j,:), tv%S(i,j,:), P_Ref, rho_col, tv%eqn_of_state) z_top = min(max(-1.0*mass_shelf(i,j)/rho_col(1),-G%bathyT(i,j)),0.) h_tmp = 0.0 @@ -1271,29 +1275,29 @@ subroutine calc_sfc_displacement(PF, G, GV, US, mass_shelf, tv, h) do k=1,nz+1 if (ei_tmp(k)tol .and. z_top>-G%bathyT(i,j) .and. iter .lt. max_iter) + do while (abs(residual) .gt. tol .and. z_top .gt. -G%bathyT(i,j) .and. iter .lt. max_iter) z_top=min(max(z_top-(residual*0.5e-3),-G%bathyT(i,j)),0.0) h_tmp = 0.0 z_col = 0.0 - ei_tmp(1:nz+1)=ei_orig(1:nz+1) + ei_tmp(1:nz+1) = ei_orig(1:nz+1) do k=1,nz+1 if (ei_tmp(k) Date: Mon, 3 Jan 2022 11:16:29 -0500 Subject: [PATCH 3/4] add appropriate intent to ice_shelf_query --- src/ice_shelf/MOM_ice_shelf.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 88f51571cd..e0cca05aa5 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -2035,8 +2035,8 @@ end subroutine update_shelf_mass subroutine ice_shelf_query(CS, G, frac_shelf_h, mass_shelf) type(ice_shelf_CS), pointer :: CS !< ice shelf control structure type(ocean_grid_type), intent(in) :: G !< A pointer to an ocean grid control structure. - real, optional, dimension(SZI_(G),SZJ_(G)) :: frac_shelf_h !< Ice shelf area fraction [nodim]. - real, optional, dimension(SZI_(G),SZJ_(G)) :: mass_shelf ! kg m-2] + real, optional, dimension(SZI_(G),SZJ_(G)), intent(out) :: frac_shelf_h !< Ice shelf area fraction [nodim]. + real, optional, dimension(SZI_(G),SZJ_(G)), intent(out) :: mass_shelf ! kg m-2] integer :: i, j From 15f98ae76b33e693832d410a8b98cebbf618f093 Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Mon, 3 Jan 2022 14:50:14 -0500 Subject: [PATCH 4/4] fix unit scaling comments --- src/ice_shelf/MOM_ice_shelf.F90 | 2 +- src/initialization/MOM_state_initialization.F90 | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index e0cca05aa5..13af5a936a 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -2036,7 +2036,7 @@ subroutine ice_shelf_query(CS, G, frac_shelf_h, mass_shelf) type(ice_shelf_CS), pointer :: CS !< ice shelf control structure type(ocean_grid_type), intent(in) :: G !< A pointer to an ocean grid control structure. real, optional, dimension(SZI_(G),SZJ_(G)), intent(out) :: frac_shelf_h !< Ice shelf area fraction [nodim]. - real, optional, dimension(SZI_(G),SZJ_(G)), intent(out) :: mass_shelf ! kg m-2] + real, optional, dimension(SZI_(G),SZJ_(G)), intent(out) :: mass_shelf ! kg m-2] integer :: i, j diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 1e6791c3d5..f95192f5f8 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -149,7 +149,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & !! by a floating ice shelf [nondim]. real, dimension(SZI_(G),SZJ_(G)), & optional, intent(in) :: mass_shelf !< The mass per unit area of the overlying - !! ice shelf [ RZ ~> kg m-2 ] + !! ice shelf [ R Z ~> kg m-2 ] ! Local variables real :: depth_tot(SZI_(G),SZJ_(G)) ! The nominal total depth of the ocean [Z ~> m] character(len=200) :: filename ! The name of an input file. @@ -1232,7 +1232,7 @@ subroutine calc_sfc_displacement(PF, G, GV, US, mass_shelf, tv, h) type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G)), & - intent(in) :: mass_shelf !< Ice shelf mass [RZ ~> kg m-2] + intent(in) :: mass_shelf !< Ice shelf mass [R Z ~> kg m-2] type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamics structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(inout) :: h !< Layer thickness [H ~> m or kg m-2] @@ -1242,9 +1242,9 @@ subroutine calc_sfc_displacement(PF, G, GV, US, mass_shelf, tv, h) eta ! The free surface height that the model should use [Z ~> m]. ! temporary arrays real, dimension(SZK_(GV)) :: rho_col ! potential density in the column for use in ice - real, dimension(SZK_(GV)) :: rho_h ! potential density multiplied by thickness [RZ ~> kg m-2 ] + real, dimension(SZK_(GV)) :: rho_h ! potential density multiplied by thickness [R Z ~> kg m-2 ] real, dimension(SZK_(GV)) :: h_tmp ! temporary storage for thicknesses [H ~> m] - real, dimension(SZK_(GV)) :: p_ref ! pressure for density [RZ ~> kg m-2] + real, dimension(SZK_(GV)) :: p_ref ! pressure for density [R Z ~> kg m-2] real, dimension(SZK_(GV)+1) :: ei_tmp, ei_orig ! temporary storage for interface positions [Z ~> m] real :: z_top, z_col, mass_disp, residual, tol integer :: is, ie, js, je, k, nz, i, j, max_iter, iter