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

Jump to bottom

Update KE diagnostic #784

Merged
merged 2 commits into from
Nov 7, 2022
Merged

Update KE diagnostic #784

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
096320c
Update KE diagnostic
apcraig Nov 4, 2022
db254c1
Add grid average NE2TA, NE2UA
apcraig Nov 5, 2022
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
156 changes: 74 additions & 82 deletions cicecore/cicedynB/analysis/ice_diagnostics.F90
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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