diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index bac5b0fce9..ac12f18c9d 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -353,8 +353,8 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) character(len=160) :: mesg ! The text of an error message integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state integer :: i, j, is, ie, js, je, ied, jed, it1, it3 - real :: vaf0, vaf0_A, vaf0_G !The previous volumes above floatation [m3] - !for all ice sheets, Antarctica only, or Greenland only [m3] + real :: vaf0, vaf0_A, vaf0_G !The previous volumes above floatation [Z L2 ~> m3] + !for all ice sheets, Antarctica only, or Greenland only [Z L2 ~> m3] if (.not. associated(CS)) call MOM_error(FATAL, "shelf_calc_flux: "// & "initialize_ice_shelf must be called before shelf_calc_flux.") @@ -874,16 +874,18 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) end subroutine shelf_calc_flux -subroutine integrate_over_ice_sheet_area(G, ISS, var, var_scale, var_out, hemisphere) +function integrate_over_ice_sheet_area(G, ISS, var, unscale, hemisphere) result(var_out) type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf. type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe the ice-shelf state real, dimension(SZI_(G),SZJ_(G)), intent(in) :: var !< Ice variable to integrate in arbitrary units [A ~> a] - real, intent(in) :: var_scale !< Dimensional scaling for variable to integrate [a A-1 ~> 1] - real, intent(out) :: var_out !< Variable integrated over the area of the ice sheet in arbitrary units [a m2] + real, intent(in) :: unscale !< Dimensional scaling for variable to integrate [a A-1 ~> 1] integer, optional, intent(in) :: hemisphere !< 0 for Antarctica only, 1 for Greenland only. Otherwise, all ice sheets + real :: var_out !< Variable integrated over the area of the ice sheet in arbitrary scaled units [A L2 ~> a m2] + + ! Local variables integer :: IS_ID ! local copy of hemisphere real, dimension(SZI_(G),SZJ_(G)) :: var_cell !< Variable integrated over the ice-sheet area of each cell - !! in arbitrary units [a m2] + !! in arbitrary units [A L2 ~> a m2] integer, dimension(SZI_(G),SZJ_(G)) :: mask ! a mask for active cells depending on hemisphere indicated integer :: i,j @@ -903,16 +905,16 @@ subroutine integrate_over_ice_sheet_area(G, ISS, var, var_scale, var_out, hemisp if (ISS%hmask(i,j)>0 .and. G%geoLatT(i,j)>0.0) mask(i,j)=1 enddo; enddo else !All ice sheets - mask(G%isc:G%iec,G%jsc:G%jec)=ISS%hmask(G%isc:G%iec,G%jsc:G%jec) + mask(G%isc:G%iec,G%jsc:G%jec) = ISS%hmask(G%isc:G%iec,G%jsc:G%jec) endif var_cell(:,:)=0.0 do j = G%jsc,G%jec; do i = G%isc,G%iec - if (mask(i,j)>0) var_cell(i,j) = (var(i,j) * var_scale) * (ISS%area_shelf_h(i,j) * G%US%L_to_m**2) + if (mask(i,j)>0) var_cell(i,j) = var(i,j) * ISS%area_shelf_h(i,j) enddo; enddo - var_out = reproducing_sum(var_cell) -end subroutine integrate_over_ice_sheet_area + var_out = reproducing_sum(var_cell, unscale=unscale*G%US%L_to_m**2) +end function integrate_over_ice_sheet_area !> Converts the ice-shelf-to-ocean calving and calving_hflx variables from the ice-shelf state (ISS) type !! to the ocean public type @@ -1252,10 +1254,8 @@ subroutine add_shelf_flux(G, US, CS, sfc_state, fluxes, time_step) do j=js,je ; do i=is,ie last_hmask(i,j) = ISS%hmask(i,j) ; last_area_shelf_h(i,j) = ISS%area_shelf_h(i,j) enddo ; enddo - call time_interp_external(CS%mass_handle, Time0, last_mass_shelf) + call time_interp_external(CS%mass_handle, Time0, last_mass_shelf, scale=US%kg_m3_to_R*US%m_to_Z) do j=js,je ; do i=is,ie - ! This should only be done if time_interp_extern did an update. - last_mass_shelf(i,j) = US%kg_m3_to_R*US%m_to_Z * last_mass_shelf(i,j) ! Rescale after time_interp last_h_shelf(i,j) = last_mass_shelf(i,j) / CS%density_ice enddo ; enddo @@ -1853,7 +1853,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init, v_desc = var_desc("tauy_shelf", "Pa", "the meridional stress on the ocean under ice shelves", & hor_grid='Cv',z_grid='1') call register_restart_pair(sfc_state%taux_shelf, sfc_state%tauy_shelf, u_desc, v_desc, & - .false., CS%restart_CSp, conversion=US%RZ_T_to_kg_m2s*US%L_T_to_m_s) + .false., CS%restart_CSp, conversion=US%RLZ_T2_to_Pa) endif endif @@ -1997,134 +1997,161 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init, 'Fric vel under shelf', 'm/s', conversion=US%Z_to_m*US%s_to_T) if (CS%active_shelf_dynamics) then CS%id_h_mask = register_diag_field('ice_shelf_model', 'h_mask', CS%diag%axesT1, CS%Time, & - 'ice shelf thickness mask', 'none') + 'ice shelf thickness mask', 'none', conversion=1.0) CS%id_shelf_sfc_mass_flux = register_diag_field('ice_shelf_model', 'sfc_mass_flux', CS%diag%axesT1, CS%Time, & 'ice shelf surface mass flux deposition from atmosphere', & 'kg m-2 s-1', conversion=US%RZ_T_to_kg_m2s) endif - !scalars (area integrated over all ice sheets) + ! Scalars (area integrated over all ice sheets) CS%id_vaf = register_scalar_field('ice_shelf_model', 'int_vaf', CS%diag%axesT1, CS%Time, & - 'Area integrated ice sheet volume above floatation', 'm3') + 'Area integrated ice sheet volume above floatation', 'm3', conversion=US%Z_to_m*US%L_to_m**2) CS%id_adott = register_scalar_field('ice_shelf_model', 'int_a', CS%diag%axesT1, CS%Time, & - 'Area integrated change in ice-sheet thickness ' //& - 'due to surface accum+melt during a DT_THERM time step', 'm3') + 'Area integrated change in ice-sheet thickness ' //& + 'due to surface accum+melt during a DT_THERM time step', 'm3', conversion=US%Z_to_m*US%L_to_m**2) CS%id_g_adott = register_scalar_field('ice_shelf_model', 'int_a_ground', CS%diag%axesT1, CS%Time, & - 'Area integrated change in grounded ice-sheet thickness ' //& - 'due to surface accum+melt during a DT_THERM time step', 'm3') + 'Area integrated change in grounded ice-sheet thickness ' //& + 'due to surface accum+melt during a DT_THERM time step', 'm3', conversion=US%Z_to_m*US%L_to_m**2) CS%id_f_adott = register_scalar_field('ice_shelf_model', 'int_a_float', CS%diag%axesT1, CS%Time, & - 'Area integrated change in floating ice-shelf thickness ' //& - 'due to surface accum+melt during a DT_THERM time step', 'm3') + 'Area integrated change in floating ice-shelf thickness ' //& + 'due to surface accum+melt during a DT_THERM time step', 'm3', conversion=US%Z_to_m*US%L_to_m**2) CS%id_bdott = register_scalar_field('ice_shelf_model', 'int_b', CS%diag%axesT1, CS%Time, & - 'Area integrated change in floating ice-shelf thickness '//& - 'due to basal accum+melt during a DT_THERM time step', 'm3') + 'Area integrated change in floating ice-shelf thickness '//& + 'due to basal accum+melt during a DT_THERM time step', 'm3', conversion=US%Z_to_m*US%L_to_m**2) CS%id_bdott_melt = register_scalar_field('ice_shelf_model', 'int_b_melt', CS%diag%axesT1, CS%Time, & - 'Area integrated basal melt over ice shelves during a DT_THERM time step', 'm3') + 'Area integrated basal melt over ice shelves during a DT_THERM time step', & + units='m3', conversion=US%Z_to_m*US%L_to_m**2) CS%id_bdott_accum = register_scalar_field('ice_shelf_model', 'int_b_accum', CS%diag%axesT1, CS%Time, & - 'Area integrated basal accumulation over ice shelves during a DT_THERM a time step', 'm3') + 'Area integrated basal accumulation over ice shelves during a DT_THERM a time step', & + units='m3', conversion=US%Z_to_m*US%L_to_m**2) CS%id_t_area = register_scalar_field('ice_shelf_model', 'tot_area', CS%diag%axesT1, CS%Time, & - 'Total ice-sheet area', 'm2') + 'Total ice-sheet area', 'm2', conversion=US%L_to_m**2) CS%id_f_area = register_scalar_field('ice_shelf_model', 'tot_area_float', CS%diag%axesT1, CS%Time, & - 'Total area of floating ice shelves', 'm2') + 'Total area of floating ice shelves', 'm2', conversion=US%L_to_m**2) CS%id_g_area = register_scalar_field('ice_shelf_model', 'tot_area_ground', CS%diag%axesT1, CS%Time, & - 'Total area of grounded ice sheets', 'm2') + 'Total area of grounded ice sheets', 'm2', conversion=US%L_to_m**2) !scalars (area integrated rates over all ice sheets) CS%id_dvafdt = register_scalar_field('ice_shelf_model', 'int_vafdot', CS%diag%axesT1, CS%Time, & - 'Area integrated rate of change in ice-sheet volume above floatation', 'm3 s-1') - CS%id_adot = register_scalar_field('ice_shelf_model', 'int_adot', CS%diag%axesT1, CS%Time, & - 'Area integrated rate of change in ice-sheet thickness due to surface accum+melt', 'm3 s-1') + 'Area integrated rate of change in ice-sheet volume above floatation', & + units='m3 s-1', conversion=US%Z_to_m*US%L_to_m**2*US%s_to_T) + CS%id_adot = register_scalar_field('ice_shelf_model', 'int_adot', CS%diag%axesT1, CS%Time, & + 'Area integrated rate of change in ice-sheet thickness due to surface accum+melt', & + units='m3 s-1', conversion=US%Z_to_m*US%L_to_m**2*US%s_to_T) CS%id_g_adot = register_scalar_field('ice_shelf_model', 'int_adot_ground', CS%diag%axesT1, CS%Time, & - 'Area integrated rate of change in grounded ice-sheet thickness due to surface accum+melt', 'm3 s-1') + 'Area integrated rate of change in grounded ice-sheet thickness due to surface accum+melt', & + units='m3 s-1', conversion=US%Z_to_m*US%L_to_m**2*US%s_to_T) CS%id_f_adot = register_scalar_field('ice_shelf_model', 'int_adot_float', CS%diag%axesT1, CS%Time, & - 'Area integrated rate of change in floating ice-shelf thickness due to surface accum+melt', 'm3 s-1') + 'Area integrated rate of change in floating ice-shelf thickness due to surface accum+melt', & + units='m3 s-1', conversion=US%Z_to_m*US%L_to_m**2*US%s_to_T) CS%id_bdot = register_scalar_field('ice_shelf_model', 'int_bdot', CS%diag%axesT1, CS%Time, & - 'Area integrated rate of change in ice-shelf thickness due to basal accum+melt', 'm3 s-1') + 'Area integrated rate of change in ice-shelf thickness due to basal accum+melt', & + units='m3 s-1', conversion=US%Z_to_m*US%L_to_m**2*US%s_to_T) CS%id_bdot_melt = register_scalar_field('ice_shelf_model', 'int_bdot_melt', CS%diag%axesT1, CS%Time, & - 'Area integrated basal melt rate over ice shelves', 'm3 s-1') + 'Area integrated basal melt rate over ice shelves', & + units='m3 s-1', conversion=US%Z_to_m*US%L_to_m**2*US%s_to_T) CS%id_bdot_accum = register_scalar_field('ice_shelf_model', 'int_bdot_accum', CS%diag%axesT1, CS%Time, & - 'Area integrated basal accumulation rate over ice shelves', 'm3 s-1') + 'Area integrated basal accumulation rate over ice shelves', & + units='m3 s-1', conversion=US%Z_to_m*US%L_to_m**2*US%s_to_T) !scalars (area integrated over the Antarctic ice sheet) CS%id_Ant_vaf = register_scalar_field('ice_shelf_model', 'int_vaf_A', CS%diag%axesT1, CS%Time, & - 'Area integrated Antarctic ice sheet volume above floatation', 'm3') + 'Area integrated Antarctic ice sheet volume above floatation', 'm3', conversion=US%Z_to_m*US%L_to_m**2) CS%id_Ant_adott = register_scalar_field('ice_shelf_model', 'int_a_A', CS%diag%axesT1, CS%Time, & - 'Area integrated (Antarctic ice sheet) change in ice-sheet thickness ' //& - 'due to surface accum+melt during a DT_THERM time step', 'm3') + 'Area integrated (Antarctic ice sheet) change in ice-sheet thickness ' //& + 'due to surface accum+melt during a DT_THERM time step', 'm3', conversion=US%Z_to_m*US%L_to_m**2) CS%id_Ant_g_adott = register_scalar_field('ice_shelf_model', 'int_a_ground_A', CS%diag%axesT1, CS%Time, & - 'Area integrated change in Antarctic grounded ice-sheet thickness ' //& - 'due to surface accum+melt during a DT_THERM time step', 'm3') + 'Area integrated change in Antarctic grounded ice-sheet thickness ' //& + 'due to surface accum+melt during a DT_THERM time step', 'm3', conversion=US%Z_to_m*US%L_to_m**2) CS%id_Ant_f_adott = register_scalar_field('ice_shelf_model', 'int_a_float_A', CS%diag%axesT1, CS%Time, & - 'Area integrated change in Antarctic floating ice-shelf thickness ' //& - 'due to surface accum+melt during a DT_THERM time step', 'm3') + 'Area integrated change in Antarctic floating ice-shelf thickness ' //& + 'due to surface accum+melt during a DT_THERM time step', 'm3', conversion=US%Z_to_m*US%L_to_m**2) CS%id_Ant_bdott = register_scalar_field('ice_shelf_model', 'int_b_A', CS%diag%axesT1, CS%Time, & - 'Area integrated change in Antarctic floating ice-shelf thickness '//& - 'due to basal accum+melt during a DT_THERM time step', 'm3') + 'Area integrated change in Antarctic floating ice-shelf thickness '//& + 'due to basal accum+melt during a DT_THERM time step', 'm3', conversion=US%Z_to_m*US%L_to_m**2) CS%id_Ant_bdott_melt = register_scalar_field('ice_shelf_model', 'int_b_melt_A', CS%diag%axesT1, CS%Time, & - 'Area integrated basal melt over Antarctic ice shelves during a DT_THERM time step', 'm3') + 'Area integrated basal melt over Antarctic ice shelves during a DT_THERM time step', & + units='m3', conversion=US%Z_to_m*US%L_to_m**2) CS%id_Ant_bdott_accum = register_scalar_field('ice_shelf_model', 'int_b_accum_A', CS%diag%axesT1, CS%Time, & - 'Area integrated basal accumulation over Antarctic ice shelves during a DT_THERM a time step', 'm3') + 'Area integrated basal accumulation over Antarctic ice shelves during a DT_THERM a time step', & + units='m3', conversion=US%Z_to_m*US%L_to_m**2) CS%id_Ant_t_area = register_scalar_field('ice_shelf_model', 'tot_area_A', CS%diag%axesT1, CS%Time, & - 'Total area of Antarctic ice sheet', 'm2') + 'Total area of Antarctic ice sheet', 'm2', conversion=US%L_to_m**2) CS%id_Ant_f_area = register_scalar_field('ice_shelf_model', 'tot_area_float_A', CS%diag%axesT1, CS%Time, & - 'Total area of Antarctic floating ice shelves', 'm2') + 'Total area of Antarctic floating ice shelves', 'm2', conversion=US%L_to_m**2) CS%id_Ant_g_area = register_scalar_field('ice_shelf_model', 'tot_area_ground_A', CS%diag%axesT1, CS%Time, & - 'Total area of Antarctic grounded ice sheet', 'm2') + 'Total area of Antarctic grounded ice sheet', 'm2', conversion=US%L_to_m**2) !scalars (area integrated rates over the Antarctic ice sheet) CS%id_Ant_dvafdt = register_scalar_field('ice_shelf_model', 'int_vafdot_A', CS%diag%axesT1, CS%Time, & - 'Area integrated rate of change in Antarctic ice-sheet volume above floatation', 'm3 s-1') - CS%id_Ant_adot = register_scalar_field('ice_shelf_model', 'int_adot_A', CS%diag%axesT1, CS%Time, & - 'Area integrated rate of change in Antarctic ice-sheet thickness due to surface accum+melt', 'm3 s-1') + 'Area integrated rate of change in Antarctic ice-sheet volume above floatation', & + units='m3 s-1', conversion=US%Z_to_m*US%L_to_m**2*US%s_to_T) + CS%id_Ant_adot = register_scalar_field('ice_shelf_model', 'int_adot_A', CS%diag%axesT1, CS%Time, & + 'Area integrated rate of change in Antarctic ice-sheet thickness due to surface accum+melt', & + units='m3 s-1', conversion=US%Z_to_m*US%L_to_m**2*US%s_to_T) CS%id_Ant_g_adot = register_scalar_field('ice_shelf_model', 'int_adot_ground_A', CS%diag%axesT1, CS%Time, & - 'Area integrated rate of change in Antarctic grounded ice-sheet thickness due to surface accum+melt', 'm3 s-1') + 'Area integrated rate of change in Antarctic grounded ice-sheet thickness due to surface accum+melt', & + units='m3 s-1', conversion=US%Z_to_m*US%L_to_m**2*US%s_to_T) CS%id_Ant_f_adot = register_scalar_field('ice_shelf_model', 'int_adot_float_A', CS%diag%axesT1, CS%Time, & - 'Area integrated rate of change in Antarctic floating ice-shelf thickness due to surface accum+melt', 'm3 s-1') + 'Area integrated rate of change in Antarctic floating ice-shelf thickness due to surface accum+melt', & + units='m3 s-1', conversion=US%Z_to_m*US%L_to_m**2*US%s_to_T) CS%id_Ant_bdot = register_scalar_field('ice_shelf_model', 'int_bdot_A', CS%diag%axesT1, CS%Time, & - 'Area integrated rate of change in Antarctic ice-shelf thickness due to basal accum+melt', 'm3 s-1') + 'Area integrated rate of change in Antarctic ice-shelf thickness due to basal accum+melt', & + units='m3 s-1', conversion=US%Z_to_m*US%L_to_m**2*US%s_to_T) CS%id_Ant_bdot_melt = register_scalar_field('ice_shelf_model', 'int_bdot_melt_A', CS%diag%axesT1, CS%Time, & - 'Area integrated basal melt rate over Antarctic ice shelves', 'm3 s-1') + 'Area integrated basal melt rate over Antarctic ice shelves', & + units='m3 s-1', conversion=US%Z_to_m*US%L_to_m**2*US%s_to_T) CS%id_Ant_bdot_accum = register_scalar_field('ice_shelf_model', 'int_bdot_accum_A', CS%diag%axesT1, CS%Time, & - 'Area integrated basal accumulation rate over Antarctic ice shelves', 'm3 s-1') + 'Area integrated basal accumulation rate over Antarctic ice shelves', & + units='m3 s-1', conversion=US%Z_to_m*US%L_to_m**2*US%s_to_T) !scalars (area integrated over the Greenland ice sheet) CS%id_Gr_vaf = register_scalar_field('ice_shelf_model', 'int_vaf_G', CS%diag%axesT1, CS%Time, & - 'Area integrated Greenland ice sheet volume above floatation', 'm3') + 'Area integrated Greenland ice sheet volume above floatation', 'm3', conversion=US%Z_to_m*US%L_to_m**2) CS%id_Gr_adott = register_scalar_field('ice_shelf_model', 'int_a_G', CS%diag%axesT1, CS%Time, & - 'Area integrated (Greenland ice sheet) change in ice-sheet thickness ' //& - 'due to surface accum+melt during a DT_THERM time step', 'm3') + 'Area integrated (Greenland ice sheet) change in ice-sheet thickness ' //& + 'due to surface accum+melt during a DT_THERM time step', 'm3', conversion=US%Z_to_m*US%L_to_m**2) CS%id_Gr_g_adott = register_scalar_field('ice_shelf_model', 'int_a_ground_G', CS%diag%axesT1, CS%Time, & - 'Area integrated change in Greenland grounded ice-sheet thickness ' //& - 'due to surface accum+melt during a DT_THERM time step', 'm3') + 'Area integrated change in Greenland grounded ice-sheet thickness ' //& + 'due to surface accum+melt during a DT_THERM time step', 'm3', conversion=US%Z_to_m*US%L_to_m**2) CS%id_Gr_f_adott = register_scalar_field('ice_shelf_model', 'int_a_float_G', CS%diag%axesT1, CS%Time, & - 'Area integrated change in Greenland floating ice-shelf thickness ' //& - 'due to surface accum+melt during a DT_THERM time step', 'm3') + 'Area integrated change in Greenland floating ice-shelf thickness ' //& + 'due to surface accum+melt during a DT_THERM time step', 'm3', conversion=US%Z_to_m*US%L_to_m**2) CS%id_Gr_bdott = register_scalar_field('ice_shelf_model', 'int_b_G', CS%diag%axesT1, CS%Time, & - 'Area integrated change in Greenland floating ice-shelf thickness '//& - 'due to basal accum+melt during a DT_THERM time step', 'm3') + 'Area integrated change in Greenland floating ice-shelf thickness '//& + 'due to basal accum+melt during a DT_THERM time step', 'm3', conversion=US%Z_to_m*US%L_to_m**2) CS%id_Gr_bdott_melt = register_scalar_field('ice_shelf_model', 'int_b_melt_G', CS%diag%axesT1, CS%Time, & - 'Area integrated basal melt over Greenland ice shelves during a DT_THERM time step', 'm3') + 'Area integrated basal melt over Greenland ice shelves during a DT_THERM time step', & + units='m3', conversion=US%Z_to_m*US%L_to_m**2) CS%id_Gr_bdott_accum = register_scalar_field('ice_shelf_model', 'int_b_accum_G', CS%diag%axesT1, CS%Time, & - 'Area integrated basal accumulation over Greenland ice shelves during a DT_THERM a time step', 'm3') + 'Area integrated basal accumulation over Greenland ice shelves during a DT_THERM a time step', & + units='m3', conversion=US%Z_to_m*US%L_to_m**2) CS%id_Gr_t_area = register_scalar_field('ice_shelf_model', 'tot_area_G', CS%diag%axesT1, CS%Time, & - 'Total area of Greenland ice sheet', 'm2') + 'Total area of Greenland ice sheet', 'm2', conversion=US%L_to_m**2) CS%id_Gr_f_area = register_scalar_field('ice_shelf_model', 'tot_area_float_G', CS%diag%axesT1, CS%Time, & - 'Total area of Greenland floating ice shelves', 'm2') + 'Total area of Greenland floating ice shelves', 'm2', conversion=US%L_to_m**2) CS%id_Gr_g_area = register_scalar_field('ice_shelf_model', 'tot_area_ground_G', CS%diag%axesT1, CS%Time, & - 'Total area of Greenland grounded ice sheet', 'm2') + 'Total area of Greenland grounded ice sheet', 'm2', conversion=US%L_to_m**2) !scalars (area integrated rates over the Greenland ice sheet) CS%id_Gr_dvafdt = register_scalar_field('ice_shelf_model', 'int_vafdot_G', CS%diag%axesT1, CS%Time, & - 'Area integrated rate of change in Greenland ice-sheet volume above floatation', 'm3 s-1') - CS%id_Gr_adot = register_scalar_field('ice_shelf_model', 'int_adot_G', CS%diag%axesT1, CS%Time, & - 'Area integrated rate of change in Greenland ice-sheet thickness due to surface accum+melt', 'm3 s-1') + 'Area integrated rate of change in Greenland ice-sheet volume above floatation', & + units='m3 s-1', conversion=US%Z_to_m*US%L_to_m**2*US%s_to_T) + CS%id_Gr_adot = register_scalar_field('ice_shelf_model', 'int_adot_G', CS%diag%axesT1, CS%Time, & + 'Area integrated rate of change in Greenland ice-sheet thickness due to surface accum+melt', & + units='m3 s-1', conversion=US%Z_to_m*US%L_to_m**2*US%s_to_T) CS%id_Gr_g_adot = register_scalar_field('ice_shelf_model', 'int_adot_ground_G', CS%diag%axesT1, CS%Time, & - 'Area integrated rate of change in Greenland grounded ice-sheet thickness due to surface accum+melt', 'm3 s-1') + 'Area integrated rate of change in Greenland grounded ice-sheet thickness due to surface accum+melt', & + units='m3 s-1', conversion=US%Z_to_m*US%L_to_m**2*US%s_to_T) CS%id_Gr_f_adot = register_scalar_field('ice_shelf_model', 'int_adot_float_G', CS%diag%axesT1, CS%Time, & - 'Area integrated rate of change in Greenland floating ice-shelf thickness due to surface accum+melt', 'm3 s-1') + 'Area integrated rate of change in Greenland floating ice-shelf thickness due to surface accum+melt', & + units='m3 s-1', conversion=US%Z_to_m*US%L_to_m**2*US%s_to_T) CS%id_Gr_bdot = register_scalar_field('ice_shelf_model', 'int_bdot_G', CS%diag%axesT1, CS%Time, & - 'Area integrated rate of change in Greenland ice-shelf thickness due to basal accum+melt', 'm3 s-1') + 'Area integrated rate of change in Greenland ice-shelf thickness due to basal accum+melt', & + units='m3 s-1', conversion=US%Z_to_m*US%L_to_m**2*US%s_to_T) CS%id_Gr_bdot_melt = register_scalar_field('ice_shelf_model', 'int_bdot_melt_G', CS%diag%axesT1, CS%Time, & - 'Area integrated basal melt rate over Greenland ice shelves', 'm3 s-1') + 'Area integrated basal melt rate over Greenland ice shelves', & + units='m3 s-1', conversion=US%Z_to_m*US%L_to_m**2*US%s_to_T) CS%id_Gr_bdot_accum = register_scalar_field('ice_shelf_model', 'int_bdot_accum_G', CS%diag%axesT1, CS%Time, & - 'Area integrated basal accumulation rate over Greenland ice shelves', 'm3 s-1') + 'Area integrated basal accumulation rate over Greenland ice shelves', & + units='m3 s-1', conversion=US%Z_to_m*US%L_to_m**2*US%s_to_T) !Flags to calculate diagnostics related to surface/basal mass balance if (CS%id_adott>0 .or. CS%id_g_adott>0 .or. CS%id_f_adott>0 .or. & @@ -2385,7 +2412,7 @@ subroutine update_shelf_mass(G, US, CS, ISS, Time) ! local variables integer :: i, j, is, ie, js, je - real, allocatable, dimension(:,:) :: tmp2d ! Temporary array for storing ice shelf input data + real, allocatable, dimension(:,:) :: tmp2d ! Temporary array for storing ice shelf input data [R Z ~> kg m-2] is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec @@ -2396,15 +2423,10 @@ subroutine update_shelf_mass(G, US, CS, ISS, Time) allocate(tmp2d(is:ie,js:je), source=0.0) endif - call time_interp_external(CS%mass_handle, Time, tmp2d) + call time_interp_external(CS%mass_handle, Time, tmp2d, scale=US%kg_m3_to_R*US%m_to_Z) call rotate_array(tmp2d, CS%turns, ISS%mass_shelf) deallocate(tmp2d) - ! This should only be done if time_interp_external did an update. - do j=js,je ; do i=is,ie - ISS%mass_shelf(i,j) = US%kg_m3_to_R*US%m_to_Z * ISS%mass_shelf(i,j) ! Rescale after time_interp - enddo ; enddo - do j=js,je ; do i=is,ie ISS%area_shelf_h(i,j) = 0.0 ISS%hmask(i,j) = 0. @@ -2525,7 +2547,7 @@ subroutine solo_step_ice_shelf(CS, time_interval, nsteps, Time, min_time_step_in ! coupled ice-ocean dynamics. integer :: is, ie, js, je, i, j real :: vaf0, vaf0_A, vaf0_G !The previous volumes above floatation - !for all ice sheets, Antarctica only, or Greenland only [m3] + !for all ice sheets, Antarctica only, or Greenland only [Z L2 ~> m3] real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: & dh_adott_sum, & ! Surface melt/accumulation over a full time step, used for diagnostics [Z ~> m] dh_adott ! Surface melt/accumulation over a partial time step, used for diagnostics [Z ~> m] @@ -2609,16 +2631,19 @@ end subroutine solo_step_ice_shelf !> Post_data calls for ice-sheet scalars subroutine process_and_post_scalar_data(CS, vaf0, vaf0_A, vaf0_G, Itime_step, dh_adott, dh_bdott) type(ice_shelf_CS), pointer :: CS !< A pointer to the ice shelf control structure - real :: vaf0 !< The previous volumes above floatation for all ice sheets [m3] - real :: vaf0_A !< The previous volumes above floatation for the Antarctic ice sheet [m3] - real :: vaf0_G !< The previous volumes above floatation for the Greenland ice sheet [m3] + real :: vaf0 !< The previous volumes above floatation for all ice sheets [Z L2 ~> m3] + real :: vaf0_A !< The previous volumes above floatation for the Antarctic ice sheet [Z L2 ~> m3] + real :: vaf0_G !< The previous volumes above floatation for the Greenland ice sheet [Z L2 ~> m3] real :: Itime_step !< Inverse of the time step [T-1 ~> s-1] real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: dh_adott !< Surface (plus basal if solo shelf mode) !! melt/accumulation over a time step [Z ~> m] real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: dh_bdott !< Surface (plus basal if solo shelf mode) !! melt/accumulation over a time step [Z ~> m] + + ! Local variables real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: tmp ! Temporary field used when calculating diagnostics [various] - real :: vaf ! The current ice-sheet volume above floatation [m3] + real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: ones ! Temporary field used when calculating diagnostics [various] + real :: vaf ! The current ice-sheet volume above floatation [Z L2 ~> m3] real :: val ! Temporary value when calculating scalar diagnostics [various] type(ocean_grid_type), pointer :: G => NULL() ! A pointer to the ocean's grid structure type(unit_scale_type), pointer :: US => NULL() ! Pointer to a structure containing various unit conversion factors @@ -2636,13 +2661,13 @@ subroutine process_and_post_scalar_data(CS, vaf0, vaf0_A, vaf0_G, Itime_step, dh if (CS%id_vaf > 0) call post_scalar_data(CS%id_vaf ,vaf ,CS%diag) !current vaf if (CS%id_dvafdt > 0) call post_scalar_data(CS%id_dvafdt,(vaf-vaf0)*Itime_step,CS%diag) !d(vaf)/dt if (CS%id_adott > 0 .or. CS%id_adot > 0) then !surface accumulation - surface melt - call integrate_over_ice_sheet_area(G, ISS, dh_adott, US%Z_to_m, val) + val = integrate_over_ice_sheet_area(G, ISS, dh_adott, unscale=US%Z_to_m) if (CS%id_adott > 0) call post_scalar_data(CS%id_adott,val ,CS%diag) if (CS%id_adot > 0) call post_scalar_data(CS%id_adot ,val*Itime_step,CS%diag) endif if (CS%id_g_adott > 0 .or. CS%id_g_adot > 0) then !grounded only: surface accumulation - surface melt call masked_var_grounded(G,CS%dCS,dh_adott,tmp) - call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val) + val = integrate_over_ice_sheet_area(G, ISS, tmp, unscale=US%Z_to_m) if (CS%id_g_adott > 0) call post_scalar_data(CS%id_g_adott,val ,CS%diag) if (CS%id_g_adot > 0) call post_scalar_data(CS%id_g_adot ,val*Itime_step,CS%diag) endif @@ -2651,12 +2676,12 @@ subroutine process_and_post_scalar_data(CS, vaf0, vaf0_A, vaf0_G, Itime_step, dh do j=js,je ; do i=is,ie tmp(i,j) = dh_adott(i,j) - tmp(i,j) enddo; enddo - call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val) + val = integrate_over_ice_sheet_area(G, ISS, tmp, unscale=US%Z_to_m) if (CS%id_f_adott > 0) call post_scalar_data(CS%id_f_adott,val ,CS%diag) if (CS%id_f_adot > 0) call post_scalar_data(CS%id_f_adot ,val*Itime_step,CS%diag) endif if (CS%id_bdott > 0 .or. CS%id_bdot > 0) then !bottom accumulation - bottom melt - call integrate_over_ice_sheet_area(G, ISS, dh_bdott, US%Z_to_m, val) + val = integrate_over_ice_sheet_area(G, ISS, dh_bdott, unscale=US%Z_to_m) if (CS%id_bdott > 0) call post_scalar_data(CS%id_bdott,val ,CS%diag) if (CS%id_bdot > 0) call post_scalar_data(CS%id_bdot ,val*Itime_step,CS%diag) endif @@ -2665,7 +2690,7 @@ subroutine process_and_post_scalar_data(CS, vaf0, vaf0_A, vaf0_G, Itime_step, dh do j=js,je ; do i=is,ie if (dh_bdott(i,j) < 0) tmp(i,j) = -dh_bdott(i,j) enddo; enddo - call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val) + val = integrate_over_ice_sheet_area(G, ISS, tmp, unscale=US%Z_to_m) if (CS%id_bdott_melt > 0) call post_scalar_data(CS%id_bdott_melt,val ,CS%diag) if (CS%id_bdot_melt > 0) call post_scalar_data(CS%id_bdot_melt ,val*Itime_step,CS%diag) endif @@ -2674,22 +2699,22 @@ subroutine process_and_post_scalar_data(CS, vaf0, vaf0_A, vaf0_G, Itime_step, dh do j=js,je ; do i=is,ie if (dh_bdott(i,j) > 0) tmp(i,j) = dh_bdott(i,j) enddo; enddo - call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val) + val = integrate_over_ice_sheet_area(G, ISS, tmp, unscale=US%Z_to_m) if (CS%id_bdott_accum > 0) call post_scalar_data(CS%id_bdott_accum,val ,CS%diag) if (CS%id_bdot_accum > 0) call post_scalar_data(CS%id_bdot_accum ,val*Itime_step,CS%diag) endif if (CS%id_t_area > 0) then !ice sheet area - tmp(:,:) = 1.0; call integrate_over_ice_sheet_area(G, ISS, tmp, 1.0, val) + tmp(:,:) = 1.0 ; val = integrate_over_ice_sheet_area(G, ISS, tmp, unscale=1.0) call post_scalar_data(CS%id_t_area,val,CS%diag) endif if (CS%id_g_area > 0 .or. CS%id_f_area > 0) then - tmp(:,:) = 1.0; call masked_var_grounded(G,CS%dCS,tmp,tmp) + ones(:,:) = 1.0 ; call masked_var_grounded(G, CS%dCS, ones, tmp) if (CS%id_g_area > 0) then !grounded only ice sheet area - call integrate_over_ice_sheet_area(G, ISS, tmp, 1.0, val) + val = integrate_over_ice_sheet_area(G, ISS, tmp, unscale=1.0) call post_scalar_data(CS%id_g_area,val,CS%diag) endif if (CS%id_f_area > 0) then !floating only ice sheet area (ice shelf area) - call integrate_over_ice_sheet_area(G, ISS, 1.0-tmp, 1.0, val) + val = integrate_over_ice_sheet_area(G, ISS, 1.0-tmp, unscale=1.0) call post_scalar_data(CS%id_f_area,val,CS%diag) endif endif @@ -2700,13 +2725,13 @@ subroutine process_and_post_scalar_data(CS, vaf0, vaf0_A, vaf0_G, Itime_step, dh if (CS%id_Ant_vaf > 0) call post_scalar_data(CS%id_Ant_vaf ,vaf ,CS%diag) !current vaf if (CS%id_Ant_dvafdt > 0) call post_scalar_data(CS%id_Ant_dvafdt,(vaf-vaf0_A)*Itime_step,CS%diag) !d(vaf)/dt if (CS%id_Ant_adott > 0 .or. CS%id_Ant_adot > 0) then !surface accumulation - surface melt - call integrate_over_ice_sheet_area(G, ISS, dh_adott, US%Z_to_m, val, hemisphere=0) + val = integrate_over_ice_sheet_area(G, ISS, dh_adott, unscale=US%Z_to_m, hemisphere=0) if (CS%id_Ant_adott > 0) call post_scalar_data(CS%id_Ant_adott,val ,CS%diag) if (CS%id_Ant_adot > 0) call post_scalar_data(CS%id_Ant_adot ,val*Itime_step,CS%diag) endif if (CS%id_Ant_g_adott > 0 .or. CS%id_Ant_g_adot > 0) then !grounded only: surface accumulation - surface melt call masked_var_grounded(G,CS%dCS,dh_adott,tmp) - call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val, hemisphere=0) + val = integrate_over_ice_sheet_area(G, ISS, tmp, unscale=US%Z_to_m, hemisphere=0) if (CS%id_Ant_g_adott > 0) call post_scalar_data(CS%id_Ant_g_adott,val ,CS%diag) if (CS%id_Ant_g_adot > 0) call post_scalar_data(CS%id_Ant_g_adot ,val*Itime_step,CS%diag) endif @@ -2715,12 +2740,12 @@ subroutine process_and_post_scalar_data(CS, vaf0, vaf0_A, vaf0_G, Itime_step, dh do j=js,je ; do i=is,ie tmp(i,j) = dh_adott(i,j) - tmp(i,j) enddo; enddo - call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val, hemisphere=0) + val = integrate_over_ice_sheet_area(G, ISS, tmp, unscale=US%Z_to_m, hemisphere=0) if (CS%id_Ant_f_adott > 0) call post_scalar_data(CS%id_Ant_f_adott,val ,CS%diag) if (CS%id_Ant_f_adot > 0) call post_scalar_data(CS%id_Ant_f_adot ,val*Itime_step,CS%diag) endif if (CS%id_Ant_bdott > 0 .or. CS%id_Ant_bdot > 0) then !bottom accumulation - bottom melt - call integrate_over_ice_sheet_area(G, ISS, dh_bdott, US%Z_to_m, val, hemisphere=0) + val = integrate_over_ice_sheet_area(G, ISS, dh_bdott, unscale=US%Z_to_m, hemisphere=0) if (CS%id_Ant_bdott > 0) call post_scalar_data(CS%id_Ant_bdott,val ,CS%diag) if (CS%id_Ant_bdot > 0) call post_scalar_data(CS%id_Ant_bdot ,val*Itime_step,CS%diag) endif @@ -2729,7 +2754,7 @@ subroutine process_and_post_scalar_data(CS, vaf0, vaf0_A, vaf0_G, Itime_step, dh do j=js,je ; do i=is,ie if (dh_bdott(i,j) < 0) tmp(i,j) = -dh_bdott(i,j) enddo; enddo - call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val, hemisphere=0) + val = integrate_over_ice_sheet_area(G, ISS, tmp, unscale=US%Z_to_m, hemisphere=0) if (CS%id_Ant_bdott_melt > 0) call post_scalar_data(CS%id_Ant_bdott_melt,val ,CS%diag) if (CS%id_Ant_bdot_melt > 0) call post_scalar_data(CS%id_Ant_bdot_melt ,val*Itime_step,CS%diag) endif @@ -2738,22 +2763,22 @@ subroutine process_and_post_scalar_data(CS, vaf0, vaf0_A, vaf0_G, Itime_step, dh do j=js,je ; do i=is,ie if (dh_bdott(i,j) > 0) tmp(i,j) = dh_bdott(i,j) enddo; enddo - call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val, hemisphere=0) + val = integrate_over_ice_sheet_area(G, ISS, tmp, unscale=US%Z_to_m, hemisphere=0) if (CS%id_Ant_bdott_accum > 0) call post_scalar_data(CS%id_Ant_bdott_accum,val ,CS%diag) if (CS%id_Ant_bdot_accum > 0) call post_scalar_data(CS%id_Ant_bdot_accum ,val*Itime_step,CS%diag) endif if (CS%id_Ant_t_area > 0) then !ice sheet area - tmp(:,:) = 1.0; call integrate_over_ice_sheet_area(G, ISS, tmp, 1.0, val, hemisphere=0) + tmp(:,:) = 1.0; val = integrate_over_ice_sheet_area(G, ISS, tmp, unscale=1.0, hemisphere=0) call post_scalar_data(CS%id_Ant_t_area,val,CS%diag) endif if (CS%id_Ant_g_area > 0 .or. CS%id_Ant_f_area > 0) then - tmp(:,:) = 1.0; call masked_var_grounded(G,CS%dCS,tmp,tmp) + ones(:,:) = 1.0 ; call masked_var_grounded(G, CS%dCS, ones, tmp) if (CS%id_Ant_g_area > 0) then !grounded only ice sheet area - call integrate_over_ice_sheet_area(G, ISS, tmp, 1.0, val, hemisphere=0) + val = integrate_over_ice_sheet_area(G, ISS, tmp, unscale=1.0, hemisphere=0) call post_scalar_data(CS%id_Ant_g_area,val,CS%diag) endif if (CS%id_Ant_f_area > 0) then !floating only ice sheet area (ice shelf area) - call integrate_over_ice_sheet_area(G, ISS, 1.0-tmp, 1.0, val, hemisphere=0) + val = integrate_over_ice_sheet_area(G, ISS, 1.0-tmp, unscale=1.0, hemisphere=0) call post_scalar_data(CS%id_Ant_f_area,val,CS%diag) endif endif @@ -2764,13 +2789,13 @@ subroutine process_and_post_scalar_data(CS, vaf0, vaf0_A, vaf0_G, Itime_step, dh if (CS%id_Gr_vaf > 0) call post_scalar_data(CS%id_Gr_vaf ,vaf ,CS%diag) !current vaf if (CS%id_Gr_dvafdt > 0) call post_scalar_data(CS%id_Gr_dvafdt,(vaf-vaf0_A)*Itime_step,CS%diag) !d(vaf)/dt if (CS%id_Gr_adott > 0 .or. CS%id_Gr_adot > 0) then !surface accumulation - surface melt - call integrate_over_ice_sheet_area(G, ISS, dh_adott, US%Z_to_m, val, hemisphere=1) + val = integrate_over_ice_sheet_area(G, ISS, dh_adott, unscale=US%Z_to_m, hemisphere=1) if (CS%id_Gr_adott > 0) call post_scalar_data(CS%id_Gr_adott,val ,CS%diag) if (CS%id_Gr_adot > 0) call post_scalar_data(CS%id_Gr_adot ,val*Itime_step,CS%diag) endif if (CS%id_Gr_g_adott > 0 .or. CS%id_Gr_g_adot > 0) then !grounded only: surface accumulation - surface melt call masked_var_grounded(G,CS%dCS,dh_adott,tmp) - call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val, hemisphere=1) + val = integrate_over_ice_sheet_area(G, ISS, tmp, unscale=US%Z_to_m, hemisphere=1) if (CS%id_Gr_g_adott > 0) call post_scalar_data(CS%id_Gr_g_adott,val ,CS%diag) if (CS%id_Gr_g_adot > 0) call post_scalar_data(CS%id_Gr_g_adot ,val*Itime_step,CS%diag) endif @@ -2779,12 +2804,12 @@ subroutine process_and_post_scalar_data(CS, vaf0, vaf0_A, vaf0_G, Itime_step, dh do j=js,je ; do i=is,ie tmp(i,j) = dh_adott(i,j) - tmp(i,j) enddo; enddo - call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val, hemisphere=1) + val = integrate_over_ice_sheet_area(G, ISS, tmp, unscale=US%Z_to_m, hemisphere=1) if (CS%id_Gr_f_adott > 0) call post_scalar_data(CS%id_Gr_f_adott,val ,CS%diag) if (CS%id_Gr_f_adot > 0) call post_scalar_data(CS%id_Gr_f_adot ,val*Itime_step,CS%diag) endif if (CS%id_Gr_bdott > 0 .or. CS%id_Gr_bdot > 0) then !bottom accumulation - bottom melt - call integrate_over_ice_sheet_area(G, ISS, dh_bdott, US%Z_to_m, val, hemisphere=1) + val = integrate_over_ice_sheet_area(G, ISS, dh_bdott, unscale=US%Z_to_m, hemisphere=1) if (CS%id_Gr_bdott > 0) call post_scalar_data(CS%id_Gr_bdott,val ,CS%diag) if (CS%id_Gr_bdot > 0) call post_scalar_data(CS%id_Gr_bdot ,val*Itime_step,CS%diag) endif @@ -2793,7 +2818,7 @@ subroutine process_and_post_scalar_data(CS, vaf0, vaf0_A, vaf0_G, Itime_step, dh do j=js,je ; do i=is,ie if (dh_bdott(i,j) < 0) tmp(i,j) = -dh_bdott(i,j) enddo; enddo - call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val, hemisphere=1) + val = integrate_over_ice_sheet_area(G, ISS, tmp, unscale=US%Z_to_m, hemisphere=1) if (CS%id_Gr_bdott_melt > 0) call post_scalar_data(CS%id_Gr_bdott_melt,val ,CS%diag) if (CS%id_Gr_bdot_melt > 0) call post_scalar_data(CS%id_Gr_bdot_melt ,val*Itime_step,CS%diag) endif @@ -2802,22 +2827,22 @@ subroutine process_and_post_scalar_data(CS, vaf0, vaf0_A, vaf0_G, Itime_step, dh do j=js,je ; do i=is,ie if (dh_bdott(i,j) > 0) tmp(i,j) = dh_bdott(i,j) enddo; enddo - call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val, hemisphere=1) + val = integrate_over_ice_sheet_area(G, ISS, tmp, unscale=US%Z_to_m, hemisphere=1) if (CS%id_Gr_bdott_accum > 0) call post_scalar_data(CS%id_Gr_bdott_accum,val ,CS%diag) if (CS%id_Gr_bdot_accum > 0) call post_scalar_data(CS%id_Gr_bdot_accum ,val*Itime_step,CS%diag) endif if (CS%id_Gr_t_area > 0) then !ice sheet area - tmp(:,:) = 1.0; call integrate_over_ice_sheet_area(G, ISS, tmp, 1.0, val, hemisphere=1) + tmp(:,:) = 1.0; val = integrate_over_ice_sheet_area(G, ISS, tmp, unscale=1.0, hemisphere=1) call post_scalar_data(CS%id_Gr_t_area,val,CS%diag) endif if (CS%id_Gr_g_area > 0 .or. CS%id_Gr_f_area > 0) then - tmp(:,:) = 1.0; call masked_var_grounded(G,CS%dCS,tmp,tmp) + ones(:,:) = 1.0 ; call masked_var_grounded(G, CS%dCS, ones, tmp) if (CS%id_Gr_g_area > 0) then !grounded only ice sheet area - call integrate_over_ice_sheet_area(G, ISS, tmp, 1.0, val, hemisphere=1) + val = integrate_over_ice_sheet_area(G, ISS, tmp, unscale=1.0, hemisphere=1) call post_scalar_data(CS%id_Gr_g_area,val,CS%diag) endif if (CS%id_Gr_f_area > 0) then !floating only ice sheet area (ice shelf area) - call integrate_over_ice_sheet_area(G, ISS, 1.0-tmp, 1.0, val, hemisphere=1) + val = integrate_over_ice_sheet_area(G, ISS, 1.0-tmp, unscale=1.0, hemisphere=1) call post_scalar_data(CS%id_Gr_f_area,val,CS%diag) endif endif diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index 9c7dda22de..09eea05127 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -108,10 +108,10 @@ module MOM_ice_shelf_dynamics real, pointer, dimension(:,:) :: basal_traction => NULL() !< The area-integrated taub_beta field !! (m2 Pa s m-1, or kg s-1) related to the nonlinear part - !! of "linearized" basal stress (Pa) [R L3 T-1 ~> kg s-1] + !! of "linearized" basal stress (Pa) [R Z L2 T-1 ~> kg s-1] !! The exact form depends on basal law exponent and/or whether flow is "hybridized" a la Goldberg 2011 real, pointer, dimension(:,:) :: C_basal_friction => NULL()!< Coefficient in sliding law tau_b = C u^(n_basal_fric), - !! units= Pa (s m-1)^(n_basal_fric) + !! units= [Pa (s m-1)^(n_basal_fric)] real, pointer, dimension(:,:) :: OD_rt => NULL() !< A running total for calculating OD_av [Z ~> m]. real, pointer, dimension(:,:) :: ground_frac_rt => NULL() !< A running total for calculating ground_frac. real, pointer, dimension(:,:) :: OD_av => NULL() !< The time average open ocean depth [Z ~> m]. @@ -169,7 +169,7 @@ module MOM_ice_shelf_dynamics !! i.e. dt <= CFL_factor * min(dx / u) [nondim] real :: min_h_shelf !< The minimum ice thickness used during ice dynamics [L ~> m]. - real :: min_basal_traction !< The minimum basal traction for grounded ice (Pa m-1 s) [R L T-1 ~> kg m-2 s-1] + real :: min_basal_traction !< The minimum basal traction for grounded ice (Pa m-1 s) [R Z T-1 ~> kg m-2 s-1] real :: max_surface_slope !< The maximum allowed ice-sheet surface slope (to ignore, set to zero) [nondim] real :: min_ice_visc !< The minimum allowed Glen's law ice viscosity (Pa s), in [R L2 T-1 ~> kg m-1 s-1]. @@ -358,7 +358,7 @@ subroutine register_ice_shelf_dyn_restarts(G, US, param_file, CS, restart_CS) allocate(CS%t_shelf(isd:ied,jsd:jed), source=T_shelf_missing) ! [C ~> degC] 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%basal_traction(isd:ied,jsd:jed), source=0.0) ! [R Z L2 T-1 ~> kg s-1] allocate(CS%C_basal_friction(isd:ied,jsd:jed), source=5.0e10) ! [Pa (s m-1)^n_sliding] allocate(CS%OD_av(isd:ied,jsd:jed), source=0.0) allocate(CS%ground_frac(isd:ied,jsd:jed), source=0.0) @@ -838,7 +838,7 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ CS%id_visc_shelf = register_diag_field('ice_shelf_model','ice_visc',CS%diag%axesT1, Time, & 'vi-viscosity', 'Pa m s', conversion=US%RL2_T2_to_Pa*US%Z_to_m*US%T_to_s) !vertically integrated viscosity CS%id_taub = register_diag_field('ice_shelf_model','taub_beta',CS%diag%axesT1, Time, & - 'taub', 'MPa s m-1', conversion=1e-6*US%RL2_T2_to_Pa/(365.0*86400.0*US%L_T_to_m_s)) + 'taub', units='MPa yr m-1', conversion=1e-6*US%RLZ_T2_to_Pa/(365.0*86400.0*US%L_T_to_m_s)) 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) @@ -1010,10 +1010,10 @@ subroutine volume_above_floatation(CS, G, ISS, vaf, hemisphere) type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf. type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe !! the ice-shelf state - real, intent(out) :: vaf !< area integrated volume above floatation [m3] + real, intent(out) :: vaf !< area integrated volume above floatation [Z L2 ~> m3] integer, optional, intent(in) :: hemisphere !< 0 for Antarctica only, 1 for Greenland only. Otherwise, all ice sheets integer :: IS_ID ! local copy of hemisphere - real, dimension(SZI_(G),SZJ_(G)) :: vaf_cell !< cell-wise volume above floatation [m3] + real, dimension(SZI_(G),SZJ_(G)) :: vaf_cell !< cell-wise volume above floatation [Z L2 ~> m3] integer, dimension(SZI_(G),SZJ_(G)) :: mask ! a mask for active cells depending on hemisphere indicated integer :: is,ie,js,je,i,j real :: rhoi_rhow, rhow_rhoi @@ -1049,16 +1049,15 @@ subroutine volume_above_floatation(CS, G, ISS, vaf, hemisphere) if (mask(i,j)>0) then if (CS%bed_elev(i,j) <= 0) then !grounded above sea level - vaf_cell(i,j)= (ISS%h_shelf(i,j) * G%US%Z_to_m) * (ISS%area_shelf_h(i,j) * G%US%L_to_m**2) + vaf_cell(i,j) = ISS%h_shelf(i,j) * ISS%area_shelf_h(i,j) else !grounded if vaf_cell(i,j) > 0 - vaf_cell(i,j) = (max(ISS%h_shelf(i,j) - rhow_rhoi * CS%bed_elev(i,j), 0.0) * G%US%Z_to_m) * & - (ISS%area_shelf_h(i,j) * G%US%L_to_m**2) + vaf_cell(i,j) = max(ISS%h_shelf(i,j) - rhow_rhoi * CS%bed_elev(i,j), 0.0) * ISS%area_shelf_h(i,j) endif endif enddo; enddo - vaf = reproducing_sum(vaf_cell) + vaf = reproducing_sum(vaf_cell, unscale=G%US%Z_to_m*G%US%L_to_m**2) end subroutine volume_above_floatation !> multiplies a variable with the ice sheet grounding fraction @@ -1193,7 +1192,8 @@ subroutine write_ice_shelf_energy(CS, G, US, mass, area, day, time_step) ! Local variables type(time_type) :: dt ! A time_type version of the timestep. real, dimension(SZDI_(G),SZDJ_(G)) :: tmp1 ! A temporary array used in reproducing sums [various] - real :: KE_tot, mass_tot, KE_scale_factor, mass_scale_factor + real :: KE_tot ! THe total kinetic energy [J] + real :: mass_tot ! The total mass [kg] integer :: is, ie, js, je, isr, ier, jsr, jer, i, j character(len=32) :: mesg_intro, time_units, day_str, n_str, date_str integer :: start_of_day, num_days @@ -1243,24 +1243,23 @@ subroutine write_ice_shelf_energy(CS, G, US, mass, area, day, time_step) isr = is - (G%isd-1) ; ier = ie - (G%isd-1) ; jsr = js - (G%jsd-1) ; jer = je - (G%jsd-1) !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) + tmp1(:,:) = 0.0 do j=js,je ; do i=is,ie - tmp1(i,j) = (KE_scale_factor * 0.03125) * (mass(i,j) * area(i,j)) * & + tmp1(i,j) = 0.03125 * (mass(i,j) * area(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) + KE_tot = US%RZL2_to_kg*US%L_T_to_m_s**2 * & + reproducing_sum(tmp1, isr, ier, jsr, jer, unscale=(US%RZL2_to_kg*US%L_T_to_m_s**2)) !calculate mass - tmp1(:,:)=0.0 - mass_scale_factor = US%L_to_m**2 * US%RZ_to_kg_m2 + tmp1(:,:) = 0.0 do j=js,je ; do i=is,ie - tmp1(i,j) = mass_scale_factor * (mass(i,j) * area(i,j)) + tmp1(i,j) = mass(i,j) * area(i,j) enddo; enddo - mass_tot = reproducing_sum(tmp1, isr, ier, jsr, jer) + mass_tot = US%RZL2_to_kg * reproducing_sum(tmp1, isr, ier, jsr, jer, unscale=US%RZL2_to_kg) if (is_root_pe()) then ! Only the root PE actually writes anything. if (day > CS%Start_time) then @@ -1449,12 +1448,15 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i 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, indicates cells containing ! the grounding line (float_cond=1) or not (float_cond=0) - real, dimension(SZDIB_(G),SZDJB_(G)) :: Normvec ! Used for convergence + real, dimension(SZDIB_(G),SZDJB_(G)) :: Normvec ! Velocities used for convergence [L2 T-2 ~> m2 s-2] character(len=160) :: mesg ! The text of an error message integer :: conv_flag, i, j, k,l, iter, nodefloat integer :: Isdq, Iedq, Jsdq, Jedq, isd, ied, jsd, jed integer :: Iscq, Iecq, Jscq, Jecq, isc, iec, jsc, jec - real :: err_max, err_tempu, err_tempv, err_init, max_vel, tempu, tempv, Norm, PrevNorm + real :: err_max, err_tempu, err_tempv, err_init ! Errors in [R L3 Z T-2 ~> kg m s-2] or [L T-1 ~> m s-1] + real :: max_vel ! The maximum velocity magnitude [L T-1 ~> m s-1] + real :: tempu, tempv ! Temporary variables with velocity magnitudes [L T-1 ~> m s-1] + real :: Norm, PrevNorm ! Velocities used to assess convergence [L T-1 ~> m s-1] real :: rhoi_rhow ! The density of ice divided by a typical water density [nondim] integer :: Is_sum, Js_sum, Ie_sum, Je_sum ! Loop bounds for global sums or arrays starting at 1. integer :: Iscq_sv, Jscq_sv ! Starting loop bound for sum_vec @@ -1553,7 +1555,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i call max_across_PEs(err_init) elseif (CS%nonlin_solve_err_mode == 3) then - Normvec=0.0 + Normvec(:,:) = 0.0 ! Determine the loop limits for sums, bearing in mind that the arrays will be starting at 1. ! Includes the edge of the tile is at the western/southern bdry (if symmetric) @@ -1570,11 +1572,10 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i Ie_sum = Iecq + (1-Isdq) ; Je_sum = Jecq + (1-Jsdq) do J=Jscq_sv,Jecq ; do I=Iscq_sv,Iecq - 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) + if (CS%umask(I,J) == 1) Normvec(I,J) = Normvec(I,J) + u_shlf(I,J)**2 + if (CS%vmask(I,J) == 1) Normvec(I,J) = Normvec(I,J) + v_shlf(I,J)**2 enddo ; enddo - Norm = reproducing_sum( Normvec, Is_sum, Ie_sum, Js_sum, Je_sum ) - Norm = sqrt(Norm) + Norm = sqrt( reproducing_sum( Normvec, Is_sum, Ie_sum, Js_sum, Je_sum, unscale=US%L_T_to_m_s**2 ) ) endif u_last(:,:) = u_shlf(:,:) ; v_last(:,:) = v_shlf(:,:) @@ -1662,14 +1663,13 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i err_init = max_vel elseif (CS%nonlin_solve_err_mode == 3) then - PrevNorm=Norm; Norm=0.0; Normvec=0.0 + PrevNorm = Norm ; Norm = 0.0 ; Normvec=0.0 do J=Jscq_sv,Jecq ; do I=Iscq_sv,Iecq - 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) + if (CS%umask(I,J) == 1) Normvec(I,J) = Normvec(I,J) + u_shlf(I,J)**2 + if (CS%vmask(I,J) == 1) Normvec(I,J) = Normvec(I,J) + v_shlf(I,J)**2 enddo; enddo - Norm = reproducing_sum( Normvec, Is_sum, Ie_sum, Js_sum, Je_sum ) - Norm = sqrt(Norm) - err_max=2.*abs(Norm-PrevNorm); err_init=Norm+PrevNorm + Norm = sqrt( reproducing_sum( Normvec, Is_sum, Ie_sum, Js_sum, Je_sum, unscale=US%L_T_to_m_s**2 ) ) + err_max = 2.*abs(Norm-PrevNorm) ; err_init = Norm+PrevNorm endif write(mesg,*) "ice_shelf_solve_outer: nonlinear fractional residual = ", err_max/err_init @@ -1797,7 +1797,7 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H call pass_vector(Au, Av, G%domain, TO_ALL, BGRID_NE, complete=.true.) 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) + resid_scale = US%s_to_T*(US%RZL2_to_kg*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 @@ -2642,7 +2642,7 @@ subroutine CG_action(CS, uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask, !! relative to sea-level [Z ~> m]. 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]. + !! part of the "linearized" basal stress [R Z L2 T-1 ~> kg s-1]. real, intent(in) :: dens_ratio !< The density of ice divided by the density !! of seawater, nondimensional @@ -2675,10 +2675,11 @@ subroutine CG_action(CS, uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask, real :: uq, vq ! Interpolated velocities [L T-1 ~> m s-1] 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, dimension(2,2,4) :: uret_qp, vret_qp - real, dimension(SZDIB_(G),SZDJB_(G),4) :: uret_b, vret_b + real, dimension(2) :: xquad ! Nondimensional quadrature ratios [nondim] + real, dimension(2,2) :: Ucell, Vcell, Usub, Vsub ! Velocities at the nodal points around the cell [L T-1 ~> m s-1] + real, dimension(2,2) :: Hcell ! Ice shelf thickness at notal (corner) points [Z ~> m] + real, dimension(2,2,4) :: uret_qp, vret_qp ! Temporary arrays in [R Z L3 T-2 ~> kg m s-2] + real, dimension(SZDIB_(G),SZDJB_(G),4) :: uret_b, vret_b ! Temporary arrays in [R Z L3 T-2 ~> kg m s-2] xquad(1) = .5 * (1-sqrt(1./3)) ; xquad(2) = .5 * (1+sqrt(1./3)) @@ -3316,9 +3317,9 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) (v_shlf(I,J-1) * CS%PhiC(4,i,j))) CS%ice_visc(i,j,1) = (G%areaT(i,j) * max(ISS%h_shelf(i,j),CS%min_h_shelf)) * & - max(0.5 * Visc_coef * & - (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),CS%min_ice_visc) + max(0.5 * Visc_coef * & + (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),CS%min_ice_visc) ! Rescale after the fractional power law. elseif (model_qp4) then !calculate viscosity at 4 quadrature points per cell @@ -3347,9 +3348,9 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) (v_shlf(I-1,J) * CS%Phi(6,2*(jq-1)+iq,i,j))) CS%ice_visc(i,j,2*(jq-1)+iq) = (G%areaT(i,j) * max(ISS%h_shelf(i,j),CS%min_h_shelf)) * & - max(0.5 * Visc_coef * & - (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),CS%min_ice_visc) + max(0.5 * Visc_coef * & + (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),CS%min_ice_visc) ! Rescale after the fractional power law. enddo; enddo endif endif @@ -3376,7 +3377,9 @@ subroutine calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) integer :: i, j, iscq, iecq, jscq, jecq, isd, jsd, ied, jed, iegq, jegq integer :: giec, gjec, gisc, gjsc, isc, jsc, iec, jec, is, js - real :: umid, vmid, unorm, eps_min ! Velocities [L T-1 ~> m s-1] + real :: umid, vmid ! Velocities [L T-1 ~> m s-1] + real :: eps_min ! A minimal strain rate used in the Glens flow law expression [T-1 ~> s-1] + real :: unorm ! The magnitude of the velocity in mks units for use with fractional powers [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 [Pa] @@ -3399,7 +3402,7 @@ subroutine calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) else alpha = 1.0 endif - fN_scale = US%R_to_kg_m3 * US%L_T_to_m_s**2 + fN_scale = US%R_to_kg_m3 * US%L_T_to_m_s**2 ! Unscaling factor for use in the fractional power law. endif do j=jsd+1,jed @@ -3417,12 +3420,12 @@ subroutine calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) 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))) * & - (US%Pa_to_RLZ_T2*US%L_T_to_m_s) + (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) ! Restore the scaling after the fractional power law. else !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) + (US%Pa_to_RLZ_T2*US%L_T_to_m_s) ! Rescale after the fractional power law. endif CS%basal_traction(i,j)=max(CS%basal_traction(i,j), CS%min_basal_traction * G%areaT(i,j)) diff --git a/src/ice_shelf/MOM_ice_shelf_initialize.F90 b/src/ice_shelf/MOM_ice_shelf_initialize.F90 index f976187c2b..3dfe1045cf 100644 --- a/src/ice_shelf/MOM_ice_shelf_initialize.F90 +++ b/src/ice_shelf/MOM_ice_shelf_initialize.F90 @@ -321,7 +321,7 @@ subroutine initialize_ice_shelf_boundary_channel(u_face_mask_bdry, v_face_mask_b call get_param(PF, mdl, "INPUT_VEL_ICE_SHELF", input_vel, & "inflow ice velocity at upstream boundary", & - units="m s-1", default=0., scale=US%m_s_to_L_T*US%m_to_Z) !### This conversion factor is wrong? + units="m s-1", default=0., scale=US%m_s_to_L_T) call get_param(PF, mdl, "INPUT_THICK_ICE_SHELF", input_thick, & "flux thickness at upstream boundary", & units="m", default=1000., scale=US%m_to_Z) @@ -556,7 +556,7 @@ end subroutine initialize_ice_shelf_boundary_from_file subroutine initialize_ice_C_basal_friction(C_basal_friction, G, US, PF) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure real, dimension(SZDI_(G),SZDJ_(G)), & - intent(inout) :: C_basal_friction !< Ice-stream basal friction + intent(inout) :: C_basal_friction !< Ice-stream basal friction [Pa (s m-1)^(n_basal_fric)] 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