From 7f4065b50961ddc2b7d4169d1ff741737c148ef7 Mon Sep 17 00:00:00 2001 From: Anton Steketee <79179784+anton-seaice@users.noreply.github.com> Date: Wed, 6 Nov 2024 03:32:05 +1100 Subject: [PATCH] Allow kmt variable name of "mask" and different size netcdf 2d reads (#985) This change allows for a single processor read of netcdf variables that are not of size nx_global by ny_global. This is to make it easier to extent the ice_grid routines. This changes allows the mask variable in the kmt netcdf file to be called "mask" instead of "kmc". Both options now work. This change splits reading the mask and grid into seperate subroutines, so the routine can be reused independently later. * Allow for varname=mask and reading different size nc files * Update ug_implementation.rst * fix gadi test submission * Tidy gadi compiler flags --- cicecore/cicedyn/infrastructure/ice_grid.F90 | 220 +++++++++++------- .../cicedyn/infrastructure/ice_read_write.F90 | 42 +++- configuration/scripts/cice.batch.csh | 11 + .../scripts/machines/Macros.gadi_intel | 4 +- configuration/scripts/tests/base_suite.ts | 18 +- doc/source/user_guide/ug_implementation.rst | 25 +- 6 files changed, 213 insertions(+), 107 deletions(-) diff --git a/cicecore/cicedyn/infrastructure/ice_grid.F90 b/cicecore/cicedyn/infrastructure/ice_grid.F90 index 54bc3ad92..de0fbb586 100644 --- a/cicecore/cicedyn/infrastructure/ice_grid.F90 +++ b/cicecore/cicedyn/infrastructure/ice_grid.F90 @@ -169,6 +169,10 @@ module ice_grid logical (kind=log_kind), private :: & l_readCenter ! If anglet exist in grid file read it otherwise calculate it + character (len=char_len), private :: & + mask_fieldname !field/var name for the mask variable (in nc files) + + interface grid_average_X2Y module procedure grid_average_X2Y_base , & grid_average_X2Y_userwghts, & @@ -291,6 +295,11 @@ end subroutine dealloc_grid subroutine init_grid1 +#ifdef USE_NETCDF + use netcdf, only: nf90_inq_varid , nf90_noerr + integer (kind=int_kind) :: status, varid +#endif + integer (kind=int_kind) :: & fid_grid, & ! file id for netCDF grid file fid_kmt ! file id for netCDF kmt file @@ -342,19 +351,31 @@ subroutine init_grid1 trim(grid_type) == 'regional' ) then if (trim(grid_format) == 'nc') then - - call ice_open_nc(grid_file,fid_grid) - call ice_open_nc(kmt_file,fid_kmt) - + fieldname='ulat' + call ice_open_nc(grid_file,fid_grid) call ice_read_global_nc(fid_grid,1,fieldname,work_g1,.true.) - fieldname='kmt' - call ice_read_global_nc(fid_kmt,1,fieldname,work_g2,.true.) - - if (my_task == master_task) then - call ice_close_nc(fid_grid) - call ice_close_nc(fid_kmt) + call ice_close_nc(fid_grid) + + ! mask variable name might be kmt or mask, check both + call ice_open_nc(kmt_file,fid_kmt) +#ifdef USE_NETCDF + if ( my_task==master_task ) then + status = nf90_inq_varid(fid_kmt, 'kmt', varid) + if (status == nf90_noerr) then + mask_fieldname = 'kmt' + else + status = nf90_inq_varid(fid_kmt, 'mask', varid) + call ice_check_nc(status, subname//' ERROR: does '//trim(kmt_file)//& + ' contain "kmt" or "mask" variable?', file=__FILE__, line=__LINE__) + mask_fieldname = 'mask' + endif endif +#endif + call broadcast_scalar(mask_fieldname, master_task) + + call ice_read_global_nc(fid_kmt,1,mask_fieldname,work_g2,.true.) + call ice_close_nc(fid_kmt) else @@ -466,8 +487,10 @@ subroutine init_grid2 trim(grid_type) == 'tripole' .or. & trim(grid_type) == 'regional' ) then if (trim(grid_format) == 'nc') then + call kmtmask_nc ! read mask from nc file call popgrid_nc ! read POP grid lengths from nc file else + call kmtmask ! read kmt directly call popgrid ! read POP grid lengths directly endif #ifdef CESMCOUPLED @@ -607,6 +630,16 @@ subroutine init_grid2 file=__FILE__, line=__LINE__) endif + if (l_readCenter) then + out_of_range = .false. + where (ANGLET < -pi .or. ANGLET > pi) out_of_range = .true. + if (count(out_of_range) > 0) then + write(nu_diag,*) subname,' angle = ',minval(ANGLET),maxval(ANGLET),count(out_of_range) + call abort_ice (subname//' ANGLET out of expected range', & + file=__FILE__, line=__LINE__) + endif + endif + !----------------------------------------------------------------- ! Compute ANGLE on T-grid !----------------------------------------------------------------- @@ -722,21 +755,10 @@ end subroutine init_grid2 !======================================================================= -! POP displaced pole grid and land mask (or tripole). -! Grid record number, field and units are: \\ -! (1) ULAT (radians) \\ -! (2) ULON (radians) \\ -! (3) HTN (cm) \\ -! (4) HTE (cm) \\ -! (5) HUS (cm) \\ -! (6) HUW (cm) \\ -! (7) ANGLE (radians) -! +! POP land mask ! Land mask record number and field is (1) KMT. -! -! author: Elizabeth C. Hunke, LANL - subroutine popgrid + subroutine kmtmask integer (kind=int_kind) :: & i, j, iblk, & @@ -744,18 +766,14 @@ subroutine popgrid logical (kind=log_kind) :: diag - real (kind=dbl_kind), dimension(:,:), allocatable :: & - work_g1 - real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: & work1 type (block) :: & this_block ! block information for current block - character(len=*), parameter :: subname = '(popgrid)' + character(len=*), parameter :: subname = '(kmtmask)' - call ice_open(nu_grid,grid_file,64) call ice_open(nu_kmt,kmt_file,32) diag = .true. ! write diagnostic info @@ -763,13 +781,13 @@ subroutine popgrid !----------------------------------------------------------------- ! topography !----------------------------------------------------------------- + kmt(:,:,:) = c0 + hm (:,:,:) = c0 - call ice_read(nu_kmt,1,work1,'ida4',diag, & + call ice_read(nu_kmt,1,kmt,'ida4',diag, & field_loc=field_loc_center, & field_type=field_type_scalar) - hm (:,:,:) = c0 - kmt(:,:,:) = c0 !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block) do iblk = 1, nblocks this_block = get_block(blocks_ice(iblk),iblk) @@ -780,12 +798,44 @@ subroutine popgrid do j = jlo, jhi do i = ilo, ihi - kmt(i,j,iblk) = work1(i,j,iblk) if (kmt(i,j,iblk) >= p5) hm(i,j,iblk) = c1 enddo enddo enddo !$OMP END PARALLEL DO + + if (my_task == master_task) then + close (nu_kmt) + endif + + end subroutine kmtmask + +!======================================================================= + +! POP displaced pole grid (or tripole). +! Grid record number, field and units are: \\ +! (1) ULAT (radians) \\ +! (2) ULON (radians) \\ +! (3) HTN (cm) \\ +! (4) HTE (cm) \\ +! (5) HUS (cm) \\ +! (6) HUW (cm) \\ +! (7) ANGLE (radians) +! +! author: Elizabeth C. Hunke, LANL + + subroutine popgrid + + logical (kind=log_kind) :: diag + + real (kind=dbl_kind), dimension(:,:), allocatable :: & + work_g1 + + character(len=*), parameter :: subname = '(popgrid)' + + call ice_open(nu_grid,grid_file,64) + + diag = .true. ! write diagnostic info !----------------------------------------------------------------- ! lat, lon, angle @@ -827,13 +877,65 @@ subroutine popgrid if (my_task == master_task) then close (nu_grid) - close (nu_kmt) endif end subroutine popgrid !======================================================================= +! POP/MOM land mask. +! Land mask field is kmt or mask, saved in mask_fieldname. + + subroutine kmtmask_nc + + integer (kind=int_kind) :: & + i, j, iblk, & + ilo,ihi,jlo,jhi, & ! beginning and end of physical domain + fid_kmt ! file id for netCDF kmt file + + logical (kind=log_kind) :: diag + + real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: & + work1 + + type (block) :: & + this_block ! block information for current block + + character(len=*), parameter :: subname = '(kmtmask_nc)' + + diag = .true. ! write diagnostic info + + hm (:,:,:) = c0 + kmt(:,:,:) = c0 + + call ice_open_nc(kmt_file,fid_kmt) + + call ice_read_nc(fid_kmt,1,mask_fieldname,kmt,diag, & + field_loc=field_loc_center, & + field_type=field_type_scalar) + + call ice_close_nc(fid_kmt) + + !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block) + do iblk = 1, nblocks + this_block = get_block(blocks_ice(iblk),iblk) + ilo = this_block%ilo + ihi = this_block%ihi + jlo = this_block%jlo + jhi = this_block%jhi + + do j = jlo, jhi + do i = ilo, ihi + if (kmt(i,j,iblk) >= c1) hm(i,j,iblk) = c1 + enddo + enddo + enddo + !$OMP END PARALLEL DO + + end subroutine kmtmask_nc + +!======================================================================= + ! POP displaced pole grid and land mask. ! Grid record number, field and units are: \\ ! (1) ULAT (radians) \\ @@ -844,8 +946,6 @@ end subroutine popgrid ! (6) HUW (cm) \\ ! (7) ANGLE (radians) ! -! Land mask record number and field is (1) KMT. -! ! author: Elizabeth C. Hunke, LANL ! Revised for netcdf input: Ann Keen, Met Office, May 2007 @@ -858,8 +958,7 @@ subroutine popgrid_nc integer (kind=int_kind) :: & i, j, iblk, & ilo,ihi,jlo,jhi, & ! beginning and end of physical domain - fid_grid, & ! file id for netCDF grid file - fid_kmt ! file id for netCDF kmt file + fid_grid ! file id for netCDF grid file logical (kind=log_kind) :: diag @@ -872,17 +971,8 @@ subroutine popgrid_nc real (kind=dbl_kind), dimension(:,:), allocatable :: & work_g1 - real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: & - work1 - - type (block) :: & - this_block ! block information for current block - integer(kind=int_kind) :: & - varid - integer (kind=int_kind) :: & - status ! status flag - + varid, status character(len=*), parameter :: subname = '(popgrid_nc)' @@ -893,37 +983,9 @@ subroutine popgrid_nc file=__FILE__, line=__LINE__) call ice_open_nc(grid_file,fid_grid) - call ice_open_nc(kmt_file,fid_kmt) diag = .true. ! write diagnostic info - !----------------------------------------------------------------- - ! topography - !----------------------------------------------------------------- - - fieldname='kmt' - call ice_read_nc(fid_kmt,1,fieldname,work1,diag, & - field_loc=field_loc_center, & - field_type=field_type_scalar) - - hm (:,:,:) = c0 - kmt(:,:,:) = c0 - !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block) - do iblk = 1, nblocks - this_block = get_block(blocks_ice(iblk),iblk) - ilo = this_block%ilo - ihi = this_block%ihi - jlo = this_block%jlo - jhi = this_block%jhi - - do j = jlo, jhi - do i = ilo, ihi - kmt(i,j,iblk) = work1(i,j,iblk) - if (kmt(i,j,iblk) >= c1) hm(i,j,iblk) = c1 - enddo - enddo - enddo - !$OMP END PARALLEL DO - + !----------------------------------------------------------------- ! lat, lon, angle !----------------------------------------------------------------- @@ -996,10 +1058,8 @@ subroutine popgrid_nc deallocate(work_g1) - if (my_task == master_task) then - call ice_close_nc(fid_grid) - call ice_close_nc(fid_kmt) - endif + call ice_close_nc(fid_grid) + #else call abort_ice(subname//' ERROR: USE_NETCDF cpp not defined', & file=__FILE__, line=__LINE__) diff --git a/cicecore/cicedyn/infrastructure/ice_read_write.F90 b/cicecore/cicedyn/infrastructure/ice_read_write.F90 index 4613843b5..b8bafc10b 100644 --- a/cicecore/cicedyn/infrastructure/ice_read_write.F90 +++ b/cicecore/cicedyn/infrastructure/ice_read_write.F90 @@ -2474,7 +2474,7 @@ subroutine ice_read_global_nc (fid, nrec, varname, work_g, diag) character (char_len), intent(in) :: & varname ! field name in netcdf file - real (kind=dbl_kind), dimension(nx_global,ny_global), intent(out) :: & + real (kind=dbl_kind), dimension(:,:), intent(out) :: & work_g ! output array (real, 8-byte) logical (kind=log_kind) :: & @@ -2485,13 +2485,15 @@ subroutine ice_read_global_nc (fid, nrec, varname, work_g, diag) character(len=*), parameter :: subname = '(ice_read_global_nc)' #ifdef USE_NETCDF -! netCDF file diagnostics: + ! netCDF file diagnostics: integer (kind=int_kind) :: & - varid, & ! netcdf id for field - status ! status output from netcdf routines -! ndim, nvar, & ! sizes of netcdf file + varid, & ! netcdf id for field + status, & ! status output from netcdf routines + ndim, & ! number of variable dimensions + dimids(NF90_MAX_VAR_DIMS) , & !ids of dimensions + dimlen ! size of dimension +! nvar, & ! sizes of netcdf file ! id, & ! dimension index -! dimlen ! size of dimension real (kind=dbl_kind) :: & amin, amax, asum ! min, max values and sum of input array @@ -2534,8 +2536,32 @@ subroutine ice_read_global_nc (fid, nrec, varname, work_g, diag) file=__FILE__, line=__LINE__) work_g=work_g3(2:nx_global+1,1:ny_global) else - status = nf90_get_var( fid, varid, work_g, & - start=(/1,1,nrec/), count=(/nx_global,ny_global,1/)) + ! Check var size : is var 2d ? + status = nf90_inquire_variable(fid, varid, ndims=ndim, dimids=dimids) + call ice_check_nc(status, subname//' ERROR: Cannot check variable '//trim(varname), & + file=__FILE__, line=__LINE__) + if ( ndim > 2 ) then + call abort_ice(subname//' ERROR: '//trim(varname)//' cannot have more than 2 dimensions', & + file=__FILE__, line=__LINE__) + endif + ! Is work_g the same size as the variable? + status = nf90_inquire_dimension(fid, dimids(1), len=dimlen) + call ice_check_nc(status, subname//' ERROR: Cannot check variable '//trim(varname), & + file=__FILE__, line=__LINE__) + if ( dimlen /= size(work_g,1) ) then + call abort_ice(subname//' ERROR: x dim of '//trim(varname)//' wrong size, check nx_global', & + file=__FILE__, line=__LINE__) + endif + status = nf90_inquire_dimension(fid, dimids(2), len=dimlen) + call ice_check_nc(status, subname//' ERROR: Cannot check variable '//trim(varname), & + file=__FILE__, line=__LINE__) + if ( dimlen /= size(work_g,2) ) then + call abort_ice(subname//' ERROR: y dim of '//trim(varname)//' wrong size, check ny_global', & + file=__FILE__, line=__LINE__) + endif + + ! Get the data + status = nf90_get_var( fid, varid, work_g, start=(/1,1,nrec/)) call ice_check_nc(status, subname//' ERROR: Cannot get variable '//trim(varname), & file=__FILE__, line=__LINE__) endif diff --git a/configuration/scripts/cice.batch.csh b/configuration/scripts/cice.batch.csh index a68852602..bb2edc17e 100755 --- a/configuration/scripts/cice.batch.csh +++ b/configuration/scripts/cice.batch.csh @@ -56,12 +56,23 @@ EOFB else if (${ICE_MACHINE} =~ gadi*) then if (${queue} =~ *sr) then #sapphire rapids @ memuse = ( $ncores * 481 / 100 ) + set corespernode = 52 else if (${queue} =~ *bw) then #broadwell @ memuse = ( $ncores * 457 / 100 ) + set corespernode = 28 else if (${queue} =~ *sl) then @ memuse = ( $ncores * 6 ) + set corespernode = 32 else #normal queues @ memuse = ( $ncores * 395 / 100 ) + set corespernode = 48 +endif +if (${ncores} > ${corespernode}) then + # need to use a whole number of nodes + @ ncores = ( ( $ncores + $corespernode - 1 ) / $corespernode ) * $corespernode +endif +if (${runlength} <= 0) then + set batchtime = "00:30:00" endif cat >> ${jobfile} << EOFB #PBS -q ${queue} diff --git a/configuration/scripts/machines/Macros.gadi_intel b/configuration/scripts/machines/Macros.gadi_intel index df7746731..3663af20d 100644 --- a/configuration/scripts/machines/Macros.gadi_intel +++ b/configuration/scripts/machines/Macros.gadi_intel @@ -13,11 +13,11 @@ NCI_INTEL_FLAGS := -r8 -i4 -traceback -w -fpe0 -ftz -convert big_endian -assume NCI_REPRO_FLAGS := -fp-model precise -fp-model source -align all ifeq ($(ICE_BLDDEBUG), true) - NCI_DEBUG_FLAGS := -g3 -O0 -debug all -check all -no-vec -assume nobuffered_io + NCI_DEBUG_FLAGS := -O0 -debug all -check all -no-vec -assume nobuffered_io FFLAGS := $(NCI_INTEL_FLAGS) $(NCI_REPRO_FLAGS) $(NCI_DEBUG_FLAGS) CPPDEFS := $(CPPDEFS) -DDEBUG=$(DEBUG) else - NCI_OPTIM_FLAGS := -g3 -O2 -axCORE-AVX2 -debug all -check none -qopt-report=5 -qopt-report-annotate -assume buffered_io + NCI_OPTIM_FLAGS := -O2 -axCORE-AVX2 -debug all -check none -qopt-report=5 -qopt-report-annotate -assume buffered_io FFLAGS := $(NCI_INTEL_FLAGS) $(NCI_REPRO_FLAGS) $(NCI_OPTIM_FLAGS) endif diff --git a/configuration/scripts/tests/base_suite.ts b/configuration/scripts/tests/base_suite.ts index 2ac855a24..f148a8f57 100644 --- a/configuration/scripts/tests/base_suite.ts +++ b/configuration/scripts/tests/base_suite.ts @@ -14,21 +14,21 @@ smoke gx3 1x8 diag1,run5day,evp1d restart gx1 40x4 droundrobin,medium restart tx1 40x4 dsectrobin,medium restart tx1 40x4 dsectrobin,medium,jra55do -restart gx3 4x4 none -restart gx3 4x4 gx3nc -restart gx3 10x4 maskhalo +restart gx3 4x4 short +restart gx3 4x4 gx3nc,short +restart gx3 10x4 maskhalo,short restart gx3 6x2 alt01 restart gx3 8x2 alt02 restart gx3 4x2 alt03 restart gx3 12x2 alt03,maskhalo,droundrobin restart gx3 4x4 alt04 -restart gx3 4x4 alt05 +restart gx3 4x4 alt05,short restart gx3 8x2 alt06 restart gx3 8x3 alt07 restart gx3 16x2 snicar restart gx3 12x2 snicartest restart gx3 8x2 congel -restart gx3 8x3 saltflux +restart gx3 8x3 saltflux,short restart gx3 18x2 debug,maskhalo restart gx3 6x2 alt01,debug,short restart gx3 8x2 alt02,debug,short @@ -48,7 +48,7 @@ restart gbox128 4x2 boxnodyn,short restart gbox128 4x2 boxnodyn,short,debug restart gbox128 2x2 boxadv,short smoke gbox128 2x2 boxadv,short,debug -restart gbox128 4x4 boxrestore,short +restart gbox128 4x4 boxrestore,medium smoke gbox128 4x4 boxrestore,short,debug restart gbox80 1x1 box2001 smoke gbox80 1x1 boxslotcyl @@ -79,13 +79,13 @@ restart gx3 8x2 isotope smoke gx3 4x1 snwitdrdg,snwgrain,icdefault,debug smoke gx3 4x1 snw30percent,icdefault,debug restart gx3 8x2 snwitdrdg,icdefault,snwgrain -restart gx3 4x4 gx3ncarbulk,iobinary -restart gx3 4x4 histall,precision8,cdf64 +restart gx3 4x4 gx3ncarbulk,iobinary,medium +restart gx3 4x4 cdf64,histall,precision8,short smoke gx3 30x1 bgcz,histall smoke gx3 14x2 fsd12,histall smoke gx3 4x1 dynpicard restart gx3 8x2 gx3ncarbulk,debug -restart gx3 4x4 gx3ncarbulk,diag1 +restart gx3 4x4 diag1,gx3ncarbulk,short smoke gx3 4x1 calcdragio restart gx3 4x2 atmbndyconstant restart gx3 4x2 atmbndymixed diff --git a/doc/source/user_guide/ug_implementation.rst b/doc/source/user_guide/ug_implementation.rst index 91909082c..0e03c2f13 100644 --- a/doc/source/user_guide/ug_implementation.rst +++ b/doc/source/user_guide/ug_implementation.rst @@ -129,16 +129,25 @@ This is shown in Figure :ref:`fig-Cgrid`. Schematic of CICE CD-grid. -The user has several ways to initialize the grid: *popgrid* reads grid -lengths and other parameters for a nonuniform grid (including tripole -and regional grids), and *rectgrid* creates a regular rectangular grid. -The input files **global_gx3.grid** and **global_gx3.kmt** contain the +The user has several ways to initialize the grid, which can be read from +files or created internally. The *rectgrid* code creates a regular rectangular +grid (use the namelist option ``grid_type='rectangular'``). The *popgrid* and *popgrid_nc* +code reads grid lengths and other parameters for a nonuniform grid (including tripole +and regional grids). +The input files **grid_gx3.bin** and **kmt_gx3.bin** contain the :math:`\left<3^\circ\right>` POP grid and land mask; -**global_gx1.grid** and **global_gx1.kmt** contain the -:math:`\left<1^\circ\right>` grid and land mask, and **global_tx1.grid** -and **global_tx1.kmt** contain the :math:`\left<1^\circ\right>` POP +**grid_gx1.bin** and **kmt_gx1.bin** contain the +:math:`\left<1^\circ\right>` grid and land mask, and **grid_tx1.bin** +and **kmt_tx1.bin** contain the :math:`\left<1^\circ\right>` POP tripole grid and land mask. These are binary unformatted, direct access, -Big Endian files. +Big Endian files. + +The are also input files in netcdf format for the **gx3** grid, +(**grid_gx3.nc** and **kmt_gx3.nc**) which can serve as a template for defining +other grids. At a minimum the grid file needs to to contain ULAT, ULON, HTN, HTE +and ANGLE variables. If the variables exist, ANGLET, TLON and TLAT will also be +read from a netcdf grid file. The kmt (mask) netcdf file needs a variable named +kmt or mask, set to 0 for land and 1 for ocean. The input grid file for the B-grid and CD-grid is identical. That file contains each cells' HTN, HTE, ULON, ULAT, and kmt value. From those