Implement box model test from 2001 JCP paper (#151)

* Add gbox80 grid configuration

* Add option files for the box2001 case (Hunke, E.C., 2001.  Viscout-plastic sea ice dynamics with the EVP model: linearization issues. J. Comp. Phys. 170, 18-38.)

* Modifications and additions to F90 files required for the box2001 case

* Update box2001 namelist file to have both N/S and E/W boundary types set to 'open'

* modifications to allow for constant coriolis

* modifications for updates to namelist file

* Update ice_in to include coriolis and x/y_spacing variables

* replace x_spacing and y_spacing with dxrect and dyrect in the namelists

* Add  to the box2001 namelist file

* Add new 'thermo' namelist variable to disable thermodynamics.  Add code to set coriolis term to zero.  Update box2001 namelist option file

* Update rectgrid subroutine to set a land mask of 2 grid points around the 80x80 box.

* Update user-guide case settings documentation to include the new namelist variables added for the box2001 test case

* Fix bug in documentation added in previous commit for box2001

* More updates to user guide case settings for box2001 test case

* Remove duplicate  namelist definition in user guide testing

* Add description of box2001 test to the testing documentation

* Typo fix in documentation

* Fix indentation

* Updates to remove the thermo namelist variable, and instead turn off thermodynamics with kdyn=-1.

* Modifications to not call init_flux_* in cice_run if the 2001 test case is run

* Add new kridge and ktransport namelist options, as well as conditionals to turn off ridging and transport

* Resolve namelist errors that were introduced during the merge conflict resolution

* Fix a bug that caused binary data to be written to the log file

* Remove unused variables from ice_step

* Remove unused variables from ice_forcing file

* initialize ice state

* replace land_override with close_boundaries.

* Replace 'default' with 'latitude' as the default coriolis value

* Update documentation to include new namelist variables kridge and ktransport

* Re-add dxrect and dyrect as variables in ice_init.  This was a bug introduced during merge conflict resolution

* correct configuration options

* miscellaneous updates, increased run length

* use strax,y (input directly) instead of strairx,yT (calculated) and set calc_strair=F

Author: mattdturner

cicecore/cicedynB/dynamics/ice_dyn_eap.F90

@@ -28,9 +28,6 @@ module ice_dyn_eap
       use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
       use icepack_intfc, only: icepack_query_parameters
       use icepack_intfc, only: icepack_ice_strength
-#ifdef CICE_IN_NEMO
-      use icepack_intfc, only: calc_strair
-#endif
 
       implicit none
       private
@@ -134,9 +131,6 @@ 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
-#ifdef CICE_IN_NEMO
-      use ice_flux, only: strax, stray
-#endif
       use ice_grid, only: tmask, umask, dxt, dyt, dxhy, dyhx, cxp, cyp, cxm, cym, &
           tarear, uarear, to_ugrid, t2ugrid_vector, u2tgrid_vector
       use ice_state, only: aice, vice, vsno, uvel, vvel, divu, shear, &
@@ -183,6 +177,8 @@ subroutine eap (dt)
       real (kind=dbl_kind), dimension(nx_block,ny_block,8):: &
          strtmp       ! stress combinations for momentum equation
 
+      logical (kind=log_kind) :: calc_strair
+
       integer (kind=int_kind), dimension (nx_block,ny_block,max_blocks) :: &
          icetmask, &  ! ice extent mask (T-cell)
          halomask     ! ice mask for halo update
@@ -259,21 +255,22 @@ subroutine eap (dt)
       call to_ugrid(tmass,umass)
       call to_ugrid(aice_init, aiu)
 
-#ifdef CICE_IN_NEMO
       !----------------------------------------------------------------
-      ! Set wind stress to values supplied via NEMO
+      ! Set wind stress to values supplied via NEMO or other forcing
       ! This wind stress is rotated on u grid and multiplied by aice
       !----------------------------------------------------------------
+      call icepack_query_parameters(calc_strair_out=calc_strair)
+      call icepack_warnings_flush(nu_diag)
+      if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
+         file=__FILE__, line=__LINE__)
+
       if (.not. calc_strair) then       
          strairx(:,:,:) = strax(:,:,:)
          strairy(:,:,:) = stray(:,:,:)
       else
-#endif
          call t2ugrid_vector(strairx)
          call t2ugrid_vector(strairy)
-#ifdef CICE_IN_NEMO
       endif
-#endif
 
 ! tcraig, tcx, turned off this threaded region, in evp, this block and 
 ! the icepack_ice_strength call seems to not be thread safe.  more

cicecore/cicedynB/dynamics/ice_dyn_evp.F90

@@ -45,9 +45,6 @@ module ice_dyn_evp
       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
-#ifdef CICE_IN_NEMO
-      use icepack_intfc, only: calc_strair
-#endif
 
       implicit none
       private
@@ -87,9 +84,6 @@ subroutine evp (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
-#ifdef CICE_IN_NEMO
-      use ice_flux, only: strax, stray
-#endif
       use ice_grid, only: tmask, umask, dxt, dyt, dxhy, dyhx, cxp, cyp, cxm, cym, &
           tarear, uarear, tinyarea, to_ugrid, t2ugrid_vector, u2tgrid_vector, &
           grid_type
@@ -134,6 +128,8 @@ subroutine evp (dt)
       real (kind=dbl_kind), dimension(nx_block,ny_block,8):: &
          strtmp       ! stress combinations for momentum equation
 
+      logical (kind=log_kind) :: calc_strair
+
       integer (kind=int_kind), dimension (nx_block,ny_block,max_blocks) :: &
          icetmask, &  ! ice extent mask (T-cell)
          halomask     ! generic halo mask
@@ -217,21 +213,22 @@ subroutine evp (dt)
       call to_ugrid(tmass,umass)
       call to_ugrid(aice_init, aiu)
 
-#ifdef CICE_IN_NEMO
       !----------------------------------------------------------------
-      ! Set wind stress to values supplied via NEMO
+      ! Set wind stress to values supplied via NEMO or other forcing
       ! This wind stress is rotated on u grid and multiplied by aice
       !----------------------------------------------------------------
+      call icepack_query_parameters(calc_strair_out=calc_strair)
+      call icepack_warnings_flush(nu_diag)
+      if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
+         file=__FILE__, line=__LINE__)
+
       if (.not. calc_strair) then       
          strairx(:,:,:) = strax(:,:,:)
          strairy(:,:,:) = stray(:,:,:)
       else
-#endif
-      call t2ugrid_vector(strairx)
-      call t2ugrid_vector(strairy)
-#ifdef CICE_IN_NEMO
+         call t2ugrid_vector(strairx)
+         call t2ugrid_vector(strairy)
       endif      
-#endif
 
 ! tcraig, tcx, threading here leads to some non-reproducbile results and failures in icepack_ice_strength
 ! need to do more debugging

cicecore/cicedynB/dynamics/ice_dyn_shared.F90

@@ -28,9 +28,14 @@ module ice_dyn_shared
       ! namelist parameters
 
       integer (kind=int_kind), public :: &
-         kdyn     , & ! type of dynamics ( 1 = evp, 2 = eap )
+         kdyn       , & ! type of dynamics ( -1, 0 = off, 1 = evp, 2 = eap )
+         kridge     , & ! set to "-1" to turn off ridging
+         ktransport , & ! set to "-1" to turn off transport
          ndte         ! number of subcycles:  ndte=dt/dte
 
+      character (len=char_len), public :: &
+         coriolis     ! 'constant', 'zero', or 'latitude'
+
       logical (kind=log_kind), public :: &
          revised_evp  ! if true, use revised evp procedure
 
@@ -148,8 +153,13 @@ 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
-         fcor_blk(i,j,iblk) = c2*omega*sin(ULAT(i,j,iblk)) ! 1/s
+         if (trim(coriolis) == 'constant') then
+            fcor_blk(i,j,iblk) = 1.46e-4_dbl_kind ! Hibler 1979, N. Hem; 1/s
+         else if (trim(coriolis) == 'zero') then
+            fcor_blk(i,j,iblk) = 0.0
+         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

cicecore/cicedynB/general/ice_forcing.F90

@@ -29,7 +29,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
@@ -240,6 +240,8 @@ subroutine init_forcing_atmo
          call oned_files
       elseif (trim(atm_data_type) == 'ISPOL') then 
          call ISPOL_files
+      elseif (trim(atm_data_type) == 'box') then
+         call box_data
       endif
 
       end subroutine init_forcing_atmo
@@ -530,6 +532,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
       else    ! default values set in init_flux
          return
       endif
@@ -4405,8 +4409,96 @@ subroutine ocn_data_ispol_init
 
       end subroutine ocn_data_ispol_init
 
+!=======================================================================
+!
+      subroutine box_data
+
+! wind and current fields as in Hunke, JCP 2001
+! authors: Elizabeth Hunke, LANL
+
+      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, wind, rhoa, strax, stray
+      use ice_fileunits, only: nu_diag, nu_forcing
+      use ice_grid, only: uvm
+
+      ! local parameters
+
+      integer (kind=int_kind) :: &
+         iblk, i,j           ! loop indices
+
+      real (kind=dbl_kind) :: &
+          secday, pi , c10, c12, c20, puny, period, pi2, tau
+      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
+
+         uocn(i,j,iblk) = uocn(i,j,iblk) * uvm(i,j,iblk)
+         vocn(i,j,iblk) = vocn(i,j,iblk) * uvm(i,j,iblk)
+
+         ! 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))
+         ! wind stress
+         wind(i,j,iblk) = sqrt(uatm(i,j,iblk)**2 + vatm(i,j,iblk)**2)
+         tau = rhoa(i,j,iblk) * 0.0012_dbl_kind * wind(i,j,iblk)
+         strax(i,j,iblk) = tau * uatm(i,j,iblk)
+         stray(i,j,iblk) = tau * vatm(i,j,iblk)
+
+! 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
 
 !=======================================================================
+