Skip to content

Commit 98f87ff

Browse files
mjrenomjreno
authored andcommitted
rework structured export proj and coords, support feet lenuni
1 parent a8d72c5 commit 98f87ff

File tree

6 files changed

+65
-35
lines changed

6 files changed

+65
-35
lines changed

doc/mf6io/mf6ivar/dfn/utl-ncf.dfn

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -10,7 +10,7 @@ shape lenbigline
1010
reader urword
1111
optional true
1212
longname CRS well-known text (WKT) string
13-
description is the coordinate reference system (CRS) well-known text (WKT) string.
13+
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.
1414

1515
block options
1616
name deflate
@@ -96,7 +96,7 @@ shape (ncpl)
9696
optional true
9797
reader readarray
9898
longname cell center latitude
99-
description cell center latitude.
99+
description cell center latitude. Only supported for NETCDF\_STRUCTURED export types.
100100

101101
block griddata
102102
name longitude
@@ -105,4 +105,4 @@ shape (ncpl)
105105
optional true
106106
reader readarray
107107
longname cell center longitude
108-
description cell center longitude.
108+
description cell center longitude. Only supported for NETCDF\_STRUCTURED export types.

src/Utilities/Export/DisNCMesh.f90

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -69,7 +69,7 @@ subroutine dis_export_init(this, modelname, modeltype, modelfname, nc_fname, &
6969

7070
! initialize base class
7171
call this%mesh_init(modelname, modeltype, modelfname, nc_fname, disenum, &
72-
nctype, iout)
72+
nctype, this%dis%lenuni, iout)
7373
end subroutine dis_export_init
7474

7575
!> @brief netcdf export dis destroy

src/Utilities/Export/DisNCStructured.f90

Lines changed: 45 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -75,7 +75,7 @@ module DisNCStructuredModule
7575
procedure :: define_dim
7676
procedure :: define_dependent
7777
procedure :: define_gridmap
78-
procedure :: define_projection
78+
procedure :: define_geocoords
7979
procedure :: add_proj_data
8080
procedure :: add_grid_data
8181
end type DisNCStructuredType
@@ -153,11 +153,24 @@ subroutine dis_export_init(this, modelname, modeltype, modelfname, nc_fname, &
153153

154154
if (latsz > 0 .and. lonsz > 0) then
155155
this%latlon = .true.
156+
if (this%wkt /= '') then
157+
write (warnmsg, '(a)') 'Ignoring user provided NetCDF wkt parameter &
158+
&as longitude and latitude arrays have been provided.'
159+
call store_warning(warnmsg)
160+
this%wkt = ''
161+
this%gridmap_name = ''
162+
end if
156163
call mem_setptr(this%latitude, 'LATITUDE', this%ncf_mempath)
157164
call mem_setptr(this%longitude, 'LONGITUDE', this%ncf_mempath)
158165
end if
159166
end if
160167

168+
if (this%dis%lenuni == 1) then
169+
this%lenunits = 'ft'
170+
else
171+
this%lenunits = 'm'
172+
end if
173+
161174
! create the netcdf file
162175
call nf_verify(nf90_create(this%nc_fname, &
163176
IOR(NF90_CLOBBER, NF90_NETCDF4), this%ncid), &
@@ -191,7 +204,7 @@ subroutine df(this)
191204
! define root group dimensions and coordinate variables
192205
call this%define_dim()
193206
! define grid projection variables
194-
call this%define_projection()
207+
call this%define_geocoords()
195208
if (isim_mode /= MVALIDATE) then
196209
! define the dependent variable
197210
call this%define_dependent()
@@ -707,8 +720,8 @@ subroutine define_dim(this)
707720
this%nc_fname)
708721
call nf_verify(nf90_def_var(this%ncid, 'y', NF90_DOUBLE, this%dim_ids%y, &
709722
this%var_ids%y), this%nc_fname)
710-
call nf_verify(nf90_put_att(this%ncid, this%var_ids%y, 'units', 'm'), &
711-
this%nc_fname)
723+
call nf_verify(nf90_put_att(this%ncid, this%var_ids%y, 'units', &
724+
this%lenunits), this%nc_fname)
712725
call nf_verify(nf90_put_att(this%ncid, this%var_ids%y, 'axis', 'Y'), &
713726
this%nc_fname)
714727
call nf_verify(nf90_put_att(this%ncid, this%var_ids%y, 'standard_name', &
@@ -730,8 +743,8 @@ subroutine define_dim(this)
730743
this%nc_fname)
731744
call nf_verify(nf90_def_var(this%ncid, 'x', NF90_DOUBLE, this%dim_ids%x, &
732745
this%var_ids%x), this%nc_fname)
733-
call nf_verify(nf90_put_att(this%ncid, this%var_ids%x, 'units', 'm'), &
734-
this%nc_fname)
746+
call nf_verify(nf90_put_att(this%ncid, this%var_ids%x, 'units', &
747+
this%lenunits), this%nc_fname)
735748
call nf_verify(nf90_put_att(this%ncid, this%var_ids%x, 'axis', 'X'), &
736749
this%nc_fname)
737750
call nf_verify(nf90_put_att(this%ncid, this%var_ids%x, 'standard_name', &
@@ -782,7 +795,7 @@ subroutine define_dependent(this)
782795

783796
! put attr
784797
call nf_verify(nf90_put_att(this%ncid, this%var_ids%dependent, &
785-
'units', 'm'), this%nc_fname)
798+
'units', this%lenunits), this%nc_fname)
786799
call nf_verify(nf90_put_att(this%ncid, this%var_ids%dependent, &
787800
'standard_name', this%annotation%stdname), &
788801
this%nc_fname)
@@ -816,9 +829,9 @@ end subroutine define_gridmap
816829

817830
!> @brief define grid projection variables
818831
!<
819-
subroutine define_projection(this)
832+
subroutine define_geocoords(this)
820833
class(DisNCStructuredType), intent(inout) :: this
821-
if (this%latlon .and. this%wkt /= '') then
834+
if (this%latlon) then
822835
! lat
823836
call nf_verify(nf90_def_var(this%ncid, 'lat', NF90_DOUBLE, &
824837
(/this%dim_ids%x, this%dim_ids%y/), &
@@ -841,13 +854,13 @@ subroutine define_projection(this)
841854
call nf_verify(nf90_put_att(this%ncid, this%var_ids%longitude, &
842855
'long_name', 'longitude'), this%nc_fname)
843856
end if
844-
end subroutine define_projection
857+
end subroutine define_geocoords
845858

846859
!> @brief add grid projection data
847860
!<
848861
subroutine add_proj_data(this)
849862
class(DisNCStructuredType), intent(inout) :: this
850-
if (this%latlon .and. this%wkt /= '') then
863+
if (this%latlon) then
851864
! lat
852865
call nf_verify(nf90_put_var(this%ncid, this%var_ids%latitude, &
853866
this%latitude, start=(/1, 1/), &
@@ -865,20 +878,30 @@ end subroutine add_proj_data
865878
!> @brief add grid coordinates
866879
!<
867880
subroutine add_grid_data(this)
881+
use ConstantsModule, only: DZERO
868882
class(DisNCStructuredType), intent(inout) :: this
869883
integer(I4B) :: ibnd, n !, k, i, j
870884
real(DP), dimension(:, :), pointer, contiguous :: dbl2d
871885
real(DP), dimension(:), allocatable :: x, y
886+
real(DP) :: xoff, yoff
887+
888+
if (this%dis%angrot /= DZERO) then
889+
xoff = DZERO
890+
yoff = DZERO
891+
else
892+
xoff = this%dis%xorigin
893+
yoff = this%dis%yorigin
894+
end if
872895

873896
allocate (x(size(this%dis%cellx)))
874897
allocate (y(size(this%dis%celly)))
875898

876899
do n = 1, size(this%dis%cellx)
877-
x(n) = this%dis%cellx(n) + this%dis%xorigin
900+
x(n) = this%dis%cellx(n) + xoff
878901
end do
879902

880903
do n = 1, size(this%dis%celly)
881-
y(n) = this%dis%celly(n) + this%dis%yorigin
904+
y(n) = this%dis%celly(n) + yoff
882905
end do
883906

884907
call nf_verify(nf90_put_var(this%ncid, this%var_ids%x, x), &
@@ -897,8 +920,8 @@ subroutine add_grid_data(this)
897920
ibnd = 1
898921
do n = 1, size(this%dis%cellx)
899922
if (ibnd == 1) then
900-
dbl2d(1, ibnd) = this%dis%xorigin
901-
dbl2d(2, ibnd) = this%dis%xorigin + this%dis%delr(ibnd)
923+
dbl2d(1, ibnd) = xoff
924+
dbl2d(2, ibnd) = xoff + this%dis%delr(ibnd)
902925
else
903926
dbl2d(1, ibnd) = dbl2d(1, ibnd - 1) + this%dis%delr(ibnd)
904927
dbl2d(2, ibnd) = dbl2d(2, ibnd - 1) + this%dis%delr(ibnd)
@@ -914,8 +937,8 @@ subroutine add_grid_data(this)
914937
ibnd = 1
915938
do n = size(this%dis%celly), 1, -1
916939
if (ibnd == 1) then
917-
dbl2d(1, ibnd) = this%dis%yorigin + sum(this%dis%delc) - this%dis%delc(n)
918-
dbl2d(2, ibnd) = this%dis%yorigin + sum(this%dis%delc)
940+
dbl2d(1, ibnd) = yoff + sum(this%dis%delc) - this%dis%delc(n)
941+
dbl2d(2, ibnd) = yoff + sum(this%dis%delc)
919942
else
920943
dbl2d(1, ibnd) = dbl2d(1, ibnd - 1) - this%dis%delc(n)
921944
dbl2d(2, ibnd) = dbl2d(2, ibnd - 1) - this%dis%delc(n)
@@ -982,15 +1005,13 @@ subroutine ncvar_gridmap(ncid, varid, gridmap_name, latlon, nc_fname)
9821005
logical(LGP), intent(in) :: latlon
9831006
character(len=*), intent(in) :: nc_fname
9841007
if (gridmap_name /= '') then
985-
if (latlon) then
986-
call nf_verify(nf90_put_att(ncid, varid, 'coordinates', 'lon lat'), &
987-
nc_fname)
988-
else
989-
call nf_verify(nf90_put_att(ncid, varid, 'coordinates', 'x y'), &
990-
nc_fname)
991-
end if
1008+
call nf_verify(nf90_put_att(ncid, varid, 'coordinates', 'x y'), &
1009+
nc_fname)
9921010
call nf_verify(nf90_put_att(ncid, varid, 'grid_mapping', gridmap_name), &
9931011
nc_fname)
1012+
else if (latlon) then
1013+
call nf_verify(nf90_put_att(ncid, varid, 'coordinates', 'lon lat'), &
1014+
nc_fname)
9941015
end if
9951016
end subroutine ncvar_gridmap
9961017

src/Utilities/Export/DisvNCMesh.f90

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -67,7 +67,7 @@ subroutine disv_export_init(this, modelname, modeltype, modelfname, nc_fname, &
6767

6868
! initialize base class
6969
call this%mesh_init(modelname, modeltype, modelfname, nc_fname, disenum, &
70-
nctype, iout)
70+
nctype, this%disv%lenuni, iout)
7171
end subroutine disv_export_init
7272

7373
!> @brief netcdf export disv destroy

src/Utilities/Export/MeshNCModel.f90

Lines changed: 13 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -96,7 +96,7 @@ subroutine nc_array_export_if(this, pkgtype, pkgname, mempath, idt)
9696
!> @brief initialize
9797
!<
9898
subroutine mesh_init(this, modelname, modeltype, modelfname, nc_fname, &
99-
disenum, nctype, iout)
99+
disenum, nctype, lenuni, iout)
100100
use MemoryManagerExtModule, only: mem_set_value
101101
class(MeshModelType), intent(inout) :: this
102102
character(len=*), intent(in) :: modelname
@@ -105,6 +105,7 @@ subroutine mesh_init(this, modelname, modeltype, modelfname, nc_fname, &
105105
character(len=*), intent(in) :: nc_fname
106106
integer(I4B), intent(in) :: disenum
107107
integer(I4B), intent(in) :: nctype
108+
integer(I4B), intent(in) :: lenuni
108109
integer(I4B), intent(in) :: iout
109110
logical(LGP) :: found
110111

@@ -131,6 +132,12 @@ subroutine mesh_init(this, modelname, modeltype, modelfname, nc_fname, &
131132
call store_warning(warnmsg)
132133
end if
133134

135+
if (lenuni == 1) then
136+
this%lenunits = 'ft'
137+
else
138+
this%lenunits = 'm'
139+
end if
140+
134141
! create the netcdf file
135142
call nf_verify(nf90_create(this%nc_fname, &
136143
IOR(NF90_CLOBBER, NF90_NETCDF4), this%ncid), &
@@ -298,7 +305,7 @@ subroutine define_dependent(this)
298305

299306
! assign variable attributes
300307
call nf_verify(nf90_put_att(this%ncid, this%var_ids%dependent(k), &
301-
'units', 'm'), this%nc_fname)
308+
'units', this%lenunits), this%nc_fname)
302309
call nf_verify(nf90_put_att(this%ncid, this%var_ids%dependent(k), &
303310
'standard_name', this%annotation%stdname), &
304311
this%nc_fname)
@@ -377,7 +384,7 @@ subroutine create_mesh(this)
377384

378385
! assign mesh x node variable attributes
379386
call nf_verify(nf90_put_att(this%ncid, this%var_ids%mesh_node_x, &
380-
'units', 'm'), this%nc_fname)
387+
'units', this%lenunits), this%nc_fname)
381388
call nf_verify(nf90_put_att(this%ncid, this%var_ids%mesh_node_x, &
382389
'standard_name', 'projection_x_coordinate'), &
383390
this%nc_fname)
@@ -398,7 +405,7 @@ subroutine create_mesh(this)
398405

399406
! assign mesh y variable attributes
400407
call nf_verify(nf90_put_att(this%ncid, this%var_ids%mesh_node_y, &
401-
'units', 'm'), this%nc_fname)
408+
'units', this%lenunits), this%nc_fname)
402409
call nf_verify(nf90_put_att(this%ncid, this%var_ids%mesh_node_y, &
403410
'standard_name', 'projection_y_coordinate'), &
404411
this%nc_fname)
@@ -419,7 +426,7 @@ subroutine create_mesh(this)
419426

420427
! assign mesh x face variable attributes
421428
call nf_verify(nf90_put_att(this%ncid, this%var_ids%mesh_face_x, &
422-
'units', 'm'), this%nc_fname)
429+
'units', this%lenunits), this%nc_fname)
423430
call nf_verify(nf90_put_att(this%ncid, this%var_ids%mesh_face_x, &
424431
'standard_name', 'projection_x_coordinate'), &
425432
this%nc_fname)
@@ -448,7 +455,7 @@ subroutine create_mesh(this)
448455

449456
! assign mesh y face variable attributes
450457
call nf_verify(nf90_put_att(this%ncid, this%var_ids%mesh_face_y, &
451-
'units', 'm'), this%nc_fname)
458+
'units', this%lenunits), this%nc_fname)
452459
call nf_verify(nf90_put_att(this%ncid, this%var_ids%mesh_face_y, &
453460
'standard_name', 'projection_y_coordinate'), &
454461
this%nc_fname)

src/Utilities/Export/NCModel.f90

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -77,6 +77,7 @@ module NCModelExportModule
7777
character(len=LENBIGLINE) :: wkt !< wkt user string
7878
character(len=LINELENGTH) :: datetime !< export file creation time
7979
character(len=LINELENGTH) :: xname !< dependent variable name
80+
character(len=LINELENGTH) :: lenunits !< unidata udunits length units
8081
type(NCExportAnnotation) :: annotation !< export file annotation
8182
real(DP), dimension(:), pointer, contiguous :: x !< dependent variable pointer
8283
integer(I4B) :: disenum !< type of discretization
@@ -296,6 +297,7 @@ subroutine export_init(this, modelname, modeltype, modelfname, nc_fname, &
296297
this%wkt = ''
297298
this%datetime = ''
298299
this%xname = ''
300+
this%lenunits = ''
299301
this%disenum = disenum
300302
this%ncid = 0
301303
this%stepcnt = 0

0 commit comments

Comments
 (0)