Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add vorticity as a diagnostic output #924

Merged
merged 7 commits into from
Jan 11, 2024
Merged

Add vorticity as a diagnostic output #924

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
a098554
Added new variable vort for vorticity output
JFLemieux73 Dec 22, 2023
005ba4a
Added calc of diag vorticity for evp, vp and eap for B, C and CD grids
JFLemieux73 Dec 22, 2023
43b6847
updated doc and ice_in file for new vorticity variable
JFLemieux73 Dec 22, 2023
9d3107f
Changed output frequency of vorticity from m to x
JFLemieux73 Jan 3, 2024
1345aa5
Added f_vort to set_nml.histall and set_nml.histdbg
JFLemieux73 Jan 5, 2024
0425cbd
Specified location of divu, shear and vort in ice_history.F90
JFLemieux73 Jan 11, 2024
baf25d5
Merge branch 'main' into vorticity
apcraig Jan 11, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 16 additions & 4 deletions cicecore/cicedyn/analysis/ice_history.F90
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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)

Comment on lines +1369 to +1375
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here the computation of vorticity is the same as in ice_dyn_shared::deformations, right ?

If so, I would prefer we compute the vorticity in a separate subroutine, and call that subroutine from ice_dyn_eap::stress_eap and ice_dyn_shared::deformations. This would reduce code duplication.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But then it is the same thing for divu and shear...the same computations are done in ice_dyn_eap and ice_dyn_shared. It has been like this for a long time...I say we keep it simple the way it is done right now.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm okay with this as it is, for now. We do need to do a massive cleanup of the whole dynamics code, to fix these kinds of things. There is a lot of duplication for historical reasons.

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
Open in desktop
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
Open in desktop
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
Loading