Skip to content

Commit

Permalink
Update KE diagnostic (CICE-Consortium#784)
Browse files Browse the repository at this point in the history 
* Update KE diagnostic

- Average U and V fields to T grid for KE calculation for B, C, and CD grid,
  raw U and V were combined with T grid variables in prior implementation.

Also
- Update max speed calculation for C grid, calc speed at both E and N points
- Refactor highfreq option u,v passed into icepack for windstress calculation
  - This is bit-for-bit
  - Create local variables uvelT_icep and vvelT_icep to store highfreq info
  - Compute [uvelT,vvelT]_icep at start of therm calc instead of end of dynamics calc
  - No longer need to recalc [uvelT,vvelT]_icep on restart
- Update some indentation in ice_diagnostics.F90

* Add grid average NE2TA, NE2UA
  • Loading branch information
@apcraig
apcraig authored Nov 7, 2022
1 parent 3820cde commit aa1e066
Show file tree
Hide file tree
Showing 5 changed files with 207 additions and 149 deletions.
156 changes: 74 additions & 82 deletions cicecore/cicedynB/analysis/ice_diagnostics.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -127,7 +127,7 @@ subroutine runtime_diags (dt)
alvdr_init, alvdf_init, alidr_init, alidf_init
use ice_flux_bgc, only: faero_atm, faero_ocn, fiso_atm, fiso_ocn
use ice_global_reductions, only: global_sum, global_sum_prod, global_maxval
use ice_grid, only: lmask_n, lmask_s, tarean, tareas, grid_ice
use ice_grid, only: lmask_n, lmask_s, tarean, tareas, grid_ice, grid_average_X2Y
use ice_state ! everything
! tcraig, this is likely to cause circular dependency because ice_prescribed_mod is high level routine
#ifdef CESMCOUPLED
Expand Down Expand Up @@ -199,7 +199,8 @@ subroutine runtime_diags (dt)
prsnwavg, prhosavg, psmicetot, psmliqtot, psmtot

real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
work1, work2
uvelT, vvelT, & ! u,v on T points
work1, work2 ! temporary

real (kind=dbl_kind), parameter :: &
maxval_spval = -0.9_dbl_kind*HUGE(0.0_dbl_kind) ! spval to detect
Expand All @@ -209,9 +210,6 @@ subroutine runtime_diags (dt)
! returns -HUGE which we want to avoid writing. The
! return value is checked against maxval_spval before writing.

! real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
! uvelT, vvelT

character(len=*), parameter :: subname = '(runtime_diags)'

call icepack_query_parameters(ktherm_out=ktherm, calc_Tsfc_out=calc_Tsfc)
Expand Down Expand Up @@ -284,51 +282,45 @@ subroutine runtime_diags (dt)
if (tr_pond_topo) then
!$OMP PARALLEL DO PRIVATE(iblk,i,j,n)
do iblk = 1, nblocks
do j = 1, ny_block
do i = 1, nx_block
work1(i,j,iblk) = c0
do n = 1, ncat
work1(i,j,iblk) = work1(i,j,iblk) &
+ aicen(i,j,n,iblk) &
* trcrn(i,j,nt_apnd,n,iblk) &
* trcrn(i,j,nt_hpnd,n,iblk)
do j = 1, ny_block
do i = 1, nx_block
work1(i,j,iblk) = c0
do n = 1, ncat
work1(i,j,iblk) = work1(i,j,iblk) &
+ aicen(i,j,n,iblk) &
* trcrn(i,j,nt_apnd,n,iblk) &
* trcrn(i,j,nt_hpnd,n,iblk)
enddo
enddo
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
ptotn = global_sum(work1, distrb_info, field_loc_center, tarean)
ptots = global_sum(work1, distrb_info, field_loc_center, tareas)
endif

! total ice-snow kinetic energy
! total ice-snow kinetic energy, on T points.
if (grid_ice == 'B') then
call grid_average_X2Y('A',uvel ,'U',uvelT,'T')
call grid_average_X2Y('A',vvel ,'U',vvelT,'T')
elseif (grid_ice == 'C') then
call grid_average_X2Y('A',uvelE,'E',uvelT,'T')
call grid_average_X2Y('A',vvelN,'N',vvelT,'T')
elseif (grid_ice == 'CD') then
call grid_average_X2Y('A',uvelE,'E',uvelN,'N',uvelT,'T')
call grid_average_X2Y('A',vvelE,'E',vvelN,'N',vvelT,'T')
endif

!$OMP PARALLEL DO PRIVATE(iblk,i,j)
do iblk = 1, nblocks
do j = 1, ny_block
do i = 1, nx_block
work1(i,j,iblk) = p5 &
* (rhos*vsno(i,j,iblk) + rhoi*vice(i,j,iblk)) &
* (uvel(i,j,iblk)**2 + vvel(i,j,iblk)**2)
work1(i,j,iblk) = p5 * (rhos*vsno(i,j,iblk) + rhoi*vice(i,j,iblk)) &
* (uvelT(i,j,iblk)**2 + vvelT(i,j,iblk)**2)
enddo
enddo
enddo
! Eventually do energy diagnostic on T points.
! if (grid_ice == 'CD') then
! !$OMP PARALLEL DO PRIVATE(iblk,i,j)
! do iblk = 1, nblocks
! do j = 1, ny_block
! do i = 1, nx_block
! call grid_average_X2Y('E2TS',uvelE,uvelT)
! call grid_average_X2Y('N2TS',vvelN,vvelT)
! work1(i,j,iblk) = p5 &
! * (rhos*vsno(i,j,iblk) + rhoi*vice(i,j,iblk)) &
! * (uvelT(i,j,iblk)*uvelT(i,j,iblk) &
! + vvelT(i,j,iblk)*vvelT(i,j,iblk))
! enddo
! enddo
! enddo
! endif
! !$OMP END PARALLEL DO
!$OMP END PARALLEL DO
ketotn = c0
ketots = c0
ketotn = global_sum(work1, distrb_info, field_loc_center, tarean)
Expand Down Expand Up @@ -420,21 +412,22 @@ subroutine runtime_diags (dt)
do iblk = 1, nblocks
do j = 1, ny_block
do i = 1, nx_block
work1(i,j,iblk) = max(sqrt(uvelE(i,j,iblk)**2 &
+ vvelE(i,j,iblk)**2), &
sqrt(uvelN(i,j,iblk)**2 &
+ vvelN(i,j,iblk)**2))
work1(i,j,iblk) = max(sqrt(uvelE(i,j,iblk)**2 + vvelE(i,j,iblk)**2), &
sqrt(uvelN(i,j,iblk)**2 + vvelN(i,j,iblk)**2))
enddo
enddo
enddo
!$OMP END PARALLEL DO
elseif (grid_ice == 'C') then
! map uvelE to N and vvelN to E then compute max on E and N
call grid_average_X2Y('A',uvelE,'E',work1,'N') ! work1 =~ uvelN
call grid_average_X2Y('A',vvelN,'N',work2,'E') ! work2 =~ vvelE
!$OMP PARALLEL DO PRIVATE(iblk,i,j)
do iblk = 1, nblocks
do j = 1, ny_block
do i = 1, nx_block
work1(i,j,iblk) = sqrt(uvelE(i,j,iblk)**2 &
+ vvelN(i,j,iblk)**2)
work1(i,j,iblk) = max(sqrt(uvelE(i,j,iblk)**2 + work2(i,j,iblk)**2), &
sqrt(work1(i,j,iblk)**2 + vvelN(i,j,iblk)**2))
enddo
enddo
enddo
Expand All @@ -444,8 +437,7 @@ subroutine runtime_diags (dt)
do iblk = 1, nblocks
do j = 1, ny_block
do i = 1, nx_block
work1(i,j,iblk) = sqrt(uvel(i,j,iblk)**2 &
+ vvel(i,j,iblk)**2)
work1(i,j,iblk) = sqrt(uvel(i,j,iblk)**2 + vvel(i,j,iblk)**2)
enddo
enddo
enddo
Expand All @@ -466,31 +458,31 @@ subroutine runtime_diags (dt)
if (umaxn > umax_stab) then
!$OMP PARALLEL DO PRIVATE(iblk,i,j)
do iblk = 1, nblocks
do j = 1, ny_block
do i = 1, nx_block
if (abs(work1(i,j,iblk) - umaxn) < puny) then
write(nu_diag,*) ' '
write(nu_diag,*) 'Warning, large ice speed'
write(nu_diag,*) 'my_task, iblk, i, j, umaxn:', &
my_task, iblk, i, j, umaxn
endif
enddo
enddo
do j = 1, ny_block
do i = 1, nx_block
if (abs(work1(i,j,iblk) - umaxn) < puny) then
write(nu_diag,*) ' '
write(nu_diag,*) 'Warning, large ice speed'
write(nu_diag,*) 'my_task, iblk, i, j, umaxn:', &
my_task, iblk, i, j, umaxn
endif
enddo
enddo
enddo
!$OMP END PARALLEL DO
elseif (umaxs > umax_stab) then
!$OMP PARALLEL DO PRIVATE(iblk,i,j)
do iblk = 1, nblocks
do j = 1, ny_block
do i = 1, nx_block
if (abs(work1(i,j,iblk) - umaxs) < puny) then
write(nu_diag,*) ' '
write(nu_diag,*) 'Warning, large ice speed'
write(nu_diag,*) 'my_task, iblk, i, j, umaxs:', &
my_task, iblk, i, j, umaxs
endif
enddo
enddo
do j = 1, ny_block
do i = 1, nx_block
if (abs(work1(i,j,iblk) - umaxs) < puny) then
write(nu_diag,*) ' '
write(nu_diag,*) 'Warning, large ice speed'
write(nu_diag,*) 'my_task, iblk, i, j, umaxs:', &
my_task, iblk, i, j, umaxs
endif
enddo
enddo
enddo
!$OMP END PARALLEL DO
endif ! umax
Expand Down Expand Up @@ -1357,14 +1349,14 @@ subroutine init_mass_diags
do n=1,n_aero
!$OMP PARALLEL DO PRIVATE(iblk,i,j)
do iblk = 1, nblocks
do j = 1, ny_block
do i = 1, nx_block
work1(i,j,iblk) = trcr(i,j,nt_aero +4*(n-1),iblk)*vsno(i,j,iblk) &
+ trcr(i,j,nt_aero+1+4*(n-1),iblk)*vsno(i,j,iblk) &
+ trcr(i,j,nt_aero+2+4*(n-1),iblk)*vice(i,j,iblk) &
+ trcr(i,j,nt_aero+3+4*(n-1),iblk)*vice(i,j,iblk)
enddo
enddo
do j = 1, ny_block
do i = 1, nx_block
work1(i,j,iblk) = trcr(i,j,nt_aero +4*(n-1),iblk)*vsno(i,j,iblk) &
+ trcr(i,j,nt_aero+1+4*(n-1),iblk)*vsno(i,j,iblk) &
+ trcr(i,j,nt_aero+2+4*(n-1),iblk)*vice(i,j,iblk) &
+ trcr(i,j,nt_aero+3+4*(n-1),iblk)*vice(i,j,iblk)
enddo
enddo
enddo
!$OMP END PARALLEL DO
totaeron(n)= global_sum(work1, distrb_info, field_loc_center, tarean)
Expand All @@ -1377,17 +1369,17 @@ subroutine init_mass_diags
totps = c0
!$OMP PARALLEL DO PRIVATE(iblk,i,j,n)
do iblk = 1, nblocks
do j = 1, ny_block
do i = 1, nx_block
work1(i,j,iblk) = c0
do n = 1, ncat
work1(i,j,iblk) = work1(i,j,iblk) &
+ aicen(i,j,n,iblk) &
* trcrn(i,j,nt_apnd,n,iblk) &
* trcrn(i,j,nt_hpnd,n,iblk)
do j = 1, ny_block
do i = 1, nx_block
work1(i,j,iblk) = c0
do n = 1, ncat
work1(i,j,iblk) = work1(i,j,iblk) &
+ aicen(i,j,n,iblk) &
* trcrn(i,j,nt_apnd,n,iblk) &
* trcrn(i,j,nt_hpnd,n,iblk)
enddo
enddo
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
totpn = global_sum(work1, distrb_info, field_loc_center, tarean)
Expand Down
4 changes: 0 additions & 4 deletions cicecore/cicedynB/general/ice_state.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -110,8 +110,6 @@ module ice_state
public :: &
uvel , & ! x-component of velocity on U grid (m/s)
vvel , & ! y-component of velocity on U grid (m/s)
uvelT , & ! x-component of velocity on T grid (m/s)
vvelT , & ! y-component of velocity on T grid (m/s)
uvelE , & ! x-component of velocity on E grid (m/s)
vvelE , & ! y-component of velocity on E grid (m/s)
uvelN , & ! x-component of velocity on N grid (m/s)
Expand Down Expand Up @@ -159,8 +157,6 @@ subroutine alloc_state
aice0 (nx_block,ny_block,max_blocks) , & ! concentration of open water
uvel (nx_block,ny_block,max_blocks) , & ! x-component of velocity on U grid (m/s)
vvel (nx_block,ny_block,max_blocks) , & ! y-component of velocity on U grid (m/s)
uvelT (nx_block,ny_block,max_blocks) , & ! x-component of velocity on T grid (m/s)
vvelT (nx_block,ny_block,max_blocks) , & ! y-component of velocity on T grid (m/s)
uvelE (nx_block,ny_block,max_blocks) , & ! x-component of velocity on E grid (m/s)
vvelE (nx_block,ny_block,max_blocks) , & ! y-component of velocity on E grid (m/s)
uvelN (nx_block,ny_block,max_blocks) , & ! x-component of velocity on N grid (m/s)
Expand Down
57 changes: 35 additions & 22 deletions cicecore/cicedynB/general/ice_step_mod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,10 @@ module ice_step_mod
step_snow, prep_radiation, step_radiation, ocean_mixed_layer, &
update_state, biogeochemistry, step_dyn_wave, step_prep

real (kind=dbl_kind), dimension (:,:,:), allocatable :: &
uvelT_icep, & ! uvel for wind stress computation in icepack
vvelT_icep ! vvel for wind stress computation in icepack

!=======================================================================

contains
Expand Down Expand Up @@ -80,6 +84,13 @@ subroutine step_prep

use ice_flux, only: uatm, vatm, uatmT, vatmT
use ice_grid, only: grid_atm_dynu, grid_atm_dynv, grid_average_X2Y
use ice_state, only: uvel, vvel

logical (kind=log_kind) :: &
highfreq ! highfreq flag

logical (kind=log_kind), save :: &
first_call = .true. ! first call flag

character(len=*), parameter :: subname = '(step_prep)'

Expand All @@ -92,6 +103,26 @@ subroutine step_prep
call grid_average_X2Y('S',uatm,grid_atm_dynu,uatmT,'T')
call grid_average_X2Y('S',vatm,grid_atm_dynv,vatmT,'T')

!-----------------------------------------------------------------
! Compute uvelT_icep, vvelT_icep
!-----------------------------------------------------------------

if (first_call) then
allocate(uvelT_icep(nx_block,ny_block,max_blocks))
allocate(vvelT_icep(nx_block,ny_block,max_blocks))
uvelT_icep = c0
vvelT_icep = c0
endif

call icepack_query_parameters(highfreq_out=highfreq)

if (highfreq) then
call grid_average_X2Y('A', uvel, 'U', uvelT_icep, 'T')
call grid_average_X2Y('A', vvel, 'U', vvelT_icep, 'T')
endif

first_call = .false.

end subroutine step_prep

!=======================================================================
Expand Down Expand Up @@ -209,7 +240,7 @@ subroutine step_therm1 (dt, iblk)
Qa_iso, Qref_iso, fiso_evap, HDO_ocn, H2_16O_ocn, H2_18O_ocn
use ice_grid, only: lmask_n, lmask_s, tmask
use ice_state, only: aice, aicen, aicen_init, vicen_init, &
vice, vicen, vsno, vsnon, trcrn, uvelT, vvelT, vsnon_init
vice, vicen, vsno, vsnon, trcrn, vsnon_init
#ifdef CICE_IN_NEMO
use ice_state, only: aice_init
#endif
Expand Down Expand Up @@ -251,8 +282,6 @@ subroutine step_therm1 (dt, iblk)
tr_pond_lvl, tr_pond_topo, calc_Tsfc, highfreq, tr_snow

real (kind=dbl_kind) :: &
uvelTij, & ! cell-centered velocity, x component (m/s)
vvelTij, & ! cell-centered velocity, y component (m/s)
puny ! a very small number

real (kind=dbl_kind), dimension(n_aero,2,ncat) :: &
Expand Down Expand Up @@ -336,14 +365,6 @@ subroutine step_therm1 (dt, iblk)
do j = jlo, jhi
do i = ilo, ihi

if (highfreq) then ! include ice velocity in calculation of wind stress
uvelTij = uvelT(i,j,iblk)
vvelTij = vvelT(i,j,iblk)
else
uvelTij = c0
vvelTij = c0
endif ! highfreq

if (tr_snow) then
do n = 1, ncat
do k = 1, nslyr
Expand Down Expand Up @@ -389,8 +410,8 @@ subroutine step_therm1 (dt, iblk)
vicen = vicen (i,j,:,iblk), &
vsno = vsno (i,j, iblk), &
vsnon = vsnon (i,j,:,iblk), &
uvel = uvelTij , &
vvel = vvelTij , &
uvel = uvelT_icep (i,j, iblk), &
vvel = vvelT_icep (i,j, iblk), &
Tsfc = trcrn (i,j,nt_Tsfc,:,iblk), &
zqsn = trcrn (i,j,nt_qsno:nt_qsno+nslyr-1,:,iblk), &
zqin = trcrn (i,j,nt_qice:nt_qice+nilyr-1,:,iblk), &
Expand Down Expand Up @@ -939,7 +960,7 @@ subroutine step_dyn_horiz (dt)
use ice_flux, only: strocnxU, strocnyU, strocnxT_iavg, strocnyT_iavg
use ice_flux, only: init_history_dyn
use ice_grid, only: grid_average_X2Y
use ice_state, only: aiU, uvel, vvel, uvelT, vvelT
use ice_state, only: aiU
use ice_transport_driver, only: advection, transport_upwind, transport_remap

real (kind=dbl_kind), intent(in) :: &
Expand Down Expand Up @@ -971,14 +992,6 @@ subroutine step_dyn_horiz (dt)
if (kdyn == 2) call eap (dt)
if (kdyn == 3) call implicit_solver (dt)

!-----------------------------------------------------------------
! Compute uvelT, vvelT
! only needed for highfreq, but compute anyway
!-----------------------------------------------------------------

call grid_average_X2Y('A', uvel, 'U', uvelT, 'T')
call grid_average_X2Y('A', vvel, 'U', vvelT, 'T')

!-----------------------------------------------------------------
! Compute strocnxT_iavg, strocnyT_iavg for thermo and coupling
!-----------------------------------------------------------------
Expand Down
Loading

0 comments on commit aa1e066

Please sign in to comment.