First round of housekeeping on ice_grid (CICE-Consortium#921)

* removal of unused variables.

* moved xav to transport. Could remove commented code. Could remove xav and yav as they are zero

* Move derived parameters and only allocate if needed

* bugfixes for cxp, cyp...

* fix index and remove commented code in ice_grid

* new version of transport_remap. xav, yav array where needed. xxav, yyav parameter

* Removed comments rom ice_transport_remap and arrays for nonuniform grids

cicecore/cicedyn/dynamics/ice_dyn_eap.F90

@@ -109,7 +109,7 @@ subroutine eap (dt)
           seabed_stress_factor_LKD, seabed_stress_factor_prob, &
           seabed_stress_method, seabed_stress, &
           stack_fields, unstack_fields, iceTmask, iceUmask, &
-          fld2, fld3, fld4
+          fld2, fld3, fld4, dxhy, dyhx, cxp, cyp, cxm, cym
       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, &
@@ -118,7 +118,7 @@ subroutine eap (dt)
           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_grid, only: tmask, umask, dxT, dyT, dxhy, dyhx, cxp, cyp, cxm, cym, &
+      use ice_grid, only: tmask, umask, dxT, dyT, &
           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, &

cicecore/cicedyn/dynamics/ice_dyn_evp.F90

@@ -84,6 +84,12 @@ module ice_dyn_evp
          emass    (:,:,:) , & ! total mass of ice and snow (E grid)
          emassdti (:,:,:)     ! mass of E-cell/dte (kg/m^2 s)
 
+      real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: &
+         ratiodxN    , & ! - dxN(i+1,j)   / dxN(i,j)
+         ratiodyE    , & ! - dyE(i  ,j+1) / dyE(i,j)
+         ratiodxNr   , & !   1 / ratiodxN
+         ratiodyEr       !   1 / ratiodyE
+
       real (kind=dbl_kind), allocatable :: &
          strengthU(:,:,:) , & ! strength averaged to U points
          divergU  (:,:,:) , & ! div array on U points, differentiate from divu
@@ -119,9 +125,10 @@ module ice_dyn_evp
 ! Elastic-viscous-plastic dynamics driver
 !
       subroutine init_evp
-      use ice_blocks, only: nx_block, ny_block, nghost
-      use ice_domain_size, only: max_blocks, nx_global, ny_global
-      use ice_grid, only: grid_ice, dyT, dxT, uarear, tmask, G_HTE, G_HTN
+      use ice_blocks, only: get_block, nx_block, ny_block, nghost, block
+      use ice_domain_size, only: max_blocks
+      use ice_domain, only: nblocks, blocks_ice
+      use ice_grid, only: grid_ice, dyT, dxT, uarear, tmask, G_HTE, G_HTN, dxN, dyE
       use ice_calendar, only: dt_dyn
       use ice_dyn_shared, only: init_dyn_shared, evp_algorithm
       use ice_dyn_evp1d, only: dyn_evp1d_init
@@ -131,6 +138,14 @@ subroutine init_evp
 
       character(len=*), parameter :: subname = '(init_evp)'
 
+     type (block) :: &
+         this_block   ! block information for current block
+
+      integer (kind=int_kind) :: &
+         i, j, iblk            , & ! block index
+         ilo,ihi,jlo,jhi           ! beginning and end of physical domain
+
+
       call init_dyn_shared(dt_dyn)
 
       if (evp_algorithm == "shared_mem_1d" ) then
@@ -197,6 +212,32 @@ subroutine init_evp
                    stat=ierr)
          if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory E evp')
 
+         allocate( ratiodxN (nx_block,ny_block,max_blocks), &
+                   ratiodyE (nx_block,ny_block,max_blocks), &
+                   ratiodxNr(nx_block,ny_block,max_blocks), &
+                   ratiodyEr(nx_block,ny_block,max_blocks), &
+                   stat=ierr)
+         if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory ratio')
+
+         !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
+         do iblk = 1, nblocks
+            this_block = get_block(blocks_ice(iblk),iblk)
+            ilo = this_block%ilo
+            ihi = this_block%ihi
+            jlo = this_block%jlo
+            jhi = this_block%jhi
+
+            do j = jlo, jhi
+            do i = ilo, ihi
+               ratiodxN (i,j,iblk) = - dxN(i+1,j  ,iblk) / dxN(i,j,iblk)
+               ratiodyE (i,j,iblk) = - dyE(i  ,j+1,iblk) / dyE(i,j,iblk)
+               ratiodxNr(i,j,iblk) =   c1 / ratiodxN(i,j,iblk)
+               ratiodyEr(i,j,iblk) =   c1 / ratiodyE(i,j,iblk)
+            enddo
+            enddo
+         enddo                     ! iblk
+         !$OMP END PARALLEL DO
+
       endif
 
       end subroutine init_evp
@@ -219,7 +260,7 @@ subroutine evp (dt)
           ice_HaloDestroy, ice_HaloUpdate_stress
       use ice_blocks, only: block, get_block, nx_block, ny_block, nghost
       use ice_domain, only: nblocks, blocks_ice, halo_info, maskhalo_dyn
-      use ice_domain_size, only: max_blocks, ncat, nx_global, ny_global
+      use ice_domain_size, only: max_blocks, ncat
       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, &
@@ -238,8 +279,6 @@ subroutine evp (dt)
           stresspU, stressmU, stress12U
       use ice_grid, only: tmask, umask, umaskCD, nmask, emask, uvm, epm, npm, &
           dxE, dxN, dxT, dxU, dyE, dyN, dyT, dyU, &
-          ratiodxN, ratiodxNr, ratiodyE, ratiodyEr, &
-          dxhy, dyhx, cxp, cyp, cxm, cym, &
           tarear, uarear, earear, narear, grid_average_X2Y, uarea, &
           grid_type, grid_ice, &
           grid_atm_dynu, grid_atm_dynv, grid_ocn_dynu, grid_ocn_dynv
@@ -250,7 +289,7 @@ subroutine evp (dt)
           ice_timer_start, ice_timer_stop, timer_evp
       use ice_dyn_shared, only: evp_algorithm, stack_fields, unstack_fields, &
           DminTarea, visc_method, deformations, deformationsC_T, deformationsCD_T, &
-          strain_rates_U, &
+          strain_rates_U, dxhy, dyhx, cxp, cyp, cxm, cym, &
           iceTmask, iceUmask, iceEmask, iceNmask, &
           dyn_haloUpdate, fld2, fld3, fld4
       use ice_dyn_evp1d, only: dyn_evp1d_run

cicecore/cicedyn/dynamics/ice_dyn_shared.F90

@@ -11,7 +11,7 @@ module ice_dyn_shared
 
       use ice_kinds_mod
       use ice_communicate, only: my_task, master_task, get_num_procs
-      use ice_constants, only: c0, c1, c2, c3, c4, c6
+      use ice_constants, only: c0, c1, c2, c3, c4, c6, c1p5
       use ice_constants, only: omega, spval_dbl, p01, p001, p5
       use ice_blocks, only: nx_block, ny_block
       use ice_domain_size, only: max_blocks
@@ -119,6 +119,14 @@ module ice_dyn_shared
       real (kind=dbl_kind), allocatable, public :: &
          DminTarea(:,:,:)    ! deltamin * tarea (m^2/s)
 
+      real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: &
+         cyp    , & ! 1.5*HTE(i,j)-0.5*HTW(i,j) = 1.5*HTE(i,j)-0.5*HTE(i-1,j)
+         cxp    , & ! 1.5*HTN(i,j)-0.5*HTS(i,j) = 1.5*HTN(i,j)-0.5*HTN(i,j-1)
+         cym    , & ! 0.5*HTE(i,j)-1.5*HTW(i,j) = 0.5*HTE(i,j)-1.5*HTE(i-1,j)
+         cxm    , & ! 0.5*HTN(i,j)-1.5*HTS(i,j) = 0.5*HTN(i,j)-1.5*HTN(i,j-1)
+         dxhy   , & ! 0.5*(HTE(i,j) - HTW(i,j)) = 0.5*(HTE(i,j) - HTE(i-1,j))
+         dyhx       ! 0.5*(HTN(i,j) - HTS(i,j)) = 0.5*(HTN(i,j) - HTN(i,j-1))
+
       ! ice isotropic tensile strength parameter
       real (kind=dbl_kind), public :: &
          Ktens               ! T=Ktens*P (tensile strength: see Konig and Holland, 2010)
@@ -192,6 +200,22 @@ subroutine alloc_dyn_shared
          stat=ierr)
       if (ierr/=0) call abort_ice(subname//': Out of memory')
 
+      allocate( &
+         cyp(nx_block,ny_block,max_blocks), & ! 1.5*HTE - 0.5*HTW
+         cxp(nx_block,ny_block,max_blocks), & ! 1.5*HTN - 0.5*HTS
+         cym(nx_block,ny_block,max_blocks), & ! 0.5*HTE - 1.5*HTW
+         cxm(nx_block,ny_block,max_blocks), & ! 0.5*HTN - 1.5*HTS
+         stat=ierr)
+      if (ierr/=0) call abort_ice(subname//': Out of memory')
+
+      if (grid_ice == 'B' .and. evp_algorithm == "standard_2d") then
+         allocate( &
+            dxhy(nx_block,ny_block,max_blocks), & ! 0.5*(HTE - HTW)
+            dyhx(nx_block,ny_block,max_blocks), & ! 0.5*(HTN - HTS)
+            stat=ierr)
+         if (ierr/=0) call abort_ice(subname//': Out of memory')
+      endif
+
       if (grid_ice == 'CD' .or. grid_ice == 'C') then
          allocate( &
             uvelE_init (nx_block,ny_block,max_blocks), & ! x-component of velocity (m/s), beginning of timestep
@@ -214,8 +238,9 @@ end subroutine alloc_dyn_shared
 
       subroutine init_dyn_shared (dt)
 
-      use ice_blocks, only: nx_block, ny_block
-      use ice_domain, only: nblocks, halo_dynbundle
+      use ice_blocks, only: block, get_block
+      use ice_boundary, only: ice_halo, ice_haloUpdate
+      use ice_domain, only: nblocks, halo_dynbundle, blocks_ice, halo_info
       use ice_domain_size, only: max_blocks
       use ice_flux, only: &
           stressp_1, stressp_2, stressp_3, stressp_4, &
@@ -224,17 +249,22 @@ subroutine init_dyn_shared (dt)
           stresspT, stressmT, stress12T, &
           stresspU, stressmU, stress12U
       use ice_state, only: uvel, vvel, uvelE, vvelE, uvelN, vvelN
-      use ice_grid, only: ULAT, NLAT, ELAT, tarea
+      use ice_grid, only: ULAT, NLAT, ELAT, tarea, HTE, HTN
+      use ice_constants, only: field_loc_center, field_type_vector
 
       real (kind=dbl_kind), intent(in) :: &
          dt      ! time step
 
       ! local variables
 
       integer (kind=int_kind) :: &
-         i, j  , &  ! indices
-         nprocs, &  ! number of processors
-         iblk       ! block index
+         i, j  ,             & ! indices
+         ilo, ihi, jlo, jhi, & !min and max index for interior of blocks
+         nprocs,             & ! number of processors
+         iblk                  ! block index
+
+      type (block) :: &
+         this_block   ! block information for current block
 
       character(len=*), parameter :: subname = '(init_dyn_shared)'
 
@@ -333,6 +363,48 @@ subroutine init_dyn_shared (dt)
       enddo                     ! iblk
       !$OMP END PARALLEL DO
 
+      if (grid_ice == 'B' .and. evp_algorithm == "standard_2d") then
+         do iblk = 1, nblocks
+            this_block = get_block(blocks_ice(iblk),iblk)
+            ilo = this_block%ilo
+            ihi = this_block%ihi
+            jlo = this_block%jlo
+            jhi = this_block%jhi
+
+            do j = jlo, jhi
+            do i = ilo, ihi
+               dxhy(i,j,iblk) = p5*(HTE(i,j,iblk) - HTE(i-1,j,iblk))
+               dyhx(i,j,iblk) = p5*(HTN(i,j,iblk) - HTN(i,j-1,iblk))
+            enddo
+            enddo
+         enddo
+
+         call ice_HaloUpdate (dxhy,               halo_info, &
+                              field_loc_center,   field_type_vector, &
+                              fillValue=c1)
+         call ice_HaloUpdate (dyhx,               halo_info, &
+                              field_loc_center,   field_type_vector, &
+                              fillValue=c1)
+
+     endif
+
+     do iblk = 1, nblocks
+        this_block = get_block(blocks_ice(iblk),iblk)
+        ilo = this_block%ilo
+        ihi = this_block%ihi
+        jlo = this_block%jlo
+        jhi = this_block%jhi
+        do j = jlo, jhi+1
+        do i = ilo, ihi+1
+           cyp(i,j,iblk) = (c1p5*HTE(i,j,iblk) - p5*HTE(i-1,j,iblk))
+           cxp(i,j,iblk) = (c1p5*HTN(i,j,iblk) - p5*HTN(i,j-1,iblk))
+            ! match order of operations in cyp, cxp for tripole grids
+           cym(i,j,iblk) = -(c1p5*HTE(i-1,j,iblk) - p5*HTE(i,j,iblk))
+           cxm(i,j,iblk) = -(c1p5*HTN(i,j-1,iblk) - p5*HTN(i,j,iblk))
+        enddo
+        enddo
+     enddo
+
   end subroutine init_dyn_shared
 
 !=======================================================================

cicecore/cicedyn/dynamics/ice_dyn_vp.F90

@@ -52,7 +52,7 @@ module ice_dyn_vp
       use ice_fileunits, only: nu_diag
       use ice_flux, only: fmU
       use ice_global_reductions, only: global_sum
-      use ice_grid, only: dxT, dyT, dxhy, dyhx, cxp, cyp, cxm, cym, uarear
+      use ice_grid, only: dxT, dyT, uarear
       use ice_exit, only: abort_ice
       use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
       use icepack_intfc, only: icepack_ice_strength, icepack_query_parameters
@@ -120,19 +120,9 @@ subroutine init_vp
       use ice_boundary, only: ice_HaloUpdate
       use ice_constants, only: c1, &
           field_loc_center, field_type_scalar
-      use ice_domain, only: blocks_ice, halo_info
+      use ice_domain, only: blocks_ice
       use ice_calendar, only: dt_dyn
       use ice_dyn_shared, only: init_dyn_shared
-!      use ice_grid, only: tarea
-
-      ! local variables
-
-      integer (kind=int_kind) :: &
-         i, j, iblk, &
-         ilo,ihi,jlo,jhi      ! beginning and end of physical domain
-
-      type (block) :: &
-         this_block           ! block information for current block
 
       call init_dyn_shared(dt_dyn)
 
@@ -167,7 +157,8 @@ 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, iceTmask, iceUmask
+      use ice_dyn_shared, only: deformations, iceTmask, iceUmask, &
+          cxp, cyp, cxm, cym
       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, &
@@ -176,7 +167,7 @@ subroutine implicit_solver (dt)
           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_grid, only: tmask, umask, dxT, dyT, cxp, cyp, cxm, cym, &
+      use ice_grid, only: tmask, umask, dxT, dyT, &
           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, &
@@ -686,9 +677,8 @@ subroutine anderson_solver (icellT  , icellU , &
       use ice_domain, only: maskhalo_dyn, halo_info
       use ice_domain_size, only: max_blocks
       use ice_flux, only:   fmU, TbU
-      use ice_grid, only: dxT, dyT, dxhy, dyhx, cxp, cyp, cxm, cym, &
-           uarear
-      use ice_dyn_shared, only: DminTarea
+      use ice_grid, only: dxT, dyT, uarear
+      use ice_dyn_shared, only: DminTarea, dxhy, dyhx, cxp, cyp, cxm, cym
       use ice_state, only: uvel, vvel, strength
       use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bound
 
@@ -2753,6 +2743,7 @@ subroutine fgmres (zetax2   , etax2   , &
       use ice_boundary, only: ice_HaloUpdate
       use ice_domain, only: maskhalo_dyn, halo_info
       use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bound
+      use ice_dyn_shared, only: dxhy, dyhx, cxp, cyp, cxm, cym
 
       real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks,4), intent(in) :: &
          zetax2   , & ! zetax2 = 2*zeta (bulk viscosity)
@@ -3154,6 +3145,7 @@ subroutine pgmres (zetax2   , etax2   , &
       use ice_boundary, only: ice_HaloUpdate
       use ice_domain, only: maskhalo_dyn, halo_info
       use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bound
+      use ice_dyn_shared, only: dyhx, dxhy, cxp, cyp, cxm, cym
 
       real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks,4), intent(in) :: &
          zetax2   , & ! zetax2 = 2*zeta (bulk viscosity)