Add grid_ocn, grid_atm feature (CICE-Consortium#47)

* - Refactor strair and strocn implementation a bit
  - Note strocnxT, strocnyT are per ice area and different units from other strocn[x,y] variables.  Added task to C-grid clean issue.
  - Remove strocnxT_f, strocnyT_f from dyn_finish interface and compute separately
  - Remove strairxU/T, strairyU/T arguments from dyn_prep1, compute strairxU, strairyU more cleanly
  - Update variable documentation in code

- Remove straxE, strayE, straxN, strayN, probably won't need

- Update grid_average_X2Y
  - Add grid_average_X2Y for 'NE2T' and 'NE2U' as an overloaded interface
  - Add grid_average_X2Y "A" implementation which is unmasked normalized weighted average.  This is like S but ignores masks
  - Add grid_average_X2Y explicit implementation (dir,work1,wght1,mask1,work2) for 'S' and 'A' averaging options
  - Eliminate "in-place" operation, not really needed, can only cause confusion, require all averaging from one variable to a different variable.

- Update gridavgchk unit test to test new grid_average_X2Y options

- Identify bug related to location of uocn,vocn variables.  Used on U grid but variables on T grid.  Added task to C-grid clean issue.

* - Add grid_atm and grid_ocn to namelist as well as grid_*_[thrm,dynu,dynv] to support flexible grid definitions
- Refactor grid_average_X2Y interfaces to better support flexibility wrt grids
- Update gridavgchk unit test

* - Update [u,v]ocn, ss_tlt[x,y] implementation to improve flexibility and use grid_atm and grid_ocn info
  - Migrate averaging and memory allocation of U, N, E fields to dynamics

- Update history capability to support grid_atm and grid_ocn values for some history variables

* - Update u/vocn and ss_tltx/y usage to support the grid_ocn value

- Fix bug in init/get_forcing_ocn, affects only cases with atm_data_type=box2001 and ocn_data_type/=box2001

- Update set_nml files for box cases to more clearly specify atm/ocn/ice_data_type and define grid_ocn where needed

* update documentation

* recover lost grid_file initialization line

* update uatm, vatm handling

* - Update strax, stray implementation to support flexible grid interpolation
- Update grid_average_X2Y internal interface names
- Update documentation

* - Rename grid_system to grid_ice
- Update documentation

* rename grid_system to grid_ice

* rename grid_system to grid_ice

* update comment for uocn, vocn

cicecore/cicedynB/analysis/ice_diagnostics.F90

@@ -127,7 +127,7 @@ subroutine runtime_diags (dt)
           alvdr_init, alvdf_init, alidr_init, alidf_init
       use ice_flux_bgc, only: faero_atm, faero_ocn, fiso_atm, fiso_ocn
       use ice_global_reductions, only: global_sum, global_sum_prod, global_maxval
-      use ice_grid, only: lmask_n, lmask_s, tarean, tareas, grid_system
+      use ice_grid, only: lmask_n, lmask_s, tarean, tareas, grid_ice
       use ice_state   ! everything
 ! tcraig, this is likely to cause circular dependency because ice_prescribed_mod is high level routine
 #ifdef CESMCOUPLED
@@ -297,7 +297,7 @@ subroutine runtime_diags (dt)
          enddo
       enddo
       ! Eventually do energy diagnostic on T points.
-!     if (grid_system == 'CD') then
+!     if (grid_ice == 'CD') then
 !     !$OMP PARALLEL DO PRIVATE(iblk,i,j)
 !     do iblk = 1, nblocks
 !        do j = 1, ny_block
@@ -403,7 +403,7 @@ subroutine runtime_diags (dt)
          enddo
       enddo
       !$OMP END PARALLEL DO
-      if (grid_system == 'CD') then
+      if (grid_ice == 'CD') then
       !$OMP PARALLEL DO PRIVATE(iblk,i,j)
       do iblk = 1, nblocks
          do j = 1, ny_block
@@ -1663,7 +1663,7 @@ end subroutine debug_ice
 
       subroutine print_state(plabel,i,j,iblk)
 
-      use ice_grid, only: grid_system
+      use ice_grid, only: grid_ice
       use ice_blocks, only: block, get_block
       use ice_domain, only: blocks_ice
       use ice_domain_size, only: ncat, nilyr, nslyr, nfsd
@@ -1790,7 +1790,7 @@ subroutine print_state(plabel,i,j,iblk)
 
       write(nu_diag,*) 'uvel(i,j)',uvel(i,j,iblk)
       write(nu_diag,*) 'vvel(i,j)',vvel(i,j,iblk)
-      if (grid_system == 'CD') then
+      if (grid_ice == 'CD') then
          write(nu_diag,*) 'uvelE(i,j)',uvelE(i,j,iblk)
          write(nu_diag,*) 'vvelE(i,j)',vvelE(i,j,iblk)
          write(nu_diag,*) 'uvelN(i,j)',uvelN(i,j,iblk)
@@ -1843,7 +1843,7 @@ end subroutine print_state
 
       subroutine print_points_state(plabel,ilabel)
 
-      use ice_grid, only: grid_system
+      use ice_grid, only: grid_ice
       use ice_blocks, only: block, get_block
       use ice_domain, only: blocks_ice
       use ice_domain_size, only: ncat, nilyr, nslyr
@@ -1940,7 +1940,7 @@ subroutine print_points_state(plabel,ilabel)
 
       write(nu_diag,*) trim(llabel),'uvel=',uvel(i,j,iblk)
       write(nu_diag,*) trim(llabel),'vvel=',vvel(i,j,iblk)
-      if (grid_system == 'CD') then
+      if (grid_ice == 'CD') then
          write(nu_diag,*) trim(llabel),'uvelE=',uvelE(i,j,iblk)
          write(nu_diag,*) trim(llabel),'vvelE=',vvelE(i,j,iblk)
          write(nu_diag,*) trim(llabel),'uvelN=',uvelN(i,j,iblk)

cicecore/cicedynB/analysis/ice_history.F90

@@ -68,7 +68,9 @@ subroutine init_hist (dt)
       use ice_domain_size, only: max_blocks, max_nstrm, nilyr, nslyr, nblyr, ncat, nfsd
       use ice_dyn_shared, only: kdyn
       use ice_flux, only: mlt_onset, frz_onset, albcnt, snwcnt
-      use ice_grid, only: grid_system
+      use ice_grid, only: grid_ice, &
+          grid_atm_thrm, grid_atm_dynu, grid_atm_dynv, &
+          grid_ocn_thrm, grid_ocn_dynu, grid_ocn_dynv
       use ice_history_shared ! everything
       use ice_history_mechred, only: init_hist_mechred_2D, init_hist_mechred_3Dc
       use ice_history_pond, only: init_hist_pond_2D, init_hist_pond_3Dc
@@ -94,9 +96,108 @@ subroutine init_hist (dt)
       integer (kind=int_kind), dimension(max_nstrm) :: &
          ntmp
       integer (kind=int_kind) :: nml_error ! namelist i/o error flag
+      character (len=25) :: &
+         str2D_gat, str2d_gau, str2d_gav, &  ! dimensions for t, u, v atm grid (ga)
+         str2D_got, str2d_gou, str2d_gov     ! dimensions for t, u, v ocn grid (go)
+      character (len=25) :: &
+         cstr_gat, cstr_gau, cstr_gav, &     ! mask area name for t, u, v atm grid (ga)
+         cstr_got, cstr_gou, cstr_gov        ! mask area name for t, u, v ocn grid (go)
       character(len=char_len) :: description
       character(len=*), parameter :: subname = '(init_hist)'
 
+      !-----------------------------------------------------------------
+      ! set atm/ocn forcing grids
+      !-----------------------------------------------------------------
+
+      !--- ATM ---
+
+      if     (grid_atm_thrm == 'T') then
+         str2D_gat = tstr2D
+         cstr_gat  = tcstr
+      elseif (grid_atm_thrm == 'U') then
+         str2D_gat = ustr2D
+         cstr_gat  = ucstr
+      elseif (grid_atm_thrm == 'N') then
+         str2D_gat = nstr2D
+         cstr_gat  = ncstr
+      elseif (grid_atm_thrm == 'E') then
+         str2D_gat = estr2D
+         cstr_gat  = ecstr
+      endif
+
+      if     (grid_atm_dynu == 'T') then
+         str2D_gau = tstr2D
+         cstr_gau  = tcstr
+      elseif (grid_atm_dynu == 'U') then
+         str2D_gau = ustr2D
+         cstr_gau  = ucstr
+      elseif (grid_atm_dynu == 'N') then
+         str2D_gau = nstr2D
+         cstr_gau  = ncstr
+      elseif (grid_atm_dynu == 'E') then
+         str2D_gau = estr2D
+         cstr_gau  = ecstr
+      endif
+
+      if     (grid_atm_dynv == 'T') then
+         str2D_gav = tstr2D
+         cstr_gav  = tcstr
+      elseif (grid_atm_dynv == 'U') then
+         str2D_gav = ustr2D
+         cstr_gav  = ucstr
+      elseif (grid_atm_dynv == 'N') then
+         str2D_gav = nstr2D
+         cstr_gav  = ncstr
+      elseif (grid_atm_dynv == 'E') then
+         str2D_gav = estr2D
+         cstr_gav  = ecstr
+      endif
+
+      !--- OCN ---
+
+      if     (grid_ocn_thrm == 'T') then
+         str2D_got = tstr2D
+         cstr_got  = tcstr
+      elseif (grid_ocn_thrm == 'U') then
+         str2D_got = ustr2D
+         cstr_got  = ucstr
+      elseif (grid_ocn_thrm == 'N') then
+         str2D_got = nstr2D
+         cstr_got  = ncstr
+      elseif (grid_ocn_thrm == 'E') then
+         str2D_got = estr2D
+         cstr_got  = ecstr
+      endif
+
+      if     (grid_ocn_dynu == 'T') then
+         str2D_gou = tstr2D
+         cstr_gou  = tcstr
+      elseif (grid_ocn_dynu == 'U') then
+         str2D_gou = ustr2D
+         cstr_gou  = ucstr
+      elseif (grid_ocn_dynu == 'N') then
+         str2D_gou = nstr2D
+         cstr_gou  = ncstr
+      elseif (grid_ocn_dynu == 'E') then
+         str2D_gou = estr2D
+         cstr_gou  = ecstr
+      endif
+
+      if     (grid_ocn_dynv == 'T') then
+         str2D_gov = tstr2D
+         cstr_gov  = tcstr
+      elseif (grid_ocn_dynv == 'U') then
+         str2D_gov = ustr2D
+         cstr_gov  = ucstr
+      elseif (grid_ocn_dynv == 'N') then
+         str2D_gov = nstr2D
+         cstr_gov  = ncstr
+      elseif (grid_ocn_dynv == 'E') then
+         str2D_gov = estr2D
+         cstr_gov  = ecstr
+      endif
+
+
       !-----------------------------------------------------------------
       ! set history dimensions
       !-----------------------------------------------------------------
@@ -279,7 +380,7 @@ subroutine init_hist (dt)
          f_sispeed = f_CMIP
       endif
 
-      if (grid_system == 'CD') then
+      if (grid_ice == 'CD') then
          f_uvelE = f_uvel
          f_vvelE = f_vvel
          f_icespdE = f_icespd
@@ -684,22 +785,22 @@ subroutine init_hist (dt)
              "vector direction - coming from", c1, c0,                 &
              ns1, f_icedir)
 
-         call define_hist_field(n_uatm,"uatm","m/s",ustr2D, ucstr,  &
+         call define_hist_field(n_uatm,"uatm","m/s",str2D_gau, cstr_gau,  &
              "atm velocity (x)",                                  &
              "positive is x direction on U grid", c1, c0,         &
              ns1, f_uatm)
 
-         call define_hist_field(n_vatm,"vatm","m/s",ustr2D, ucstr,  &
+         call define_hist_field(n_vatm,"vatm","m/s",str2D_gav, cstr_gav,  &
              "atm velocity (y)",                                  &
              "positive is y direction on U grid", c1, c0,         &
              ns1, f_vatm)
 
-         call define_hist_field(n_atmspd,"atmspd","m/s",ustr2D, ucstr, &
+         call define_hist_field(n_atmspd,"atmspd","m/s",str2D_gau, cstr_gau, &
              "atmosphere wind speed",                                  &
              "vector magnitude", c1, c0,                               &
              ns1, f_atmspd)
 
-         call define_hist_field(n_atmdir,"atmdir","deg",ustr2D, ucstr, &
+         call define_hist_field(n_atmdir,"atmdir","deg",str2D_gau, cstr_gau, &
              "atmosphere wind direction",                              &
              "vector direction - coming from", c1, c0,                 &
              ns1, f_atmdir)
@@ -754,22 +855,22 @@ subroutine init_hist (dt)
              "none", c1, c0,                                   &
              ns1, f_sss)
 
-         call define_hist_field(n_uocn,"uocn","m/s",ustr2D, ucstr, &
+         call define_hist_field(n_uocn,"uocn","m/s",str2D_gou, cstr_gou, &
              "ocean current (x)",                                &
              "positive is x direction on U grid", c1, c0,        &
              ns1, f_uocn)
 
-         call define_hist_field(n_vocn,"vocn","m/s",ustr2D, ucstr, &
+         call define_hist_field(n_vocn,"vocn","m/s",str2D_gov, cstr_gov, &
              "ocean current (y)",                                &
              "positive is y direction on U grid", c1, c0,        &
              ns1, f_vocn)
 
-         call define_hist_field(n_ocnspd,"ocnspd","m/s",ustr2D, ucstr, &
+         call define_hist_field(n_ocnspd,"ocnspd","m/s",str2D_gou, cstr_gou, &
              "ocean current speed",                                    &
              "vector magnitude", c1, c0,                               &
              ns1, f_ocnspd)
 
-         call define_hist_field(n_ocndir,"ocndir","deg",ustr2D, ucstr, &
+         call define_hist_field(n_ocndir,"ocndir","deg",str2D_gou, cstr_gou, &
              "ocean current direction",                                &
              "vector direction - going to", c1, c0,                    &
              ns1, f_ocndir)
@@ -1200,7 +1301,7 @@ subroutine init_hist (dt)
              "none", secday*c100, c0,                                &
              ns1, f_shear)
 
-         select case (grid_system)
+         select case (grid_ice)
          case('B')
             description = ", on U grid  (NE corner values)"
          case ('CD')
@@ -1987,7 +2088,7 @@ subroutine accum_hist (dt)
       use ice_blocks, only: block, get_block, nx_block, ny_block
       use ice_domain, only: blocks_ice, nblocks
       use ice_domain_size, only: nfsd
-      use ice_grid, only: tmask, lmask_n, lmask_s, dxu, dyu, grid_system
+      use ice_grid, only: tmask, lmask_n, lmask_s, dxu, dyu, grid_ice
       use ice_calendar, only: new_year, write_history, &
                               write_ic, timesecs, histfreq, nstreams, mmonth, &
                               new_month
@@ -4297,7 +4398,7 @@ subroutine accum_hist (dt)
       !---------------------------------------------------------------
 
             ! compute sig1 and sig2
-            select case (grid_system)
+            select case (grid_ice)
             case('B')
                call principal_stress (nx_block,  ny_block,  &
                                       stressp_1 (:,:,iblk), &