Skip to content

Commit

Permalink
Add vorticity as a diagnostic output (CICE-Consortium#924)
Browse files Browse the repository at this point in the history 
* 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 <[email protected]>
  • Loading branch information
@JFLemieux73 @apcraig @DeniseWorthen
2 people authored and DeniseWorthen committed Apr 7, 2024
1 parent f514034 commit bfbb221
Show file tree
Hide file tree
Showing 15 changed files with 85 additions and 21 deletions.
20 changes: 16 additions & 4 deletions cicecore/cicedyn/analysis/ice_history.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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).
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)"
Expand Down Expand Up @@ -2623,14 +2630,16 @@ 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).

! if (f_divu (1:1) /= 'x') &
! 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') &
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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) = &
Expand Down
9 changes: 5 additions & 4 deletions cicecore/cicedyn/analysis/ice_history_shared.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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).
Expand Down Expand Up @@ -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', &
Expand Down Expand Up @@ -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, &
Expand Down Expand Up @@ -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 , &
Expand Down
20 changes: 18 additions & 2 deletions cicecore/cicedyn/dynamics/ice_dyn_eap.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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), &
Expand All @@ -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), &
Expand Down Expand Up @@ -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, &
Expand All @@ -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, &
Expand Down Expand Up @@ -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
Expand All @@ -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 , & !
Expand Down Expand Up @@ -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 , &
Expand Down Expand Up @@ -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)
Expand Down
10 changes: 6 additions & 4 deletions cicecore/cicedyn/dynamics/ice_dyn_evp.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
28 changes: 25 additions & 3 deletions cicecore/cicedyn/dynamics/ice_dyn_shared.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 )

Expand All @@ -1732,13 +1733,16 @@ 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
cxm , & ! 0.5*HTN - 1.5*HTS
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)
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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 )

Expand All @@ -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)
Expand Down Expand Up @@ -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

Expand All @@ -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 )

Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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

Expand Down
Loading

0 comments on commit bfbb221

Please sign in to comment.