From e3112ebb854d2e7f6b660f5823af9705dd9da602 Mon Sep 17 00:00:00 2001
From: "David A. Bailey" <dbailey@ucar.edu>
Date: Thu, 18 Nov 2021 13:29:36 -0700
Subject: [PATCH] Add some variables for CD discretization (#21)

* Prototype variables for CD grid discretization

* Fix syntax error

* Change variable names to ratio and only need 4

* Needed minus signs and fix alignment

* Make the ratios consistent with the doc
---
 cicecore/cicedynB/infrastructure/ice_grid.F90 | 27 +++++++++++++++++++
 1 file changed, 27 insertions(+)

diff --git a/cicecore/cicedynB/infrastructure/ice_grid.F90 b/cicecore/cicedynB/infrastructure/ice_grid.F90
index 6e554b001..6a3cf1171 100644
--- a/cicecore/cicedynB/infrastructure/ice_grid.F90
+++ b/cicecore/cicedynB/infrastructure/ice_grid.F90
@@ -107,6 +107,12 @@ module ice_grid
          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)) 
 
+      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
+
       ! grid dimensions for rectangular grid
       real (kind=dbl_kind), public ::  &
          dxrect, & !  user_specified spacing (cm) in x-direction (uniform HTN)
@@ -255,6 +261,16 @@ subroutine alloc_grid
          stat=ierr)
       if (ierr/=0) call abort_ice(subname//'ERROR: Out of memory')
 
+      if (grid_system == 'CD') then
+         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')
+      endif
+
       if (pgl_global_ext) then
          allocate( &
             G_HTE(nx_global+2*nghost, ny_global+2*nghost), & ! length of eastern edge of T-cell (global ext.)
@@ -502,6 +518,17 @@ subroutine init_grid2
          enddo
          enddo
 
+         if (grid_system == 'CD') then
+            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
+         endif
+
       enddo                     ! iblk
       !$OMP END PARALLEL DO