Merge branch 'njeffery_debugZbgc3D' of https://github.com/njeffery/cice…

… into njeffery_debugZbgc3D

cice.setup

@@ -499,6 +499,13 @@ foreach compiler ( $ncompilers )
     set casescr = "${casedir}/casescripts"
     if !( -d ${casescr}) mkdir ${casescr}
 
+    # set default test output as failure
+    if (${docase} == 0) then
+      echo "#---" >! test_output
+      echo "FAIL ${testname_noid} build" >> test_output
+      echo "FAIL ${testname_noid} run" >> test_output
+    endif
+
     # from basic script dir to case
     foreach file (cice.build cice.settings Makefile ice_in makdep.c setup_run_dirs.csh)
       if !(-e ${ICE_SCRIPTS}/$file) then
@@ -703,6 +710,7 @@ EOF1
       # from test options to casescr in case any test time changes are applied
       if (-e ${ICE_SCRIPTS}/tests/test_${test}.files) then
         cp -f -p ${ICE_SCRIPTS}/tests/test_${test}.files ${casescr}
+        cp -f -p ${ICE_SCRIPTS}/tests/comparebfb.csh ${casescr}
         foreach file (`cat ${casescr}/test_${test}.files`)
           if (-e ${ICE_SCRIPTS}/options/$file) then
             cp -f -p ${ICE_SCRIPTS}/options/$file ${casescr}
@@ -758,6 +766,10 @@ EOF2
 
     ${casescr}/parse_settings.sh cice.settings ${fsmods}
     ${casescr}/parse_namelist.sh ice_in ${fimods}
+    if ($status != 0) then
+      echo "${0}: ERROR, parse_namelist.sh aborted"
+      exit -1
+    endif
     source ./cice.settings
     source ./env.${machcomp} -nomodules || exit 2
     ${casescr}/parse_namelist_from_env.sh ice_in
@@ -785,9 +797,9 @@ EOF2
         exit -1
       endif
 
-      # Initial test_output file
-      echo "#---" >! test_output
-      echo "PEND ${testname_noid} " >> test_output
+#      # Initial test_output file
+#      echo "#---" >! test_output
+#      echo "PEND ${testname_noid} " >> test_output
     endif
 
     #------------------------------------------------------------
@@ -837,7 +849,7 @@ cat >> ./${tsdir}/results.csh << EOF
 cat ./results.log
 set pends = \`cat ./results.log | grep PEND | wc -l\`
 set failures = \`cat ./results.log | grep FAIL | wc -l\`
-set success = \`cat ./results.log | grep PASS | wc -l\`
+set success = \`cat ./results.log | grep 'PASS\|COPY' | wc -l\`
 set comments = \`cat ./results.log | grep "#" | wc -l\`
 set alltotal = \`cat ./results.log | wc -l\`
 @ total = \$alltotal - \$comments

cicecore/cicedynB/dynamics/ice_dyn_eap.F90

@@ -19,6 +19,7 @@ module ice_dyn_eap
 
       use ice_kinds_mod
       use ice_blocks, only: nx_block, ny_block
+      use ice_communicate, only: my_task, master_task
       use ice_domain_size, only: max_blocks, ncat
       use ice_constants, only: c0, c1, c2, c3, c12, p1, p2, p5, &
           p001, p027, p055, p111, p166, p222, p25, p333
@@ -27,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
@@ -133,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, &
@@ -182,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
@@ -258,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
@@ -416,7 +414,7 @@ subroutine eap (dt)
       ! stress tensor equation, total surface stress
       !-----------------------------------------------------------------
 
-         !$OMP PARALLEL DO PRIVATE(iblk,strtmp)
+         !$TCXOMP PARALLEL DO PRIVATE(iblk,strtmp)
          do iblk = 1, nblocks
 
 !      call ice_timer_start(timer_tmp1) ! dynamics
@@ -502,7 +500,7 @@ subroutine eap (dt)
             endif
 !      call ice_timer_stop(timer_tmp3) ! dynamics
          enddo
-         !$OMP END PARALLEL DO
+         !$TCXOMP END PARALLEL DO
 
          call ice_timer_start(timer_bound)
          if (maskhalo_dyn) then
@@ -1624,15 +1622,15 @@ subroutine update_stress_rdg (ksub, ndte, divu, tension, &
          gamma, alpha, x, y, dx, dy, da, &
          invdx, invdy, invda, invsin, &
          invleng, dtemp1, dtemp2, atempprime, &
-         puny, pi, pi2, piq
+         puny, pi, pi2, piq, pih
 
       real (kind=dbl_kind), parameter :: &
          kfriction = 0.45_dbl_kind
 
       character(len=*), parameter :: subname = '(update_stress_rdg)'
 
          call icepack_query_parameters(puny_out=puny, &
-            pi_out=pi, pi2_out=pi2, piq_out=piq)
+            pi_out=pi, pi2_out=pi2, piq_out=piq, pih_out=pih)
          call icepack_warnings_flush(nu_diag)
          if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
             file=__FILE__, line=__LINE__)
@@ -1649,8 +1647,12 @@ subroutine update_stress_rdg (ksub, ndte, divu, tension, &
 
          a22 = c1-a11
 
-! gamma: angle between general coordiantes and principal axis 
-         gamma = p5*atan2((c2*a12),(a11 - a22))
+! gamma: angle between general coordiantes and principal axis
+         if ((a11-a22) == c0) then
+           gamma = p5*(pih)
+         else
+           gamma = p5*atan2((c2*a12),(a11 - a22))
+         endif
 
 ! rotational tensor from general coordinates into principal axis
          Q11 = cos(gamma)
@@ -1673,7 +1675,11 @@ subroutine update_stress_rdg (ksub, ndte, divu, tension, &
          dtemp22 = p5*(divu - tension)
 
 ! alpha: angle between general coordiantes and principal axis
-         alpha = p5*atan2((c2*dtemp12),(dtemp11 - dtemp22))
+         if ((dtemp11-dtemp22) == c0) then
+           alpha = p5*(pih)
+         else
+           alpha = p5*atan2((c2*dtemp12),(dtemp11 - dtemp22))
+         endif
 
 ! y: angle between major principal axis of strain rate and structure tensor
 ! to make sure y between 0 and pi/2
@@ -1699,7 +1705,11 @@ subroutine update_stress_rdg (ksub, ndte, divu, tension, &
            invleng = c1/sqrt(dtemp1*dtemp1 + dtemp2*dtemp2)
            dtemp1 = dtemp1*invleng
            dtemp2 = dtemp2*invleng
-           x = atan2(dtemp2,dtemp1)
+           if ((dtemp1) == c0) then
+             x = (pih)
+           else
+             x = atan2(dtemp2,dtemp1)
+           endif
          endif
 
 !echmod to ensure the angle lies between pi/4 and 9 pi/4 
@@ -1928,7 +1938,7 @@ subroutine calc_ffrac (blockno, stressp, stressm, &
 
       real (kind=dbl_kind) :: &
          sigma11, sigma12, sigma22, &
-         gamma, sigma_1, sigma_2, &
+         gamma, sigma_1, sigma_2, pih, &
          Q11, Q12, Q11Q11, Q11Q12, Q12Q12
 
       real (kind=dbl_kind), parameter :: &
@@ -1937,11 +1947,20 @@ subroutine calc_ffrac (blockno, stressp, stressm, &
 
       character(len=*), parameter :: subname = '(calc_ffrac)'
 
+       call icepack_query_parameters(pih_out=pih)
+       call icepack_warnings_flush(nu_diag)
+       if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
+          file=__FILE__, line=__LINE__)
+
        sigma11 = p5*(stressp+stressm) 
        sigma12 = stress12 
        sigma22 = p5*(stressp-stressm) 
 
-       gamma = p5*atan2((c2*sigma12),(sigma11-sigma22))
+       if ((sigma11-sigma22) == c0) then
+         gamma = p5*(pih)
+       else
+         gamma = p5*atan2((c2*sigma12),(sigma11-sigma22))
+       endif
 
 ! rotate tensor to get into sigma principal axis
 
@@ -2023,7 +2042,6 @@ subroutine read_restart_eap()
 
       use ice_blocks, only: nghost
       use ice_boundary, only: ice_HaloUpdate_stress
-      use ice_communicate, only: my_task, master_task
       use ice_constants, only:  &
           field_loc_center, field_type_scalar
       use ice_domain, only: nblocks, halo_info

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
@@ -353,7 +350,7 @@ subroutine evp (dt)
       ! stress tensor equation, total surface stress
       !-----------------------------------------------------------------
 
-         !$OMP PARALLEL DO PRIVATE(iblk,strtmp)
+         !$TCXOMP PARALLEL DO PRIVATE(iblk,strtmp)
          do iblk = 1, nblocks
 
 !            if (trim(yield_curve) == 'ellipse') then
@@ -403,7 +400,7 @@ subroutine evp (dt)
             fld2(:,:,1,iblk) = uvel(:,:,iblk)
             fld2(:,:,2,iblk) = vvel(:,:,iblk)
          enddo
-         !$OMP END PARALLEL DO
+         !$TCXOMP END PARALLEL DO
 
          call ice_timer_start(timer_bound)
          if (maskhalo_dyn) then

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