Added dxgrow, dygrow to facilitate variable spaced grid. Modified rec…

cicecore/cicedynB/general/ice_init.F90

@@ -110,7 +110,7 @@ subroutine input_data
                           grid_ice, grid_ice_thrm, grid_ice_dynu, grid_ice_dynv, &
                           grid_ocn, grid_ocn_thrm, grid_ocn_dynu, grid_ocn_dynv, & 
                           grid_atm, grid_atm_thrm, grid_atm_dynu, grid_atm_dynv, &
-                          dxrect, dyrect, &
+                          dxrect, dyrect, dxscale, dyscale, scale_dxdy, &
                           pgl_global_ext
       use ice_dyn_shared, only: ndte, kdyn, revised_evp, yield_curve, &
                                 evp_algorithm, visc_method,     &
@@ -208,6 +208,7 @@ subroutine input_data
         bathymetry_file, use_bathymetry, nfsd,          bathymetry_format, &
         ncat,           nilyr,           nslyr,         nblyr,          &
         kcatbound,      gridcpl_file,    dxrect,        dyrect,         &
+        dxscale,        dyscale,         scale_dxdy,                    &
         close_boundaries, orca_halogrid, grid_ice,      kmt_type,       &
         grid_atm,       grid_ocn
 
@@ -394,6 +395,9 @@ subroutine input_data
       ksno = 0.3_dbl_kind     ! snow thermal conductivity
       dxrect = 0.0_dbl_kind   ! user defined grid spacing in cm in x direction
       dyrect = 0.0_dbl_kind   ! user defined grid spacing in cm in y direction
+      scale_dxdy = .false.    ! apply dxscale, dyscale to rectgrid
+      dxscale = 1.0_dbl_kind   ! user defined rectgrid x-grid scale factor (e.g., 1.02)
+      dyscale = 1.0_dbl_kind   ! user defined rectgrid y-grid scale factor (e.g., 1.02)
       close_boundaries = .false.   ! true = set land on edges of grid
       seabed_stress= .false.  ! if true, seabed stress for landfast is on
       seabed_stress_method  = 'LKD'! LKD = Lemieux et al 2015, probabilistic = Dupont et al. in prep
@@ -847,6 +851,9 @@ subroutine input_data
       call broadcast_scalar(grid_format,          master_task)
       call broadcast_scalar(dxrect,               master_task)
       call broadcast_scalar(dyrect,               master_task)
+      call broadcast_scalar(scale_dxdy,           master_task)
+      call broadcast_scalar(dxscale,              master_task)
+      call broadcast_scalar(dyscale,              master_task)
       call broadcast_scalar(close_boundaries,     master_task)
       call broadcast_scalar(grid_type,            master_task)
       call broadcast_scalar(grid_ice,             master_task)

cicecore/cicedynB/infrastructure/ice_grid.F90

@@ -129,6 +129,12 @@ module ice_grid
          dxrect, & !  user_specified spacing (cm) in x-direction (uniform HTN)
          dyrect    !  user_specified spacing (cm) in y-direction (uniform HTE)
 
+      ! growth factor for variable spaced grid
+      real (kind=dbl_kind), public ::  &
+         dxscale, & !  scale factor for grid spacing in x direction (e.g., 1.02)
+         dyscale    !  scale factor for gird spacing in y direction (e.g., 1.02)
+
+
       ! Corners of grid boxes for history output
       real (kind=dbl_kind), dimension (:,:,:,:), allocatable, public :: &
          lont_bounds, & ! longitude of gridbox corners for T point
@@ -171,7 +177,8 @@ module ice_grid
 
       logical (kind=log_kind), public :: &
          use_bathymetry, & ! flag for reading in bathymetry_file
-         pgl_global_ext    ! flag for init primary grid lengths (global ext.)
+         pgl_global_ext, & ! flag for init primary grid lengths (global ext.)
+         scale_dxdy        ! flag to apply scale factor to vary dx/dy in rectgrid
 
       logical (kind=log_kind), dimension (:,:,:), allocatable, public :: &
          tmask  , & ! land/boundary mask, thickness (T-cell)
@@ -1360,8 +1367,12 @@ subroutine rectgrid
          i, j, iblk, &
          imid, jmid
 
+      integer (kind=int_kind) :: &
+           center1, center2 ! array centers for expanding dx, dy
+
       real (kind=dbl_kind) :: &
-         length, &
+         length,  &
+         latref, lonref, & ! reference lat, lon
          rad_to_deg
 
       real (kind=dbl_kind), dimension(:,:), allocatable :: &
@@ -1386,65 +1397,150 @@ subroutine rectgrid
 
       ! Weddell Sea
       ! lower left corner of grid is 55W, 75S
+      !lonref = -55.0_dbl_kind
+      !latref = -75.0_dbl_kind
 
       ! Barrow AK
       ! lower left corner of grid is 156.5W, 71.35N
-
+      lonref = -156.5_dbl_kind
+      latref = 71.35_dbl_kind
+
+      ! determine dx spacing
+      ! strategy: initialize with dxrect.
+      ! if want to scale the grid, work from center outwards,
+      ! multplying neighbor cell by scale factor.
+      ! this assumes dx varies in x direction only. 
+      !(i.e, dx is the same across same y location)
       if (my_task == master_task) then
-         work_g1 = c0
-         length = dxrect*cm_to_m/radius*rad_to_deg
 
-!         work_g1(1,:) = -55._dbl_kind   ! Weddell Sea
-         work_g1(1,:) = -156.5_dbl_kind ! Barrow AK
+         ! initialize with initial dxrect
+         work_g1(:,:) = dxrect
+
+         ! scale if desired
+         if (scale_dxdy) then
+
+            ! check if nx is even or odd
+            ! if even, middle 2 columns are center
+            ! of odd,  middle 1 column is center
+            if (mod(nx_global,2) == 0) then ! nx_global is even
+
+               ! with even number of x locatons, 
+               ! the center two y columns are center
+               center1 = nx_global/2  ! integer math
+               center2 = center1 + 1  ! integer math
+
+            else ! nx_global = odd
+               ! only one center index. set center2=center1
+               center1 = ceiling(real(nx_global/2),int_kind) 
+               center2 = center1
+            endif
+
+            ! note loop over only half the x grid points (center1)-1 
+            ! working from the center outward.
+            do j = 1, ny_global
+            do i = 1, center1-1 
+               ! work from center1 to left
+               work_g1(center1-i,j) = dxscale*work_g1(center1-i+1,j)
+
+               ! work from center2 to right
+               work_g1(center2+i,j) = dxscale*work_g1(center2+i-1,j)
+            enddo ! i
+            enddo ! j 
+         endif    ! scale_dxdy
+      endif       ! my_task == master_task
+
+      ! note work_g1 is converted to meters in primary_grid_lengths_HTN
+      call primary_grid_lengths_HTN(work_g1)  ! dxU, dxT, dxN, dxE
+
+      ! make ULON array
+      if (my_task == master_task) then
+
+         ! make first column reference lon in radians.
+         ! the remaining work_g1 is still dx in meters
+         work_g1(1,:) = lonref/rad_to_deg ! radians
 
+         ! loop over remaining points and add spacing to successive 
+         ! x locations
          do j = 1, ny_global
-         do i = 2, nx_global
+         do i = 2, nx_global ! start from i=2. i=1 is lonref
+            length = work_g1(i,j)/radius             ! grid spacing in radians
             work_g1(i,j) = work_g1(i-1,j) + length   ! ULON
-         enddo
-         enddo
-         work_g1(:,:) = work_g1(:,:) / rad_to_deg
-      endif
+         enddo ! i
+         enddo ! j
+      endif    ! mytask == master_task
       call scatter_global(ULON, work_g1, master_task, distrb_info, &
                           field_loc_NEcorner, field_type_scalar)
       call ice_HaloExtrapolate(ULON, distrb_info, &
                                ew_boundary_type, ns_boundary_type)
 
+      ! determine dy spacing
+      ! strategy: initialize with dyrect.
+      ! if want to scale the grid, work from center outwards,
+      ! multplying neighbor cell by scale factor.
+      ! this assumes dy varies in y direction only. 
+      ! (i.e, dy is the same across same x location)
       if (my_task == master_task) then
-         work_g1 = c0
-         length = dyrect*cm_to_m/radius*rad_to_deg
 
-!         work_g1(:,1) = -75._dbl_kind ! Weddell Sea
-         work_g1(:,1) = 71.35_dbl_kind ! Barrow AK
+         ! initialize with initial dxrect
+         work_g1(:,:) = dyrect
+
+         ! scale if desired
+         if (scale_dxdy) then
+
+            ! check if ny is even or odd
+            ! if even, middle 2 rows are center
+            ! of odd,  middle 1 row is center
+            if (mod(ny_global,2) == 0) then ! ny_global is even
+
+               ! with even number of x locatons, 
+               ! the center two y columns are center
+               center1 = ny_global/2  ! integer math
+               center2 = center1 + 1  ! integer math
+
+            else ! ny_global = odd
+               ! only one center index. set center2=center1
+               center1 = ceiling(real(ny_global/2),int_kind) 
+               center2 = center1
+            endif
+
+            ! note loop over only half the y grid points (center1)-1 
+            ! working from the center outward.
+            do i = 1, nx_global
+            do j = 1, center1-1 
+               ! work from center1 to bottom
+               work_g1(i,center1-j) = dyscale*work_g1(i,center1-j+1)
+
+               ! work from center2 to top
+               work_g1(i,center2+j) = dyscale*work_g1(i,center2+j-1)
+            enddo ! i
+            enddo ! j 
+         endif    ! scale_dxdy
+      endif       ! my_task == master_task
+
+      ! note work_g1 is converted to meters primary_grid_lengths_HTE
+      call primary_grid_lengths_HTE(work_g1)  ! dyU, dyT, dyN, dyE
+
+      ! make ULAT array
+      if (my_task == master_task) then
+
+         ! make first row reference lat in radians.
+         ! the remaining work_g1 is still dy in meters
+         work_g1(:,1) = latref/rad_to_deg ! radians
 
+         ! loop over remaining points and add spacing to successive 
+         ! x locations
+         do j = 2, ny_global ! start from j=2. j=1 is latref
          do i = 1, nx_global
-         do j = 2, ny_global
+            length = work_g1(i,j)/radius             ! grid spacing in radians
             work_g1(i,j) = work_g1(i,j-1) + length   ! ULAT
-         enddo
-         enddo
-         work_g1(:,:) = work_g1(:,:) / rad_to_deg
-      endif
+         enddo ! i
+         enddo ! j
+      endif    ! mytask == master_task
       call scatter_global(ULAT, work_g1, master_task, distrb_info, &
                           field_loc_NEcorner, field_type_scalar)
       call ice_HaloExtrapolate(ULAT, distrb_info, &
                                ew_boundary_type, ns_boundary_type)
 
-      if (my_task == master_task) then
-         do j = 1, ny_global
-         do i = 1, nx_global
-            work_g1(i,j) = dxrect             ! HTN
-         enddo
-         enddo
-      endif
-      call primary_grid_lengths_HTN(work_g1)  ! dxU, dxT, dxN, dxE
-
-      if (my_task == master_task) then
-         do j = 1, ny_global
-         do i = 1, nx_global
-            work_g1(i,j) = dyrect             ! HTE
-         enddo
-         enddo
-      endif
-      call primary_grid_lengths_HTE(work_g1)  ! dyU, dyT, dyN, dyE
 
       !-----------------------------------------------------------------
       ! Construct T-cell land mask

configuration/scripts/ice_in

@@ -75,6 +75,9 @@
     kcatbound    = 0
     dxrect         = 30.e5
     dyrect         = 30.e5
+    scale_dxdy     = .false.
+    dxscale        = 1.0
+    dyscale        = 1.0
     close_boundaries = .false.
     ncat         = 5
     nfsd         = 1

configuration/scripts/options/set_nml.vargrid

@@ -0,0 +1,13 @@
+scale_dxdy = .true.
+dxscale    = 1.02
+dyscale    = 1.01
+f_dxt      = .true.
+f_dyt      = .true.
+f_dxu      = .true.
+f_dyu      = .true.
+f_dxe      = .true.
+f_dye      = .true.
+f_dxn      = .true.
+f_dyn      = .true.
+f_HTN      = .true.
+f_HTE      = .true.

configuration/scripts/tests/gridsys_suite.ts

@@ -10,9 +10,11 @@ smoke          gbox80  1x1        boxslotcyl
 smoke          gbox80  2x4        boxnodyn
 #smoke	       gbox80  2x2        boxsymn,run1day
 smoke	       gbox80  4x2        boxsyme,run1day
+smoke	       gbox80  4x2        boxsyme,run1day,vargrid
 #smoke	       gbox80  4x1        boxsymne,run1day
 #smoke	       gbox80  2x2        boxsymn,run1day,kmtislands
 smoke	       gbox80  4x1        boxsyme,run1day,kmtislands
+smoke	       gbox80  4x1        boxsyme,run1day,kmtislands,vargrid
 #smoke	       gbox80  4x2        boxsymne,run1day,kmtislands
 #smoke	       gbox80  8x1        boxislandsn,run1day
 smoke	       gbox80  4x2        boxislandse,run1day

configuration/scripts/tests/vargrid_suite.ts

@@ -0,0 +1,6 @@
+# Test         Grid    PEs        Sets    BFB-compare
+smoke	       gbox80  4x2        boxsyme,run1day,vargrid
+smoke	       gbox80  4x2        boxsyme,run1day
+smoke          gbox80  4x1        boxsyme,run1day,kmtislands,vargrid
+smoke          gbox80  4x1        boxsyme,run1day,kmtislands
+

doc/source/user_guide/ug_case_settings.rst

@@ -248,6 +248,9 @@ grid_nml
    "``close_boundaries``", "logical", "force two gridcell wide land mask on boundaries for rectangular grids", "``.false.``"
    "``dxrect``", "real", "x-direction grid spacing for rectangular grid in cm", "0.0"
    "``dyrect``", "real", "y-direction grid spacing for rectangular grid in cm", "0.0"
+   "``scale_dxdy``", "logical", "apply dxscale, dyscale to rectgrid", "``false``"
+   "``dxscale``", "real", "user defined rectgrid x-grid scale factor", "1.0"
+   "``dyscale``", "real", "user defined rectgrid y-grid scale factor", "1.0"
    "``gridcpl_file``", "string", "input file for coupling grid info", "'unknown_gridcpl_file'"
    "``grid_atm``", "``A``", "atm forcing/coupling grid, all fields on T grid", "``A``"
    "", "``B``", "atm forcing/coupling grid, thermo fields on T grid, dyn fields on U grid", ""