Skip to content

Commit

Permalink
Update computation of cdn_ocn for use in dynamics (CICE-Consortium#771)
Browse files Browse the repository at this point in the history 
- Compute and used cdn_ocn on U, E, and N cell locations as needed for dynamics.
- Add halo updates in dynamics before calling grid_average.  This doesn't
  change answers, but is the safe thing to do in general.  Performance does
  not seem to be affected.
  • Loading branch information
@apcraig
apcraig authored Oct 14, 2022
1 parent 6a62a11 commit 0447b9e
Show file tree
Hide file tree
Showing 4 changed files with 80 additions and 29 deletions.
27 changes: 21 additions & 6 deletions cicecore/cicedynB/dynamics/ice_dyn_eap.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -177,6 +177,7 @@ subroutine eap (dt)
vocnU , & ! j ocean current (m/s)
ss_tltxU , & ! sea surface slope, x-direction (m/m)
ss_tltyU , & ! sea surface slope, y-direction (m/m)
cdn_ocnU , & ! ocn drag coefficient
tmass , & ! total mass of ice and snow (kg/m^2)
waterxU , & ! for ocean stress calculation, x (m/s)
wateryU , & ! for ocean stress calculation, y (m/s)
Expand All @@ -186,7 +187,9 @@ subroutine eap (dt)
umassdti ! mass of U-cell/dte (kg/m^2 s)

real (kind=dbl_kind), allocatable :: &
fld2(:,:,:,:) ! temporary for stacking fields for halo update
fld2(:,:,:,:), & ! temporary for stacking fields for halo update
fld3(:,:,:,:), & ! temporary for stacking fields for halo update
fld4(:,:,:,:) ! temporary for stacking fields for halo update

real (kind=dbl_kind), dimension(nx_block,ny_block,8):: &
strtmp ! stress combinations for momentum equation
Expand All @@ -213,6 +216,8 @@ subroutine eap (dt)
!-----------------------------------------------------------------

allocate(fld2(nx_block,ny_block,2,max_blocks))
allocate(fld3(nx_block,ny_block,3,max_blocks))
allocate(fld4(nx_block,ny_block,4,max_blocks))

! This call is needed only if dt changes during runtime.
! call set_evp_parameters (dt)
Expand Down Expand Up @@ -265,8 +270,18 @@ subroutine eap (dt)
! convert fields from T to U grid
!-----------------------------------------------------------------

call grid_average_X2Y('S', tmass , 'T' , umass, 'U')
call grid_average_X2Y('S', aice_init, 'T' , aiU , 'U')
call stack_fields(tmass, aice_init, cdn_ocn, fld3)
call ice_HaloUpdate (fld3, halo_info, &
field_loc_center, field_type_scalar)
call stack_fields(uocn, vocn, ss_tltx, ss_tlty, fld4)
call ice_HaloUpdate (fld4, halo_info, &
field_loc_center, field_type_vector)
call unstack_fields(fld3, tmass, aice_init, cdn_ocn)
call unstack_fields(fld4, uocn, vocn, ss_tltx, ss_tlty)

call grid_average_X2Y('S', tmass , 'T' , umass , 'U')
call grid_average_X2Y('S', aice_init, 'T' , aiU , 'U')
call grid_average_X2Y('S', cdn_ocn , 'T' , cdn_ocnU, 'U')
call grid_average_X2Y('S', uocn , grid_ocn_dynu, uocnU , 'U')
call grid_average_X2Y('S', vocn , grid_ocn_dynv, vocnU , 'U')
call grid_average_X2Y('S', ss_tltx , grid_ocn_dynu, ss_tltxU, 'U')
Expand Down Expand Up @@ -486,7 +501,7 @@ subroutine eap (dt)

! call ice_timer_start(timer_tmp2,iblk)
call stepu (nx_block, ny_block, &
icellu (iblk), Cdn_ocn (:,:,iblk), &
icellu (iblk), Cdn_ocnU (:,:,iblk), &
indxui (:,iblk), indxuj (:,iblk), &
aiU (:,:,iblk), strtmp (:,:,:), &
uocnU (:,:,iblk), vocnU (:,:,iblk), &
Expand Down Expand Up @@ -540,7 +555,7 @@ subroutine eap (dt)

enddo ! subcycling

deallocate(fld2)
deallocate(fld2,fld3,fld4)
if (maskhalo_dyn) call ice_HaloDestroy(halo_info_mask)

!-----------------------------------------------------------------
Expand All @@ -552,7 +567,7 @@ subroutine eap (dt)

call dyn_finish &
(nx_block, ny_block, &
icellu (iblk), Cdn_ocn (:,:,iblk), &
icellu (iblk), Cdn_ocnU(:,:,iblk), &
indxui (:,iblk), indxuj (:,iblk), &
uvel (:,:,iblk), vvel (:,:,iblk), &
uocnU (:,:,iblk), vocnU (:,:,iblk), &
Expand Down
44 changes: 32 additions & 12 deletions cicecore/cicedynB/dynamics/ice_dyn_evp.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -150,6 +150,7 @@ subroutine evp (dt)
vocnU , & ! j ocean current (m/s)
ss_tltxU , & ! sea surface slope, x-direction (m/m)
ss_tltyU , & ! sea surface slope, y-direction (m/m)
cdn_ocnU , & ! ocn drag coefficient
tmass , & ! total mass of ice and snow (kg/m^2)
waterxU , & ! for ocean stress calculation, x (m/s)
wateryU , & ! for ocean stress calculation, y (m/s)
Expand All @@ -163,6 +164,7 @@ subroutine evp (dt)
vocnN , & ! j ocean current (m/s)
ss_tltxN , & ! sea surface slope, x-direction (m/m)
ss_tltyN , & ! sea surface slope, y-direction (m/m)
cdn_ocnN , & ! ocn drag coefficient
waterxN , & ! for ocean stress calculation, x (m/s)
wateryN , & ! for ocean stress calculation, y (m/s)
forcexN , & ! work array: combined atm stress and ocn tilt, x
Expand All @@ -176,6 +178,7 @@ subroutine evp (dt)
vocnE , & ! j ocean current (m/s)
ss_tltxE , & ! sea surface slope, x-direction (m/m)
ss_tltyE , & ! sea surface slope, y-direction (m/m)
cdn_ocnE , & ! ocn drag coefficient
waterxE , & ! for ocean stress calculation, x (m/s)
wateryE , & ! for ocean stress calculation, y (m/s)
forcexE , & ! work array: combined atm stress and ocn tilt, x
Expand All @@ -185,9 +188,9 @@ subroutine evp (dt)
emassdti ! mass of E-cell/dte (kg/m^2 s)

real (kind=dbl_kind), allocatable :: &
fld2(:,:,:,:) , & ! 2 bundled fields
fld3(:,:,:,:) , & ! 3 bundled fields
fld4(:,:,:,:) ! 4 bundled fields
fld2(:,:,:,:), & ! 2 bundled fields
fld3(:,:,:,:), & ! 3 bundled fields
fld4(:,:,:,:) ! 4 bundled fields

real (kind=dbl_kind), allocatable :: &
strengthU(:,:,:), & ! strength averaged to U points
Expand Down Expand Up @@ -312,8 +315,18 @@ subroutine evp (dt)
! convert fields from T to U grid
!-----------------------------------------------------------------

call stack_fields(tmass, aice_init, cdn_ocn, fld3)
call ice_HaloUpdate (fld3, halo_info, &
field_loc_center, field_type_scalar)
call stack_fields(uocn, vocn, ss_tltx, ss_tlty, fld4)
call ice_HaloUpdate (fld4, halo_info, &
field_loc_center, field_type_vector)
call unstack_fields(fld3, tmass, aice_init, cdn_ocn)
call unstack_fields(fld4, uocn, vocn, ss_tltx, ss_tlty)

call grid_average_X2Y('S', tmass , 'T' , umass , 'U')
call grid_average_X2Y('S', aice_init, 'T' , aiU , 'U')
call grid_average_X2Y('S', cdn_ocn , 'T' , cdn_ocnU, 'U')
call grid_average_X2Y('S', uocn , grid_ocn_dynu, uocnU , 'U')
call grid_average_X2Y('S', vocn , grid_ocn_dynv, vocnU , 'U')
call grid_average_X2Y('S', ss_tltx , grid_ocn_dynu, ss_tltxU, 'U')
Expand All @@ -322,12 +335,14 @@ subroutine evp (dt)
if (grid_ice == 'CD' .or. grid_ice == 'C') then
call grid_average_X2Y('S', tmass , 'T' , emass , 'E')
call grid_average_X2Y('S', aice_init, 'T' , aie , 'E')
call grid_average_X2Y('S', cdn_ocn , 'T' , cdn_ocnE, 'E')
call grid_average_X2Y('S', uocn , grid_ocn_dynu, uocnE , 'E')
call grid_average_X2Y('S', vocn , grid_ocn_dynv, vocnE , 'E')
call grid_average_X2Y('S', ss_tltx , grid_ocn_dynu, ss_tltxE, 'E')
call grid_average_X2Y('S', ss_tlty , grid_ocn_dynv, ss_tltyE, 'E')
call grid_average_X2Y('S', tmass , 'T' , nmass , 'N')
call grid_average_X2Y('S', aice_init, 'T' , ain , 'N')
call grid_average_X2Y('S', cdn_ocn , 'T' , cdn_ocnN, 'N')
call grid_average_X2Y('S', uocn , grid_ocn_dynu, uocnN , 'N')
call grid_average_X2Y('S', vocn , grid_ocn_dynv, vocnN , 'N')
call grid_average_X2Y('S', ss_tltx , grid_ocn_dynu, ss_tltxN, 'N')
Expand Down Expand Up @@ -720,7 +735,7 @@ subroutine evp (dt)
call ice_dyn_evp_1d_copyin( &
nx_block,ny_block,nblocks,nx_global+2*nghost,ny_global+2*nghost, &
icetmask, iceumask, &
cdn_ocn,aiU,uocnU,vocnU,forcexU,forceyU,TbU, &
Cdn_ocnU,aiU,uocnU,vocnU,forcexU,forceyU,TbU, &
umassdti,fmU,uarear,tarear,strintxU,strintyU,uvel_init,vvel_init,&
strength,uvel,vvel,dxT,dyT, &
stressp_1 ,stressp_2, stressp_3, stressp_4, &
Expand Down Expand Up @@ -773,7 +788,7 @@ subroutine evp (dt)
! momentum equation
!-----------------------------------------------------------------
call stepu (nx_block , ny_block , &
icellu (iblk), Cdn_ocn (:,:,iblk), &
icellu (iblk), Cdn_ocnU(:,:,iblk), &
indxui (:,iblk), indxuj (:,iblk), &
aiU (:,:,iblk), strtmp (:,:,:), &
uocnU (:,:,iblk), vocnU (:,:,iblk), &
Expand Down Expand Up @@ -924,7 +939,7 @@ subroutine evp (dt)
do iblk = 1, nblocks

call stepu_C (nx_block , ny_block , & ! u, E point
icelle (iblk), Cdn_ocn (:,:,iblk), &
icelle (iblk), Cdn_ocnE (:,:,iblk), &
indxei (:,iblk), indxej (:,iblk), &
aiE (:,:,iblk), &
uocnE (:,:,iblk), vocnE (:,:,iblk), &
Expand All @@ -936,7 +951,7 @@ subroutine evp (dt)
TbE (:,:,iblk))

call stepv_C (nx_block, ny_block, & ! v, N point
icelln (iblk), Cdn_ocn (:,:,iblk), &
icelln (iblk), Cdn_ocnN (:,:,iblk), &
indxni (:,iblk), indxnj (:,iblk), &
aiN (:,:,iblk), &
uocnN (:,:,iblk), vocnN (:,:,iblk), &
Expand Down Expand Up @@ -1124,7 +1139,7 @@ subroutine evp (dt)
do iblk = 1, nblocks

call stepuv_CD (nx_block , ny_block , & ! E point
icelle (iblk), Cdn_ocn (:,:,iblk), &
icelle (iblk), Cdn_ocnE (:,:,iblk), &
indxei (:,iblk), indxej (:,iblk), &
aiE (:,:,iblk), &
uocnE (:,:,iblk), vocnE (:,:,iblk), &
Expand All @@ -1138,7 +1153,7 @@ subroutine evp (dt)
TbE (:,:,iblk))

call stepuv_CD (nx_block , ny_block , & ! N point
icelln (iblk), Cdn_ocn (:,:,iblk), &
icelln (iblk), Cdn_ocnN (:,:,iblk), &
indxni (:,iblk), indxnj (:,iblk), &
aiN (:,:,iblk), &
uocnN (:,:,iblk), vocnN (:,:,iblk), &
Expand Down Expand Up @@ -1284,7 +1299,7 @@ subroutine evp (dt)
do iblk = 1, nblocks
call dyn_finish &
(nx_block , ny_block , &
icellu (iblk), Cdn_ocn (:,:,iblk), &
icellu (iblk), Cdn_ocnU(:,:,iblk), &
indxui (:,iblk), indxuj (:,iblk), &
uvel (:,:,iblk), vvel (:,:,iblk), &
uocnU (:,:,iblk), vocnU (:,:,iblk), &
Expand All @@ -1300,7 +1315,7 @@ subroutine evp (dt)

call dyn_finish &
(nx_block , ny_block , &
icelln (iblk), Cdn_ocn (:,:,iblk), &
icelln (iblk), Cdn_ocnN(:,:,iblk), &
indxni (:,iblk), indxnj (:,iblk), &
uvelN (:,:,iblk), vvelN (:,:,iblk), &
uocnN (:,:,iblk), vocnN (:,:,iblk), &
Expand All @@ -1309,7 +1324,7 @@ subroutine evp (dt)

call dyn_finish &
(nx_block , ny_block , &
icelle (iblk), Cdn_ocn (:,:,iblk), &
icelle (iblk), Cdn_ocnE(:,:,iblk), &
indxei (:,iblk), indxej (:,iblk), &
uvelE (:,:,iblk), vvelE (:,:,iblk), &
uocnE (:,:,iblk), vocnE (:,:,iblk), &
Expand All @@ -1322,6 +1337,11 @@ subroutine evp (dt)
endif

if (grid_ice == 'CD' .or. grid_ice == 'C') then
call ice_HaloUpdate (strintxE, halo_info, &
field_loc_Eface, field_type_vector)
call ice_HaloUpdate (strintyN, halo_info, &
field_loc_Nface, field_type_vector)
call ice_timer_stop(timer_bound)
call grid_average_X2Y('S', strintxE, 'E', strintxU, 'U') ! diagnostic
call grid_average_X2Y('S', strintyN, 'N', strintyU, 'U') ! diagnostic
endif
Expand Down
37 changes: 26 additions & 11 deletions cicecore/cicedynB/dynamics/ice_dyn_vp.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -106,7 +106,9 @@ module ice_dyn_vp
indxuj(:,:) ! compressed index in j-direction

real (kind=dbl_kind), allocatable :: &
fld2(:,:,:,:) ! work array for boundary updates
fld2(:,:,:,:), & ! work array for boundary updates
fld3(:,:,:,:), & ! work array for boundary updates
fld4(:,:,:,:) ! work array for boundary updates

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

Expand Down Expand Up @@ -142,6 +144,8 @@ subroutine init_vp
indxui(nx_block*ny_block, max_blocks), &
indxuj(nx_block*ny_block, max_blocks))
allocate(fld2(nx_block,ny_block,2,max_blocks))
allocate(fld3(nx_block,ny_block,3,max_blocks))
allocate(fld4(nx_block,ny_block,4,max_blocks))

end subroutine init_vp

Expand Down Expand Up @@ -200,6 +204,7 @@ subroutine implicit_solver (dt)
vocnU , & ! j ocean current (m/s)
ss_tltxU , & ! sea surface slope, x-direction (m/m)
ss_tltyU , & ! sea surface slope, y-direction (m/m)
cdn_ocnU , & ! ocn drag coefficient
tmass , & ! total mass of ice and snow (kg/m^2)
waterxU , & ! for ocean stress calculation, x (m/s)
wateryU , & ! for ocean stress calculation, y (m/s)
Expand Down Expand Up @@ -299,12 +304,22 @@ subroutine implicit_solver (dt)
! convert fields from T to U grid
!-----------------------------------------------------------------

call grid_average_X2Y('S',tmass , 'T', umass, 'U')
call grid_average_X2Y('S',aice_init, 'T', aiU , 'U')
call grid_average_X2Y('S',uocn , grid_ocn_dynu, uocnU , 'U')
call grid_average_X2Y('S',vocn , grid_ocn_dynv, vocnU , 'U')
call grid_average_X2Y('S',ss_tltx, grid_ocn_dynu, ss_tltxU, 'U')
call grid_average_X2Y('S',ss_tlty, grid_ocn_dynv, ss_tltyU, 'U')
call stack_fields(tmass, aice_init, cdn_ocn, fld3)
call ice_HaloUpdate (fld3, halo_info, &
field_loc_center, field_type_scalar)
call stack_fields(uocn, vocn, ss_tltx, ss_tlty, fld4)
call ice_HaloUpdate (fld4, halo_info, &
field_loc_center, field_type_vector)
call unstack_fields(fld3, tmass, aice_init, cdn_ocn)
call unstack_fields(fld4, uocn, vocn, ss_tltx, ss_tlty)

call grid_average_X2Y('S',tmass , 'T' , umass , 'U')
call grid_average_X2Y('S',aice_init, 'T' , aiU , 'U')
call grid_average_X2Y('S',cdn_ocn , 'T' , cdn_ocnU, 'U')
call grid_average_X2Y('S',uocn , grid_ocn_dynu, uocnU , 'U')
call grid_average_X2Y('S',vocn , grid_ocn_dynv, vocnU , 'U')
call grid_average_X2Y('S',ss_tltx , grid_ocn_dynu, ss_tltxU, 'U')
call grid_average_X2Y('S',ss_tlty , grid_ocn_dynv, ss_tltyU, 'U')

!----------------------------------------------------------------
! Set wind stress to values supplied via NEMO or other forcing
Expand Down Expand Up @@ -479,7 +494,7 @@ subroutine implicit_solver (dt)
umassdti, sol , &
fpresx , fpresy , &
zetax2 , etax2 , &
rep_prs , &
rep_prs , cdn_ocnU,&
Cb, halo_info_mask)
!-----------------------------------------------------------------
! End of nonlinear iteration
Expand Down Expand Up @@ -623,7 +638,7 @@ subroutine implicit_solver (dt)

call dyn_finish &
(nx_block, ny_block, &
icellu (iblk), Cdn_ocn (:,:,iblk), &
icellu (iblk), Cdn_ocnU(:,:,iblk), &
indxui (:,iblk), indxuj (:,iblk), &
uvel (:,:,iblk), vvel (:,:,iblk), &
uocnU (:,:,iblk), vocnU (:,:,iblk), &
Expand Down Expand Up @@ -669,10 +684,9 @@ subroutine anderson_solver (icellt , icellu , &
umassdti, sol , &
fpresx , fpresy , &
zetax2 , etax2 , &
rep_prs , &
rep_prs , cdn_ocn, &
Cb, halo_info_mask)

use ice_arrays_column, only: Cdn_ocn
use ice_blocks, only: nx_block, ny_block
use ice_boundary, only: ice_HaloUpdate
use ice_constants, only: c1
Expand Down Expand Up @@ -702,6 +716,7 @@ subroutine anderson_solver (icellt , icellu , &
aiU , & ! ice fraction on u-grid
uocn , & ! i ocean current (m/s)
vocn , & ! j ocean current (m/s)
cdn_ocn , & ! ocn drag coefficient
waterxU , & ! for ocean stress calculation, x (m/s)
wateryU , & ! for ocean stress calculation, y (m/s)
bxfix , & ! part of bx that is constant during Picard
Expand Down
1 change: 1 addition & 0 deletions cicecore/shared/ice_arrays_column.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@ module ice_arrays_column
public :: alloc_arrays_column

! icepack_atmo.F90
! Cdn variables on the T-grid
real (kind=dbl_kind), public, &
dimension (:,:,:), allocatable :: &
Cdn_atm , & ! atm drag coefficient
Expand Down

0 comments on commit 0447b9e

Please sign in to comment.