Skip to content

Commit

Permalink
Allow kmt variable name of "mask" and different size netcdf 2d reads (C…
Browse files Browse the repository at this point in the history 
…ICE-Consortium#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
  • Loading branch information
@anton-seaice
anton-seaice authored Nov 5, 2024
1 parent 3ebdef7 commit 7f4065b
Show file tree
Hide file tree
Showing 6 changed files with 213 additions and 107 deletions.
220 changes: 140 additions & 80 deletions cicecore/cicedyn/infrastructure/ice_grid.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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, &
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
!-----------------------------------------------------------------
Expand Down Expand Up @@ -722,54 +755,39 @@ 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, &
ilo,ihi,jlo,jhi ! beginning and end of physical domain

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

!-----------------------------------------------------------------
! 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)
Expand All @@ -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
Expand Down Expand Up @@ -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) \\
Expand All @@ -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

Expand All @@ -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

Expand All @@ -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)'

Expand All @@ -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
!-----------------------------------------------------------------
Expand Down Expand Up @@ -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__)
Expand Down
Loading

0 comments on commit 7f4065b

Please sign in to comment.