Advection test - slotted cylinder

cicecore/cicedynB/general/ice_forcing.F90

@@ -119,6 +119,7 @@ module ice_forcing
          bgc_data_type, & ! 'default', 'clim'
          ocn_data_type, & ! 'default', 'clim', 'ncar', 'oned',
                           ! 'hadgem_sst' or 'hadgem_sst_uvocn'
+         ice_data_type, & ! 'default', 'box2001', 'boxslotcyl'
          precip_units     ! 'mm_per_month', 'mm_per_sec', 'mks','m_per_sec'
 
       character(char_len_long), public :: & 

cicecore/cicedynB/general/ice_init.F90

@@ -84,7 +84,8 @@ subroutine input_data
           atm_data_format, ocn_data_format, &
           bgc_data_type, &
           ocn_data_type, ocn_data_dir,      &
-          oceanmixed_file, restore_ocn,   trestore
+          oceanmixed_file, restore_ocn,   trestore, & 
+          ice_data_type
       use ice_arrays_column, only: bgc_data_dir, fe_data_type
       use ice_grid, only: grid_file, gridcpl_file, kmt_file, &
                           bathymetry_file, use_bathymetry, &
@@ -191,6 +192,7 @@ subroutine input_data
         oceanmixed_ice, restore_ice,     restore_ocn,   trestore,       &
         precip_units,   default_season,                                 &
         atm_data_type,  ocn_data_type,   bgc_data_type, fe_data_type,   &
+        ice_data_type,                                                  &
         fyear_init,     ycycle,                                         &
         atm_data_dir,   ocn_data_dir,    bgc_data_dir,                  &
         atm_data_format, ocn_data_format,                               &
@@ -347,6 +349,7 @@ subroutine input_data
       ocn_data_format = 'bin'     ! file format ('bin'=binary or 'nc'=netcdf)
       bgc_data_type   = 'default'
       fe_data_type    = 'default'
+      ice_data_type   = 'default' ! used by some tests to initialize ice state (concentration, velocities)
       bgc_data_dir    = 'unknown_bgc_data_dir'
       ocn_data_type   = 'default'
       ocn_data_dir    = 'unknown_ocn_data_dir'
@@ -602,6 +605,7 @@ subroutine input_data
       call broadcast_scalar(ocn_data_format,    master_task)
       call broadcast_scalar(bgc_data_type,      master_task)
       call broadcast_scalar(fe_data_type,       master_task)
+      call broadcast_scalar(ice_data_type,      master_task)
       call broadcast_scalar(bgc_data_dir,       master_task)
       call broadcast_scalar(ocn_data_type,      master_task)
       call broadcast_scalar(ocn_data_dir,       master_task)
@@ -1091,6 +1095,8 @@ subroutine input_data
                                trim(bgc_data_type)
          write(nu_diag,*)    ' fe_data_type              = ', &
                                trim(fe_data_type)
+         write(nu_diag,*)    ' ice_data_type             = ', &
+                               trim(ice_data_type)
          write(nu_diag,*)    ' bgc_data_dir              = ', &
                                trim(bgc_data_dir)
          write(nu_diag,*)    ' ocn_data_type             = ', &
@@ -1220,7 +1226,7 @@ subroutine init_state
       use ice_grid, only: tmask, ULON, TLAT
       use ice_state, only: trcr_depend, aicen, trcrn, vicen, vsnon, &
           aice0, aice, vice, vsno, trcr, aice_init, bound_state, &
-          n_trcr_strata, nt_strata, trcr_base
+          n_trcr_strata, nt_strata, trcr_base, uvel, vvel
 
       integer (kind=int_kind) :: &
          ilo, ihi    , & ! physical domain indices
@@ -1410,7 +1416,8 @@ subroutine init_state
                              Tf   (:,:,    iblk),                      &
                              salinz(:,:,:, iblk), Tmltz(:,:,:,  iblk), &
                              aicen(:,:,  :,iblk), trcrn(:,:,:,:,iblk), &
-                             vicen(:,:,  :,iblk), vsnon(:,:,  :,iblk))
+                             vicen(:,:,  :,iblk), vsnon(:,:,  :,iblk), &
+                             uvel (:,:,    iblk), vvel (:,:,    iblk))
 
       enddo                     ! iblk
       !$OMP END PARALLEL DO
@@ -1487,12 +1494,13 @@ subroutine set_state_var (nx_block, ny_block, &
                                 Tf,       &
                                 salinz,   Tmltz, &
                                 aicen,    trcrn, &
-                                vicen,    vsnon)
+                                vicen,    vsnon, &
+                                uvel,     vvel)
 
       use ice_arrays_column, only: hin_max
       use ice_domain_size, only: nilyr, nslyr, nx_global, ny_global, ncat
       use ice_grid, only: grid_type
-      use ice_forcing, only: atm_data_type
+      use ice_forcing, only: ice_data_type
 
       integer (kind=int_kind), intent(in) :: &
          nx_block, ny_block, & ! block dimensions
@@ -1529,6 +1537,10 @@ subroutine set_state_var (nx_block, ny_block, &
          trcrn     ! ice tracers
                    ! 1: surface temperature of ice/snow (C)
 
+      real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out) :: &
+         uvel    , & ! ice velocity
+         vvel        ! 
+
       ! local variables
 
       integer (kind=int_kind) :: &
@@ -1619,18 +1631,30 @@ subroutine set_state_var (nx_block, ny_block, &
       !       extend to the prescribed edges.
       !-----------------------------------------------------------------
 
-         if (trim(atm_data_type) == 'box2001') then
+         if (trim(ice_data_type) == 'box2001') then
 
             hbar = c2  ! initial ice thickness
             do n = 1, ncat
                hinit(n) = c0
                ainit(n) = c0
                if (hbar > hin_max(n-1) .and. hbar < hin_max(n)) then
                   hinit(n) = hbar
-                  ainit(n) = 0.50 !echmod symm
+                  ainit(n) = p5 !echmod symm
                endif
             enddo
 
+         elseif (trim(ice_data_type) == 'boxslotcyl') then
+
+            hbar = c1  ! initial ice thickness (1 m)
+            do n = 1, ncat
+               hinit(n) = c0
+               ainit(n) = c0
+               if (hbar > hin_max(n-1) .and. hbar < hin_max(n)) then
+                  hinit(n) = hbar
+                  ainit(n) = c1 !echmod symm
+               endif
+            enddo
+
          else
 
       ! initial category areas in cells with ice
@@ -1653,7 +1677,7 @@ subroutine set_state_var (nx_block, ny_block, &
             ainit(n) = ainit(n) / (sum + puny/ncat) ! normalize
          enddo
 
-         endif ! atm_data_type
+         endif ! ice_data_type
 
          if (trim(grid_type) == 'rectangular') then
 
@@ -1704,7 +1728,7 @@ subroutine set_state_var (nx_block, ny_block, &
 
                aicen(i,j,n) = ainit(n)
 
-               if (trim(atm_data_type) == 'box2001') then
+               if (trim(ice_data_type) == 'box2001') then
                   if (hinit(n) > c0) then
 !                  ! constant slope from 0 to 1 in x direction
                      aicen(i,j,n) = (real(iglob(i), kind=dbl_kind)-p5) &
@@ -1725,6 +1749,15 @@ subroutine set_state_var (nx_block, ny_block, &
 !                                         / (real(ny_global,kind=dbl_kind)) * p5)
                   endif
                   vicen(i,j,n) = hinit(n) * aicen(i,j,n) ! m
+               elseif (trim(ice_data_type) == 'boxslotcyl') then
+                  if (hinit(n) > c0) then
+                   ! slotted cylinder
+                   call boxslotcyl_data_aice(aicen, i, j,        &
+                                             nx_block, ny_block, &
+                                             n,        ainit,    &
+                                             iglob,    jglob)
+                  endif
+                  vicen(i,j,n) = hinit(n) * aicen(i,j,n) ! m
                else
                   vicen(i,j,n) = hinit(n) * ainit(n) ! m
                endif
@@ -1752,6 +1785,18 @@ subroutine set_state_var (nx_block, ny_block, &
 
             enddo               ! ij
          enddo                  ! ncat
+
+         ! velocity initialization for special tests
+         if (trim(ice_data_type) == 'boxslotcyl') then
+            do j = 1, ny_block
+            do i = 1, nx_block
+               call boxslotcyl_data_vel(i,        j,       &
+                                        nx_block, ny_block, &
+                                        iglob,    jglob,   &
+                                        uvel,     vvel)
+            enddo               ! j
+            enddo               ! i
+         endif
       endif                     ! ice_ic
 
       call icepack_warnings_flush(nu_diag)
@@ -1760,6 +1805,133 @@ subroutine set_state_var (nx_block, ny_block, &
 
       end subroutine set_state_var
 
+!=======================================================================
+
+! Set ice concentration for slotted cylinder advection test
+!
+! author: Philippe Blain (ECCC)
+
+      subroutine boxslotcyl_data_aice(aicen, i, j,        &
+                                      nx_block, ny_block, &
+                                      n,        ainit,    &
+                                      iglob,    jglob)
+
+      use ice_constants, only: c0, c2, c5, p3, p166, p75, p5
+      use ice_domain_size, only: nx_global, ny_global, ncat
+      use ice_grid, only: dxrect, dyrect
+
+      integer (kind=int_kind), intent(in) :: &
+         i, j              , & ! local indices
+         nx_block, ny_block, & ! block dimensions
+         iglob(nx_block)   , & ! global indices
+         jglob(ny_block)   , &
+         n                     ! thickness category index
+
+      real (kind=dbl_kind), dimension(ncat) :: &
+         ainit ! initial area
+
+      real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), intent(out) :: &
+         aicen ! concentration of ice
+
+      ! local variables
+
+      logical :: in_slot, in_cyl, in_slotted_cyl
+
+      real (kind=dbl_kind), dimension (2) :: &
+         slot_x, &  ! geometric limits of the slot
+         slot_y
+
+      real (kind=dbl_kind) :: & 
+         diam    , & ! cylinder diameter
+         radius  , & ! cylinder radius
+         center_x, & ! cylinder center
+         center_y, &
+         width   , & ! slot width
+         length      ! slot height
+
+      character(len=*), parameter :: subname = '(boxslotcyl_data_aice)'
+
+      ! Geometric configuration of the slotted cylinder
+      diam     = p3 *dxrect*(nx_global-1)
+      center_x = p5 *dxrect*(nx_global-1)
+      center_y = p75*dxrect*(ny_global-1)
+      radius   = p5*diam
+      width    = p166*diam
+      length   = c5*p166*diam
+
+      slot_x(1) = center_x - width/c2
+      slot_x(2) = center_x + width/c2
+      slot_y(1) = center_y - radius
+      slot_y(2) = center_y + (length - radius)
+
+      ! check if grid point is inside slotted cylinder
+      in_slot = (dxrect*real(iglob(i)-1, kind=dbl_kind) >= slot_x(1)) .and. &
+                (dxrect*real(iglob(i)-1, kind=dbl_kind) <= slot_x(2)) .and. & 
+                (dxrect*real(jglob(j)-1, kind=dbl_kind) >= slot_y(1)) .and. &
+                (dxrect*real(jglob(j)-1, kind=dbl_kind) <= slot_y(2))
+
+      in_cyl  = sqrt((dxrect*real(iglob(i)-1, kind=dbl_kind) - center_x)**c2 + &
+                     (dxrect*real(jglob(j)-1, kind=dbl_kind) - center_y)**c2) <= radius
+
+      in_slotted_cyl = in_cyl .and. .not. in_slot
+
+      if (in_slotted_cyl) then
+         aicen(i,j,n) = ainit(n)
+      else
+         aicen(i,j,n) = c0
+      endif
+
+
+      end subroutine boxslotcyl_data_aice
+
+!=======================================================================
+
+! Set ice velocity for slotted cylinder advection test
+!
+! author: Philippe Blain (ECCC)
+
+      subroutine boxslotcyl_data_vel(i,        j,       &
+                                     nx_block, ny_block, &
+                                     iglob,    jglob,   &
+                                     uvel,     vvel)
+
+      use ice_constants, only: c1, c4, c2, c12, p5, cm_to_m
+      use ice_domain_size, only: nx_global, ny_global
+      use ice_grid, only: dxrect
+
+      integer (kind=int_kind), intent(in) :: &
+         i, j,               & ! local indices
+         nx_block, ny_block, & ! block dimensions
+         iglob(nx_block),    & ! global indices
+         jglob(ny_block)
+
+      real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out) :: &
+         uvel, vvel            ! ice velocity
+
+      ! local variables
+
+      real (kind=dbl_kind) :: &
+         max_vel        , & ! max velocity
+         domain_length  , & ! physical domain length
+         period             ! rotational period
+
+      real (kind=dbl_kind), parameter :: &
+         pi        = c4*atan(c1), & ! pi
+         days_to_s = 86400_dbl_kind
+
+      character(len=*), parameter :: subname = '(boxslotcyl_data_vel)'
+
+      domain_length = dxrect*cm_to_m*nx_global
+      period        = c12*days_to_s            ! 12 days rotational period
+      max_vel       = pi*domain_length/period
+
+      uvel(i,j) =  c2*max_vel*(real(jglob(j), kind=dbl_kind) - p5) &
+                    / real(ny_global - 1, kind=dbl_kind) - max_vel
+      vvel(i,j) = -c2*max_vel*(real(iglob(i), kind=dbl_kind) - p5) &
+                    / real(nx_global - 1, kind=dbl_kind) + max_vel
+
+      end subroutine boxslotcyl_data_vel
+
 !=======================================================================
 
       end module ice_init

cicecore/cicedynB/infrastructure/ice_grid.F90

@@ -77,9 +77,11 @@ module ice_grid
          cxm    , & ! 0.5*HTN - 1.5*HTN
          dxhy   , & ! 0.5*(HTE - HTE)
          dyhx       ! 0.5*(HTN - HTN)
+
+      ! grid dimensions for rectangular grid
       real (kind=dbl_kind), public ::  &
-         dxrect, & !  user_specified spacing (m) in x-direction
-         dyrect    !  user_specified spacing (m) in y-direction
+         dxrect, & !  user_specified spacing (cm) in x-direction (uniform HTN)
+         dyrect    !  user_specified spacing (cm) in y-direction (uniform HTE)
 
       ! Corners of grid boxes for history output
       real (kind=dbl_kind), dimension (:,:,:,:), allocatable, public :: &
@@ -125,13 +127,6 @@ module ice_grid
          lmask_n, & ! northern hemisphere mask
          lmask_s    ! southern hemisphere mask
 
-      ! grid dimensions for rectangular grid
-!      real (kind=dbl_kind), public ::  &
-!         dxrect = 30.e5_dbl_kind   ,&! uniform HTN (cm)
-!         dyrect = 30.e5_dbl_kind     ! uniform HTE (cm)
-!         dxrect = 16.e5_dbl_kind   ,&! uniform HTN (cm)
-!         dyrect = 16.e5_dbl_kind     ! uniform HTE (cm)
-
       real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: &
          rndex_global       ! global index for local subdomain (dbl)
 

cicecore/shared/ice_constants.F90

@@ -71,6 +71,7 @@ module ice_constants
         p15  = 0.15_dbl_kind, &
         p2   = 0.2_dbl_kind, &
         p25  = 0.25_dbl_kind, &
+        p3   = 0.3_dbl_kind, &
         p4   = 0.4_dbl_kind, &
         p5   = 0.5_dbl_kind, &
         p6   = 0.6_dbl_kind, &

configuration/scripts/ice_in

@@ -169,6 +169,7 @@
     ocn_data_type   = 'default'
     bgc_data_type   = 'default'
     fe_data_type    = 'default'
+    ice_data_type   = 'default'
     fyear_init      = 1997
     ycycle          = 1
     atm_data_format = 'bin'

configuration/scripts/machines/Macros.badger_intel

@@ -54,9 +54,9 @@ SCC:= icc
 SFC:= ifort 
 
 ifeq ($(ICE_THREADED), true) 
-   LDFLAGS += -openmp
-   CFLAGS += -openmp
-   FFLAGS += -openmp
+   LDFLAGS += -qopenmp
+   CFLAGS += -qopenmp
+   FFLAGS += -qopenmp
 endif
 
 ### if using parallel I/O, load all 3 libraries.  PIO must be first!

configuration/scripts/machines/env.badger_intel

@@ -33,7 +33,7 @@ setenv ICE_MACHINE_WKDIR /net/scratch3/$user/CICE_RUNS
 setenv ICE_MACHINE_INPUTDATA /usr/projects/climate/eclare/DATA/Consortium
 setenv ICE_MACHINE_BASELINE /net/scratch3/$user/CICE_BASELINE
 setenv ICE_MACHINE_SUBMIT "sbatch "
-setenv ICE_MACHINE_ACCT climateacme
+setenv ICE_MACHINE_ACCT e3sm
 setenv ICE_MACHINE_QUEUE "default"
 setenv ICE_MACHINE_TPNODE 16
 setenv ICE_MACHINE_BLDTHRDS 1

configuration/scripts/options/set_nml.box2001

@@ -23,6 +23,7 @@ ktransport     = -1
 coriolis       = 'constant'
 atmbndy        = 'constant'
 atm_data_type  = 'box2001'
+ice_data_type  = 'box2001'
 calc_strair    = .false.
 restore_ice    = .true.
 f_aice         = 'd'