Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Allow kmt variable name of "mask" and different size netcdf 2d reads #985

Open
wants to merge 8 commits into
base: main
Choose a base branch
from
Open

Allow kmt variable name of "mask" and different size netcdf 2d reads #985

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
d0e0595
Allow for varname=mask and reading different size nc files
anton-seaice Oct 18, 2024
f498d2b
clean and remove a little duplication
anton-seaice Oct 18, 2024
d5eeebc
Update ug_implementation.rst
anton-seaice Oct 18, 2024
359a2d1
fix gadi test submission
anton-seaice Oct 22, 2024
5c5ff05
Tidy compiler flags
anton-seaice Oct 22, 2024
1cc7450
Merge branch 'main' into kmt_read_tidy
anton-seaice Oct 22, 2024
a106246
extend a few tests
anton-seaice Oct 23, 2024
58bc70b
extend one more
anton-seaice Oct 23, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
220 changes: 140 additions & 80 deletions cicecore/cicedyn/infrastructure/ice_grid.F90
View file
Open in desktop
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
Loading