ncmodelexportmodule Module Reference

This module contains the NCModelExportModule. More...

Data Types
type	exportpackagetype
type	ncexportannotation
	netcdf export attribute annotations More...
type	ncmodelexporttype
	base class for an export model More...
type	ncbasemodelexporttype
	abstract type for model netcdf export type More...
interface	model_define
	abstract interfaces for model netcdf export type More...
interface	model_step
interface	package_export
interface	package_export_ilayer

Enumerations
enum
	netcdf export types enumerator More...

Functions/Subroutines
subroutine	epkg_init (this, mf6_input, mshape, naux, param_names, nparam)
	initialize dynamic package export object More...
subroutine	epkg_destroy (this)
	destroy dynamic package export object More...
subroutine	set (this, modelname, modeltype, modelfname, nctype, disenum)
	set netcdf file scoped attributes More...
subroutine	export_init (this, modelname, modeltype, modelfname, nc_fname, disenum, nctype, iout)
	initialization of model netcdf export More...
class(exportpackagetype) function, pointer	export_get (this, idx)
	retrieve dynamic export object from package list More...
character(len=linelength) function	input_attribute (this, pkgname, idt)
	build modflow_input attribute string More...
integer(i4b) function	istp (this)
	step index for timeseries data More...
character(len=linelength) function, public	export_varname (pkgname, tagname, mempath, layer, iaux)
	build netcdf variable name More...
character(len=linelength) function, public	export_longname (longname, pkgname, tagname, mempath, layer, iaux, component_type, subcomponent_type)
	build netcdf variable longname More...
subroutine	export_input (this)
	netcdf dynamic package period export More...
character(len=linelength) function, public	wkt_to_cf_gridmapping (wkt)
	map a WKT string to a CF-1.13 grid_mapping_name More...
real(dp) function	wkt_find_param (wkt, name, found)
	find the first number following a quoted parameter name More...
subroutine	wkt_find_ellipsoid (wkt, semi_major_axis, inverse_flattening, found)
	find semi-major axis and inverse flattening from a WKT string More...
subroutine, public	wkt_transverse_mercator_params (wkt, longitude_of_central_meridian, latitude_of_projection_origin, scale_factor_at_central_meridian, false_easting, false_northing, semi_major_axis, inverse_flattening, found)
	extract CF transverse_mercator grid_mapping parameters More...
subroutine, public	wkt_lambert_conformal_conic_params (wkt, standard_parallel_1, standard_parallel_2, longitude_of_central_meridian, latitude_of_projection_origin, false_easting, false_northing, semi_major_axis, inverse_flattening, is_2sp, found)
	extract CF lambert_conformal_conic grid_mapping parameters More...
subroutine, public	wkt_albers_conical_equal_area_params (wkt, standard_parallel_1, standard_parallel_2, longitude_of_central_meridian, latitude_of_projection_origin, false_easting, false_northing, semi_major_axis, inverse_flattening, found)
	extract CF albers_conical_equal_area grid_mapping parameters More...
character(len=lenbigline) function, public	wrap_rotated_crs (wkt, xorigin, yorigin, angrot)
	wrap a WKT2 PROJCRS in a derived CRS encoding grid rotation More...
subroutine	export_destroy (this)
	destroy model netcdf export object More...

Variables
@, public	netcdf_undef = 0
	undefined netcdf export type More...
@, public	netcdf_structured = 1
	netcdf structrured export More...
@, public	netcdf_mesh2d = 2
	netcdf ugrid layered mesh export More...

Detailed Description

This module defines a model export and base type for supported netcdf files and is not dependent on netcdf libraries.

Enumeration Type Documentation

anonymous enum

anonymous enum

private

Definition at line 35 of file NCModel.f90.

Function/Subroutine Documentation

epkg_destroy()

subroutine ncmodelexportmodule::epkg_destroy

(

class(exportpackagetype), intent(inout)

this

)

Definition at line 198 of file NCModel.f90.

    use inputdefinitionmodule, only: inputparamdefinitiontype
    class(ExportPackageType), intent(inout) :: this
    if (allocated(this%param_names)) deallocate (this%param_names)

epkg_init()

subroutine ncmodelexportmodule::epkg_init ( class(exportpackagetype), intent(inout)  this, type(modflowinputtype), intent(in)  mf6_input, integer(i4b), dimension(:), intent(in), pointer, contiguous  mshape, integer(i4b), intent(in)  naux, character(len=linelength), dimension(:), intent(in), allocatable  param_names, integer(i4b), intent(in)  nparam  )

private

Parameters

[in]

mshape

model shape

Definition at line 150 of file NCModel.f90.

    use simvariablesmodule, only: idm_context
    use memorymanagermodule, only: mem_setptr
    use memorymanagerextmodule, only: mem_set_value
    use memoryhelpermodule, only: create_mem_path
    class(ExportPackageType), intent(inout) :: this
    type(ModflowInputType), intent(in) :: mf6_input
    integer(I4B), dimension(:), pointer, contiguous, intent(in) :: mshape !< model shape
    integer(I4B), intent(in) :: naux
    character(len=LINELENGTH), dimension(:), allocatable, &
      intent(in) :: param_names
    integer(I4B), intent(in) :: nparam
    integer(I4B) :: n
    character(len=LENVARNAME) :: rs_varname
    character(len=LENMEMPATH) :: input_mempath
    integer(I4B), pointer :: rsvar
 
    this%mf6_input = mf6_input
    this%mshape => mshape
    this%nparam = nparam
    this%naux = naux
    this%eper = 0
 
    input_mempath = create_mem_path(component=mf6_input%component_name, &
                                    subcomponent=mf6_input%subcomponent_name, &
                                    context=idm_context)
 
    ! allocate param arrays
    allocate (this%param_names(nparam))
    allocate (this%param_reads(nparam))
    allocate (this%varids_param(nparam, mshape(1)))
    allocate (this%varids_aux(naux, mshape(1)))
 
    ! set param arrays
    do n = 1, nparam
      this%param_names(n) = param_names(n)
      rs_varname = rsv_name(param_names(n))
      call mem_setptr(rsvar, rs_varname, mf6_input%mempath)
      this%param_reads(n)%invar => rsvar
    end do
 
    ! set pointer to loaded input period
    call mem_setptr(this%iper, 'IPER', mf6_input%mempath)

Here is the call graph for this function:

export_destroy()

subroutine ncmodelexportmodule::export_destroy

(

class(ncmodelexporttype), intent(inout)

this

)

Definition at line 1012 of file NCModel.f90.

    use memorymanagerextmodule, only: memorystore_remove
    use simvariablesmodule, only: idm_context
    class(NCModelExportType), intent(inout) :: this
    ! override in derived class
    deallocate (this%deflate)
    deallocate (this%shuffle)
    deallocate (this%input_attr)
    deallocate (this%chunk_time)
    ! Deallocate idm memory
    if (this%ncf_mempath /= '') then
      call memorystore_remove(this%modelname, 'NCF', idm_context)
    end if

Here is the call graph for this function:

export_get()

class(exportpackagetype) function, pointer ncmodelexportmodule::export_get	(	class(ncmodelexporttype), intent(inout)	this,
		integer(i4b), intent(in)	idx
	)

Definition at line 403 of file NCModel.f90.

    use listmodule, only: listtype
    class(NCModelExportType), intent(inout) :: this
    integer(I4B), intent(in) :: idx
    class(ExportPackageType), pointer :: res
    class(*), pointer :: obj
    nullify (res)
    obj => this%pkglist%GetItem(idx)
    if (associated(obj)) then
      select type (obj)
      class is (exportpackagetype)
        res => obj
      end select
    end if

export_init()

subroutine ncmodelexportmodule::export_init	(	class(ncmodelexporttype), intent(inout)	this,
		character(len=*), intent(in)	modelname,
		character(len=*), intent(in)	modeltype,
		character(len=*), intent(in)	modelfname,
		character(len=*), intent(in)	nc_fname,
		integer(i4b), intent(in)	disenum,
		integer(i4b), intent(in)	nctype,
		integer(i4b), intent(in)	iout
	)

Definition at line 281 of file NCModel.f90.

    use tdismodule, only: datetime0, nper, nstp
    use memorymanagermodule, only: mem_setptr
    use memoryhelpermodule, only: create_mem_path
    use memorymanagerextmodule, only: mem_set_value
    use inputoutputmodule, only: lowcase
    use sourcecommonmodule, only: filein_fname
    use utlncfinputmodule, only: utlncfparamfoundtype
    class(NCModelExportType), intent(inout) :: this
    character(len=*), intent(in) :: modelname
    character(len=*), intent(in) :: modeltype
    character(len=*), intent(in) :: modelfname
    character(len=*), intent(in) :: nc_fname
    integer(I4B), intent(in) :: disenum
    integer(I4B), intent(in) :: nctype
    integer(I4B), intent(in) :: iout
    character(len=LENMEMPATH) :: model_mempath
    type(UtlNcfParamFoundType) :: ncf_found
 
    ! allocate
    allocate (this%deflate)
    allocate (this%shuffle)
    allocate (this%input_attr)
    allocate (this%chunk_time)
 
    ! initialize
    this%modelname = modelname
    this%modeltype = modeltype
    this%modelfname = modelfname
    this%nc_fname = nc_fname
    this%gridmap_name = ''
    this%ncf_mempath = ''
    this%wkt = ''
    this%crs_wkt = ''
    this%datetime = ''
    this%xname = ''
    this%lenunits = ''
    this%disenum = disenum
    this%ncid = 0
    this%totnstp = 0
    this%deflate = -1
    this%shuffle = 0
    this%input_attr = 1
    this%chunk_time = -1
    this%iout = iout
    this%chunking_active = .false.
 
    ! set file scoped attributes
    call this%annotation%set(modelname, modeltype, modelfname, nctype, disenum)
 
    ! set dependent variable basename
    select case (modeltype)
    case ('GWF')
      this%xname = 'head'
    case ('GWT')
      this%xname = 'concentration'
    case ('GWE')
      this%xname = 'temperature'
    case default
      errmsg = trim(modeltype)//' models not supported for NetCDF export.'
      call store_error(errmsg)
      call store_error_filename(modelfname)
    end select
 
    ! set discretization input mempath
    if (disenum == dis) then
      this%dis_mempath = create_mem_path(modelname, 'DIS', idm_context)
    else if (disenum == disu) then
      this%dis_mempath = create_mem_path(modelname, 'DISU', idm_context)
    else if (disenum == disv) then
      this%dis_mempath = create_mem_path(modelname, 'DISV', idm_context)
    end if
 
    ! set dependent variable pointer
    model_mempath = create_mem_path(component=modelname)
    call mem_setptr(this%x, 'X', model_mempath)
 
    ! set ncf_mempath if provided
    if (filein_fname(this%ncf_mempath, 'NCF6_MEMPATH', this%dis_mempath, &
                     modelfname)) then
      call mem_set_value(this%wkt, 'WKT', this%ncf_mempath, &
                         ncf_found%wkt)
      call mem_set_value(this%crs_wkt, 'CRS_WKT', this%ncf_mempath, &
                         ncf_found%crs_wkt)
      call mem_set_value(this%deflate, 'DEFLATE', this%ncf_mempath, &
                         ncf_found%deflate)
      call mem_set_value(this%shuffle, 'SHUFFLE', this%ncf_mempath, &
                         ncf_found%shuffle)
      call mem_set_value(this%input_attr, 'ATTR_OFF', this%ncf_mempath, &
                         ncf_found%attr_off)
      call mem_set_value(this%chunk_time, 'CHUNK_TIME', this%ncf_mempath, &
                         ncf_found%chunk_time)
    end if
 
    if (ncf_found%wkt .or. ncf_found%crs_wkt) then
      this%gridmap_name = 'projection'
    end if
 
    ! ATTR_OFF turns off modflow 6 input attributes
    if (ncf_found%attr_off) then
      this%input_attr = 0
    end if
 
    ! set datetime string
    if (datetime0 /= '') then
      this%datetime = 'days since '//trim(datetime0)
    else
      ! January 1, 1970 at 00:00:00 UTC
      this%datetime = 'days since 1970-01-01T00:00:00'
    end if
 
    ! set total nstp
    if (isim_mode == mvalidate) then
      this%totnstp = nper
    else
      this%totnstp = sum(nstp)
    end if

Here is the call graph for this function:

export_input()

subroutine ncmodelexportmodule::export_input

(

class(ncbasemodelexporttype), intent(inout)

this

)

Definition at line 548 of file NCModel.f90.

    use tdismodule, only: kper
    class(NCBaseModelExportType), intent(inout) :: this
    integer(I4B) :: idx
    class(ExportPackageType), pointer :: export_pkg
    do idx = 1, this%pkglist%Count()
      export_pkg => this%get(idx)
      ! period input already exported
      if (export_pkg%eper >= kper) cycle
      ! update export package
      call this%package_step(export_pkg)
      ! update exported iper
      export_pkg%eper = kper
    end do

export_longname()

character(len=linelength) function, public ncmodelexportmodule::export_longname	(	character(len=*), intent(in)	longname,
		character(len=*), intent(in)	pkgname,
		character(len=*), intent(in)	tagname,
		character(len=*), intent(in)	mempath,
		integer(i4b), intent(in), optional	layer,
		integer(i4b), intent(in), optional	iaux,
		character(len=*), intent(in), optional	component_type,
		character(len=*), intent(in), optional	subcomponent_type
	)

Definition at line 494 of file NCModel.f90.

    use memorymanagermodule, only: mem_setptr
    use characterstringmodule, only: characterstringtype
    use inputoutputmodule, only: lowcase
    use idmdfnselectormodule, only: idm_multi_package
    character(len=*), intent(in) :: longname
    character(len=*), intent(in) :: pkgname
    character(len=*), intent(in) :: tagname
    character(len=*), intent(in) :: mempath
    integer(I4B), optional, intent(in) :: layer
    integer(I4B), optional, intent(in) :: iaux
    character(len=*), optional, intent(in) :: component_type
    character(len=*), optional, intent(in) :: subcomponent_type
    character(len=LINELENGTH) :: lname
    type(CharacterStringType), dimension(:), pointer, &
      contiguous :: auxnames
    character(len=LINELENGTH) :: pname, vname, auxname
    pname = pkgname
    vname = tagname
    call lowcase(pname)
    call lowcase(vname)
    if (longname == '') then
      lname = trim(pname)//' '//trim(vname)
    else
      lname = longname
      if (present(component_type) .and. present(subcomponent_type)) then
        if (idm_multi_package(component_type, subcomponent_type)) then
          lname = trim(pname)//' '//trim(lname)
        end if
      end if
    end if
 
    if (present(iaux)) then
      if (iaux > 0) then
        if (tagname == 'AUX') then
          ! reset vname to auxiliary variable name
          call mem_setptr(auxnames, 'AUXILIARY', mempath)
          auxname = auxnames(iaux)
          call lowcase(auxname)
          lname = trim(lname)//' '//trim(auxname)
        end if
      end if
    end if
 
    if (present(layer)) then
      if (layer > 0) then
        write (lname, '(a,i0)') trim(lname)//' layer ', layer
      end if
    end if

Here is the call graph for this function:

Here is the caller graph for this function:

export_varname()

character(len=linelength) function, public ncmodelexportmodule::export_varname	(	character(len=*), intent(in)	pkgname,
		character(len=*), intent(in)	tagname,
		character(len=*), intent(in)	mempath,
		integer(i4b), intent(in), optional	layer,
		integer(i4b), intent(in), optional	iaux
	)

Definition at line 453 of file NCModel.f90.

    use memorymanagermodule, only: mem_setptr
    use characterstringmodule, only: characterstringtype
    use inputoutputmodule, only: lowcase
    character(len=*), intent(in) :: pkgname
    character(len=*), intent(in) :: tagname
    character(len=*), intent(in) :: mempath
    integer(I4B), optional, intent(in) :: layer
    integer(I4B), optional, intent(in) :: iaux
    character(len=LINELENGTH) :: varname
    type(CharacterStringType), dimension(:), pointer, &
      contiguous :: auxnames
    character(len=LINELENGTH) :: pname, vname
    vname = tagname
    pname = pkgname
 
    if (present(iaux)) then
      if (iaux > 0) then
        if (tagname == 'AUX') then
          ! reset vname to auxiliary variable name
          call mem_setptr(auxnames, 'AUXILIARY', mempath)
          vname = auxnames(iaux)
        end if
      end if
    end if
 
    call lowcase(vname)
    call lowcase(pname)
    varname = trim(pname)//'_'//trim(vname)
 
    if (present(layer)) then
      if (layer > 0) then
        !write (varname, '(a,i0)') trim(varname)//'_L', layer
        write (varname, '(a,i0)') trim(varname)//'_l', layer
      end if
    end if

Here is the call graph for this function:

Here is the caller graph for this function:

input_attribute()

character(len=linelength) function ncmodelexportmodule::input_attribute	(	class(ncmodelexporttype), intent(inout)	this,
		character(len=*), intent(in)	pkgname,
		type(inputparamdefinitiontype), intent(in), pointer	idt
	)

Definition at line 421 of file NCModel.f90.

    use inputoutputmodule, only: lowcase
    use memoryhelpermodule, only: mempathseparator
    use inputdefinitionmodule, only: inputparamdefinitiontype
    class(NCModelExportType), intent(inout) :: this
    character(len=*), intent(in) :: pkgname
    type(InputParamDefinitionType), pointer, intent(in) :: idt
    character(len=LINELENGTH) :: attr
    attr = ''
    if (this%input_attr > 0) then
      attr = trim(this%modelname)//mempathseparator//trim(pkgname)// &
             mempathseparator//trim(idt%tagname)
      call lowcase(attr)
    end if

Here is the call graph for this function:

istp()

integer(i4b) function ncmodelexportmodule::istp

(

class(ncmodelexporttype), intent(inout)

this

)

Definition at line 439 of file NCModel.f90.

    use tdismodule, only: kstp, kper, nstp
    class(NCModelExportType), intent(inout) :: this
    integer(I4B) :: n, istp
    istp = kstp
    if (kper > 1) then
      do n = 1, kper - 1
        istp = istp + nstp(n)
      end do
    end if

set()

subroutine ncmodelexportmodule::set	(	class(ncexportannotation), intent(inout)	this,
		character(len=*), intent(in)	modelname,
		character(len=*), intent(in)	modeltype,
		character(len=*), intent(in)	modelfname,
		integer(i4b), intent(in)	nctype,
		integer(i4b), intent(in)	disenum
	)

Definition at line 206 of file NCModel.f90.

    use versionmodule, only: fullversion
    use inputoutputmodule, only: lowcase
    class(NCExportAnnotation), intent(inout) :: this
    character(len=*), intent(in) :: modelname
    character(len=*), intent(in) :: modeltype
    character(len=*), intent(in) :: modelfname
    integer(I4B), intent(in) :: nctype
    integer(I4B), intent(in) :: disenum
    integer :: values(8)
 
    this%title = ''
    this%model = ''
    this%mesh = ''
    this%grid = ''
    this%history = ''
    this%source = ''
    this%conventions = ''
    this%stdname = ''
    this%longname = ''
 
    ! set file conventions
    this%conventions = 'CF-1.13'
    if (nctype == netcdf_mesh2d) this%conventions = &
      trim(this%conventions)//' UGRID-1.0'
 
    ! set model specific attributes
    select case (modeltype)
    case ('GWF')
      this%title = trim(modelname)//' hydraulic head'
      this%longname = 'head'
    case ('GWT')
      this%title = trim(modelname)//' concentration'
      this%longname = 'concentration'
    case ('GWE')
      this%title = trim(modelname)//' temperature'
      this%longname = 'temperature'
    case default
      errmsg = trim(modeltype)//' models not supported for NetCDF export.'
      call store_error(errmsg)
      call store_error_filename(modelfname)
    end select
 
    if (isim_mode == mvalidate) then
      this%title = trim(this%title)//' array input'
    end if
 
    ! set mesh type
    if (nctype == netcdf_mesh2d) then
      this%mesh = 'layered'
    end if
 
    ! set grid type
    if (disenum == dis) then
      this%grid = 'structured'
    else if (disenum == disv) then
      this%grid = 'vertex'
    end if
 
    ! model description string (lowercase, no version suffix, matching CF convention)
    this%model = trim(modeltype)//': '//trim(modelname)
    call lowcase(this%model)
 
    ! modflow6 version string
    this%source = 'MODFLOW 6 '//trim(adjustl(fullversion))
 
    ! create timestamp
    call date_and_time(values=values)
    write (this%history, '(a,i0,a,i0,a,i0,a,i0,a,i0,a,i0,a,i0)') &
      'first created ', values(1), '/', values(2), '/', values(3), ' ', &
      values(5), ':', values(6), ':', values(7), '.', values(8)

Here is the call graph for this function:

wkt_albers_conical_equal_area_params()

subroutine, public ncmodelexportmodule::wkt_albers_conical_equal_area_params	(	character(len=*), intent(in)	wkt,
		real(dp), intent(out)	standard_parallel_1,
		real(dp), intent(out)	standard_parallel_2,
		real(dp), intent(out)	longitude_of_central_meridian,
		real(dp), intent(out)	latitude_of_projection_origin,
		real(dp), intent(out)	false_easting,
		real(dp), intent(out)	false_northing,
		real(dp), intent(out)	semi_major_axis,
		real(dp), intent(out)	inverse_flattening,
		logical(lgp), intent(out)	found
	)

Best-effort extraction from a WKT1 or WKT2 string – MF6 has no CRS parsing library. found=.false. if any parameter cannot be located.

Definition at line 867 of file NCModel.f90.

    character(len=*), intent(in) :: wkt
    real(DP), intent(out) :: standard_parallel_1
    real(DP), intent(out) :: standard_parallel_2
    real(DP), intent(out) :: longitude_of_central_meridian
    real(DP), intent(out) :: latitude_of_projection_origin
    real(DP), intent(out) :: false_easting
    real(DP), intent(out) :: false_northing
    real(DP), intent(out) :: semi_major_axis
    real(DP), intent(out) :: inverse_flattening
    logical(LGP), intent(out) :: found
    logical(LGP) :: f1, f2, f3, f4, f5, f6, f7
 
    standard_parallel_1 = wkt_find_param(wkt, 'standard_parallel_1', f1)
    if (.not. f1) standard_parallel_1 = &
      wkt_find_param(wkt, 'Latitude of 1st standard parallel', f1)
 
    standard_parallel_2 = wkt_find_param(wkt, 'standard_parallel_2', f2)
    if (.not. f2) standard_parallel_2 = &
      wkt_find_param(wkt, 'Latitude of 2nd standard parallel', f2)
 
    longitude_of_central_meridian = wkt_find_param(wkt, 'longitude_of_center', f3)
    if (.not. f3) longitude_of_central_meridian = &
      wkt_find_param(wkt, 'Longitude of false origin', f3)
 
    latitude_of_projection_origin = wkt_find_param(wkt, 'latitude_of_center', f4)
    if (.not. f4) latitude_of_projection_origin = &
      wkt_find_param(wkt, 'Latitude of false origin', f4)
 
    false_easting = wkt_find_param(wkt, 'false_easting', f5)
    if (.not. f5) false_easting = &
      wkt_find_param(wkt, 'Easting at false origin', f5)
 
    false_northing = wkt_find_param(wkt, 'false_northing', f6)
    if (.not. f6) false_northing = &
      wkt_find_param(wkt, 'Northing at false origin', f6)
 
    call wkt_find_ellipsoid(wkt, semi_major_axis, inverse_flattening, f7)
 
    found = f1 .and. f2 .and. f3 .and. f4 .and. f5 .and. f6 .and. f7

Here is the call graph for this function:

Here is the caller graph for this function:

wkt_find_ellipsoid()

subroutine ncmodelexportmodule::wkt_find_ellipsoid	(	character(len=*), intent(in)	wkt,
		real(dp), intent(out)	semi_major_axis,
		real(dp), intent(out)	inverse_flattening,
		logical(lgp), intent(out)	found
	)

Matches SPHEROID["name",A,B,...] (WKT1) or ELLIPSOID["name",A,B,...] (WKT2). Best-effort only. found=.false. if neither keyword is present or the two numbers cannot be read.

Definition at line 698 of file NCModel.f90.

    use inputoutputmodule, only: upcase
    character(len=*), intent(in) :: wkt
    real(DP), intent(out) :: semi_major_axis
    real(DP), intent(out) :: inverse_flattening
    logical(LGP), intent(out) :: found
    character(len=LENBIGLINE) :: wkt_upper
    integer :: ipos, i, iend, ios, wkt_len, nquote
 
    found = .false.
    semi_major_axis = dzero
    inverse_flattening = dzero
    wkt_len = len_trim(wkt)
    wkt_upper = wkt
    call upcase(wkt_upper)
 
    ipos = index(wkt_upper(1:wkt_len), 'SPHEROID[')
    if (ipos == 0) ipos = index(wkt_upper(1:wkt_len), 'ELLIPSOID[')
    if (ipos == 0) return
 
    ! skip past the quoted ellipsoid name (its closing quote)
    i = ipos
    nquote = 0
    do while (i <= wkt_len)
      if (wkt(i:i) == '"') nquote = nquote + 1
      if (nquote == 2) exit
      i = i + 1
    end do
    if (nquote < 2) return
    i = i + 1
    if (wkt(i:i) == ',') i = i + 1
 
    iend = i
    do while (iend <= wkt_len)
      if (wkt(iend:iend) == ',') exit
      iend = iend + 1
    end do
    if (iend <= i) return
    read (wkt(i:iend - 1), *, iostat=ios) semi_major_axis
    if (ios /= 0) return
 
    i = iend + 1
    iend = i
    do while (iend <= wkt_len)
      if (wkt(iend:iend) == ',' .or. wkt(iend:iend) == ']') exit
      iend = iend + 1
    end do
    if (iend <= i) return
    read (wkt(i:iend - 1), *, iostat=ios) inverse_flattening
    if (ios == 0) found = .true.

Here is the call graph for this function:

Here is the caller graph for this function:

wkt_find_param()

real(dp) function ncmodelexportmodule::wkt_find_param	(	character(len=*), intent(in)	wkt,
		character(len=*), intent(in)	name,
		logical(lgp), intent(out)	found
	)

Matches e.g. PARAMETER["central_meridian",-117] (WKT1) or PARAMETER["Longitude of natural origin",-117,...] (WKT2). Best-effort only – MF6 has no CRS parsing library. found=.false. if name is not present or no number follows it.

Definition at line 656 of file NCModel.f90.

    use inputoutputmodule, only: upcase
    character(len=*), intent(in) :: wkt
    character(len=*), intent(in) :: name
    logical(LGP), intent(out) :: found
    real(DP) :: value
    character(len=LENBIGLINE) :: wkt_upper
    character(len=LINELENGTH) :: target_upper
    integer :: ipos, istart, iend, ios, wkt_len
 
    found = .false.
    value = dzero
    wkt_len = len_trim(wkt)
    target_upper = '"'//trim(name)//'"'
    wkt_upper = wkt
    call upcase(wkt_upper)
    call upcase(target_upper)
    ipos = index(wkt_upper(1:wkt_len), trim(target_upper))
    if (ipos == 0) return
 
    istart = ipos + len_trim(target_upper)
    do while (istart <= wkt_len)
      if (wkt(istart:istart) == ',') exit
      istart = istart + 1
    end do
    istart = istart + 1
    iend = istart
    do while (iend <= wkt_len)
      if (wkt(iend:iend) == ',' .or. wkt(iend:iend) == ']') exit
      iend = iend + 1
    end do
    if (iend <= istart) return
    read (wkt(istart:iend - 1), *, iostat=ios) value
    if (ios == 0) found = .true.

Here is the call graph for this function:

Here is the caller graph for this function:

wkt_lambert_conformal_conic_params()

subroutine, public ncmodelexportmodule::wkt_lambert_conformal_conic_params	(	character(len=*), intent(in)	wkt,
		real(dp), intent(out)	standard_parallel_1,
		real(dp), intent(out)	standard_parallel_2,
		real(dp), intent(out)	longitude_of_central_meridian,
		real(dp), intent(out)	latitude_of_projection_origin,
		real(dp), intent(out)	false_easting,
		real(dp), intent(out)	false_northing,
		real(dp), intent(out)	semi_major_axis,
		real(dp), intent(out)	inverse_flattening,
		logical(lgp), intent(out)	is_2sp,
		logical(lgp), intent(out)	found
	)

Best-effort extraction from a WKT1 or WKT2 string – MF6 has no CRS parsing library. 2SP (two standard parallels) only: CF's lambert_conformal_conic has no scale_factor attribute, so a 1SP (scale-factor-based) WKT cannot be exactly represented without trigonometric conversion MF6 has no library for. is_2sp=.false. if a second standard parallel is not present (i.e. likely 1SP); found=.false. if is_2sp but any other parameter cannot be located.

Definition at line 807 of file NCModel.f90.

    character(len=*), intent(in) :: wkt
    real(DP), intent(out) :: standard_parallel_1
    real(DP), intent(out) :: standard_parallel_2
    real(DP), intent(out) :: longitude_of_central_meridian
    real(DP), intent(out) :: latitude_of_projection_origin
    real(DP), intent(out) :: false_easting
    real(DP), intent(out) :: false_northing
    real(DP), intent(out) :: semi_major_axis
    real(DP), intent(out) :: inverse_flattening
    logical(LGP), intent(out) :: is_2sp
    logical(LGP), intent(out) :: found
    logical(LGP) :: f1, f2, f3, f4, f5, f6, f7
 
    standard_parallel_2 = wkt_find_param(wkt, 'standard_parallel_2', f2)
    if (.not. f2) standard_parallel_2 = &
      wkt_find_param(wkt, 'Latitude of 2nd standard parallel', f2)
    is_2sp = f2
    if (.not. is_2sp) then
      found = .false.
      return
    end if
 
    standard_parallel_1 = wkt_find_param(wkt, 'standard_parallel_1', f1)
    if (.not. f1) standard_parallel_1 = &
      wkt_find_param(wkt, 'Latitude of 1st standard parallel', f1)
 
    longitude_of_central_meridian = wkt_find_param(wkt, 'central_meridian', f3)
    if (.not. f3) longitude_of_central_meridian = &
      wkt_find_param(wkt, 'Longitude of false origin', f3)
 
    latitude_of_projection_origin = wkt_find_param(wkt, 'latitude_of_origin', f4)
    if (.not. f4) latitude_of_projection_origin = &
      wkt_find_param(wkt, 'Latitude of false origin', f4)
 
    false_easting = wkt_find_param(wkt, 'false_easting', f5)
    if (.not. f5) false_easting = &
      wkt_find_param(wkt, 'Easting at false origin', f5)
 
    false_northing = wkt_find_param(wkt, 'false_northing', f6)
    if (.not. f6) false_northing = &
      wkt_find_param(wkt, 'Northing at false origin', f6)
 
    call wkt_find_ellipsoid(wkt, semi_major_axis, inverse_flattening, f7)
 
    found = f1 .and. f3 .and. f4 .and. f5 .and. f6 .and. f7

Here is the call graph for this function:

Here is the caller graph for this function:

wkt_to_cf_gridmapping()

character(len=linelength) function, public ncmodelexportmodule::wkt_to_cf_gridmapping

(

character(len=*), intent(in)

wkt

)

Accepts either WKT1 (OGC 01-009, PROJECTION["name"]) or WKT2 (ISO 19162:2019, METHOD["name"]) format. WKT1 is tried first; if PROJECTION[ is absent the WKT2 METHOD[ path is attempted. Returns '' for geographic CRS or unrecognised projection names. Only common groundwater projections included. Matching is case-insensitive.

Definition at line 573 of file NCModel.f90.

    use inputoutputmodule, only: upcase
    character(len=*), intent(in) :: wkt
    character(len=LINELENGTH) :: gmname
    character(len=LINELENGTH) :: proj_name
    character(len=LENBIGLINE) :: wkt_upper
    integer :: istart, iend
 
    gmname = ''
    proj_name = ''
 
    ! uppercase copy for case-insensitive keyword search
    wkt_upper = wkt
    call upcase(wkt_upper)
 
    ! --- WKT1 path: PROJECTION["name"] ---
    istart = index(wkt_upper, 'PROJECTION[')
    if (istart /= 0) then
      istart = istart + len('PROJECTION[')
      do while (istart <= len(wkt))
        if (wkt(istart:istart) == '"') exit
        istart = istart + 1
      end do
      if (istart > len(wkt)) return
      istart = istart + 1
      iend = index(wkt(istart:), '"')
      if (iend == 0) return
      iend = istart + iend - 2
      if (iend < istart) return
      proj_name = wkt(istart:iend)
      call upcase(proj_name)
      select case (trim(proj_name))
      case ('TRANSVERSE_MERCATOR')
        gmname = 'transverse_mercator'
      case ('LAMBERT_CONFORMAL_CONIC_2SP', 'LAMBERT_CONFORMAL_CONIC_1SP')
        gmname = 'lambert_conformal_conic'
      case ('ALBERS_CONIC_EQUAL_AREA')
        gmname = 'albers_conical_equal_area'
      case ('MERCATOR_1SP', 'MERCATOR_2SP')
        gmname = 'mercator'
      case ('POLAR_STEREOGRAPHIC')
        gmname = 'polar_stereographic'
      end select
      return
    end if
 
    ! --- WKT2 path: METHOD["name"] ---
    istart = index(wkt_upper, 'METHOD[')
    if (istart == 0) return
    istart = istart + len('METHOD[')
    do while (istart <= len(wkt))
      if (wkt(istart:istart) == '"') exit
      istart = istart + 1
    end do
    if (istart > len(wkt)) return
    istart = istart + 1
    iend = index(wkt(istart:), '"')
    if (iend == 0) return
    iend = istart + iend - 2
    if (iend < istart) return
    proj_name = wkt(istart:iend)
    call upcase(proj_name)
    select case (trim(proj_name))
    case ('TRANSVERSE MERCATOR')
      gmname = 'transverse_mercator'
    case ('LAMBERT CONIC CONFORMAL (2SP)', 'LAMBERT CONIC CONFORMAL (1SP)')
      gmname = 'lambert_conformal_conic'
    case ('ALBERS EQUAL AREA')
      gmname = 'albers_conical_equal_area'
    case ('MERCATOR (VARIANT A)', 'MERCATOR (VARIANT B)')
      gmname = 'mercator'
    case ('POLAR STEREOGRAPHIC (VARIANT A)', 'POLAR STEREOGRAPHIC (VARIANT B)')
      gmname = 'polar_stereographic'
    end select

Here is the call graph for this function:

Here is the caller graph for this function:

wkt_transverse_mercator_params()

subroutine, public ncmodelexportmodule::wkt_transverse_mercator_params	(	character(len=*), intent(in)	wkt,
		real(dp), intent(out)	longitude_of_central_meridian,
		real(dp), intent(out)	latitude_of_projection_origin,
		real(dp), intent(out)	scale_factor_at_central_meridian,
		real(dp), intent(out)	false_easting,
		real(dp), intent(out)	false_northing,
		real(dp), intent(out)	semi_major_axis,
		real(dp), intent(out)	inverse_flattening,
		logical(lgp), intent(out)	found
	)

Best-effort extraction from a WKT1 or WKT2 string – MF6 has no CRS parsing library. Tries WKT1 parameter names first, falling back to WKT2 names. found=.false. if any parameter cannot be located.

Definition at line 757 of file NCModel.f90.

    character(len=*), intent(in) :: wkt
    real(DP), intent(out) :: longitude_of_central_meridian
    real(DP), intent(out) :: latitude_of_projection_origin
    real(DP), intent(out) :: scale_factor_at_central_meridian
    real(DP), intent(out) :: false_easting
    real(DP), intent(out) :: false_northing
    real(DP), intent(out) :: semi_major_axis
    real(DP), intent(out) :: inverse_flattening
    logical(LGP), intent(out) :: found
    logical(LGP) :: f1, f2, f3, f4, f5, f6
 
    longitude_of_central_meridian = wkt_find_param(wkt, 'central_meridian', f1)
    if (.not. f1) longitude_of_central_meridian = &
      wkt_find_param(wkt, 'Longitude of natural origin', f1)
 
    latitude_of_projection_origin = wkt_find_param(wkt, 'latitude_of_origin', f2)
    if (.not. f2) latitude_of_projection_origin = &
      wkt_find_param(wkt, 'Latitude of natural origin', f2)
 
    scale_factor_at_central_meridian = wkt_find_param(wkt, 'scale_factor', f3)
    if (.not. f3) scale_factor_at_central_meridian = &
      wkt_find_param(wkt, 'Scale factor at natural origin', f3)
 
    false_easting = wkt_find_param(wkt, 'false_easting', f4)
    if (.not. f4) false_easting = wkt_find_param(wkt, 'False easting', f4)
 
    false_northing = wkt_find_param(wkt, 'false_northing', f5)
    if (.not. f5) false_northing = wkt_find_param(wkt, 'False northing', f5)
 
    call wkt_find_ellipsoid(wkt, semi_major_axis, inverse_flattening, f6)
 
    found = f1 .and. f2 .and. f3 .and. f4 .and. f5 .and. f6

Here is the call graph for this function:

Here is the caller graph for this function:

wrap_rotated_crs()

character(len=lenbigline) function, public ncmodelexportmodule::wrap_rotated_crs	(	character(len=*), intent(in)	wkt,
		real(dp), intent(in)	xorigin,
		real(dp), intent(in)	yorigin,
		real(dp), intent(in)	angrot
	)

Builds a DerivedProjectedCRS wrapping the given WKT2 PROJCRS in an EPSG:9624 (affine parametric transformation) deriving conversion parameterized from xorigin/yorigin/angrot (degrees). Per ISO 19111, a deriving conversion is directed base->derived, so the parameters encode the world->local (inverse) rotation; a CRS-aware consumer applies the inverse to resolve true position from local coordinates. Returns '' if wkt is not recognized as a WKT2 PROJCRS – a best-effort structural check only, as MF6 has no CRS parsing library and does not validate the WKT beyond this.

Definition at line 928 of file NCModel.f90.

    use inputoutputmodule, only: upcase
    character(len=*), intent(in) :: wkt
    real(DP), intent(in) :: xorigin, yorigin, angrot
    character(len=LENBIGLINE) :: derived_wkt
    character(len=LENBIGLINE) :: wkt_trim, wkt_upper, base_body, conversion
    character(len=30) :: a0s, a1s, a2s, b0s, b1s, b2s
    integer :: nstart, i, depth, closures, iend, wkt_len
    real(DP) :: ang, a0, a1, a2, b0, b1, b2
 
    derived_wkt = ''
 
    ! best-effort check that wkt is a WKT2 PROJCRS
    wkt_trim = trim(adjustl(wkt))
    wkt_len = len_trim(wkt_trim)
    wkt_upper = wkt_trim
    call upcase(wkt_upper)
    if (wkt_len < 9 .or. wkt_upper(1:8) /= 'PROJCRS[') return
 
    ! locate the end of PROJCRS's 2nd top-level child (the mandatory
    ! BASEGEOGCRS/BASEGEODCRS and CONVERSION elements, in that order),
    ! by tracking bracket depth from just inside the opening '['
    nstart = 9
    depth = 0
    closures = 0
    iend = 0
    do i = nstart, wkt_len
      if (wkt_trim(i:i) == '[') then
        depth = depth + 1
      else if (wkt_trim(i:i) == ']') then
        depth = depth - 1
        if (depth == 0) then
          closures = closures + 1
          if (closures == 2) then
            iend = i
            exit
          end if
        end if
      end if
    end do
    if (iend == 0) return
    base_body = wkt_trim(nstart:iend)
 
    ! EPSG:9624 affine parameters, base(world) -> derived(local grid)
    ang = angrot * dpio180
    a0 = -(xorigin * cos(ang) + yorigin * sin(ang))
    a1 = cos(ang)
    a2 = sin(ang)
    b0 = xorigin * sin(ang) - yorigin * cos(ang)
    b1 = -sin(ang)
    b2 = cos(ang)
    write (a0s, '(es22.15)') a0
    write (a1s, '(es22.15)') a1
    write (a2s, '(es22.15)') a2
    write (b0s, '(es22.15)') b0
    write (b1s, '(es22.15)') b1
    write (b2s, '(es22.15)') b2
 
    conversion = 'DERIVINGCONVERSION["MODFLOW 6 grid rotation",'// &
                 'METHOD["Affine parametric transformation",'// &
                 'ID["EPSG",9624]],'// &
                 'PARAMETER["A0",'//trim(adjustl(a0s))// &
                 ',LENGTHUNIT["metre",1]],'// &
                 'PARAMETER["A1",'//trim(adjustl(a1s))// &
                 ',SCALEUNIT["unity",1]],'// &
                 'PARAMETER["A2",'//trim(adjustl(a2s))// &
                 ',SCALEUNIT["unity",1]],'// &
                 'PARAMETER["B0",'//trim(adjustl(b0s))// &
                 ',LENGTHUNIT["metre",1]],'// &
                 'PARAMETER["B1",'//trim(adjustl(b1s))// &
                 ',SCALEUNIT["unity",1]],'// &
                 'PARAMETER["B2",'//trim(adjustl(b2s))// &
                 ',SCALEUNIT["unity",1]]]'
 
    derived_wkt = 'DERIVEDPROJCRS["MODFLOW 6 rotated grid CRS",'// &
                  'BASEPROJCRS['//trim(base_body)//'],'// &
                  trim(conversion)//','// &
                  'CS[Cartesian,2],'// &
                  'AXIS["easting (X)",east,ORDER[1],LENGTHUNIT["metre",1]],'// &
                  'AXIS["northing (Y)",north,ORDER[2],LENGTHUNIT["metre",1]]]'

Here is the call graph for this function:

Here is the caller graph for this function:

Variable Documentation

netcdf_mesh2d

@, public ncmodelexportmodule::netcdf_mesh2d = 2

Definition at line 38 of file NCModel.f90.

    ENUMERATOR :: NETCDF_MESH2D = 2 !< netcdf ugrid layered mesh export

netcdf_structured

@, public ncmodelexportmodule::netcdf_structured = 1

Definition at line 37 of file NCModel.f90.

    ENUMERATOR :: NETCDF_STRUCTURED = 1 !< netcdf structrured export

netcdf_undef

@, public ncmodelexportmodule::netcdf_undef = 0

Definition at line 36 of file NCModel.f90.

    ENUMERATOR :: NETCDF_UNDEF = 0 !< undefined netcdf export type