Implement box model test from 2001 JCP paper

cicecore/cicedynB/dynamics/ice_dyn_shared.F90

@@ -30,6 +30,9 @@ module ice_dyn_shared
          kdyn     , & ! type of dynamics ( 1 = evp, 2 = eap )
          ndte         ! number of subcycles:  ndte=dt/dte
 
+         character (len=char_len), public :: &
+         coriolis     ! 'constant' or 'default'
+
       logical (kind=log_kind), public :: &
          revised_evp  ! if true, use revised evp procedure
 
@@ -130,8 +133,11 @@ subroutine init_evp (dt)
          rdg_shear(i,j,iblk) = c0
 
          ! Coriolis parameter
-!!         fcor_blk(i,j,iblk) = 1.46e-4_dbl_kind ! Hibler 1979, N. Hem; 1/s
+         if (trim(coriolis) == 'constant') then
+         fcor_blk(i,j,iblk) = 1.46e-4_dbl_kind ! Hibler 1979, N. Hem; 1/s
+         else
          fcor_blk(i,j,iblk) = c2*omega*sin(ULAT(i,j,iblk)) ! 1/s
+         endif
 
          ! stress tensor,  kg/s^2
          stressp_1 (i,j,iblk) = c0
@@ -154,6 +160,7 @@ subroutine init_evp (dt)
       enddo                     ! j
       enddo                     ! iblk
       !$OMP END PARALLEL DO
+         write (nu_diag,*) 'fcor_blk(1,1,1)=', fcor_blk(1,1,1)
 
       end subroutine init_evp
 

cicecore/cicedynB/general/ice_forcing.F90

@@ -30,7 +30,7 @@ module ice_forcing
       use ice_timers, only: ice_timer_start, ice_timer_stop, timer_readwrite, &
                             timer_bound
       use ice_arrays_column, only: oceanmixed_ice, restore_bgc
-      use ice_constants, only: c0, c1, c2, c4, c10, c12, c20, &
+      use ice_constants, only: c0, c1, c2, c3, c4, c5, c10, c12, c20, &
                                c180, c365, c1000, c3600
       use ice_constants, only: p001, p01, p1, p25, p5, p6
       use ice_constants, only: cm_to_m
@@ -121,7 +121,7 @@ module ice_forcing
          atm_data_format, & ! 'bin'=binary or 'nc'=netcdf
          ocn_data_format, & ! 'bin'=binary or 'nc'=netcdf
          atm_data_type, & ! 'default', 'monthly', 'ncar', 
-                          ! 'LYq' or 'hadgem' or 'oned'
+                          ! 'LYq' , 'hadgem' , 'oned', or 'box'
          sss_data_type, & ! 'default', 'clim', 'ncar', 'oned'
          sst_data_type, & ! 'default', 'clim', 'ncar', 'oned',
                           ! 'hadgem_sst' or 'hadgem_sst_uvocn'
@@ -492,6 +492,8 @@ subroutine get_forcing_atmo
          call monthly_data
       elseif (trim(atm_data_type) == 'oned') then
          call oned_data
+      elseif (trim(atm_data_type) == 'box') then
+         call box_data(fyear)
       else    ! default values set in init_flux
          return
       endif
@@ -4440,6 +4442,112 @@ end subroutine ocn_data_ispol_init
 
 !=======================================================================
 
-      end module ice_forcing
+!      end module ice_forcing
 
 !=======================================================================
+!=======================================================================
+!
+!BOP
+!
+! !IROUTINE: box_data - define atmospheric data fields 
+!
+! !INTERFACE:
+!
+      subroutine box_data (yr)
+!
+! !DESCRIPTION:
+!
+! wind and current fields as in Hunke, JCP 2001
+!
+! !REVISION HISTORY:
+!
+! authors: Elizabeth Hunke, LANL
+!
+! !USES:
+!
+      use ice_domain, only: nblocks
+      use ice_constants, only: c0, c1, c2, c3, c4, c5, p2
+      use ice_blocks, only: nx_block, ny_block, nghost
+      use ice_flux, only: uocn, vocn, uatm, vatm
+      use ice_fileunits, only: nu_diag, nu_forcing
+
+!
+! !INPUT/OUTPUT PARAMETERS:
+! 
+! EOP
+! 
+!local parameters
+
+      integer (kind=int_kind) :: &
+         iblk, i,j           ! loop indices
+
+      integer (kind=int_kind), intent(in) :: &
+           yr                   ! current forcing year
+      real (kind=dbl_kind) :: &
+          secday, pi , c10, c12, c20, puny, period, pi2
+      call icepack_query_parameters(pi_out=pi, pi2_out=pi2, puny_out=puny)
+      call icepack_query_parameters(secday_out=secday)
+        period = c4*secday 
+
+      do iblk = 1, nblocks
+         do j = 1, ny_block   
+         do i = 1, nx_block   
+
+         ! ocean current
+         ! constant in time---could be initialized in ice_flux.F90
+         uocn(i,j,iblk) =  p2*real(j-nghost, kind=dbl_kind) &
+                            / real(nx_global,kind=dbl_kind) - p1
+         vocn(i,j,iblk) = -p2*real(i-nghost, kind=dbl_kind) &
+                            / real(ny_global,kind=dbl_kind) + p1
+
+         ! wind components
+         uatm(i,j,iblk) = c5 + (sin(pi2*time/period)-c3) &
+                              * sin(pi2*real(i-nghost, kind=dbl_kind)  &
+                                       /real(nx_global,kind=dbl_kind)) &
+                              * sin(pi *real(j-nghost, kind=dbl_kind)  &
+                                       /real(ny_global,kind=dbl_kind))
+         vatm(i,j,iblk) = c5 + (sin(pi2*time/period)-c3) &
+                              * sin(pi *real(i-nghost, kind=dbl_kind)  &
+                                       /real(nx_global,kind=dbl_kind)) &
+                              * sin(pi2*real(j-nghost, kind=dbl_kind)  &
+                                       /real(ny_global,kind=dbl_kind))
+!echmod symm
+!         uocn(i,j,iblk) = c0
+!         vocn(i,j,iblk) = c0
+        ! uatm(i,j,iblk) = p1
+        ! vatm(i,j,iblk) = c0
+!echmod symm
+
+! initialization test
+       ! Diagonal wind vectors 1
+         !uatm(i,j,iblk) = c1 *real(j-nghost, kind=dbl_kind) &
+         !                   / real(ny_global,kind=dbl_kind)
+         !vatm(i,j,iblk) = c1 *real(j-nghost, kind=dbl_kind) &
+         !                   / real(ny_global,kind=dbl_kind)
+
+       ! Diagonal wind vectors 2
+         !uatm(i,j,iblk) = c1 *real(i-nghost, kind=dbl_kind) &
+         !                   / real(nx_global,kind=dbl_kind)
+         !vatm(i,j,iblk) = -c1 *real(i-nghost, kind=dbl_kind) &
+         !                   / real(nx_global,kind=dbl_kind)
+
+       ! Wind in x direction
+        ! uatm(i,j,iblk) = c1 *real(i-nghost, kind=dbl_kind) &
+        !                    / real(nx_global,kind=dbl_kind)
+        ! vatm(i,j,iblk) = c0
+
+       ! Wind in y direction
+        ! uatm(i,j,iblk) = c0
+        ! vatm(i,j,iblk) = c1 *real(j-nghost, kind=dbl_kind) &
+        !                    / real(ny_global,kind=dbl_kind)
+! initialization test
+
+
+         enddo  
+         enddo  
+      enddo ! nblocks
+
+      end subroutine box_data
+
+      end module ice_forcing
+

cicecore/cicedynB/general/ice_init.F90

@@ -76,9 +76,10 @@ subroutine input_data
           atm_data_format, ocn_data_format, &
           sss_data_type,   sst_data_type, ocn_data_dir, &
           oceanmixed_file, restore_sst,   trestore
-      use ice_grid, only: grid_file, gridcpl_file, kmt_file, grid_type, grid_format
+      use ice_grid, only: grid_file, gridcpl_file, kmt_file, grid_type, grid_format, &
+                          grid_spacing, x_spacing, y_spacing
       use ice_dyn_shared, only: ndte, kdyn, revised_evp, yield_curve, &
-                                basalstress, Ktens, e_ratio
+                                basalstress, Ktens, e_ratio, coriolis
       use ice_transport_driver, only: advection
       use ice_restoring, only: restore_ice
 #ifdef CESMCOUPLED
@@ -141,7 +142,8 @@ subroutine input_data
 
       namelist /grid_nml/ &
         grid_format,    grid_type,       grid_file,     kmt_file,       &
-        kcatbound,      gridcpl_file
+        kcatbound,      gridcpl_file, grid_spacing,     x_spacing,      &
+        y_spacing
 
       namelist /thermo_nml/ &
         kitd,           ktherm,          conduct,                       &
@@ -150,7 +152,7 @@ subroutine input_data
 
       namelist /dynamics_nml/ &
         kdyn,           ndte,           revised_evp,    yield_curve,    &
-        advection,                                                      &
+        advection,      coriolis,                                       &
         kstrength,      krdg_partic,    krdg_redist,    mu_rdg,         &
         e_ratio,        Ktens,          Cf,             basalstress
 
@@ -262,6 +264,7 @@ subroutine input_data
       albedo_type = 'default'! or 'constant'
       ktherm = 1             ! 0 = 0-layer, 1 = BL99, 2 = mushy thermo
       conduct = 'bubbly'     ! 'MU71' or 'bubbly' (Pringle et al 2007)
+      coriolis = 'default'   ! latitude dependent, or 'constant'
       calc_Tsfc = .true.     ! calculate surface temperature
       update_ocn_f = .false. ! include fresh water and salt fluxes for frazil
       ustar_min = 0.005      ! minimum friction velocity for ocean heat flux (m/s)
@@ -473,6 +476,9 @@ subroutine input_data
       call broadcast_scalar(pointer_file,       master_task)
       call broadcast_scalar(ice_ic,             master_task)
       call broadcast_scalar(grid_format,        master_task)
+      call broadcast_scalar(grid_spacing,       master_task)
+      call broadcast_scalar(x_spacing,          master_task)
+      call broadcast_scalar(y_spacing,          master_task)
       call broadcast_scalar(grid_type,          master_task)
       call broadcast_scalar(grid_file,          master_task)
       call broadcast_scalar(gridcpl_file,       master_task)
@@ -496,6 +502,7 @@ subroutine input_data
       call broadcast_scalar(shortwave,          master_task)
       call broadcast_scalar(albedo_type,        master_task)
       call broadcast_scalar(ktherm,             master_task)
+      call broadcast_scalar(coriolis,           master_task)
       call broadcast_scalar(conduct,            master_task)
       call broadcast_scalar(R_ice,              master_task)
       call broadcast_scalar(R_pnd,              master_task)
@@ -866,6 +873,7 @@ subroutine input_data
          write(nu_diag,*)    ' yield_curve               = ', &
                                trim(yield_curve)
          write(nu_diag,1020) ' kstrength                 = ', kstrength
+         write(nu_diag,1030) ' coriolis                 = ', coriolis
          write(nu_diag,1020) ' krdg_partic               = ', &
                                krdg_partic
          write(nu_diag,1020) ' krdg_redist               = ', &
@@ -1681,8 +1689,8 @@ subroutine set_state_var (nx_block, ny_block, &
                if (trim(atm_data_type) == 'box') 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) &
-!                                  / (real(nx_global,kind=dbl_kind))
+                     aicen(i,j,n) = (real(iglob(i), kind=dbl_kind)-p5) &
+                                  / (real(nx_global,kind=dbl_kind))
 !                  ! constant slope from 0 to 0.5 in x direction
 !                     aicen(i,j,n) = (real(iglob(i), kind=dbl_kind)-p5) &
 !                                  / (real(nx_global,kind=dbl_kind)) * p5
@@ -1691,12 +1699,12 @@ subroutine set_state_var (nx_block, ny_block, &
 !                                         / (real(nx_global,kind=dbl_kind)) &
 !                                         * (real(jglob(j), kind=dbl_kind)-p5) &
 !                                         / (real(ny_global,kind=dbl_kind)) * p5)
-                     aicen(i,j,n) = max(c0,(real(nx_global, kind=dbl_kind) &
-                                         -  real(iglob(i), kind=dbl_kind)-p5) &
-                                         / (real(nx_global,kind=dbl_kind)) &
-                                         * (real(ny_global, kind=dbl_kind) &
-                                         -  real(jglob(j), kind=dbl_kind)-p5) &
-                                         / (real(ny_global,kind=dbl_kind)) * p5)
+!                     aicen(i,j,n) = max(c0,(real(nx_global, kind=dbl_kind) &
+!                                         -  real(iglob(i), kind=dbl_kind)-p5) &
+!                                         / (real(nx_global,kind=dbl_kind)) &
+!                                         * (real(ny_global, kind=dbl_kind) &
+!                                         -  real(jglob(j), kind=dbl_kind)-p5) &
+!                                         / (real(ny_global,kind=dbl_kind)) * p5)
                   endif
                   vicen(i,j,n) = hinit(n) * aicen(i,j,n) ! m
                else

cicecore/cicedynB/infrastructure/ice_grid.F90

@@ -42,6 +42,7 @@ module ice_grid
          gridcpl_file , & !  input file for POP coupling grid info
          grid_file    , & !  input file for POP grid info
          kmt_file     , & !  input file for POP grid info
+         grid_spacing , & !  default of 30.e3m or set by user in namelist 
          grid_type        !  current options are rectangular (default),
                           !  displaced_pole, tripole, regional
 
@@ -76,6 +77,9 @@ module ice_grid
          cxm    , & ! 0.5*HTN - 1.5*HTN
          dxhy   , & ! 0.5*(HTE - HTE)
          dyhx       ! 0.5*(HTN - HTN)
+      real (kind=dbl_kind), public ::  &
+         x_spacing, & !  user_specified spacing (m) in x-direction
+         y_spacing    !  user_specified spacing (m) in y-direction
 
       ! Corners of grid boxes for history output
       real (kind=dbl_kind), dimension (4,nx_block,ny_block,max_blocks), public :: &
@@ -119,9 +123,11 @@ module ice_grid
          lmask_s    ! southern hemisphere mask
 
       ! grid dimensions for rectangular grid
-      real (kind=dbl_kind), parameter ::  &
+      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 (nx_block,ny_block,max_blocks), public :: &
          rndex_global       ! global index for local subdomain (dbl)
@@ -230,6 +236,12 @@ subroutine init_grid1
         write(nu_diag,'(a26,i6)') '  Block size:  nx_block = ',nx_block
         write(nu_diag,'(a26,i6)') '               ny_block = ',ny_block
       endif
+      if (trim(grid_spacing) == 'user_specified') then
+       dxrect = x_spacing
+       dyrect = y_spacing
+      endif     
+        write(nu_diag,*), 'dxrect= ', dxrect
+        write(nu_diag,*), 'dyrect= ', dyrect
 
       end subroutine init_grid1
 

configuration/scripts/cice_decomp.csh

@@ -44,6 +44,17 @@ else if (${grid} == 'gbox128') then
     set blckx = 8; set blcky = 8
   endif
 
+else if (${grid} == 'gbox80') then
+  set nxglob = 80
+  set nyglob = 80
+  if (${cicepes} <= 1) then
+    set blckx = 80; set blcky = 80
+  else if (${cicepes} <= 8) then
+    set blckx = 20; set blcky = 20
+  else
+    set blckx = 8; set blcky = 8
+  endif
+
 else if (${grid} == 'gx3') then
   set nxglob = 100
   set nyglob = 116

configuration/scripts/ice_in

@@ -48,6 +48,8 @@
     kmt_file     = 'kmt'
     gridcpl_file = 'unknown_gridcpl_file'
     kcatbound    = 0
+    x_spacing    = 30.e3
+    y_spacing    = 30.e3
 /
 
 &domain_nml
@@ -104,6 +106,7 @@
     Ktens           = 0.
     e_ratio         = 2.
     basalstress     = .false.
+    coriolis        = 'default'
 /
 
 &shortwave_nml

configuration/scripts/options/set_env.box2001

@@ -0,0 +1 @@
+setenv NICELYR   1

configuration/scripts/options/set_nml.box2001

@@ -0,0 +1,35 @@
+days_per_year = 360
+npt = 72
+ice_ic = 'default'
+restart = .false.
+restart_ext = .true.
+restart_format = 'bin'
+dumpfreq_n = 5
+histfreq       = 'd','x','x','x','x'
+grid_type    = 'rectangular'
+grid_file = 'grid'
+kmt_file = 'kmt'
+ew_boundary_type  = 'open'
+ns_boundary_type  = 'open'
+tr_iage      = .false.
+tr_FY        = .false.
+tr_lvl       = .false.
+tr_pond_lvl  = .false.
+kitd              = 0
+ktherm            = 0
+kstrength       = 0
+atmbndy         = 'constant'
+atm_data_type   = 'box'
+restore_ice     = .true.
+f_aice         = 'd'
+f_hi           = 'd'
+f_hs           = 'd'
+f_uvel         = 'd'
+f_vvel         = 'd'
+f_uatm         = 'd'
+f_vatm         = 'd'
+f_uocn         = 'd'
+f_vocn         = 'd'
+coriolis       = 'constant'
+x_spacing      = 16.e3
+y_spacing      = 16.e3

configuration/scripts/options/set_nml.gbox80

@@ -0,0 +1,3 @@
+ice_ic = 'default'
+grid_type = 'rectangular'
+atm_data_type = 'box'