diff --git a/cicecore/cicedyn/analysis/ice_history.F90 b/cicecore/cicedyn/analysis/ice_history.F90 index f19158f6a..54b6ce934 100644 --- a/cicecore/cicedyn/analysis/ice_history.F90 +++ b/cicecore/cicedyn/analysis/ice_history.F90 @@ -1341,7 +1341,7 @@ subroutine init_hist (dt) "sig2 is instantaneous" // trim(description), c1, c0, & ns1, f_sig2) - call define_hist_field(n_sigP,"sigP","1",gridstr2d, gridstr, & + call define_hist_field(n_sigP,"sigP","N/m",gridstr2d, gridstr, & "ice pressure", & "sigP is instantaneous" // trim(description), c1, c0, & ns1, f_sigP) diff --git a/cicecore/cicedyn/general/ice_flux.F90 b/cicecore/cicedyn/general/ice_flux.F90 index 8d190753e..5145fec66 100644 --- a/cicecore/cicedyn/general/ice_flux.F90 +++ b/cicecore/cicedyn/general/ice_flux.F90 @@ -366,6 +366,7 @@ module ice_flux vatmT , & ! vatm on T grid (m/s) rside , & ! fraction of ice that melts laterally fside , & ! lateral heat flux (W/m^2) + wlat , & ! lateral heat rate (m/s) fsw , & ! incoming shortwave radiation (W/m^2) coszen , & ! cosine solar zenith angle, < 0 for sun below horizon rdg_conv, & ! convergence term for ridging (1/s) @@ -540,7 +541,8 @@ subroutine alloc_flux uatmT (nx_block,ny_block,max_blocks), & ! uatm on T grid vatmT (nx_block,ny_block,max_blocks), & ! vatm on T grid rside (nx_block,ny_block,max_blocks), & ! fraction of ice that melts laterally - fside (nx_block,ny_block,max_blocks), & ! lateral melt rate (W/m^2) + fside (nx_block,ny_block,max_blocks), & ! lateral melt flux (W/m^2) + wlat (nx_block,ny_block,max_blocks), & ! lateral melt rate (m/s) fsw (nx_block,ny_block,max_blocks), & ! incoming shortwave radiation (W/m^2) coszen (nx_block,ny_block,max_blocks), & ! cosine solar zenith angle, < 0 for sun below horizon rdg_conv (nx_block,ny_block,max_blocks), & ! convergence term for ridging (1/s) diff --git a/cicecore/cicedyn/general/ice_step_mod.F90 b/cicecore/cicedyn/general/ice_step_mod.F90 index 8dd6fe49a..56510c247 100644 --- a/cicecore/cicedyn/general/ice_step_mod.F90 +++ b/cicecore/cicedyn/general/ice_step_mod.F90 @@ -227,7 +227,7 @@ subroutine step_therm1 (dt, iblk) use ice_calendar, only: yday use ice_domain_size, only: ncat, nilyr, nslyr, n_iso, n_aero use ice_flux, only: frzmlt, sst, Tf, strocnxT_iavg, strocnyT_iavg, rside, fbot, Tbot, Tsnice, & - meltsn, melttn, meltbn, congeln, snoicen, uatmT, vatmT, fside, & + meltsn, melttn, meltbn, congeln, snoicen, uatmT, vatmT, fside, wlat, & wind, rhoa, potT, Qa, zlvl, zlvs, strax, stray, flatn, fsensn, fsurfn, fcondtopn, & flw, fsnow, fpond, sss, mlt_onset, frz_onset, fcondbotn, fcondbot, fsloss, & frain, Tair, strairxT, strairyT, fsurf, fcondtop, fsens, & @@ -469,6 +469,7 @@ subroutine step_therm1 (dt, iblk) frzmlt = frzmlt (i,j, iblk), & rside = rside (i,j, iblk), & fside = fside (i,j, iblk), & + wlat = wlat (i,j, iblk), & fsnow = fsnow (i,j, iblk), & frain = frain (i,j, iblk), & fpond = fpond (i,j, iblk), & @@ -618,7 +619,7 @@ subroutine step_therm2 (dt, iblk) use ice_calendar, only: yday use ice_domain_size, only: ncat, nilyr, nslyr, nblyr, nfsd use ice_flux, only: fresh, frain, fpond, frzmlt, frazil, frz_onset, & - update_ocn_f, fsalt, Tf, sss, salinz, fhocn, rside, fside, & + update_ocn_f, fsalt, Tf, sss, salinz, fhocn, rside, fside, wlat, & meltl, frazil_diag use ice_flux_bgc, only: flux_bio, faero_ocn, & fiso_ocn, HDO_ocn, H2_16O_ocn, H2_18O_ocn @@ -702,6 +703,7 @@ subroutine step_therm2 (dt, iblk) rside = rside (i,j, iblk), & meltl = meltl (i,j, iblk), & fside = fside (i,j, iblk), & + wlat = wlat (i,j, iblk), & frzmlt = frzmlt (i,j, iblk), & frazil = frazil (i,j, iblk), & frain = frain (i,j, iblk), & diff --git a/cicecore/cicedyn/infrastructure/comm/mpi/ice_boundary.F90 b/cicecore/cicedyn/infrastructure/comm/mpi/ice_boundary.F90 index 68436cd0f..2a7d68c11 100644 --- a/cicecore/cicedyn/infrastructure/comm/mpi/ice_boundary.F90 +++ b/cicecore/cicedyn/infrastructure/comm/mpi/ice_boundary.F90 @@ -12,14 +12,53 @@ module ice_boundary ! fixes for non-existent blocks ! 2008-01-28: Elizabeth Hunke replaced old routines with new POP ! infrastructure +! +!----------------------------------------------------------------------- +! +! Some notes on tripole, A-H below are gridpoints at i = 1:nx_global +! where nx_global=8. The schematics below show the general layout of the center +! points on the tripole fold. More complex pictures are needed to show +! relative orientation and offsets of east, north, and northeast points +! across the fold. See also appendix E of the NEMO_manual, +! https://zenodo.org/record/6334656#.YiYirhPMLXQ. Note the NFtype=T +! is the tripole u-fold grid with T-grid=center, U-grid=east, V-grid=north, +! and F-grid=northeast points in CICE. NFtype=F is similar to tripoleT +! except for the treatment of the poles. The CICE implementation also +! averages all degenerate points, NEMO's strategy seems to be to copy +! data from one side of the tripole to the other for degenerate points. +! +! tripole: u-fold, fold is on north edge of ny_global +! north and northeast points on the fold are degenerate and averaged +! A,H,D,and E are pole points +! +! ny_global+2 H G F E D C B A @ny_global-1 +! ny_global+1 H G F E D C B A @ny_global +! ny_global A B C D E F G H +! ny_global-1 A B C D E F G H +! +! tripoleT: t-fold, fold is thru center of ny_global +! center and east points at ny_global are degenerate and averaged +! north and northeast point at ny_global are not prognostic, they are halos +! A and E are pole points +! +! ny_global+2 H G F E D C B A @ny_global-2 +! ny_global+1 H G F E D C B A @ny_global-1 +! ny_global A BH CG DF E FD GC HB A +! ny_global-1 A B C D E F G H +! ny_global-2 A B C D E F G H +! +!----------------------------------------------------------------------- + use mpi ! MPI Fortran module use ice_kinds_mod use ice_communicate, only: my_task, mpiR4, mpiR8, mpitagHalo use ice_constants, only: field_type_scalar, & field_type_vector, field_type_angle, & + field_type_unknown, field_type_noupdate, & field_loc_center, field_loc_NEcorner, & - field_loc_Nface, field_loc_Eface + field_loc_Nface, field_loc_Eface, & + field_loc_unknown, field_loc_noupdate use ice_global_reductions, only: global_maxval use ice_exit, only: abort_ice use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted @@ -1229,6 +1268,23 @@ subroutine ice_HaloUpdate2DR8(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate2DR8)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -1552,7 +1608,7 @@ subroutine ice_HaloUpdate2DR8(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -1653,6 +1709,23 @@ subroutine ice_HaloUpdate2DR4(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate2DR4)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -1950,7 +2023,7 @@ subroutine ice_HaloUpdate2DR4(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -2051,6 +2124,23 @@ subroutine ice_HaloUpdate2DI4(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate2DI4)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -2348,7 +2438,7 @@ subroutine ice_HaloUpdate2DI4(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -2426,6 +2516,23 @@ subroutine ice_HaloUpdate2DL1(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate2DL1)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! copy logical into integer array and call haloupdate on integer array @@ -2519,6 +2626,23 @@ subroutine ice_HaloUpdate3DR8(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate3DR8)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -2858,7 +2982,7 @@ subroutine ice_HaloUpdate3DR8(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -2976,6 +3100,23 @@ subroutine ice_HaloUpdate3DR4(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate3DR4)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -3315,7 +3456,7 @@ subroutine ice_HaloUpdate3DR4(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -3433,6 +3574,23 @@ subroutine ice_HaloUpdate3DI4(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate3DI4)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -3772,7 +3930,7 @@ subroutine ice_HaloUpdate3DI4(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -3890,6 +4048,23 @@ subroutine ice_HaloUpdate4DR8(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate4DR8)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -4251,7 +4426,7 @@ subroutine ice_HaloUpdate4DR8(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -4371,6 +4546,23 @@ subroutine ice_HaloUpdate4DR4(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate4DR4)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -4732,7 +4924,7 @@ subroutine ice_HaloUpdate4DR4(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -4852,6 +5044,23 @@ subroutine ice_HaloUpdate4DI4(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate4DI4)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -5213,7 +5422,7 @@ subroutine ice_HaloUpdate4DI4(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -5265,6 +5474,7 @@ end subroutine ice_HaloUpdate4DI4 !*********************************************************************** ! This routine updates ghost cells for an input array using ! a second array as needed by the stress fields. +! This is just like 2DR8 except no averaging and only on tripole subroutine ice_HaloUpdate_stress(array1, array2, halo, & fieldLoc, fieldKind, & @@ -5319,6 +5529,23 @@ subroutine ice_HaloUpdate_stress(array1, array2, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate_stress)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -5485,30 +5712,61 @@ subroutine ice_HaloUpdate_stress(array1, array2, halo, & call abort_ice(subname//'ERROR: Unknown field kind') end select - select case (fieldLoc) - case (field_loc_center) ! cell center location + if (halo%tripoleTFlag) then + + select case (fieldLoc) + case (field_loc_center) ! cell center location + + ioffset = -1 + joffset = 0 - ioffset = 0 - joffset = 0 + case (field_loc_NEcorner) ! cell corner location - case (field_loc_NEcorner) ! cell corner location + ioffset = 0 + joffset = 1 - ioffset = 1 - joffset = 1 + case (field_loc_Eface) ! cell center location - case (field_loc_Eface) + ioffset = 0 + joffset = 0 - ioffset = 1 - joffset = 0 + case (field_loc_Nface) ! cell corner (velocity) location - case (field_loc_Nface) + ioffset = -1 + joffset = 1 - ioffset = 0 - joffset = 1 + case default + call abort_ice(subname//'ERROR: Unknown field location') + end select - case default - call abort_ice(subname//'ERROR: Unknown field location') - end select + else ! tripole u-fold + + select case (fieldLoc) + case (field_loc_center) ! cell center location + + ioffset = 0 + joffset = 0 + + case (field_loc_NEcorner) ! cell corner location + + ioffset = 1 + joffset = 1 + + case (field_loc_Eface) + + ioffset = 1 + joffset = 0 + + case (field_loc_Nface) + + ioffset = 0 + joffset = 1 + + case default + call abort_ice(subname//'ERROR: Unknown field location') + end select + + endif !*** copy out of global tripole buffer into local !*** ghost cells @@ -5531,14 +5789,15 @@ subroutine ice_HaloUpdate_stress(array1, array2, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal + if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface, do not need to replace !*** top row of physical domain, so jSrc should be !*** out of range and skipped !*** otherwise do the copy - if (jSrc <= nghost+1 .AND. jDst /= -1 ) then + if (jSrc <= halo%tripoleRows .and. jSrc>0 .and. jDst>0) then array1(iDst,jDst,dstBlock) = isign*bufTripoleR8(iSrc,jSrc) endif diff --git a/cicecore/cicedyn/infrastructure/comm/serial/ice_boundary.F90 b/cicecore/cicedyn/infrastructure/comm/serial/ice_boundary.F90 index 2b81c4441..aaebcfaad 100644 --- a/cicecore/cicedyn/infrastructure/comm/serial/ice_boundary.F90 +++ b/cicecore/cicedyn/infrastructure/comm/serial/ice_boundary.F90 @@ -12,13 +12,20 @@ module ice_boundary ! fixes for non-existent blocks ! 2008-01-28: Elizabeth Hunke replaced old routines with new POP ! infrastructure +! 2023-03-09: Tony Craig updated the implementation to fix bug in +! tripoleT and reduce number of copies in tripole overall. +! Because all blocks are local, can fill the tripole +! buffer from "north" copies. This is not true for +! the MPI version. use ice_kinds_mod use ice_communicate, only: my_task use ice_constants, only: field_type_scalar, & field_type_vector, field_type_angle, & + field_type_unknown, field_type_noupdate, & field_loc_center, field_loc_NEcorner, & - field_loc_Nface, field_loc_Eface + field_loc_Nface, field_loc_Eface, & + field_loc_unknown, field_loc_noupdate use ice_global_reductions, only: global_maxval use ice_exit, only: abort_ice use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted @@ -310,10 +317,11 @@ function ice_HaloCreate(dist, nsBoundaryType, ewBoundaryType, & !*** for that !echmod if (tripoleBlock .and. dstProc /= srcProc) then - if (tripoleBlock) then - call ice_HaloIncrementMsgCount(sendCount, recvCount, & - srcProc, dstProc, northMsgSize) - endif +! tcx,tcraig, 3/2023, this is not needed +! if (tripoleBlock) then +! call ice_HaloIncrementMsgCount(sendCount, recvCount, & +! srcProc, dstProc, northMsgSize) +! endif !*** find west neighbor block and add to message count @@ -336,10 +344,11 @@ function ice_HaloCreate(dist, nsBoundaryType, ewBoundaryType, & !*** for that !echmod if (tripoleBlock .and. dstProc /= srcProc) then - if (tripoleBlock) then - call ice_HaloIncrementMsgCount(sendCount, recvCount, & - srcProc, dstProc, northMsgSize) - endif +! tcx,tcraig, 3/2023, this is not needed +! if (tripoleBlock) then +! call ice_HaloIncrementMsgCount(sendCount, recvCount, & +! srcProc, dstProc, northMsgSize) +! endif !*** find northeast neighbor block and add to message count @@ -352,11 +361,12 @@ function ice_HaloCreate(dist, nsBoundaryType, ewBoundaryType, & call ice_distributionGetBlockLoc(dist, neBlock, dstProc, & dstLocalID) - else if (neBlock < 0) then ! tripole north row - msgSize = northMsgSize ! tripole needs whole top row of block - - call ice_distributionGetBlockLoc(dist, abs(neBlock), dstProc, & - dstLocalID) +! tcx,tcraig, 3/2023, this is not needed +! else if (neBlock < 0) then ! tripole north row +! msgSize = northMsgSize ! tripole needs whole top row of block +! +! call ice_distributionGetBlockLoc(dist, abs(neBlock), dstProc, & +! dstLocalID) else dstProc = 0 dstLocalID = 0 @@ -376,11 +386,12 @@ function ice_HaloCreate(dist, nsBoundaryType, ewBoundaryType, & call ice_distributionGetBlockLoc(dist, nwBlock, dstProc, & dstLocalID) - else if (nwBlock < 0) then ! tripole north row, count block - msgSize = northMsgSize ! tripole NE corner update - entire row needed - - call ice_distributionGetBlockLoc(dist, abs(nwBlock), dstProc, & - dstLocalID) +! tcx,tcraig, 3/2023, this is not needed +! else if (nwBlock < 0) then ! tripole north row, count block +! msgSize = northMsgSize ! tripole NE corner update - entire row needed +! +! call ice_distributionGetBlockLoc(dist, abs(nwBlock), dstProc, & +! dstLocalID) else dstProc = 0 @@ -482,8 +493,6 @@ function ice_HaloCreate(dist, nsBoundaryType, ewBoundaryType, & northBlock = ice_blocksGetNbrID(iblock, ice_blocksNorth, & ewBoundaryType, nsBoundaryType) - call ice_HaloMsgCreate(halo, dist, iblock, northBlock, 'north') - !*** set tripole flag and add two copies for inserting !*** and extracting info from the tripole buffer @@ -493,6 +502,7 @@ function ice_HaloCreate(dist, nsBoundaryType, ewBoundaryType, & call ice_HaloMsgCreate(halo, dist, -iblock, iblock, 'north') else tripoleBlock = .false. + call ice_HaloMsgCreate(halo, dist, iblock, northBlock, 'north') endif !*** find south neighbor block @@ -513,9 +523,10 @@ function ice_HaloCreate(dist, nsBoundaryType, ewBoundaryType, & !*** the east block to make sure enough information is !*** available for tripole manipulations - if (tripoleBlock) then - call ice_HaloMsgCreate(halo, dist, iblock, -eastBlock, 'north') - endif +! tcx,tcraig, 3/2023, this is not needed +! if (tripoleBlock) then +! call ice_HaloMsgCreate(halo, dist, iblock, -eastBlock, 'north') +! endif !*** find west neighbor block @@ -528,9 +539,10 @@ function ice_HaloCreate(dist, nsBoundaryType, ewBoundaryType, & !*** the west block to make sure enough information is !*** available for tripole manipulations - if (tripoleBlock) then - call ice_HaloMsgCreate(halo, dist, iblock, -westBlock, 'north') - endif +! tcx,tcraig, 3/2023, this is not needed +! if (tripoleBlock) then +! call ice_HaloMsgCreate(halo, dist, iblock, -westBlock, 'north') +! endif !*** find northeast neighbor block @@ -699,6 +711,23 @@ subroutine ice_HaloUpdate2DR8(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate2DR8)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -936,7 +965,7 @@ subroutine ice_HaloUpdate2DR8(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -1011,6 +1040,23 @@ subroutine ice_HaloUpdate2DR4(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate2DR4)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -1228,7 +1274,7 @@ subroutine ice_HaloUpdate2DR4(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -1303,6 +1349,23 @@ subroutine ice_HaloUpdate2DI4(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate2DI4)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -1520,7 +1583,7 @@ subroutine ice_HaloUpdate2DI4(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -1583,6 +1646,23 @@ subroutine ice_HaloUpdate2DL1(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate2DL1)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! copy logical into integer array and call haloupdate on integer array @@ -1662,6 +1742,23 @@ subroutine ice_HaloUpdate3DR8(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate3DR8)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -1896,7 +1993,7 @@ subroutine ice_HaloUpdate3DR8(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -1980,6 +2077,23 @@ subroutine ice_HaloUpdate3DR4(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate3DR4)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -2214,7 +2328,7 @@ subroutine ice_HaloUpdate3DR4(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -2298,6 +2412,23 @@ subroutine ice_HaloUpdate3DI4(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate3DI4)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -2532,7 +2663,7 @@ subroutine ice_HaloUpdate3DI4(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -2616,6 +2747,23 @@ subroutine ice_HaloUpdate4DR8(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate4DR8)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -2865,7 +3013,7 @@ subroutine ice_HaloUpdate4DR8(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -2951,6 +3099,23 @@ subroutine ice_HaloUpdate4DR4(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate4DR4)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -3200,7 +3365,7 @@ subroutine ice_HaloUpdate4DR4(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -3286,6 +3451,23 @@ subroutine ice_HaloUpdate4DI4(array, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate4DI4)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -3535,7 +3717,7 @@ subroutine ice_HaloUpdate4DI4(array, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface on u-fold, and NE corner and Nface @@ -3567,6 +3749,7 @@ end subroutine ice_HaloUpdate4DI4 !*********************************************************************** ! This routine updates ghost cells for an input array using ! a second array as needed by the stress fields. +! This is just like 2DR8 except no averaging and only on tripole subroutine ice_HaloUpdate_stress(array1, array2, halo, & fieldLoc, fieldKind, & @@ -3610,6 +3793,23 @@ subroutine ice_HaloUpdate_stress(array1, array2, halo, & character(len=*), parameter :: subname = '(ice_HaloUpdate_stress)' +!----------------------------------------------------------------------- +! +! abort or return on unknown or noupdate field_loc or field_type +! +!----------------------------------------------------------------------- + + if (fieldLoc == field_loc_unknown .or. & + fieldKind == field_type_unknown) then + call abort_ice(subname//'ERROR: use of field_loc/type_unknown not allowed') + return + endif + + if (fieldLoc == field_loc_noupdate .or. & + fieldKind == field_type_noupdate) then + return + endif + !----------------------------------------------------------------------- ! ! initialize error code and fill value @@ -3697,30 +3897,61 @@ subroutine ice_HaloUpdate_stress(array1, array2, halo, & call abort_ice(subname//'ERROR: Unknown field kind') end select - select case (fieldLoc) - case (field_loc_center) ! cell center location + if (halo%tripoleTFlag) then - ioffset = 0 - joffset = 0 + select case (fieldLoc) + case (field_loc_center) ! cell center location - case (field_loc_NEcorner) ! cell corner location + ioffset = -1 + joffset = 0 - ioffset = 1 - joffset = 1 + case (field_loc_NEcorner) ! cell corner location - case (field_loc_Eface) + ioffset = 0 + joffset = 1 - ioffset = 1 - joffset = 0 + case (field_loc_Eface) ! cell center location - case (field_loc_Nface) + ioffset = 0 + joffset = 0 - ioffset = 0 - joffset = 1 + case (field_loc_Nface) ! cell corner (velocity) location - case default - call abort_ice(subname//'ERROR: Unknown field location') - end select + ioffset = -1 + joffset = 1 + + case default + call abort_ice(subname//'ERROR: Unknown field location') + end select + + else ! tripole u-fold + + select case (fieldLoc) + case (field_loc_center) ! cell center location + + ioffset = 0 + joffset = 0 + + case (field_loc_NEcorner) ! cell corner location + + ioffset = 1 + joffset = 1 + + case (field_loc_Eface) + + ioffset = 1 + joffset = 0 + + case (field_loc_Nface) + + ioffset = 0 + joffset = 1 + + case default + call abort_ice(subname//'ERROR: Unknown field location') + end select + + endif !*** copy out of global tripole buffer into local !*** ghost cells @@ -3743,14 +3974,15 @@ subroutine ice_HaloUpdate_stress(array1, array2, halo, & !*** correct for offsets iSrc = iSrc - ioffset jSrc = jSrc - joffset - if (iSrc == 0) iSrc = nxGlobal + if (iSrc < 1 ) iSrc = iSrc + nxGlobal + if (iSrc > nxGlobal) iSrc = iSrc - nxGlobal !*** for center and Eface, do not need to replace !*** top row of physical domain, so jSrc should be !*** out of range and skipped !*** otherwise do the copy - if (jSrc <= nghost+1) then + if (jSrc <= halo%tripoleRows .and. jSrc>0 .and. jDst>0) then array1(iDst,jDst,dstBlock) = isign*bufTripoleR8(iSrc,jSrc) endif @@ -4159,36 +4391,37 @@ subroutine ice_HaloMsgCreate(halo, dist, srcBlock, dstBlock, direction) halo%numLocalCopies = msgIndx - else - - !*** tripole grid - copy entire top halo+1 - !*** rows into global buffer at src location - - msgIndx = halo%numLocalCopies - - do j=1,nghost+1 - do i=1,ieSrc-ibSrc+1 - - msgIndx = msgIndx + 1 - - halo%srcLocalAddr(1,msgIndx) = ibSrc + i - 1 - halo%srcLocalAddr(2,msgIndx) = jeSrc-1-nghost+j - halo%srcLocalAddr(3,msgIndx) = srcLocalID - - halo%dstLocalAddr(1,msgIndx) = iGlobal(ibSrc + i - 1) - halo%dstLocalAddr(2,msgIndx) = j - halo%dstLocalAddr(3,msgIndx) = -dstLocalID - - end do - end do - - halo%numLocalCopies = msgIndx +! tcx,tcraig, 3/2023, this is not needed +! else +! +! !*** tripole grid - copy entire top halo+1 +! !*** rows into global buffer at src location +! +! msgIndx = halo%numLocalCopies +! +! do j=1,nghost+1 +! do i=1,ieSrc-ibSrc+1 +! +! msgIndx = msgIndx + 1 +! +! halo%srcLocalAddr(1,msgIndx) = ibSrc + i - 1 +! halo%srcLocalAddr(2,msgIndx) = jeSrc-1-nghost+j +! halo%srcLocalAddr(3,msgIndx) = srcLocalID +! +! halo%dstLocalAddr(1,msgIndx) = iGlobal(ibSrc + i - 1) +! halo%dstLocalAddr(2,msgIndx) = j +! halo%dstLocalAddr(3,msgIndx) = -dstLocalID +! +! end do +! end do +! +! halo%numLocalCopies = msgIndx endif case ('northwest') - !*** normal northeast boundary - just copy NW corner + !*** normal northwest boundary - just copy NW corner !*** of physical domain into SE halo of NW nbr block if (dstBlock > 0) then @@ -4213,30 +4446,31 @@ subroutine ice_HaloMsgCreate(halo, dist, srcBlock, dstBlock, direction) halo%numLocalCopies = msgIndx - else - - !*** tripole grid - copy entire top halo+1 - !*** rows into global buffer at src location - - msgIndx = halo%numLocalCopies - - do j=1,nghost+1 - do i=1,ieSrc-ibSrc+1 - - msgIndx = msgIndx + 1 - - halo%srcLocalAddr(1,msgIndx) = ibSrc + i - 1 - halo%srcLocalAddr(2,msgIndx) = jeSrc-1-nghost+j - halo%srcLocalAddr(3,msgIndx) = srcLocalID - - halo%dstLocalAddr(1,msgIndx) = iGlobal(ibSrc + i - 1) - halo%dstLocalAddr(2,msgIndx) = j - halo%dstLocalAddr(3,msgIndx) = -dstLocalID - - end do - end do - - halo%numLocalCopies = msgIndx +! tcx,tcraig, 3/2023, this is not needed +! else +! +! !*** tripole grid - copy entire top halo+1 +! !*** rows into global buffer at src location +! +! msgIndx = halo%numLocalCopies +! +! do j=1,nghost+1 +! do i=1,ieSrc-ibSrc+1 +! +! msgIndx = msgIndx + 1 +! +! halo%srcLocalAddr(1,msgIndx) = ibSrc + i - 1 +! halo%srcLocalAddr(2,msgIndx) = jeSrc-1-nghost+j +! halo%srcLocalAddr(3,msgIndx) = srcLocalID +! +! halo%dstLocalAddr(1,msgIndx) = iGlobal(ibSrc + i - 1) +! halo%dstLocalAddr(2,msgIndx) = j +! halo%dstLocalAddr(3,msgIndx) = -dstLocalID +! +! end do +! end do +! +! halo%numLocalCopies = msgIndx endif @@ -4444,7 +4678,7 @@ subroutine ice_HaloMsgCreate(halo, dist, srcBlock, dstBlock, direction) case ('northwest') - !*** normal northeast boundary - just copy NW corner + !*** normal northwest boundary - just copy NW corner !*** of physical domain into SE halo of NW nbr block if (dstBlock > 0) then diff --git a/cicecore/cicedyn/infrastructure/ice_grid.F90 b/cicecore/cicedyn/infrastructure/ice_grid.F90 index 0d56d3400..770ee9ed9 100644 --- a/cicecore/cicedyn/infrastructure/ice_grid.F90 +++ b/cicecore/cicedyn/infrastructure/ice_grid.F90 @@ -360,6 +360,11 @@ subroutine init_grid1 file=__FILE__, line=__LINE__) endif + if (grid_type == 'tripole' .and. (mod(nx_global,2)/=0)) then + call abort_ice(subname//'ERROR: grid_type tripole requires even nx_global number', & + file=__FILE__, line=__LINE__) + endif + if (trim(grid_type) == 'displaced_pole' .or. & trim(grid_type) == 'tripole' .or. & trim(grid_type) == 'regional' ) then diff --git a/cicecore/cicedyn/infrastructure/io/io_binary/ice_history_write.F90 b/cicecore/cicedyn/infrastructure/io/io_binary/ice_history_write.F90 index 2a3f042c3..9df51635d 100644 --- a/cicecore/cicedyn/infrastructure/io/io_binary/ice_history_write.F90 +++ b/cicecore/cicedyn/infrastructure/io/io_binary/ice_history_write.F90 @@ -158,6 +158,7 @@ subroutine ice_write_hist(ns) trim(avail_hist_fields(n)%vcomment) if (histfreq(ns) == '1' .or. .not. hist_avg & + .or. write_ic & .or. n==n_divu(ns) .or. n==n_shear(ns) & ! snapshots .or. n==n_sig1(ns) .or. n==n_sig2(ns) & .or. n==n_sigP(ns) .or. n==n_trsig(ns) & @@ -186,7 +187,7 @@ subroutine ice_write_hist(ns) write (nu_hdr, 994) nrec,trim(avail_hist_fields(n)%vname), & trim(avail_hist_fields(n)%vdesc),trim(avail_hist_fields(n)%vunit),nn - if (histfreq(ns) == '1' .or. .not. hist_avg) then + if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then write (nu_hdr, 996) nrec,trim(avail_hist_fields(n)%vname), & 'time_rep','instantaneous' else @@ -210,7 +211,7 @@ subroutine ice_write_hist(ns) write (nu_hdr, 993) nrec,trim(avail_hist_fields(n)%vname), & trim(avail_hist_fields(n)%vdesc),trim(avail_hist_fields(n)%vunit),k - if (histfreq(ns) == '1' .or. .not. hist_avg) then + if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then write (nu_hdr, 996) nrec,trim(avail_hist_fields(n)%vname), & 'time_rep','instantaneous' else @@ -234,7 +235,7 @@ subroutine ice_write_hist(ns) write (nu_hdr, 993) nrec,trim(avail_hist_fields(n)%vname), & trim(avail_hist_fields(n)%vdesc),trim(avail_hist_fields(n)%vunit),nn,k - if (histfreq(ns) == '1' .or. .not. hist_avg) then + if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then write (nu_hdr, 996) nrec,trim(avail_hist_fields(n)%vname), & 'time_rep','instantaneous' else @@ -258,7 +259,7 @@ subroutine ice_write_hist(ns) write (nu_hdr, 993) nrec,trim(avail_hist_fields(n)%vname), & trim(avail_hist_fields(n)%vdesc),trim(avail_hist_fields(n)%vunit),nn,k - if (histfreq(ns) == '1' .or. .not. hist_avg) then + if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then write (nu_hdr, 996) nrec,trim(avail_hist_fields(n)%vname), & 'time_rep','instantaneous' else @@ -282,7 +283,7 @@ subroutine ice_write_hist(ns) write (nu_hdr, 993) nrec,trim(avail_hist_fields(n)%vname), & trim(avail_hist_fields(n)%vdesc),trim(avail_hist_fields(n)%vunit),nn,k - if (histfreq(ns) == '1' .or. .not. hist_avg) then + if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then write (nu_hdr, 996) nrec,trim(avail_hist_fields(n)%vname), & 'time_rep','instantaneous' else @@ -307,7 +308,7 @@ subroutine ice_write_hist(ns) write (nu_hdr, 993) nrec,trim(avail_hist_fields(n)%vname), & trim(avail_hist_fields(n)%vdesc),trim(avail_hist_fields(n)%vunit),nn,k - if (histfreq(ns) == '1' .or. .not. hist_avg) then + if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then write (nu_hdr, 996) nrec,trim(avail_hist_fields(n)%vname), & 'time_rep','instantaneous' else @@ -333,7 +334,7 @@ subroutine ice_write_hist(ns) write (nu_hdr, 993) nrec,trim(avail_hist_fields(n)%vname), & trim(avail_hist_fields(n)%vdesc),trim(avail_hist_fields(n)%vunit),nn,k - if (histfreq(ns) == '1' .or. .not. hist_avg) then + if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then write (nu_hdr, 996) nrec,trim(avail_hist_fields(n)%vname), & 'time_rep','instantaneous' else @@ -359,7 +360,7 @@ subroutine ice_write_hist(ns) write (nu_hdr, 993) nrec,trim(avail_hist_fields(n)%vname), & trim(avail_hist_fields(n)%vdesc),trim(avail_hist_fields(n)%vunit),nn,k - if (histfreq(ns) == '1' .or. .not. hist_avg) then + if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then write (nu_hdr, 996) nrec,trim(avail_hist_fields(n)%vname), & 'time_rep','instantaneous' else diff --git a/cicecore/cicedyn/infrastructure/io/io_netcdf/ice_history_write.F90 b/cicecore/cicedyn/infrastructure/io/io_netcdf/ice_history_write.F90 index d85ec5e3c..10d750300 100644 --- a/cicecore/cicedyn/infrastructure/io/io_netcdf/ice_history_write.F90 +++ b/cicecore/cicedyn/infrastructure/io/io_netcdf/ice_history_write.F90 @@ -159,7 +159,7 @@ subroutine ice_write_hist (ns) ! define dimensions !----------------------------------------------------------------- - if (hist_avg) then + if (hist_avg .and. .not. write_ic) then status = nf90_def_dim(ncid,'d2',2,boundid) if (status /= nf90_noerr) call abort_ice(subname// & 'ERROR: defining dim d2') @@ -241,7 +241,7 @@ subroutine ice_write_hist (ns) call abort_ice(subname//'ERROR: invalid calendar settings') endif - if (hist_avg) then + if (hist_avg .and. .not. write_ic) then status = nf90_put_att(ncid,varid,'bounds','time_bounds') if (status /= nf90_noerr) call abort_ice(subname// & 'ERROR: time bounds') @@ -251,7 +251,7 @@ subroutine ice_write_hist (ns) ! Define attributes for time bounds if hist_avg is true !----------------------------------------------------------------- - if (hist_avg) then + if (hist_avg .and. .not. write_ic) then dimid(1) = boundid dimid(2) = timid status = nf90_def_var(ncid,'time_bounds',lprecision,dimid(1:2),varid) @@ -745,7 +745,7 @@ subroutine ice_write_hist (ns) ! write time_bounds info !----------------------------------------------------------------- - if (hist_avg) then + if (hist_avg .and. .not. write_ic) then status = nf90_inq_varid(ncid,'time_bounds',varid) if (status /= nf90_noerr) call abort_ice(subname// & 'ERROR: getting time_bounds id') @@ -1236,7 +1236,7 @@ end subroutine ice_write_hist subroutine ice_write_hist_attrs(ncid, varid, hfield, ns) use ice_kinds_mod - use ice_calendar, only: histfreq, histfreq_n + use ice_calendar, only: histfreq, histfreq_n, write_ic use ice_history_shared, only: ice_hist_field, history_precision, & hist_avg #ifdef USE_NETCDF @@ -1279,7 +1279,7 @@ subroutine ice_write_hist_attrs(ncid, varid, hfield, ns) call ice_write_hist_fill(ncid,varid,hfield%vname,history_precision) ! Add cell_methods attribute to variables if averaged - if (hist_avg) then + if (hist_avg .and. .not. write_ic) then if (TRIM(hfield%vname(1:4))/='sig1' & .and.TRIM(hfield%vname(1:4))/='sig2' & .and.TRIM(hfield%vname(1:9))/='sistreave' & @@ -1293,6 +1293,7 @@ subroutine ice_write_hist_attrs(ncid, varid, hfield, ns) if ((histfreq(ns) == '1' .and. histfreq_n(ns) == 1) & .or..not. hist_avg & + .or. write_ic & .or.TRIM(hfield%vname(1:4))=='divu' & .or.TRIM(hfield%vname(1:5))=='shear' & .or.TRIM(hfield%vname(1:4))=='sig1' & diff --git a/cicecore/cicedyn/infrastructure/io/io_pio2/ice_history_write.F90 b/cicecore/cicedyn/infrastructure/io/io_pio2/ice_history_write.F90 index a697a98d5..25f9850ce 100644 --- a/cicecore/cicedyn/infrastructure/io/io_pio2/ice_history_write.F90 +++ b/cicecore/cicedyn/infrastructure/io/io_pio2/ice_history_write.F90 @@ -195,7 +195,7 @@ subroutine ice_write_hist (ns) ! define dimensions !----------------------------------------------------------------- - if (hist_avg) then + if (hist_avg .and. .not. write_ic) then status = pio_def_dim(File,'d2',2,boundid) endif @@ -233,12 +233,12 @@ subroutine ice_write_hist (ns) call abort_ice(subname//'ERROR: invalid calendar settings') endif - if (hist_avg) then + if (hist_avg .and. .not. write_ic) then status = pio_put_att(File,varid,'bounds','time_bounds') endif ! Define attributes for time_bounds if hist_avg is true - if (hist_avg) then + if (hist_avg .and. .not. write_ic) then dimid2(1) = boundid dimid2(2) = timid status = pio_def_var(File,'time_bounds',pio_double,dimid2,varid) @@ -702,7 +702,7 @@ subroutine ice_write_hist (ns) ! write time_bounds info !----------------------------------------------------------------- - if (hist_avg) then + if (hist_avg .and. .not. write_ic) then status = pio_inq_varid(File,'time_bounds',varid) time_bounds=(/time_beg(ns),time_end(ns)/) bnd_start = (/1,1/) @@ -1219,7 +1219,7 @@ end subroutine ice_write_hist subroutine ice_write_hist_attrs(File, varid, hfield, ns) use ice_kinds_mod - use ice_calendar, only: histfreq, histfreq_n + use ice_calendar, only: histfreq, histfreq_n, write_ic use ice_history_shared, only: ice_hist_field, history_precision, & hist_avg use ice_pio @@ -1250,7 +1250,7 @@ subroutine ice_write_hist_attrs(File, varid, hfield, ns) call ice_write_hist_fill(File,varid,hfield%vname,history_precision) ! Add cell_methods attribute to variables if averaged - if (hist_avg) then + if (hist_avg .and. .not. write_ic) then if (TRIM(hfield%vname(1:4))/='sig1' & .and.TRIM(hfield%vname(1:4))/='sig2' & .and.TRIM(hfield%vname(1:9))/='sistreave' & @@ -1262,6 +1262,7 @@ subroutine ice_write_hist_attrs(File, varid, hfield, ns) if ((histfreq(ns) == '1' .and. histfreq_n(ns) == 1) & .or..not. hist_avg & + .or. write_ic & .or.TRIM(hfield%vname(1:4))=='divu' & .or.TRIM(hfield%vname(1:5))=='shear' & .or.TRIM(hfield%vname(1:4))=='sig1' & diff --git a/cicecore/drivers/direct/nemo_concepts/CICE_InitMod.F90 b/cicecore/drivers/direct/nemo_concepts/CICE_InitMod.F90 index 85050d8c9..3a8f5e33d 100644 --- a/cicecore/drivers/direct/nemo_concepts/CICE_InitMod.F90 +++ b/cicecore/drivers/direct/nemo_concepts/CICE_InitMod.F90 @@ -185,6 +185,8 @@ subroutine cice_init if (trim(runtype) == 'continue' .or. restart) & call init_shortwave ! initialize radiative transfer + if (write_ic) call accum_hist(dt) ! write initial conditions + ! determine the time and date at the end of the first timestep call advance_timestep() @@ -215,8 +217,6 @@ subroutine cice_init call init_flux_atm ! initialize atmosphere fluxes sent to coupler call init_flux_ocn ! initialize ocean fluxes sent to coupler - if (write_ic) call accum_hist(dt) ! write initial conditions - end subroutine cice_init !======================================================================= diff --git a/cicecore/drivers/standalone/cice/CICE_InitMod.F90 b/cicecore/drivers/standalone/cice/CICE_InitMod.F90 index 8de05a121..0371c7f38 100644 --- a/cicecore/drivers/standalone/cice/CICE_InitMod.F90 +++ b/cicecore/drivers/standalone/cice/CICE_InitMod.F90 @@ -199,6 +199,8 @@ subroutine cice_init if (trim(runtype) == 'continue' .or. restart) & call init_shortwave ! initialize radiative transfer + if (write_ic) call accum_hist(dt) ! write initial conditions + ! tcraig, use advance_timestep here ! istep = istep + 1 ! update time step counters ! istep1 = istep1 + 1 @@ -243,8 +245,6 @@ subroutine cice_init call init_flux_atm ! initialize atmosphere fluxes sent to coupler call init_flux_ocn ! initialize ocean fluxes sent to coupler - if (write_ic) call accum_hist(dt) ! write initial conditions - if (my_task == master_task) then call ice_memusage_print(nu_diag,subname//':end') endif diff --git a/cicecore/drivers/unittest/halochk/CICE_InitMod.F90 b/cicecore/drivers/unittest/halochk/CICE_InitMod.F90 new file mode 100644 index 000000000..9ed1c5cbc --- /dev/null +++ b/cicecore/drivers/unittest/halochk/CICE_InitMod.F90 @@ -0,0 +1,472 @@ +!======================================================================= +! +! This module contains the CICE initialization routine that sets model +! parameters and initializes the grid and CICE state variables. +! +! authors Elizabeth C. Hunke, LANL +! William H. Lipscomb, LANL +! Philip W. Jones, LANL +! +! 2006: Converted to free form source (F90) by Elizabeth Hunke +! 2008: E. Hunke moved ESMF code to its own driver + + module CICE_InitMod + + use ice_kinds_mod + use ice_exit, only: abort_ice + use ice_fileunits, only: init_fileunits, nu_diag + use icepack_intfc, only: icepack_aggregate + use icepack_intfc, only: icepack_init_itd, icepack_init_itd_hist + use icepack_intfc, only: icepack_init_fsd_bounds, icepack_init_wave + use icepack_intfc, only: icepack_configure + use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted + use icepack_intfc, only: icepack_query_parameters, icepack_query_tracer_flags, & + icepack_query_tracer_indices, icepack_query_tracer_sizes + + implicit none + private + public :: CICE_Initialize, cice_init + +!======================================================================= + + contains + +!======================================================================= + +! Initialize the basic state, grid and all necessary parameters for +! running the CICE model. Return the initial state in routine +! export state. +! Note: This initialization driver is designed for standalone and +! CESM-coupled applications. For other +! applications (e.g., standalone CAM), this driver would be +! replaced by a different driver that calls subroutine cice_init, +! where most of the work is done. + + subroutine CICE_Initialize + + character(len=*), parameter :: subname='(CICE_Initialize)' + !-------------------------------------------------------------------- + ! model initialization + !-------------------------------------------------------------------- + + call cice_init + + end subroutine CICE_Initialize + +!======================================================================= +! +! Initialize CICE model. + + subroutine cice_init + + use ice_arrays_column, only: hin_max, c_hi_range, alloc_arrays_column + use ice_arrays_column, only: floe_rad_l, floe_rad_c, & + floe_binwidth, c_fsd_range + use ice_state, only: alloc_state + use ice_flux_bgc, only: alloc_flux_bgc + use ice_calendar, only: dt, dt_dyn, istep, istep1, write_ic, & + init_calendar, advance_timestep, calc_timesteps + use ice_communicate, only: init_communicate, my_task, master_task + use ice_diagnostics, only: init_diags + use ice_domain, only: init_domain_blocks + use ice_domain_size, only: ncat, nfsd + use ice_dyn_eap, only: init_eap + use ice_dyn_evp, only: init_evp + use ice_dyn_vp, only: init_vp + use ice_dyn_shared, only: kdyn + use ice_flux, only: init_coupler_flux, init_history_therm, & + init_history_dyn, init_flux_atm, init_flux_ocn, alloc_flux + use ice_forcing, only: init_forcing_ocn, init_forcing_atmo, & + get_forcing_atmo, get_forcing_ocn, get_wave_spec + use ice_forcing_bgc, only: get_forcing_bgc, get_atm_bgc, & + faero_default, faero_optics, alloc_forcing_bgc, fiso_default + use ice_grid, only: init_grid1, init_grid2, alloc_grid + use ice_history, only: init_hist, accum_hist + use ice_restart_shared, only: restart, runtype + use ice_init, only: input_data, init_state + use ice_init_column, only: init_thermo_vertical, init_shortwave, init_zbgc, input_zbgc, count_tracers + use ice_kinds_mod + use ice_restoring, only: ice_HaloRestore_init + use ice_timers, only: timer_total, init_ice_timers, ice_timer_start + use ice_transport_driver, only: init_transport + + logical(kind=log_kind) :: tr_aero, tr_zaero, skl_bgc, z_tracers, & + tr_iso, tr_fsd, wave_spec + character(len=*), parameter :: subname = '(cice_init)' + + call init_communicate ! initial setup for message passing + call init_fileunits ! unit numbers + + ! tcx debug, this will create a different logfile for each pe + ! if (my_task /= master_task) nu_diag = 100+my_task + + call icepack_configure() ! initialize icepack + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(trim(subname), & + file=__FILE__,line= __LINE__) + + call input_data ! namelist variables + call input_zbgc ! vertical biogeochemistry namelist + call count_tracers ! count tracers + + call init_domain_blocks ! set up block decomposition + call init_grid1 ! domain distribution + call alloc_grid ! allocate grid arrays + call alloc_arrays_column ! allocate column arrays + call alloc_state ! allocate state arrays + call alloc_flux_bgc ! allocate flux_bgc arrays + call alloc_flux ! allocate flux arrays + call init_ice_timers ! initialize all timers + call ice_timer_start(timer_total) ! start timing entire run + call init_grid2 ! grid variables + call init_zbgc ! vertical biogeochemistry initialization + call init_calendar ! initialize some calendar stuff + call init_hist (dt) ! initialize output history file + + if (kdyn == 1) then + call init_evp + else if (kdyn == 2) then + call init_eap ! define eap dynamics parameters, variables + else if (kdyn == 3) then + call init_vp ! define vp dynamics parameters, variables + endif + + call init_coupler_flux ! initialize fluxes exchanged with coupler + + call init_thermo_vertical ! initialize vertical thermodynamics + + call icepack_init_itd(ncat=ncat, hin_max=hin_max) ! ice thickness distribution + if (my_task == master_task) then + call icepack_init_itd_hist(ncat=ncat, hin_max=hin_max, c_hi_range=c_hi_range) ! output + endif + + call icepack_query_tracer_flags(tr_fsd_out=tr_fsd) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(trim(subname), & + file=__FILE__,line= __LINE__) + + if (tr_fsd) call icepack_init_fsd_bounds (nfsd, & ! floe size distribution + floe_rad_l, & ! fsd size lower bound in m (radius) + floe_rad_c, & ! fsd size bin centre in m (radius) + floe_binwidth, & ! fsd size bin width in m (radius) + c_fsd_range, & ! string for history output + write_diags=(my_task == master_task)) ! write diag on master only + + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) + + call init_forcing_ocn(dt) ! initialize sss and sst from data + call init_state ! initialize the ice state + call init_transport ! initialize horizontal transport + call ice_HaloRestore_init ! restored boundary conditions + + call icepack_query_parameters(skl_bgc_out=skl_bgc, z_tracers_out=z_tracers, & + wave_spec_out=wave_spec) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(trim(subname), & + file=__FILE__,line= __LINE__) + + if (skl_bgc .or. z_tracers) call alloc_forcing_bgc ! allocate biogeochemistry arrays + + call init_restart ! initialize restart variables + call init_diags ! initialize diagnostic output points + call init_history_therm ! initialize thermo history variables + call init_history_dyn ! initialize dynamic history variables + call calc_timesteps ! update timestep counter if not using npt_unit="1" + + call icepack_query_tracer_flags(tr_aero_out=tr_aero, tr_zaero_out=tr_zaero) + call icepack_query_tracer_flags(tr_iso_out=tr_iso) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(trim(subname), & + file=__FILE__,line= __LINE__) + + if (tr_aero .or. tr_zaero) call faero_optics !initialize aerosol optical + !property tables + + ! Initialize shortwave components using swdn from previous timestep + ! if restarting. These components will be scaled to current forcing + ! in prep_radiation. + if (trim(runtype) == 'continue' .or. restart) & + call init_shortwave ! initialize radiative transfer + +! tcraig, use advance_timestep here +! istep = istep + 1 ! update time step counters +! istep1 = istep1 + 1 +! time = time + dt ! determine the time and date +! call calendar(time) ! at the end of the first timestep + call advance_timestep() + + !-------------------------------------------------------------------- + ! coupler communication or forcing data initialization + !-------------------------------------------------------------------- + + call init_forcing_atmo ! initialize atmospheric forcing (standalone) + + if (tr_fsd .and. wave_spec) call get_wave_spec ! wave spectrum in ice + call get_forcing_atmo ! atmospheric forcing from data + call get_forcing_ocn(dt) ! ocean forcing from data + + ! isotopes + if (tr_iso) call fiso_default ! default values + ! aerosols + ! if (tr_aero) call faero_data ! data file + ! if (tr_zaero) call fzaero_data ! data file (gx1) + if (tr_aero .or. tr_zaero) call faero_default ! default values + if (skl_bgc .or. z_tracers) call get_forcing_bgc ! biogeochemistry + if (z_tracers) call get_atm_bgc ! biogeochemistry + + if (runtype == 'initial' .and. .not. restart) & + call init_shortwave ! initialize radiative transfer using current swdn + + call init_flux_atm ! initialize atmosphere fluxes sent to coupler + call init_flux_ocn ! initialize ocean fluxes sent to coupler + + if (write_ic) call accum_hist(dt) ! write initial conditions + + end subroutine cice_init + +!======================================================================= + + subroutine init_restart + + use ice_arrays_column, only: dhsn + use ice_blocks, only: nx_block, ny_block + use ice_calendar, only: calendar + use ice_constants, only: c0 + use ice_domain, only: nblocks + use ice_domain_size, only: ncat, n_iso, n_aero, nfsd + use ice_dyn_eap, only: read_restart_eap + use ice_dyn_shared, only: kdyn + use ice_grid, only: tmask + use ice_init, only: ice_ic + use ice_init_column, only: init_age, init_FY, init_lvl, & + init_meltponds_lvl, init_meltponds_topo, & + init_isotope, init_aerosol, init_hbrine, init_bgc, init_fsd + use ice_restart_column, only: restart_age, read_restart_age, & + restart_FY, read_restart_FY, restart_lvl, read_restart_lvl, & + restart_pond_lvl, read_restart_pond_lvl, & + restart_pond_topo, read_restart_pond_topo, & + restart_fsd, read_restart_fsd, & + restart_iso, read_restart_iso, & + restart_aero, read_restart_aero, & + restart_hbrine, read_restart_hbrine, & + restart_zsal, restart_bgc + use ice_restart_driver, only: restartfile + use ice_restart_shared, only: runtype, restart + use ice_state ! almost everything + + integer(kind=int_kind) :: & + i, j , & ! horizontal indices + iblk ! block index + logical(kind=log_kind) :: & + tr_iage, tr_FY, tr_lvl, tr_pond_lvl, & + tr_pond_topo, tr_fsd, tr_iso, tr_aero, tr_brine, & + skl_bgc, z_tracers, solve_zsal + integer(kind=int_kind) :: & + ntrcr + integer(kind=int_kind) :: & + nt_alvl, nt_vlvl, nt_apnd, nt_hpnd, nt_ipnd, & + nt_iage, nt_FY, nt_aero, nt_fsd, nt_isosno, nt_isoice + + character(len=*), parameter :: subname = '(init_restart)' + + call icepack_query_tracer_sizes(ntrcr_out=ntrcr) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) + + call icepack_query_parameters(skl_bgc_out=skl_bgc, & + z_tracers_out=z_tracers, solve_zsal_out=solve_zsal) + call icepack_query_tracer_flags(tr_iage_out=tr_iage, tr_FY_out=tr_FY, & + tr_lvl_out=tr_lvl, tr_pond_lvl_out=tr_pond_lvl, & + tr_pond_topo_out=tr_pond_topo, tr_aero_out=tr_aero, tr_brine_out=tr_brine, & + tr_fsd_out=tr_fsd, tr_iso_out=tr_iso) + call icepack_query_tracer_indices(nt_alvl_out=nt_alvl, nt_vlvl_out=nt_vlvl, & + nt_apnd_out=nt_apnd, nt_hpnd_out=nt_hpnd, nt_ipnd_out=nt_ipnd, & + nt_iage_out=nt_iage, nt_FY_out=nt_FY, nt_aero_out=nt_aero, nt_fsd_out=nt_fsd, & + nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) + + if (trim(runtype) == 'continue') then + ! start from core restart file + call restartfile() ! given by pointer in ice_in + call calendar() ! update time parameters + if (kdyn == 2) call read_restart_eap ! EAP + else if (restart) then ! ice_ic = core restart file + call restartfile (ice_ic) ! or 'default' or 'none' + !!! uncomment to create netcdf + ! call restartfile_v4 (ice_ic) ! CICE v4.1 binary restart file + !!! uncomment if EAP restart data exists + ! if (kdyn == 2) call read_restart_eap + endif + + ! tracers + ! ice age tracer + if (tr_iage) then + if (trim(runtype) == 'continue') & + restart_age = .true. + if (restart_age) then + call read_restart_age + else + do iblk = 1, nblocks + call init_age(trcrn(:,:,nt_iage,:,iblk)) + enddo ! iblk + endif + endif + ! first-year area tracer + if (tr_FY) then + if (trim(runtype) == 'continue') restart_FY = .true. + if (restart_FY) then + call read_restart_FY + else + do iblk = 1, nblocks + call init_FY(trcrn(:,:,nt_FY,:,iblk)) + enddo ! iblk + endif + endif + ! level ice tracer + if (tr_lvl) then + if (trim(runtype) == 'continue') restart_lvl = .true. + if (restart_lvl) then + call read_restart_lvl + else + do iblk = 1, nblocks + call init_lvl(iblk,trcrn(:,:,nt_alvl,:,iblk), & + trcrn(:,:,nt_vlvl,:,iblk)) + enddo ! iblk + endif + endif + ! level-ice melt ponds + if (tr_pond_lvl) then + if (trim(runtype) == 'continue') & + restart_pond_lvl = .true. + if (restart_pond_lvl) then + call read_restart_pond_lvl + else + do iblk = 1, nblocks + call init_meltponds_lvl(trcrn(:,:,nt_apnd,:,iblk), & + trcrn(:,:,nt_hpnd,:,iblk), & + trcrn(:,:,nt_ipnd,:,iblk), & + dhsn(:,:,:,iblk)) + enddo ! iblk + endif + endif + ! topographic melt ponds + if (tr_pond_topo) then + if (trim(runtype) == 'continue') & + restart_pond_topo = .true. + if (restart_pond_topo) then + call read_restart_pond_topo + else + do iblk = 1, nblocks + call init_meltponds_topo(trcrn(:,:,nt_apnd,:,iblk), & + trcrn(:,:,nt_hpnd,:,iblk), & + trcrn(:,:,nt_ipnd,:,iblk)) + enddo ! iblk + endif ! .not. restart_pond + endif + ! floe size distribution + if (tr_fsd) then + if (trim(runtype) == 'continue') restart_fsd = .true. + if (restart_fsd) then + call read_restart_fsd + else + call init_fsd(trcrn(:,:,nt_fsd:nt_fsd+nfsd-1,:,:)) + endif + endif + + ! isotopes + if (tr_iso) then + if (trim(runtype) == 'continue') restart_iso = .true. + if (restart_iso) then + call read_restart_iso + else + do iblk = 1, nblocks + call init_isotope(trcrn(:,:,nt_isosno:nt_isosno+n_iso-1,:,iblk), & + trcrn(:,:,nt_isoice:nt_isoice+n_iso-1,:,iblk)) + enddo ! iblk + endif + endif + + if (tr_aero) then ! ice aerosol + if (trim(runtype) == 'continue') restart_aero = .true. + if (restart_aero) then + call read_restart_aero + else + do iblk = 1, nblocks + call init_aerosol(trcrn(:,:,nt_aero:nt_aero+4*n_aero-1,:,iblk)) + enddo ! iblk + endif ! .not. restart_aero + endif + + if (trim(runtype) == 'continue') then + if (tr_brine) & + restart_hbrine = .true. + if (solve_zsal) & + restart_zsal = .true. + if (skl_bgc .or. z_tracers) & + restart_bgc = .true. + endif + + if (tr_brine .or. skl_bgc) then ! brine height tracer + call init_hbrine + if (tr_brine .and. restart_hbrine) call read_restart_hbrine + endif + + if (solve_zsal .or. skl_bgc .or. z_tracers) then ! biogeochemistry + if (tr_fsd) then + write (nu_diag,*) 'FSD implementation incomplete for use with BGC' + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) + endif + call init_bgc + endif + + !----------------------------------------------------------------- + ! aggregate tracers + !----------------------------------------------------------------- + + !$OMP PARALLEL DO PRIVATE(iblk) + do iblk = 1, nblocks + do j = 1, ny_block + do i = 1, nx_block + if (tmask(i,j,iblk)) then + call icepack_aggregate(ncat = ncat, & + aicen = aicen(i,j,:,iblk), & + trcrn = trcrn(i,j,:,:,iblk), & + vicen = vicen(i,j,:,iblk), & + vsnon = vsnon(i,j,:,iblk), & + aice = aice (i,j, iblk), & + trcr = trcr (i,j,:,iblk), & + vice = vice (i,j, iblk), & + vsno = vsno (i,j, iblk), & + aice0 = aice0(i,j, iblk), & + ntrcr = ntrcr, & + trcr_depend = trcr_depend, & + trcr_base = trcr_base, & + n_trcr_strata = n_trcr_strata, & + nt_strata = nt_strata) + else + ! tcraig, reset all tracer values on land to zero + trcrn(i,j,:,:,iblk) = c0 + endif + enddo + enddo + enddo + !$OMP END PARALLEL DO + + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call abort_ice(error_message=subname, & + file=__FILE__, line=__LINE__) + + end subroutine init_restart + +!======================================================================= + + end module CICE_InitMod + +!======================================================================= diff --git a/cicecore/drivers/unittest/halochk/halochk.F90 b/cicecore/drivers/unittest/halochk/halochk.F90 new file mode 100644 index 000000000..29eaa8150 --- /dev/null +++ b/cicecore/drivers/unittest/halochk/halochk.F90 @@ -0,0 +1,811 @@ + + module halochk_data + + use CICE_InitMod + use ice_kinds_mod, only: int_kind, dbl_kind, real_kind, log_kind + use ice_blocks, only: block, get_block, nx_block, ny_block, nblocks_tot, nghost + use ice_boundary, only: ice_HaloUpdate, ice_HaloUpdate_stress + use ice_constants, only: c0, c1, p5, & + field_loc_unknown, field_loc_noupdate, & + field_loc_center, field_loc_NEcorner, & + field_loc_Nface, field_loc_Eface, & + field_type_unknown, field_type_noupdate, & + field_type_scalar, field_type_vector, field_type_angle + use ice_communicate, only: my_task, master_task, get_num_procs, MPI_COMM_ICE + use ice_distribution, only: ice_distributionGetBlockID, ice_distributionGet + use ice_domain_size, only: nx_global, ny_global, & + block_size_x, block_size_y, max_blocks + use ice_domain, only: distrb_info, halo_info, & + ew_boundary_type, ns_boundary_type + use ice_exit, only: abort_ice, end_run + use ice_global_reductions, only: global_minval, global_maxval, global_sum + + implicit none + + integer(int_kind), parameter :: & + passflag = 0, & + failflag = 1 + + end module halochk_data + +!======================================================================= + + program halochk + + ! This tests the CICE halo update methods by + ! using CICE_InitMod (from the standalone model) to read/initialize + ! a CICE grid/configuration. + + use halochk_data + + implicit none + + integer(int_kind) :: nn, nl, nt, i, j, k1, k2, n, ib, ie, jb, je + integer(int_kind) :: iblock, itrip, ioffset, joffset + integer(int_kind) :: blockID, numBlocks, jtrip + type (block) :: this_block + + real(dbl_kind) , allocatable :: darrayi1(:,:,:) , darrayj1(:,:,:) + real(dbl_kind) , allocatable :: darrayi2(:,:,:,:) , darrayj2(:,:,:,:) + real(dbl_kind) , allocatable :: darrayi3(:,:,:,:,:), darrayj3(:,:,:,:,:) + real(real_kind) , allocatable :: rarrayi1(:,:,:) , rarrayj1(:,:,:) + real(real_kind) , allocatable :: rarrayi2(:,:,:,:) , rarrayj2(:,:,:,:) + real(real_kind) , allocatable :: rarrayi3(:,:,:,:,:), rarrayj3(:,:,:,:,:) + integer(int_kind), allocatable :: iarrayi1(:,:,:) , iarrayj1(:,:,:) + integer(int_kind), allocatable :: iarrayi2(:,:,:,:) , iarrayj2(:,:,:,:) + integer(int_kind), allocatable :: iarrayi3(:,:,:,:,:), iarrayj3(:,:,:,:,:) + logical(log_kind), allocatable :: larrayi1(:,:,:) , larrayj1(:,:,:) + real(dbl_kind) , allocatable :: darrayi1str(:,:,:) , darrayj1str(:,:,:) + real(dbl_kind) , allocatable :: darrayi10(:,:,:) , darrayj10(:,:,:) + + real(dbl_kind), allocatable :: cidata_bas(:,:,:,:,:),cjdata_bas(:,:,:,:,:) + real(dbl_kind), allocatable :: cidata_nup(:,:,:,:,:),cjdata_nup(:,:,:,:,:) + real(dbl_kind), allocatable :: cidata_std(:,:,:,:,:),cjdata_std(:,:,:,:,:) + + integer(int_kind), parameter :: maxtests = 11 + integer(int_kind), parameter :: maxtypes = 4 + integer(int_kind), parameter :: maxlocs = 5 + integer(int_kind), parameter :: nz1 = 3 + integer(int_kind), parameter :: nz2 = 4 + real(dbl_kind) :: aichk,ajchk,cichk,cjchk,rival,rjval,rsign + character(len=16) :: locs_name(maxlocs), types_name(maxtypes) + integer(int_kind) :: field_loc(maxlocs), field_type(maxtypes) + integer(int_kind) :: npes, ierr, ntask, testcnt, tottest, tpcnt, tfcnt + integer(int_kind) :: errorflag0, gflag, k1m, k2m, ptcntsum, failcntsum + integer(int_kind), allocatable :: errorflag(:) + integer(int_kind), allocatable :: ptcnt(:), failcnt(:) + character(len=128), allocatable :: teststring(:) + character(len=32) :: halofld + logical :: tripole_average, tripole_pole, spvalL1 + logical :: first_call = .true. + + real(dbl_kind) , parameter :: fillval = -88888.0_dbl_kind + real(dbl_kind) , parameter :: dhalofillval = -999.0_dbl_kind + real(real_kind) , parameter :: rhalofillval = -999.0_real_kind + integer(int_kind), parameter :: ihalofillval = -999 + character(len=*) , parameter :: subname='(halochk)' + + !----------------------------------------------------------------- + ! Initialize CICE + !----------------------------------------------------------------- + + call CICE_Initialize + npes = get_num_procs() + + locs_name (:) = 'unknown' + types_name(:) = 'unknown' + field_type(:) = field_type_unknown + field_loc (:) = field_loc_unknown + + types_name(1) = 'scalar' + field_type(1) = field_type_scalar + types_name(2) = 'vector' + field_type(2) = field_type_vector + types_name(3) = 'angle' + field_type(3) = field_type_angle + types_name(4) = 'none' + field_type(4) = field_type_noupdate +! types_name(5) = 'unknown' +! field_type(5) = field_type_unknown ! aborts in CICE, as expected + + locs_name (1) = 'center' + field_loc (1) = field_loc_center + locs_name (2) = 'NEcorner' + field_loc (2) = field_loc_NEcorner + locs_name (3) = 'Nface' + field_loc (3) = field_loc_Nface + locs_name (4) = 'Eface' + field_loc (4) = field_loc_Eface + locs_name (5) = 'none' + field_loc (5) = field_loc_noupdate +! locs_name (6) = 'unknown' +! field_loc (6) = field_loc_unknown ! aborts in CICE, as expected + + tottest = maxtests * maxlocs * maxtypes + allocate(errorflag(tottest)) + allocate(teststring(tottest)) + allocate(ptcnt(tottest)) + allocate(failcnt(tottest)) + ptcnt(:) = 0 + failcnt(:) = 0 + + !----------------------------------------------------------------- + ! Testing + !----------------------------------------------------------------- + + if (my_task == master_task) then + write(6,*) ' ' + write(6,*) '==========================================================' + write(6,*) ' ' + write(6,*) 'RunningUnitTest HALOCHK' + write(6,*) ' ' + write(6,*) ' npes = ',npes + write(6,*) ' my_task = ',my_task + write(6,*) ' nx_global = ',nx_global + write(6,*) ' ny_global = ',ny_global + write(6,*) ' block_size_x = ',block_size_x + write(6,*) ' block_size_y = ',block_size_y + write(6,*) ' nblocks_tot = ',nblocks_tot + write(6,*) ' tottest = ',tottest + write(6,*) ' ' + endif + + errorflag0 = passflag + errorflag(:) = passflag + teststring(:) = ' ' + + ! --------------------------- + ! TEST HALO UPDATE + ! --------------------------- + + if (my_task == master_task) write(6,*) ' ' + + allocate(darrayi1 (nx_block,ny_block,max_blocks)) + allocate(darrayj1 (nx_block,ny_block,max_blocks)) + allocate(darrayi2 (nx_block,ny_block,nz1,max_blocks)) + allocate(darrayj2 (nx_block,ny_block,nz1,max_blocks)) + allocate(darrayi3 (nx_block,ny_block,nz1,nz2,max_blocks)) + allocate(darrayj3 (nx_block,ny_block,nz1,nz2,max_blocks)) + allocate(rarrayi1 (nx_block,ny_block,max_blocks)) + allocate(rarrayj1 (nx_block,ny_block,max_blocks)) + allocate(rarrayi2 (nx_block,ny_block,nz1,max_blocks)) + allocate(rarrayj2 (nx_block,ny_block,nz1,max_blocks)) + allocate(rarrayi3 (nx_block,ny_block,nz1,nz2,max_blocks)) + allocate(rarrayj3 (nx_block,ny_block,nz1,nz2,max_blocks)) + allocate(iarrayi1 (nx_block,ny_block,max_blocks)) + allocate(iarrayj1 (nx_block,ny_block,max_blocks)) + allocate(iarrayi2 (nx_block,ny_block,nz1,max_blocks)) + allocate(iarrayj2 (nx_block,ny_block,nz1,max_blocks)) + allocate(iarrayi3 (nx_block,ny_block,nz1,nz2,max_blocks)) + allocate(iarrayj3 (nx_block,ny_block,nz1,nz2,max_blocks)) + allocate(larrayi1 (nx_block,ny_block,max_blocks)) + allocate(larrayj1 (nx_block,ny_block,max_blocks)) + allocate(darrayi1str(nx_block,ny_block,max_blocks)) + allocate(darrayj1str(nx_block,ny_block,max_blocks)) + allocate(darrayi10 (nx_block,ny_block,max_blocks)) + allocate(darrayj10 (nx_block,ny_block,max_blocks)) + + allocate(cidata_bas(nx_block,ny_block,nz1,nz2,max_blocks)) + allocate(cjdata_bas(nx_block,ny_block,nz1,nz2,max_blocks)) + allocate(cidata_std(nx_block,ny_block,nz1,nz2,max_blocks)) + allocate(cjdata_std(nx_block,ny_block,nz1,nz2,max_blocks)) + allocate(cidata_nup(nx_block,ny_block,nz1,nz2,max_blocks)) + allocate(cjdata_nup(nx_block,ny_block,nz1,nz2,max_blocks)) + + darrayi1 = fillval + darrayj1 = fillval + darrayi2 = fillval + darrayj2 = fillval + darrayi3 = fillval + darrayj3 = fillval + rarrayi1 = fillval + rarrayj1 = fillval + rarrayi2 = fillval + rarrayj2 = fillval + rarrayi3 = fillval + rarrayj3 = fillval + iarrayi1 = fillval + iarrayj1 = fillval + iarrayi2 = fillval + iarrayj2 = fillval + iarrayi3 = fillval + iarrayj3 = fillval + larrayi1 = .false. + larrayj1 = .true. + darrayi1str = fillval + darrayj1str = fillval + darrayi10 = fillval + darrayj10 = fillval + cidata_bas = fillval + cjdata_bas = fillval + cidata_std = fillval + cjdata_std = fillval + cidata_nup = fillval + cjdata_nup = fillval + + call ice_distributionGet(distrb_info, numLocalBlocks = numBlocks) + + !--- baseline data --- + + do iblock = 1,numBlocks + call ice_distributionGetBlockID(distrb_info, iblock, blockID) + this_block = get_block(blockID, blockID) + do k2 = 1,nz2 + do k1 = 1,nz1 + do j = 1,ny_block + do i = 1,nx_block + cidata_bas(i,j,k1,k2,iblock) = real(this_block%i_glob(i),kind=dbl_kind) + & + real(k1,kind=dbl_kind)*1000._dbl_kind + real(k2,kind=dbl_kind)*10000._dbl_kind + cjdata_bas(i,j,k1,k2,iblock) = real(this_block%j_glob(j),kind=dbl_kind) + & + real(k1,kind=dbl_kind)*1000._dbl_kind + real(k2,kind=dbl_kind)*10000._dbl_kind + enddo + enddo + enddo + enddo + enddo + + !--- setup nup (noupdate) solution, set halo/pad will fillval --- + + cidata_nup(:,:,:,:,:) = cidata_bas(:,:,:,:,:) + cjdata_nup(:,:,:,:,:) = cjdata_bas(:,:,:,:,:) + + do iblock = 1,numBlocks + call ice_distributionGetBlockID(distrb_info, iblock, blockID) + this_block = get_block(blockID, blockID) + ib = this_block%ilo + ie = this_block%ihi + jb = this_block%jlo + je = this_block%jhi + cidata_nup(1:ib-1 ,: ,:,:,iblock) = fillval + cjdata_nup(1:ib-1 ,: ,:,:,iblock) = fillval + cidata_nup(ie+1:nx_block,: ,:,:,iblock) = fillval + cjdata_nup(ie+1:nx_block,: ,:,:,iblock) = fillval + cidata_nup(: ,1:jb-1 ,:,:,iblock) = fillval + cjdata_nup(: ,1:jb-1 ,:,:,iblock) = fillval + cidata_nup(: ,je+1:ny_block,:,:,iblock) = fillval + cjdata_nup(: ,je+1:ny_block,:,:,iblock) = fillval + enddo + + !--- setup std solution for cyclic, closed, open, tripole solution --- + + cidata_std(:,:,:,:,:) = cidata_bas(:,:,:,:,:) + cjdata_std(:,:,:,:,:) = cjdata_bas(:,:,:,:,:) + + !--- halo off on east and west boundary --- + if (ew_boundary_type == 'closed' .or. & + ew_boundary_type == 'open' ) then + do iblock = 1,numBlocks + call ice_distributionGetBlockID(distrb_info, iblock, blockID) + this_block = get_block(blockID, blockID) + ib = this_block%ilo + ie = this_block%ihi + jb = this_block%jlo + je = this_block%jhi + if (this_block%i_glob(ib) == 1) then + cidata_std(1:ib-1 ,:,:,:,iblock) = dhalofillval + cjdata_std(1:ib-1 ,:,:,:,iblock) = dhalofillval + endif + if (this_block%i_glob(ie) == nx_global) then + cidata_std(ie+1:nx_block,:,:,:,iblock) = dhalofillval + cjdata_std(ie+1:nx_block,:,:,:,iblock) = dhalofillval + endif + enddo + endif + + !--- halo off on south boundary --- + if (ns_boundary_type == 'closed' .or. & + ns_boundary_type == 'open' .or. & + ns_boundary_type == 'tripole' .or. & + ns_boundary_type == 'tripoleT' ) then + do iblock = 1,numBlocks + call ice_distributionGetBlockID(distrb_info, iblock, blockID) + this_block = get_block(blockID, blockID) + ib = this_block%ilo + ie = this_block%ihi + jb = this_block%jlo + je = this_block%jhi + if (this_block%j_glob(jb) == 1) then + cidata_std(:,1:jb-1,:,:,iblock) = dhalofillval + cjdata_std(:,1:jb-1,:,:,iblock) = dhalofillval + endif + enddo + endif + + !--- halo off on north boundary, tripole handled later --- + if (ns_boundary_type == 'closed' .or. & + ns_boundary_type == 'open' .or. & + ns_boundary_type == 'tripole' .or. & + ns_boundary_type == 'tripoleT' ) then + do iblock = 1,numBlocks + call ice_distributionGetBlockID(distrb_info, iblock, blockID) + this_block = get_block(blockID, blockID) + ib = this_block%ilo + ie = this_block%ihi + jb = this_block%jlo + je = this_block%jhi + if (this_block%j_glob(je) == ny_global) then + cidata_std(:,je+1:ny_block,:,:,iblock) = dhalofillval + cjdata_std(:,je+1:ny_block,:,:,iblock) = dhalofillval + endif + enddo + endif + + !--------------------------------------------------------------- + + testcnt = 0 + do nn = 1, maxtests + do nl = 1, maxlocs + do nt = 1, maxtypes + + !--- setup test --- + first_call = .true. + testcnt = testcnt + 1 + if (testcnt > tottest) then + if (my_task == master_task) then + write(6,*) ' ' + write(6,*) 'HALOCHK FAILED' + write(6,*) ' ' + endif + call abort_ice(subname//' testcnt > tottest',file=__FILE__,line=__LINE__) + endif + + !--- fill arrays --- + darrayi1(:,:,:) = fillval + darrayj1(:,:,:) = fillval + darrayi2(:,:,:,:) = fillval + darrayj2(:,:,:,:) = fillval + darrayi3(:,:,:,:,:) = fillval + darrayj3(:,:,:,:,:) = fillval + darrayi1str(:,:,:) = fillval + darrayj1str(:,:,:) = fillval + do iblock = 1,numBlocks + call ice_distributionGetBlockID(distrb_info, iblock, blockID) + this_block = get_block(blockID, blockID) + ib = this_block%ilo + ie = this_block%ihi + jb = this_block%jlo + je = this_block%jhi + do j = jb,je + do i = ib,ie + k1 = 1 + k2 = 1 + darrayi1(i,j,iblock) = cidata_bas(i,j,k1,k2,iblock) + darrayj1(i,j,iblock) = cjdata_bas(i,j,k1,k2,iblock) + do k1 = 1,nz1 + k2 = 1 + darrayi2(i,j,k1,iblock) = cidata_bas(i,j,k1,k2,iblock) + darrayj2(i,j,k1,iblock) = cjdata_bas(i,j,k1,k2,iblock) + do k2 = 1,nz2 + darrayi3(i,j,k1,k2,iblock) = cidata_bas(i,j,k1,k2,iblock) + darrayj3(i,j,k1,k2,iblock) = cjdata_bas(i,j,k1,k2,iblock) + enddo + enddo + enddo + enddo + enddo + + ! copy original darray1 for "stress" compare + darrayi10 = darrayi1 + darrayj10 = darrayj1 + + + !--- halo update --- + + if (nn == 1) then + k1m = 1 + k2m = 1 + halofld = '2DR8' + call ice_haloUpdate(darrayi1, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) + call ice_haloUpdate(darrayj1, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) + elseif (nn == 2) then + k1m = nz1 + k2m = 1 + halofld = '3DR8' + call ice_haloUpdate(darrayi2, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) + call ice_haloUpdate(darrayj2, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) + elseif (nn == 3) then + k1m = nz1 + k2m = nz2 + halofld = '4DR8' + call ice_haloUpdate(darrayi3, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) + call ice_haloUpdate(darrayj3, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) + elseif (nn == 4) then + k1m = 1 + k2m = 1 + halofld = '2DR4' + rarrayi1 = real(darrayi1,kind=real_kind) + rarrayj1 = real(darrayj1,kind=real_kind) + call ice_haloUpdate(rarrayi1, halo_info, field_loc(nl), field_type(nt), fillvalue=rhalofillval) + call ice_haloUpdate(rarrayj1, halo_info, field_loc(nl), field_type(nt), fillvalue=rhalofillval) + darrayi1 = real(rarrayi1,kind=dbl_kind) + darrayj1 = real(rarrayj1,kind=dbl_kind) + elseif (nn == 5) then + k1m = nz1 + k2m = 1 + halofld = '3DR4' + rarrayi2 = real(darrayi2,kind=real_kind) + rarrayj2 = real(darrayj2,kind=real_kind) + call ice_haloUpdate(rarrayi2, halo_info, field_loc(nl), field_type(nt), fillvalue=rhalofillval) + call ice_haloUpdate(rarrayj2, halo_info, field_loc(nl), field_type(nt), fillvalue=rhalofillval) + darrayi2 = real(rarrayi2,kind=dbl_kind) + darrayj2 = real(rarrayj2,kind=dbl_kind) + elseif (nn == 6) then + k1m = nz1 + k2m = nz2 + halofld = '4DR4' + rarrayi3 = real(darrayi3,kind=real_kind) + rarrayj3 = real(darrayj3,kind=real_kind) + call ice_haloUpdate(rarrayi3, halo_info, field_loc(nl), field_type(nt), fillvalue=rhalofillval) + call ice_haloUpdate(rarrayj3, halo_info, field_loc(nl), field_type(nt), fillvalue=rhalofillval) + darrayi3 = real(rarrayi3,kind=dbl_kind) + darrayj3 = real(rarrayj3,kind=dbl_kind) + elseif (nn == 7) then + k1m = 1 + k2m = 1 + halofld = '2DI4' + iarrayi1 = nint(darrayi1) + iarrayj1 = nint(darrayj1) + call ice_haloUpdate(iarrayi1, halo_info, field_loc(nl), field_type(nt), fillvalue=ihalofillval) + call ice_haloUpdate(iarrayj1, halo_info, field_loc(nl), field_type(nt), fillvalue=ihalofillval) + darrayi1 = real(iarrayi1,kind=dbl_kind) + darrayj1 = real(iarrayj1,kind=dbl_kind) + elseif (nn == 8) then + k1m = nz1 + k2m = 1 + halofld = '3DI4' + iarrayi2 = nint(darrayi2) + iarrayj2 = nint(darrayj2) + call ice_haloUpdate(iarrayi2, halo_info, field_loc(nl), field_type(nt), fillvalue=ihalofillval) + call ice_haloUpdate(iarrayj2, halo_info, field_loc(nl), field_type(nt), fillvalue=ihalofillval) + darrayi2 = real(iarrayi2,kind=dbl_kind) + darrayj2 = real(iarrayj2,kind=dbl_kind) + elseif (nn == 9) then + k1m = nz1 + k2m = nz2 + halofld = '4DI4' + iarrayi3 = nint(darrayi3) + iarrayj3 = nint(darrayj3) + call ice_haloUpdate(iarrayi3, halo_info, field_loc(nl), field_type(nt), fillvalue=ihalofillval) + call ice_haloUpdate(iarrayj3, halo_info, field_loc(nl), field_type(nt), fillvalue=ihalofillval) + darrayi3 = real(iarrayi3,kind=dbl_kind) + darrayj3 = real(iarrayj3,kind=dbl_kind) + elseif (nn == 10) then + k1m = 1 + k2m = 1 + halofld = '2DL1' + larrayi1 = .true. + where (darrayi1 == fillval) larrayi1 = .false. + larrayj1 = .false. + where (darrayj1 == fillval) larrayj1 = .true. + call ice_haloUpdate(larrayi1, halo_info, field_loc(nl), field_type(nt), fillvalue=0) + call ice_haloUpdate(larrayj1, halo_info, field_loc(nl), field_type(nt), fillvalue=1) + darrayi1 = c0 + where (larrayi1) darrayi1 = c1 + darrayj1 = c0 + where (larrayj1) darrayj1 = c1 + elseif (nn == 11) then + k1m = 1 + k2m = 1 + halofld = 'STRESS' + darrayi1str = -darrayi1 ! flip sign for testing + darrayj1str = -darrayj1 + call ice_haloUpdate_stress(darrayi1, darrayi1str, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) + call ice_haloUpdate_stress(darrayj1, darrayj1str, halo_info, field_loc(nl), field_type(nt), fillvalue=dhalofillval) + endif + + write(teststring(testcnt),'(5a10)') trim(halofld),trim(locs_name(nl)),trim(types_name(nt)), & + trim(ew_boundary_type),trim(ns_boundary_type) + + do iblock = 1,numBlocks + call ice_distributionGetBlockID(distrb_info, iblock, blockID) + this_block = get_block(blockID, blockID) + ib = this_block%ilo + ie = this_block%ihi + jb = this_block%jlo + je = this_block%jhi + ! just check non-padded gridcells +! do j = 1,ny_block +! do i = 1,nx_block + do j = jb-nghost, je+nghost + do i = ib-nghost, ie+nghost + do k1 = 1,k1m + do k2 = 1,k2m + tripole_average = .false. + tripole_pole = .false. + spvalL1 = .false. + if (index(halofld,'2D') > 0) then + aichk = darrayi1(i,j,iblock) + ajchk = darrayj1(i,j,iblock) + elseif (index(halofld,'STRESS') > 0) then + aichk = darrayi1(i,j,iblock) + ajchk = darrayj1(i,j,iblock) + elseif (index(halofld,'3D') > 0) then + aichk = darrayi2(i,j,k1,iblock) + ajchk = darrayj2(i,j,k1,iblock) + elseif (index(halofld,'4D') > 0) then + aichk = darrayi3(i,j,k1,k2,iblock) + ajchk = darrayj3(i,j,k1,k2,iblock) + else + if (my_task == master_task) then + write(6,*) ' ' + write(6,*) 'HALOCHK FAILED' + write(6,*) ' ' + endif + call abort_ice(subname//' halofld not matched '//trim(halofld),file=__FILE__,line=__LINE__) + endif + + + if (field_loc (nl) == field_loc_noupdate .or. & + field_type(nt) == field_type_noupdate) then + cichk = cidata_nup(i,j,k1,k2,iblock) + cjchk = cjdata_nup(i,j,k1,k2,iblock) + else + cichk = cidata_std(i,j,k1,k2,iblock) + cjchk = cjdata_std(i,j,k1,k2,iblock) + + if (index(halofld,'STRESS') > 0) then + ! only updates on tripole zipper for tripole grids + ! darrayi10 is copy of darrayi1 before halo call + cichk = darrayi10(i,j,iblock) + cjchk = darrayj10(i,j,iblock) + endif + + !--- tripole on north boundary, need to hardcode --- + !--- tripole and tripoleT slightly different --- + !--- establish special set of points here --- + if ((this_block%j_glob(je) == ny_global) .and. & + ((ns_boundary_type == 'tripole' .and. & + (j > je .or. & + (j == je .and. (field_loc(nl) == field_loc_Nface .or. field_loc(nl) == field_loc_NEcorner)))) .or. & + (ns_boundary_type == 'tripoleT' .and. & + (j >= je)))) then + + ! flip sign for vector/angle + if (field_type(nt) == field_type_vector .or. field_type(nt) == field_type_angle ) then + rsign = -c1 + else + rsign = c1 + endif + + ! for tripole + if (ns_boundary_type == 'tripole') then + + ! compute itrip and jtrip, these are the location where the halo values are defined for i,j + ! for j=je averaging, itrip and jtrip are the 2nd gridpoint associated with averaging + + ! standard center tripole u-fold + itrip = nx_global-this_block%i_glob(i)+1 + jtrip = max(je - (j-je) + 1 , je) + ioffset = 0 + joffset = 0 + + if (field_loc(nl) == field_loc_NEcorner .or. field_loc(nl) == field_loc_Nface) then + ! need j offset + joffset = -1 + if (j == je) then + tripole_average = .true. + endif + endif + + if (field_loc(nl) == field_loc_NEcorner .or. field_loc(nl) == field_loc_Eface) then + ! fold plus cell offset + ioffset = -1 + ! CICE treats j=ny_global tripole edge points incorrectly + ! should do edge wraparound and average + ! CICE does not update those points, assumes it's "land" + if (j == je) then + if (this_block%i_glob(i) == nx_global/2) then + tripole_pole = .true. + elseif (this_block%i_glob(i) == nx_global ) then + tripole_pole = .true. + endif + endif + endif + + ! for tripoleT + elseif (ns_boundary_type == 'tripoleT') then + + ! compute itrip and jtrip, these are the location where the halo values are defined for i,j + ! for j=je averaging, itrip and jtrip are the 2nd gridpoint associated with averaging + + ! standard center tripoleT t-fold + itrip = nx_global-this_block%i_glob(i)+2 + jtrip = je - (j-je) + ioffset = 0 + joffset = 0 + + if (field_loc(nl) == field_loc_NEcorner .or. field_loc(nl) == field_loc_Eface) then + ! fold plus cell offset + ioffset = -1 + endif + + if (field_loc(nl) == field_loc_NEcorner .or. field_loc(nl) == field_loc_Nface) then + ! need j offset + joffset = -1 + endif + + if (field_loc(nl) == field_loc_Center .or. field_loc(nl) == field_loc_Eface) then + if (j == je) then + tripole_average = .true. + endif + endif + + ! center point poles need to be treated special + if (field_loc(nl) == field_loc_Center) then + if (j == je .and. & + (this_block%i_glob(i) == 1 .or. this_block%i_glob(i) == nx_global/2+1)) then + tripole_pole = .true. + endif + endif + + endif + + itrip = mod(itrip + ioffset + nx_global-1,nx_global)+1 + jtrip = jtrip + joffset + + rival = (real(itrip,kind=dbl_kind) + & + real(k1,kind=dbl_kind)*1000._dbl_kind + real(k2,kind=dbl_kind)*10000._dbl_kind) + rjval = (real(this_block%j_glob(jtrip),kind=dbl_kind) + & + real(k1,kind=dbl_kind)*1000._dbl_kind + real(k2,kind=dbl_kind)*10000._dbl_kind) + + if (index(halofld,'STRESS') > 0) then + ! only updates on tripole zipper for tripole grids, not tripoleT + if (tripole_pole) then + ! flip sign due to sign of darrayi1str + ! ends of tripole seam not averaged in CICE + cichk = -rsign * cidata_std(i,j,k1,k2,iblock) + cjchk = -rsign * cjdata_std(i,j,k1,k2,iblock) + else + cichk = -rsign * rival + cjchk = -rsign * rjval + endif + elseif (index(halofld,'L1') > 0 .and. j == je) then + ! force cichk and cjchk to match on tripole average index, calc not well defined + spvalL1 = .true. + cichk = aichk + cjchk = ajchk + elseif (tripole_pole) then + ! ends of tripole seam not averaged in CICE + cichk = rsign * cidata_std(i,j,k1,k2,iblock) + cjchk = rsign * cjdata_std(i,j,k1,k2,iblock) + elseif (tripole_average) then + ! tripole average + cichk = p5 * (cidata_std(i,j,k1,k2,iblock) + rsign * rival) + cjchk = p5 * (cjdata_std(i,j,k1,k2,iblock) + rsign * rjval) + else + ! standard tripole fold + cichk = rsign * rival + cjchk = rsign * rjval + endif + +! if (testcnt == 6 .and. j == 61 .and. i < 3) then +! if (testcnt == 186 .and. j == 61 .and. i<4) then +! if (testcnt == 13 .and. j > 61 .and. (i < 3 .or. i > 89)) then +! if (testcnt == 5 .and. j >= 61 .and. (i < 3 .or. i > 90)) then +! write(100+my_task,'(a,5i6,2l3,f6.2,i6)') 'tcx1 ',i,j,iblock,itrip,jtrip, & +! tripole_average,tripole_pole,rsign,this_block%i_glob(i) +! write(100+my_task,'(a,4f12.2)') 'tcx2 ',cidata_std(i,j,k1,k2,iblock),rival,cichk,aichk +! write(100+my_task,'(a,4f12.2)') 'tcx3 ',cjdata_std(i,j,k1,k2,iblock),rjval,cjchk,ajchk +! endif + endif ! tripole or tripoleT + endif + + if (index(halofld,'I4') > 0) then + cichk = real(nint(cichk),kind=dbl_kind) + cjchk = real(nint(cjchk),kind=dbl_kind) + endif + + if (index(halofld,'L1') > 0 .and. .not.spvalL1) then + if (cichk == dhalofillval .or. cichk == fillval) then + cichk = c0 + else + cichk = c1 + endif + if (cjchk == dhalofillval .or. cjchk == fillval) then + cjchk = c1 + else + cjchk = c0 + endif + endif + + ptcnt(testcnt) = ptcnt(testcnt) + 1 + call chkresults(aichk,cichk,errorflag(testcnt),testcnt,failcnt(testcnt), & + i,j,k1,k2,iblock,first_call,teststring(testcnt),trim(halofld)//'_I') + call chkresults(ajchk,cjchk,errorflag(testcnt),testcnt,failcnt(testcnt), & + i,j,k1,k2,iblock,first_call,teststring(testcnt),trim(halofld)//'_J') + enddo ! k2 + enddo ! k1 + enddo ! i + enddo ! j + enddo ! iblock + + enddo ! maxtypes + enddo ! maxlocs + enddo ! maxtests + + ! --------------------------- + ! SUMMARY + ! --------------------------- + + do n = 1,tottest + gflag = global_maxval(errorflag(n), MPI_COMM_ICE) + errorflag(n) = gflag + ptcntsum = global_sum(ptcnt(n),distrb_info) + ptcnt(n) = ptcntsum + failcntsum = global_sum(failcnt(n),distrb_info) + failcnt(n) = failcntsum + enddo + errorflag0 = maxval(errorflag(:)) + + if (my_task == master_task) then + write(6,*) ' ' + write(6,*) 'HALOCHK COMPLETED SUCCESSFULLY' + write(6,*) ' ' + tpcnt = 0 + tfcnt = 0 + do n = 1,tottest + if (errorflag(n) == passflag) then + tpcnt = tpcnt + 1 + write(6,*) 'PASS ',trim(teststring(n)),ptcnt(n),failcnt(n) + else + tfcnt = tfcnt + 1 + write(6,*) 'FAIL ',trim(teststring(n)),ptcnt(n),failcnt(n) + endif + enddo + write(6,*) ' ' + write(6,*) ' total pass = ',tpcnt + write(6,*) ' total fail = ',tfcnt + write(6,*) ' ' + if (errorflag0 == passflag) then + write(6,*) 'HALOCHK TEST COMPLETED SUCCESSFULLY' + else + write(6,*) 'HALOCHK TEST FAILED' + endif + write(6,*) ' ' + write(6,*) '==========================================================' + write(6,*) ' ' + endif + + + !----------------------------------------------------------------- + ! Gracefully end + !----------------------------------------------------------------- + + call end_run() + + end program halochk + +!======================================================================= + + subroutine chkresults(a1,r1,errorflag,testcnt,failcnt,i,j,k1,k2,iblock,first_call,teststring,halofld) + + use halochk_data + + implicit none + + real(dbl_kind) , intent(in) :: a1,r1 + integer(int_kind), intent(inout) :: errorflag, failcnt + integer(int_kind), intent(in) :: i,j,k1,k2,iblock,testcnt + logical , intent(inout) :: first_call + character(len=*) , intent(in) :: teststring,halofld + + logical,parameter :: print_always = .false. + character(len=*) , parameter :: subname='(chkresults)' + + if (a1 /= r1 .or. print_always) then + errorflag = failflag + failcnt = failcnt + 1 + if (first_call) then + write(100+my_task,*) ' ' + write(100+my_task,'(a,i4,2a)') '------- TEST = ',testcnt,' ',trim(teststring) + write(100+my_task,*) ' ' + write(100+my_task,'(a)') ' test task i j k1 k2 iblock expected halocomp diff' + first_call = .false. + endif + write(100+my_task,1001) trim(halofld),testcnt,my_task,i,j,k1,k2,iblock,r1,a1,r1-a1 + endif + + 1001 format(a8,7i6,3f12.3) + + end subroutine chkresults +!======================================================================= diff --git a/cicecore/shared/ice_constants.F90 b/cicecore/shared/ice_constants.F90 index f2da2ef9d..6656213be 100644 --- a/cicecore/shared/ice_constants.F90 +++ b/cicecore/shared/ice_constants.F90 @@ -97,8 +97,7 @@ module ice_constants field_loc_center = 1, & field_loc_NEcorner = 2, & field_loc_Nface = 3, & - field_loc_Eface = 4, & - field_loc_Wface = 5 + field_loc_Eface = 4 !----------------------------------------------------------------- ! field type attribute - necessary for handling diff --git a/configuration/scripts/Makefile b/configuration/scripts/Makefile index a2f17256f..872f426ad 100644 --- a/configuration/scripts/Makefile +++ b/configuration/scripts/Makefile @@ -74,7 +74,7 @@ AR := ar .SUFFIXES: -.PHONY: all cice libcice targets target db_files db_flags clean realclean helloworld calchk sumchk bcstchk gridavgchk optargs +.PHONY: all cice libcice targets target db_files db_flags clean realclean helloworld calchk sumchk bcstchk gridavgchk halochk optargs all: $(EXEC) cice: $(EXEC) @@ -93,7 +93,7 @@ targets: @echo " " @echo "Supported Makefile Targets are: cice, libcice, makdep, depends, clean, realclean" @echo " Diagnostics: targets, db_files, db_flags" - @echo " Unit Tests : helloworld, calchk, sumchk, bcstchk, gridavgchk, optargs" + @echo " Unit Tests : helloworld, calchk, sumchk, bcstchk, gridavgchk, halochk, optargs" target: targets db_files: @@ -151,6 +151,8 @@ bcstchk: $(EXEC) gridavgchk: $(EXEC) +halochk: $(EXEC) + # this builds just a subset of source code specified explicitly and requires a separate target HWOBJS := helloworld.o diff --git a/configuration/scripts/cice_decomp.csh b/configuration/scripts/cice_decomp.csh index 0c6715f3b..bcf27beee 100755 --- a/configuration/scripts/cice_decomp.csh +++ b/configuration/scripts/cice_decomp.csh @@ -156,8 +156,10 @@ if (${ICE_DECOMP_MXBLCKS} > 0) set mxblcks = ${ICE_DECOMP_MXBLCKS} set decomp = 'cartesian' set dshape = 'slenderX2' -if (${nxglob} % ${cicepes} != 0) set decomp = 'roundrobin' -if (${mxblcks} * ${blcky} * 2 < ${nyglob}) set decomp = 'roundrobin' +if (${cicepes} % 2 != 0) set decomp = 'roundrobin' +if (${nyglob} % (${blcky} * 2) != 0) set decomp = 'roundrobin' +if (${nxglob} % ${blckx} != 0) set decomp = 'roundrobin' +if (((${nxglob} * 2) % (${cicepes} * ${blckx})) != 0) set decomp = 'roundrobin' #--- outputs --- diff --git a/configuration/scripts/options/set_env.halochk b/configuration/scripts/options/set_env.halochk new file mode 100644 index 000000000..d09166d2f --- /dev/null +++ b/configuration/scripts/options/set_env.halochk @@ -0,0 +1,2 @@ +setenv ICE_DRVOPT unittest/halochk +setenv ICE_TARGET halochk diff --git a/configuration/scripts/options/set_nml.cyclic b/configuration/scripts/options/set_nml.cyclic new file mode 100644 index 000000000..3a5ae1a7b --- /dev/null +++ b/configuration/scripts/options/set_nml.cyclic @@ -0,0 +1,3 @@ +ew_boundary_type = 'cyclic' +ns_boundary_type = 'cyclic' + diff --git a/configuration/scripts/options/set_nml.open b/configuration/scripts/options/set_nml.open new file mode 100644 index 000000000..0e2d5f388 --- /dev/null +++ b/configuration/scripts/options/set_nml.open @@ -0,0 +1,3 @@ +ew_boundary_type = 'open' +ns_boundary_type = 'open' + diff --git a/configuration/scripts/options/set_nml.tripole b/configuration/scripts/options/set_nml.tripole new file mode 100644 index 000000000..7904b8134 --- /dev/null +++ b/configuration/scripts/options/set_nml.tripole @@ -0,0 +1,3 @@ +grid_type = 'tripole' +ew_boundary_type = 'cyclic' +ns_boundary_type = 'tripole' diff --git a/configuration/scripts/options/set_nml.tripolet b/configuration/scripts/options/set_nml.tripolet new file mode 100644 index 000000000..4bb63dc17 --- /dev/null +++ b/configuration/scripts/options/set_nml.tripolet @@ -0,0 +1,3 @@ +grid_type = 'tripole' +ew_boundary_type = 'cyclic' +ns_boundary_type = 'tripoleT' diff --git a/configuration/scripts/tests/unittest_suite.ts b/configuration/scripts/tests/unittest_suite.ts index 319c91aa6..7486e87aa 100644 --- a/configuration/scripts/tests/unittest_suite.ts +++ b/configuration/scripts/tests/unittest_suite.ts @@ -1,14 +1,29 @@ # Test Grid PEs Sets BFB-compare -unittest gx3 1x1 helloworld -unittest gx3 1x1 optargs -unittest gx3 1x1 calchk,short -unittest gx3 4x1x25x29x4 sumchk -unittest gx3 1x1x25x29x16 sumchk -unittest tx1 8x1 sumchk -unittest gx3 4x1 bcstchk -unittest gx3 1x1 bcstchk -unittest gx3 8x2 gridavgchk,dwblockall -unittest gx3 12x1 gridavgchk -unittest gx1 28x1 gridavgchk,dwblockall -unittest gx1 16x2 gridavgchk -unittest gbox128 8x2 gridavgchk +unittest gx3 1x1 helloworld +unittest gx3 1x1 optargs +unittest gx3 1x1 calchk,short +unittest gx3 4x1x25x29x4 sumchk +unittest gx3 1x1x25x29x16 sumchk +unittest tx1 8x1 sumchk +unittest gx3 4x1 bcstchk +unittest gx3 1x1 bcstchk +unittest gx3 8x2 gridavgchk,dwblockall +unittest gx3 12x1 gridavgchk +unittest gx1 28x1 gridavgchk,dwblockall +unittest gx1 16x2 gridavgchk +unittest gbox128 8x2 gridavgchk +unittest gbox80 1x1x10x10x80 halochk,cyclic,debug +unittest gbox80 1x1x24x23x16 halochk +unittest gbox80 1x1x23x24x16 halochk,cyclic +unittest gbox80 1x1x23x23x16 halochk,open +unittest tx1 1x1x90x60x16 halochk,dwblockall +unittest tx1 1x1x90x60x16 halochk,dwblockall,tripolet +unittest tx1 1x1x95x65x16 halochk,dwblockall +unittest tx1 1x1x95x65x16 halochk,dwblockall,tripolet +unittest gx3 4x2 halochk,dwblockall,debug +unittest gx3 8x2x16x12x10 halochk,cyclic,dwblockall +unittest gx3 17x1x16x12x10 halochk,open,dwblockall +unittest tx1 4x2 halochk,dwblockall +unittest tx1 4x2 halochk,dwblockall,tripolet +unittest tx1 4x2x65x45x10 halochk,dwblockall +unittest tx1 4x2x57x43x12 halochk,dwblockall,tripolet diff --git a/doc/source/cice_index.rst b/doc/source/cice_index.rst index 0e9d21517..000004bb9 100644 --- a/doc/source/cice_index.rst +++ b/doc/source/cice_index.rst @@ -620,7 +620,8 @@ either Celsius or Kelvin units). Deprecated parameters are listed at the end. "shear", "strain rate II component", "1/s" "shlat", "southern latitude of artificial mask edge", "30\ :math:`^\circ`\ N" "shortwave", "flag for shortwave parameterization (‘ccsm3’ or ‘dEdd’)", "" - "sig1(2)", "principal stress components (diagnostic)", "" + "sig1(2)", "principal stress components :math:`\sigma_{n,1}`, :math:`\sigma_{n,2}` (diagnostic)", "" + "sigP", "internal ice pressure", "N/m" "sil", "silicate concentration", "mmol/m\ :math:`^3`" "sinw", "sine of the turning angle in water", "0." "Sinz", "ice salinity profile", "ppt" @@ -660,8 +661,8 @@ either Celsius or Kelvin units). Deprecated parameters are listed at the end. "strax(y)", "wind stress components from data", "N/m\ :math:`^2`" "strength", "ice strength", "N/m" "stress12", "internal ice stress, :math:`\sigma_{12}`", "N/m" - "stressm", "internal ice stress, :math:`\sigma_{11}-\sigma_{22}`", "N/m" - "stressp", "internal ice stress, :math:`\sigma_{11}+\sigma_{22}`", "N/m" + "stressm", "internal ice stress, :math:`\sigma_{11}-\sigma_{22}` (:math:`\sigma_2` in the doc)", "N/m" + "stressp", "internal ice stress, :math:`\sigma_{11}+\sigma_{22}` (:math:`\sigma_1` in the doc)", "N/m" "strintx(y)U", "divergence of internal ice stress, x(y)", "N/m\ :math:`^2`" "strocnx(y)U", "ice–ocean stress in the x(y)-direction (U-cell)", "N/m\ :math:`^2`" "strocnx(y)T", "ice–ocean stress, x(y)-dir. (T-cell)", "N/m\ :math:`^2`" diff --git a/doc/source/conf.py b/doc/source/conf.py index 88b98bc09..7d79f7b43 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -38,6 +38,9 @@ 'sphinxcontrib.bibtex', ] +# Name of the bibliography file for sphinxcontrib.bibtex. +bibtex_bibfiles = ['master_list.bib'] + # Add any paths that contain templates here, relative to this directory. templates_path = ['_templates'] diff --git a/doc/source/science_guide/sg_dynamics.rst b/doc/source/science_guide/sg_dynamics.rst index 6b269d453..585c18616 100644 --- a/doc/source/science_guide/sg_dynamics.rst +++ b/doc/source/science_guide/sg_dynamics.rst @@ -237,15 +237,15 @@ In the VP approach, equation :eq:`momsys` is discretized implicitly using a Back and stresses are not computed explicitly: .. math:: - \begin{align} + \begin{aligned} m\frac{(u^{n}-u^{n-1})}{\Delta t} &= \frac{\partial \sigma_{1j}^n}{\partial x_j} - \tau_{w,x}^n + \tau_{b,x}^n + mfv^n + r_{x}^n, \\ m\frac{(v^{n}-v^{n-1})}{\Delta t} &= \frac{\partial \sigma^{n} _{2j}}{\partial x_j} - - \tau_{w,y}^n + \tau_{b,y}^n -mfu^{n} + - \tau_{w,y}^n + \tau_{b,y}^n - mfu^{n} + r_{y}^n - \end{align} + \end{aligned} :label: u_sit where :math:`r = (r_x,r_y)` contains all terms that do not depend on the velocities :math:`u^n, v^n` (namely the sea surface tilt and the wind stress). @@ -297,7 +297,7 @@ The Hibler-Bryan form for the ice-ocean stress :cite:`Hibler87` is included in **ice\_dyn\_shared.F90** but is currently commented out, pending further testing. -.. _seabed-stress: +.. _seabedstress: Seabed stress ~~~~~~~~~~~~~ @@ -449,8 +449,8 @@ Rheology ******** For convenience we formulate the stress tensor :math:`\bf \sigma` in -terms of :math:`\sigma_1=\sigma_{11}+\sigma_{22}`, -:math:`\sigma_2=\sigma_{11}-\sigma_{22}`, and introduce the +terms of :math:`\sigma_1=\sigma_{11}+\sigma_{22}` (``stressp``), +:math:`\sigma_2=\sigma_{11}-\sigma_{22}` (``stressm``), and introduce the divergence, :math:`D_D`, and the horizontal tension and shearing strain rates, :math:`D_T` and :math:`D_S` respectively: @@ -468,8 +468,16 @@ where .. math:: \dot{\epsilon}_{ij} = {1\over 2}\left({{\partial u_i}\over{\partial x_j}} + {{\partial u_j}\over{\partial x_i}}\right) -CICE can output the internal ice pressure which is an important field to support navigation in ice-infested water. -The internal ice pressure (``sigP``) is the average of the normal stresses multiplied by :math:`-1` and +Note that :math:`\sigma_1` and :math:`\sigma_2` are not to be confused with the normalized principal stresses, +:math:`\sigma_{n,1}` and :math:`\sigma_{n,2}` (``sig1`` and ``sig2``), which are defined as: + +.. math:: + \sigma_{n,1}, \sigma_{n,2} = \frac{1}{P} \left( \frac{\sigma_1}{2} \pm \sqrt{\left(\frac{\sigma_2}{2}\right)^2 + \sigma_{12}^2} \right) + +where :math:`P` is the ice strength. + +In addition to the normalized principal stresses, CICE can output the internal ice pressure which is an important field to support navigation in ice-infested water. +The internal ice pressure (``sigP``) is the average of the normal stresses (:math:`\sigma_{11}`, :math:`\sigma_{22}`) multiplied by :math:`-1` and is therefore simply equal to :math:`-\sigma_1/2`. .. _stress-vp: diff --git a/doc/source/user_guide/ug_implementation.rst b/doc/source/user_guide/ug_implementation.rst index 5ed2092c0..d9ea07a02 100644 --- a/doc/source/user_guide/ug_implementation.rst +++ b/doc/source/user_guide/ug_implementation.rst @@ -275,32 +275,74 @@ the namelist variable ``ns_boundary_type``, ‘tripole’ for the U-fold and ‘tripoleT’ for the T-fold grid. In the U-fold tripole grid, the poles have U-index -:math:`{\tt nx\_global}/2` and ``nx_global`` on the top U-row of the -physical grid, and points with U-index i and :math:`{\tt nx\_global-i}` +:math:`nx\_global/2` and :math:`nx\_global` on the top U-row of the +physical grid, and points with U-index :math:`i` and :math:`nx\_global-i` are coincident. Let the fold have U-row index :math:`n` on the global grid; this will also be the T-row index of the T-row to the south of the fold. There are ghost (halo) T- and U-rows to the north, beyond the fold, on the logical grid. The point with index i along the ghost T-row of index :math:`n+1` physically coincides with point -:math:`{\tt nx\_global}-{\tt i}+1` on the T-row of index :math:`n`. The +:math:`nx\_global-i+1` on the T-row of index :math:`n`. The ghost U-row of index :math:`n+1` physically coincides with the U-row of -index :math:`n-1`. +index :math:`n-1`. In the schematics below, symbols A-H represent +grid points from 1:nx_global at a given j index and the setup of the +tripole seam is depicted within a few rows of the seam. -In the T-fold tripole grid, the poles have T-index 1 and and -:math:`{\tt nx\_global}/2+1` on the top T-row of the physical grid, and -points with T-index i and :math:`{\tt nx\_global}-{\tt i}+2` are +.. _tab-tripole: + +.. table:: Tripole (u-fold) Grid Schematic + :align: center + + +--------------+---------------------------------------+--------------+ + | global j | | global j | + | index | grid point IDs (i index) | index source | + +==============+====+====+====+====+====+====+====+====+==============+ + | ny_global+2 | H | G | F | E | D | C | B | A | ny_global-1 | + +--------------+----+----+----+----+----+----+----+----+--------------+ + | ny_global+1 | H | G | F | E | D | C | B | A | ny_global | + +--------------+----+----+----+----+----+----+----+----+--------------+ + | ny_global | A | B | C | D | E | F | G | H | | + +--------------+----+----+----+----+----+----+----+----+--------------+ + | ny_global-1 | A | B | C | D | E | F | G | H | | + +--------------+----+----+----+----+----+----+----+----+--------------+ + + +In the T-fold tripole grid, the poles have T-index :math:`1` and and +:math:`nx\_global/2+1` on the top T-row of the physical grid, and +points with T-index :math:`i` and :math:`nx\_global-i+2` are coincident. Let the fold have T-row index :math:`n` on the global grid. It is usual for the northernmost row of the physical domain to be a U-row, but in the case of the T-fold, the U-row of index :math:`n` is “beyond” the fold; although it is not a ghost row, it is not physically independent, because it coincides with U-row :math:`n-1`, and it therefore has to be treated like a ghost row. Points i on U-row -:math:`n` coincides with :math:`{\tt nx\_global}-{\tt i}+1` on U-row +:math:`n` coincides with :math:`nx\_global-i+1` on U-row :math:`n-1`. There are still ghost T- and U-rows :math:`n+1` to the north of U-row :math:`n`. Ghost T-row :math:`n+1` coincides with T-row :math:`n-1`, and ghost U-row :math:`n+1` coincides with U-row :math:`n-2`. +.. _tab-tripoleT: + +.. table:: TripoleT (t-fold) Grid Schematic + :align: center + + +--------------+--------------------------------------------+--------------+ + | global j | | global j | + | index | grid point IDs (i index) | index source | + +==============+====+====+====+====+====+====+====+====+====+==============+ + | ny_global+2 | | H | G | F | E | D | C | B | A | ny_global-2 | + +--------------+----+----+----+----+----+----+----+----+----+--------------+ + | ny_global+1 | | H | G | F | E | D | C | B | A | ny_global-1 | + +--------------+----+----+----+----+----+----+----+----+----+--------------+ + | ny_global | A | BH | CG | DF | E | FD | GC | HB | | | + +--------------+----+----+----+----+----+----+----+----+----+--------------+ + | ny_global-1 | A | B | C | D | E | F | G | H | | | + +--------------+----+----+----+----+----+----+----+----+----+--------------+ + | ny_global-2 | A | B | C | D | E | F | G | H | | | + +--------------+----+----+----+----+----+----+----+----+----+--------------+ + + The tripole grid thus requires two special kinds of treatment for certain rows, arranged by the halo-update routines. First, within rows along the fold, coincident points must always have the same value. This @@ -310,7 +352,8 @@ the coincident physical rows. Both operations involve the tripole buffer, which is used to assemble the data for the affected rows. Special treatment is also required in the scattering routine, and when computing global sums one of each pair of coincident points has to be -excluded. +excluded. Halos of center, east, north, and northeast points are supported, +and each requires slightly different halo indexing across the tripole seam. ***************** Rectangular grids @@ -1207,7 +1250,7 @@ directory in **iceh\_ic.[timeID].nc(da)**. Several history variables are hard-coded for instantaneous output regardless of the ``hist_avg`` averaging flag, at the frequency given by their namelist flag. -The normalized principal components of internal ice stress are computed +The normalized principal components of internal ice stress (``sig1``, ``sig2``) are computed in *principal\_stress* and written to the history file. This calculation is not necessary for the simulation; principal stresses are merely computed for diagnostic purposes and included here for the user’s diff --git a/doc/source/user_guide/ug_testing.rst b/doc/source/user_guide/ug_testing.rst index 289f626a9..606ae1397 100644 --- a/doc/source/user_guide/ug_testing.rst +++ b/doc/source/user_guide/ug_testing.rst @@ -736,6 +736,7 @@ The following are brief descriptions of some of the current unit tests, - **calchk** is a unit test that exercises the CICE calendar over 100,000 years and verifies correctness. This test does not depend on the CICE initialization. - **gridavgchk** is a unit test that exercises the CICE grid_average_X2Y methods and verifies results. + - **halochk** is a unit test that exercises the CICE haloUpdate methods and verifies results. - **helloworld** is a simple test that writes out helloworld and uses no CICE infrastructure. This tests exists to demonstrate how to build a unit test by specifying the object files directly in the Makefile diff --git a/icepack b/icepack index 8f96707a9..a4779cc71 160000 --- a/icepack +++ b/icepack @@ -1 +1 @@ -Subproject commit 8f96707a90132ca119d81ed84e5a62ca0ff3ed96 +Subproject commit a4779cc71125b40a7db3a4da8512247cbf2b0955