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

refactor(netcdf): updates focused on flopy parity #2163

Open
wants to merge 7 commits into
base: develop
Choose a base branch
from
Open

refactor(netcdf): updates focused on flopy parity #2163

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
8bb85f3
fix dis mesh_face_x array
Jan 24, 2025
3e2a886
rotate mesh x and y coordinates
Feb 3, 2025
520dd09
name adjustment for internal global attrs
Feb 3, 2025
07029af
rework structured export proj and coords, support feet lenuni
Feb 15, 2025
fe11b5c
add warning when wkt provided for rotated grid with structured export
Feb 19, 2025
55fa6f8
add filename to warnings
Feb 20, 2025
491df74
make error string clear
Feb 28, 2025
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
6 changes: 3 additions & 3 deletions doc/mf6io/mf6ivar/dfn/utl-ncf.dfn
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ shape lenbigline
reader urword
optional true
longname CRS well-known text (WKT) string
description is the coordinate reference system (CRS) well-known text (WKT) string.
description is the coordinate reference system (CRS) well-known text (WKT) string. Ignored if latitude and longitude griddata arrays have been provided for NETCDF\_STRUCTURED export type.

block options
name deflate
Expand Down Expand Up @@ -96,7 +96,7 @@ shape (ncpl)
optional true
reader readarray
longname cell center latitude
description cell center latitude.
description cell center latitude. Only supported for NETCDF\_STRUCTURED export type.

block griddata
name longitude
Expand All @@ -105,4 +105,4 @@ shape (ncpl)
optional true
reader readarray
longname cell center longitude
description cell center longitude.
description cell center longitude. Only supported for NETCDF\_STRUCTURED export type.
32 changes: 21 additions & 11 deletions src/Utilities/Export/DisNCMesh.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,7 @@ subroutine dis_export_init(this, modelname, modeltype, modelfname, nc_fname, &

! initialize base class
call this%mesh_init(modelname, modeltype, modelfname, nc_fname, disenum, &
nctype, iout)
nctype, this%dis%lenuni, iout)
end subroutine dis_export_init

!> @brief netcdf export dis destroy
Expand Down Expand Up @@ -457,12 +457,13 @@ end subroutine define_dim
!> @brief netcdf export add mesh information
!<
subroutine add_mesh_data(this)
use BaseDisModule, only: dis_transform_xy
class(Mesh2dDisExportType), intent(inout) :: this
integer(I4B) :: cnt, maxvert, m
integer(I4B), dimension(:), allocatable :: verts
real(DP), dimension(:), allocatable :: bnds
integer(I4B) :: i, j
real(DP) :: x, y
real(DP) :: x, y, x_transform, y_transform
real(DP), dimension(:), allocatable :: node_x, node_y
real(DP), dimension(:), allocatable :: cell_x, cell_y

Expand All @@ -485,13 +486,18 @@ subroutine add_mesh_data(this)
cnt = 0
node_x = NF90_FILL_DOUBLE
node_y = NF90_FILL_DOUBLE
y = this%dis%yorigin + sum(this%dis%delc)
y = sum(this%dis%delc)
do j = this%dis%nrow, 0, -1
x = this%dis%xorigin
x = 0
do i = this%dis%ncol, 0, -1
cnt = cnt + 1
node_x(cnt) = x
node_y(cnt) = y
call dis_transform_xy(x, y, &
this%dis%xorigin, &
this%dis%yorigin, &
this%dis%angrot, &
x_transform, y_transform)
node_x(cnt) = x_transform
node_y(cnt) = y_transform
if (i > 0) x = x + this%dis%delr(i)
end do
if (j > 0) y = y - this%dis%delc(j)
Expand All @@ -508,13 +514,17 @@ subroutine add_mesh_data(this)
cell_x = NF90_FILL_DOUBLE
cell_y = NF90_FILL_DOUBLE
do j = 1, this%dis%nrow
x = this%dis%xorigin
y = this%dis%celly(j) + this%dis%yorigin
y = this%dis%celly(j)
do i = 1, this%dis%ncol
cell_x(cnt) = x
cell_y(cnt) = y
x = this%dis%cellx(i)
call dis_transform_xy(x, y, &
this%dis%xorigin, &
this%dis%yorigin, &
this%dis%angrot, &
x_transform, y_transform)
cell_x(cnt) = x_transform
cell_y(cnt) = y_transform
cnt = cnt + 1
x = this%dis%cellx(i) + this%dis%xorigin
end do
end do

Expand Down
89 changes: 59 additions & 30 deletions src/Utilities/Export/DisNCStructured.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ module DisNCStructuredModule

use KindModule, only: DP, I4B, LGP
use ConstantsModule, only: LINELENGTH, LENBIGLINE, LENCOMPONENTNAME, &
LENMEMPATH, LENVARNAME
LENMEMPATH, LENVARNAME, DNODATA, DZERO
use SimVariablesModule, only: errmsg, warnmsg
use SimModule, only: store_error, store_warning, store_error_filename
use MemoryManagerModule, only: mem_setptr
Expand Down Expand Up @@ -75,7 +75,7 @@ module DisNCStructuredModule
procedure :: define_dim
procedure :: define_dependent
procedure :: define_gridmap
procedure :: define_projection
procedure :: define_geocoords
procedure :: add_proj_data
procedure :: add_grid_data
end type DisNCStructuredType
Expand Down Expand Up @@ -127,6 +127,7 @@ subroutine dis_export_init(this, modelname, modeltype, modelfname, nc_fname, &
! initialize base class
call this%NCModelExportType%init(modelname, modeltype, modelfname, nc_fname, &
disenum, nctype, iout)

! update values from input context
if (this%ncf_mempath /= '') then
call mem_set_value(this%chunk_z, 'CHUNK_Z', this%ncf_mempath, found)
Expand All @@ -144,7 +145,7 @@ subroutine dis_export_init(this, modelname, modeltype, modelfname, nc_fname, &
this%chunk_x = -1
write (warnmsg, '(a)') 'Ignoring user provided NetCDF chunking &
&parameters. Define chunk_time, chunk_x, chunk_y and chunk_z input &
&parameters to see an effect.'
&parameters to see an effect in file "'//trim(nc_fname)//'".'
call store_warning(warnmsg)
end if

Expand All @@ -153,9 +154,32 @@ subroutine dis_export_init(this, modelname, modeltype, modelfname, nc_fname, &

if (latsz > 0 .and. lonsz > 0) then
this%latlon = .true.
if (this%wkt /= '') then
write (warnmsg, '(a)') 'Ignoring user provided NetCDF wkt parameter &
&as longitude and latitude arrays have been provided. &
&Applies to file "'//trim(nc_fname)//'".'
call store_warning(warnmsg)
this%wkt = ''
this%gridmap_name = ''
end if
call mem_setptr(this%latitude, 'LATITUDE', this%ncf_mempath)
call mem_setptr(this%longitude, 'LONGITUDE', this%ncf_mempath)
end if

if (this%wkt /= '') then
if (this%dis%angrot /= DZERO) then
write (warnmsg, '(a)') 'WKT parameter set with structured rotated &
&grid. Projected coordinates will have grid local values. &
&Applies to file "'//trim(nc_fname)//'".'
call store_warning(warnmsg)
end if
end if
end if

if (this%dis%lenuni == 1) then
this%lenunits = 'ft'
else
this%lenunits = 'm'
end if

! create the netcdf file
Expand Down Expand Up @@ -191,7 +215,7 @@ subroutine df(this)
! define root group dimensions and coordinate variables
call this%define_dim()
! define grid projection variables
call this%define_projection()
call this%define_geocoords()
if (isim_mode /= MVALIDATE) then
! define the dependent variable
call this%define_dependent()
Expand Down Expand Up @@ -367,7 +391,6 @@ end subroutine export_input_arrays
!> @brief netcdf export package dynamic input with ilayer index variable
!<
subroutine package_step_ilayer(this, export_pkg, ilayer_varname, ilayer)
use ConstantsModule, only: DNODATA, DZERO
use TdisModule, only: kper
use NCModelExportModule, only: ExportPackageType
use DefinitionSelectModule, only: get_param_definition_type
Expand Down Expand Up @@ -515,7 +538,6 @@ end subroutine package_step
!<
subroutine export_layer_3d(this, export_pkg, idt, ilayer_read, ialayer, &
dbl1d, nc_varname, input_attr, iaux)
use ConstantsModule, only: DNODATA, DZERO
use TdisModule, only: kper
use NCModelExportModule, only: ExportPackageType
class(DisNCStructuredType), intent(inout) :: this
Expand Down Expand Up @@ -640,10 +662,10 @@ subroutine add_global_att(this)
call nf_verify(nf90_put_att(this%ncid, NF90_GLOBAL, 'source', &
this%annotation%source), this%nc_fname)
! export type (MODFLOW 6)
call nf_verify(nf90_put_att(this%ncid, NF90_GLOBAL, 'modflow6_grid', &
call nf_verify(nf90_put_att(this%ncid, NF90_GLOBAL, 'modflow_grid', &
this%annotation%grid), this%nc_fname)
! MODFLOW 6 model type
call nf_verify(nf90_put_att(this%ncid, NF90_GLOBAL, 'modflow6_model', &
call nf_verify(nf90_put_att(this%ncid, NF90_GLOBAL, 'modflow_model', &
this%annotation%model), this%nc_fname)
! generation datetime
call nf_verify(nf90_put_att(this%ncid, NF90_GLOBAL, 'history', &
Expand Down Expand Up @@ -707,8 +729,8 @@ subroutine define_dim(this)
this%nc_fname)
call nf_verify(nf90_def_var(this%ncid, 'y', NF90_DOUBLE, this%dim_ids%y, &
this%var_ids%y), this%nc_fname)
call nf_verify(nf90_put_att(this%ncid, this%var_ids%y, 'units', 'm'), &
this%nc_fname)
call nf_verify(nf90_put_att(this%ncid, this%var_ids%y, 'units', &
this%lenunits), this%nc_fname)
call nf_verify(nf90_put_att(this%ncid, this%var_ids%y, 'axis', 'Y'), &
this%nc_fname)
call nf_verify(nf90_put_att(this%ncid, this%var_ids%y, 'standard_name', &
Expand All @@ -730,8 +752,8 @@ subroutine define_dim(this)
this%nc_fname)
call nf_verify(nf90_def_var(this%ncid, 'x', NF90_DOUBLE, this%dim_ids%x, &
this%var_ids%x), this%nc_fname)
call nf_verify(nf90_put_att(this%ncid, this%var_ids%x, 'units', 'm'), &
this%nc_fname)
call nf_verify(nf90_put_att(this%ncid, this%var_ids%x, 'units', &
this%lenunits), this%nc_fname)
call nf_verify(nf90_put_att(this%ncid, this%var_ids%x, 'axis', 'X'), &
this%nc_fname)
call nf_verify(nf90_put_att(this%ncid, this%var_ids%x, 'standard_name', &
Expand Down Expand Up @@ -782,7 +804,7 @@ subroutine define_dependent(this)

! put attr
call nf_verify(nf90_put_att(this%ncid, this%var_ids%dependent, &
'units', 'm'), this%nc_fname)
'units', this%lenunits), this%nc_fname)
call nf_verify(nf90_put_att(this%ncid, this%var_ids%dependent, &
'standard_name', this%annotation%stdname), &
this%nc_fname)
Expand Down Expand Up @@ -816,9 +838,9 @@ end subroutine define_gridmap

!> @brief define grid projection variables
!<
subroutine define_projection(this)
subroutine define_geocoords(this)
class(DisNCStructuredType), intent(inout) :: this
if (this%latlon .and. this%wkt /= '') then
if (this%latlon) then
! lat
call nf_verify(nf90_def_var(this%ncid, 'lat', NF90_DOUBLE, &
(/this%dim_ids%x, this%dim_ids%y/), &
Expand All @@ -841,13 +863,13 @@ subroutine define_projection(this)
call nf_verify(nf90_put_att(this%ncid, this%var_ids%longitude, &
'long_name', 'longitude'), this%nc_fname)
end if
end subroutine define_projection
end subroutine define_geocoords

!> @brief add grid projection data
!<
subroutine add_proj_data(this)
class(DisNCStructuredType), intent(inout) :: this
if (this%latlon .and. this%wkt /= '') then
if (this%latlon) then
! lat
call nf_verify(nf90_put_var(this%ncid, this%var_ids%latitude, &
this%latitude, start=(/1, 1/), &
Expand All @@ -869,16 +891,25 @@ subroutine add_grid_data(this)
integer(I4B) :: ibnd, n !, k, i, j
real(DP), dimension(:, :), pointer, contiguous :: dbl2d
real(DP), dimension(:), allocatable :: x, y
real(DP) :: xoff, yoff

if (this%dis%angrot /= DZERO) then
xoff = DZERO
yoff = DZERO
else
xoff = this%dis%xorigin
yoff = this%dis%yorigin
end if

allocate (x(size(this%dis%cellx)))
allocate (y(size(this%dis%celly)))

do n = 1, size(this%dis%cellx)
x(n) = this%dis%cellx(n) + this%dis%xorigin
x(n) = this%dis%cellx(n) + xoff
end do

do n = 1, size(this%dis%celly)
y(n) = this%dis%celly(n) + this%dis%yorigin
y(n) = this%dis%celly(n) + yoff
end do

call nf_verify(nf90_put_var(this%ncid, this%var_ids%x, x), &
Expand All @@ -897,8 +928,8 @@ subroutine add_grid_data(this)
ibnd = 1
do n = 1, size(this%dis%cellx)
if (ibnd == 1) then
dbl2d(1, ibnd) = this%dis%xorigin
dbl2d(2, ibnd) = this%dis%xorigin + this%dis%delr(ibnd)
dbl2d(1, ibnd) = xoff
dbl2d(2, ibnd) = xoff + this%dis%delr(ibnd)
else
dbl2d(1, ibnd) = dbl2d(1, ibnd - 1) + this%dis%delr(ibnd)
dbl2d(2, ibnd) = dbl2d(2, ibnd - 1) + this%dis%delr(ibnd)
Expand All @@ -914,8 +945,8 @@ subroutine add_grid_data(this)
ibnd = 1
do n = size(this%dis%celly), 1, -1
if (ibnd == 1) then
dbl2d(1, ibnd) = this%dis%yorigin + sum(this%dis%delc) - this%dis%delc(n)
dbl2d(2, ibnd) = this%dis%yorigin + sum(this%dis%delc)
dbl2d(1, ibnd) = yoff + sum(this%dis%delc) - this%dis%delc(n)
dbl2d(2, ibnd) = yoff + sum(this%dis%delc)
else
dbl2d(1, ibnd) = dbl2d(1, ibnd - 1) - this%dis%delc(n)
dbl2d(2, ibnd) = dbl2d(2, ibnd - 1) - this%dis%delc(n)
Expand Down Expand Up @@ -982,15 +1013,13 @@ subroutine ncvar_gridmap(ncid, varid, gridmap_name, latlon, nc_fname)
logical(LGP), intent(in) :: latlon
character(len=*), intent(in) :: nc_fname
if (gridmap_name /= '') then
if (latlon) then
call nf_verify(nf90_put_att(ncid, varid, 'coordinates', 'lon lat'), &
nc_fname)
else
call nf_verify(nf90_put_att(ncid, varid, 'coordinates', 'x y'), &
nc_fname)
end if
call nf_verify(nf90_put_att(ncid, varid, 'coordinates', 'x y'), &
nc_fname)
call nf_verify(nf90_put_att(ncid, varid, 'grid_mapping', gridmap_name), &
nc_fname)
else if (latlon) then
call nf_verify(nf90_put_att(ncid, varid, 'coordinates', 'lon lat'), &
nc_fname)
end if
end subroutine ncvar_gridmap

Expand Down
Loading
Loading