From f01afd99eae4012da682d02d63ae73a28b601fc7 Mon Sep 17 00:00:00 2001 From: apcraig Date: Wed, 12 Oct 2022 10:20:13 -0600 Subject: [PATCH] Update computation of cdn_ocn for use in dynamics - 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. --- cicecore/cicedynB/dynamics/ice_dyn_eap.F90 | 27 ++++++++++--- cicecore/cicedynB/dynamics/ice_dyn_evp.F90 | 44 ++++++++++++++++------ cicecore/cicedynB/dynamics/ice_dyn_vp.F90 | 37 ++++++++++++------ cicecore/shared/ice_arrays_column.F90 | 1 + 4 files changed, 80 insertions(+), 29 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 b/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 index 854fd061b..ea254cbc0 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_eap.F90 @@ -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) @@ -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 @@ -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) @@ -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') @@ -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), & @@ -540,7 +555,7 @@ subroutine eap (dt) enddo ! subcycling - deallocate(fld2) + deallocate(fld2,fld3,fld4) if (maskhalo_dyn) call ice_HaloDestroy(halo_info_mask) !----------------------------------------------------------------- @@ -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), & diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index 3bbfc01bc..d1264bae7 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -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) @@ -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 @@ -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 @@ -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 @@ -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') @@ -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') @@ -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, & @@ -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), & @@ -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), & @@ -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), & @@ -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), & @@ -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), & @@ -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), & @@ -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), & @@ -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), & @@ -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 diff --git a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 index 631192587..2a0dbc4b3 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 @@ -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 !======================================================================= @@ -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 @@ -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) @@ -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 @@ -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 @@ -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), & @@ -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 @@ -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 diff --git a/cicecore/shared/ice_arrays_column.F90 b/cicecore/shared/ice_arrays_column.F90 index c9e8be8db..b4727d3fd 100644 --- a/cicecore/shared/ice_arrays_column.F90 +++ b/cicecore/shared/ice_arrays_column.F90 @@ -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