From 6449f40c41aa1a5c00096696202d7bd7ebd2a69a Mon Sep 17 00:00:00 2001
From: JFLemieux73 <31927797+JFLemieux73@users.noreply.github.com>
Date: Thu, 11 Jan 2024 19:17:20 +0000
Subject: [PATCH] Add vorticity as a diagnostic output (#924)

* Added new variable vort for vorticity output

* Added calc of diag vorticity for evp, vp and eap for B, C and CD grids

* updated doc and ice_in file for new vorticity variable

* Changed output frequency of vorticity from m to x

* Added f_vort to set_nml.histall and set_nml.histdbg

* Specified location of divu, shear and vort in ice_history.F90

---------

Co-authored-by: Tony Craig <apcraig@users.noreply.github.com>
---
 cicecore/cicedyn/analysis/ice_history.F90     | 20 ++++++++++---
 .../cicedyn/analysis/ice_history_shared.F90   |  9 +++---
 cicecore/cicedyn/dynamics/ice_dyn_eap.F90     | 20 +++++++++++--
 cicecore/cicedyn/dynamics/ice_dyn_evp.F90     | 10 ++++---
 cicecore/cicedyn/dynamics/ice_dyn_shared.F90  | 28 +++++++++++++++++--
 cicecore/cicedyn/dynamics/ice_dyn_vp.F90      |  7 +++--
 cicecore/cicedyn/general/ice_flux.F90         |  3 +-
 cicecore/cicedyn/general/ice_state.F90        |  2 ++
 .../io/io_binary/ice_history_write.F90        |  1 +
 .../io/io_netcdf/ice_history_write.F90        |  1 +
 .../io/io_pio2/ice_history_write.F90          |  1 +
 configuration/scripts/ice_in                  |  1 +
 configuration/scripts/options/set_nml.histall |  1 +
 configuration/scripts/options/set_nml.histdbg |  1 +
 doc/source/cice_index.rst                     |  1 +
 15 files changed, 85 insertions(+), 21 deletions(-)

diff --git a/cicecore/cicedyn/analysis/ice_history.F90 b/cicecore/cicedyn/analysis/ice_history.F90
index 34e5a9131..0243d9861 100644
--- a/cicecore/cicedyn/analysis/ice_history.F90
+++ b/cicecore/cicedyn/analysis/ice_history.F90
@@ -3,7 +3,8 @@
 !
 ! The following variables are currently hard-wired as snapshots
 !   (instantaneous rather than time-averages):
-!   divu, shear, sig1, sig2, sigP, trsig, mlt_onset, frz_onset, hisnap, aisnap
+!   divu, shear, vort, sig1, sig2, sigP, trsig, mlt_onset, 
+!   frz_onset, hisnap, aisnap
 !
 ! Options for histfreq: '1','h','d','m','y','x', where x means that
 !   output stream will not be used (recommended for efficiency).
@@ -597,6 +598,7 @@ subroutine init_hist (dt)
       call broadcast_scalar (f_strength, master_task)
       call broadcast_scalar (f_divu, master_task)
       call broadcast_scalar (f_shear, master_task)
+      call broadcast_scalar (f_vort, master_task)
       call broadcast_scalar (f_sig1, master_task)
       call broadcast_scalar (f_sig2, master_task)
       call broadcast_scalar (f_sigP, master_task)
@@ -1312,14 +1314,19 @@ subroutine init_hist (dt)
 
          call define_hist_field(n_divu,"divu","%/day",tstr2D, tcstr, &
              "strain rate (divergence)",                           &
-             "none", secday*c100, c0,                              &
+             "divu is instantaneous, on T grid", secday*c100, c0,                              &
              ns1, f_divu)
 
          call define_hist_field(n_shear,"shear","%/day",tstr2D, tcstr, &
              "strain rate (shear)",                                  &
-             "none", secday*c100, c0,                                &
+             "shear is instantaneous, on T grid", secday*c100, c0,                                &
              ns1, f_shear)
 
+         call define_hist_field(n_vort,"vort","%/day",tstr2D, tcstr, &
+             "strain rate (vorticity)",                                  &
+             "vort is instantaneous, on T grid", secday*c100, c0,                                &
+             ns1, f_vort)
+
          select case (grid_ice)
          case('B')
             description = ", on U grid  (NE corner values)"
@@ -2623,7 +2630,7 @@ subroutine accum_hist (dt)
          if (f_strength(1:1)/= 'x') &
              call accum_hist_field(n_strength,iblk, strength(:,:,iblk), a2D)
 
-! The following fields (divu, shear, sig1, and sig2) will be smeared
+! The following fields (divu, shear, vort, sig1, and sig2) will be smeared
 !  if averaged over more than a few days.
 ! Snapshots may be more useful (see below).
 
@@ -2631,6 +2638,8 @@ subroutine accum_hist (dt)
 !             call accum_hist_field(n_divu,    iblk, divu(:,:,iblk), a2D)
 !        if (f_shear  (1:1) /= 'x') &
 !             call accum_hist_field(n_shear,   iblk, shear(:,:,iblk), a2D)
+!        if (f_vort  (1:1) /= 'x') &
+!             call accum_hist_field(n_vort,    iblk, vort(:,:,iblk), a2D)
 !        if (f_sig1   (1:1) /= 'x') &
 !             call accum_hist_field(n_sig1,    iblk, sig1(:,:,iblk), a2D)
 !        if (f_sig2   (1:1) /= 'x') &
@@ -3967,6 +3976,7 @@ subroutine accum_hist (dt)
               if (.not. tmask(i,j,iblk)) then ! mask out land points
                  if (n_divu     (ns) /= 0) a2D(i,j,n_divu(ns),     iblk) = spval_dbl
                  if (n_shear    (ns) /= 0) a2D(i,j,n_shear(ns),    iblk) = spval_dbl
+                 if (n_vort     (ns) /= 0) a2D(i,j,n_vort(ns),     iblk) = spval_dbl
                  if (n_sig1     (ns) /= 0) a2D(i,j,n_sig1(ns),     iblk) = spval_dbl
                  if (n_sig2     (ns) /= 0) a2D(i,j,n_sig2(ns),     iblk) = spval_dbl
                  if (n_sigP     (ns) /= 0) a2D(i,j,n_sigP(ns),     iblk) = spval_dbl
@@ -3996,6 +4006,8 @@ subroutine accum_hist (dt)
                        divu (i,j,iblk)*avail_hist_fields(n_divu(ns))%cona
                  if (n_shear    (ns) /= 0) a2D(i,j,n_shear(ns),iblk)     = &
                        shear(i,j,iblk)*avail_hist_fields(n_shear(ns))%cona
+                 if (n_vort     (ns) /= 0) a2D(i,j,n_vort(ns),iblk)      = &
+                       vort(i,j,iblk)*avail_hist_fields(n_vort(ns))%cona
                  if (n_sig1     (ns) /= 0) a2D(i,j,n_sig1(ns),iblk)      = &
                        sig1 (i,j,iblk)*avail_hist_fields(n_sig1(ns))%cona
                  if (n_sig2     (ns) /= 0) a2D(i,j,n_sig2(ns),iblk)      = &
diff --git a/cicecore/cicedyn/analysis/ice_history_shared.F90 b/cicecore/cicedyn/analysis/ice_history_shared.F90
index 16a153c93..2c64bb021 100644
--- a/cicecore/cicedyn/analysis/ice_history_shared.F90
+++ b/cicecore/cicedyn/analysis/ice_history_shared.F90
@@ -4,7 +4,8 @@
 !
 ! The following variables are currently hard-wired as snapshots
 !   (instantaneous rather than time-averages):
-!   divu, shear, sig1, sig2, sigP, trsig, mlt_onset, frz_onset, hisnap, aisnap
+!   divu, shear, vort, sig1, sig2, sigP, trsig, mlt_onset, 
+!   frz_onset, hisnap, aisnap
 !
 ! Options for histfreq: '1','h','d','m','y','x', where x means that
 !   output stream will not be used (recommended for efficiency).
@@ -267,7 +268,7 @@ module ice_history_shared
            f_strocnxE  = 'x', f_strocnyE   = 'x', &
            f_strintxE  = 'x', f_strintyE   = 'x', &
            f_taubxE    = 'x', f_taubyE     = 'x', &
-           f_strength  = 'm', &
+           f_strength  = 'm', f_vort       = 'm', &
            f_divu      = 'm', f_shear      = 'm', &
            f_sig1      = 'm', f_sig2       = 'm', &
            f_sigP      = 'm', &
@@ -434,7 +435,7 @@ module ice_history_shared
 !          f_strocnxE,  f_strocnyE , &
 !          f_strintxE,  f_strintyE , &
 !          f_taubxE,    f_taubyE   , &
-           f_strength,  &
+           f_strength,  f_vort     , &
            f_divu,      f_shear    , &
            f_sig1,      f_sig2     , &
            f_sigP,      &
@@ -626,7 +627,7 @@ module ice_history_shared
            n_strocnxE   , n_strocnyE   , &
            n_strintxE   , n_strintyE   , &
            n_taubxE     , n_taubyE     , &
-           n_strength   , &
+           n_strength   , n_vort       , &
            n_divu       , n_shear      , &
            n_sig1       , n_sig2       , &
            n_sigP       , &
diff --git a/cicecore/cicedyn/dynamics/ice_dyn_eap.F90 b/cicecore/cicedyn/dynamics/ice_dyn_eap.F90
index a8ac62797..cc85d8ab6 100644
--- a/cicecore/cicedyn/dynamics/ice_dyn_eap.F90
+++ b/cicecore/cicedyn/dynamics/ice_dyn_eap.F90
@@ -118,10 +118,10 @@ 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
-      use ice_grid, only: tmask, umask, dxT, dyT, &
+      use ice_grid, only: tmask, umask, dxT, dyT, dxU, dyU, &
           tarear, uarear, grid_average_X2Y, &
           grid_atm_dynu, grid_atm_dynv, grid_ocn_dynu, grid_ocn_dynv
-      use ice_state, only: aice, aiU, vice, vsno, uvel, vvel, divu, shear, &
+      use ice_state, only: aice, aiU, vice, vsno, uvel, vvel, divu, shear, vort, &
           aice_init, aice0, aicen, vicen, strength
       use ice_timers, only: timer_dynamics, timer_bound, &
           ice_timer_start, ice_timer_stop
@@ -195,6 +195,7 @@ subroutine eap (dt)
             rdg_shear(i,j,iblk) = c0 ! always zero. Could be moved
             divu (i,j,iblk) = c0
             shear(i,j,iblk) = c0
+            vort(i,j,iblk)  = c0
             e11(i,j,iblk) = c0
             e12(i,j,iblk) = c0
             e22(i,j,iblk) = c0
@@ -433,6 +434,7 @@ subroutine eap (dt)
                               arlx1i,               denom1,               &
                               uvel      (:,:,iblk), vvel      (:,:,iblk), &
                               dxT       (:,:,iblk), dyT       (:,:,iblk), &
+                              dxU       (:,:,iblk), dyU       (:,:,iblk), &
                               dxhy      (:,:,iblk), dyhx      (:,:,iblk), &
                               cxp       (:,:,iblk), cyp       (:,:,iblk), &
                               cxm       (:,:,iblk), cym       (:,:,iblk), &
@@ -448,6 +450,7 @@ subroutine eap (dt)
                               stress12_1(:,:,iblk), stress12_2(:,:,iblk), &
                               stress12_3(:,:,iblk), stress12_4(:,:,iblk), &
                               shear     (:,:,iblk), divu      (:,:,iblk), &
+                              vort      (:,:,iblk),                       &
                               e11       (:,:,iblk), e12       (:,:,iblk), &
                               e22       (:,:,iblk),                       &
                               s11       (:,:,iblk), s12       (:,:,iblk), &
@@ -1162,6 +1165,7 @@ subroutine stress_eap  (nx_block,   ny_block,       &
                               arlx1i,     denom1,         &
                               uvel,       vvel,           &
                               dxT,        dyT,            &
+                              dxU,        dyU,            &
                               dxhy,       dyhx,           &
                               cxp,        cyp,            &
                               cxm,        cym,            &
@@ -1175,6 +1179,7 @@ subroutine stress_eap  (nx_block,   ny_block,       &
                               stress12_1, stress12_2,     &
                               stress12_3, stress12_4,     &
                               shear,      divu,           &
+                              vort,                       &
                               e11,        e12,            &
                               e22,                        &
                               s11,        s12,            &
@@ -1206,6 +1211,8 @@ subroutine stress_eap  (nx_block,   ny_block,       &
          vvel     , & ! y-component of velocity (m/s)
          dxT      , & ! width of T-cell through the middle (m)
          dyT      , & ! height of T-cell through the middle (m)
+         dxU      , & ! width of U-cell through the middle (m)
+         dyU      , & ! height of U-cell through the middle (m)
          dxhy     , & ! 0.5*(HTE - HTW)
          dyhx     , & ! 0.5*(HTN - HTS)
          cyp      , & ! 1.5*HTE - 0.5*HTW
@@ -1226,6 +1233,7 @@ subroutine stress_eap  (nx_block,   ny_block,       &
       real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
          shear    , & ! strain rate II component (1/s)
          divu     , & ! strain rate I component, velocity divergence (1/s)
+         vort     , & ! vorticity (1/s)
          e11      , & ! components of strain rate tensor (1/s)
          e12      , & !
          e22      , & !
@@ -1255,6 +1263,7 @@ subroutine stress_eap  (nx_block,   ny_block,       &
          divune, divunw, divuse, divusw            , & ! divergence
          tensionne, tensionnw, tensionse, tensionsw, & ! tension
          shearne, shearnw, shearse, shearsw        , & ! shearing
+         dvdxn, dvdxs, dudye, dudyw                , & ! for vorticity calc
          ssigpn, ssigps, ssigpe, ssigpw            , &
          ssigmn, ssigms, ssigme, ssigmw            , &
          ssig12n, ssig12s, ssig12e, ssig12w        , &
@@ -1357,6 +1366,13 @@ subroutine stress_eap  (nx_block,   ny_block,       &
                  (tensionne + tensionnw + tensionse + tensionsw)**2 &
                 +  (shearne +   shearnw +   shearse +   shearsw)**2)
 
+            ! vorticity
+            dvdxn = dyU(i,j)*vvel(i,j) - dyU(i-1,j)*vvel(i-1,j)
+            dvdxs = dyU(i,j-1)*vvel(i,j-1) - dyU(i-1,j-1)*vvel(i-1,j-1)
+            dudye = dxU(i,j)*uvel(i,j) - dxU(i,j-1)*uvel(i,j-1)
+            dudyw = dxU(i-1,j)*uvel(i-1,j) - dxU(i-1,j-1)*uvel(i-1,j-1)
+            vort(i,j) = p5*tarear(i,j)*(dvdxn + dvdxs - dudye - dudyw)
+
             divu(i,j) = p25*(divune + divunw + divuse + divusw) * tarear(i,j)
             rdg_conv(i,j)  = -min(p25*(alpharne + alpharnw &
                                      + alpharsw + alpharse),c0) * tarear(i,j)
diff --git a/cicecore/cicedyn/dynamics/ice_dyn_evp.F90 b/cicecore/cicedyn/dynamics/ice_dyn_evp.F90
index 6a71d6a14..301a89916 100644
--- a/cicecore/cicedyn/dynamics/ice_dyn_evp.F90
+++ b/cicecore/cicedyn/dynamics/ice_dyn_evp.F90
@@ -283,7 +283,7 @@ subroutine evp (dt)
           grid_type, grid_ice, &
           grid_atm_dynu, grid_atm_dynv, grid_ocn_dynu, grid_ocn_dynv
       use ice_state, only: aice, aiU, vice, vsno, uvel, vvel, uvelN, vvelN, &
-          uvelE, vvelE, divu, shear, &
+          uvelE, vvelE, divu, shear, vort, &
           aice_init, aice0, aicen, vicen, strength
       use ice_timers, only: timer_dynamics, timer_bound, &
           ice_timer_start, ice_timer_stop, timer_evp
@@ -385,6 +385,7 @@ subroutine evp (dt)
             rdg_shear(i,j,iblk) = c0
             divu (i,j,iblk) = c0
             shear(i,j,iblk) = c0
+            vort (i,j,iblk) = c0
          enddo
          enddo
 
@@ -921,9 +922,10 @@ subroutine evp (dt)
                                indxTi    (:,iblk), indxTj     (:,iblk), &
                                uvel    (:,:,iblk), vvel     (:,:,iblk), &
                                dxT     (:,:,iblk), dyT      (:,:,iblk), &
+                               dxU     (:,:,iblk), dyU      (:,:,iblk), &
                                cxp     (:,:,iblk), cyp      (:,:,iblk), &
                                cxm     (:,:,iblk), cym      (:,:,iblk), &
-                               tarear  (:,:,iblk),                      &
+                               tarear  (:,:,iblk), vort     (:,:,iblk), &
                                shear   (:,:,iblk), divu     (:,:,iblk), &
                                rdg_conv(:,:,iblk), rdg_shear(:,:,iblk) )
          enddo
@@ -1109,7 +1111,7 @@ subroutine evp (dt)
                                   dxN     (:,:,iblk), dyE      (:,:,iblk), &
                                   dxT     (:,:,iblk), dyT      (:,:,iblk), &
                                   tarear  (:,:,iblk), uarea    (:,:,iblk), &
-                                  shearU    (:,:,iblk),                    &
+                                  shearU  (:,:,iblk), vort     (:,:,iblk), &
                                   shear   (:,:,iblk), divu     (:,:,iblk), &
                                   rdg_conv(:,:,iblk), rdg_shear(:,:,iblk))
          enddo
@@ -1299,7 +1301,7 @@ subroutine evp (dt)
                                    uvelN   (:,:,iblk), vvelN    (:,:,iblk), &
                                    dxN     (:,:,iblk), dyE      (:,:,iblk), &
                                    dxT     (:,:,iblk), dyT      (:,:,iblk), &
-                                   tarear  (:,:,iblk),                      &
+                                   tarear  (:,:,iblk), vort     (:,:,iblk), &
                                    shear   (:,:,iblk), divu     (:,:,iblk), &
                                    rdg_conv(:,:,iblk), rdg_shear(:,:,iblk))
          enddo
diff --git a/cicecore/cicedyn/dynamics/ice_dyn_shared.F90 b/cicecore/cicedyn/dynamics/ice_dyn_shared.F90
index d3f819f20..015c925a6 100644
--- a/cicecore/cicedyn/dynamics/ice_dyn_shared.F90
+++ b/cicecore/cicedyn/dynamics/ice_dyn_shared.F90
@@ -1711,9 +1711,10 @@ subroutine deformations (nx_block,   ny_block,   &
                                indxTi,     indxTj,     &
                                uvel,       vvel,       &
                                dxT,        dyT,        &
+                               dxU,        dyU,        &
                                cxp,        cyp,        &
                                cxm,        cym,        &
-                               tarear,                 &
+                               tarear,     vort,       &
                                shear,      divu,       &
                                rdg_conv,   rdg_shear )
 
@@ -1732,6 +1733,8 @@ subroutine deformations (nx_block,   ny_block,   &
          vvel     , & ! y-component of velocity (m/s)
          dxT      , & ! width of T-cell through the middle (m)
          dyT      , & ! height of T-cell through the middle (m)
+         dxU      , & ! width of U-cell through the middle (m)
+         dyU      , & ! height of U-cell through the middle (m)
          cyp      , & ! 1.5*HTE - 0.5*HTW
          cxp      , & ! 1.5*HTN - 0.5*HTS
          cym      , & ! 0.5*HTE - 1.5*HTW
@@ -1739,6 +1742,7 @@ subroutine deformations (nx_block,   ny_block,   &
          tarear       ! 1/tarea
 
       real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
+         vort     , & ! vorticity (1/s)   
          shear    , & ! strain rate II component (1/s)
          divu     , & ! strain rate I component, velocity divergence (1/s)
          rdg_conv , & ! convergence term for ridging (1/s)
@@ -1756,6 +1760,9 @@ subroutine deformations (nx_block,   ny_block,   &
         Deltane, Deltanw, Deltase, Deltasw        , & ! Delta
         tmp                                           ! useful combination
 
+      real (kind=dbl_kind) :: &                       ! at edges for vorticity calc :
+         dvdxn, dvdxs, dudye, dudyw                   ! dvdx and dudy terms on edges
+
       character(len=*), parameter :: subname = '(deformations)'
 
       do ij = 1, icellT
@@ -1794,6 +1801,13 @@ subroutine deformations (nx_block,   ny_block,   &
                       (tensionne + tensionnw + tensionse + tensionsw)**2 + &
                       (shearne   + shearnw   + shearse   + shearsw  )**2)
 
+         ! vorticity
+         dvdxn = dyU(i,j)*vvel(i,j) - dyU(i-1,j)*vvel(i-1,j)
+         dvdxs = dyU(i,j-1)*vvel(i,j-1) - dyU(i-1,j-1)*vvel(i-1,j-1)
+         dudye = dxU(i,j)*uvel(i,j) - dxU(i,j-1)*uvel(i,j-1)
+         dudyw = dxU(i-1,j)*uvel(i-1,j) - dxU(i-1,j-1)*uvel(i-1,j-1)
+         vort(i,j) = p5*tarear(i,j)*(dvdxn + dvdxs - dudye - dudyw)
+
       enddo                     ! ij
 
       end subroutine deformations
@@ -1811,7 +1825,7 @@ subroutine deformationsCD_T (nx_block,   ny_block,   &
                                    uvelN,      vvelN,      &
                                    dxN,        dyE,        &
                                    dxT,        dyT,        &
-                                   tarear,                 &
+                                   tarear,     vort,       &
                                    shear,      divu,       &
                                    rdg_conv,   rdg_shear )
 
@@ -1837,6 +1851,7 @@ subroutine deformationsCD_T (nx_block,   ny_block,   &
          tarear       ! 1/tarea
 
       real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
+         vort     , & ! vorticity (1/s)
          shear    , & ! strain rate II component (1/s)
          divu     , & ! strain rate I component, velocity divergence (1/s)
          rdg_conv , & ! convergence term for ridging (1/s)
@@ -1888,6 +1903,9 @@ subroutine deformationsCD_T (nx_block,   ny_block,   &
          ! diagnostic only
          ! shear = sqrt(tension**2 + shearing**2)
          shear(i,j) = tarear(i,j)*sqrt( tensionT(i,j)**2 + shearT(i,j)**2 )
+         ! vorticity
+         vort (i,j) = tarear(i,j)*( ( dyE(i,j)*vvelE(i,j) - dyE(i-1,j)*vvelE(i-1,j) ) &
+                                  - ( dxN(i,j)*uvelN(i,j) - dxN(i,j-1)*uvelN(i,j-1)) )
 
       enddo                     ! ij
 
@@ -1908,7 +1926,7 @@ subroutine deformationsC_T (nx_block,   ny_block,   &
                                 dxN,        dyE,        &
                                 dxT,        dyT,        &
                                 tarear,     uarea,      &
-                                shearU,                 &
+                                shearU,     vort,       &
                                 shear,      divu,       &
                                 rdg_conv,   rdg_shear )
 
@@ -1936,6 +1954,7 @@ subroutine deformationsC_T (nx_block,   ny_block,   &
          shearU       ! shearU
 
       real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
+         vort     , & ! vorticity (1/s)
          shear    , & ! strain rate II component (1/s)
          divu     , & ! strain rate I component, velocity divergence (1/s)
          rdg_conv , & ! convergence term for ridging (1/s)
@@ -1999,6 +2018,9 @@ subroutine deformationsC_T (nx_block,   ny_block,   &
          ! diagnostic only...maybe we dont want to use shearTsqr here????
          ! shear = sqrt(tension**2 + shearing**2)
          shear(i,j) = tarear(i,j)*sqrt( tensionT(i,j)**2 + shearT(i,j)**2 )
+         ! vorticity
+         vort (i,j) = tarear(i,j)*( ( dyE(i,j)*vvelE(i,j) - dyE(i-1,j)*vvelE(i-1,j) ) &
+                                  - ( dxN(i,j)*uvelN(i,j) - dxN(i,j-1)*uvelN(i,j-1)) )
 
       enddo                     ! ij
 
diff --git a/cicecore/cicedyn/dynamics/ice_dyn_vp.F90 b/cicecore/cicedyn/dynamics/ice_dyn_vp.F90
index 477d19515..0d04bf974 100644
--- a/cicecore/cicedyn/dynamics/ice_dyn_vp.F90
+++ b/cicecore/cicedyn/dynamics/ice_dyn_vp.F90
@@ -167,10 +167,10 @@ subroutine implicit_solver (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
-      use ice_grid, only: tmask, umask, dxT, dyT, &
+      use ice_grid, only: tmask, umask, dxT, dyT, dxU, dyU, &
           tarear, grid_type, grid_average_X2Y, &
           grid_atm_dynu, grid_atm_dynv, grid_ocn_dynu, grid_ocn_dynv
-      use ice_state, only: aice, aiU, vice, vsno, uvel, vvel, divu, shear, &
+      use ice_state, only: aice, aiU, vice, vsno, uvel, vvel, divu, shear, vort, &
           aice_init, aice0, aicen, vicen, strength
       use ice_timers, only: timer_dynamics, timer_bound, &
           ice_timer_start, ice_timer_stop
@@ -521,9 +521,10 @@ subroutine implicit_solver (dt)
                             indxTi      (:,iblk), indxTj      (:,iblk), &
                             uvel      (:,:,iblk), vvel      (:,:,iblk), &
                             dxT       (:,:,iblk), dyT       (:,:,iblk), &
+                            dxU       (:,:,iblk), dyU       (:,:,iblk), &
                             cxp       (:,:,iblk), cyp       (:,:,iblk), &
                             cxm       (:,:,iblk), cym       (:,:,iblk), &
-                            tarear    (:,:,iblk),                       &
+                            tarear    (:,:,iblk), vort      (:,:,iblk), &
                             shear     (:,:,iblk), divu      (:,:,iblk), &
                             rdg_conv  (:,:,iblk), rdg_shear (:,:,iblk))
       enddo
diff --git a/cicecore/cicedyn/general/ice_flux.F90 b/cicecore/cicedyn/general/ice_flux.F90
index 4c37a0696..2d61bf642 100644
--- a/cicecore/cicedyn/general/ice_flux.F90
+++ b/cicecore/cicedyn/general/ice_flux.F90
@@ -1022,7 +1022,7 @@ end subroutine init_history_therm
 
       subroutine init_history_dyn
 
-      use ice_state, only: aice, vice, trcr, strength, divu, shear
+      use ice_state, only: aice, vice, trcr, strength, divu, shear, vort
       use ice_grid,  only: grid_ice
 
       logical (kind=log_kind) :: &
@@ -1043,6 +1043,7 @@ subroutine init_history_dyn
       sig2    (:,:,:) = c0
       divu    (:,:,:) = c0
       shear   (:,:,:) = c0
+      vort    (:,:,:) = c0
       taubxU  (:,:,:) = c0
       taubyU  (:,:,:) = c0
       strength (:,:,:) = c0
diff --git a/cicecore/cicedyn/general/ice_state.F90 b/cicecore/cicedyn/general/ice_state.F90
index 862f0a8bc..21ddf562c 100644
--- a/cicecore/cicedyn/general/ice_state.F90
+++ b/cicecore/cicedyn/general/ice_state.F90
@@ -116,6 +116,7 @@ module ice_state
          vvelN    , & ! y-component of velocity on N grid (m/s)
          divu     , & ! strain rate I component, velocity divergence (1/s)
          shear    , & ! strain rate II component (1/s)
+         vort     , & ! vorticity (1/s)
          strength     ! ice strength (N/m)
 
       !-----------------------------------------------------------------
@@ -163,6 +164,7 @@ subroutine alloc_state
          vvelN     (nx_block,ny_block,max_blocks) , & ! y-component of velocity on N grid (m/s)
          divu      (nx_block,ny_block,max_blocks) , & ! strain rate I component, velocity divergence (1/s)
          shear     (nx_block,ny_block,max_blocks) , & ! strain rate II component (1/s)
+         vort      (nx_block,ny_block,max_blocks) , & ! vorticity (1/s)
          strength  (nx_block,ny_block,max_blocks) , & ! ice strength (N/m)
          aice_init (nx_block,ny_block,max_blocks) , & ! initial concentration of ice, for diagnostics
          aicen     (nx_block,ny_block,ncat,max_blocks) , & ! concentration of ice
diff --git a/cicecore/cicedyn/infrastructure/io/io_binary/ice_history_write.F90 b/cicecore/cicedyn/infrastructure/io/io_binary/ice_history_write.F90
index 526d0d96d..b16d00f07 100644
--- a/cicecore/cicedyn/infrastructure/io/io_binary/ice_history_write.F90
+++ b/cicecore/cicedyn/infrastructure/io/io_binary/ice_history_write.F90
@@ -160,6 +160,7 @@ subroutine ice_write_hist(ns)
             if (histfreq(ns) == '1' .or. .not. hist_avg(ns)     &
                 .or. write_ic                                   &
                 .or. n==n_divu(ns)      .or. n==n_shear(ns)     &  ! snapshots
+                .or. n==n_vort(ns)                              &  ! snapshots
                 .or. n==n_sig1(ns)      .or. n==n_sig2(ns)      &
                 .or. n==n_sigP(ns)      .or. n==n_trsig(ns)     &
                 .or. n==n_sistreave(ns) .or. n==n_sistremax(ns) &
diff --git a/cicecore/cicedyn/infrastructure/io/io_netcdf/ice_history_write.F90 b/cicecore/cicedyn/infrastructure/io/io_netcdf/ice_history_write.F90
index 51d76a6f4..4a0e86233 100644
--- a/cicecore/cicedyn/infrastructure/io/io_netcdf/ice_history_write.F90
+++ b/cicecore/cicedyn/infrastructure/io/io_netcdf/ice_history_write.F90
@@ -1325,6 +1325,7 @@ subroutine ice_write_hist_attrs(ncid, varid, hfield, ns)
           .or. write_ic                                   &
           .or.TRIM(hfield%vname(1:4))=='divu' &
           .or.TRIM(hfield%vname(1:5))=='shear' &
+          .or.TRIM(hfield%vname(1:4))=='vort' &
           .or.TRIM(hfield%vname(1:4))=='sig1' &
           .or.TRIM(hfield%vname(1:4))=='sig2' &
           .or.TRIM(hfield%vname(1:4))=='sigP' &
diff --git a/cicecore/cicedyn/infrastructure/io/io_pio2/ice_history_write.F90 b/cicecore/cicedyn/infrastructure/io/io_pio2/ice_history_write.F90
index cf2f40521..650005a83 100644
--- a/cicecore/cicedyn/infrastructure/io/io_pio2/ice_history_write.F90
+++ b/cicecore/cicedyn/infrastructure/io/io_pio2/ice_history_write.F90
@@ -1286,6 +1286,7 @@ subroutine ice_write_hist_attrs(File, varid, hfield, ns)
           .or. write_ic                                   &
           .or.TRIM(hfield%vname(1:4))=='divu' &
           .or.TRIM(hfield%vname(1:5))=='shear' &
+          .or.TRIM(hfield%vname(1:4))=='vort' &
           .or.TRIM(hfield%vname(1:4))=='sig1' &
           .or.TRIM(hfield%vname(1:4))=='sig2' &
           .or.TRIM(hfield%vname(1:4))=='sigP' &
diff --git a/configuration/scripts/ice_in b/configuration/scripts/ice_in
index 85f502683..e33d16c18 100644
--- a/configuration/scripts/ice_in
+++ b/configuration/scripts/ice_in
@@ -571,6 +571,7 @@
     f_strength     = 'm'
     f_divu         = 'm'
     f_shear        = 'm'
+    f_vort         = 'x'
     f_sig1         = 'm' 
     f_sig2         = 'm' 
     f_sigP         = 'm' 
diff --git a/configuration/scripts/options/set_nml.histall b/configuration/scripts/options/set_nml.histall
index 78932cba8..83421aca0 100644
--- a/configuration/scripts/options/set_nml.histall
+++ b/configuration/scripts/options/set_nml.histall
@@ -106,6 +106,7 @@
     f_strength     = 'md'
     f_divu         = 'md'
     f_shear        = 'md'
+    f_vort         = 'md'			
     f_sig1         = 'md' 
     f_sig2         = 'md' 
     f_sigP         = 'md' 
diff --git a/configuration/scripts/options/set_nml.histdbg b/configuration/scripts/options/set_nml.histdbg
index 43ae8e566..a70e734e5 100644
--- a/configuration/scripts/options/set_nml.histdbg
+++ b/configuration/scripts/options/set_nml.histdbg
@@ -106,6 +106,7 @@
     f_strength     = 'md1'
     f_divu         = 'md1'
     f_shear        = 'md1'
+    f_vort         = 'md1'
     f_sig1         = 'md1' 
     f_sig2         = 'md1' 
     f_sigP         = 'md1' 
diff --git a/doc/source/cice_index.rst b/doc/source/cice_index.rst
index dae10eda4..1249feb08 100644
--- a/doc/source/cice_index.rst
+++ b/doc/source/cice_index.rst
@@ -745,6 +745,7 @@ section :ref:`tabnamelist`.
    "visc_method", "method for calculating viscosities (‘avg_strength’ or ‘avg_zeta’)", "avg_zeta"
    "vocn", "ocean current in the y-direction", "m/s"
    "vonkar", "von Karman constant", "0.4"
+   "vort", "vorticity", "1/s"
    "vraftn", "volume of rafted ice", "m"
    "vrdgn", "volume of ridged ice", "m"
    "vredistrn", "redistribution function: fraction of new ridge volume", ""