Skip to content

Commit

Permalink
Update mct/cesm1 driver, fix ss_tlt implementation (CICE-Consortium#726)
Browse files Browse the repository at this point in the history 
* Update ss_tltx, ss_tlty averaging/usage in evp, eap, vp
Update mct/cesm1 driver for snow mods and other recent changes
Remove USE_NETCDF cpp in io_pio2/ice_restart.F90, not needed or correct

* Add tripoleOnly flag to ice_HaloUpdate2DR8.  This is needed in very
limited situations where we want to ensure bit-for-bit across the seam,
but do NOT want to update the halo in general.  In particular, this is
needed for UAREA during initialization.  UAREA is explicitly computed
in the halo to avoid missing data.  In some cases, uarea computed
independently along the seam in the halo points results in roundoff
level differences.  This will produce a small error in the simulation
but more important, breaks the symmetry across the seam.  In that case,
we just want to sync up the halo, but we do NOT want to update the
independent UAREA halo computations elsewhere because of land block
elimination and other situations where blocks have no neighbors.
  • Loading branch information
@apcraig
apcraig authored Jun 3, 2022
1 parent c334aee commit 7705e13
Show file tree
Hide file tree
Showing 10 changed files with 337 additions and 166 deletions.
10 changes: 7 additions & 3 deletions cicecore/cicedynB/dynamics/ice_dyn_eap.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -175,6 +175,8 @@ subroutine eap (dt)
real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
uocnU , & ! i ocean current (m/s)
vocnU , & ! j ocean current (m/s)
ss_tltxU , & ! sea surface slope, x-direction (m/m)
ss_tltyU , & ! sea surface slope, y-direction (m/m)
tmass , & ! total mass of ice and snow (kg/m^2)
waterx , & ! for ocean stress calculation, x (m/s)
watery , & ! for ocean stress calculation, y (m/s)
Expand Down Expand Up @@ -270,8 +272,10 @@ subroutine eap (dt)

call grid_average_X2Y('F', tmass , 'T' , umass, 'U')
call grid_average_X2Y('F', 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', 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 @@ -318,7 +322,7 @@ subroutine eap (dt)
umask (:,:,iblk), &
uocnU (:,:,iblk), vocnU (:,:,iblk), &
strairx (:,:,iblk), strairy (:,:,iblk), &
ss_tltx (:,:,iblk), ss_tlty (:,:,iblk), &
ss_tltxU (:,:,iblk), ss_tltyU (:,:,iblk), &
icetmask (:,:,iblk), iceumask (:,:,iblk), &
fm (:,:,iblk), dt, &
strtltx (:,:,iblk), strtlty (:,:,iblk), &
Expand Down
20 changes: 12 additions & 8 deletions cicecore/cicedynB/dynamics/ice_dyn_evp.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -324,14 +324,18 @@ subroutine evp (dt)
call grid_average_X2Y('S', ss_tlty , grid_ocn_dynv, ss_tltyU, 'U')

if (grid_ice == 'CD' .or. grid_ice == 'C') then
call grid_average_X2Y('F', tmass , 'T' , emass, 'E')
call grid_average_X2Y('F', aice_init, 'T' , aie , '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('F', tmass , 'T' , nmass, 'N')
call grid_average_X2Y('F', aice_init, 'T' , ain , '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('F', tmass , 'T' , emass , 'E')
call grid_average_X2Y('F', aice_init, 'T' , aie , '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('F', tmass , 'T' , nmass , 'N')
call grid_average_X2Y('F', aice_init, 'T' , ain , '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')
call grid_average_X2Y('S', ss_tlty , grid_ocn_dynv, ss_tltyN, 'N')
endif
!----------------------------------------------------------------
! Set wind stress to values supplied via NEMO or other forcing
Expand Down
14 changes: 9 additions & 5 deletions cicecore/cicedynB/dynamics/ice_dyn_vp.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -196,6 +196,8 @@ subroutine implicit_solver (dt)
real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
uocnU , & ! i ocean current (m/s)
vocnU , & ! j ocean current (m/s)
ss_tltxU , & ! sea surface slope, x-direction (m/m)
ss_tltyU , & ! sea surface slope, y-direction (m/m)
tmass , & ! total mass of ice and snow (kg/m^2)
waterx , & ! for ocean stress calculation, x (m/s)
watery , & ! for ocean stress calculation, y (m/s)
Expand Down Expand Up @@ -300,10 +302,12 @@ subroutine implicit_solver (dt)
! convert fields from T to U grid
!-----------------------------------------------------------------

call grid_average_X2Y('F',tmass,'T',umass,'U')
call grid_average_X2Y('F',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('F',tmass , 'T', umass, 'U')
call grid_average_X2Y('F',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')

!----------------------------------------------------------------
! Set wind stress to values supplied via NEMO or other forcing
Expand Down Expand Up @@ -351,7 +355,7 @@ subroutine implicit_solver (dt)
umask (:,:,iblk), &
uocnU (:,:,iblk), vocnU (:,:,iblk), &
strairx (:,:,iblk), strairy (:,:,iblk), &
ss_tltx (:,:,iblk), ss_tlty (:,:,iblk), &
ss_tltxU (:,:,iblk), ss_tltyU (:,:,iblk), &
icetmask (:,:,iblk), iceumask (:,:,iblk), &
fm (:,:,iblk), dt, &
strtltx (:,:,iblk), strtlty (:,:,iblk), &
Expand Down
87 changes: 61 additions & 26 deletions cicecore/cicedynB/infrastructure/comm/mpi/ice_boundary.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1159,7 +1159,7 @@ end subroutine ice_HaloMask

subroutine ice_HaloUpdate2DR8(array, halo, &
fieldLoc, fieldKind, &
fillValue)
fillValue, tripoleOnly)

! This routine updates ghost cells for an input array and is a
! member of a group of routines under the generic interface
Expand All @@ -1181,6 +1181,12 @@ subroutine ice_HaloUpdate2DR8(array, halo, &
! (e.g. eliminated land blocks or
! closed boundaries)

logical (log_kind), intent(in), optional :: &
tripoleOnly ! optional flag to execute halo only across tripole seam.
! this is required for a few fields where we just want to
! ensure the tripole seam is synced up to preserve symmetry.
! Added June, 2022 by tcraig. Only added to 2DR8 for now.

real (dbl_kind), dimension(:,:,:), intent(inout) :: &
array ! array containing field for which halo
! needs to be updated
Expand Down Expand Up @@ -1215,6 +1221,9 @@ subroutine ice_HaloUpdate2DR8(array, halo, &
fill, &! value to use for unknown points
x1,x2,xavg ! scalars for enforcing symmetry at U pts

logical (log_kind) :: &
ltripoleOnly ! local tripoleOnly value

integer (int_kind) :: len ! length of messages

character(len=*), parameter :: subname = '(ice_HaloUpdate2DR8)'
Expand All @@ -1231,6 +1240,12 @@ subroutine ice_HaloUpdate2DR8(array, halo, &
fill = 0.0_dbl_kind
endif

if (present(tripoleOnly)) then
ltripoleOnly = tripoleOnly
else
ltripoleOnly = .false.
endif

nxGlobal = 0
if (allocated(bufTripoleR8)) then
nxGlobal = size(bufTripoleR8,dim=1)
Expand Down Expand Up @@ -1302,19 +1317,24 @@ subroutine ice_HaloUpdate2DR8(array, halo, &
!
!-----------------------------------------------------------------------

do iblk = 1, halo%numLocalBlocks
call get_block_parameter(halo%blockGlobalID(iblk), &
ilo=ilo, ihi=ihi, &
jlo=jlo, jhi=jhi)
do j = 1,nghost
array(1:nx_block, jlo-j,iblk) = fill
array(1:nx_block, jhi+j,iblk) = fill
enddo
do i = 1,nghost
array(ilo-i, 1:ny_block,iblk) = fill
array(ihi+i, 1:ny_block,iblk) = fill
if (ltripoleOnly) then
! skip fill, not needed since tripole seam always exists if running
! on tripole grid and set tripoleOnly flag
else
do iblk = 1, halo%numLocalBlocks
call get_block_parameter(halo%blockGlobalID(iblk), &
ilo=ilo, ihi=ihi, &
jlo=jlo, jhi=jhi)
do j = 1,nghost
array(1:nx_block, jlo-j,iblk) = fill
array(1:nx_block, jhi+j,iblk) = fill
enddo
do i = 1,nghost
array(ilo-i, 1:ny_block,iblk) = fill
array(ihi+i, 1:ny_block,iblk) = fill
enddo
enddo
enddo
endif

!-----------------------------------------------------------------------
!
Expand All @@ -1334,16 +1354,25 @@ subroutine ice_HaloUpdate2DR8(array, halo, &
jDst = halo%dstLocalAddr(2,nmsg)
dstBlock = halo%dstLocalAddr(3,nmsg)

if (srcBlock > 0) then
if (dstBlock > 0) then
array(iDst,jDst,dstBlock) = &
array(iSrc,jSrc,srcBlock)
else if (dstBlock < 0) then ! tripole copy into buffer
bufTripoleR8(iDst,jDst) = &
array(iSrc,jSrc,srcBlock)
if (ltripoleOnly) then
if (srcBlock > 0) then
if (dstBlock < 0) then ! tripole copy into buffer
bufTripoleR8(iDst,jDst) = &
array(iSrc,jSrc,srcBlock)
endif
endif
else
if (srcBlock > 0) then
if (dstBlock > 0) then
array(iDst,jDst,dstBlock) = &
array(iSrc,jSrc,srcBlock)
else if (dstBlock < 0) then ! tripole copy into buffer
bufTripoleR8(iDst,jDst) = &
array(iSrc,jSrc,srcBlock)
endif
else if (srcBlock == 0) then
array(iDst,jDst,dstBlock) = fill
endif
else if (srcBlock == 0) then
array(iDst,jDst,dstBlock) = fill
endif
end do

Expand All @@ -1362,10 +1391,16 @@ subroutine ice_HaloUpdate2DR8(array, halo, &
jDst = halo%recvAddr(2,n,nmsg)
dstBlock = halo%recvAddr(3,n,nmsg)

if (dstBlock > 0) then
array(iDst,jDst,dstBlock) = bufRecvR8(n,nmsg)
else if (dstBlock < 0) then !tripole
bufTripoleR8(iDst,jDst) = bufRecvR8(n,nmsg)
if (ltripoleOnly) then
if (dstBlock < 0) then !tripole
bufTripoleR8(iDst,jDst) = bufRecvR8(n,nmsg)
endif
else
if (dstBlock > 0) then
array(iDst,jDst,dstBlock) = bufRecvR8(n,nmsg)
else if (dstBlock < 0) then !tripole
bufTripoleR8(iDst,jDst) = bufRecvR8(n,nmsg)
endif
endif
end do
end do
Expand Down
70 changes: 48 additions & 22 deletions cicecore/cicedynB/infrastructure/comm/serial/ice_boundary.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -643,7 +643,7 @@ end subroutine ice_HaloMask

subroutine ice_HaloUpdate2DR8(array, halo, &
fieldLoc, fieldKind, &
fillValue)
fillValue, tripoleOnly)

! This routine updates ghost cells for an input array and is a
! member of a group of routines under the generic interface
Expand All @@ -665,6 +665,9 @@ subroutine ice_HaloUpdate2DR8(array, halo, &
! (e.g. eliminated land blocks or
! closed boundaries)

logical (log_kind), intent(in), optional :: &
tripoleOnly ! optional flag to execute halo only across tripole seam

real (dbl_kind), dimension(:,:,:), intent(inout) :: &
array ! array containing field for which halo
! needs to be updated
Expand All @@ -690,6 +693,9 @@ subroutine ice_HaloUpdate2DR8(array, halo, &
fill, &! value to use for unknown points
x1,x2,xavg ! scalars for enforcing symmetry at U pts

logical (log_kind) :: &
ltripoleOnly ! local tripoleOnly value

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

!-----------------------------------------------------------------------
Expand All @@ -704,6 +710,12 @@ subroutine ice_HaloUpdate2DR8(array, halo, &
fill = 0.0_dbl_kind
endif

if (present(tripoleOnly)) then
ltripoleOnly = tripoleOnly
else
ltripoleOnly = .false.
endif

nxGlobal = 0
if (allocated(bufTripoleR8)) then
nxGlobal = size(bufTripoleR8,dim=1)
Expand All @@ -718,19 +730,24 @@ subroutine ice_HaloUpdate2DR8(array, halo, &
!
!-----------------------------------------------------------------------

do iblk = 1, halo%numLocalBlocks
call get_block_parameter(halo%blockGlobalID(iblk), &
ilo=ilo, ihi=ihi, &
jlo=jlo, jhi=jhi)
do j = 1,nghost
array(1:nx_block, jlo-j,iblk) = fill
array(1:nx_block, jhi+j,iblk) = fill
enddo
do i = 1,nghost
array(ilo-i, 1:ny_block,iblk) = fill
array(ihi+i, 1:ny_block,iblk) = fill
if (ltripoleOnly) then
! skip fill, not needed since tripole seam always exists if running
! on tripole grid and set tripoleOnly flag
else
do iblk = 1, halo%numLocalBlocks
call get_block_parameter(halo%blockGlobalID(iblk), &
ilo=ilo, ihi=ihi, &
jlo=jlo, jhi=jhi)
do j = 1,nghost
array(1:nx_block, jlo-j,iblk) = fill
array(1:nx_block, jhi+j,iblk) = fill
enddo
do i = 1,nghost
array(ilo-i, 1:ny_block,iblk) = fill
array(ihi+i, 1:ny_block,iblk) = fill
enddo
enddo
enddo
endif

!-----------------------------------------------------------------------
!
Expand All @@ -750,16 +767,25 @@ subroutine ice_HaloUpdate2DR8(array, halo, &
jDst = halo%dstLocalAddr(2,nmsg)
dstBlock = halo%dstLocalAddr(3,nmsg)

if (srcBlock > 0) then
if (dstBlock > 0) then
array(iDst,jDst,dstBlock) = &
array(iSrc,jSrc,srcBlock)
else if (dstBlock < 0) then ! tripole copy into buffer
bufTripoleR8(iDst,jDst) = &
array(iSrc,jSrc,srcBlock)
if (ltripoleOnly) then
if (srcBlock > 0) then
if (dstBlock < 0) then ! tripole copy into buffer
bufTripoleR8(iDst,jDst) = &
array(iSrc,jSrc,srcBlock)
endif
endif
else
if (srcBlock > 0) then
if (dstBlock > 0) then
array(iDst,jDst,dstBlock) = &
array(iSrc,jSrc,srcBlock)
else if (dstBlock < 0) then ! tripole copy into buffer
bufTripoleR8(iDst,jDst) = &
array(iSrc,jSrc,srcBlock)
endif
else if (srcBlock == 0) then
array(iDst,jDst,dstBlock) = fill
endif
else if (srcBlock == 0) then
array(iDst,jDst,dstBlock) = fill
endif
end do

Expand Down
16 changes: 16 additions & 0 deletions cicecore/cicedynB/infrastructure/ice_grid.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -604,12 +604,28 @@ subroutine init_grid2
!-----------------------------------------------------------------

call ice_timer_start(timer_bound)

call ice_HaloUpdate (dxhy, halo_info, &
field_loc_center, field_type_vector, &
fillValue=c1)
call ice_HaloUpdate (dyhx, halo_info, &
field_loc_center, field_type_vector, &
fillValue=c1)

! Update just on the tripole seam to ensure bit-for-bit symmetry across seam
call ice_HaloUpdate (tarea, halo_info, &
field_loc_center, field_type_scalar, &
fillValue=c1, tripoleOnly=.true.)
call ice_HaloUpdate (uarea, halo_info, &
field_loc_NEcorner, field_type_scalar, &
fillValue=c1, tripoleOnly=.true.)
call ice_HaloUpdate (tarear, halo_info, &
field_loc_center, field_type_scalar, &
fillValue=c1, tripoleOnly=.true.)
call ice_HaloUpdate (uarear, halo_info, &
field_loc_NEcorner, field_type_scalar, &
fillValue=c1, tripoleOnly=.true.)

call ice_timer_stop(timer_bound)

!-----------------------------------------------------------------
Expand Down
6 changes: 1 addition & 5 deletions cicecore/cicedynB/infrastructure/io/io_pio2/ice_restart.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -946,16 +946,12 @@ logical function query_field(nu,vname)
character(len=*), parameter :: subname = '(query_field)'

query_field = .false.
#ifdef USE_NETCDF

if (my_task == master_task) then
status = pio_inq_varid(File,trim(vname),vardesc)
if (status == PIO_noerr) query_field = .true.
endif
call broadcast_scalar(query_field,master_task)
#else
call abort_ice(subname//'ERROR: USE_NETCDF cpp not defined for '//trim(ice_ic), &
file=__FILE__, line=__LINE__)
#endif

end function query_field

Expand Down
Loading

0 comments on commit 7705e13

Please sign in to comment.