Merge remote-tracking branch 'Consortium/main' into feature/updcice

.github/workflows/test-cice.yml

@@ -74,6 +74,20 @@ jobs:
            ln -s ${GITHUB_WORKSPACE}/../CICE ${HOME}/cice
 #           ls -al ${HOME}/
 #           ls -al ${GITHUB_WORKSPACE}/
+      - name: check for tabs
+        run: |
+          cd $HOME/cice/cicecore
+          set cnt = 0
+          set ffiles = `find -P . -iname "*.f*"`
+          set cfiles = `find -P . -iname "*.c*"`
+          foreach file ($ffiles $cfiles)
+            set fcnt = `sed -n '/\t/p' $file | wc -l`
+            @ cnt = $cnt + $fcnt
+            if ($fcnt > 0) then
+               echo "TAB found: $fcnt $file"
+            endif
+          end
+          exit $cnt
       - name: setup conda env
         shell: /bin/bash {0}
         run: |

cice.setup

@@ -455,7 +455,22 @@ if ( ${dosuite} == 0 ) then
   set sets = ""
 
 else
-  set tarrays = `echo ${testsuite} | sed 's/,/ /g' | fmt -1 | sort -u`
+  # generate unique set of suites in tarrays in order they are set
+  set tarrays0 = `echo ${testsuite} | sed 's/,/ /g' | fmt -1 `
+  #echo "${0}: tarrays0 = ${tarrays0}"
+  set tarrays = $tarrays0[1]
+  foreach t1 ( ${tarrays0} )
+    set found = 0
+    foreach t2 ( ${tarrays} )
+      if ( ${t1} == ${t2} ) then
+         set found = 1
+      endif
+    end
+    if ( ${found} == 0 ) then
+       set tarrays = ( ${tarrays} ${t1} )
+    endif
+  end
+  #echo "${0}: tarrays = ${tarrays}"
   set testsuitecnt = 0
   foreach tarray ( ${tarrays} )
     @ testsuitecnt = ${testsuitecnt} + 1

cicecore/cicedynB/analysis/ice_history.F90 → cicecore/cicedyn/analysis/ice_history.F90

@@ -102,6 +102,8 @@ subroutine init_hist (dt)
       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=25) :: &
+         gridstr2D, gridstr                  ! temporary string names
       character(len=char_len) :: description
 
       character(len=*), parameter :: subname = '(init_hist)'
@@ -418,10 +420,6 @@ subroutine init_hist (dt)
          f_taubyE = f_tauby
       endif
 
-#ifndef ncdf
-      f_bounds = .false.
-#endif
-
       ! write dimensions for 3D or 4D history variables
       ! note: list of variables checked here is incomplete
       if (f_aicen(1:1) /= 'x' .or. f_vicen(1:1) /= 'x' .or. &
@@ -1307,21 +1305,25 @@ subroutine init_hist (dt)
          select case (grid_ice)
          case('B')
             description = ", on U grid  (NE corner values)"
+            gridstr2d = trim(ustr2D)
+            gridstr   = trim(ucstr)
          case ('CD','C')
             description = ", on T grid"
+            gridstr2d = trim(tstr2D)
+            gridstr   = trim(tcstr)
          end select
 
-         call define_hist_field(n_sig1,"sig1","1",ustr2D, ucstr, &
+         call define_hist_field(n_sig1,"sig1","1",gridstr2d, gridstr, &
              "norm. principal stress 1",                       &
              "sig1 is instantaneous" // trim(description), c1, c0, &
              ns1, f_sig1)
 
-         call define_hist_field(n_sig2,"sig2","1",ustr2D, ucstr, &
+         call define_hist_field(n_sig2,"sig2","1",gridstr2d, gridstr, &
              "norm. principal stress 2",                       &
              "sig2 is instantaneous" // trim(description), c1, c0, &
              ns1, f_sig2)
 
-         call define_hist_field(n_sigP,"sigP","1",ustr2D, ucstr, &
+         call define_hist_field(n_sigP,"sigP","1",gridstr2d, gridstr, &
              "ice pressure",                       &
              "sigP is instantaneous" // trim(description), c1, c0, &
              ns1, f_sigP)

cicecore/cicedynB/dynamics/ice_dyn_eap.F90 → cicecore/cicedyn/dynamics/ice_dyn_eap.F90

@@ -25,6 +25,7 @@ module ice_dyn_eap
           p001, p027, p055, p111, p166, p222, p25, p333
       use ice_fileunits, only: nu_diag, nu_dump_eap, nu_restart_eap
       use ice_exit, only: abort_ice
+      use ice_flux, only: rdg_shear
 !      use ice_timers, only:  &
 !          ice_timer_start, ice_timer_stop, &
 !          timer_tmp1, timer_tmp2, timer_tmp3, timer_tmp4, &
@@ -36,8 +37,7 @@ module ice_dyn_eap
 
       implicit none
       private
-      public :: eap, init_eap, write_restart_eap, read_restart_eap, &
-                alloc_dyn_eap
+      public :: eap, init_eap, write_restart_eap, read_restart_eap
 
       ! Look-up table needed for calculating structure tensor
       integer (int_kind), parameter :: &
@@ -71,42 +71,16 @@ module ice_dyn_eap
       real (kind=dbl_kind) :: &
          puny, pi, pi2, piq, pih
 
-!=======================================================================
-
-      contains
-
-!=======================================================================
-! Allocate space for all variables
-!
-      subroutine alloc_dyn_eap
+      real (kind=dbl_kind), parameter :: &
+         kfriction = 0.45_dbl_kind
 
-      integer (int_kind) :: ierr
+      real (kind=dbl_kind), save :: &
+         invdx, invdy, invda, invsin
 
-      character(len=*), parameter :: subname = '(alloc_dyn_eap)'
 
-      allocate( a11_1        (nx_block,ny_block,max_blocks), &
-                a11_2        (nx_block,ny_block,max_blocks), &
-                a11_3        (nx_block,ny_block,max_blocks), &
-                a11_4        (nx_block,ny_block,max_blocks), &
-                a12_1        (nx_block,ny_block,max_blocks), &
-                a12_2        (nx_block,ny_block,max_blocks), &
-                a12_3        (nx_block,ny_block,max_blocks), &
-                a12_4        (nx_block,ny_block,max_blocks), &
-                e11          (nx_block,ny_block,max_blocks), &
-                e12          (nx_block,ny_block,max_blocks), &
-                e22          (nx_block,ny_block,max_blocks), &
-                yieldstress11(nx_block,ny_block,max_blocks), &
-                yieldstress12(nx_block,ny_block,max_blocks), &
-                yieldstress22(nx_block,ny_block,max_blocks), &
-                s11          (nx_block,ny_block,max_blocks), &
-                s12          (nx_block,ny_block,max_blocks), &
-                s22          (nx_block,ny_block,max_blocks), &
-                a11          (nx_block,ny_block,max_blocks), &
-                a12          (nx_block,ny_block,max_blocks), &
-                stat=ierr)
-      if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory')
+!=======================================================================
 
-      end subroutine alloc_dyn_eap
+      contains
 
 !=======================================================================
 ! Elastic-anisotropic-plastic dynamics driver
@@ -134,7 +108,8 @@ subroutine eap (dt)
           dyn_prep1, dyn_prep2, stepu, dyn_finish, &
           seabed_stress_factor_LKD, seabed_stress_factor_prob, &
           seabed_stress_method, seabed_stress, &
-          stack_fields, unstack_fields, iceTmask, iceUmask
+          stack_fields, unstack_fields, iceTmask, iceUmask, &
+          fld2, fld3, fld4
       use ice_flux, only: rdg_conv, strairxT, strairyT, &
           strairxU, strairyU, uocn, vocn, ss_tltx, ss_tlty, fmU, &
           strtltxU, strtltyU, strocnxU, strocnyU, strintxU, strintyU, taubxU, taubyU, &
@@ -186,11 +161,6 @@ subroutine eap (dt)
          umass      , & ! total mass of ice and snow (u grid)
          umassdti       ! mass of U-cell/dte (kg/m^2 s)
 
-      real (kind=dbl_kind), allocatable :: &
-         fld2(:,:,:,:), & ! temporary for stacking fields for halo update
-         fld3(:,:,:,:), & ! temporary for stacking fields for halo update
-         fld4(:,:,:,:)    ! temporary for stacking fields for halo update
-
       real (kind=dbl_kind), dimension(nx_block,ny_block,8):: &
          strtmp         ! stress combinations for momentum equation
 
@@ -214,10 +184,6 @@ subroutine eap (dt)
       ! Initialize
       !-----------------------------------------------------------------
 
-      allocate(fld2(nx_block,ny_block,2,max_blocks))
-      allocate(fld3(nx_block,ny_block,3,max_blocks))
-      allocate(fld4(nx_block,ny_block,4,max_blocks))
-
        ! This call is needed only if dt changes during runtime.
 !      call set_evp_parameters (dt)
 
@@ -226,7 +192,7 @@ subroutine eap (dt)
          do j = 1, ny_block
          do i = 1, nx_block
             rdg_conv (i,j,iblk) = c0
-!           rdg_shear(i,j,iblk) = c0
+            rdg_shear(i,j,iblk) = c0 ! always zero. Could be moved
             divu (i,j,iblk) = c0
             shear(i,j,iblk) = c0
             e11(i,j,iblk) = c0
@@ -554,7 +520,6 @@ subroutine eap (dt)
 
       enddo                     ! subcycling
 
-      deallocate(fld2,fld3,fld4)
       if (maskhalo_dyn) call ice_HaloDestroy(halo_info_mask)
 
       !-----------------------------------------------------------------
@@ -588,6 +553,8 @@ subroutine init_eap
 
       use ice_blocks, only: nx_block, ny_block
       use ice_domain, only: nblocks
+      use ice_calendar, only: dt_dyn
+      use ice_dyn_shared, only: init_dyn_shared
 
       ! local variables
 
@@ -599,7 +566,7 @@ subroutine init_eap
          eps6 = 1.0e-6_dbl_kind
 
       integer (kind=int_kind) :: &
-         ix, iy, iz, ia
+         ix, iy, iz, ia, ierr
 
       integer (kind=int_kind), parameter :: &
          nz = 100
@@ -609,6 +576,8 @@ subroutine init_eap
          da, dx, dy, dz, &
          phi
 
+      real (kind=dbl_kind) :: invstressconviso
+
       character(len=*), parameter :: subname = '(init_eap)'
 
       call icepack_query_parameters(puny_out=puny, &
@@ -619,6 +588,31 @@ subroutine init_eap
 
       phi = pi/c12 ! diamond shaped floe smaller angle (default phi = 30 deg)
 
+      call init_dyn_shared(dt_dyn)
+
+      allocate( a11_1        (nx_block,ny_block,max_blocks), &
+                a11_2        (nx_block,ny_block,max_blocks), &
+                a11_3        (nx_block,ny_block,max_blocks), &
+                a11_4        (nx_block,ny_block,max_blocks), &
+                a12_1        (nx_block,ny_block,max_blocks), &
+                a12_2        (nx_block,ny_block,max_blocks), &
+                a12_3        (nx_block,ny_block,max_blocks), &
+                a12_4        (nx_block,ny_block,max_blocks), &
+                e11          (nx_block,ny_block,max_blocks), &
+                e12          (nx_block,ny_block,max_blocks), &
+                e22          (nx_block,ny_block,max_blocks), &
+                yieldstress11(nx_block,ny_block,max_blocks), &
+                yieldstress12(nx_block,ny_block,max_blocks), &
+                yieldstress22(nx_block,ny_block,max_blocks), &
+                s11          (nx_block,ny_block,max_blocks), &
+                s12          (nx_block,ny_block,max_blocks), &
+                s22          (nx_block,ny_block,max_blocks), &
+                a11          (nx_block,ny_block,max_blocks), &
+                a12          (nx_block,ny_block,max_blocks), &
+                stat=ierr)
+      if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory')
+
+
       !$OMP PARALLEL DO PRIVATE(iblk,i,j) SCHEDULE(runtime)
       do iblk = 1, nblocks
       do j = 1, ny_block
@@ -640,6 +634,7 @@ subroutine init_eap
          a12_2        (i,j,iblk) = c0
          a12_3        (i,j,iblk) = c0
          a12_4        (i,j,iblk) = c0
+         rdg_shear    (i,j,iblk) = c0
       enddo                     ! i
       enddo                     ! j
       enddo                     ! iblk
@@ -657,6 +652,9 @@ subroutine init_eap
       zinit = -pih
       dy = pi/real(ny_yield-1,kind=dbl_kind)
       yinit = -dy
+      invdx = c1/dx
+      invdy = c1/dy
+      invda = c1/da
 
       do ia=1,na_yield
       do ix=1,nx_yield
@@ -712,6 +710,12 @@ subroutine init_eap
       enddo
       enddo
 
+      ! Factor to maintain the same stress as in EVP (see Section 3)
+      ! Can be set to 1 otherwise
+
+      invstressconviso = c1/(c1+kfriction*kfriction)
+      invsin = c1/sin(pi2/c12) * invstressconviso
+
       end subroutine init_eap
 
 !=======================================================================
@@ -1590,37 +1594,19 @@ subroutine update_stress_rdg (ksub, ndte, divu, tension, &
          rotstemp11s, rotstemp12s, rotstemp22s, &
          sig11, sig12, sig22, &
          sgprm11, sgprm12, sgprm22, &
-         invstressconviso, &
          Angle_denom_gamma,  Angle_denom_alpha, &
          Tany_1, Tany_2, &
          x, y, dx, dy, da, &
          dtemp1, dtemp2, atempprime, &
          kxw, kyw, kaw
 
-      real (kind=dbl_kind), save :: &
-         invdx, invdy, invda, invsin
-
-      logical (kind=log_kind), save :: &
-         first_call = .true.
-
-      real (kind=dbl_kind), parameter :: &
-         kfriction = 0.45_dbl_kind
-
       ! tcraig, temporary, should be moved to namelist
       ! turns on interpolation in stress_rdg
       logical(kind=log_kind), parameter :: &
          interpolate_stress_rdg = .false.
 
       character(len=*), parameter :: subname = '(update_stress_rdg)'
 
-      ! Factor to maintain the same stress as in EVP (see Section 3)
-      ! Can be set to 1 otherwise
-
-      if (first_call) then
-         invstressconviso = c1/(c1+kfriction*kfriction)
-         invsin = c1/sin(pi2/c12) * invstressconviso
-      endif
-
       ! compute eigenvalues, eigenvectors and angles for structure tensor, strain rates
 
       ! 1) structure tensor
@@ -1717,17 +1703,6 @@ subroutine update_stress_rdg (ksub, ndte, divu, tension, &
       if (y > pi) y = y - pi
       if (y < 0)  y = y + pi
 
-      ! Now calculate updated stress tensor
-
-      if (first_call) then
-         dx   = pi/real(nx_yield-1,kind=dbl_kind)
-         dy   = pi/real(ny_yield-1,kind=dbl_kind)
-         da   = p5/real(na_yield-1,kind=dbl_kind)
-         invdx = c1/dx
-         invdy = c1/dy
-         invda = c1/da
-      endif
-
       if (interpolate_stress_rdg) then
 
          ! Interpolated lookup
@@ -1869,8 +1844,6 @@ subroutine update_stress_rdg (ksub, ndte, divu, tension, &
                  + rotstemp22s*dtemp22
       endif
 
-      first_call = .false.
-
       end subroutine update_stress_rdg
 
 !=======================================================================