diff --git a/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 b/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 index ad99f8482..bc5823be8 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 @@ -122,7 +122,8 @@ subroutine eap (dt) use ice_dyn_shared, only: fcor_blk, ndte, dtei, & denom1, uvel_init, vvel_init, arlx1i, & dyn_prep1, dyn_prep2, stepu, dyn_finish, & - basal_stress_coeff, basalstress, ice_HaloUpdate_vel + basal_stress_coeff, basalstress, & + stack_velocity_field, unstack_velocity_field use ice_flux, only: rdg_conv, strairxT, strairyT, & strairx, strairy, uocn, vocn, ss_tltx, ss_tlty, iceumask, fm, & strtltx, strtlty, strocnx, strocny, strintx, strinty, taubx, tauby, & @@ -172,6 +173,8 @@ subroutine eap (dt) umass , & ! total mass of ice and snow (u grid) umassdti ! mass of U-cell/dte (kg/m^2 s) + real (kind=dbl_kind), allocatable :: fld2(:,:,:,:) + real (kind=dbl_kind), dimension(nx_block,ny_block,8):: & strtmp ! stress combinations for momentum equation @@ -195,6 +198,8 @@ subroutine eap (dt) ! Initialize !----------------------------------------------------------------- + allocate(fld2(nx_block,ny_block,2,max_blocks)) + ! This call is needed only if dt changes during runtime. ! call set_evp_parameters (dt) @@ -373,7 +378,17 @@ subroutine eap (dt) endif ! velocities may have changed in dyn_prep2 - call ice_HaloUpdate_vel(uvel, vvel, halo_info_mask) + call stack_velocity_field(uvel, vvel, fld2) + call ice_timer_start(timer_bound) + if (maskhalo_dyn) then + call ice_HaloUpdate (fld2, halo_info_mask, & + field_loc_NEcorner, field_type_vector) + else + call ice_HaloUpdate (fld2, halo_info, & + field_loc_NEcorner, field_type_vector) + endif + call ice_timer_stop(timer_bound) + call unstack_velocity_field(fld2, uvel, vvel) !----------------------------------------------------------------- ! basal stress coefficients (landfast ice) @@ -480,10 +495,21 @@ subroutine eap (dt) enddo !$TCXOMP END PARALLEL DO - call ice_HaloUpdate_vel(uvel, vvel, halo_info_mask) + call stack_velocity_field(uvel, vvel, fld2) + call ice_timer_start(timer_bound) + if (maskhalo_dyn) then + call ice_HaloUpdate (fld2, halo_info_mask, & + field_loc_NEcorner, field_type_vector) + else + call ice_HaloUpdate (fld2, halo_info, & + field_loc_NEcorner, field_type_vector) + endif + call ice_timer_stop(timer_bound) + call unstack_velocity_field(fld2, uvel, vvel) enddo ! subcycling + deallocate(fld2) if (maskhalo_dyn) call ice_HaloDestroy(halo_info_mask) !----------------------------------------------------------------- diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index af06e5d70..c4312c325 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -94,7 +94,7 @@ subroutine evp (dt) ice_timer_start, ice_timer_stop, timer_evp_1d, timer_evp_2d use ice_dyn_evp_1d, only: ice_dyn_evp_1d_copyin, ice_dyn_evp_1d_kernel, & ice_dyn_evp_1d_copyout - use ice_dyn_shared, only: kevp_kernel, ice_HaloUpdate_vel + use ice_dyn_shared, only: kevp_kernel, stack_velocity_field, unstack_velocity_field real (kind=dbl_kind), intent(in) :: & dt ! time step @@ -127,6 +127,8 @@ subroutine evp (dt) umass , & ! total mass of ice and snow (u grid) umassdti ! mass of U-cell/dte (kg/m^2 s) + real (kind=dbl_kind), allocatable :: fld2(:,:,:,:) + real (kind=dbl_kind), dimension(nx_block,ny_block,8):: & strtmp ! stress combinations for momentum equation @@ -152,6 +154,8 @@ subroutine evp (dt) ! Initialize !----------------------------------------------------------------- + allocate(fld2(nx_block,ny_block,2,max_blocks)) + ! This call is needed only if dt changes during runtime. ! call set_evp_parameters (dt) @@ -316,7 +320,17 @@ subroutine evp (dt) endif ! velocities may have changed in dyn_prep2 - call ice_HaloUpdate_vel(uvel, vvel, halo_info_mask) + call stack_velocity_field(uvel, vvel, fld2) + call ice_timer_start(timer_bound) + if (maskhalo_dyn) then + call ice_HaloUpdate (fld2, halo_info_mask, & + field_loc_NEcorner, field_type_vector) + else + call ice_HaloUpdate (fld2, halo_info, & + field_loc_NEcorner, field_type_vector) + endif + call ice_timer_stop(timer_bound) + call unstack_velocity_field(fld2, uvel, vvel) !----------------------------------------------------------------- ! basal stress coefficients (landfast ice) @@ -429,12 +443,23 @@ subroutine evp (dt) enddo !$TCXOMP END PARALLEL DO - call ice_HaloUpdate_vel(uvel, vvel, halo_info_mask) + call stack_velocity_field(uvel, vvel, fld2) + call ice_timer_start(timer_bound) + if (maskhalo_dyn) then + call ice_HaloUpdate (fld2, halo_info_mask, & + field_loc_NEcorner, field_type_vector) + else + call ice_HaloUpdate (fld2, halo_info, & + field_loc_NEcorner, field_type_vector) + endif + call ice_timer_stop(timer_bound) + call unstack_velocity_field(fld2, uvel, vvel) enddo ! subcycling endif ! kevp_kernel call ice_timer_stop(timer_evp_2d) + deallocate(fld2) if (maskhalo_dyn) call ice_HaloDestroy(halo_info_mask) ! Force symmetry across the tripole seam diff --git a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 index 202df51bf..56bd54b00 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -24,7 +24,8 @@ module ice_dyn_shared private public :: init_evp, set_evp_parameters, stepu, principal_stress, & dyn_prep1, dyn_prep2, dyn_finish, basal_stress_coeff, & - alloc_dyn_shared, deformations, strain_rates, ice_HaloUpdate_vel + alloc_dyn_shared, deformations, strain_rates, & + stack_velocity_field, unstack_velocity_field ! namelist parameters @@ -94,9 +95,6 @@ module ice_dyn_shared ! see keel data from Amundrud et al. 2004 (JGR) u0 = 5e-5_dbl_kind ! residual velocity for basal stress (m/s) - real (kind=dbl_kind), allocatable :: & - fld2(:,:,:,:) ! work array for boundary updates - !======================================================================= contains @@ -112,7 +110,6 @@ subroutine alloc_dyn_shared allocate( & uvel_init (nx_block,ny_block,max_blocks), & ! x-component of velocity (m/s), beginning of timestep vvel_init (nx_block,ny_block,max_blocks), & ! y-component of velocity (m/s), beginning of timestep - fld2 (nx_block,ny_block,2,max_blocks), & ! work array for boundary updates stat=ierr) if (ierr/=0) call abort_ice('(alloc_dyn_shared): Out of memory') @@ -1183,29 +1180,25 @@ end subroutine strain_rates !======================================================================= -! Perform a halo update for the velocity field -! author: Philippe Blain, ECCC - - subroutine ice_HaloUpdate_vel(uvel, vvel, halo_info_mask) +! Load velocity components into array for boundary updates - use ice_boundary, only: ice_HaloUpdate, ice_halo - use ice_constants, only: field_loc_NEcorner, field_type_vector - use ice_domain, only: halo_info, maskhalo_dyn, nblocks - use ice_timers, only: timer_bound, ice_timer_start, ice_timer_stop + subroutine stack_velocity_field(uvel, vvel, fld2) - type (ice_halo), intent(in) :: & - halo_info_mask ! ghost cell update info for masked halo + use ice_domain, only: nblocks - real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), intent(inout) :: & + real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), intent(in) :: & uvel , & ! u components of velocity vector vvel ! v components of velocity vector + real (kind=dbl_kind), dimension (nx_block,ny_block,2,max_blocks), intent(out) :: & + fld2 ! work array for boundary updates + ! local variables integer (kind=int_kind) :: & iblk ! block index - character(len=*), parameter :: subname = '(ice_HaloUpdate_vel)' + character(len=*), parameter :: subname = '(stack_velocity_field)' ! load velocity into array for boundary updates !$OMP PARALLEL DO PRIVATE(iblk) @@ -1215,17 +1208,31 @@ subroutine ice_HaloUpdate_vel(uvel, vvel, halo_info_mask) enddo !$OMP END PARALLEL DO - call ice_timer_start(timer_bound) - if (maskhalo_dyn) then - call ice_HaloUpdate (fld2, halo_info_mask, & - field_loc_NEcorner, field_type_vector) - else - call ice_HaloUpdate (fld2, halo_info, & - field_loc_NEcorner, field_type_vector) - endif - call ice_timer_stop(timer_bound) + end subroutine stack_velocity_field + +!======================================================================= + +! Unload velocity components from array after boundary updates + + subroutine unstack_velocity_field(fld2, uvel, vvel) + + use ice_domain, only: nblocks + + real (kind=dbl_kind), dimension (nx_block,ny_block,2,max_blocks), intent(in) :: & + fld2 ! work array for boundary updates - ! Unload + real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), intent(out) :: & + uvel , & ! u components of velocity vector + vvel ! v components of velocity vector + + ! local variables + + integer (kind=int_kind) :: & + iblk ! block index + + character(len=*), parameter :: subname = '(unstack_velocity_field)' + + ! Unload velocity from array after boundary updates !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks uvel(:,:,iblk) = fld2(:,:,1,iblk) @@ -1233,10 +1240,10 @@ subroutine ice_HaloUpdate_vel(uvel, vvel, halo_info_mask) enddo !$OMP END PARALLEL DO - end subroutine ice_HaloUpdate_vel + end subroutine unstack_velocity_field !======================================================================= - + end module ice_dyn_shared !======================================================================= diff --git a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 index a448eb6e3..6f43b2fe1 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 @@ -47,7 +47,8 @@ module ice_dyn_vp use ice_domain_size, only: max_blocks use ice_dyn_shared, only: dyn_prep1, dyn_prep2, dyn_finish, & ecci, cosw, sinw, fcor_blk, uvel_init, & - vvel_init, basal_stress_coeff, basalstress, Ktens, ice_HaloUpdate_vel + vvel_init, basal_stress_coeff, basalstress, Ktens, & + stack_velocity_field, unstack_velocity_field use ice_fileunits, only: nu_diag use ice_flux, only: fm use ice_global_reductions, only: global_sum, global_allreduce_sum @@ -102,6 +103,9 @@ module ice_dyn_vp indxui(:,:) , & ! compressed index in i-direction indxuj(:,:) ! compressed index in j-direction + real (kind=dbl_kind), allocatable :: & + fld2(:,:,:,:) ! work array for boundary updates + !======================================================================= contains @@ -145,6 +149,7 @@ subroutine init_vp (dt) indxtj(nx_block*ny_block, max_blocks), & indxui(nx_block*ny_block, max_blocks), & indxuj(nx_block*ny_block, max_blocks)) + allocate(fld2(nx_block,ny_block,2,max_blocks)) ! Redefine tinyarea using min_strain_rate @@ -433,7 +438,17 @@ subroutine implicit_solver (dt) endif ! velocities may have changed in dyn_prep2 - call ice_HaloUpdate_vel(uvel, vvel, halo_info_mask) + call stack_velocity_field(uvel, vvel, fld2) + call ice_timer_start(timer_bound) + if (maskhalo_dyn) then + call ice_HaloUpdate (fld2, halo_info_mask, & + field_loc_NEcorner, field_type_vector) + else + call ice_HaloUpdate (fld2, halo_info, & + field_loc_NEcorner, field_type_vector) + endif + call ice_timer_stop(timer_bound) + call unstack_velocity_field(fld2, uvel, vvel) !----------------------------------------------------------------- ! basal stress coefficients (landfast ice) @@ -664,13 +679,13 @@ subroutine anderson_solver (icellt , icellu, & use ice_blocks, only: nx_block, ny_block use ice_boundary, only: ice_HaloUpdate use ice_constants, only: c1 - use ice_domain, only: maskhalo_dyn + use ice_domain, only: maskhalo_dyn, halo_info use ice_domain_size, only: max_blocks use ice_flux, only: uocn, vocn, fm, Tbu use ice_grid, only: dxt, dyt, dxhy, dyhx, cxp, cyp, cxm, cym, & uarear, tinyarea use ice_state, only: uvel, vvel, strength - use ice_timers, only: ice_timer_start, ice_timer_stop + use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bound integer (kind=int_kind), intent(in) :: & ntot ! size of problem for Anderson @@ -1069,7 +1084,17 @@ subroutine anderson_solver (icellt , icellu, & uvel (:,:,:), vvel (:,:,:)) ! Do halo update so that halo cells contain up to date info for advection - call ice_HaloUpdate_vel(uvel, vvel, halo_info_mask) + call stack_velocity_field(uvel, vvel, fld2) + call ice_timer_start(timer_bound) + if (maskhalo_dyn) then + call ice_HaloUpdate (fld2, halo_info_mask, & + field_loc_NEcorner, field_type_vector) + else + call ice_HaloUpdate (fld2, halo_info, & + field_loc_NEcorner, field_type_vector) + endif + call ice_timer_stop(timer_bound) + call unstack_velocity_field(fld2, uvel, vvel) ! Compute "progress" residual norm !$OMP PARALLEL DO PRIVATE(iblk) @@ -2644,6 +2669,10 @@ subroutine fgmres (zetaD , & solx , soly , & nbiter) + use ice_boundary, only: ice_HaloUpdate + use ice_domain, only: maskhalo_dyn, halo_info + use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bound + real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks,4), intent(in) :: & zetaD ! zetaD = 2*zeta (viscous coefficient) @@ -2839,8 +2868,18 @@ subroutine fgmres (zetaD , & orig_basis_y(:,:,:,initer) = workspace_y ! Update workspace with boundary values - call ice_HaloUpdate_vel(workspace_x, workspace_y, & - halo_info_mask) + call stack_velocity_field(workspace_x, workspace_y, fld2) + call ice_timer_start(timer_bound) + if (maskhalo_dyn) then + call ice_HaloUpdate (fld2, halo_info_mask, & + field_loc_NEcorner, field_type_vector) + else + call ice_HaloUpdate (fld2, halo_info, & + field_loc_NEcorner, field_type_vector) + endif + call ice_timer_stop(timer_bound) + call unstack_velocity_field(fld2, workspace_x, workspace_y) + !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks call matvec (nx_block , ny_block , &