From 2435fa7824b88903acb7d2d3a616d5975dee8164 Mon Sep 17 00:00:00 2001 From: Tony Craig Date: Mon, 17 Oct 2022 06:42:33 -0700 Subject: [PATCH] Change icetmask to logical consistent with iceumask, icenmask, iceemask (#773) * Change icetmask to logical consistent with iceumask, icenmask, iceemask - Add icetmask as logical array to ice_grid.F90, was integer array - Update use of icetmask in code for consistency with new type - Add ice_HaloUpdate2DL1 to support halo updates for logical fields in both mpi and serial ice_boundary.F90 - Modify some capital T,U,N,E in ice_dyn_shared.F90 to t,u,n,e for better consistency in code * Update cicecore/cicedynB/infrastructure/comm/mpi/ice_boundary.F90 * Update comment in code * Revert changes to T,U,N,E in ice_dyn_shared.F90, working toward additional changes * Move ice[T,U,N,E}mask from ice_grid to ice_dyn_shared * rename icell[t,u,n,e] to icell[T,U,N,E], rename indx[t,u,n,e] to indx[T,U,N,E] * remove ice[t,u,n,e]grid from ice_grid --- cicecore/cicedynB/dynamics/ice_dyn_eap.F90 | 97 ++-- cicecore/cicedynB/dynamics/ice_dyn_evp.F90 | 303 ++++++----- cicecore/cicedynB/dynamics/ice_dyn_evp_1d.F90 | 54 +- cicecore/cicedynB/dynamics/ice_dyn_shared.F90 | 48 +- cicecore/cicedynB/dynamics/ice_dyn_vp.F90 | 509 +++++++++--------- .../infrastructure/comm/mpi/ice_boundary.F90 | 64 +++ .../comm/serial/ice_boundary.F90 | 64 +++ cicecore/cicedynB/infrastructure/ice_grid.F90 | 8 - .../infrastructure/ice_restart_driver.F90 | 32 +- doc/source/cice_index.rst | 8 +- doc/source/user_guide/ug_implementation.rst | 6 +- 11 files changed, 657 insertions(+), 536 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 b/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 index ea254cbc0..28a047c4e 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 @@ -134,7 +134,7 @@ subroutine eap (dt) dyn_prep1, dyn_prep2, stepu, dyn_finish, & seabed_stress_factor_LKD, seabed_stress_factor_prob, & seabed_stress_method, seabed_stress, & - stack_fields, unstack_fields + stack_fields, unstack_fields, iceTmask, iceUmask use ice_flux, only: rdg_conv, strairxT, strairyT, & strairxU, strairyU, uocn, vocn, ss_tltx, ss_tlty, fmU, & strtltxU, strtltyU, strocnxU, strocnyU, strintxU, strintyU, taubxU, taubyU, & @@ -144,7 +144,7 @@ subroutine eap (dt) stressm_1, stressm_2, stressm_3, stressm_4, & stress12_1, stress12_2, stress12_3, stress12_4 use ice_grid, only: tmask, umask, dxT, dyT, dxhy, dyhx, cxp, cyp, cxm, cym, & - tarear, uarear, grid_average_X2Y, iceumask, & + tarear, uarear, grid_average_X2Y, & grid_atm_dynu, grid_atm_dynv, grid_ocn_dynu, grid_ocn_dynv use ice_state, only: aice, aiU, vice, vsno, uvel, vvel, divu, shear, & aice_init, aice0, aicen, vicen, strength @@ -163,14 +163,14 @@ subroutine eap (dt) i, j, ij integer (kind=int_kind), dimension(max_blocks) :: & - icellt , & ! no. of cells where icetmask = 1 - icellu ! no. of cells where iceumask = 1 + icellT , & ! no. of cells where iceTmask = .true. + icellU ! no. of cells where iceUmask = .true. integer (kind=int_kind), dimension (nx_block*ny_block, max_blocks) :: & - indxti , & ! compressed index in i-direction - indxtj , & ! compressed index in j-direction - indxui , & ! compressed index in i-direction - indxuj ! compressed index in j-direction + indxTi , & ! compressed index in i-direction + indxTj , & ! compressed index in j-direction + indxUi , & ! compressed index in i-direction + indxUj ! compressed index in j-direction real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: & uocnU , & ! i ocean current (m/s) @@ -198,7 +198,6 @@ subroutine eap (dt) calc_strair integer (kind=int_kind), dimension (nx_block,ny_block,max_blocks) :: & - icetmask , & ! ice extent mask (T-cell) halomask ! ice mask for halo update type (ice_halo) :: & @@ -256,13 +255,13 @@ subroutine eap (dt) ilo, ihi, jlo, jhi, & aice (:,:,iblk), vice (:,:,iblk), & vsno (:,:,iblk), tmask (:,:,iblk), & - tmass (:,:,iblk), icetmask(:,:,iblk)) + tmass (:,:,iblk), iceTmask(:,:,iblk)) enddo ! iblk !$OMP END PARALLEL DO call ice_timer_start(timer_bound) - call ice_HaloUpdate (icetmask, halo_info, & + call ice_HaloUpdate (iceTmask, halo_info, & field_loc_center, field_type_scalar) call ice_timer_stop(timer_bound) @@ -324,16 +323,16 @@ subroutine eap (dt) call dyn_prep2 (nx_block, ny_block, & ilo, ihi, jlo, jhi, & - icellt (iblk), icellu (iblk), & - indxti (:,iblk), indxtj (:,iblk), & - indxui (:,iblk), indxuj (:,iblk), & + icellT (iblk), icellU (iblk), & + indxTi (:,iblk), indxTj (:,iblk), & + indxUi (:,iblk), indxUj (:,iblk), & aiU (:,:,iblk), umass (:,:,iblk), & umassdti (:,:,iblk), fcor_blk (:,:,iblk), & umask (:,:,iblk), & uocnU (:,:,iblk), vocnU (:,:,iblk), & strairxU (:,:,iblk), strairyU (:,:,iblk), & ss_tltxU (:,:,iblk), ss_tltyU (:,:,iblk), & - icetmask (:,:,iblk), iceumask (:,:,iblk), & + iceTmask (:,:,iblk), iceUmask (:,:,iblk), & fmU (:,:,iblk), dt, & strtltxU (:,:,iblk), strtltyU (:,:,iblk), & strocnxU (:,:,iblk), strocnyU (:,:,iblk), & @@ -357,7 +356,7 @@ subroutine eap (dt) do j = 1, ny_block do i = 1, nx_block - if (icetmask(i,j,iblk)==0) then + if (.not.iceTmask(i,j,iblk)) then if (tmask(i,j,iblk)) then ! structure tensor a11_1(i,j,iblk) = p5 @@ -374,7 +373,7 @@ subroutine eap (dt) a12_2(i,j,iblk) = c0 a12_3(i,j,iblk) = c0 a12_4(i,j,iblk) = c0 - endif ! icetmask + endif ! iceTmask enddo ! i enddo ! j @@ -384,9 +383,9 @@ subroutine eap (dt) !----------------------------------------------------------------- strength(:,:,iblk) = c0 ! initialize - do ij = 1, icellt(iblk) - i = indxti(ij, iblk) - j = indxtj(ij, iblk) + do ij = 1, icellT(iblk) + i = indxTi(ij, iblk) + j = indxTj(ij, iblk) call icepack_ice_strength(ncat=ncat, & aice = aice (i,j, iblk), & vice = vice (i,j, iblk), & @@ -415,7 +414,7 @@ subroutine eap (dt) if (maskhalo_dyn) then call ice_timer_start(timer_bound) halomask = 0 - where (iceumask) halomask = 1 + where (iceUmask) halomask = 1 call ice_HaloUpdate (halomask, halo_info, & field_loc_center, field_type_scalar) call ice_timer_stop(timer_bound) @@ -431,8 +430,8 @@ subroutine eap (dt) !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime) do iblk = 1, nblocks call seabed_stress_factor_LKD (nx_block , ny_block , & - icellu (iblk), & - indxui (:,iblk), indxuj(:,iblk), & + icellU (iblk), & + indxUi (:,iblk), indxUj(:,iblk), & vice (:,:,iblk), aice(:,:,iblk), & hwater(:,:,iblk), TbU (:,:,iblk)) enddo @@ -442,8 +441,8 @@ subroutine eap (dt) !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime) do iblk = 1, nblocks call seabed_stress_factor_prob (nx_block , ny_block , & - icellt(iblk), indxti(:,iblk), indxtj(:,iblk), & - icellu(iblk), indxui(:,iblk), indxuj(:,iblk), & + icellT(iblk), indxTi(:,iblk), indxTj(:,iblk), & + icellU(iblk), indxUi(:,iblk), indxUj(:,iblk), & aicen(:,:,:,iblk), vicen(:,:,:,iblk), & hwater (:,:,iblk), TbU (:,:,iblk)) enddo @@ -463,8 +462,8 @@ subroutine eap (dt) ! call ice_timer_start(timer_tmp1,iblk) call stress_eap (nx_block, ny_block, & ksub, ndte, & - icellt (iblk), & - indxti (:,iblk), indxtj (:,iblk), & + icellT (iblk), & + indxTi (:,iblk), indxTj (:,iblk), & arlx1i, denom1, & uvel (:,:,iblk), vvel (:,:,iblk), & dxT (:,:,iblk), dyT (:,:,iblk), & @@ -501,8 +500,8 @@ subroutine eap (dt) ! call ice_timer_start(timer_tmp2,iblk) call stepu (nx_block, ny_block, & - icellu (iblk), Cdn_ocnU (:,:,iblk), & - indxui (:,iblk), indxuj (:,iblk), & + icellU (iblk), Cdn_ocnU (:,:,iblk), & + indxUi (:,iblk), indxUj (:,iblk), & aiU (:,:,iblk), strtmp (:,:,:), & uocnU (:,:,iblk), vocnU (:,:,iblk), & waterxU (:,:,iblk), wateryU (:,:,iblk), & @@ -523,8 +522,8 @@ subroutine eap (dt) ! call ice_timer_start(timer_tmp3,iblk) if (mod(ksub,10) == 1) then ! only called every 10th timestep call stepa (nx_block , ny_block , & - dtei , icellt (iblk), & - indxti (:,iblk), indxtj (:,iblk), & + dtei , icellT (iblk), & + indxTi (:,iblk), indxTj (:,iblk), & a11 (:,:,iblk), a12 (:,:,iblk), & a11_1 (:,:,iblk), a11_2 (:,:,iblk), & a11_3 (:,:,iblk), a11_4 (:,:,iblk), & @@ -567,8 +566,8 @@ subroutine eap (dt) call dyn_finish & (nx_block, ny_block, & - icellu (iblk), Cdn_ocnU(:,:,iblk), & - indxui (:,iblk), indxuj (:,iblk), & + icellU (iblk), Cdn_ocnU(:,:,iblk), & + indxUi (:,iblk), indxUj (:,iblk), & uvel (:,:,iblk), vvel (:,:,iblk), & uocnU (:,:,iblk), vocnU (:,:,iblk), & aiU (:,:,iblk), fmU (:,:,iblk), & @@ -1154,8 +1153,8 @@ end FUNCTION s22ks subroutine stress_eap (nx_block, ny_block, & ksub, ndte, & - icellt, & - indxti, indxtj, & + icellT, & + indxTi, indxTj, & arlx1i, denom1, & uvel, vvel, & dxT, dyT, & @@ -1187,11 +1186,11 @@ subroutine stress_eap (nx_block, ny_block, & nx_block, ny_block, & ! block dimensions ksub , & ! subcycling step ndte , & ! number of subcycles - icellt ! no. of cells where icetmask = 1 + icellT ! no. of cells where iceTmask = .true. integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: & - indxti , & ! compressed index in i-direction - indxtj ! compressed index in j-direction + indxTi , & ! compressed index in i-direction + indxTj ! compressed index in j-direction real (kind=dbl_kind), intent(in) :: & arlx1i , & ! dte/2T (original) or 1/alpha1 (revised) @@ -1274,9 +1273,9 @@ subroutine stress_eap (nx_block, ny_block, & strtmp(:,:,:) = c0 - do ij = 1, icellt - i = indxti(ij) - j = indxtj(ij) + do ij = 1, icellT + i = indxTi(ij) + j = indxTj(ij) !----------------------------------------------------------------- ! strain rates @@ -1878,8 +1877,8 @@ end subroutine update_stress_rdg ! Solves evolution equation for structure tensor (A19, A20) subroutine stepa (nx_block, ny_block, & - dtei, icellt, & - indxti, indxtj, & + dtei, icellT, & + indxTi, indxTj, & a11, a12, & a11_1, a11_2, a11_3, a11_4, & a12_1, a12_2, a12_3, a12_4, & @@ -1892,14 +1891,14 @@ subroutine stepa (nx_block, ny_block, & integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions - icellt ! no. of cells where icetmask = 1 + icellT ! no. of cells where iceTmask = .true. real (kind=dbl_kind), intent(in) :: & dtei ! 1/dte, where dte is subcycling timestep (1/s) integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: & - indxti , & ! compressed index in i-direction - indxtj ! compressed index in j-direction + indxTi , & ! compressed index in i-direction + indxTj ! compressed index in j-direction real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & ! ice stress tensor (kg/s^2) in each corner of T cell @@ -1929,9 +1928,9 @@ subroutine stepa (nx_block, ny_block, & dteikth = c1 / (dtei + kth) p5kth = p5 * kth - do ij = 1, icellt - i = indxti(ij) - j = indxtj(ij) + do ij = 1, icellT + i = indxTi(ij) + j = indxTj(ij) ! ne call calc_ffrac(stressp_1(i,j), stressm_1(i,j), & diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index d1264bae7..8eab5e260 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -99,7 +99,6 @@ subroutine evp (dt) stresspT, stressmT, stress12T, & stresspU, stressmU, stress12U use ice_grid, only: tmask, umask, umaskCD, nmask, emask, uvm, epm, npm, & - iceumask, iceemask, icenmask, & dxE, dxN, dxT, dxU, dyE, dyN, dyT, dyU, & ratiodxN, ratiodxNr, ratiodyE, ratiodyEr, & dxhy, dyhx, cxp, cyp, cxm, cym, & @@ -116,6 +115,7 @@ subroutine evp (dt) use ice_dyn_shared, only: evp_algorithm, stack_fields, unstack_fields, & DminTarea, visc_method, deformations, deformationsC_T, deformationsCD_T, & strain_rates_U, & + iceTmask, iceUmask, iceEmask, iceNmask, & dyn_haloUpdate real (kind=dbl_kind), intent(in) :: & @@ -130,20 +130,20 @@ subroutine evp (dt) i, j, ij ! local indices integer (kind=int_kind), dimension(max_blocks) :: & - icellt , & ! no. of cells where icetmask = 1 - icelln , & ! no. of cells where icenmask = .true. - icelle , & ! no. of cells where iceemask = .true. - icellu ! no. of cells where iceumask = .true. + icellT , & ! no. of cells where iceTmask = .true. + icellN , & ! no. of cells where iceNmask = .true. + icellE , & ! no. of cells where iceEmask = .true. + icellU ! no. of cells where iceUmask = .true. integer (kind=int_kind), dimension (nx_block*ny_block, max_blocks) :: & - indxti , & ! compressed index in i-direction - indxtj , & ! compressed index in j-direction - indxei , & ! compressed index in i-direction - indxej , & ! compressed index in j-direction - indxni , & ! compressed index in i-direction - indxnj , & ! compressed index in j-direction - indxui , & ! compressed index in i-direction - indxuj ! compressed index in j-direction + indxTi , & ! compressed index in i-direction + indxTj , & ! compressed index in j-direction + indxEi , & ! compressed index in i-direction + indxEj , & ! compressed index in j-direction + indxNi , & ! compressed index in i-direction + indxNj , & ! compressed index in j-direction + indxUi , & ! compressed index in i-direction + indxUj ! compressed index in j-direction real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: & uocnU , & ! i ocean current (m/s) @@ -210,7 +210,6 @@ subroutine evp (dt) calc_strair ! calculate air/ice stress integer (kind=int_kind), dimension (nx_block,ny_block,max_blocks) :: & - icetmask, & ! ice extent mask (T-cell) halomask ! generic halo mask type (ice_halo) :: & @@ -301,13 +300,13 @@ subroutine evp (dt) ilo, ihi, jlo, jhi, & aice (:,:,iblk), vice (:,:,iblk), & vsno (:,:,iblk), tmask (:,:,iblk), & - tmass (:,:,iblk), icetmask(:,:,iblk)) + tmass (:,:,iblk), iceTmask(:,:,iblk)) enddo ! iblk !$OMP END PARALLEL DO call ice_timer_start(timer_bound) - call ice_HaloUpdate (icetmask, halo_info, & + call ice_HaloUpdate (iceTmask, halo_info, & field_loc_center, field_type_scalar) call ice_timer_stop(timer_bound) @@ -400,16 +399,16 @@ subroutine evp (dt) if (trim(grid_ice) == 'B') then call dyn_prep2 (nx_block, ny_block, & ilo, ihi, jlo, jhi, & - icellt (iblk), icellu (iblk), & - indxti (:,iblk), indxtj (:,iblk), & - indxui (:,iblk), indxuj (:,iblk), & + icellT (iblk), icellU (iblk), & + indxTi (:,iblk), indxTj (:,iblk), & + indxUi (:,iblk), indxUj (:,iblk), & aiU (:,:,iblk), umass (:,:,iblk), & umassdti (:,:,iblk), fcor_blk (:,:,iblk), & umask (:,:,iblk), & uocnU (:,:,iblk), vocnU (:,:,iblk), & strairxU (:,:,iblk), strairyU (:,:,iblk), & ss_tltxU (:,:,iblk), ss_tltyU (:,:,iblk), & - icetmask (:,:,iblk), iceumask (:,:,iblk), & + iceTmask (:,:,iblk), iceUmask (:,:,iblk), & fmU (:,:,iblk), dt, & strtltxU (:,:,iblk), strtltyU (:,:,iblk), & strocnxU (:,:,iblk), strocnyU (:,:,iblk), & @@ -430,16 +429,16 @@ subroutine evp (dt) elseif (trim(grid_ice) == 'CD' .or. grid_ice == 'C') then call dyn_prep2 (nx_block, ny_block, & ilo, ihi, jlo, jhi, & - icellt (iblk), icellu (iblk), & - indxti (:,iblk), indxtj (:,iblk), & - indxui (:,iblk), indxuj (:,iblk), & + icellT (iblk), icellU (iblk), & + indxTi (:,iblk), indxTj (:,iblk), & + indxUi (:,iblk), indxUj (:,iblk), & aiU (:,:,iblk), umass (:,:,iblk), & umassdti (:,:,iblk), fcor_blk (:,:,iblk), & umaskCD (:,:,iblk), & uocnU (:,:,iblk), vocnU (:,:,iblk), & strairxU (:,:,iblk), strairyU (:,:,iblk), & ss_tltxU (:,:,iblk), ss_tltyU (:,:,iblk), & - icetmask (:,:,iblk), iceumask (:,:,iblk), & + iceTmask (:,:,iblk), iceUmask (:,:,iblk), & fmU (:,:,iblk), dt, & strtltxU (:,:,iblk), strtltyU (:,:,iblk), & strocnxU (:,:,iblk), strocnyU (:,:,iblk), & @@ -463,9 +462,9 @@ subroutine evp (dt) !----------------------------------------------------------------- strength(:,:,iblk) = c0 ! initialize - do ij = 1, icellt(iblk) - i = indxti(ij, iblk) - j = indxtj(ij, iblk) + do ij = 1, icellT(iblk) + i = indxTi(ij, iblk) + j = indxTj(ij, iblk) call icepack_ice_strength(ncat = ncat, & aice = aice (i,j, iblk), & vice = vice (i,j, iblk), & @@ -495,16 +494,16 @@ subroutine evp (dt) call dyn_prep2 (nx_block, ny_block, & ilo, ihi, jlo, jhi, & - icellt (iblk), icelln (iblk), & - indxti (:,iblk), indxtj (:,iblk), & - indxni (:,iblk), indxnj (:,iblk), & + icellT (iblk), icellN (iblk), & + indxTi (:,iblk), indxTj (:,iblk), & + indxNi (:,iblk), indxNj (:,iblk), & aiN (:,:,iblk), nmass (:,:,iblk), & nmassdti (:,:,iblk), fcorN_blk (:,:,iblk), & nmask (:,:,iblk), & uocnN (:,:,iblk), vocnN (:,:,iblk), & strairxN (:,:,iblk), strairyN (:,:,iblk), & ss_tltxN (:,:,iblk), ss_tltyN (:,:,iblk), & - icetmask (:,:,iblk), icenmask (:,:,iblk), & + iceTmask (:,:,iblk), iceNmask (:,:,iblk), & fmN (:,:,iblk), dt, & strtltxN (:,:,iblk), strtltyN (:,:,iblk), & strocnxN (:,:,iblk), strocnyN (:,:,iblk), & @@ -528,16 +527,16 @@ subroutine evp (dt) call dyn_prep2 (nx_block, ny_block, & ilo, ihi, jlo, jhi, & - icellt (iblk), icelle (iblk), & - indxti (:,iblk), indxtj (:,iblk), & - indxei (:,iblk), indxej (:,iblk), & + icellT (iblk), icellE (iblk), & + indxTi (:,iblk), indxTj (:,iblk), & + indxEi (:,iblk), indxEj (:,iblk), & aiE (:,:,iblk), emass (:,:,iblk), & emassdti (:,:,iblk), fcorE_blk (:,:,iblk), & emask (:,:,iblk), & uocnE (:,:,iblk), vocnE (:,:,iblk), & strairxE (:,:,iblk), strairyE (:,:,iblk), & ss_tltxE (:,:,iblk), ss_tltyE (:,:,iblk), & - icetmask (:,:,iblk), iceemask (:,:,iblk), & + iceTmask (:,:,iblk), iceEmask (:,:,iblk), & fmE (:,:,iblk), dt, & strtltxE (:,:,iblk), strtltyE (:,:,iblk), & strocnxE (:,:,iblk), strocnyE (:,:,iblk), & @@ -558,12 +557,12 @@ subroutine evp (dt) do i=1,nx_block do j=1,ny_block - if (.not.iceumask(i,j,iblk)) then + if (.not.iceUmask(i,j,iblk)) then stresspU (i,j,iblk) = c0 stressmU (i,j,iblk) = c0 stress12U(i,j,iblk) = c0 endif - if (icetmask(i,j,iblk) == 0) then + if (.not.iceTmask(i,j,iblk)) then stresspT (i,j,iblk) = c0 stressmT (i,j,iblk) = c0 stress12T(i,j,iblk) = c0 @@ -622,7 +621,7 @@ subroutine evp (dt) if (maskhalo_dyn) then halomask = 0 if (grid_ice == 'B') then - where (iceumask) halomask = 1 + where (iceUmask) halomask = 1 elseif (grid_ice == 'C' .or. grid_ice == 'CD') then !$OMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block,i,j) SCHEDULE(runtime) do iblk = 1, nblocks @@ -633,11 +632,11 @@ subroutine evp (dt) jhi = this_block%jhi do j = jlo,jhi do i = ilo,ihi - if (icetmask(i ,j ,iblk) /= 0 .or. & - icetmask(i-1,j ,iblk) /= 0 .or. & - icetmask(i+1,j ,iblk) /= 0 .or. & - icetmask(i ,j-1,iblk) /= 0 .or. & - icetmask(i ,j+1,iblk) /= 0) then + if (iceTmask(i ,j ,iblk) .or. & + iceTmask(i-1,j ,iblk) .or. & + iceTmask(i+1,j ,iblk) .or. & + iceTmask(i ,j-1,iblk) .or. & + iceTmask(i ,j+1,iblk)) then halomask(i,j,iblk) = 1 endif enddo @@ -664,8 +663,8 @@ subroutine evp (dt) !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime) do iblk = 1, nblocks call seabed_stress_factor_LKD (nx_block , ny_block , & - icellu (iblk), & - indxui (:,iblk), indxuj(:,iblk), & + icellU (iblk), & + indxUi (:,iblk), indxUj(:,iblk), & vice (:,:,iblk), aice(:,:,iblk), & hwater(:,:,iblk), TbU (:,:,iblk)) enddo @@ -675,8 +674,8 @@ subroutine evp (dt) !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime) do iblk = 1, nblocks call seabed_stress_factor_prob (nx_block , ny_block , & - icellt(iblk), indxti(:,iblk), indxtj(:,iblk), & - icellu(iblk), indxui(:,iblk), indxuj(:,iblk), & + icellT(iblk), indxTi(:,iblk), indxTj(:,iblk), & + icellU(iblk), indxUi(:,iblk), indxUj(:,iblk), & aicen(:,:,:,iblk), vicen(:,:,:,iblk) , & hwater (:,:,iblk), TbU (:,:,iblk)) enddo @@ -689,13 +688,13 @@ subroutine evp (dt) !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime) do iblk = 1, nblocks call seabed_stress_factor_LKD (nx_block , ny_block, & - icelle (iblk), & - indxei (:,iblk), indxej(:,iblk), & + icellE (iblk), & + indxEi (:,iblk), indxEj(:,iblk), & vice (:,:,iblk), aice(:,:,iblk), & hwater(:,:,iblk), TbE (:,:,iblk)) call seabed_stress_factor_LKD (nx_block , ny_block, & - icelln (iblk), & - indxni (:,iblk), indxnj(:,iblk), & + icellN (iblk), & + indxNi (:,iblk), indxNj(:,iblk), & vice (:,:,iblk), aice(:,:,iblk), & hwater(:,:,iblk), TbN (:,:,iblk)) enddo @@ -705,13 +704,13 @@ subroutine evp (dt) !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime) do iblk = 1, nblocks call seabed_stress_factor_prob (nx_block , ny_block , & - icellt(iblk), indxti(:,iblk), indxtj(:,iblk), & - icellu(iblk), indxui(:,iblk), indxuj(:,iblk), & + icellT(iblk), indxTi(:,iblk), indxTj(:,iblk), & + icellU(iblk), indxUi(:,iblk), indxUj(:,iblk), & aicen(:,:,:,iblk), vicen(:,:,:,iblk) , & hwater (:,:,iblk), TbU (:,:,iblk) , & TbE (:,:,iblk), TbN (:,:,iblk) , & - icelle(iblk), indxei(:,iblk), indxej(:,iblk), & - icelln(iblk), indxni(:,iblk), indxnj(:,iblk) ) + icellE(iblk), indxEi(:,iblk), indxEj(:,iblk), & + icellN(iblk), indxNi(:,iblk), indxNj(:,iblk) ) enddo !$OMP END PARALLEL DO endif @@ -734,8 +733,8 @@ subroutine evp (dt) call ice_timer_start(timer_evp_1d) call ice_dyn_evp_1d_copyin( & nx_block,ny_block,nblocks,nx_global+2*nghost,ny_global+2*nghost, & - icetmask, iceumask, & - Cdn_ocnU,aiU,uocnU,vocnU,forcexU,forceyU,TbU, & + iceTmask, iceUmask, & + cdn_ocnU,aiU,uocnU,vocnU,forcexU,forceyU,TbU, & umassdti,fmU,uarear,tarear,strintxU,strintyU,uvel_init,vvel_init,& strength,uvel,vvel,dxT,dyT, & stressp_1 ,stressp_2, stressp_3, stressp_4, & @@ -767,8 +766,8 @@ subroutine evp (dt) ! stress tensor equation, total surface stress !----------------------------------------------------------------- call stress (nx_block , ny_block , & - icellt (iblk), & - indxti (:,iblk), indxtj (:,iblk), & + icellT (iblk), & + indxTi (:,iblk), indxTj (:,iblk), & uvel (:,:,iblk), vvel (:,:,iblk), & dxT (:,:,iblk), dyT (:,:,iblk), & dxhy (:,:,iblk), dyhx (:,:,iblk), & @@ -788,8 +787,8 @@ subroutine evp (dt) ! momentum equation !----------------------------------------------------------------- call stepu (nx_block , ny_block , & - icellu (iblk), Cdn_ocnU(:,:,iblk), & - indxui (:,iblk), indxuj (:,iblk), & + icellU (iblk), Cdn_ocnU(:,:,iblk), & + indxUi (:,iblk), indxUj (:,iblk), & aiU (:,:,iblk), strtmp (:,:,:), & uocnU (:,:,iblk), vocnU (:,:,iblk), & waterxU (:,:,iblk), wateryU (:,:,iblk), & @@ -820,8 +819,8 @@ subroutine evp (dt) !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime) do iblk = 1, nblocks call deformations (nx_block , ny_block , & - icellt (iblk), & - indxti (:,iblk), indxtj (:,iblk), & + icellT (iblk), & + indxTi (:,iblk), indxTj (:,iblk), & uvel (:,:,iblk), vvel (:,:,iblk), & dxT (:,:,iblk), dyT (:,:,iblk), & cxp (:,:,iblk), cyp (:,:,iblk), & @@ -845,8 +844,8 @@ subroutine evp (dt) ! NOTE these are actually strain rates * area (m^2/s) !----------------------------------------------------------------- call strain_rates_U (nx_block , ny_block , & - icellu (iblk), & - indxui (:,iblk), indxuj (:,iblk), & + icellU (iblk), & + indxUi (:,iblk), indxUj (:,iblk), & uvelE (:,:,iblk), vvelE (:,:,iblk), & uvelN (:,:,iblk), vvelN (:,:,iblk), & uvel (:,:,iblk), vvel (:,:,iblk), & @@ -869,8 +868,8 @@ subroutine evp (dt) !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks call stressC_T (nx_block , ny_block , & - icellt (iblk), & - indxti (:,iblk), indxtj (:,iblk), & + icellT (iblk), & + indxTi (:,iblk), indxTj (:,iblk), & uvelE (:,:,iblk), vvelE (:,:,iblk), & uvelN (:,:,iblk), vvelN (:,:,iblk), & dxN (:,:,iblk), dyE (:,:,iblk), & @@ -897,8 +896,8 @@ subroutine evp (dt) !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks call stressC_U (nx_block , ny_block , & - icellu (iblk), & - indxui (:,iblk), indxuj (:,iblk), & + icellU (iblk), & + indxUi (:,iblk), indxUj (:,iblk), & uarea (:,:,iblk), & etax2U (:,:,iblk), deltaU (:,:,iblk), & strengthU (:,:,iblk), shearU (:,:,iblk), & @@ -915,8 +914,8 @@ subroutine evp (dt) do iblk = 1, nblocks call div_stress_Ex (nx_block , ny_block , & - icelle (iblk), & - indxei (:,iblk), indxej (:,iblk), & + icellE (iblk), & + indxEi (:,iblk), indxEj (:,iblk), & dxE (:,:,iblk), dyE (:,:,iblk), & dxU (:,:,iblk), dyT (:,:,iblk), & earear (:,:,iblk) , & @@ -924,8 +923,8 @@ subroutine evp (dt) stress12U (:,:,iblk), strintxE (:,:,iblk) ) call div_stress_Ny (nx_block , ny_block , & - icelln (iblk), & - indxni (:,iblk), indxnj (:,iblk), & + icellN (iblk), & + indxNi (:,iblk), indxNj (:,iblk), & dxN (:,:,iblk), dyN (:,:,iblk), & dxT (:,:,iblk), dyU (:,:,iblk), & narear (:,:,iblk) , & @@ -939,8 +938,8 @@ subroutine evp (dt) do iblk = 1, nblocks call stepu_C (nx_block , ny_block , & ! u, E point - icelle (iblk), Cdn_ocnE (:,:,iblk), & - indxei (:,iblk), indxej (:,iblk), & + icellE (iblk), Cdn_ocnE (:,:,iblk), & + indxEi (:,iblk), indxEj (:,iblk), & aiE (:,:,iblk), & uocnE (:,:,iblk), vocnE (:,:,iblk), & waterxE (:,:,iblk), forcexE (:,:,iblk), & @@ -951,8 +950,8 @@ subroutine evp (dt) TbE (:,:,iblk)) call stepv_C (nx_block, ny_block, & ! v, N point - icelln (iblk), Cdn_ocnN (:,:,iblk), & - indxni (:,iblk), indxnj (:,iblk), & + icellN (iblk), Cdn_ocnN (:,:,iblk), & + indxNi (:,iblk), indxNj (:,iblk), & aiN (:,:,iblk), & uocnN (:,:,iblk), vocnN (:,:,iblk), & wateryN (:,:,iblk), forceyN (:,:,iblk), & @@ -1005,8 +1004,8 @@ subroutine evp (dt) !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime) do iblk = 1, nblocks call deformationsC_T (nx_block , ny_block , & - icellt (iblk), & - indxti (:,iblk), indxtj (:,iblk), & + icellT (iblk), & + indxTi (:,iblk), indxTj (:,iblk), & uvelE (:,:,iblk), vvelE (:,:,iblk), & uvelN (:,:,iblk), vvelN (:,:,iblk), & dxN (:,:,iblk), dyE (:,:,iblk), & @@ -1025,8 +1024,8 @@ subroutine evp (dt) !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks call stressCD_T (nx_block , ny_block , & - icellt (iblk), & - indxti (:,iblk), indxtj (:,iblk), & + icellT (iblk), & + indxTi (:,iblk), indxTj (:,iblk), & uvelE (:,:,iblk), vvelE (:,:,iblk), & uvelN (:,:,iblk), vvelN (:,:,iblk), & dxN (:,:,iblk), dyE (:,:,iblk), & @@ -1059,8 +1058,8 @@ subroutine evp (dt) ! NOTE these are actually strain rates * area (m^2/s) !----------------------------------------------------------------- call strain_rates_U (nx_block , ny_block , & - icellu (iblk), & - indxui (:,iblk), indxuj (:,iblk), & + icellU (iblk), & + indxUi (:,iblk), indxUj (:,iblk), & uvelE (:,:,iblk), vvelE (:,:,iblk), & uvelN (:,:,iblk), vvelN (:,:,iblk), & uvel (:,:,iblk), vvel (:,:,iblk), & @@ -1073,8 +1072,8 @@ subroutine evp (dt) shearU (:,:,iblk), DeltaU (:,:,iblk) ) call stressCD_U (nx_block , ny_block , & - icellu (iblk), & - indxui (:,iblk), indxuj (:,iblk), & + icellU (iblk), & + indxUi (:,iblk), indxUj (:,iblk), & uarea (:,:,iblk), & zetax2U (:,:,iblk), etax2U (:,:,iblk), & strengthU(:,:,iblk), & @@ -1097,8 +1096,8 @@ subroutine evp (dt) do iblk = 1, nblocks call div_stress_Ex (nx_block , ny_block , & - icelle (iblk), & - indxei (:,iblk), indxej (:,iblk), & + icellE (iblk), & + indxEi (:,iblk), indxEj (:,iblk), & dxE (:,:,iblk), dyE (:,:,iblk), & dxU (:,:,iblk), dyT (:,:,iblk), & earear (:,:,iblk) , & @@ -1106,8 +1105,8 @@ subroutine evp (dt) stress12U (:,:,iblk), strintxE (:,:,iblk) ) call div_stress_Ey (nx_block , ny_block , & - icelle (iblk), & - indxei (:,iblk), indxej (:,iblk), & + icellE (iblk), & + indxEi (:,iblk), indxEj (:,iblk), & dxE (:,:,iblk), dyE (:,:,iblk), & dxU (:,:,iblk), dyT (:,:,iblk), & earear (:,:,iblk) , & @@ -1115,8 +1114,8 @@ subroutine evp (dt) stress12T (:,:,iblk), strintyE (:,:,iblk) ) call div_stress_Nx (nx_block , ny_block , & - icelln (iblk), & - indxni (:,iblk), indxnj (:,iblk), & + icellN (iblk), & + indxNi (:,iblk), indxNj (:,iblk), & dxN (:,:,iblk), dyN (:,:,iblk), & dxT (:,:,iblk), dyU (:,:,iblk), & narear (:,:,iblk) , & @@ -1124,8 +1123,8 @@ subroutine evp (dt) stress12T (:,:,iblk), strintxN (:,:,iblk) ) call div_stress_Ny (nx_block , ny_block , & - icelln (iblk), & - indxni (:,iblk), indxnj (:,iblk), & + icellN (iblk), & + indxNi (:,iblk), indxNj (:,iblk), & dxN (:,:,iblk), dyN (:,:,iblk), & dxT (:,:,iblk), dyU (:,:,iblk), & narear (:,:,iblk) , & @@ -1139,8 +1138,8 @@ subroutine evp (dt) do iblk = 1, nblocks call stepuv_CD (nx_block , ny_block , & ! E point - icelle (iblk), Cdn_ocnE (:,:,iblk), & - indxei (:,iblk), indxej (:,iblk), & + icellE (iblk), Cdn_ocnE (:,:,iblk), & + indxEi (:,iblk), indxEj (:,iblk), & aiE (:,:,iblk), & uocnE (:,:,iblk), vocnE (:,:,iblk), & waterxE (:,:,iblk), wateryE (:,:,iblk), & @@ -1153,8 +1152,8 @@ subroutine evp (dt) TbE (:,:,iblk)) call stepuv_CD (nx_block , ny_block , & ! N point - icelln (iblk), Cdn_ocnN (:,:,iblk), & - indxni (:,iblk), indxnj (:,iblk), & + icellN (iblk), Cdn_ocnN (:,:,iblk), & + indxNi (:,iblk), indxNj (:,iblk), & aiN (:,:,iblk), & uocnN (:,:,iblk), vocnN (:,:,iblk), & waterxN (:,:,iblk), wateryN (:,:,iblk), & @@ -1196,8 +1195,8 @@ subroutine evp (dt) !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks call deformationsCD_T (nx_block , ny_block , & - icellt (iblk), & - indxti (:,iblk), indxtj (:,iblk), & + icellT (iblk), & + indxTi (:,iblk), indxTj (:,iblk), & uvelE (:,:,iblk), vvelE (:,:,iblk), & uvelN (:,:,iblk), vvelN (:,:,iblk), & dxN (:,:,iblk), dyE (:,:,iblk), & @@ -1299,8 +1298,8 @@ subroutine evp (dt) do iblk = 1, nblocks call dyn_finish & (nx_block , ny_block , & - icellu (iblk), Cdn_ocnU(:,:,iblk), & - indxui (:,iblk), indxuj (:,iblk), & + icellU (iblk), Cdn_ocnU(:,:,iblk), & + indxUi (:,iblk), indxUj (:,iblk), & uvel (:,:,iblk), vvel (:,:,iblk), & uocnU (:,:,iblk), vocnU (:,:,iblk), & aiU (:,:,iblk), fmU (:,:,iblk), & @@ -1315,8 +1314,8 @@ subroutine evp (dt) call dyn_finish & (nx_block , ny_block , & - icelln (iblk), Cdn_ocnN(:,:,iblk), & - indxni (:,iblk), indxnj (:,iblk), & + icellN (iblk), Cdn_ocnN(:,:,iblk), & + indxNi (:,iblk), indxNj (:,iblk), & uvelN (:,:,iblk), vvelN (:,:,iblk), & uocnN (:,:,iblk), vocnN (:,:,iblk), & aiN (:,:,iblk), fmN (:,:,iblk), & @@ -1324,8 +1323,8 @@ subroutine evp (dt) call dyn_finish & (nx_block , ny_block , & - icelle (iblk), Cdn_ocnE(:,:,iblk), & - indxei (:,iblk), indxej (:,iblk), & + icellE (iblk), Cdn_ocnE(:,:,iblk), & + indxEi (:,iblk), indxEj (:,iblk), & uvelE (:,:,iblk), vvelE (:,:,iblk), & uocnE (:,:,iblk), vocnE (:,:,iblk), & aiE (:,:,iblk), fmE (:,:,iblk), & @@ -1358,8 +1357,8 @@ end subroutine evp ! author: Elizabeth C. Hunke, LANL subroutine stress (nx_block, ny_block, & - icellt, & - indxti, indxtj, & + icellT, & + indxTi, indxTj, & uvel, vvel, & dxT, dyT, & dxhy, dyhx, & @@ -1379,11 +1378,11 @@ subroutine stress (nx_block, ny_block, & integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions - icellt ! no. of cells where icetmask = 1 + icellT ! no. of cells where iceTmask = .true. integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: & - indxti , & ! compressed index in i-direction - indxtj ! compressed index in j-direction + indxTi , & ! compressed index in i-direction + indxTj ! compressed index in j-direction real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & strength , & ! ice strength (N/m) @@ -1439,9 +1438,9 @@ subroutine stress (nx_block, ny_block, & str(:,:,:) = c0 - do ij = 1, icellt - i = indxti(ij) - j = indxtj(ij) + do ij = 1, icellT + i = indxTi(ij) + j = indxTj(ij) !----------------------------------------------------------------- ! strain rates @@ -1660,8 +1659,8 @@ end subroutine stress ! for solving the sea ice momentum equation. Ocean Model., 101, 59-67. subroutine stressC_T (nx_block, ny_block , & - icellt , & - indxti , indxtj , & + icellT , & + indxTi , indxTj , & uvelE , vvelE , & uvelN , vvelN , & dxN , dyE , & @@ -1676,11 +1675,11 @@ subroutine stressC_T (nx_block, ny_block , & integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions - icellt ! no. of cells where icetmask = 1 + icellT ! no. of cells where iceTmask = .true. integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: & - indxti , & ! compressed index in i-direction - indxtj ! compressed index in j-direction + indxTi , & ! compressed index in i-direction + indxTj ! compressed index in j-direction real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & uvelE , & ! x-component of velocity (m/s) at the E point @@ -1723,17 +1722,17 @@ subroutine stressC_T (nx_block, ny_block , & !----------------------------------------------------------------- call strain_rates_T (nx_block , ny_block , & - icellt , & - indxti(:) , indxtj (:) , & + icellT , & + indxTi(:) , indxTj (:) , & uvelE (:,:), vvelE (:,:), & uvelN (:,:), vvelN (:,:), & dxN (:,:), dyE (:,:), & dxT (:,:), dyT (:,:), & divT (:,:), tensionT(:,:) ) - do ij = 1, icellt - i = indxti(ij) - j = indxtj(ij) + do ij = 1, icellT + i = indxTi(ij) + j = indxTj(ij) !----------------------------------------------------------------- ! Square of shear strain rate at T obtained from interpolation of @@ -1785,8 +1784,8 @@ end subroutine stressC_T ! for solving the sea ice momentum equation. Ocean Model., 101, 59-67. subroutine stressC_U (nx_block , ny_block, & - icellu, & - indxui , indxuj, & + icellU, & + indxUi , indxUj, & uarea , & etax2U , deltaU, & strengthU, shearU, & @@ -1797,11 +1796,11 @@ subroutine stressC_U (nx_block , ny_block, & integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions - icellu ! no. of cells where iceumask = 1 + icellU ! no. of cells where iceUmask = 1 integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: & - indxui , & ! compressed index in i-direction - indxuj ! compressed index in j-direction + indxUi , & ! compressed index in i-direction + indxUj ! compressed index in j-direction real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & uarea , & ! area of U point @@ -1834,17 +1833,17 @@ subroutine stressC_U (nx_block , ny_block, & !----------------------------------------------------------------- if (visc_method == 'avg_zeta') then - do ij = 1, icellu - i = indxui(ij) - j = indxuj(ij) + do ij = 1, icellU + i = indxUi(ij) + j = indxUj(ij) stress12(i,j) = (stress12(i,j)*(c1-arlx1i*revp) & + arlx1i*p5*etax2U(i,j)*shearU(i,j)) * denom1 enddo elseif (visc_method == 'avg_strength') then - do ij = 1, icellu - i = indxui(ij) - j = indxuj(ij) + do ij = 1, icellU + i = indxUi(ij) + j = indxUj(ij) DminUarea = deltaminEVP*uarea(i,j) ! only need etax2U here, but other terms are calculated with etax2U ! minimal extra calculations here even though it seems like there is @@ -1865,8 +1864,8 @@ end subroutine stressC_U ! Nov 2021 subroutine stressCD_T (nx_block, ny_block, & - icellt, & - indxti, indxtj, & + icellT, & + indxTi, indxTj, & uvelE, vvelE, & uvelN, vvelN, & dxN, dyE, & @@ -1882,11 +1881,11 @@ subroutine stressCD_T (nx_block, ny_block, & integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions - icellt ! no. of cells where icetmask = 1 + icellT ! no. of cells where iceTmask = .true. integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: & - indxti , & ! compressed index in i-direction - indxtj ! compressed index in j-direction + indxTi , & ! compressed index in i-direction + indxTj ! compressed index in j-direction real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & uvelE , & ! x-component of velocity (m/s) at the E point @@ -1929,8 +1928,8 @@ subroutine stressCD_T (nx_block, ny_block, & !----------------------------------------------------------------- call strain_rates_T (nx_block , ny_block , & - icellt , & - indxti(:) , indxtj (:) , & + icellT , & + indxTi(:) , indxTj (:) , & uvelE (:,:), vvelE (:,:), & uvelN (:,:), vvelN (:,:), & dxN (:,:), dyE (:,:), & @@ -1938,9 +1937,9 @@ subroutine stressCD_T (nx_block, ny_block, & divT (:,:), tensionT(:,:), & shearT(:,:), DeltaT (:,:) ) - do ij = 1, icellt - i = indxti(ij) - j = indxtj(ij) + do ij = 1, icellT + i = indxTi(ij) + j = indxTj(ij) !----------------------------------------------------------------- ! viscosities and replacement pressure at T point @@ -1975,8 +1974,8 @@ end subroutine stressCD_T ! Nov 2021 subroutine stressCD_U (nx_block, ny_block, & - icellu, & - indxui, indxuj, & + icellU, & + indxUi, indxUj, & uarea, & zetax2U, etax2U, & strengthU, & @@ -1991,11 +1990,11 @@ subroutine stressCD_U (nx_block, ny_block, & integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions - icellu ! no. of cells where iceumask = 1 + icellU ! no. of cells where iceUmask = 1 integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: & - indxui , & ! compressed index in i-direction - indxuj ! compressed index in j-direction + indxUi , & ! compressed index in i-direction + indxUj ! compressed index in j-direction real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & uarea , & ! area of U-cell (m^2) @@ -2025,9 +2024,9 @@ subroutine stressCD_U (nx_block, ny_block, & character(len=*), parameter :: subname = '(stressCD_U)' - do ij = 1, icellu - i = indxui(ij) - j = indxuj(ij) + do ij = 1, icellU + i = indxUi(ij) + j = indxUj(ij) !----------------------------------------------------------------- ! viscosities and replacement pressure at U point diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp_1d.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp_1d.F90 index fe04a3d63..e874611bd 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp_1d.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp_1d.F90 @@ -1017,7 +1017,7 @@ end subroutine dealloc1d !======================================================================= subroutine ice_dyn_evp_1d_copyin(nx, ny, nblk, nx_glob, ny_glob, & - I_icetmask, I_iceumask, I_cdn_ocn, I_aiu, I_uocn, I_vocn, & + I_iceTmask, I_iceUmask, I_cdn_ocn, I_aiu, I_uocn, I_vocn, & I_forcex, I_forcey, I_Tbu, I_umassdti, I_fm, I_uarear, I_tarear, & I_strintx, I_strinty, I_uvel_init, I_vvel_init, I_strength, & I_uvel, I_vvel, I_dxT, I_dyT, I_stressp_1, I_stressp_2, & @@ -1035,9 +1035,7 @@ subroutine ice_dyn_evp_1d_copyin(nx, ny, nblk, nx_glob, ny_glob, & integer(int_kind), intent(in) :: nx, ny, nblk, nx_glob, ny_glob logical(kind=log_kind), dimension(nx, ny, nblk), intent(in) :: & - I_iceumask - integer(kind=int_kind), dimension(nx, ny, nblk), intent(in) :: & - I_icetmask + I_iceTmask, I_iceUmask real(kind=dbl_kind), dimension(nx, ny, nblk), intent(in) :: & I_cdn_ocn, I_aiu, I_uocn, I_vocn, I_forcex, I_forcey, I_Tbu, & I_umassdti, I_fm, I_uarear, I_tarear, I_strintx, I_strinty, & @@ -1049,9 +1047,7 @@ subroutine ice_dyn_evp_1d_copyin(nx, ny, nblk, nx_glob, ny_glob, & ! local variables logical(kind=log_kind), dimension(nx_glob, ny_glob) :: & - G_iceumask - integer(kind=int_kind), dimension(nx_glob, ny_glob) :: & - G_icetmask + G_iceTmask, G_iceUmask real(kind=dbl_kind), dimension(nx_glob, ny_glob) :: & G_cdn_ocn, G_aiu, G_uocn, G_vocn, G_forcex, G_forcey, G_Tbu, & G_umassdti, G_fm, G_uarear, G_tarear, G_strintx, G_strinty, & @@ -1063,8 +1059,8 @@ subroutine ice_dyn_evp_1d_copyin(nx, ny, nblk, nx_glob, ny_glob, & character(len=*), parameter :: & subname = '(ice_dyn_evp_1d_copyin)' - call gather_global_ext(G_icetmask, I_icetmask, master_task, distrb_info ) - call gather_global_ext(G_iceumask, I_iceumask, master_task, distrb_info ) + call gather_global_ext(G_iceTmask, I_iceTmask, master_task, distrb_info ) + call gather_global_ext(G_iceUmask, I_iceUmask, master_task, distrb_info ) call gather_global_ext(G_cdn_ocn, I_cdn_ocn, master_task, distrb_info ) call gather_global_ext(G_aiu, I_aiu, master_task, distrb_info ) call gather_global_ext(G_uocn, I_uocn, master_task, distrb_info ) @@ -1101,9 +1097,9 @@ subroutine ice_dyn_evp_1d_copyin(nx, ny, nblk, nx_glob, ny_glob, & ! all calculations id done on master task if (my_task == master_task) then ! find number of active points and allocate 1D vectors - call calc_na(nx_glob, ny_glob, NA_len, G_icetmask, G_iceumask) + call calc_na(nx_glob, ny_glob, NA_len, G_iceTmask, G_iceUmask) call alloc1d(NA_len) - call calc_2d_indices(nx_glob, ny_glob, NA_len, G_icetmask, G_iceumask) + call calc_2d_indices(nx_glob, ny_glob, NA_len, G_iceTmask, G_iceUmask) call calc_navel(nx_glob, ny_glob, NA_len, NAVEL_len) call alloc1d_navel(NAVEL_len) ! initialize OpenMP. FIXME: ought to be called from main @@ -1120,7 +1116,7 @@ subroutine ice_dyn_evp_1d_copyin(nx, ny, nblk, nx_glob, ny_glob, & G_stressp_2, G_stressp_3, G_stressp_4, G_stressm_1, & G_stressm_2, G_stressm_3, G_stressm_4, G_stress12_1, & G_stress12_2, G_stress12_3, G_stress12_4) - call calc_halo_parent(nx_glob, ny_glob, NA_len, NAVEL_len, G_icetmask) + call calc_halo_parent(nx_glob, ny_glob, NA_len, NAVEL_len, G_iceTmask) end if end subroutine ice_dyn_evp_1d_copyin @@ -1336,7 +1332,7 @@ end subroutine ice_dyn_evp_1d_kernel !======================================================================= - subroutine calc_na(nx, ny, na, icetmask, iceumask) + subroutine calc_na(nx, ny, na, iceTmask, iceUmask) ! Calculate number of active points use ice_blocks, only : nghost @@ -1344,10 +1340,8 @@ subroutine calc_na(nx, ny, na, icetmask, iceumask) implicit none integer(kind=int_kind), intent(in) :: nx, ny - integer(kind=int_kind), dimension(nx, ny), intent(in) :: & - icetmask logical(kind=log_kind), dimension(nx, ny), intent(in) :: & - iceumask + iceTmask, iceUmask integer(kind=int_kind), intent(out) :: na ! local variables @@ -1360,7 +1354,7 @@ subroutine calc_na(nx, ny, na, icetmask, iceumask) ! NOTE: T mask includes northern and eastern ghost cells do j = 1 + nghost, ny do i = 1 + nghost, nx - if (icetmask(i,j) == 1 .or. iceumask(i,j)) na = na + 1 + if (iceTmask(i,j) .or. iceUmask(i,j)) na = na + 1 end do end do @@ -1368,17 +1362,15 @@ end subroutine calc_na !======================================================================= - subroutine calc_2d_indices(nx, ny, na, icetmask, iceumask) + subroutine calc_2d_indices(nx, ny, na, iceTmask, iceUmask) use ice_blocks, only : nghost implicit none integer(kind=int_kind), intent(in) :: nx, ny, na - integer(kind=int_kind), dimension(nx, ny), intent(in) :: & - icetmask logical(kind=log_kind), dimension(nx, ny), intent(in) :: & - iceumask + iceTmask, iceUmask ! local variables @@ -1394,12 +1386,12 @@ subroutine calc_2d_indices(nx, ny, na, icetmask, iceumask) ! NOTE: T mask includes northern and eastern ghost cells do j = 1 + nghost, ny do i = 1 + nghost, nx - if (icetmask(i,j) == 1 .or. iceumask(i,j)) then + if (iceTmask(i,j) .or. iceUmask(i,j)) then Nmaskt = Nmaskt + 1 indi(Nmaskt) = i indj(Nmaskt) = j - if (icetmask(i,j) /= 1) skiptcell(Nmaskt) = .true. - if (.not. iceumask(i,j)) skipucell(Nmaskt) = .true. + if (.not. iceTmask(i,j)) skiptcell(Nmaskt) = .true. + if (.not. iceUmask(i,j)) skipucell(Nmaskt) = .true. ! NOTE: U mask does not include northern and eastern ! ghost cells. Skip northern and eastern ghost cells if (i == nx) skipucell(Nmaskt) = .true. @@ -1588,13 +1580,13 @@ end subroutine convert_2d_1d !======================================================================= - subroutine calc_halo_parent(nx, ny, na, navel, I_icetmask) + subroutine calc_halo_parent(nx, ny, na, navel, I_iceTmask) implicit none integer(kind=int_kind), intent(in) :: nx, ny, na, navel - integer(kind=int_kind), dimension(nx, ny), intent(in) :: & - I_icetmask + logical(kind=log_kind), dimension(nx, ny), intent(in) :: & + I_iceTmask ! local variables @@ -1619,10 +1611,10 @@ subroutine calc_halo_parent(nx, ny, na, navel, I_icetmask) j = int((indij(iw) - 1) / (nx)) + 1 i = indij(iw) - (j - 1) * nx ! if within ghost zone - if (i == nx .and. I_icetmask(2, j) == 1) Ihalo(iw) = 2 + (j - 1) * nx - if (i == 1 .and. I_icetmask(nx - 1, j) == 1) Ihalo(iw) = (nx - 1) + (j - 1) * nx - if (j == ny .and. I_icetmask(i, 2) == 1) Ihalo(iw) = i + nx - if (j == 1 .and. I_icetmask(i, ny - 1) == 1) Ihalo(iw) = i + (ny - 2) * nx + if (i == nx .and. I_iceTmask(2, j) ) Ihalo(iw) = 2 + (j - 1) * nx + if (i == 1 .and. I_iceTmask(nx - 1, j) ) Ihalo(iw) = (nx - 1) + (j - 1) * nx + if (j == ny .and. I_iceTmask(i, 2) ) Ihalo(iw) = i + nx + if (j == 1 .and. I_iceTmask(i, ny - 1) ) Ihalo(iw) = i + (ny - 2) * nx end do ! relate halo indices to indij vector diff --git a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 index 95d2eedb1..187ec55cc 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -106,6 +106,12 @@ module ice_dyn_shared uvelE_init , & ! x-component of velocity (m/s), beginning of timestep vvelE_init ! y-component of velocity (m/s), beginning of timestep + logical (kind=log_kind), dimension (:,:,:), allocatable, public :: & + iceTmask, & ! ice extent mask (T-cell) + iceUmask, & ! ice extent mask (U-cell) + iceNmask, & ! ice extent mask (N-cell) + iceEmask ! ice extent mask (E-cell) + real (kind=dbl_kind), allocatable, public :: & DminTarea(:,:,:) ! deltamin * tarea (m^2/s) @@ -168,6 +174,8 @@ 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 + iceTmask (nx_block,ny_block,max_blocks), & ! T mask for dynamics + iceUmask (nx_block,ny_block,max_blocks), & ! U mask for dynamics stat=ierr) if (ierr/=0) call abort_ice(subname//': Out of memory') @@ -177,6 +185,8 @@ subroutine alloc_dyn_shared vvelE_init (nx_block,ny_block,max_blocks), & ! y-component of velocity (m/s), beginning of timestep uvelN_init (nx_block,ny_block,max_blocks), & ! x-component of velocity (m/s), beginning of timestep vvelN_init (nx_block,ny_block,max_blocks), & ! y-component of velocity (m/s), beginning of timestep + iceEmask (nx_block,ny_block,max_blocks), & ! T mask for dynamics + iceNmask (nx_block,ny_block,max_blocks), & ! U mask for dynamics stat=ierr) if (ierr/=0) call abort_ice(subname//': Out of memory') endif @@ -199,7 +209,7 @@ subroutine init_dyn (dt) stresspT, stressmT, stress12T, & stresspU, stressmU, stress12U use ice_state, only: uvel, vvel, uvelE, vvelE, uvelN, vvelN, divu, shear - use ice_grid, only: ULAT, NLAT, ELAT, tarea, iceumask, iceemask, icenmask + use ice_grid, only: ULAT, NLAT, ELAT, tarea real (kind=dbl_kind), intent(in) :: & dt ! time step @@ -310,10 +320,10 @@ subroutine init_dyn (dt) endif ! ice extent mask on velocity points - iceumask(i,j,iblk) = .false. + iceUmask(i,j,iblk) = .false. if (grid_ice == 'CD' .or. grid_ice == 'C') then - iceemask(i,j,iblk) = .false. - icenmask(i,j,iblk) = .false. + iceEmask(i,j,iblk) = .false. + iceNmask(i,j,iblk) = .false. end if enddo ! i enddo ! j @@ -394,7 +404,7 @@ subroutine dyn_prep1 (nx_block, ny_block, & real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out) :: & Tmass ! total mass of ice and snow (kg/m^2) - integer (kind=int_kind), dimension (nx_block,ny_block), intent(out) :: & + logical (kind=log_kind), dimension (nx_block,ny_block), intent(out) :: & iceTmask ! ice extent mask (T-cell) ! local variables @@ -438,7 +448,7 @@ subroutine dyn_prep1 (nx_block, ny_block, & !----------------------------------------------------------------- ! augmented mask (land + open ocean) !----------------------------------------------------------------- - iceTmask (i,j) = 0 + iceTmask (i,j) = .false. enddo enddo @@ -450,10 +460,10 @@ subroutine dyn_prep1 (nx_block, ny_block, & if (tmphm(i-1,j+1) .or. tmphm(i,j+1) .or. tmphm(i+1,j+1) .or. & tmphm(i-1,j) .or. tmphm(i,j) .or. tmphm(i+1,j) .or. & tmphm(i-1,j-1) .or. tmphm(i,j-1) .or. tmphm(i+1,j-1) ) then - iceTmask(i,j) = 1 + iceTmask(i,j) = .true. endif - if (.not.Tmask(i,j)) iceTmask(i,j) = 0 + if (.not.Tmask(i,j)) iceTmask(i,j) = .false. enddo enddo @@ -504,8 +514,8 @@ subroutine dyn_prep2 (nx_block, ny_block, & ilo,ihi,jlo,jhi ! beginning and end of physical domain integer (kind=int_kind), intent(out) :: & - icellT , & ! no. of cells where iceTmask = 1 - icellX ! no. of cells where iceXmask = 1 + icellT , & ! no. of cells where iceTmask = .true. + icellX ! no. of cells where iceXmask = .true. integer (kind=int_kind), dimension (nx_block*ny_block), intent(out) :: & indxTi , & ! compressed index in i-direction on T grid @@ -516,7 +526,7 @@ subroutine dyn_prep2 (nx_block, ny_block, & logical (kind=log_kind), dimension (nx_block,ny_block), intent(in) :: & Xmask ! land/boundary mask, thickness (X-grid-cell) - integer (kind=int_kind), dimension (nx_block,ny_block), intent(in) :: & + logical (kind=log_kind), dimension (nx_block,ny_block), intent(in) :: & iceTmask ! ice extent mask (T-cell) logical (kind=log_kind), dimension (nx_block,ny_block), intent(inout) :: & @@ -590,7 +600,7 @@ subroutine dyn_prep2 (nx_block, ny_block, & taubx (i,j) = c0 tauby (i,j) = c0 - if (iceTmask(i,j)==0) then + if (.not.iceTmask(i,j)) then stressp_1 (i,j) = c0 stressp_2 (i,j) = c0 stressp_3 (i,j) = c0 @@ -608,7 +618,7 @@ subroutine dyn_prep2 (nx_block, ny_block, & enddo ! j !----------------------------------------------------------------- - ! Identify cells where iceTmask = 1 + ! Identify cells where iceTmask = .true. ! Note: The icellT mask includes north and east ghost cells ! where stresses are needed. !----------------------------------------------------------------- @@ -616,7 +626,7 @@ subroutine dyn_prep2 (nx_block, ny_block, & icellT = 0 do j = jlo, jhi+1 do i = ilo, ihi+1 - if (iceTmask(i,j) == 1) then + if (iceTmask(i,j)) then icellT = icellT + 1 indxTi(icellT) = i indxTj(icellT) = j @@ -729,7 +739,7 @@ subroutine stepu (nx_block, ny_block, & integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions - icellU ! total count when iceumask is true + icellU ! total count when iceUmask is true integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: & indxUi , & ! compressed index in i-direction @@ -1166,7 +1176,7 @@ subroutine dyn_finish (nx_block, ny_block, & integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions - icellU ! total count when iceumask is true + icellU ! total count when iceUmask is true integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: & indxUi , & ! compressed index in i-direction @@ -1635,7 +1645,7 @@ subroutine deformations (nx_block, ny_block, & integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions - icellT ! no. of cells where iceTmask = 1 + icellT ! no. of cells where iceTmask = .true. integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: & indxTi , & ! compressed index in i-direction @@ -1733,7 +1743,7 @@ subroutine deformationsCD_T (nx_block, ny_block, & integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions - icellT ! no. of cells where iceTmask = 1 + icellT ! no. of cells where iceTmask = .true. integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: & indxTi , & ! compressed index in i-direction @@ -1830,7 +1840,7 @@ subroutine deformationsC_T (nx_block, ny_block, & integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions - icellT ! no. of cells where iceTmask = 1 + icellT ! no. of cells where iceTmask = .true. integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: & indxTi , & ! compressed index in i-direction diff --git a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 index 2a0dbc4b3..caedecc1e 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 @@ -96,14 +96,14 @@ module ice_dyn_vp ! module variables integer (kind=int_kind), allocatable :: & - icellt(:) , & ! no. of cells where icetmask = 1 - icellu(:) ! no. of cells where iceumask = 1 + icellT(:) , & ! no. of cells where iceTmask = .true. + icellU(:) ! no. of cells where iceUmask = .true. integer (kind=int_kind), allocatable :: & - indxti(:,:) , & ! compressed index in i-direction - indxtj(:,:) , & ! compressed index in j-direction - indxui(:,:) , & ! compressed index in i-direction - indxuj(:,:) ! compressed index in j-direction + indxTi(:,:) , & ! compressed index in i-direction + indxTj(:,:) , & ! compressed index in j-direction + indxUi(:,:) , & ! compressed index in i-direction + indxUj(:,:) ! compressed index in j-direction real (kind=dbl_kind), allocatable :: & fld2(:,:,:,:), & ! work array for boundary updates @@ -138,11 +138,11 @@ subroutine init_vp this_block ! block information for current block ! Initialize module variables - allocate(icellt(max_blocks), icellu(max_blocks)) - allocate(indxti(nx_block*ny_block, max_blocks), & - indxtj(nx_block*ny_block, max_blocks), & - indxui(nx_block*ny_block, max_blocks), & - indxuj(nx_block*ny_block, max_blocks)) + allocate(icellT(max_blocks), icellU(max_blocks)) + allocate(indxTi(nx_block*ny_block, max_blocks), & + 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)) allocate(fld3(nx_block,ny_block,3,max_blocks)) allocate(fld4(nx_block,ny_block,4,max_blocks)) @@ -171,7 +171,7 @@ subroutine implicit_solver (dt) use ice_blocks, only: block, get_block, nx_block, ny_block use ice_domain, only: blocks_ice, halo_info, maskhalo_dyn use ice_domain_size, only: max_blocks, ncat - use ice_dyn_shared, only: deformations + use ice_dyn_shared, only: deformations, iceTmask, iceUmask use ice_flux, only: rdg_conv, rdg_shear, strairxT, strairyT, & strairxU, strairyU, uocn, vocn, ss_tltx, ss_tlty, fmU, & strtltxU, strtltyU, strocnxU, strocnyU, strintxU, strintyU, taubxU, taubyU, & @@ -181,7 +181,7 @@ subroutine implicit_solver (dt) stressm_1, stressm_2, stressm_3, stressm_4, & stress12_1, stress12_2, stress12_3, stress12_4 use ice_grid, only: tmask, umask, dxT, dyT, cxp, cyp, cxm, cym, & - tarear, grid_type, grid_average_X2Y, iceumask, & + tarear, grid_type, grid_average_X2Y, & grid_atm_dynu, grid_atm_dynv, grid_ocn_dynu, grid_ocn_dynv use ice_state, only: aice, aiU, vice, vsno, uvel, vvel, divu, shear, & aice_init, aice0, aicen, vicen, strength @@ -226,7 +226,6 @@ subroutine implicit_solver (dt) logical (kind=log_kind) :: calc_strair integer (kind=int_kind), dimension (nx_block,ny_block,max_blocks) :: & - icetmask, & ! ice extent mask (T-cell) halomask ! generic halo mask type (ice_halo) :: & @@ -290,13 +289,13 @@ subroutine implicit_solver (dt) ilo, ihi, jlo, jhi, & aice (:,:,iblk), vice (:,:,iblk), & vsno (:,:,iblk), tmask (:,:,iblk), & - tmass (:,:,iblk), icetmask(:,:,iblk)) + tmass (:,:,iblk), iceTmask(:,:,iblk)) enddo ! iblk !$OMP END PARALLEL DO call ice_timer_start(timer_bound) - call ice_HaloUpdate (icetmask, halo_info, & + call ice_HaloUpdate (iceTmask, halo_info, & field_loc_center, field_type_scalar) call ice_timer_stop(timer_bound) @@ -359,16 +358,16 @@ subroutine implicit_solver (dt) call dyn_prep2 (nx_block, ny_block, & ilo, ihi, jlo, jhi, & - icellt(iblk), icellu(iblk), & - indxti (:,iblk), indxtj (:,iblk), & - indxui (:,iblk), indxuj (:,iblk), & + icellT(iblk), icellU(iblk), & + indxTi (:,iblk), indxTj (:,iblk), & + indxUi (:,iblk), indxUj (:,iblk), & aiU (:,:,iblk), umass (:,:,iblk), & umassdti (:,:,iblk), fcor_blk (:,:,iblk), & umask (:,:,iblk), & uocnU (:,:,iblk), vocnU (:,:,iblk), & strairxU (:,:,iblk), strairyU (:,:,iblk), & ss_tltxU (:,:,iblk), ss_tltyU (:,:,iblk), & - icetmask (:,:,iblk), iceumask (:,:,iblk), & + iceTmask (:,:,iblk), iceUmask (:,:,iblk), & fmU (:,:,iblk), dt, & strtltxU (:,:,iblk), strtltyU (:,:,iblk), & strocnxU (:,:,iblk), strocnyU (:,:,iblk), & @@ -387,8 +386,8 @@ subroutine implicit_solver (dt) TbU (:,:,iblk)) call calc_bfix (nx_block , ny_block , & - icellu(iblk) , & - indxui (:,iblk), indxuj (:,iblk), & + icellU(iblk) , & + indxUi (:,iblk), indxUj (:,iblk), & umassdti (:,:,iblk), & forcexU (:,:,iblk), forceyU (:,:,iblk), & uvel_init (:,:,iblk), vvel_init (:,:,iblk), & @@ -399,9 +398,9 @@ subroutine implicit_solver (dt) !----------------------------------------------------------------- strength(:,:,iblk) = c0 ! initialize - do ij = 1, icellt(iblk) - i = indxti(ij, iblk) - j = indxtj(ij, iblk) + do ij = 1, icellT(iblk) + i = indxTi(ij, iblk) + j = indxTj(ij, iblk) call icepack_ice_strength (ncat, & aice (i,j, iblk), & vice (i,j, iblk), & @@ -431,7 +430,7 @@ subroutine implicit_solver (dt) if (maskhalo_dyn) then call ice_timer_start(timer_bound) halomask = 0 - where (iceumask) halomask = 1 + where (iceUmask) halomask = 1 call ice_HaloUpdate (halomask, halo_info, & field_loc_center, field_type_scalar) call ice_timer_stop(timer_bound) @@ -446,8 +445,8 @@ subroutine implicit_solver (dt) !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks call seabed_stress_factor_LKD (nx_block, ny_block, & - icellu (iblk), & - indxui(:,iblk), indxuj(:,iblk), & + icellU (iblk), & + indxUi(:,iblk), indxUj(:,iblk), & vice(:,:,iblk), aice(:,:,iblk), & hwater(:,:,iblk), TbU(:,:,iblk)) enddo @@ -458,8 +457,8 @@ subroutine implicit_solver (dt) do iblk = 1, nblocks call seabed_stress_factor_prob (nx_block, ny_block, & - icellt(iblk), indxti(:,iblk), indxtj(:,iblk), & - icellu(iblk), indxui(:,iblk), indxuj(:,iblk), & + icellT(iblk), indxTi(:,iblk), indxTj(:,iblk), & + icellU(iblk), indxUi(:,iblk), indxUj(:,iblk), & aicen(:,:,:,iblk), vicen(:,:,:,iblk), & hwater(:,:,iblk), TbU(:,:,iblk)) enddo @@ -475,7 +474,7 @@ subroutine implicit_solver (dt) ntot = 0 do iblk = 1, nblocks - ntot = ntot + icellu(iblk) + ntot = ntot + icellU(iblk) enddo ntot = 2 * ntot ! times 2 because of u and v @@ -484,9 +483,9 @@ subroutine implicit_solver (dt) !----------------------------------------------------------------- ! Start of nonlinear iteration !----------------------------------------------------------------- - call anderson_solver (icellt , icellu , & - indxti , indxtj , & - indxui , indxuj , & + call anderson_solver (icellT , icellU , & + indxTi , indxTj , & + indxUi , indxUj , & aiU , ntot , & uocnU , vocnU , & waterxU , wateryU, & @@ -510,8 +509,8 @@ subroutine implicit_solver (dt) !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks call stress_vp (nx_block , ny_block , & - icellt(iblk) , & - indxti (:,iblk), indxtj (:,iblk), & + icellT(iblk) , & + indxTi (:,iblk), indxTj (:,iblk), & uvel (:,:,iblk), vvel (:,:,iblk), & dxT (:,:,iblk), dyT (:,:,iblk), & cxp (:,:,iblk), cyp (:,:,iblk), & @@ -533,8 +532,8 @@ subroutine implicit_solver (dt) !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks call deformations (nx_block , ny_block , & - icellt(iblk) , & - indxti (:,iblk), indxtj (:,iblk), & + icellT(iblk) , & + indxTi (:,iblk), indxTj (:,iblk), & uvel (:,:,iblk), vvel (:,:,iblk), & dxT (:,:,iblk), dyT (:,:,iblk), & cxp (:,:,iblk), cyp (:,:,iblk), & @@ -552,8 +551,8 @@ subroutine implicit_solver (dt) !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks call calc_seabed_stress (nx_block , ny_block , & - icellu(iblk) , & - indxui (:,iblk), indxuj (:,iblk), & + icellU(iblk) , & + indxUi (:,iblk), indxUj (:,iblk), & uvel (:,:,iblk), vvel (:,:,iblk), & Cb (:,:,iblk), & taubxU (:,:,iblk), taubyU (:,:,iblk)) @@ -638,8 +637,8 @@ subroutine implicit_solver (dt) call dyn_finish & (nx_block, ny_block, & - icellu (iblk), Cdn_ocnU(:,:,iblk), & - indxui (:,iblk), indxuj (:,iblk), & + icellU (iblk), Cdn_ocnU(:,:,iblk), & + indxUi (:,iblk), indxUj (:,iblk), & uvel (:,:,iblk), vvel (:,:,iblk), & uocnU (:,:,iblk), vocnU (:,:,iblk), & aiU (:,:,iblk), fmU (:,:,iblk), & @@ -674,9 +673,9 @@ end subroutine implicit_solver ! H. F. Walker, “Anderson Acceleration: Algorithms and Implementations” ! [Online]. Available: https://users.wpi.edu/~walker/Papers/anderson_accn_algs_imps.pdf - subroutine anderson_solver (icellt , icellu , & - indxti , indxtj , & - indxui , indxuj , & + subroutine anderson_solver (icellT , icellU , & + indxTi , indxTj , & + indxUi , indxUj , & aiU , ntot , & uocn , vocn , & waterxU , wateryU, & @@ -703,14 +702,14 @@ subroutine anderson_solver (icellt , icellu , & ntot ! size of problem for Anderson integer (kind=int_kind), dimension(max_blocks), intent(in) :: & - icellt , & ! no. of cells where icetmask = 1 - icellu ! no. of cells where iceumask = 1 + icellT , & ! no. of cells where iceTmask = .true. + icellU ! no. of cells where iceUmask = .true. integer (kind=int_kind), dimension (nx_block*ny_block, max_blocks), intent(in) :: & - indxti , & ! compressed index in i-direction - indxtj , & ! compressed index in j-direction - indxui , & ! compressed index in i-direction - indxuj ! compressed index in j-direction + indxTi , & ! compressed index in i-direction + indxTj , & ! compressed index in j-direction + indxUi , & ! compressed index in i-direction + indxUj ! compressed index in j-direction real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), intent(in) :: & aiU , & ! ice fraction on u-grid @@ -839,8 +838,8 @@ subroutine anderson_solver (icellt , icellu , & vprev_k(:,:,iblk) = vvel(:,:,iblk) call calc_zeta_dPr (nx_block , ny_block , & - icellt (iblk), & - indxti (:,iblk), indxtj (:,iblk), & + icellT (iblk), & + indxTi (:,iblk), indxTj (:,iblk), & uprev_k (:,:,iblk), vprev_k (:,:,iblk), & dxT (:,:,iblk), dyT (:,:,iblk), & dxhy (:,:,iblk), dyhx (:,:,iblk), & @@ -851,8 +850,8 @@ subroutine anderson_solver (icellt , icellu , & rep_prs(:,:,iblk,:), stress_Pr (:,:,:)) call calc_vrel_Cb (nx_block , ny_block , & - icellu (iblk), Cdn_ocn (:,:,iblk), & - indxui (:,iblk), indxuj (:,iblk), & + icellU (iblk), Cdn_ocn (:,:,iblk), & + indxUi (:,iblk), indxUj (:,iblk), & aiU (:,:,iblk), TbU (:,:,iblk), & uocn (:,:,iblk), vocn (:,:,iblk), & ulin (:,:,iblk), vlin (:,:,iblk), & @@ -860,8 +859,8 @@ subroutine anderson_solver (icellt , icellu , & ! prepare b vector (RHS) call calc_bvec (nx_block , ny_block , & - icellu (iblk), & - indxui (:,iblk), indxuj (:,iblk), & + icellU (iblk), & + indxUi (:,iblk), indxUj (:,iblk), & stress_Pr (:,:,:), uarear (:,:,iblk), & waterxU (:,:,iblk), wateryU (:,:,iblk), & bxfix (:,:,iblk), byfix (:,:,iblk), & @@ -870,9 +869,9 @@ subroutine anderson_solver (icellt , icellu , & ! Compute nonlinear residual norm (PDE residual) call matvec (nx_block , ny_block , & - icellu (iblk) , icellt (iblk), & - indxui (:,iblk) , indxuj (:,iblk), & - indxti (:,iblk) , indxtj (:,iblk), & + icellU (iblk) , icellT (iblk), & + indxUi (:,iblk) , indxUj (:,iblk), & + indxTi (:,iblk) , indxTj (:,iblk), & dxT (:,:,iblk) , dyT (:,:,iblk), & dxhy (:,:,iblk) , dyhx (:,:,iblk), & cxp (:,:,iblk) , cyp (:,:,iblk), & @@ -884,8 +883,8 @@ subroutine anderson_solver (icellt , icellu , & uarear (:,:,iblk) , & Au (:,:,iblk) , Av (:,:,iblk)) call residual_vec (nx_block , ny_block , & - icellu (iblk), & - indxui (:,iblk), indxuj (:,iblk), & + icellU (iblk), & + indxUi (:,iblk), indxUj (:,iblk), & bx (:,:,iblk), by (:,:,iblk), & Au (:,:,iblk), Av (:,:,iblk), & Fx (:,:,iblk), Fy (:,:,iblk), & @@ -912,8 +911,8 @@ subroutine anderson_solver (icellt , icellu , & soly = vprev_k call arrays_to_vec (nx_block , ny_block , & nblocks , max_blocks , & - icellu (:), ntot , & - indxui (:,:), indxuj (:,:), & + icellU (:), ntot , & + indxUi (:,:), indxUj (:,:), & uprev_k (:,:,:), vprev_k (:,:,:), & sol (:)) @@ -927,8 +926,8 @@ subroutine anderson_solver (icellt , icellu , & do iblk = 1, nblocks ! first compute diagonal contributions due to rheology term call formDiag_step1 (nx_block , ny_block , & - icellu (iblk) , & - indxui (:,iblk) , indxuj(:,iblk), & + icellU (iblk) , & + indxUi (:,iblk) , indxUj(:,iblk), & dxT (:,:,iblk) , dyT (:,:,iblk), & dxhy (:,:,iblk) , dyhx(:,:,iblk), & cxp (:,:,iblk) , cyp (:,:,iblk), & @@ -937,8 +936,8 @@ subroutine anderson_solver (icellt , icellu , & diag_rheo(:,:,:)) ! second compute the full diagonal call formDiag_step2 (nx_block , ny_block , & - icellu (iblk), & - indxui (:,iblk), indxuj (:,iblk), & + icellU (iblk), & + indxUi (:,iblk), indxUj (:,iblk), & diag_rheo (:,:,:), vrel (:,:,iblk), & umassdti (:,:,iblk), & uarear (:,:,iblk), Cb (:,:,iblk), & @@ -961,8 +960,8 @@ subroutine anderson_solver (icellt , icellu , & ! Put FGMRES solution solx,soly in fpfunc vector (needed for Anderson) call arrays_to_vec (nx_block , ny_block , & nblocks , max_blocks , & - icellu (:), ntot , & - indxui (:,:), indxuj (:,:), & + icellU (:), ntot , & + indxUi (:,:), indxUj (:,:), & solx (:,:,:), soly (:,:,:), & fpfunc (:)) elseif (fpfunc_andacc == 2) then @@ -978,15 +977,15 @@ subroutine anderson_solver (icellt , icellu , & #else call vec_to_arrays (nx_block , ny_block , & nblocks , max_blocks , & - icellu (:), ntot , & - indxui (:,:), indxuj(:,:) , & + icellU (:), ntot , & + indxUi (:,:), indxUj(:,:) , & res (:), & fpresx (:,:,:), fpresy (:,:,:)) !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks call calc_L2norm_squared (nx_block , ny_block , & - icellu (iblk), & - indxui (:,iblk), indxuj (:,iblk), & + icellU (iblk), & + indxUi (:,iblk), indxUj (:,iblk), & fpresx(:,:,iblk), fpresy(:,:,iblk), & L2norm (iblk)) enddo @@ -1097,8 +1096,8 @@ subroutine anderson_solver (icellt , icellu , & !----------------------------------------------------------------------- call vec_to_arrays (nx_block , ny_block , & nblocks , max_blocks , & - icellu (:), ntot , & - indxui (:,:), indxuj (:,:), & + icellU (:), ntot , & + indxUi (:,:), indxUj (:,:), & sol (:), & uvel (:,:,:), vvel (:,:,:)) @@ -1121,8 +1120,8 @@ subroutine anderson_solver (icellt , icellu , & fpresx(:,:,iblk) = uvel(:,:,iblk) - uprev_k(:,:,iblk) fpresy(:,:,iblk) = vvel(:,:,iblk) - vprev_k(:,:,iblk) call calc_L2norm_squared (nx_block , ny_block , & - icellu (iblk), & - indxui (:,iblk), indxuj (:,iblk), & + icellU (iblk), & + indxUi (:,iblk), indxUj (:,iblk), & fpresx(:,:,iblk), fpresy(:,:,iblk), & L2norm (iblk)) enddo @@ -1142,8 +1141,8 @@ end subroutine anderson_solver ! Computes the viscosities and dPr/dx, dPr/dy subroutine calc_zeta_dPr (nx_block, ny_block, & - icellt , & - indxti , indxtj , & + icellT , & + indxTi , indxTj , & uvel , vvel , & dxT , dyT , & dxhy , dyhx , & @@ -1158,11 +1157,11 @@ subroutine calc_zeta_dPr (nx_block, ny_block, & integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions - icellt ! no. of cells where icetmask = 1 + icellT ! no. of cells where iceTmask = .true. integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: & - indxti , & ! compressed index in i-direction - indxtj ! compressed index in j-direction + indxTi , & ! compressed index in i-direction + indxTj ! compressed index in j-direction real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & strength , & ! ice strength (N/m) @@ -1204,14 +1203,14 @@ subroutine calc_zeta_dPr (nx_block, ny_block, & character(len=*), parameter :: subname = '(calc_zeta_dPr)' ! Initialize stPr, zetax2 and etax2 to zero - ! (for cells where icetmask is false) + ! (for cells where iceTmask is false) stPr = c0 zetax2 = c0 etax2 = c0 - do ij = 1, icellt - i = indxti(ij) - j = indxtj(ij) + do ij = 1, icellT + i = indxTi(ij) + j = indxTj(ij) !----------------------------------------------------------------- ! strain rates @@ -1337,8 +1336,8 @@ end subroutine calc_zeta_dPr ! viscous-plastic sea ice stresses, Geosci. Model Dev., 13, 1763–1769, subroutine stress_vp (nx_block , ny_block , & - icellt , & - indxti , indxtj , & + icellT , & + indxTi , indxTj , & uvel , vvel , & dxT , dyT , & cxp , cyp , & @@ -1356,11 +1355,11 @@ subroutine stress_vp (nx_block , ny_block , & integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions - icellt ! no. of cells where icetmask = 1 + icellT ! no. of cells where iceTmask = .true. integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: & - indxti , & ! compressed index in i-direction - indxtj ! compressed index in j-direction + indxTi , & ! compressed index in i-direction + indxTj ! compressed index in j-direction real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & uvel , & ! x-component of velocity (m/s) @@ -1395,9 +1394,9 @@ subroutine stress_vp (nx_block , ny_block , & character(len=*), parameter :: subname = '(stress_vp)' - do ij = 1, icellt - i = indxti(ij) - j = indxtj(ij) + do ij = 1, icellT + i = indxTi(ij) + j = indxTj(ij) !----------------------------------------------------------------- ! strain rates @@ -1447,8 +1446,8 @@ end subroutine stress_vp ! Compute vrel and seabed stress coefficients subroutine calc_vrel_Cb (nx_block, ny_block, & - icellu , Cw , & - indxui , indxuj , & + icellU , Cw , & + indxUi , indxUj , & aiU , TbU , & uocn , vocn , & uvel , vvel , & @@ -1458,11 +1457,11 @@ subroutine calc_vrel_Cb (nx_block, ny_block, & integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions - icellu ! total count when iceumask is true + icellU ! total count when iceUmask = .true. integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: & - indxui , & ! compressed index in i-direction - indxuj ! compressed index in j-direction + indxUi , & ! compressed index in i-direction + indxUj ! compressed index in j-direction real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & TbU, & ! seabed stress factor (N/m^2) @@ -1494,9 +1493,9 @@ subroutine calc_vrel_Cb (nx_block, ny_block, & if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) - do ij = 1, icellu - i = indxui(ij) - j = indxuj(ij) + do ij = 1, icellU + i = indxUi(ij) + j = indxUj(ij) ! (magnitude of relative ocean current)*rhow*drag*aice vrel(i,j) = aiU(i,j)*rhow*Cw(i,j)*sqrt((uocn(i,j) - uvel(i,j))**2 + & @@ -1512,19 +1511,19 @@ end subroutine calc_vrel_Cb ! Compute seabed stress (diagnostic) subroutine calc_seabed_stress (nx_block, ny_block, & - icellu , & - indxui , indxuj , & + icellU , & + indxUi , indxUj , & uvel , vvel , & Cb , & taubxU , taubyU) integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions - icellu ! total count when iceumask is true + icellU ! total count when iceUmask = .true. integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: & - indxui , & ! compressed index in i-direction - indxuj ! compressed index in j-direction + indxUi , & ! compressed index in i-direction + indxUj ! compressed index in j-direction real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & uvel , & ! x-component of velocity (m/s) @@ -1542,9 +1541,9 @@ subroutine calc_seabed_stress (nx_block, ny_block, & character(len=*), parameter :: subname = '(calc_seabed_stress)' - do ij = 1, icellu - i = indxui(ij) - j = indxuj(ij) + do ij = 1, icellU + i = indxUi(ij) + j = indxUj(ij) taubxU(i,j) = -uvel(i,j)*Cb(i,j) taubyU(i,j) = -vvel(i,j)*Cb(i,j) @@ -1559,9 +1558,9 @@ end subroutine calc_seabed_stress ! Av = A(u,v)_[y] * vvel (y components of A(u,v) * (u,v)) subroutine matvec (nx_block, ny_block, & - icellu , icellt , & - indxui , indxuj , & - indxti , indxtj , & + icellU , icellT , & + indxUi , indxUj , & + indxTi , indxTj , & dxT , dyT , & dxhy , dyhx , & cxp , cyp , & @@ -1577,14 +1576,14 @@ subroutine matvec (nx_block, ny_block, & integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions - icellu, & ! total count when iceumask is true - icellt ! no. of cells where icetmask = 1 + icellU, & ! total count when iceUmask = .true. + icellT ! total count when iceTmask = .true. integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: & - indxui , & ! compressed index in i-direction - indxuj , & ! compressed index in j-direction - indxti , & ! compressed index in i-direction - indxtj ! compressed index in j-direction + indxUi , & ! compressed index in i-direction + indxUj , & ! compressed index in j-direction + indxTi , & ! compressed index in i-direction + indxTj ! compressed index in j-direction real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & dxT , & ! width of T-cell through the middle (m) @@ -1653,9 +1652,9 @@ subroutine matvec (nx_block, ny_block, & str(:,:,:) = c0 - do ij = 1, icellt - i = indxti(ij) - j = indxtj(ij) + do ij = 1, icellT + i = indxTi(ij) + j = indxTj(ij) !----------------------------------------------------------------- ! strain rates @@ -1796,15 +1795,15 @@ subroutine matvec (nx_block, ny_block, & str(i,j,8) = strp_tmp - strm_tmp + str12sn & - dyhx(i,j)*(csigpsw + csigmsw) + dxhy(i,j)*csig12sw - enddo ! ij - icellt + enddo ! ij - icellT !----------------------------------------------------------------- ! Form Au and Av !----------------------------------------------------------------- - do ij = 1, icellu - i = indxui(ij) - j = indxuj(ij) + do ij = 1, icellU + i = indxUi(ij) + j = indxUj(ij) ccaimp = umassdti(i,j) + vrel(i,j) * cosw + Cb(i,j) ! kg/m^2 s @@ -1818,7 +1817,7 @@ subroutine matvec (nx_block, ny_block, & Au(i,j) = ccaimp*uvel(i,j) - ccb*vvel(i,j) - strintx Av(i,j) = ccaimp*vvel(i,j) + ccb*uvel(i,j) - strinty - enddo ! ij - icellu + enddo ! ij - icellU end subroutine matvec @@ -1828,8 +1827,8 @@ end subroutine matvec ! does not depend on (u,v) and thus do not change during the nonlinear iteration subroutine calc_bfix (nx_block , ny_block , & - icellu , & - indxui , indxuj , & + icellU , & + indxUi , indxUj , & umassdti , & forcexU , forceyU , & uvel_init, vvel_init, & @@ -1837,11 +1836,11 @@ subroutine calc_bfix (nx_block , ny_block , & integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions - icellu ! no. of cells where iceumask = 1 + icellU ! no. of cells where iceUmask = .true. integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: & - indxui , & ! compressed index in i-direction - indxuj ! compressed index in j-direction + indxUi , & ! compressed index in i-direction + indxUj ! compressed index in j-direction real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & uvel_init,& ! x-component of velocity (m/s), beginning of time step @@ -1861,9 +1860,9 @@ subroutine calc_bfix (nx_block , ny_block , & character(len=*), parameter :: subname = '(calc_bfix)' - do ij = 1, icellu - i = indxui(ij) - j = indxuj(ij) + do ij = 1, icellU + i = indxUi(ij) + j = indxUj(ij) bxfix(i,j) = umassdti(i,j)*uvel_init(i,j) + forcexU(i,j) byfix(i,j) = umassdti(i,j)*vvel_init(i,j) + forceyU(i,j) @@ -1878,8 +1877,8 @@ end subroutine calc_bfix ! depending on (u,v) subroutine calc_bvec (nx_block, ny_block, & - icellu , & - indxui , indxuj , & + icellU , & + indxUi , indxUj , & stPr , uarear , & waterxU , wateryU , & bxfix , byfix , & @@ -1888,11 +1887,11 @@ subroutine calc_bvec (nx_block, ny_block, & integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions - icellu ! total count when iceumask is true + icellU ! total count when iceUmask = .true. integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: & - indxui , & ! compressed index in i-direction - indxuj ! compressed index in j-direction + indxUi , & ! compressed index in i-direction + indxUj ! compressed index in j-direction real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & uarear , & ! 1/uarea @@ -1930,9 +1929,9 @@ subroutine calc_bvec (nx_block, ny_block, & if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & file=__FILE__, line=__LINE__) - do ij = 1, icellu - i = indxui(ij) - j = indxuj(ij) + do ij = 1, icellU + i = indxUi(ij) + j = indxUj(ij) ! ice/ocean stress taux = vrel(i,j)*waterxU(i,j) ! NOTE this is not the entire @@ -1958,8 +1957,8 @@ end subroutine calc_bvec ! Av = A(u,v)_[y] * vvel (y components of A(u,v) * (u,v)) subroutine residual_vec (nx_block , ny_block, & - icellu , & - indxui , indxuj , & + icellU , & + indxUi , indxUj , & bx , by , & Au , Av , & Fx , Fy , & @@ -1967,11 +1966,11 @@ subroutine residual_vec (nx_block , ny_block, & integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions - icellu ! total count when iceumask is true + icellU ! total count when iceUmask = .true. integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: & - indxui , & ! compressed index in i-direction - indxuj ! compressed index in j-direction + indxUi , & ! compressed index in i-direction + indxUj ! compressed index in j-direction real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & bx , & ! b vector, bx = taux + bxfix (N/m^2) @@ -2001,9 +2000,9 @@ subroutine residual_vec (nx_block , ny_block, & sum_squared = c0 endif - do ij = 1, icellu - i = indxui(ij) - j = indxuj(ij) + do ij = 1, icellU + i = indxUi(ij) + j = indxUj(ij) Fx(i,j) = bx(i,j) - Au(i,j) Fy(i,j) = by(i,j) - Av(i,j) @@ -2020,8 +2019,8 @@ end subroutine residual_vec ! Part 1: compute the contributions to the diagonal from the rheology term subroutine formDiag_step1 (nx_block, ny_block, & - icellu , & - indxui , indxuj , & + icellU , & + indxUi , indxUj , & dxT , dyT , & dxhy , dyhx , & cxp , cyp , & @@ -2031,11 +2030,11 @@ subroutine formDiag_step1 (nx_block, ny_block, & integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions - icellu ! no. of cells where icetmask = 1 + icellU ! no. of cells where iceUmask = .true. integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: & - indxui , & ! compressed index in i-direction - indxuj ! compressed index in j-direction + indxUi , & ! compressed index in i-direction + indxUj ! compressed index in j-direction real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & dxT , & ! width of T-cell through the middle (m) @@ -2191,10 +2190,10 @@ subroutine formDiag_step1 (nx_block, ny_block, & dj = 1 endif - do ij = 1, icellu + do ij = 1, icellU - iu = indxui(ij) - ju = indxuj(ij) + iu = indxUi(ij) + ju = indxUj(ij) i = iu + di j = ju + dj @@ -2395,8 +2394,8 @@ end subroutine formDiag_step1 ! Part 2: compute diagonal subroutine formDiag_step2 (nx_block, ny_block, & - icellu , & - indxui , indxuj , & + icellU , & + indxUi , indxUj , & Drheo , vrel , & umassdti, & uarear , Cb , & @@ -2404,11 +2403,11 @@ subroutine formDiag_step2 (nx_block, ny_block, & integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions - icellu ! total count when iceumask is true + icellU ! total count when iceUmask = .true. integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: & - indxui , & ! compressed index in i-direction - indxuj ! compressed index in j-direction + indxUi , & ! compressed index in i-direction + indxUj ! compressed index in j-direction real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & vrel, & ! coefficient for tauw @@ -2454,9 +2453,9 @@ subroutine formDiag_step2 (nx_block, ny_block, & ! Drheo(i,j,7) corresponds to str(i+1,j,7) ! Drheo(i,j,8) corresponds to str(i+1,j+1,8)) - do ij = 1, icellu - i = indxui(ij) - j = indxuj(ij) + do ij = 1, icellU + i = indxUi(ij) + j = indxUj(ij) ccaimp = umassdti(i,j) + vrel(i,j) * cosw + Cb(i,j) ! kg/m^2 s @@ -2476,18 +2475,18 @@ end subroutine formDiag_step2 ! Compute squared l^2 norm of a grid function (tpu,tpv) subroutine calc_L2norm_squared (nx_block, ny_block, & - icellu , & - indxui , indxuj , & + icellU , & + indxUi , indxUj , & tpu , tpv , & L2norm) integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions - icellu ! total count when iceumask is true + icellU ! total count when iceUmask = .true. integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: & - indxui , & ! compressed index in i-direction - indxuj ! compressed index in j-direction + indxUi , & ! compressed index in i-direction + indxUj ! compressed index in j-direction real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: & tpu , & ! x-component of vector grid function @@ -2509,9 +2508,9 @@ subroutine calc_L2norm_squared (nx_block, ny_block, & L2norm = c0 - do ij = 1, icellu - i = indxui(ij) - j = indxuj(ij) + do ij = 1, icellU + i = indxUi(ij) + j = indxUj(ij) L2norm = L2norm + tpu(i,j)**2 + tpv(i,j)**2 enddo ! ij @@ -2524,8 +2523,8 @@ end subroutine calc_L2norm_squared subroutine arrays_to_vec (nx_block, ny_block , & nblocks , max_blocks, & - icellu , ntot , & - indxui , indxuj , & + icellU , ntot , & + indxUi , indxUj , & tpu , tpv , & outvec) @@ -2536,11 +2535,11 @@ subroutine arrays_to_vec (nx_block, ny_block , & ntot ! size of problem for Anderson integer (kind=int_kind), dimension (max_blocks), intent(in) :: & - icellu + icellU integer (kind=int_kind), dimension (nx_block*ny_block, max_blocks), intent(in) :: & - indxui , & ! compressed index in i-direction - indxuj ! compressed index in j-direction + indxUi , & ! compressed index in i-direction + indxUj ! compressed index in j-direction real (kind=dbl_kind), dimension (nx_block,ny_block, max_blocks), intent(in) :: & tpu , & ! x-component of vector @@ -2564,9 +2563,9 @@ subroutine arrays_to_vec (nx_block, ny_block , & tot = 0 do iblk = 1, nblocks - do ij = 1, icellu(iblk) - i = indxui(ij, iblk) - j = indxuj(ij, iblk) + do ij = 1, icellU(iblk) + i = indxUi(ij, iblk) + j = indxUj(ij, iblk) tot = tot + 1 outvec(tot) = tpu(i, j, iblk) tot = tot + 1 @@ -2582,8 +2581,8 @@ end subroutine arrays_to_vec subroutine vec_to_arrays (nx_block, ny_block , & nblocks , max_blocks, & - icellu , ntot , & - indxui , indxuj , & + icellU , ntot , & + indxUi , indxUj , & invec , & tpu , tpv) @@ -2594,11 +2593,11 @@ subroutine vec_to_arrays (nx_block, ny_block , & ntot ! size of problem for Anderson integer (kind=int_kind), dimension (max_blocks), intent(in) :: & - icellu + icellU integer (kind=int_kind), dimension (nx_block*ny_block, max_blocks), intent(in) :: & - indxui , & ! compressed index in i-direction - indxuj ! compressed index in j-direction + indxUi , & ! compressed index in i-direction + indxUj ! compressed index in j-direction real (kind=dbl_kind), dimension (ntot), intent(in) :: & invec ! input 1D vector @@ -2623,9 +2622,9 @@ subroutine vec_to_arrays (nx_block, ny_block , & tot = 0 do iblk = 1, nblocks - do ij = 1, icellu(iblk) - i = indxui(ij, iblk) - j = indxuj(ij, iblk) + do ij = 1, icellU(iblk) + i = indxUi(ij, iblk) + j = indxUj(ij, iblk) tot = tot + 1 tpu(i, j, iblk) = invec(tot) tot = tot + 1 @@ -2813,9 +2812,9 @@ subroutine fgmres (zetax2 , etax2 , & !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks call matvec (nx_block , ny_block , & - icellu (iblk) , icellt (iblk), & - indxui (:,iblk) , indxuj (:,iblk), & - indxti (:,iblk) , indxtj (:,iblk), & + icellU (iblk) , icellT (iblk), & + indxUi (:,iblk) , indxUj (:,iblk), & + indxTi (:,iblk) , indxTj (:,iblk), & dxT (:,:,iblk) , dyT (:,:,iblk), & dxhy (:,:,iblk) , dyhx (:,:,iblk), & cxp (:,:,iblk) , cyp (:,:,iblk), & @@ -2827,8 +2826,8 @@ subroutine fgmres (zetax2 , etax2 , & uarear (:,:,iblk) , & workspace_x(:,:,iblk) , workspace_y(:,:,iblk)) call residual_vec (nx_block , ny_block , & - icellu (iblk), & - indxui (:,iblk), indxuj (:,iblk), & + icellU (iblk), & + indxUi (:,iblk), indxUj (:,iblk), & bx (:,:,iblk), by (:,:,iblk), & workspace_x(:,:,iblk), workspace_y(:,:,iblk), & arnoldi_basis_x (:,:,iblk, 1), & @@ -2842,8 +2841,8 @@ subroutine fgmres (zetax2 , etax2 , & !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks call calc_L2norm_squared(nx_block , ny_block , & - icellu (iblk), & - indxui (:,iblk), indxuj(:,iblk), & + icellU (iblk), & + indxUi (:,iblk), indxUj(:,iblk), & arnoldi_basis_x(:,:,iblk, 1) , & arnoldi_basis_y(:,:,iblk, 1) , & norm_squared(iblk)) @@ -2866,9 +2865,9 @@ subroutine fgmres (zetax2 , etax2 , & inverse_norm = c1 / norm_residual !$OMP PARALLEL DO PRIVATE(iblk,ij,i,j) do iblk = 1, nblocks - do ij = 1, icellu(iblk) - i = indxui(ij, iblk) - j = indxuj(ij, iblk) + do ij = 1, icellU(iblk) + i = indxUi(ij, iblk) + j = indxUj(ij, iblk) arnoldi_basis_x(i, j, iblk, 1) = arnoldi_basis_x(i, j, iblk, 1) * inverse_norm arnoldi_basis_y(i, j, iblk, 1) = arnoldi_basis_y(i, j, iblk, 1) * inverse_norm @@ -2920,9 +2919,9 @@ subroutine fgmres (zetax2 , etax2 , & !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks call matvec (nx_block , ny_block , & - icellu (iblk) , icellt (iblk), & - indxui (:,iblk) , indxuj (:,iblk), & - indxti (:,iblk) , indxtj (:,iblk), & + icellU (iblk) , icellT (iblk), & + indxUi (:,iblk) , indxUj (:,iblk), & + indxTi (:,iblk) , indxTj (:,iblk), & dxT (:,:,iblk) , dyT (:,:,iblk), & dxhy (:,:,iblk) , dyhx (:,:,iblk), & cxp (:,:,iblk) , cyp (:,:,iblk), & @@ -2947,8 +2946,8 @@ subroutine fgmres (zetax2 , etax2 , & !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks call calc_L2norm_squared(nx_block , ny_block , & - icellu (iblk), & - indxui (:,iblk), indxuj(:, iblk) , & + icellU (iblk), & + indxUi (:,iblk), indxUj(:, iblk) , & arnoldi_basis_x(:,:,iblk, nextit), & arnoldi_basis_y(:,:,iblk, nextit), & norm_squared(iblk)) @@ -2962,9 +2961,9 @@ subroutine fgmres (zetax2 , etax2 , & inverse_norm = c1 / hessenberg(nextit,initer) !$OMP PARALLEL DO PRIVATE(iblk,ij,i,j) do iblk = 1, nblocks - do ij = 1, icellu(iblk) - i = indxui(ij, iblk) - j = indxuj(ij, iblk) + do ij = 1, icellU(iblk) + i = indxUi(ij, iblk) + j = indxUj(ij, iblk) arnoldi_basis_x(i, j, iblk, nextit) = arnoldi_basis_x(i, j, iblk, nextit)*inverse_norm arnoldi_basis_y(i, j, iblk, nextit) = arnoldi_basis_y(i, j, iblk, nextit)*inverse_norm @@ -3028,9 +3027,9 @@ subroutine fgmres (zetax2 , etax2 , & t = rhs_hess(it) !$OMP PARALLEL DO PRIVATE(iblk,ij,i,j) do iblk = 1, nblocks - do ij = 1, icellu(iblk) - i = indxui(ij, iblk) - j = indxuj(ij, iblk) + do ij = 1, icellU(iblk) + i = indxUi(ij, iblk) + j = indxUj(ij, iblk) solx(i, j, iblk) = solx(i, j, iblk) + t * orig_basis_x(i, j, iblk, it) soly(i, j, iblk) = soly(i, j, iblk) + t * orig_basis_y(i, j, iblk, it) @@ -3072,9 +3071,9 @@ subroutine fgmres (zetax2 , etax2 , & do it = 1, nextit !$OMP PARALLEL DO PRIVATE(iblk,ij,i,j) do iblk = 1, nblocks - do ij = 1, icellu(iblk) - i = indxui(ij, iblk) - j = indxuj(ij, iblk) + do ij = 1, icellU(iblk) + i = indxUi(ij, iblk) + j = indxUj(ij, iblk) workspace_x(i, j, iblk) = workspace_x(i, j, iblk) + rhs_hess(it) * arnoldi_basis_x(i, j, iblk, it) workspace_y(i, j, iblk) = workspace_x(i, j, iblk) + rhs_hess(it) * arnoldi_basis_y(i, j, iblk, it) @@ -3206,9 +3205,9 @@ subroutine pgmres (zetax2 , etax2 , & !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks call matvec (nx_block , ny_block , & - icellu (iblk) , icellt (iblk), & - indxui (:,iblk) , indxuj (:,iblk), & - indxti (:,iblk) , indxtj (:,iblk), & + icellU (iblk) , icellT (iblk), & + indxUi (:,iblk) , indxUj (:,iblk), & + indxTi (:,iblk) , indxTj (:,iblk), & dxT (:,:,iblk) , dyT (:,:,iblk), & dxhy (:,:,iblk) , dyhx (:,:,iblk), & cxp (:,:,iblk) , cyp (:,:,iblk), & @@ -3220,8 +3219,8 @@ subroutine pgmres (zetax2 , etax2 , & uarear (:,:,iblk) , & workspace_x(:,:,iblk) , workspace_y(:,:,iblk)) call residual_vec (nx_block , ny_block , & - icellu (iblk), & - indxui (:,iblk), indxuj (:,iblk), & + icellU (iblk), & + indxUi (:,iblk), indxUj (:,iblk), & bx (:,:,iblk), by (:,:,iblk), & workspace_x(:,:,iblk), workspace_y(:,:,iblk), & arnoldi_basis_x (:,:,iblk, 1), & @@ -3235,8 +3234,8 @@ subroutine pgmres (zetax2 , etax2 , & !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks call calc_L2norm_squared(nx_block , ny_block , & - icellu (iblk), & - indxui (:,iblk), indxuj(:, iblk), & + icellU (iblk), & + indxUi (:,iblk), indxUj(:, iblk), & arnoldi_basis_x(:,:,iblk, 1), & arnoldi_basis_y(:,:,iblk, 1), & norm_squared(iblk)) @@ -3259,9 +3258,9 @@ subroutine pgmres (zetax2 , etax2 , & inverse_norm = c1 / norm_residual !$OMP PARALLEL DO PRIVATE(iblk,ij,i,j) do iblk = 1, nblocks - do ij = 1, icellu(iblk) - i = indxui(ij, iblk) - j = indxuj(ij, iblk) + do ij = 1, icellU(iblk) + i = indxUi(ij, iblk) + j = indxUj(ij, iblk) arnoldi_basis_x(i, j, iblk, 1) = arnoldi_basis_x(i, j, iblk, 1) * inverse_norm arnoldi_basis_y(i, j, iblk, 1) = arnoldi_basis_y(i, j, iblk, 1) * inverse_norm @@ -3302,9 +3301,9 @@ subroutine pgmres (zetax2 , etax2 , & !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks call matvec (nx_block , ny_block , & - icellu (iblk) , icellt (iblk), & - indxui (:,iblk) , indxuj (:,iblk), & - indxti (:,iblk) , indxtj (:,iblk), & + icellU (iblk) , icellT (iblk), & + indxUi (:,iblk) , indxUj (:,iblk), & + indxTi (:,iblk) , indxTj (:,iblk), & dxT (:,:,iblk) , dyT (:,:,iblk), & dxhy (:,:,iblk) , dyhx (:,:,iblk), & cxp (:,:,iblk) , cyp (:,:,iblk), & @@ -3329,8 +3328,8 @@ subroutine pgmres (zetax2 , etax2 , & !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks call calc_L2norm_squared(nx_block , ny_block , & - icellu (iblk), & - indxui (:,iblk), indxuj(:, iblk) , & + icellU (iblk), & + indxUi (:,iblk), indxUj(:, iblk) , & arnoldi_basis_x(:,:,iblk, nextit), & arnoldi_basis_y(:,:,iblk, nextit), & norm_squared(iblk)) @@ -3344,9 +3343,9 @@ subroutine pgmres (zetax2 , etax2 , & inverse_norm = c1 / hessenberg(nextit,initer) !$OMP PARALLEL DO PRIVATE(iblk,ij,i,j) do iblk = 1, nblocks - do ij = 1, icellu(iblk) - i = indxui(ij, iblk) - j = indxuj(ij, iblk) + do ij = 1, icellU(iblk) + i = indxUi(ij, iblk) + j = indxUj(ij, iblk) arnoldi_basis_x(i, j, iblk, nextit) = arnoldi_basis_x(i, j, iblk, nextit)*inverse_norm arnoldi_basis_y(i, j, iblk, nextit) = arnoldi_basis_y(i, j, iblk, nextit)*inverse_norm @@ -3412,9 +3411,9 @@ subroutine pgmres (zetax2 , etax2 , & t = rhs_hess(it) !$OMP PARALLEL DO PRIVATE(iblk,ij,i,j) do iblk = 1, nblocks - do ij = 1, icellu(iblk) - i = indxui(ij, iblk) - j = indxuj(ij, iblk) + do ij = 1, icellU(iblk) + i = indxUi(ij, iblk) + j = indxUj(ij, iblk) workspace_x(i, j, iblk) = workspace_x(i, j, iblk) + t * arnoldi_basis_x(i, j, iblk, it) workspace_y(i, j, iblk) = workspace_y(i, j, iblk) + t * arnoldi_basis_y(i, j, iblk, it) @@ -3468,9 +3467,9 @@ subroutine pgmres (zetax2 , etax2 , & do it = 1, nextit !$OMP PARALLEL DO PRIVATE(iblk,ij,i,j) do iblk = 1, nblocks - do ij = 1, icellu(iblk) - i = indxui(ij, iblk) - j = indxuj(ij, iblk) + do ij = 1, icellU(iblk) + i = indxUi(ij, iblk) + j = indxUj(ij, iblk) workspace_x(i, j, iblk) = workspace_x(i, j, iblk) + rhs_hess(it) * arnoldi_basis_x(i, j, iblk, it) workspace_y(i, j, iblk) = workspace_x(i, j, iblk) + rhs_hess(it) * arnoldi_basis_y(i, j, iblk, it) @@ -3549,9 +3548,9 @@ subroutine precondition(zetax2 , etax2, & elseif (precond_type == 'diag') then ! Jacobi preconditioner (diagonal) !$OMP PARALLEL DO PRIVATE(iblk,ij,i,j) do iblk = 1, nblocks - do ij = 1, icellu(iblk) - i = indxui(ij, iblk) - j = indxuj(ij, iblk) + do ij = 1, icellU(iblk) + i = indxUi(ij, iblk) + j = indxUj(ij, iblk) wx(i,j,iblk) = vx(i,j,iblk) / diagx(i,j,iblk) wy(i,j,iblk) = vy(i,j,iblk) / diagy(i,j,iblk) @@ -3632,9 +3631,9 @@ subroutine orthogonalize(ortho_type , initer , & !$OMP PARALLEL DO PRIVATE(iblk,ij,i,j) do iblk = 1, nblocks - do ij = 1, icellu(iblk) - i = indxui(ij, iblk) - j = indxuj(ij, iblk) + do ij = 1, icellU(iblk) + i = indxUi(ij, iblk) + j = indxUj(ij, iblk) local_dot(iblk) = local_dot(iblk) + & (arnoldi_basis_x(i, j, iblk, it) * arnoldi_basis_x(i, j, iblk, nextit)) + & @@ -3652,9 +3651,9 @@ subroutine orthogonalize(ortho_type , initer , & do it = 1, initer !$OMP PARALLEL DO PRIVATE(iblk,ij,i,j) do iblk = 1, nblocks - do ij = 1, icellu(iblk) - i = indxui(ij, iblk) - j = indxuj(ij, iblk) + do ij = 1, icellU(iblk) + i = indxUi(ij, iblk) + j = indxUj(ij, iblk) arnoldi_basis_x(i, j, iblk, nextit) = arnoldi_basis_x(i, j, iblk, nextit) & - hessenberg(it,initer) * arnoldi_basis_x(i, j, iblk, it) @@ -3671,9 +3670,9 @@ subroutine orthogonalize(ortho_type , initer , & !$OMP PARALLEL DO PRIVATE(iblk,ij,i,j) do iblk = 1, nblocks - do ij = 1, icellu(iblk) - i = indxui(ij, iblk) - j = indxuj(ij, iblk) + do ij = 1, icellU(iblk) + i = indxUi(ij, iblk) + j = indxUj(ij, iblk) local_dot(iblk) = local_dot(iblk) + & (arnoldi_basis_x(i, j, iblk, it) * arnoldi_basis_x(i, j, iblk, nextit)) + & @@ -3686,9 +3685,9 @@ subroutine orthogonalize(ortho_type , initer , & !$OMP PARALLEL DO PRIVATE(iblk,ij,i,j) do iblk = 1, nblocks - do ij = 1, icellu(iblk) - i = indxui(ij, iblk) - j = indxuj(ij, iblk) + do ij = 1, icellU(iblk) + i = indxUi(ij, iblk) + j = indxUj(ij, iblk) arnoldi_basis_x(i, j, iblk, nextit) = arnoldi_basis_x(i, j, iblk, nextit) & - hessenberg(it,initer) * arnoldi_basis_x(i, j, iblk, it) diff --git a/cicecore/cicedynB/infrastructure/comm/mpi/ice_boundary.F90 b/cicecore/cicedynB/infrastructure/comm/mpi/ice_boundary.F90 index 2b64f8932..9fda67dad 100644 --- a/cicecore/cicedynB/infrastructure/comm/mpi/ice_boundary.F90 +++ b/cicecore/cicedynB/infrastructure/comm/mpi/ice_boundary.F90 @@ -81,6 +81,7 @@ module ice_boundary module procedure ice_HaloUpdate2DR8, & ice_HaloUpdate2DR4, & ice_HaloUpdate2DI4, & + ice_HaloUpdate2DL1, & ice_HaloUpdate3DR8, & ice_HaloUpdate3DR4, & ice_HaloUpdate3DI4, & @@ -2384,6 +2385,69 @@ subroutine ice_HaloUpdate2DI4(array, halo, & end subroutine ice_HaloUpdate2DI4 +!*********************************************************************** + + subroutine ice_HaloUpdate2DL1(array, halo, & + fieldLoc, fieldKind, & + fillValue) + +! This routine updates ghost cells for an input array and is a +! member of a group of routines under the generic interface +! ice\_HaloUpdate. This routine is the specific interface +! for 2d horizontal logical arrays. + + type (ice_halo), intent(in) :: & + halo ! precomputed halo structure containing all + ! information needed for halo update + + integer (int_kind), intent(in) :: & + fieldKind, &! id for type of field (scalar, vector, angle) + fieldLoc ! id for location on horizontal grid + ! (center, NEcorner, Nface, Eface) + + integer (int_kind), intent(in), optional :: & + fillValue ! optional value to put in ghost cells + ! where neighbor points are unknown + ! (e.g. eliminated land blocks or + ! closed boundaries) + + logical (log_kind), dimension(:,:,:), intent(inout) :: & + array ! array containing field for which halo + ! needs to be updated + +!----------------------------------------------------------------------- +! +! local variables +! +!----------------------------------------------------------------------- + + integer (int_kind), dimension(:,:,:), allocatable :: & + iarray ! integer array for logical + + character(len=*), parameter :: subname = '(ice_HaloUpdate2DL1)' + +!----------------------------------------------------------------------- +! +! copy logical into integer array and call haloupdate on integer array +! +!----------------------------------------------------------------------- + + allocate(iarray(size(array,dim=1),size(array,dim=2),size(array,dim=3))) + iarray(:,:,:) = 0 + where (array) iarray = 1 + + call ice_HaloUpdate(iarray, halo, & + fieldLoc, fieldKind, & + fillValue) + + array = .false. + where (iarray /= 0) array = .true. + deallocate(iarray) + +!----------------------------------------------------------------------- + + end subroutine ice_HaloUpdate2DL1 + !*********************************************************************** subroutine ice_HaloUpdate3DR8(array, halo, & diff --git a/cicecore/cicedynB/infrastructure/comm/serial/ice_boundary.F90 b/cicecore/cicedynB/infrastructure/comm/serial/ice_boundary.F90 index cafe4dc05..f10a9f432 100644 --- a/cicecore/cicedynB/infrastructure/comm/serial/ice_boundary.F90 +++ b/cicecore/cicedynB/infrastructure/comm/serial/ice_boundary.F90 @@ -68,6 +68,7 @@ module ice_boundary module procedure ice_HaloUpdate2DR8, & ice_HaloUpdate2DR4, & ice_HaloUpdate2DI4, & + ice_HaloUpdate2DL1, & ice_HaloUpdate3DR8, & ice_HaloUpdate3DR4, & ice_HaloUpdate3DI4, & @@ -1541,6 +1542,69 @@ subroutine ice_HaloUpdate2DI4(array, halo, & end subroutine ice_HaloUpdate2DI4 +!*********************************************************************** + + subroutine ice_HaloUpdate2DL1(array, halo, & + fieldLoc, fieldKind, & + fillValue) + +! This routine updates ghost cells for an input array and is a +! member of a group of routines under the generic interface +! ice\_HaloUpdate. This routine is the specific interface +! for 2d horizontal logical arrays. + + type (ice_halo), intent(in) :: & + halo ! precomputed halo structure containing all + ! information needed for halo update + + integer (int_kind), intent(in) :: & + fieldKind, &! id for type of field (scalar, vector, angle) + fieldLoc ! id for location on horizontal grid + ! (center, NEcorner, Nface, Eface) + + integer (int_kind), intent(in), optional :: & + fillValue ! optional value to put in ghost cells + ! where neighbor points are unknown + ! (e.g. eliminated land blocks or + ! closed boundaries) + + logical (log_kind), dimension(:,:,:), intent(inout) :: & + array ! array containing field for which halo + ! needs to be updated + +!----------------------------------------------------------------------- +! +! local variables +! +!----------------------------------------------------------------------- + + integer (int_kind), dimension(:,:,:), allocatable :: & + iarray ! integer array for logical + + character(len=*), parameter :: subname = '(ice_HaloUpdate2DL1)' + +!----------------------------------------------------------------------- +! +! copy logical into integer array and call haloupdate on integer array +! +!----------------------------------------------------------------------- + + allocate(iarray(size(array,dim=1),size(array,dim=2),size(array,dim=3))) + iarray(:,:,:) = 0 + where (array) iarray = 1 + + call ice_HaloUpdate(iarray, halo, & + fieldLoc, fieldKind, & + fillValue) + + array = .false. + where (iarray /= 0) array = .true. + deallocate(iarray) + +!----------------------------------------------------------------------- + + end subroutine ice_HaloUpdate2DL1 + !*********************************************************************** subroutine ice_HaloUpdate3DR8(array, halo, & diff --git a/cicecore/cicedynB/infrastructure/ice_grid.F90 b/cicecore/cicedynB/infrastructure/ice_grid.F90 index a7d8c2357..b7082bf93 100644 --- a/cicecore/cicedynB/infrastructure/ice_grid.F90 +++ b/cicecore/cicedynB/infrastructure/ice_grid.F90 @@ -192,11 +192,6 @@ module ice_grid lmask_n, & ! northern hemisphere mask lmask_s ! southern hemisphere mask - logical (kind=log_kind), dimension (:,:,:), allocatable, public :: & - iceumask, & ! ice extent mask (U-cell) - icenmask, & ! ice extent mask (N-cell) - iceemask ! ice extent mask (E-cell) - real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: & rndex_global ! global index for local subdomain (dbl) @@ -277,7 +272,6 @@ subroutine alloc_grid umaskCD (nx_block,ny_block,max_blocks), & ! land/boundary mask, velocity (U-cell) nmask (nx_block,ny_block,max_blocks), & ! land/boundary mask (N-cell) emask (nx_block,ny_block,max_blocks), & ! land/boundary mask (E-cell) - iceumask (nx_block,ny_block,max_blocks), & ! u mask for dynamics lmask_n (nx_block,ny_block,max_blocks), & ! northern hemisphere mask lmask_s (nx_block,ny_block,max_blocks), & ! southern hemisphere mask rndex_global(nx_block,ny_block,max_blocks), & ! global index for local subdomain (dbl) @@ -298,8 +292,6 @@ subroutine alloc_grid if (grid_ice == 'CD' .or. grid_ice == 'C') then allocate( & - iceemask (nx_block,ny_block,max_blocks), & ! e mask for dynamics - icenmask (nx_block,ny_block,max_blocks), & ! n mask for dynamics ratiodxN (nx_block,ny_block,max_blocks), & ratiodyE (nx_block,ny_block,max_blocks), & ratiodxNr(nx_block,ny_block,max_blocks), & diff --git a/cicecore/cicedynB/infrastructure/ice_restart_driver.F90 b/cicecore/cicedynB/infrastructure/ice_restart_driver.F90 index cfc44d987..2c7d3d63c 100644 --- a/cicecore/cicedynB/infrastructure/ice_restart_driver.F90 +++ b/cicecore/cicedynB/infrastructure/ice_restart_driver.F90 @@ -55,6 +55,7 @@ subroutine dumpfile(filename_spec) use ice_blocks, only: nx_block, ny_block use ice_domain, only: nblocks use ice_domain_size, only: nilyr, nslyr, ncat, max_blocks + use ice_dyn_shared, only: iceUmask, iceEmask, iceNmask use ice_flux, only: scale_factor, swvdr, swvdf, swidr, swidf, & strocnxT_iavg, strocnyT_iavg, sst, frzmlt, & stressp_1, stressp_2, stressp_3, stressp_4, & @@ -63,7 +64,7 @@ subroutine dumpfile(filename_spec) stresspT, stressmT, stress12T, & stresspU, stressmU, stress12U use ice_flux, only: coszen - use ice_grid, only: grid_ice, tmask, iceumask, iceemask, icenmask + use ice_grid, only: grid_ice, tmask use ice_state, only: aicen, vicen, vsnon, trcrn, uvel, vvel, & uvelE, vvelE, uvelN, vvelN @@ -220,7 +221,7 @@ subroutine dumpfile(filename_spec) do j = 1, ny_block do i = 1, nx_block work1(i,j,iblk) = c0 - if (iceumask(i,j,iblk)) work1(i,j,iblk) = c1 + if (iceUmask(i,j,iblk)) work1(i,j,iblk) = c1 enddo enddo enddo @@ -234,7 +235,7 @@ subroutine dumpfile(filename_spec) do j = 1, ny_block do i = 1, nx_block work1(i,j,iblk) = c0 - if (icenmask(i,j,iblk)) work1(i,j,iblk) = c1 + if (iceNmask(i,j,iblk)) work1(i,j,iblk) = c1 enddo enddo enddo @@ -246,7 +247,7 @@ subroutine dumpfile(filename_spec) do j = 1, ny_block do i = 1, nx_block work1(i,j,iblk) = c0 - if (iceemask(i,j,iblk)) work1(i,j,iblk) = c1 + if (iceEmask(i,j,iblk)) work1(i,j,iblk) = c1 enddo enddo enddo @@ -276,6 +277,7 @@ subroutine restartfile (ice_ic) use ice_domain, only: nblocks, halo_info use ice_domain_size, only: nilyr, nslyr, ncat, & max_blocks + use ice_dyn_shared, only: iceUmask, iceEmask, iceNmask use ice_flux, only: scale_factor, swvdr, swvdf, swidr, swidf, & strocnxT_iavg, strocnyT_iavg, sst, frzmlt, & stressp_1, stressp_2, stressp_3, stressp_4, & @@ -284,8 +286,7 @@ subroutine restartfile (ice_ic) stresspT, stressmT, stress12T, & stresspU, stressmU, stress12U use ice_flux, only: coszen - use ice_grid, only: tmask, grid_type, grid_ice, & - iceumask, iceemask, icenmask, grid_average_X2Y + use ice_grid, only: tmask, grid_type, grid_ice, grid_average_X2Y use ice_state, only: trcr_depend, aice, vice, vsno, trcr, & aice0, aicen, vicen, vsnon, trcrn, aice_init, uvel, vvel, & uvelE, vvelE, uvelN, vvelN, uvelT, vvelT, & @@ -532,12 +533,12 @@ subroutine restartfile (ice_ic) call read_restart_field(nu_restart,0,work1,'ruf8', & 'iceumask',1,diag,field_loc_center, field_type_scalar) - iceumask(:,:,:) = .false. + iceUmask(:,:,:) = .false. !$OMP PARALLEL DO PRIVATE(iblk,i,j) do iblk = 1, nblocks do j = 1, ny_block do i = 1, nx_block - if (work1(i,j,iblk) > p5) iceumask(i,j,iblk) = .true. + if (work1(i,j,iblk) > p5) iceUmask(i,j,iblk) = .true. enddo enddo enddo @@ -549,12 +550,12 @@ subroutine restartfile (ice_ic) call read_restart_field(nu_restart,0,work1,'ruf8', & 'icenmask',1,diag,field_loc_center, field_type_scalar) - icenmask(:,:,:) = .false. + iceNmask(:,:,:) = .false. !$OMP PARALLEL DO PRIVATE(iblk,i,j) do iblk = 1, nblocks do j = 1, ny_block do i = 1, nx_block - if (work1(i,j,iblk) > p5) icenmask(i,j,iblk) = .true. + if (work1(i,j,iblk) > p5) iceNmask(i,j,iblk) = .true. enddo enddo enddo @@ -565,12 +566,12 @@ subroutine restartfile (ice_ic) call read_restart_field(nu_restart,0,work1,'ruf8', & 'iceemask',1,diag,field_loc_center, field_type_scalar) - iceemask(:,:,:) = .false. + iceEmask(:,:,:) = .false. !$OMP PARALLEL DO PRIVATE(iblk,i,j) do iblk = 1, nblocks do j = 1, ny_block do i = 1, nx_block - if (work1(i,j,iblk) > p5) iceemask(i,j,iblk) = .true. + if (work1(i,j,iblk) > p5) iceEmask(i,j,iblk) = .true. enddo enddo enddo @@ -710,13 +711,14 @@ subroutine restartfile_v4 (ice_ic) use ice_domain, only: nblocks, distrb_info use ice_domain_size, only: nilyr, nslyr, ncat, nx_global, ny_global, & max_blocks + use ice_dyn_shared, only: iceUmask use ice_flux, only: scale_factor, swvdr, swvdf, swidr, swidf, & strocnxT_iavg, strocnyT_iavg, sst, frzmlt, & stressp_1, stressp_2, stressp_3, stressp_4, & stressm_1, stressm_2, stressm_3, stressm_4, & stress12_1, stress12_2, stress12_3, stress12_4 use ice_gather_scatter, only: scatter_global_stress - use ice_grid, only: tmask, iceumask + use ice_grid, only: tmask use ice_read_write, only: ice_open, ice_read, ice_read_global use ice_state, only: trcr_depend, aice, vice, vsno, trcr, & aice0, aicen, vicen, vsnon, trcrn, aice_init, uvel, vvel, & @@ -945,12 +947,12 @@ subroutine restartfile_v4 (ice_ic) call ice_read(nu_restart,0,work1,'ruf8',diag, & field_loc_center, field_type_scalar) - iceumask(:,:,:) = .false. + iceUmask(:,:,:) = .false. !$OMP PARALLEL DO PRIVATE(iblk,i,j) do iblk = 1, nblocks do j = 1, ny_block do i = 1, nx_block - if (work1(i,j,iblk) > p5) iceumask(i,j,iblk) = .true. + if (work1(i,j,iblk) > p5) iceUmask(i,j,iblk) = .true. enddo enddo enddo diff --git a/doc/source/cice_index.rst b/doc/source/cice_index.rst index 99679e791..9dc2a06c4 100644 --- a/doc/source/cice_index.rst +++ b/doc/source/cice_index.rst @@ -354,10 +354,10 @@ either Celsius or Kelvin units). "icells", "number of grid cells with specified property (for vectorization)", "" "iceruf", "ice surface roughness at atmosphere interface", "5.\ :math:`\times`\ 10\ :math:`^{-4}` m" "iceruf_ocn", "under-ice roughness (at ocean interface)", "0.03 m" - "iceemask", "ice extent mask (E-cell)", "" - "icenmask", "ice extent mask (N-cell)", "" - "icetmask", "ice extent mask (T-cell)", "" - "iceumask", "ice extent mask (U-cell)", "" + "iceEmask", "dynamics ice extent mask (E-cell)", "" + "iceNmask", "dynamics ice extent mask (N-cell)", "" + "iceTmask", "dynamics ice extent mask (T-cell)", "" + "iceUmask", "dynamics ice extent mask (U-cell)", "" "idate", "the date at the end of the current time step (yyyymmdd)", "" "idate0", "initial date", "" "ierr", "general-use error flag", "" diff --git a/doc/source/user_guide/ug_implementation.rst b/doc/source/user_guide/ug_implementation.rst index 4bcfe1ede..a7cc66948 100644 --- a/doc/source/user_guide/ug_implementation.rst +++ b/doc/source/user_guide/ug_implementation.rst @@ -426,16 +426,16 @@ respectively) are useful in conditional statements. In addition to the land masks, two other masks are implemented in *dyn\_prep* in order to reduce the dynamics component’s work on a global -grid. At each time step the logical masks ``icetmask`` and ``iceumask`` are +grid. At each time step the logical masks ``iceTmask`` and ``iceUmask`` are determined from the current ice extent, such that they have the value “true” wherever ice exists. They also include a border of cells around the ice pack for numerical purposes. These masks are used in the dynamics component to prevent unnecessary calculations on grid points where there is no ice. They are not used in the thermodynamics component, so that ice may form in previously ice-free cells. Like the -land masks ``hm`` and ``uvm``, the ice extent masks ``icetmask`` and ``iceumask`` +land masks ``hm`` and ``uvm``, the ice extent masks ``iceTmask`` and ``iceUmask`` are for T-cells and U-cells, respectively. Note that the ice extent masks -``iceemask`` and ``icenmask`` are also defined when using the C or CD grid. +``iceEmask`` and ``iceNmask`` are also defined when using the C or CD grid. Improved parallel performance may result from utilizing halo masks for boundary updates of the full ice state, incremental remapping transport,