diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index bac5b0fce9..582518d4f1 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -874,13 +874,15 @@ 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 unscaled units [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] @@ -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 = unscale*G%US%L_to_m**2 * 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 @@ -2385,7 +2385,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 +2396,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. @@ -2618,6 +2613,7 @@ subroutine process_and_post_scalar_data(CS, vaf0, vaf0_A, vaf0_G, Itime_step, dh 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] real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: tmp ! Temporary field used when calculating diagnostics [various] + 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 [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 @@ -2636,13 +2632,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 +2647,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 +2661,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 +2670,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 +2696,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 +2711,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 +2725,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 +2734,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 +2760,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 +2775,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 +2789,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 +2798,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..140b9746d8 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -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 = G%US%Z_to_m*G%US%L_to_m**2 * 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 @@ -3316,9 +3316,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 +3347,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 +3376,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 +3401,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 +3419,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))