text:
      subroutine leq1s (a,n,b,m,ib,ijob,wa,ier)
!-------------------------------------------------------------
!
!   latest revision     - april 19, 1993
!
!   purpose             - linear equation solution -
!                           indefinite matrix - symmetric
!                           storage mode - space economizer
!                           solution
!
!   usage               - call leq1s (a,n,b,m,ib,ijob,wa,ier)
!
!   arguments    a      - input/output vector of dimension
!                           n*(n+1)/2. see parameter ijob.
!                n      - order of matrix a and the number of
!                           rows in b. (input)
!                b      - input/output matrix of dimension n
!                           by m. on input, b contains the m
!                           right-hand sides of the equation
!                           ax = b. on output, the solution
!                           matrix x replaces b. if ijob = 1,
!                           b is not used.
!                m      - number of right hand sides (columns
!                           in b). (input)
!                ib     - column dimension of b exactly as
!                           specified in the dimension
!                           statement in the calling program.
!                           (input)
!                ijob   - input option parameter. ijob = i
!                           implies when
!                           i = 0, factor the matrix a and
!                             solve the equation ax = b. on
!                             input, a contains the coefficient
!                             matrix of the equation ax = b,
!                             where a is assumed to be an
!                             n by n symmetric matrix. a is
!                             stored in symmetric storage mode
!                             and therefore has dimension
!                             n*(n+1)/2. on output, a is
!                             replaced by its factorized form.
!                           i = 1, factor the matrix a.
!                             a contains the same input/output
!                             information as if ijob = 0. b is
!                             not used.
!                           i = 2, solve the equation ax = b.
!                             this option implies that leq1s
!                             has already been called using
!                             ijob = 0 or 1 so that the matrix
!                             a has already been factored.
!                             in this case, output matrix a
!                             must have been saved for reuse
!                             in the call to leq1s.
!                wa     - work area of length n+1,
!                           on output wa(n+1) contains the
!                           number of eigenvalues below zero
!                ier    - error parameter. (output)
!                         warning error (with fix)
!                           ier = 65 indicates that matrix a
!                             is algorithmically singular.
!                         terminal error
!                           ier = 129 indicates that matrix a
!                             is exactly singular on one step
!                             of procedure.
!
!-------------------------------------------------------------
      implicit real (16) (a-h,o-z)
      dimension          a(*),b(ib,*),wa(*)
      real (16)     :: zero=0.0q0
      real (16)     :: alpha=0.6403882q0
      real (16)     :: sixtn=16.0q0
      save               zero,alpha,sixtn
!
!                                  first executable statement
!                                    initialize ier
!
      ier = 0
!
!   ****  first stage  ****        decompose a into the
!                                    product m*d*m-transpose
!                                    where m is unit lower
!                                    triangular and d is block
!                                    diagonal with blocks of
!                                    order 1 or 2.
!                                    m and d are written over a.
!
      if (n.le.0) goto 9005
      if (ijob.eq.2) goto 80
!
!                                  calculate equilibration
!                                    factors
!
      rn = sixtn*n
      wa(n+1) = 0
      ll = 1
      do 10 i=1,n
         wa(i) = zero
         l = ll
         do 5 j=1,n
            temp = abs(a(l))
            if (temp.gt.wa(i)) wa(i) = temp
            if (j.lt.i) then
               l = l+1
            else
               l = l+j
            end if
    5    continue
         wa(i) = wa(i)*rn
         ll = ll+i
   10 continue
      i = n
      if (i.eq.1) goto 75
!
!                                  ****  begin of the loop  ****
!
   15 im1 = i-1
      im2 = i-2
      ix = i*im1/2
      ixi = ix+i
!
!                                  calculate maximum off
!                                    diagonal element in row i
!
      aii = abs(a(ixi))
      save = zero
      do 20 l=1,im1
         temp = abs(a(ix+l))
         if (temp.gt.save) then
            save = temp
            j = l
         end if
   20 continue
      ichang = i
      if (aii.ge.alpha*save) goto 35
!
!                                  calculate maximum off
!                                    diagonal element in row j
!
      sigma = save
      jx = j*(j-1)/2
      jxj = jx+j
      jxt = jx+1
      ajj = abs(a(jxj))
      do 25 l=1,im1
         temp = abs(a(jxt))
         if (temp.gt.sigma.and.l.ne.j) sigma = temp
         if (l.lt.j) then
            jxt = jxt+1
         else
            jxt = jxt+l
         end if
   25 continue
!
!                                  choose a strategy for the
!                                    pivoting
!
      if (aii*sigma.ge.alpha*save*save) goto 35
!
      if (ajj.lt.alpha*sigma) goto 50
!
!                                  interchange rows i and j
!
      jxt = jx+1
      do 30 l=1,im1
         temp = a(ix+l)
         a(ix+l) = a(jxt)
         a(jxt) = temp
         if (l.lt.j) then
            jxt = jxt+1
         else
            jxt = jxt+l
         end if
   30 continue
      temp = a(jxj)
      a(jxj) = a(ixi)
      a(ixi) = a(ix+j)
      a(ix+j) = temp
      temp = wa(i)
      wa(i) = wa(j)
      wa(j) = temp
      aii = ajj
      ichang = j
!
!                                  we use a 1 by 1 pivot
!
   35 if (wa(i)+aii.le.wa(i)) ier = 65
      if (aii.eq.zero) goto 9005
      wa(i) = ichang
!
      aii = a(ixi)
      if (aii.lt.zero) wa(n+1) = wa(n+1)+1
      ixt = ix+1
      do 45 j=im1,1,-1
         save = -a(ix+j)/aii
         ixl = ixt-1
         ixt = ixt-j
         ish = ix-ixt+1
         mm = mod(j,5)
         if (mm.ne.0) then
            do 39 k=ixt,ixt+mm-1
               a(k) = a(k)+save*a(ish+k)
   39       continue
         end if
         ixf = ixt+mm
         do 40 k=ixf,ixl,5
            a(k) = a(k)+save*a(ish+k)
            a(k+1) = a(k+1)+save*a(ish+k+1)
            a(k+2) = a(k+2)+save*a(ish+k+2)
            a(k+3) = a(k+3)+save*a(ish+k+3)
            a(k+4) = a(k+4)+save*a(ish+k+4)
   40    continue
         a(ix+j) = save
   45 continue
      i = im1
      goto 75
!
!                                  we use a 2 by 2 pivot
!
   50 ix1 = ix-im1
      wa(n+1) = wa(n+1)+1
      if (j.eq.im1) goto 60
!
!                                  interchange rows i-1 and j
!
      jxt = jx+1
      do 55 l=1,im2
         temp = a(ix1+l)
         a(ix1+l) = a(jxt)
         a(jxt) = temp
         if (l.lt.j) then
            jxt = jxt+1
         else
            jxt = jxt+l
         end if
   55 continue
      temp = a(jxj)
      a(jxj) = a(ix)
      a(ix) = a(ix1+j)
      a(ix1+j) = temp
      temp = a(ix+j)
      a(ix+j) = a(ix+im1)
      a(ix+im1) = temp
      temp = wa(im1)
      wa(im1) = wa(j)
      wa(j) = temp
!
!                                  det must be negative
!
   60 det = a(ixi)*a(ix)-a(ixi-1)**2
!
      temp = max(wa(i),wa(im1))
      if (temp+abs(det).le.temp) ier = 65
      if (det.ge.zero) goto 9005
      wa(i) = j
      wa(im1) = det
!
      aim1i = a(ixi-1)/det
      aii = a(ixi)/det
      aim1 = a(ix)/det
      ixt = ix1
      do 70 j=im2,1,-1
         save = aim1i*a(ix1+j)-aim1*a(ix+j)
         temp = aim1i*a(ix+j)-aii*a(ix1+j)
         ixt = ixt-j
         do 65 k=1,j
            a(ixt+k) = a(ixt+k)+a(ix+k)*save+a(ix1+k)*temp
   65    continue
         a(ix+j) = save
         a(ix1+j) = temp
   70 continue
      i = im2
!
!                                  *** end of the loop ***
!
   75 if (i.gt.1) goto 15
      if (i.eq.1) then
         if (wa(1)+abs(a(1)).le.wa(1)) ier = 65
         if (a(1).eq.zero) goto 9005
         if (a(1).lt.zero) wa(n+1) = wa(n+1)+1
         wa(1) = 1
      end if
!
!   ****  second stage  ****       solve m*d*mt*x = b where
!                                    mt = m-transpose
!
   80 if (ijob.eq.1) go to 9000
!
      do 120 jc=1,m
!
!                                  solve m*d*y = b and store
!                                    y in b
!
         i = n
   85    if (i.le.1) goto 100
         im1 = i-1
         ix = i*im1/2
         ixi = ix+i
         ichang = wa(i)
         if (ichang.le.0) then
            write (*,'(2x,a)') 'error in routine leq1s |',              &
     &                         'array wa is wrong |'
            return
         end if
         save = b(ichang,jc)
         if (wa(im1).gt.zero) then
!
!                                  we use a 1 by 1 pivot
!
            b(ichang,jc) = b(i,jc)
            b(i,jc) = save/a(ixi)
            do 90 j=1,im1
               b(j,jc) = b(j,jc)+a(ix+j)*save
   90       continue
            i = im1
         else
!
!                                  we use a 2 by 2 pivot
!
            temp = b(i,jc)
            b(ichang,jc) = b(im1,jc)
            det = wa(im1)
            ix1 = ix-im1
            im2 = i-2
            b(i,jc) = (temp*a(ix)-save*a(ix+im1))/det
            b(im1,jc) = (save*a(ixi)-temp*a(ix+im1))/det
            do 95 j=1,im2
               b(j,jc) = b(j,jc)+a(ix+j)*temp+a(ix1+j)*save
   95       continue
            i = im2
         end if
         goto 85
  100    if (i.eq.1) then
            b(1,jc) = b(1,jc)/a(1)
            i = 2
         else
            i = 3
         end if
!
!                                  solve mt*x = y and store
!                                    x in b
!
  105    if (i.gt.n) goto 120
            ii = i
            if (wa(i).le.zero) ii = i+1
            im1 = i-1
            ix = i*im1/2
            do 115 k=i,ii
               save = b(k,jc)
               do 110 j=1,im1
                  save = save+a(ix+j)*b(j,jc)
  110          continue
               b(k,jc) = save
               ix = ix+k
  115       continue
            ichang = wa(ii)
            save = b(i,jc)
            b(i,jc) = b(ichang,jc)
            b(ichang,jc) = save
            i = ii+1
         goto 105
  120 continue
!
!                                  exit
!
 9000 if (ier.eq.0) return
      if (ier.eq.65) write (*,9065)
 9065 format(/5x,'warning error  ier=65  in  *** leq1s ***'/            &
     &       7x,'matrix a is algorithmically singular.'/)
      return
 9005 ier = 129
      write (*,9129)
 9129 format(/5x,'terminal error  ier=129  in  *** leq1s ***'/          &
     &       7x,'matrix a is exactly singular.'/)
      return
      end

max_stars_repo_name: DentonW/Ps-H-Scattering

max_issues_repo_name: DentonW/Ps-H-Scattering

max_forks_repo_name: DentonW/Ps-H-Scattering


text:
module read_raster_mod
    !! Read GDAL supported rasters. Both a type-based interface and a basic procedural interface are provided.
    !!
    !! The type-based interface (using the multi_raster_type) is more flexible and should
    !! be preferred. Given a set of raster filenames (in order of preference) it can
    !! extract raster data at a given set of coordinates. At each coordinates, it will use
    !! the highest-preference raster that has non-missing data. Constant and bilinear interpolation
    !! are supported. The code does not fully read the files into memory.
    !!
    !! The procedural interface is not recommended anymore, because it does not separate "opening and closing
    !! the files" from "reading the data". But see 'read_gdal_raster' and 'get_raster_dimensions'.
    !!
    use global_mod, only: charlen, ip, dp
    use logging_mod, only: log_output_unit
    use iso_c_binding

    implicit none

    private

    ! object based interface
    public:: gdal_raster_dataset_type, multi_raster_type
    ! procedural interface
    public:: read_gdal_raster, get_raster_dimensions
    ! tests
    public:: test_read_raster1

    ! Flag to use bilinear interpolation (1) or not (0) by default
    integer(C_INT), parameter :: use_bilinear_default = 1_C_INT

    type :: gdal_raster_dataset_type
        !!
        !! Type to work with a single raster dataset.
        !!

        type(c_ptr) :: hDataset !! Pointer to GDAL raster data structure
        character(charlen):: inputFile !! Filename
        logical:: isOpen = .FALSE. !! Record if the file is open
        integer(C_INT) :: xydim(2) !! x-y dimensions of the raster
        real(C_DOUBLE) :: lowerleft(2) !! Lower left (corner) of raster
        real(C_DOUBLE) :: upperright(2) !! Upper right (corner) of raster
        real(C_DOUBLE) :: adfGeoTransform(6) !! GDAL's geotransform information
        real(C_DOUBLE) :: dx(2) !! Raster cell size
        real(C_DOUBLE) :: nodata_value !! The raster's nodata value

        contains
        procedure:: initialise => initialise_gdal_raster_dataset !! Initialise the object
        procedure:: get_xy => get_xy_values !! Interpolate at coordinates
        procedure:: finalise => close_gdal_raster !! Close the object
        procedure:: print => print_summary_gdal_raster_dataset !! Print summary information

    end type

    type :: multi_raster_type
        !!
        !! Type to interpolate data from multiple rasters with different resoltuions and spatial extents.
        !!
        !! Often in applications we have multiple raster datasets for interpolation, and have an order in which
        !! we would like to use them. This type will manage that.
        !!
        !! **The raster_files should all have the same projection, but they can have different resolutions,
        !! extents, and they can overlap. Where the raster overlap, their preference order (for interpolation)
        !! corresponds to the order of the input file names.**
        !!

        character(len=charlen), allocatable :: raster_files(:) !! Raster filenames in order of preference
        type(gdal_raster_dataset_type), allocatable :: raster_datasets(:) !! Objects corresponding to each file

        real(C_DOUBLE) :: lowerleft(2) !! minimum lower left (corner) of all rasters combined
        real(C_DOUBLE) :: upperright(2) !! maximum upper right (corner) of all rasters combined

        contains
        procedure:: initialise => initialise_multi_raster !! Initialise the object
        procedure:: get_xy => get_xy_values_multi_raster !! Interpolate at coordinates
        procedure:: finalise => close_multi_raster !! Close the object

    end type

    interface

    !
    ! Here we interface many C routines
    !

    subroutine open_gdal_raster_Cfun(inputFile, hDataset) bind(C, name='open_gdal_raster')
        !! Open raster via GDAL
        use iso_c_binding
        implicit none
        character(kind=C_CHAR), intent(in) :: inputFile(*) !! Input filename
        type(c_ptr), intent(out) :: hDataset !! GDAL raster object
    end subroutine

    subroutine close_gdal_raster_Cfun(hDataset) bind(C, name='close_gdal_raster')
        !! Close the raster via GDAL
        use iso_c_binding
        implicit none
        type(c_ptr), intent(inout) :: hDataset !! GDAL raster object
    end subroutine

    subroutine get_gdal_raster_dimensions_Cfun(hDataset, xydim, lowerleft, upperright, &
        adfGeoTransform, dx, nodata_value) bind(C, name='get_gdal_raster_dimensions')
        !! Store the raster dimensions, lower-left and upper-right coordinates, cellsize, nodata value, and some GDAL geo-transform
        !! information
        use iso_c_binding
        implicit none
        type(C_PTR), intent(in) :: hDataset !! GDAL raster object
        integer(C_INT), intent(inout):: xydim(2) !! x/y dimensions
        real(C_DOUBLE), intent(inout):: lowerleft(2), upperright(2), adfGeoTransform(6), dx(2), nodata_value
    end subroutine

    subroutine get_values_at_xy_Cfun(hDataset, adfGeoTransform, x, y, z, N, verbose, bilinear, &
        band) bind(C, name='get_values_at_xy')
        !! Get the raster values at the given x/y coordinates
        use iso_c_binding
        implicit none
        type(C_PTR), intent(in) :: hDataset !! GDAL raster object
        integer(C_INT), value, intent(in):: N !! Number of x/y/z points
        real(C_DOUBLE) :: adfGeoTransform(6) !! GDAL geo-transform data
        real(C_DOUBLE) :: x(N), y(N), z(N)
        integer(C_INT), value, intent(in):: verbose, bilinear, band
    end subroutine

    subroutine read_gdal_raster_Cfun(inputFile, x, y, z, N, verbose, bilinear) bind(C, name='read_gdal_raster')
        !! Procedural interface for reading the raster. This routine opens the inputFile, reads the data and closes the file before
        !! returning.
        use iso_c_binding
        implicit none
        character(kind = C_CHAR), intent(in) :: inputFile(*) !! Raster file name
        integer(C_INT), value, intent(in) :: N, verbose, bilinear
        real(C_DOUBLE), intent(in) :: x(N), y(N)
        real(C_DOUBLE), intent(out) :: z(N)

    end subroutine

    subroutine get_raster_dimensions_Cfun(inputFile, xydim, lowerleft, upperright) bind(C, name='get_raster_dimensions')
        use iso_c_binding
        !! Procedural interface for reading raster dimensions, lower-left and upper-right coordinates.
        implicit none
        character(kind = C_CHAR), intent(in) :: inputFile(*)
        real(C_DOUBLE), intent(out):: lowerleft(2), upperright(2)
        integer(C_INT), intent(out):: xydim(2)

    end subroutine

    end interface

    interface read_gdal_raster
        !! read_gdal_raster is a procedural interface. Each call opens the file, extracts the data, and closes the file,
        !! whereas the object based interface (multi_raster_type & gdal_raster_dataset_type) can separate opening & closing from reading.
        module procedure read_gdal_raster_C_DOUBLE, read_gdal_raster_C_FLOAT, &
            read_gdal_raster_C_DOUBLE_2D, read_gdal_raster_C_FLOAT_2D
    end interface

    contains

    subroutine open_gdal_raster(inputFile, gdal_raster_dataset)
        character(len=charlen), intent(in) :: inputFile
        type(gdal_raster_dataset_type), intent(inout):: gdal_raster_dataset

        character(kind = C_CHAR, len=len_trim(inputFile) + 1) :: inputFile_c

        inputFile_c = trim(inputFile) // C_NULL_CHAR

        gdal_raster_dataset%inputFile = inputFile

        call open_gdal_raster_Cfun(inputFile_c, gdal_raster_dataset%hDataset)

        gdal_raster_dataset%isOpen = .TRUE.

    end subroutine

    subroutine close_gdal_raster(gdal_raster_dataset)
        !! Finalise routine for gdal_raster_dataset_type
        class(gdal_raster_dataset_type), intent(inout):: gdal_raster_dataset

        call close_gdal_raster_Cfun(gdal_raster_dataset%hDataset)
        gdal_raster_dataset%isOpen = .FALSE.

    end subroutine

    subroutine get_gdal_raster_dimensions(gdal_raster_dataset)
        type(gdal_raster_dataset_type), intent(inout):: gdal_raster_dataset

        call get_gdal_raster_dimensions_Cfun(&
            gdal_raster_dataset%hDataset,&
            gdal_raster_dataset%xydim, &
            gdal_raster_dataset%lowerleft, gdal_raster_dataset%upperright, &
            gdal_raster_dataset%adfGeoTransform, gdal_raster_dataset%dx, &
            gdal_raster_dataset%nodata_value)

    end subroutine

    subroutine initialise_gdal_raster_dataset(gdal_raster_dataset, inputFile)
        !! This is the main initialiser of gdal_raster_dataset_type
        class(gdal_raster_dataset_type), intent(inout):: gdal_raster_dataset
        character(len=charlen), intent(in) :: inputFile

        call open_gdal_raster(inputFile, gdal_raster_dataset)
        call get_gdal_raster_dimensions(gdal_raster_dataset)

    end subroutine

    subroutine get_xy_values(gdal_raster_dataset, x, y, z, N, verbose, bilinear, band)
        !! This is the key useful routine for gdal_raster_dataset_type
        class(gdal_raster_dataset_type), intent(in):: gdal_raster_dataset
        integer(ip), intent(in):: N
        real(dp), intent(in):: x(n), y(n)
        real(dp), intent(out):: z(n)
        integer(ip), intent(in), optional:: verbose
        integer(ip), intent(in), optional:: bilinear
        integer(ip), intent(in), optional:: band

        real(C_DOUBLE):: x_c(N), y_c(N), z_c(N)
        integer(C_INT):: verbose_c, bilinear_c, band_c, N_c

        ! Convert to C_DOUBLE
        x_c = real(x, C_DOUBLE)
        y_c = real(y, C_DOUBLE)
        z_c = x*0.0_C_DOUBLE
        N_c = N

        if(present(verbose)) then
            verbose_c = verbose
        else
            verbose_c = 0_C_INT
        end if

        if(present(bilinear)) then
            bilinear_c = bilinear
        else
            bilinear_c = use_bilinear_default
        end if

        if(present(band)) then
            band_c = band
        else
            band_c = 1_C_INT
        end if

        call get_values_at_xy_Cfun(gdal_raster_dataset%hDataset, &
            gdal_raster_dataset%adfGeoTransform, x_c, y_c, z_c, N_c, &
            verbose_c, bilinear_c, band_c)

        z = real(z_c, dp)

    end subroutine

    subroutine get_raster_dimensions(inputFile, xydim, lowerleft, upperright)
        !! Non-object-based interface
        !! Avoid this
        character(len=charlen), intent(in) :: inputFile
        integer(ip), intent(out):: xydim(2)
        real(dp), intent(out):: lowerleft(2), upperright(2)

        real(C_DOUBLE):: lowerleft_c(2), upperright_c(2)
        integer(C_INT):: xydim_c(2)
        character(kind = C_CHAR, len = len_trim(inputFile) + 1) :: inputFile_c

        inputFile_c = trim(inputFile) // C_NULL_CHAR

        call get_raster_dimensions_Cfun(inputFile_c, xydim_c, lowerleft_c, upperright_c)

        xydim = xydim_c
        lowerleft = real(lowerleft_c, dp)
        upperright = real(upperright_c, dp)

    end subroutine

    ! Read gdal raster, assuming x,y,z are rank 1 arrays with N entries of type C_DOUBLE
    subroutine read_gdal_raster_C_DOUBLE(inputFile, x, y, z, N, verbose, bilinear)
        character(len=charlen), intent(in) :: inputFile
        integer(ip), intent(in) :: N, verbose
        integer(ip), optional, intent(in) :: bilinear
        real(C_DOUBLE), intent(in) :: x(N), y(N)
        real(C_DOUBLE), intent(out) :: z(N)

        character(kind = C_CHAR, len = len_trim(inputFile) + 1) :: inputFile_c
        integer(C_INT) :: N_c, verbose_c, use_bilinear_local

        if(present(bilinear)) then
            use_bilinear_local = bilinear
        else
            use_bilinear_local = use_bilinear_default
        end if

        ! Use C types
        inputFile_c = trim(inputFile) // C_NULL_CHAR
        N_c = size(x, kind=C_INT)
        verbose_c = 1_C_INT * verbose

        call read_gdal_raster_Cfun(inputFile_c, x, y, z, N_c, verbose_c, use_bilinear_local)

    end subroutine

    ! Read gdal raster, assuming x,y,z are rank 2 arrays with N entries of type C_DOUBLE
    subroutine read_gdal_raster_C_DOUBLE_2D(inputFile, x, y, z, N, verbose, bilinear)
        character(len=charlen), intent(in) :: inputFile
        integer(ip), intent(in) :: N, verbose
        integer(ip), optional, intent(in) :: bilinear
        real(C_DOUBLE), intent(in) :: x(:,:), y(:,:)
        real(C_DOUBLE), intent(out) :: z(:,:)

        integer(C_INT) :: use_bilinear_local

        if(present(bilinear)) then
            use_bilinear_local = bilinear
        else
            use_bilinear_local = use_bilinear_default
        end if

        call read_gdal_raster_C_DOUBLE(inputFile, x, y, z, N, verbose, use_bilinear_local)

    end subroutine

    ! Read gdal raster, assuming x,y,z are rank 1 arrays with N entries of type C_FLOAT
    subroutine read_gdal_raster_C_FLOAT(inputFile, x, y, z, N, verbose, bilinear)
        character(len=charlen), intent(in) :: inputFile
        integer(ip), intent(in) :: N, verbose
        integer(ip), optional, intent(in) :: bilinear
        real(C_FLOAT), intent(in) :: x(N), y(N)
        real(C_FLOAT), intent(out) :: z(N)

        character(kind = C_CHAR, len = len_trim(inputFile) + 1) :: inputFile_c
        integer(C_INT) :: N_c, verbose_c, use_bilinear_local
        ! Copy x,y,z into these C_DOUBLE arrays to allow the routine to work
        ! with either single or double precision
        real(C_DOUBLE), allocatable :: x0(:), y0(:), z0(:)

        if(present(bilinear)) then
            use_bilinear_local = bilinear
        else
            use_bilinear_local = use_bilinear_default
        end if

        ! Use C types
        inputFile_c = trim(inputFile) // C_NULL_CHAR
        N_c = size(x, kind=C_INT)
        verbose_c = 1_C_INT * verbose

        ! If x,y,z are not C_DOUBLE, we copy C_DOUBLE versions
        allocate(x0(N_c), y0(N_c), z0(N_c))
        x0 = real(x, C_DOUBLE)
        y0 = real(y, C_DOUBLE)
        z0 = real(z, C_DOUBLE)
        call read_gdal_raster_Cfun(inputFile_c, x0, y0, z0, N_c, verbose_c, use_bilinear_local)
        z = real(z0, C_FLOAT)
        deallocate(x0, y0, z0)

    end subroutine

    ! Read gdal raster, assuming x,y,z are rank 2 arrays with N entries of type C_FLOAT
    subroutine read_gdal_raster_C_FLOAT_2D(inputFile, x, y, z, N, verbose, bilinear)
        character(len=charlen), intent(in) :: inputFile
        integer(ip), intent(in) :: N, verbose
        integer(ip), optional, intent(in) :: bilinear
        real(C_FLOAT), intent(in) :: x(:,:), y(:,:)
        real(C_FLOAT), intent(out) :: z(:,:)
        integer(C_INT) :: use_bilinear_local

        if(present(bilinear)) then
            use_bilinear_local = bilinear
        else
            use_bilinear_local = use_bilinear_default
        end if

        call read_gdal_raster_C_FLOAT(inputFile, x, y, z, N, verbose, use_bilinear_local)

    end subroutine

    subroutine initialise_multi_raster(multi_raster, raster_files)
        !!
        !! Read a number of rasters into a single multi_raster object. The
        !! rasters should all have the same projection, but can have different resolutions
        !! and extents.
        !!
        class(multi_raster_type), intent(inout) :: multi_raster
        character(len=charlen) :: raster_files(:)

        integer(ip) :: n, i

        if(allocated(multi_raster%raster_datasets)) then
            print*, 'multi_raster already allocated'
            stop
        end if

        n = size(raster_files, kind=ip)
        allocate(multi_raster%raster_files(n))
        multi_raster%raster_files = raster_files

        allocate(multi_raster%raster_datasets(n))

        do i = 1, n
            call multi_raster%raster_datasets(i)%initialise(raster_files(i))
        end do

        !! Make the lowerleft/upperright of the multiraster correspond to the extremes of the individual rasters
        multi_raster%lowerleft = [HUGE(1.0_C_DOUBLE) , HUGE(1.0_C_DOUBLE)]
        multi_raster%upperright = [-HUGE(1.0_C_DOUBLE), -HUGE(1.0_C_DOUBLE)]
        do i = 1, n
            multi_raster%lowerleft  = min(multi_raster%lowerleft , multi_raster%raster_datasets(i)%lowerleft )
            multi_raster%upperright = max(multi_raster%upperright, multi_raster%raster_datasets(i)%upperright)
        end do

    end subroutine

    subroutine close_multi_raster(multi_raster)
        !!
        !! Close/cleanup a multi raster object
        !!
        class(multi_raster_type), intent(inout) :: multi_raster

        integer(ip) :: i, n

        do i = 1, size(multi_raster%raster_datasets, kind=ip)
            call multi_raster%raster_datasets(i)%finalise()
        end do

        deallocate(multi_raster%raster_datasets)
        deallocate(multi_raster%raster_files)

    end subroutine

    subroutine get_xy_values_multi_raster(multi_raster, x, y, z, N, verbose, bilinear, band, na_below_limit)
        !!
        !! Interpolation from multi_raster
        !!
        class(multi_raster_type), intent(in) :: multi_raster !! a multi_raster object
        integer(ip), intent(in) :: N !! Number of x/y/z coordinates
        real(dp), intent(in) :: x(N), y(N) !! Input coordinates
        real(dp), intent(inout):: z(N) !! Output raster values at x/y
        integer(ip), optional, intent(in) :: verbose !! Print information
        integer(ip), optional, intent(in) :: bilinear !! Use bilinear interpolation (1) or nearest-cell interpolation (0)
        integer(ip), optional, intent(in) :: band !! Integer giving the raster band to interpolate from
        real(dp), optional, intent(in) :: na_below_limit
        !! Treat raster values below this number as NA. This can be useful if nodata values are not preserved exactly (e.g. due to
        !! changes in precision). Often nodata values are large negative numbers, in which case we can be sure that all numbers below
        !! some threshold (e.g. `na_below_limit = -1.0e+10`) should be treated as NA.

        real(dp) :: empty_value, ll(2), ur(2), lower_limit_l
        integer(ip) :: i, j, verbose_l, bilinear_l, band_l

        if(present(verbose)) then
            verbose_l = verbose
        else
            verbose_l = 0
        end if

        if(present(bilinear)) then
            bilinear_l = bilinear
        else
            bilinear_l = use_bilinear_default
        end if

        if(present(band)) then
            band_l = band
        else
            band_l = 1
        end if

        if(present(na_below_limit)) then
            lower_limit_l = na_below_limit
        else
            lower_limit_l = -HUGE(1.0_dp)
        end if

        ! Flag for unset 'z'
        empty_value = -huge(1.0_dp)
        z = empty_value

        ! Read rasters until no values are missing
        do j = 1, size(multi_raster%raster_datasets, kind=ip)

            ll = real(multi_raster%raster_datasets(j)%lowerleft, dp)
            ur = real(multi_raster%raster_datasets(j)%upperright, dp)

            do i = 1, N

                ! If z is populated with a value, we are done
                if(z(i) /= empty_value) cycle

                ! Read values inside the raster extent
                if(x(i) >= ll(1) .and. x(i) <= ur(1) .and. y(i) >= ll(2) .and. y(i) <= ur(2)) then
                    call multi_raster%raster_datasets(j)%get_xy(x(i), y(i), z(i), 1, verbose_l, bilinear_l, band_l)
                end if

                ! Set 'nodata' values back to empty values
                if(z(i) == real(multi_raster%raster_datasets(j)%nodata_value, dp)) z(i) = empty_value
                if(z(i) < lower_limit_l) z(i) = empty_value

            end do

        end do

    end subroutine

    subroutine print_summary_gdal_raster_dataset(raster_data)
        !! Print summary information about the raster data
        class(gdal_raster_dataset_type) :: raster_data

        write(log_output_unit, *) 'inputFile   : ', trim(raster_data%inputFile)
        write(log_output_unit, *) '  isOpen     : ', raster_data%isOpen
        write(log_output_unit, *) '  xydim      : ', raster_data%xydim
        write(log_output_unit, *) '  lowerleft  : ', raster_data%lowerleft
        write(log_output_unit, *) '  upperright : ', raster_data%upperright
        write(log_output_unit, *) '  dx         : ', raster_data%dx
        write(log_output_unit, *) '  nodata     : ', raster_data%nodata_value

    end subroutine

    subroutine test_read_raster1()
        !! Unit tests for read_raster_mod
        implicit none

        character(len=charlen) :: inputFile, inputFiles(2)
        integer(ip), parameter:: N = 10, verbose = 0
        integer(ip) :: i
        integer, parameter:: cdp = C_DOUBLE, csp = C_FLOAT
        real(cdp), allocatable :: x(:), y(:), z(:), real_z(:), real_z_bl(:)
        real(csp), allocatable :: xf(:), yf(:), zf(:), real_zf(:), real_zf_bl(:)

        integer(ip):: xydim(2)
        real(dp):: lowerleft(2), upperright(2)
        type(gdal_raster_dataset_type):: test_file
        type(multi_raster_type) :: test_multi_raster

        real(dp), allocatable:: x_dp(:), y_dp(:), z_dp(:)

        ! This requires that the file exists in the subdirectory 'data' where
        ! the program is run
        inputFile = "./data/test_rast.tif"

        allocate(x(N), y(N), z(N), real_z(N), real_z_bl(N), x_dp(N), y_dp(N), z_dp(N))
        allocate(xf(N), yf(N), zf(N), real_zf(N), real_zf_bl(N))

        ! Points in the raster we will look up
        x = [502523.122040_cdp, 506073.391951_cdp, 504937.667827_cdp, 501782.285213_cdp, 502901.704270_cdp, &
             505457.970585_cdp, 500256.751295_cdp, 502305.114871_cdp, 501575.432335_cdp, 502512.685946_cdp ]

        y = [1617967.96301_cdp, 1614589.30853_cdp, 1618488.11768_cdp, 1609423.47632_cdp, 1618325.86200_cdp, &
            1613200.33784_cdp, 1616185.67459_cdp, 1613888.10362_cdp, 1618757.31551_cdp, 1617615.75308_cdp]

        ! Values with 'nearest cell' interpolation
        real_z = [3725100, 3735700, 3729300, 3731900, 3725500, 3735700, 3722500, 3728700, 3722300, 3725300]

        ! Values with bilinear interpolation
        real_z_bl = [3725078.28107_cdp, 3735557.47537_cdp, 3729387.21797_cdp, 3732141.09411_cdp, 3725477.54654_cdp, &
                     3735715.60333_cdp, 3722327.82800_cdp, 3728722.12613_cdp, 3722393.54915_cdp, 3725409.61881_cdp]

        x_dp = real(x, dp)
        y_dp = real(y, dp)
        z_dp = real(x, dp)*0

        ! Single precision variants ('float')
        xf = real(x, csp)
        yf = real(y, csp)
        real_zf = real(real_z, csp)
        real_zf_bl = real(real_z_bl, csp)

        !
        !
        ! Test gdal raster dataset class
        !
        !

        call test_file%initialise(inputFile)

        if((test_file%inputFile == inputFile) .AND. &
            c_associated(test_file%hDataset) .AND. &
            test_file%isOpen) then
            print*, 'PASS'
        else
            print*, 'FAIL'

        end if

        if( (test_file%dx(1) == 200.0_C_DOUBLE) .AND. (test_file%dx(2) == -200.0_C_DOUBLE) .AND. &
            (test_file%xydim(1) == 50_C_INT) .AND. (test_file%xydim(2) == 50_C_INT) .AND. &
            (test_file%lowerleft(1) == 499000.0_C_DOUBLE) .AND. (test_file%lowerleft(2) == 1609000.0_C_DOUBLE) .AND. &
            (test_file%upperright(1) == 509000.0_C_DOUBLE) .AND. (test_file%upperright(2) == 1619000.0_C_DOUBLE)) then
            print*, 'PASS'
        else
            print*, 'FAIL'
            print*, test_file%dx
            print*, test_file%xydim
            print*, test_file%lowerleft
            print*, test_file%upperright
        end if

        call test_file%get_xy(x_dp, y_dp, z_dp, size(z_dp), verbose=0, bilinear=0, band=1)
        if(all(abs(z_dp - real_z) < (1.0e-6_dp * real_z))) then
            print*, 'PASS'
        else
            print*, 'FAIL'
        end if

        call test_file%get_xy(x_dp, y_dp, z_dp, size(z_dp), verbose=0, bilinear=1, band=1)
        if(all(abs(z_dp - real_z_bl) < (1.0e-6_dp*real_z_bl))) then
            print*, 'PASS'
        else
            print*, 'FAIL'
            print*, z_dp
            print*, real_z_bl
            print*, abs(z_dp - real_z_bl)
        end if

        call test_file%finalise()

        if( (test_file%isOpen .eqv. .FALSE.) .AND. &
            (c_associated(test_file%hDataset) .eqv. .FALSE.)) then
            print*, 'PASS'
        else
            print*, 'FAIL'
        end if

        !
        !
        ! Test reading of raster summary data
        !
        !
        call get_raster_dimensions(inputFile, xydim, lowerleft, upperright)
        if( (xydim(1) == 50) .and. (xydim(2) == 50) ) then
            print*, 'PASS'
        else
            print*, 'FAIL'
        end if

        if( (lowerleft(1) == 499000.0_dp) .and. (lowerleft(2) == 1609000.0_dp)) then
            print*, 'PASS'
        else
            print*, 'FAIL'
        end if

        if( (upperright(1) == 509000.0_dp) .and. (upperright(2) == 1619000.0_dp)) then
            print*, 'PASS'
        else
            print*, 'FAIL'
        end if

        !
        !
        ! Test of raster IO and bilinear interpolation
        !
        !

        ! Test double precision version
        call read_gdal_raster(inputFile, x, y, z, N, verbose, bilinear=0_ip)

        ! Got these values from R's extract
        if(all(abs(z - real_z) < 1.0e-10_cdp*real_z)) then
            print*, 'PASS'
        else
            print*, 'FAIL'
            do i = 1, N
                print*, real_z(i), z(i)
            end do
        end if

        ! Test bilinear version -- got these values from R's extract with
        ! method='bilinear'.
        call read_gdal_raster(inputFile, x, y, z, N, verbose, bilinear=1_ip)
        if(all(abs(z - real_z_bl) < 1.0e-10_cdp*real_z_bl)) then
            print*, 'PASS'
        else
            print*, 'FAIL'
            do i = 1, N
                print*, real_z_bl(i), z(i)
            end do
        end if

        ! Now try single precision version
        call read_gdal_raster(inputFile, xf, yf, zf, N, verbose, bilinear=0_ip)

        if(all(abs(zf - real_zf) < 1.0e-6_csp*real_zf)) then
            print*, 'PASS'
        else
            print*, 'FAIL'
            do i = 1, N
                print*, real_zf(i), zf(i)
            end do
        end if

        ! Test bilinaer single precision version
        CALL read_gdal_raster(inputFile, xf, yf, zf, N, verbose, bilinear=1_ip)

        if(all(abs(zf - real_zf_bl) < 1.0e-6_csp*real_zf_bl)) then
            print*, 'PASS'
        else
            print*, 'FAIL'
            do i = 1, N
                print*, real_zf_bl(i), zf(i)
            end do
        end if

        !
        ! Test reading a raster with NaN values
        !

        inputFile = 'data/test_rast_nans.tif'
        call test_file%initialise(inputFile)
        !print*, test_file%nodata_value
        if( abs(test_file%nodata_value + 1.7e+308_C_DOUBLE) < 1.0e-06*(1.7e+308_C_DOUBLE)) then
            print*, 'PASS'
        else
            print*, 'FAIL'
        end if

        ! Point right in a nodata cell. Check it is identified as nodata
        x_dp(1) = 508000.0_dp
        y_dp(1) = 1614000.0_dp
        call test_file%get_xy(x_dp(1:1), y_dp(1:1), z_dp(1:1), size(z_dp(1:1)), verbose=0, bilinear=1, band=1)
        if(z_dp(1) == real(test_file%nodata_value, dp)) then
            print*, 'PASS'
        else
            print*, 'FAIL'
        end if

        ! Point not exactly in a nodata cell. Check that it is still
        ! correctly set to 'nodata', rather than interpolated.
        x_dp(1) = 500901.0_dp
        y_dp(1) = 1617111.0_dp
        call test_file%get_xy(x_dp(1:1), y_dp(1:1), z_dp(1:1), size(z_dp(1:1)), verbose=0, bilinear=1, band=1)
        if(z_dp(1) == real(test_file%nodata_value, dp)) then
            print*, 'PASS'
        else
            print*, 'FAIL'
        end if

        call test_file%finalise()

        !
        ! Test of multi-raster. Make the first raster contain NaN values, which
        ! 'pass through' to the second raster
        !
        inputFiles(1) = 'data/test_rast_nans.tif'
        inputFiles(2) = 'data/test_rast.tif'

        ! Test values
        x_dp = real(x, dp)
        y_dp = real(y, dp)
        z_dp = real(x, dp)*0
        ! ... with one point in an nan region
        x_dp(1) = 500901.0_dp
        y_dp(1) = 1617111.0_dp
        real_z_bl(1) = 3722691.0_dp

        call test_multi_raster%initialise(inputFiles)
        call test_multi_raster%get_xy(x_dp, y_dp, z_dp, size(z_dp, kind=ip), verbose=0, bilinear=1)

        if(all(abs(z_dp - real_z_bl) < (1.0e-6_dp*real_z_bl))) then
            print*, 'PASS'
        else
            print*, 'FAIL'
            print*, z_dp
            print*, real_z_bl
            print*, abs(z_dp - real_z_bl)
        end if

        call test_multi_raster%finalise()

        if(allocated(test_multi_raster%raster_files) .or. &
            allocated(test_multi_raster%raster_datasets)) then
            print*, 'FAIL'
        else
            print*, 'PASS'
        end if

    end subroutine
end module read_raster_mod

max_stars_repo_name: KayaGeoAus/ptha

max_issues_repo_name: KayaGeoAus/ptha

max_forks_repo_name: KayaGeoAus/ptha


text:
cdis    Forecast Systems Laboratory
cdis    NOAA/OAR/ERL/FSL
cdis    325 Broadway
cdis    Boulder, CO     80303
cdis
cdis    Forecast Research Division
cdis    Local Analysis and Prediction Branch
cdis    LAPS
cdis
cdis    This software and its documentation are in the public domain and
cdis    are furnished "as is."  The United States government, its
cdis    instrumentalities, officers, employees, and agents make no
cdis    warranty, express or implied, as to the usefulness of the software
cdis    and documentation for any purpose.  They assume no responsibility
cdis    (1) for the use of the software and documentation; or (2) to provide
cdis     technical support to users.
cdis
cdis    Permission to use, copy, modify, and distribute this software is
cdis    hereby granted, provided that the entire disclaimer notice appears
cdis    in all copies.  All modifications to this software must be clearly
cdis    documented, and are solely the responsibility of the agent making
cdis    the modifications.  If significant modifications or enhancements
cdis    are made to this software, the FSL Software Policy Manager
cdis    (softwaremgr@fsl.noaa.gov) should be notified.
cdis
cdis
cdis
cdis
cdis
cdis
cdis

      program diff_test

      integer imax,jmax,kmax,nlvl

      include 'lapsparms.cmn'
c
       parameter (imax=125)
       parameter (jmax=105)
c      parameter (imax=NX_L_MAX)
c      parameter (jmax=NY_L_MAX)
      parameter (kmax=200)
      parameter (nlvl=200)

      INTEGER	I4TIME,		!I4time of data
     1		KDIM,		!K dimension of DATA array
     1		LVL(nlvl),        !Level of each field (4 digit max)
     1		LVL_AVAIL(nlvl),        !Level of each field (4 digit max)
     1		FLAG,
     1		ERROR(2),
     1		INDX,
     1		I,J,K,start,
     1		ISTATUS
C
      REAL	DATA1(imax,jmax,kmax)	!Raw data to be written
      REAL	DATA2(imax,jmax,kmax)	!Raw data to be written
C
      CHARACTER*200	DIR_in  !Directory to be written to
      CHARACTER*200	DIR_out !Directory to be written to
      CHARACTER*31	EXT		!File name ext (up to 31 chars)
      CHARACTER*3	VAR(nlvl) 	        !3 letter ID of each field
      CHARACTER*3	LAPS_VAR_AVAIL(nlvl) 	!3 letter ID of each field
      character*19      VAR_AVAIL(nlvl)
      CHARACTER*4	LVL_COORD(nlvl)	!Vertical coordinate for each field
      CHARACTER*10	UNITS(nlvl)	!units of each field
      CHARACTER*125	COMMENT1(nlvl)	!Comments for each field
      CHARACTER*125	COMMENT2(nlvl)	!Comments for each field
      CHARACTER*9	GTIME
      CHARACTER*91	FILE_NAME
      CHARACTER*17	ASCTIME

      CHARACTER*4	VERSION
      CHARACTER*131	MODEL 		!Meteorological model in file
      CHARACTER*131	ORIGIN		!Location where file was created
      CHARACTER*11    LAPS_DOM_FILE   !Name of domain file e.g. NEST7GRID
      LOGICAL         l_packed_data

      logical l_pass

      character*20 filename
      character*9 a9_time
      character*3 var_last

      var_last = '   '
      idiff_msg_flag = 0
      diff_max_all = 0.
      diff_max_all_rel = 0.
      n_files = 0
      ndiff_all = 0
      l_pass = .true.

      write(6,*)' Enter 1 if comparing different machines, 0 if same'
      read(5,*)machine_factor

      write(6,*)' filename?'
5     read(5,1)filename
1     format(a)

      if(filename(1:3) .eq. 'end')then
!         write(6,*)' end'
          goto 999
      endif

      read(5,1,err=999)dir_in
      read(5,1,err=999)dir_out

      if(filename(1:6) .eq. 'static')then
          khmax = 4
          var(1) = 'LAT'
          var(2) = 'LON'
          var(3) = 'AVG'
          var(4) = 'LDF'
          ext = 'nest7grid'
          call rd_laps_static(dir_in,ext,imax,jmax,khmax,var,units
     1                     ,comment,data1,grid_spacing_m,istatus)

          if(istatus .ne. 1)then
              write(6,*)' Bad status reading 1st static file'
              stop
          endif

          call rd_laps_static(dir_out,ext,imax,jmax,khmax,var,units
     1                     ,comment,data2,grid_spacing_m,istatus)

          if(istatus .ne. 1)then
              write(6,*)' Bad status reading 2nd static file'
              stop
          endif

          thresh_write_pair = .01
          thresh_count_diff = .01

          ihmax = imax
          jhmax = jmax

      else
          n1 = index(filename,'.')
          n2 = index(filename,' ')

          a9_time = filename(1:9)
          ext = filename(n1+1:n2-1)

          call downcase(ext,ext)
          call cv_asc_i4time(a9_time,i4time)

          if(dir_in(1:3) .eq. 'end')then
!             write(6,*)' end'
              goto 999
          endif

          do i = 1,nlvl
              var(i) = '   '
          enddo

          call READ_LAPS_HEADER(I4TIME,DIR_IN,EXT,IHMAX,JHMAX,KHMAX,
     1                         LAPS_DOM_FILE,ASCTIME,VERSION,
     1                         MODEL,ORIGIN,VAR,LVL,NUM_VARIABLES,
     1                         VAR_AVAIL,LAPS_VAR_AVAIL,NUM_LEVELS,
     1                         LVL_AVAIL,LVL_COORD,UNITS,
     1                         COMMENT1,L_PACKED_DATA,ISTATUS)

!         For this extension, set default values of:
!         thresh_write_pair, thresh_count_diff, or num_diff_field_thresh

          thresh_write_pair = 1e-05
          thresh_count_diff = 0.
          num_diff_field_thresh = 10

!         For this particular extension, set values of:
!         thresh_write_pair, thresh_count_diff, or num_diff_field_thresh

          if(ext(1:3) .eq. 'lc3')then
              do i = 1,42
                  var(i) = 'LC3'
                  lvl(i) = i
              enddo
              thresh_write_pair = .01
              thresh_count_diff = .01

          elseif(ext(1:3) .eq. 'lcp')then
              thresh_write_pair = .01
              thresh_count_diff = .01

          elseif(ext(1:3) .eq. 'lcv')then
              thresh_write_pair = .01
              thresh_count_diff = .01

          elseif(ext(1:3) .eq. 'lwc')then
              thresh_write_pair = .01
              thresh_count_diff = .01

          elseif(ext(1:3) .eq. 'lw3')then
              thresh_write_pair = .01
              thresh_count_diff = .1

          elseif(ext(1:3) .eq. 'lwm')then
              thresh_write_pair = .01
              thresh_count_diff = .1

          elseif(ext(1:3) .eq. 'lps')then
              thresh_write_pair = .01
              thresh_count_diff = .1

          elseif(ext(1:3) .eq. 'lvd')then
              thresh_write_pair = 1.
              thresh_count_diff = 1.

          elseif(ext(1:3) .eq. 'vrc')then
              thresh_write_pair = 0.1
              thresh_count_diff = 0.1

          elseif(ext(1:3) .eq. 'l1s')then
              do i = 1,4
                  lvl(i) = 0
              enddo
              var(1) = 'R01'
              var(2) = 'RTO'
              var(3) = 'S01'
              var(4) = 'STO'
              thresh_write_pair = .0001
              thresh_count_diff = .0001

          elseif(ext(1:3) .eq. 'lps')then
              do i = 1,21
                  lvl(i) = 1150 - 50 * i
              enddo

          elseif(ext(1:3) .eq. 'lsx')then
              thresh_write_pair = .01
              thresh_count_diff = .01
              num_diff_field_thresh = imax * jmax

          endif

!         Adjust values for this extension (dependent on machine) of:
!         thresh_count_diff and num_diff_field_thresh

          if(machine_factor .eq. 0)then
              thresh_count_diff = 0
              num_diff_field_thresh = 0
          endif

          if(var(1) .eq. 'RH')var(1) = 'LHE'

          write(6,*)' Reading: ',dir_in,a9_time,'.',ext(1:3)

          call read_laps_data(i4time,dir_in,ext,ihmax,jhmax,khmax,
     1             khmax,var,lvl,lvl_coord,units,comment1,data1,
     1             istatus)

          if(istatus .ne. 1)then
            if(ext(1:3) .ne. 'lsx' .and. ext(1:3) .ne. 'liw')then
              if(l_pass)then
                  write(6,*)' READ ERROR: OVERALL CRITERIA FAILURE'
                  l_pass = .false.
              endif
            endif
          endif

          write(6,*)' Reading: ',dir_out,a9_time,'.',ext(1:3)

          call read_laps_data(i4time,dir_out,ext,ihmax,jhmax,khmax,
     1             khmax,var,lvl,lvl_coord,units,comment2,data2,
     1             istatus)

          if(istatus .ne. 1)then
            if(ext(1:3) .ne. 'lsx' .and. ext(1:3) .ne. 'liw')then
              if(l_pass)then
                  write(6,*)' READ ERROR: OVERALL CRITERIA FAILURE'
                  l_pass = .false.
              endif
            endif
          endif

      endif ! static file

!       write(6,*)' Hit RETURN to CONTINUE'
!	read(5,*)

        thresh_write_pair = thresh_write_pair * machine_factor

        diff_max_file_rel = 0.
        diff_max_file = 0.
        diff_max_var = 0.
        ndiff_file = 0
        nvar = 1
        n_levels = 0

        do k = 1,khmax

!       Test whether we should switch variables within this file
        if(k .gt. 1)then
            if(var(k) .ne. var(k-1))then
                if(n_levels .gt. 1)then
                    write(6,*)' Max diff for variable ',var(k-1)(1:3)
     1                      ,' =',diff_max_var
                    write(6,*)
                    diff_max_var = 0.
                    nvar = nvar + 1
                endif
                n_levels = 0
            endif
        endif

!       Test if new variable
        if(var(k)(1:3) .ne. var_last)then

            write(6,*)
            write(6,*)' New var = ',var(k)(1:3)

!           For this variable, set values of:
!           thresh_write_pair, thresh_count_diff, or num_diff_field_thresh

            if    (ext(1:3) .eq. 'lvd')thresh_write_pair = 1.0

            if    (ext(1:3) .eq. 'lt1' .and. var(k)(1:2) .eq. 'HT')then
                thresh_count_diff = 2.0 * machine_factor
            elseif(ext(1:3) .eq. 'lt1' .and. var(k)(1:2) .eq. 'T3')then
                thresh_count_diff = .01 * machine_factor
            elseif(ext(1:3) .eq. 'lmt' .and. var(k)(1:3) .eq. 'LMT')then
                thresh_count_diff = 2.0 * machine_factor
            elseif(ext(1:3) .eq. 'lmt' .and. var(k)(1:3) .eq. 'LLR')then
                thresh_count_diff = 2.0 * machine_factor
            elseif(ext(1:3) .eq. 'lvd' .and. var(k)(1:3) .eq. 'SVN')then
                thresh_count_diff = 10.0 * machine_factor
            elseif(ext(1:3) .eq. 'lvd' .and. var(k)(1:3) .eq. 'ALB')then
                thresh_count_diff = 0.1 * machine_factor
                thresh_write_pair = 0.1                                   !JSmart addition 11-1-96
!           elseif(ext(1:3) .eq. 'lw3' .and. var(k)(1:2) .ne. 'OM')then
!               thresh_count_diff = .1 * machine_factor
            endif

            write(6,*)
     1      ' Threshold to write (first ten) grid point pairs = '
     1                                               ,thresh_write_pair
            write(6,*)
     1      ' Threshold to count lvl grid point differences   = '
     1                                               ,thresh_count_diff
            write(6,*)
     1      ' Max allowed count of lvl grid point differences = '
     1                                           ,num_diff_field_thresh
            write(6,*)

        endif

        n_levels = n_levels + 1

        diff_max_field = 0.
        diff_max_field_rel = 0.
        abs_value_max = 0.
        imaxd = 0
        jmaxd = 0
        iwrite = 0
        ndiff = 0
        inan = 0
        ndiff_msg = 0

	do i = 1,ihmax
           do j = 1,jhmax

c          if(
c!    1       data1(i,j,k) .le. r_min_normal()   .or.
c     1       data1(i,j,k) .ge. r_max_normal()   .or.
c!    1       data2(i,j,k) .le. r_min_normal()   .or.
c     1       data2(i,j,k) .ge. r_max_normal()
c     1                                                    )then
              if(isnan(data1(i,j,k)).ne.0 .or. isnan(data2(i,j,k)).ne.0)
     +             then
                 iwrite = iwrite + 1
                 if(iwrite .le. 10)then
                    write(6,21)i,j,k,' Nan'
                 endif
                 inan = inan + 1
              else
                 diff     = abs(data1(i,j,k)-data2(i,j,k))

!           Test if one of the points is missing and the other isn't
            if(   (data1(i,j,k) .eq. r_missing_data .or.
     1             data2(i,j,k) .eq. r_missing_data      )
     1                          .AND.
     1                     diff .gt. 0.                     )then

                ndiff_msg = ndiff_msg + 1
                idiff_msg_flag = 1

            else ! Both data points are non-missing

                diff_max_file = max(diff_max_file,diff)
                diff_max_var = max(diff_max_var,diff)
                if(diff .gt. diff_max_field)then
                    diff_max_field = diff
                    imaxd = i
                    jmaxd = j
                endif

            endif

            if(data1(i,j,k) .ne. r_missing_data)then
                abs_value_max = max(abs_value_max,abs(data1(i,j,k)))
            endif

            if(data2(i,j,k) .ne. r_missing_data)then
                abs_value_max = max(abs_value_max,abs(data2(i,j,k)))
            endif

            if(diff .gt. thresh_count_diff)then
                ndiff = ndiff + 1
                ndiff_file = ndiff_file + 1
                ndiff_all = ndiff_all + 1
            endif

            if(diff .gt. thresh_write_pair)then
                iwrite = iwrite + 1
                if(iwrite .le. 10)then
                    write(6,21,err=22)i,j,k,data1(i,j,k),data2(i,j,k)
     1                                                  ,diff
21                  format(1x,3i5,2f12.6,f12.6)
22              endif
            endif
          endif ! Nan test
        enddo ! j
        enddo ! i

        if(comment1(k) .ne. comment2(k))then
            write(6,*)comment1(k)(1:80)
            write(6,*)comment2(k)(1:80)
        endif
        if(inan .gt. 0)write(6,*)' # of Nans = ',inan

        if(abs_value_max .gt. 0.)then
!          if(ndiff_msg .eq. 0)then
               diff_max_field_rel = diff_max_field / abs_value_max
!          endif
        else
           diff_max_field_rel = 0.
        endif

        diff_max_file_rel  = max(diff_max_file_rel,diff_max_field_rel)

 	write(6,*)' df_mx - fld #',k,' ',var(k)
     1  ,lvl(k),' abs/rel/#',diff_max_field,diff_max_field_rel
     1  ,ndiff,imaxd,jmaxd

        if(k .eq. khmax .and. nvar .gt. 1
     1                  .and. n_levels .gt. 1)then
                write(6,*)
                write(6,*)' Max diff for variable ',var(k)(1:3),' ='
     1                      ,diff_max_var
                nvar = nvar + 1
        endif

        if(ndiff_msg .gt. 0)then
            write(6,*)
            write(6,*)' WARNING: # OF POINTS DIFFERING '
     1                ,'WRT MISSING DATA = ',
     1                  ndiff_msg
        endif

        if(ndiff + ndiff_msg .gt. num_diff_field_thresh)then
          if(l_pass)then
            write(6,*)' OVERALL CRITERIA FAILURE'
     1                           ,ndiff,ndiff_msg,num_diff_field_thresh
            l_pass = .false.
          endif
        endif

        write(6,*)

        var_last = var(k)(1:3)

        enddo ! k

      if (istatus .ne. 1) write (6,*)'Error in readlapsdata'

	   write(6,*)' OVERALL FILE diff_max (',ext(1:3)
     1      ,') [abs/rel/#] = ',diff_max_file,diff_max_file_rel
     1               ,ndiff_file
           write(6,*)
           n_files = n_files + 1
           diff_max_all     = max(diff_max_all,diff_max_file)
           diff_max_all_rel = max(diff_max_all_rel,diff_max_file_rel)
c           goto 5

999     continue

        if(n_files .gt. 1)then
          if(l_pass)then
             write(6,*)' MAX difference (all files)  [abs/rel/#] = '
     1              ,diff_max_all,diff_max_all_rel
     1              ,ndiff_all,' PASSED'
          else
             write(6,*)' MAX difference (all files)  [abs/rel/#] = '
     1              ,diff_max_all,diff_max_all_rel
     1              ,ndiff_all,' FAILED'
          endif
           write(6,*)
        else
          if(l_pass)then
            write(6,*)' PASSED'
            write(6,*)
          else
            write(6,*)' FAILED'
            write(6,*)
          endif
        endif

        if(idiff_msg_flag .eq. 1)then
            write(6,*)
            write(6,*)' WARNING: DIFFERENCES WRT MISSING DATA DETECTED'
            write(6,*)
        endif

        stop

      end

ext: f

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


module NumericalExchangeModule

  use KindModule,            only: DP, I4B
  use SimVariablesModule,    only: errmsg
  use BaseExchangeModule,    only: BaseExchangeType
  use NumericalModelModule,  only: NumericalModelType
  use BaseExchangeModule,    only: BaseExchangeType, AddBaseExchangeToList
  use ConstantsModule,       only: LINELENGTH, LENAUXNAME, DZERO
  use ListModule,            only: ListType
  use BlockParserModule,     only: BlockParserType

  implicit none

  private
  public :: NumericalExchangeType,                                             &
            AddNumericalExchangeToList, GetNumericalExchangeFromList

  type, extends(BaseExchangeType) :: NumericalExchangeType
    character(len=LINELENGTH), pointer              :: filename  => null()       !name of the input file
    character(len=7), pointer                       :: typename  => null()       !name of the type (e.g., 'NM-NM')
    logical, pointer                                :: implicit  => null()       !logical flag to indicate implicit or explict exchange
    integer(I4B), pointer                           :: iprpak    => null()       !print input flag
    integer(I4B), pointer                           :: iprflow   => null()       !print flag for cell by cell flows
    integer(I4B), pointer                           :: ipakcb    => null()       !save flag for cell by cell flows
    integer(I4B), pointer                           :: nexg      => null()       !number of exchanges
    integer(I4B), dimension(:), pointer, contiguous :: nodem1    => null()       !node numbers in model 1
    integer(I4B), dimension(:), pointer, contiguous :: nodem2    => null()       !node numbers in model 2
    real(DP), dimension(:), pointer, contiguous     :: cond      => null()       !conductance
    integer(I4B), dimension(:), pointer, contiguous :: idxglo    => null()       !pointer to solution amat for each connection
    integer(I4B), dimension(:), pointer, contiguous :: idxsymglo => null()       !pointer to symmetric amat position for each connection
    class(NumericalModelType), pointer              :: m1        => null()       !pointer to model 1
    class(NumericalModelType), pointer              :: m2        => null()       !pointer to model 2
    integer(I4B), pointer                           :: naux      => null()       !number of auxiliary variables
    character(len=LENAUXNAME), dimension(:), pointer,                         &
                                 contiguous :: auxname => null()                 !vector of auxname
    real(DP), dimension(:, :), pointer, contiguous  :: auxvar    => null()       !array of auxiliary variable values
    type(BlockParserType)                           :: parser                    !block parser
  contains
    procedure :: exg_df
    procedure :: exg_ac
    procedure :: exg_mc
    procedure :: exg_ar
    !procedure :: exg_rp (not needed yet; base exg_rp does nothing)
    procedure :: exg_ad
    procedure :: exg_cf
    procedure :: exg_fc
    procedure :: exg_nr
    procedure :: exg_cc
    procedure :: exg_cq
    procedure :: exg_bd
    procedure :: exg_ot
    procedure :: exg_da
    procedure :: allocate_scalars
    procedure :: allocate_arrays
    procedure :: read_options
    procedure :: read_dimensions
    procedure :: get_iasym
  end type NumericalExchangeType

contains

  subroutine exg_df(this)
! ******************************************************************************
! exg_df -- define the exchange
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    use BaseModelModule, only: BaseModelType
    use InputOutputModule, only: getunit, openfile
    ! -- dummy
    class(NumericalExchangeType) :: this
    ! -- local
! ------------------------------------------------------------------------------
    !
    !
    ! -- return
    return
  end subroutine exg_df

  subroutine exg_ac(this, sparse)
! ******************************************************************************
! exg_ac -- If an implicit exchange then add connections to sparse
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    use SparseModule, only:sparsematrix
    ! -- dummy
    class(NumericalExchangeType) :: this
    type(sparsematrix), intent(inout) :: sparse
    ! -- local
    integer(I4B) :: n, iglo, jglo
! ------------------------------------------------------------------------------
    !
    if(this%implicit) then
      do n = 1, this%nexg
        iglo = this%nodem1(n) + this%m1%moffset
        jglo = this%nodem2(n) + this%m2%moffset
        call sparse%addconnection(iglo, jglo, 1)
        call sparse%addconnection(jglo, iglo, 1)
      enddo
    endif
    !
    ! -- return
    return
  end subroutine exg_ac

  subroutine exg_mc(this, iasln, jasln)
! ******************************************************************************
! exg_mc -- Map the connections in the global matrix
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- module
    use SparseModule, only:sparsematrix
    ! -- dummy
    class(NumericalExchangeType) :: this
    integer(I4B), dimension(:), intent(in) :: iasln
    integer(I4B), dimension(:), intent(in) :: jasln
    ! -- local
    integer(I4B) :: n, iglo, jglo, ipos
! ------------------------------------------------------------------------------
    !
    if(this%implicit) then
      do n = 1, this%nexg
        iglo = this%nodem1(n)+this%m1%moffset
        jglo = this%nodem2(n)+this%m2%moffset
        ! -- find jglobal value in row iglo and store in idxglo
        do ipos = iasln(iglo), iasln(iglo + 1) - 1
          if(jglo == jasln(ipos)) then
            this%idxglo(n) = ipos
            exit
          endif
        enddo
        ! -- find and store symmetric location
        do ipos = iasln(jglo), iasln(jglo + 1) - 1
          if(iglo == jasln(ipos)) then
            this%idxsymglo(n) = ipos
            exit
          endif
        enddo
      enddo
    endif
    !
    ! -- Return
    return
  end subroutine exg_mc

  subroutine exg_ar(this)
! ******************************************************************************
! exg_ar -- Allocate and read
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    class(NumericalExchangeType) :: this
! ------------------------------------------------------------------------------
    !
    ! -- return
    return
  end subroutine exg_ar

  subroutine exg_ad(this)
! ******************************************************************************
! exg_ad -- Advance
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- dummy
    class(NumericalExchangeType) :: this
    ! -- local
! ------------------------------------------------------------------------------
    !
    !
    ! -- return
    return
  end subroutine exg_ad

  subroutine exg_cf(this, kiter)
! ******************************************************************************
! exg_cf -- Calculate conductance, and for explicit exchanges, set the
!   conductance in the boundary package.
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    ! -- dummy
    class(NumericalExchangeType) :: this
    integer(I4B),intent(in) :: kiter
    ! -- local
! ------------------------------------------------------------------------------
    !
    !
    ! -- return
    return
  end subroutine exg_cf

  subroutine exg_fc(this, kiter, iasln, amatsln, inwtflag)
! ******************************************************************************
! exg_fc -- Fill the matrix
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- dummy
    class(NumericalExchangeType) :: this
    integer(I4B), intent(in) :: kiter
    integer(I4B), dimension(:), intent(in) :: iasln
    real(DP), dimension(:), intent(inout) :: amatsln
    integer(I4B), optional, intent(in) :: inwtflag
    ! -- local
    integer(I4B) :: i, nodem1sln, nodem2sln, idiagsln
! ------------------------------------------------------------------------------
    !
    if(this%implicit) then
      do i = 1, this%nexg
        amatsln(this%idxglo(i)) = this%cond(i)
        amatsln(this%idxsymglo(i)) = this%cond(i)
        nodem1sln = this%nodem1(i) + this%m1%moffset
        nodem2sln = this%nodem2(i) + this%m2%moffset
        idiagsln = iasln(nodem1sln)
        amatsln(idiagsln) = amatsln(idiagsln) - this%cond(i)
        idiagsln = iasln(nodem2sln)
        amatsln(idiagsln) = amatsln(idiagsln) - this%cond(i)
      enddo
    else
      ! -- nothing to do here
    endif
    !
    ! -- return
    return
  end subroutine exg_fc

  subroutine exg_nr(this, kiter, iasln, amatsln, inwtflag)
! ******************************************************************************
! exg_nr -- Add Newton-Raphson terms to the solution
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- dummy
    class(NumericalExchangeType) :: this
    integer(I4B), intent(in) :: kiter
    integer(I4B), dimension(:), intent(in) :: iasln
    real(DP), dimension(:), intent(inout) :: amatsln
    integer(I4B), optional, intent(in) :: inwtflag
    ! -- local
! ------------------------------------------------------------------------------
    !
    ! -- return
    return
  end subroutine exg_nr

  subroutine exg_cc(this, icnvg)
! ******************************************************************************
! exg_cc -- Additional convergence check
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- dummy
    class(NumericalExchangeType) :: this
    integer(I4B), intent(inout) :: icnvg
    ! -- local
! ------------------------------------------------------------------------------
    !
    ! -- return
    return
  end subroutine exg_cc

  subroutine exg_cq(this, icnvg, isuppress_output, isolnid)
! ******************************************************************************
! exg_cq -- Calculate flow
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    use ConstantsModule, only: LENBUDTXT
    ! -- dummy
    class(NumericalExchangeType) :: this
    integer(I4B), intent(inout) :: icnvg
    integer(I4B), intent(in) :: isuppress_output
    integer(I4B), intent(in) :: isolnid
    ! -- local
! ------------------------------------------------------------------------------
    !
    ! -- return
    return
  end subroutine exg_cq

  subroutine exg_bd(this, icnvg, isuppress_output, isolnid)
! ******************************************************************************
! exg_bd -- Exchange budget
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    use ConstantsModule, only: LENBUDTXT
    ! -- dummy
    class(NumericalExchangeType) :: this
    integer(I4B), intent(inout) :: icnvg
    integer(I4B), intent(in) :: isuppress_output
    integer(I4B), intent(in) :: isolnid
    ! -- local
! ------------------------------------------------------------------------------
    !
    ! -- return
    return
  end subroutine exg_bd

  subroutine exg_ot(this)
! ******************************************************************************
! exg_ot
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    use SimVariablesModule, only: iout
    ! -- dummy
    class(NumericalExchangeType) :: this
    ! -- local
    integer(I4B) :: iexg, n1, n2
    real(DP) :: flow
    character(len=LINELENGTH) :: node1str, node2str
    ! -- format
    character(len=*), parameter :: fmtheader =                                 &
     "(/1x, 'SUMMARY OF EXCHANGE RATES FOR EXCHANGE ', a, ' WITH ID ', i0, /,  &
       &2a16, 4a16, /, 96('-'))"
    character(len=*), parameter :: fmtdata =                                   &
     "(2a16, 4(1pg16.6))"
! ------------------------------------------------------------------------------
    !
    ! -- Write a table of exchanges
    if(this%iprflow /= 0) then
      write(iout, fmtheader) trim(adjustl(this%name)), this%id, 'NODEM1',      &
                             'NODEM2', 'COND', 'X_M1', 'X_M2', 'FLOW'
      do iexg = 1, this%nexg
        n1 = this%nodem1(iexg)
        n2 = this%nodem2(iexg)
        flow = this%cond(iexg) * (this%m2%x(n2) - this%m1%x(n1))
        call this%m1%dis%noder_to_string(n1, node1str)
        call this%m2%dis%noder_to_string(n2, node2str)
        write(iout, fmtdata) trim(adjustl(node1str)), trim(adjustl(node2str)), &
                             this%cond(iexg), this%m1%x(n1), this%m2%x(n2),    &
                             flow
      enddo
    endif
    !
    ! -- return
    return
  end subroutine exg_ot

  subroutine allocate_scalars(this)
! ******************************************************************************
! allocate_scalars
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    use MemoryManagerModule, only: mem_allocate
    use ConstantsModule, only: LENORIGIN
    ! -- dummy
    class(NumericalExchangeType) :: this
    ! -- local
! ------------------------------------------------------------------------------
    !
    allocate(this%filename)
    allocate(this%typename)
    call mem_allocate(this%implicit, 'IMPLICIT', this%name)
    call mem_allocate(this%iprpak, 'IPRPAK', this%name)
    call mem_allocate(this%iprflow, 'IPRFLOW', this%name)
    call mem_allocate(this%ipakcb, 'IPAKCB', this%name)
    call mem_allocate(this%nexg, 'NEXG', this%name)
    call mem_allocate(this%naux, 'NAUX', this%name)
    this%filename = ''
    this%typename = ''
    this%implicit = .false.
    this%iprpak = 0
    this%iprflow = 0
    this%ipakcb = 0
    this%nexg = 0
    this%naux = 0
    !
    ! -- return
    return
  end subroutine allocate_scalars

  subroutine allocate_arrays(this)
! ******************************************************************************
! allocate_arrays
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    use MemoryManagerModule, only: mem_allocate
    use ConstantsModule, only: LENORIGIN
    ! -- dummy
    class(NumericalExchangeType) :: this
    ! -- local
    character(len=LENORIGIN) :: origin
! ------------------------------------------------------------------------------
    !
    ! -- create the origin name
    origin = trim(this%name)
    !
    call mem_allocate(this%nodem1, this%nexg, 'NODEM1', origin)
    call mem_allocate(this%nodem2, this%nexg, 'NODEM2', origin)
    call mem_allocate(this%cond, this%nexg, 'COND', origin)
    call mem_allocate(this%idxglo, this%nexg, 'IDXGLO', origin)
    call mem_allocate(this%idxsymglo, this%nexg, 'IDXSYMGLO', origin)
    call mem_allocate(this%auxvar, this%naux, this%nexg, 'AUXVAR', origin)
    !
    ! -- return
    return
  end subroutine allocate_arrays

  subroutine exg_da(this)
! ******************************************************************************
! exg_da
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    use MemoryManagerModule, only: mem_deallocate
    use ConstantsModule, only: LENORIGIN
    ! -- dummy
    class(NumericalExchangeType) :: this
    ! -- local
! ------------------------------------------------------------------------------
    !
    ! -- scalars
    deallocate(this%filename)
    deallocate(this%typename)
    call mem_deallocate(this%implicit)
    call mem_deallocate(this%iprpak)
    call mem_deallocate(this%iprflow)
    call mem_deallocate(this%ipakcb)
    call mem_deallocate(this%nexg)
    call mem_deallocate(this%naux)
    call mem_deallocate(this%auxname, 'AUXNAME', trim(this%name))
    !
    ! -- arrays
    call mem_deallocate(this%nodem1)
    call mem_deallocate(this%nodem2)
    call mem_deallocate(this%cond)
    call mem_deallocate(this%idxglo)
    call mem_deallocate(this%idxsymglo)
    call mem_deallocate(this%auxvar)
    !
    ! -- return
    return
  end subroutine exg_da

  subroutine read_options(this, iout)
! ******************************************************************************
! read_options -- Read Options
! Subroutine: (1) read options from input file
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    use ConstantsModule, only: LINELENGTH, LENAUXNAME
    use MemoryManagerModule, only: mem_allocate
    use SimModule, only: store_error, ustop
    use InputOutputModule, only: urdaux
    use ArrayHandlersModule, only: expandarray
    ! -- dummy
    class(NumericalExchangeType) :: this
    integer(I4B), intent(in) :: iout
    ! -- local
    character(len=:), allocatable :: line
    character(len=LINELENGTH) :: keyword
    character(len=LENAUXNAME), dimension(:), allocatable :: caux
    logical :: isfound
    logical :: endOfBlock
    integer(I4B) :: istart
    integer(I4B) :: istop
    integer(I4B) :: lloc
    integer(I4B) :: ierr
    integer(I4B) :: n
! ------------------------------------------------------------------------------
    !
    ! -- get options block
    call this%parser%GetBlock('OPTIONS', isfound, ierr,                        &
      supportOpenClose=.true., blockRequired=.false.)
    !
    ! -- parse options block if detected
    if (isfound) then
      write(iout,'(1x,a)')'PROCESSING EXCHANGE OPTIONS'
      do
        call this%parser%GetNextLine(endOfBlock)
        if (endOfBlock) then
          exit
        end if
        call this%parser%GetStringCaps(keyword)
        select case (keyword)
          case('AUX', 'AUXILIARY')
            call this%parser%GetRemainingLine(line)
            lloc = 1
            call urdaux(this%naux, this%parser%iuactive, iout, lloc, istart,     &
                        istop, caux, line, 'NM_NM_Exchange')
            call mem_allocate(this%auxname, LENAUXNAME, this%naux,               &
                                'AUXNAME', trim(this%name))
            do n = 1, this%naux
              this%auxname(n) = caux(n)
            end do
            deallocate(caux)
          case ('PRINT_INPUT')
            this%iprpak = 1
            write(iout,'(4x,a)') &
              'THE LIST OF EXCHANGES WILL BE PRINTED.'
          case ('PRINT_FLOWS')
            this%iprflow = 1
            write(iout,'(4x,a)') &
              'EXCHANGE FLOWS WILL BE PRINTED TO LIST FILES.'
          case default
            errmsg = "Unknown exchange option '" // trim(keyword) // "'."
            call store_error(errmsg)
            call this%parser%StoreErrorUnit()
            call ustop()
        end select
      end do
      write(iout,'(1x,a)') 'END OF EXCHANGE OPTIONS'
    end if
    !
    ! -- return
    return
  end subroutine read_options

  subroutine read_dimensions(this, iout)
! ******************************************************************************
! read_dimensions -- Read Dimensions
! Subroutine: (1) read dimensions (size of exchange list) from input file
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    use ConstantsModule, only: LINELENGTH
    use SimModule, only: store_error, ustop
    implicit none
    ! -- dummy
    class(NumericalExchangeType) :: this
    integer(I4B), intent(in) :: iout
    ! -- local
    character(len=LINELENGTH) :: keyword
    integer(I4B) :: ierr
    logical :: isfound, endOfBlock
! ------------------------------------------------------------------------------
    !
    ! -- get options block
    call this%parser%GetBlock('DIMENSIONS', isfound, ierr,                     &
      supportOpenClose=.true.)
    !
    ! -- parse options block if detected
    if (isfound) then
      write(iout,'(1x,a)') 'PROCESSING EXCHANGE DIMENSIONS'
      do
        call this%parser%GetNextLine(endOfBlock)
        if (endOfBlock) exit
        call this%parser%GetStringCaps(keyword)
        select case (keyword)
          case ('NEXG')
            this%nexg = this%parser%GetInteger()
            write(iout,'(4x,a,i0)') 'NEXG = ', this%nexg
          case default
            errmsg = "Unknown dimension '" // trim(keyword) // "'."
            call store_error(errmsg)
            call this%parser%StoreErrorUnit()
            call ustop()
        end select
      end do
      write(iout,'(1x,a)') 'END OF EXCHANGE DIMENSIONS'
    else
      call store_error('Required dimensions block not found.')
      call this%parser%StoreErrorUnit()
      call ustop()
    end if
    !
    ! -- return
    return
  end subroutine read_dimensions

  function get_iasym(this) result (iasym)
    class(NumericalExchangeType) :: this
    integer(I4B) :: iasym
    iasym = 0
  end function get_iasym

  function CastAsNumericalExchangeClass(obj) result (res)
    implicit none
    class(*), pointer, intent(inout) :: obj
    class(NumericalExchangeType), pointer :: res
    !
    res => null()
    if (.not. associated(obj)) return
    !
    select type (obj)
    class is (NumericalExchangeType)
      res => obj
    end select
    return
  end function CastAsNumericalExchangeClass

  subroutine AddNumericalExchangeToList(list, exchange)
    implicit none
    ! -- dummy
    type(ListType),       intent(inout) :: list
    class(NumericalExchangeType), pointer, intent(in) :: exchange
    ! -- local
    class(*), pointer :: obj
    !
    obj => exchange
    call list%Add(obj)
    !
    return
  end subroutine AddNumericalExchangeToList

  function GetNumericalExchangeFromList(list, idx) result (res)
    implicit none
    ! -- dummy
    type(ListType),            intent(inout) :: list
    integer(I4B),                   intent(in)    :: idx
    class(NumericalExchangeType), pointer    :: res
    ! -- local
    class(*), pointer :: obj
    !
    obj => list%GetItem(idx)
    res => CastAsNumericalExchangeClass(obj)
    !
    return
  end function GetNumericalExchangeFromList

end module NumericalExchangeModule