text:
module m_KdTree
  !!# KdTree
  !!Build and search k-dimensional trees in 2, 3, and K dimensions.
  !!This KdTree is balanced, in that splits are made along the dimension with the largest variance.
  !!A quickselect is used to quickly find the median in each splitting dimension as the splitting value.
  !!The ends of each branch contain multiple leaves to prevent tail recursion.
  !!An in-depth example is given below on how to use all the aspects of the KdTree and KdTreeSearch classes.
  !!
  !!Important: Once a tree has been built with a set, do not change their values.  The KdTree does NOT make
  !!a copy of the input values used to build it.
  !!Important: Generating the tree does not modify the values used to build it.
  !!
  !!## Building the KdTree
  !!The KdTree object can be initialized on assignment by entering point co-ordinates.
  !!To build a 2-D tree, you can use two 1-D arrays as the x, y, co-ordinates, and optionally a third 1-D array
  !!to build a 3-D tree.
  !!```fortran
  !!use m_KdTree
  !!type(KdTree) :: tree
  !!tree = KdTree(x, y, [z])
  !!```
  !!Or you can build a k-dimensional tree using a 2-D array, where the first dimension is the number of items,
  !!and the second is the number of dimensions, k.
  !!```fortran
  !!use m_KdTree
  !!type(KdTree) :: tree
  !!tree = KdTree(D)
  !!```
  !!
  !!## Querying the KdTree
  !!After the tree is initialized, a search class can be used to perform search for the nearest neighbour,
  !!the k nearest neighbours, all neighbours within a radius, k nearest within a radius, and items within
  !!upper and lower bounds in each dimension.
  !!The searches are thread safe and can be used in a parallel region.
  !!
  !!After the KdTree is built, various queries can be carried out.  Searches that return multiple values
  !!are called using the argDynamicArray within coretran. The [[dArgDynamicArray_Class]] contains an integer index
  !!containing the indices into the co-ordinates that are closest, and a double precision that contains the
  !!distance from the query point to those points.
  !!
  !!```fortran
  !!use dArgDynamicArray_Class
  !!use m_KdTree
  !!type(KdTreeSearch) :: search
  !!integer(i32) :: i
  !!type(dArgDynamicArray) :: da
  !!! Nearest Neighbour to (0, 0), for 3D add z, and zQuery
  !!i = search%nearest(tree, x, y, [z], xQuery = 0.d0, yQuery = 0.d0, [zQuery = 0.d0])
  !!! K-Nearest to (0, 0), for 3D add z, and zQuery
  !!da = search%kNearest(tree, x, y, [z], xQuery = 0.d0, yQuery = 0.d0, [zQuery = 0.d0], k = 10)
  !!! Search for all points within a given distance
  !!da = search%kNearest(tree, x, y, [z], xQuery = 0.d0, yQuery = 0.d0, [zQuery = 0.d0], radius = 10.d0)
  !!! Search for all k points within a given distance
  !!da = search%kNearest(tree, x, y, [z], xQuery = 0.d0, yQuery = 0.d0, [zQuery = 0.d0], k = 10, radius = 10.d0)
  !!```
  !!
  !!## Full Example
  !!```fortran
  !!program kdTree_test
  !!use variableKind, only: i32, r64
  !!use m_allocate, only: allocate
  !!use m_deallocate, only: deallocate
  !!use m_random, only: rngNormal
  !!use m_KdTree, only: KdTree, KdTreeSearch
  !!use m_dArgDynamicArray, only: dArgDynamicArray
  !!use m_string, only: str
  !!implicit none
  !!real(r64), allocatable :: x(:), y(:), z(:), D(:,:)
  !!integer(i32) :: ia, N
  !!type(KdTree) :: tree
  !!type(KdSearch) :: search
  !!type(dArgDynamicArray) :: da
  !!
  !!!====================================================================!
  !!! 2D KdTree example
  !!!====================================================================!
  !!! Create some random points in space
  !!N = 1d6
  !!call allocate(x, N)
  !!call allocate(y, N)
  !!call rngNormal(x)
  !!call rngNormal(y)
  !!! Build the tree
  !!tree = KdTree(x, y)
  !!! Get the nearest neighbour to (0, 0)
  !!ia = search%kNearest(tree, x, y, xQuery = 0.d0, yQuery = 0.d0)
  !!write(*,'(a)') 'Nearest point to the query location: '//str(x(ia))//str(y(ia))
  !!! Get the 10 nearest neighbours to the query
  !!da = search%kNearest(tree, x, y, xQuery = 0.d0, yQuery = 0.d0, k = 10)
  !!write(*,'(a)') 'The 10 nearest neighbour indices and distances:'
  !!call da%print()
  !!! Find all the points within a 1.d0
  !!da = search%kNearest(tree, x, y, xQuery = 0.d0, yQuery = 0.d0, radius = 1.d0)
  !!write(*,'(a)') 'The points within a distance of 1.d0'
  !!call da%print()
  !!Deallocate any tree memory
  !!call tree%deallocate()
  !!
  !!!====================================================================!
  !!! 3D KdTree example
  !!!====================================================================!
  !!! Create the third dimension
  !!call allocate(z, N)
  !!call rngNormal(z)
  !!! Build the tree
  !!tree = KdTree(x, y, z)
  !!! Get the nearest neighbour to (0, 0, 0)
  !!ia = search%kNearest(tree, x, y, z, xQuery = 0.d0, yQuery = 0.d0, zQuery = 0.d0)
  !!write(*,'(a)')'Nearest point to the query location: '//str(x(ia))//str(y(ia))//str(z(ia))
  !!! Get the 10 nearest neighbours to the query
  !!da = search%kNearest(tree, x, y, z, xQuery = 0.d0, yQuery = 0.d0, zQuery = 0.d0, k = 10)
  !!write(*,'(a)') 'The 10 nearest neighbour indices and distances:'
  !!call da%print()
  !!! Find all the points within a 1.d0
  !!da = search%kNearest(tree, x, y, z, xQuery = 0.d0, yQuery = 0.d0, zQuery = 0.d0, radius = 1.d0)
  !!write(*,'(a)') 'The points within a distance of 1.d0'
  !!call da%print()
  !!Deallocate any tree memory
  !!call tree%deallocate()
  !!
  !!!====================================================================!
  !!! KD KdTree example
  !!!====================================================================!
  !!call allocate(D, [N, 3])
  !!D(:,1) = x
  !!D(:,2) = y
  !!D(:,3) = z
  !!! Build the tree
  !!tree = KdTree(D)
  !!! Get the nearest neighbour to (0, 0, 0)
  !!ia = search%kNearest(tree, D, query = [0.d0, 0.d0, 0.d0])
  !!write(*,'(a)')'Nearest point to the query location: '//str(D(ia, 1))//str(D(ia, 2))//str(D(ia, 3))
  !!! Get the 10 nearest neighbours to the query
  !!da = search%kNearest(tree, D, query = [0.d0, 0.d0, 0.d0], k = 10)
  !!write(*,'(a)') 'The 10 nearest neighbour indices and distances:'
  !!call da%print()
  !!! Find all the points within a 1.d0
  !!da = search%kNearest(tree, D,  query = [0.d0, 0.d0, 0.d0], radius = 1.d0)
  !!write(*,'(a)') 'The points within a distance of 1.d0'
  !!call da%print()
  !!Deallocate any tree memory
  !!call tree%deallocate()
  !!call deallocate(x)
  !!call deallocate(y)
  !!call deallocate(z)
  !!call deallocate(D)
  !!end program
  !!```
  use variableKind, only: i32, r64
  use m_allocate, only: allocate
  use m_array1D, only: arange
  use m_deallocate, only: deallocate
  use m_errors, only: eMsg
  use m_maths, only: variance
  use iDynamicArray_Class, only: iDynamicArray
  use dArgDynamicArray_Class, only: dArgDynamicArray
  use m_select, only: argSelect
  use m_strings, only: str

  implicit none

  private

  public :: KdTree
  public :: KdTreeSearch

  !====================================================================!
  ! Type Definitions
  !====================================================================!
  !====================================================================!
  type, private :: KdTreebranch
    private
    !! A branch of the tree, could be a leaf, or branch
    integer(i32) :: splitAlong  !index of median in dat d
    real(r64) :: median  !median value in dat x
    integer(i32) :: left, right ! left and right indices
    type(KdTreeBranch), pointer :: buds(:) => null()
  contains
    procedure :: init => init_branch
      !! KdTreeBranch%init() - Initialize the class
    procedure :: deallocate => deallocate_branch
      !! KdTreeBranch%deallocate() - deallocate the branch
  end type
  !====================================================================!
  !====================================================================!
  type :: KdTree
    !!KdTree in 2, 3, or N dimensions.  See [[m_KdTree]] for more information on how to use this class.
    private
    type(KdTreebranch) :: trunk
    integer(i32),allocatable :: indx(:)
    integer(i32) :: nDims
    integer(i32) :: N
    logical(i32) :: set =.false.
  contains
    procedure, public :: deallocate => deallocate_KdTree
      !! kdTree%deallocate() - deallocate the recursive pointers
  end type
  !====================================================================!
  !====================================================================!
  type :: KdTreeSearch
    !!Class to search a KdTree.  See [[m_KdTree]] for more information on how to use this class.
  contains
    generic, public :: nearest => nearest2D, nearest3D, nearestKD
      !! KdTreeSearch%nearest() - Perform a nearest neighbour search
    procedure, private :: nearest2D => nearest_2D
      !! Overloaded type bound procedure with KdTreeSearch%nearest()
    procedure, private :: nearest3D => nearest_3D
      !! Overloaded type bound procedure with KdTreeSearch%nearest()
    procedure, private :: nearestKD => nearest_KD
      !! Overloaded type bound procedure with KdTreeSearch%nearest()
    generic, public :: kNearest => kNearest2D, kNearest3D, kNearestKD
      !! KdTreeSearch%kNearest() - Perform a k nearest neighbour search or a radius search.
    procedure, private :: kNearest2D => kNearest_2D
      !! Overloaded type bound procedure with KdTreeSearch%kNearest()
    procedure, private :: kNearest3D => kNearest_3D
      !! Overloaded type bound procedure with KdTreeSearch%kNearest()
    procedure, private :: kNearestKD => kNearest_KD
      !! Overloaded type bound procedure with KdTreeSearch%kNearest()
    generic, public :: rangeSearch => rangeSearch2D, rangeSearch3D, rangeSearchKD
      !! KdTreeSearch%rangeSearch() - Find all points within axis aligned lower and upper bounds
    procedure, private :: rangeSearch2D => rangeSearch_2D
      !! Overloaded type bound procedure with KdTreeSearch%rangeSearch()
    procedure, private :: rangeSearch3D => rangeSearch_3D
      !! Overloaded type bound procedure with KdTreeSearch%rangeSearch()
    procedure, private :: rangeSearchKD => rangeSearch_KD
      !! Overloaded type bound procedure with KdTreeSearch%rangeSearch()
  end type
  !====================================================================!

  !====================================================================!
  ! Interface definitions
  !====================================================================!
  interface
    !====================================================================!
    module subroutine init_branch(this, left,right)
      !! Overloaded Type bound procedure KdTreeBranch%init()
    !====================================================================!
    class(KdTreebranch) :: this
      !! KdTreeBranch class
    integer(i32), intent(in) :: left
      !! Left index, first call normally = 1
    integer(i32), intent(in) :: right
      !! Right index, first call normally = size of coordinates
    end subroutine
    !====================================================================!
    !====================================================================!
    module recursive subroutine deallocate_branch(this)
    !====================================================================!
    class(KdTreebranch) :: this
      !! KdTreeBranch class
    end subroutine
    !====================================================================!
  end interface

  interface KdTree
    !!Overloaded Initializer for a KdTree.
    !!
    !!Can be used to create a 2D, 3D, or ND, KdTree class.
    !!
    !!See [[m_KdTree]] for more information on how to use this class
    !====================================================================!
    module function init2D_KdTree(x, y) result(this)
      !! Overloaded by interface [[KdTree(type)]]
    !====================================================================!
    type(kdTree) :: this
      !! KdTree Class
    real(r64),intent(in) :: x(:)
      !! x-coordinates of the points
    real(r64),intent(in) :: y(:)
      !! y-coordinates of the points
    end function
    !====================================================================!
    !====================================================================!
    module function init3D_KdTree(x, y, z) result(this)
      !! Overloaded by interface [[KdTree(type)]]
    !====================================================================!
    type(kdTree) :: this
      !! KdTree Class
    real(r64),intent(in) :: x(:)
      !! x-coordinates of the points
    real(r64),intent(in) :: y(:)
      !! y-coordinates of the points
    real(r64),intent(in) :: z(:)
      !! z-coordinates of the points
    end function
    !====================================================================!
    !====================================================================!
    module function initKD_KdTree(D) result(this)
      !! Overloaded by interface [[KdTree(type)]]
    !====================================================================!
    type(kdTree) :: this
      !! KdTree Class
    real(r64),intent(in) :: D(:,:)
      !! Coordinates of the points, the k columns contain the k dimensional values.
    end function
    !====================================================================!
  end interface
  interface
    !====================================================================!
    module subroutine deallocate_KdTree(this)
      !! Overloaded Type bound procedure KdTree%deallocate()
    !====================================================================!
    class(KdTree) :: this
      !! KdTree Class
    end subroutine
    !====================================================================!
  end interface

interface
    !====================================================================!
    module function nearest_2D(search, tree, x, y, xQuery, yQuery) result(nearest)
      !! Overloaded Type bound procedure KdTreeSearch%nearest()
    !====================================================================!
    class(KdTreeSearch), intent(inout) :: search
      !! KdTreeSearch class
    class(kdTree),intent(in) :: tree
      !! KdTree class
    real(r64),intent(in) :: x(:)
      !! x co-ordinates of the points
    real(r64),intent(in) :: y(:)
      !! y co-ordinates of the points
    real(r64),intent(in) :: xQuery
      !! x co-ordinate of the query location
    real(r64),intent(in) :: yQuery
      !! y co-ordinate of the query location
    integer(i32) :: nearest
      !! Index of the nearest point to the query location
    end function
    !====================================================================!
    !====================================================================!
    module function nearest_3D(search, tree, x, y, z, xQuery, yQuery, zQuery) result(nearest)
      !! Overloaded Type bound procedure KdTreeSearch%nearest()
    !====================================================================!
    class(KdTreeSearch), intent(inout) :: search
      !! KdTreeSearch class
    class(kdTree),intent(in) :: tree
      !! KdTree class
    real(r64),intent(in) :: x(:)
      !! x co-ordinates of the points
    real(r64),intent(in) :: y(:)
      !! y co-ordinates of the points
    real(r64),intent(in) :: z(:)
      !! z co-ordinates of the points
    real(r64),intent(in) :: xQuery
      !! x co-ordinate of the query location
    real(r64),intent(in) :: yQuery
      !! y co-ordinate of the query location
    real(r64),intent(in) :: zQuery
      !! z co-ordinate of the query location
    integer(i32) :: nearest
      !! Index of the nearest point to the query location
    end function
    !====================================================================!
    !====================================================================!
    module function nearest_KD(search, tree, D, query) result(nearest)
      !! Overloaded Type bound procedure KdTreeSearch%nearest()
    !====================================================================!
    class(KdTreeSearch), intent(inout) :: search
      !! KdTreeSearch class
    class(kdTree),intent(in) :: tree
      !! KdTree class
    real(r64),intent(in) :: D(:,:)
      !! Co-ordinates of the points, the second dimension contains the k dimensions
    real(r64),intent(in) :: query(:)
      !! C-ordinate of the query location
    integer(i32) :: nearest
      !! Index of the nearest point to the query location
    end function
    !====================================================================!
  end interface

  interface
    !====================================================================!
    module function kNearest_2D(search, tree, x, y, xQuery, yQuery, k, radius) result(kNearest)
      !! Overloaded Type bound procedure KdTreeSearch%kNearest()
    !====================================================================!
    class(KdTreeSearch), intent(inout) :: search
      !! KdTreeSearch class
    class(kdTree),intent(in) :: tree
      !! KdTree class
    real(r64),intent(in) :: x(:)
      !! x co-ordinates of the points
    real(r64),intent(in) :: y(:)
      !! y co-ordinates of the points
    real(r64),intent(in) :: xQuery
      !! x co-ordinate of the query location
    real(r64),intent(in) :: yQuery
      !! y co-ordinate of the query location
    integer(i32), intent(in), optional :: k
      !! Number of points to find that are closest to the query
    real(r64), intent(in), optional :: radius
      !! Only find neighbours within this distance from the query
    type(dArgDynamicArray) :: kNearest
      !! Indices of the nearest points to the query location
    end function
    !====================================================================!
    !====================================================================!
    module function kNearest_3D(search, tree, x, y, z, xQuery, yQuery, zQuery, k, radius) result(kNearest)
      !! Overloaded Type bound procedure KdTreeSearch%kNearest()
    !====================================================================!
    class(KdTreeSearch), intent(inout) :: search
      !! KdTreeSearch class
    class(kdTree),intent(in) :: tree
      !! KdTree class
    real(r64),intent(in) :: x(:)
      !! x co-ordinates of the points
    real(r64),intent(in) :: y(:)
      !! y co-ordinates of the points
    real(r64),intent(in) :: z(:)
      !! z co-ordinates of the points
    real(r64),intent(in) :: xQuery
      !! x co-ordinate of the query location
    real(r64),intent(in) :: yQuery
      !! y co-ordinate of the query location
    real(r64),intent(in) :: zQuery
      !! z co-ordinate of the query location
    integer(i32), intent(in), optional :: k
      !! Number of points to find that are closest to the query
    real(r64), intent(in), optional :: radius
      !! Only find neighbours within this distance from the query
    type(dArgDynamicArray) :: kNearest
      !! Indices of the nearest points to the query location
    end function
    !====================================================================!
    !====================================================================!
    module function kNearest_KD(search, tree, D, query, k, radius) result(kNearest)
      !! Overloaded Type bound procedure KdTreeSearch%kNearest()
    !====================================================================!
    class(KdTreeSearch), intent(inout) :: search
      !! KdTreeSearch class
    class(kdTree),intent(in) :: tree
      !! KdTree class
    real(r64),intent(in) :: D(:,:)
      !! Co-ordinates of the points, the second dimension contains the k dimensions
    real(r64),intent(in) :: query(:)
      !! Co-ordinate of the query location
    integer(i32), intent(in), optional :: k
      !! Number of points to find that are closest to the query
    real(r64), intent(in), optional :: radius
      !! Only find neighbours within this distance from the query
    type(dArgDynamicArray) :: kNearest
      !! Indices of the nearest points to the query location
    end function
    !====================================================================!
  end interface

  interface
    !====================================================================!
    module function rangeSearch_2D(search, tree, x, y, lowerBound, upperBound) result(iPoints)
      !! Overloaded Type bound procedure KdTreeSearch%rangeSearch()
    !====================================================================!
    class(KdTreeSearch), intent(inout) :: search
      !! KdTreeSearch class
    class(kdTree),intent(in) :: tree
      !! KdTree class
    real(r64),intent(in) :: x(:)
      !! x co-ordinates of the points.
    real(r64),intent(in) :: y(:)
      !! y co-ordinates of the points.
    real(r64),intent(in) :: lowerBound(2)
      !! Lower bounds in x and y for the range search [xLow, yLow].
    real(r64),intent(in) :: upperBound(2)
      !! Upper bounds in x and y for the range search [xHigh, yHigh].
    type(iDynamicArray) :: iPoints
      !! Indices of the points inside the axis aligned range.
    end function
    !====================================================================!
    !====================================================================!
    module function rangeSearch_3D(search, tree, x, y, z, lowerBound, upperBound) result(iPoints)
      !! Overloaded Type bound procedure KdTreeSearch%rangeSearch()
    !====================================================================!
    class(KdTreeSearch), intent(inout) :: search
      !! KdTreeSearch class
    class(kdTree),intent(in) :: tree
      !! KdTree class
    real(r64),intent(in) :: x(:)
      !! x co-ordinates of the points.
    real(r64),intent(in) :: y(:)
      !! y co-ordinates of the points.
    real(r64),intent(in) :: z(:)
      !! y co-ordinates of the points.
    real(r64),intent(in) :: lowerBound(3)
      !! Lower bounds in x and y for the range search [xLow, yLow].
    real(r64),intent(in) :: upperBound(3)
      !! Upper bounds in x and y for the range search [xHigh, yHigh].
    type(iDynamicArray) :: iPoints
      !! Indices of the points inside the axis aligned range.
    end function
    !====================================================================!
    !====================================================================!
    module function rangeSearch_KD(search, tree, D, lowerBound, upperBound) result(iPoints)
      !! Overloaded Type bound procedure KdTreeSearch%rangeSearch()
    !====================================================================!
    class(KdTreeSearch), intent(inout) :: search
      !! KdTreeSearch class
    class(kdTree),intent(in) :: tree
      !! KdTree class
    real(r64),intent(in) :: D(:,:)
      !! Co-ordinates of the points, the second dimension contains the k dimensions
    real(r64),intent(in) :: lowerBound(size(D,2))
      !! Lower bounds in each dimension for the range search [Low1, Low2, ..., LowK].
    real(r64),intent(in) :: upperBound(size(D,2))
      !! Upper bounds in each dimension for the range search [High1, High2, ..., HighK].
    type(iDynamicArray) :: iPoints
      !! Indices of the points inside the axis aligned range.
    end function
    !====================================================================!
  end interface
end module

max_stars_repo_name: leonfoks/coretran

max_issues_repo_name: leonfoks/coretran

max_forks_repo_name: leonfoks/coretran


text:
! This program is used to convert surface files from SUmb into binary
! tecplot files. This script optionally strips off the symmetry and
! far field zones which are often not needed.

program cgns_surf_convert

  implicit none
  include 'cgnslib_f.h'

  integer, dimension(3*3) :: isize
  integer, dimension(6) :: rinde

  ! File handles for cgns files
  integer :: cg_current

  ! CGNS counters for base
  integer current_base

  ! CGNS Names for base, zone and iterative data
  character*32 basename
  character*1024 input_name
  character*1024 output_name
  character*1024 max_name

  character*32 fieldname
  character*36 zonename
  ! CGNS Data type
  integer datatype

  ! CGNS Names
  character*32 gridName,solName
  character*32, dimension(:),allocatable :: fieldNames
  character*32, dimension(:),allocatable :: var_names

  ! CGNS cell dimension and physical dimension
  integer  CellDim, PhysDim

  ! CGNS number of bases and zones in current file
  integer nbases, nzones

  ! CGNS error flag
  integer ier

  ! Integer Counters etc
  integer :: N,i,j,ii,jj,kk,tt
  integer :: nx,ny,nz,il,jl,kl
  integer :: nsols,sol,nfields,field
  integer :: temp_shape(6),rmax(3)
  integer :: farStart, symStart, cgnsStart

  ! Character arrays for names and I/O
  character*32, dimension(3) :: coordNames

  ! Misc Data Arrays
  double precision, allocatable,dimension(:,:,:  ) :: data3d
  double precision, allocatable,dimension(:,:,:,:) :: gridTemp
  real(kind=4), allocatable, dimension(:,:,:,:) :: gridTempSingle
  double precision, allocatable,dimension(:,:,:  ) :: fieldData
  real(kind=4), allocatable, dimension(:, :, :) :: data3dSingle
  double precision :: xmin(3), xmax(3)

  ! Tecplot Data Stuff
  ! Working
  ! These Never Change
  integer VIsDouble                /1/
  integer debug                    /0/
  integer FileType                 /0/
  INTEGER ZoneType                 /0/
  INTEGER TECINI112,tecdat112,teczne112,tecend112

  INTEGER ICellMax                 / 0/
  INTEGER JCellMax                 / 0/
  INTEGER KCellMax                 / 0/
  INTEGER DIsDouble                / 1/
  INTEGER StrandID                 / 0/      !/* StaticZone */
  INTEGER ParentZn                 / 0/
  INTEGER IsBlock                  / 1/      !/* Block */
  INTEGER NFConns                  / 0/
  INTEGER FNMode                   / 0/
  INTEGER TotalNumFaceNodes        / 0/
  INTEGER TotalNumBndryFaces       / 0/
  INTEGER TotalNumBndryConnections / 0/
  INTEGER ShrConn                  / 0/
  INTEGER,dimension(:),allocatable :: ValueLocation,passiveVarList,ShareVarFromZone

  ! Start of Code
  coordNames(1) = "CoordinateX"
  coordNames(2) = "CoordinateY"
  coordNames(3) = "CoordinateZ"
  xmin(:) = 0d0
  xmax(:) = 0d0
  N = IARGC ()

  if (.not. ( N == 1 )) then
     print *,'Error: cgns_surf_convert must be called with ONE argument:'
     print *,'./cgns_surf_convert <infile.cgns>'
     print *,'Output file will be called <infile.plt>'
     stop
  end if
  CALL GETARG(1 , input_name)

  cgnsStart = INDEX (trim(input_name),'.cgns')
  if (cgnsStart == 0) then
     print *,'Error: no .cgns extension found'
     stop
  end if
 10 format(A,A)
  write(output_name,10),input_name(1:cgnsStart),'plt'
  write(max_name,10),input_name(1:cgnsStart-1),'_max.txt'
  print *,'Converting file:',trim(input_name)
  ! ------------------------------------------------------------------
  ! Open a file to get the number of zones:

  call cg_open_f(input_name,MODE_READ,cg_current,ier)
  if (ier .eq. CG_ERROR) call cg_error_exit_f

  call cg_nbases_f(cg_current, nbases, ier)
  if (ier .eq. CG_ERROR) call cg_error_exit_f

  if (nbases .lt. 1) then
     print *, 'Error: No bases found in first CGNS file'
     stop
  end if

  ! We will only deal with 1 base in each zone so ALWAYS only use 1 base
  current_base = 1

  ! Get the cell and phys dim
  call cg_base_read_f(cg_current, current_base, basename, &
       CellDim, PhysDim, ier)
  if (ier .eq. CG_ERROR) call cg_error_exit_f

  call cg_nzones_f(cg_current, current_base, nzones, ier)
  if (ier .eq. CG_ERROR) call cg_error_exit_f

  call cg_nsols_f(cg_current, current_base, 1, nsols , ier )
  if (ier .eq. CG_ERROR) call cg_error_exit_f

  if (nsols >= 1) then
     call cg_nfields_f(cg_current,current_base,1,1,nfields,ier)
     if (ier .eq. CG_ERROR) call cg_error_exit_f
  else
     nfields = 0
  end if

  ! Allocate some names
  allocate(fieldNames(nfields))
  allocate(var_names(nFields+4))
  allocate(valueLocation(nFields+physDim),&
       passiveVarList(nFields+physDim),&
       ShareVarFromZone(nFields+physDim))

  ! Loop over to get field names
  do j=1,nFields
     call cg_field_info_f(cg_current,current_base,1,1,j,&
          datatype , fieldname , ier )
     fieldNames(j) = fieldName
  end do

  ! No Passive Var List or Share from Zone
  passiveVarList(:) = 0
  ShareVarFromZone(:) = 0

  ! All Data is at Nodes
  valueLocation(:) = 1

  ! Accumulate Names
  do j=1,physDim
     var_names(j) = coordNames(j)
  end do

  do j=1,nFields
     var_names(j+physDim) = fieldNames(j)
  end do
  var_names(nFields+3+1) = char(0)

  ! Open Ouput Tecplot File
  ier = TECINI112("CGNS Surf Data"//char(0),var_names,trim(output_name)//char(0),"."//char(0),&
       FileType,Debug,VIsDouble)

  ! ------------------------------------------------------------------
  ! Loop over the zones to get the required size info:
  do i=1,nzones

     call cg_zone_read_f(cg_current,current_base,i,&
          zoneName,isize,ier)
     if (ier .eq. CG_ERROR) call cg_error_exit_f

     farStart = INDEX (trim(zoneName),'Far')
     symStart = INDEX (trim(zoneName),'Sym')

     if (farStart == 0 .and. symStart == 0) then

        ! Load the zone sizes:
        call cg_zone_read_f(cg_current,current_base,i,zoneName,isize,ier)
        if (ier .eq. CG_ERROR) call cg_error_exit_f

        if (cellDim == 2) then
           nx = isize(1)
           ny = isize(2)
           nz = 1
           il = isize(3)
           jl = isize(4)
           kl = 1
        else
           nx = isize(1)
           ny = isize(2)
           nz = isize(3)
           il = isize(4)
           jl = isize(5)
           kl = isize(6)
        end if

        ! Write Tecplot Zone Data
        ier = TECZNE112(trim(zoneName),ZoneType,nx,ny,nz,&
             ICellMax,JCellMax,KCellMax, &
             0.0,0,ParentZn,IsBlock,NFConns,FNMode,&
             TotalNumFaceNodes,TotalNumBndryFaces,TotalNumBndryConnections,&
             PassiveVarList, valueLocation, ShareVarFromZone,ShrConn)

        !-------------------------------------------------------------
        !                         Grid Coodinates
        !-------------------------------------------------------------

        ! Find out the datatype
        call cg_coord_info_f(cg_current,current_base,i,&
             1,DataType,coordNames(1),ier)

        allocate(gridTemp(nx,ny,nz,physDim))

        if (dataType == RealSingle) then
           allocate(gridTempSingle(nx,ny,nz,physDim))
           do j=1,physDim
              call cg_coord_read_f(cg_current,current_base,i,&
                   coordNames(j),dataType,(/1,1,1/),(/nx,ny,nz/),&
                   gridTemp(:,:,:,j),ier)
              if (ier .eq. CG_ERROR) call cg_error_exit_f
           end do

           do ii=1,nx
              do jj=1,ny
                 do kk = 1,nz
                    do tt = 1,physDim
                       gridTemp(ii,jj,kk,tt) = gridTempSingle(ii,jj,kk,tt)
                    end do
                 end do
              end do
           end do

           deallocate(gridTempSingle)
        else

           do j=1,physDim
              call cg_coord_read_f(cg_current,current_base,i,&
                   coordNames(j),dataType,(/1,1,1/),(/nx,ny,nz/),&
                   gridTemp(:,:,:,j),ier)
              if (ier .eq. CG_ERROR) call cg_error_exit_f
           end do
        end if

        ! Get the min and max vals:
        if (i == 1) then
           do j=1,physDim
              xmin(j) = minval(gridTemp(:,:,:,j))
              xmax(j) = maxval(gridTemp(:,:,:,j))
           end do
        else
           do j=1,physDim
              xmin(j) = minval((/minval(gridTemp(:,:,:,j)),xmin(j)/))
              xmax(j) = maxval((/maxval(gridTemp(:,:,:,j)),xmax(j)/))
           end do
        end if

        ! Write Grid Coordinates
        do j=1,physDim
           ier   = TECDAT112(nx*ny*nz,gridTemp(:,:,:,j), DIsDouble)
        end do

        deallocate(gridTemp)

        !-------------------------------------------------------------
        !                     Field Values
        !-------------------------------------------------------------
        if (nfields > 0) then
           ! Read the Rind Data
           call cg_goto_f(cg_current,current_base,ier,'Zone_t',i, &
                'FlowSolution_t',1,'end')
           rinde(:) = 0
           call cg_rind_read_f(rinde , ier )
           if (ier .eq. CG_ERROR) call cg_error_exit_f

           ! Load the Variables
           do j=1,nFields
              ! Allocate node based field data
              allocate(fieldData(nx,ny,nz))

              call cg_field_info_f(cg_current,current_base,i,1,j,&
                   datatype , fieldname , ier )
              if (ier .eq. CG_ERROR) call cg_error_exit_f
              temp_shape(:) = 0
              temp_shape = (/1 - rinde(1), il + rinde(2), &
                   1 - rinde(3), jl + rinde(4), &
                   1 - rinde(5), kl + rinde(6)/)

              ! Allocate Temporary Array for Cell Data, possibly with rind cells
              allocate(data3d(temp_shape(1):temp_shape(2),&
                   temp_shape(3):temp_shape(4),&
                   temp_shape(5):temp_shape(6)))

              rmax(1) = il+rinde(1)+rinde(2)
              rmax(2) = jl+rinde(3)+rinde(4)
              rmax(3) = kl+rinde(5)+rinde(6)

              if (dataType == RealSingle) then
                 allocate(data3dSingle(temp_shape(1):temp_shape(2),&
                      temp_shape(3):temp_shape(4),&
                      temp_shape(5):temp_shape(6)))

                 call cg_field_read_f(cg_current,current_base,i,1, &
                      fieldName, datatype, (/1,1,1/),&
                      rmax,data3dSingle,ier)
                 if (ier .eq. CG_ERROR) call cg_error_exit_f

                 do kk=temp_shape(5),temp_shape(6)
                    do jj=temp_shape(3),temp_shape(4)
                       do ii=temp_shape(1),temp_shape(2)
                          data3d(ii,jj,kk) = data3dSingle(ii,jj,kk)
                       end do
                    end do
                 end do
                 deallocate(data3dSingle)
              else
                 call cg_field_read_f(cg_current,current_base,i,1, &
                      fieldName, datatype, (/1,1,1/),&
                      rmax,data3d,ier)
                 if (ier .eq. CG_ERROR) call cg_error_exit_f
              end if

              ! Call the interpolate function to reconstruct node data
              call interpolate(data3d,temp_shape,&
                   fieldData,nx,ny,nz)

              ! Deallocate Temporary Data Array
              deallocate(data3d)

              ! Write out Field Values
              ier   = TECDAT112(nx*ny*nz,fieldData, DIsDouble)
              deallocate(fieldData)
           end do ! Field Loop
        end if ! If we have fields
     end if
  end do ! Zone Loop

  call cg_close_f(cg_current,ier)
  if (ier .eq. CG_ERROR) call cg_error_exit_f

  ! Close Tecplot File
  ier = TECEND112()

  ! Open the max file and write the 6 values
  OPEN(UNIT=11, FILE=max_name, STATUS="replace")
  write(11, '(20F15.10)') xmin(1)
  write(11, '(20F15.10)') xmin(2)
  write(11, '(20F15.10)') xmin(3)
  write(11, '(20F15.10)') xmax(1)
  write(11, '(20F15.10)') xmax(2)
  write(11, '(20F15.10)') xmax(3)
  close(11)

end program cgns_surf_convert

subroutine interpolate(cell_data,cell_shape,node_data,inode,jnode,knode)

  ! Take cell centered data which possibly contain rind cells and
  ! perform a linear reconstruction back to the nodes. cell_shape
  ! contains the start and end values for each index in
  ! cell_data. Retruned data comes back in node_data. It may not be
  ! the most efficient possible, but it works.

  implicit none

  integer, intent(in) :: cell_shape(6),inode,jnode,knode

  double precision, intent(in ), dimension(cell_shape(1):cell_shape(2),&
       cell_shape(3):cell_shape(4),&
       cell_shape(5):cell_shape(6)) :: cell_data
  double precision, intent(out), dimension(inode,jnode,knode) :: node_data
  integer :: ic_s,ic_e,jc_s,jc_e,kc_s,kc_e
  integer, dimension(inode,jnode,knode) :: node_cnts
  integer :: i,j,k

  ic_s = cell_shape(1)
  ic_e = cell_shape(2)
  jc_s = cell_shape(3)
  jc_e = cell_shape(4)
  kc_s = cell_shape(5)
  kc_e = cell_shape(6)

  node_data(:,:,:) = 0.0
  node_cnts(:,:,:) = 0

  do k=kc_s,kc_e
     do j=jc_s,jc_e
        do i=ic_s,ic_e

           if (i >= 1 .and. i <= inode .and. j >= 1 .and. j <= jnode .and. k >=1 .and. k <= knode) then
              node_data(i  ,j  ,k  ) = node_data(i  ,j  ,k  ) + cell_data(i,j,k)
              node_cnts(i  ,j  ,k  ) = node_cnts(i  ,j  ,k  ) + 1
           end if

           if (i+1 >= 1 .and. i+1 <= inode .and. j >= 1 .and. j <= jnode .and. k >=1 .and. k <= knode) then
              node_data(i+1,j  ,k  ) = node_data(i+1,j  ,k  ) + cell_data(i,j,k)
              node_cnts(i+1,j  ,k  ) = node_cnts(i+1,j  ,k  ) + 1
           end if

           if (i >= 1 .and. i <= inode .and. j+1 >= 1 .and. j+1 <= jnode .and. k >=1 .and. k <= knode) then
              node_data(i  ,j+1,k  ) = node_data(i  ,j+1,k  ) + cell_data(i,j,k)
              node_cnts(i  ,j+1,k  ) = node_cnts(i  ,j+1,k  ) + 1
           end if

           if (i+1 >= 1 .and. i+1 <= inode .and. j+1 >= 1 .and. j+1 <= jnode .and. k >=1 .and. k <= knode) then
              node_data(i+1,j+1,k  ) = node_data(i+1,j+1,k  ) + cell_data(i,j,k)
              node_cnts(i+1,j+1,k  ) = node_cnts(i+1,j+1,k  ) + 1
           end if

           if (i >= 1 .and. i <= inode .and. j >= 1 .and. j <= jnode .and. k+1 >=1 .and. k+1 <= knode) then
              node_data(i  ,j  ,k+1) = node_data(i  ,j  ,k+1) + cell_data(i,j,k)
              node_cnts(i  ,j  ,k+1) = node_cnts(i  ,j  ,k+1) + 1
           end if

           if (i+1 >= 1 .and. i+1 <= inode .and. j >= 1 .and. j <= jnode .and. k+1>=1 .and. k+1 <= knode) then
              node_data(i+1,j  ,k+1) = node_data(i+1,j  ,k+1) + cell_data(i,j,k)
              node_cnts(i+1,j  ,k+1) = node_cnts(i+1,j  ,k+1) + 1
           end if

           if (i >= 1 .and. i <= inode .and. j+1 >= 1 .and. j+1 <= jnode .and. k+1>=1 .and. k+1<= knode) then
              node_data(i  ,j+1,k+1) = node_data(i  ,j+1,k+1) + cell_data(i,j,k)
              node_cnts(i  ,j+1,k+1) = node_cnts(i  ,j+1,k+1) + 1
           end if

           if (i+1 >= 1 .and. i+1 <= inode .and. j+1 >= 1 .and. j+1 <= jnode .and. k+1 >=1 .and. k+1 <= knode) then
              node_data(i+1,j+1,k+1) = node_data(i+1,j+1,k+1) + cell_data(i,j,k)
              node_cnts(i+1,j+1,k+1) = node_cnts(i+1,j+1,k+1) + 1
           end if
        end do
     end do
  end do

  ! Divide through by the number of times a value was added to a node

  do k=1,knode
     do j=1,jnode
        do i=1,inode
           node_data(i,j,k) = node_data(i,j,k) / dble(node_cnts(i,j,k))
        end do
     end do
  end do

end subroutine interpolate

max_stars_repo_name: GuardZ/cgnsutilities

max_issues_repo_name: GuardZ/cgnsutilities

max_forks_repo_name: GuardZ/cgnsutilities


text:
This editor can edit this entry and tell us a bit about themselves by clicking the Edit icon.

20080906 17:42:37 nbsp Welcome to the Wiki Howdy, Ran, and Welcome to the Wiki!   You might want to check out the importance of using your RealName, just so we can get to know you (or not: its your choice, but people were pretty friendly here).  My names Evan, pleased to meet you!  Thanks for the reviews  you sure seem cheerful and know what you like! Users/JabberWokky Evan JabberWokky Edwards

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      DOUBLE PRECISION FUNCTION DT_RNDM(Vdummy)
C***********************************************************************
C                                                                      *
C                 7) Random number generator package                   *
C                                                                      *
C    THIS IS A PACKAGE CONTAINING A RANDOM NUMBER GENERATOR AND        *
C    SERVICE ROUTINES.                                                 *
C    THE ALGORITHM IS FROM                                             *
C      'TOWARD A UNVERSAL RANDOM NUMBER GENERATOR'                     *
C      G.MARSAGLIA, A.ZAMAN ;  FSU-SCRI-87-50                          *
C    IMPLEMENTATION BY K. HAHN  DEC. 88,                               *
C    THIS GENERATOR SHOULD NOT DEPEND ON THE HARD WARE ( IF A REAL HAS *
C    AT LEAST 24 SIGNIFICANT BITS IN INTERNAL REPRESENTATION ),        *
C    THE PERIOD IS ABOUT 2**144,                                       *
C    TIME FOR ONE CALL AT IBM-XT IS ABOUT 0.7 MILLISECONDS,            *
C    THE PACKAGE CONTAINS                                              *
C      FUNCTION DT_RNDM(I)                  : GENERATOR                *
C      SUBROUTINE DT_RNDMST(NA1,NA2,NA3,NB4): INITIALIZATION           *
C      SUBROUTINE DT_RNDMIN(U,C,CD,CM,I,J)  : PUT SEED TO GENERATOR    *
C      SUBROUTINE DT_RNDMOU(U,C,CD,CM,I,J)  : TAKE SEED FROM GENERATOR *
C      SUBROUTINE DT_RNDMTE(IO)             : TEST OF GENERATOR        *
C---                                                                   *
C    FUNCTION DT_RNDM(I)                                               *
C       GIVES UNIFORMLY DISTRIBUTED RANDOM NUMBERS  IN (0..1)          *
C       I  - DUMMY VARIABLE, NOT USED                                  *
C    SUBROUTINE DT_RNDMST(NA1,NA2,NA3,NB1)                             *
C       INITIALIZES THE GENERATOR, MUST BE CALLED BEFORE USING DT_RNDM *
C       NA1,NA2,NA3,NB1  - VALUES FOR INITIALIZING THE GENERATOR       *
C                          NA? MUST BE IN 1..178 AND NOT ALL 1         *
C                          12,34,56  ARE THE STANDARD VALUES           *
C                          NB1 MUST BE IN 1..168                       *
C                          78  IS THE STANDARD VALUE                   *
C    SUBROUTINE DT_RNDMIN(U,C,CD,CM,I,J)                               *
C       PUTS SEED TO GENERATOR ( BRINGS GENERATOR IN THE SAME STATUS   *
C       AS AFTER THE LAST DT_RNDMOU CALL )                             *
C       U(97),C,CD,CM,I,J  - SEED VALUES AS TAKEN FROM DT_RNDMOU       *
C    SUBROUTINE DT_RNDMOU(U,C,CD,CM,I,J)                               *
C       TAKES SEED FROM GENERATOR                                      *
C       U(97),C,CD,CM,I,J  - SEED VALUES                               *
C    SUBROUTINE DT_RNDMTE(IO)                                          *
C       TEST OF THE GENERATOR                                          *
C       IO     - DEFINES OUTPUT                                        *
C                  = 0  OUTPUT ONLY IF AN ERROR IS DETECTED            *
C                  = 1  OUTPUT INDEPENDEND ON AN ERROR                 *
C       DT_RNDMTE USES DT_RNDMIN AND DT_RNDMOU TO BRING GENERATOR TO   *
C       SAME STATUS                                                    *
C       AS BEFORE CALL OF DT_RNDMTE                                    *
C***********************************************************************

      IMPLICIT NONE
      DOUBLE PRECISION C , CD , CM , FLRNDM , U , Vdummy
      INTEGER I , J
      SAVE

C random number generator
#ifndef FOR_FLUKA
      COMMON /DTRAND/ U(97) , C , CD , CM , I , J
#endif
C counter of calls to random number generator
C uncomment if needed
C     COMMON /DTRNCT/ IRNCT0,IRNCT1
C     LOGICAL LFIRST
C     DATA LFIRST /.TRUE./

C counter of calls to random number generator
C uncomment if needed
C     IF (LFIRST) THEN
C        IRNCT0 = 0
C        IRNCT1 = 0
C        LFIRST = .FALSE.
C     ENDIF
#ifdef FOR_FLUKA
      DT_RNDM = FLRNDM(Vdummy)
#else
 100  DT_RNDM = U(I) - U(J)
      IF ( DT_RNDM.LT.0.0D0 ) DT_RNDM = DT_RNDM + 1.0D0
      U(I) = DT_RNDM
      I = I - 1
      IF ( I.EQ.0 ) I = 97
      J = J - 1
      IF ( J.EQ.0 ) J = 97
      C = C - CD
      IF ( C.LT.0.0D0 ) C = C + CM
      DT_RNDM = DT_RNDM - C
      IF ( DT_RNDM.LT.0.0D0 ) DT_RNDM = DT_RNDM + 1.0D0

      IF ( (DT_RNDM.EQ.0.D0) .OR. (DT_RNDM.EQ.1.D0) ) GOTO 100
#endif
C counter of calls to random number generator
C uncomment if needed
C     IRNCT0 = IRNCT0+1

      END FUNCTION

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


text:
!    Generate sigma coordinates for ivcor=1 (VQS from Dukhovskoy) using
!    multiple masters.
!    WARNING: most variables/arrays in this program use '1' as surface, and nvrt*
!    as bottom!!!

!    Inputs: (1) consts. inside the code (for max. flexibility); (2) hgrid.gr3 (in map projection or lon/lat);
!            (3) screen; (4) transect.bp (optional transect bp file; depths can be seg #s);
!            (5) v_regions.gr3 (depths define which master grid is used
!            locally. In this example, '2' is Black Sea, '1' is rest,
!            and the two are smoothly stitched in middle of Bosphorus St where
!            the depth <=100m, because the first 3 master grids are identical
!            for h<=100m)
!    Outputs: vgrid.in; vgrid_master*.out;  transect1.out; debug outputs (fort*)
!    Use plot_VQS.m to viz vgrid_master*.out, transect1.out
!    ifort -O2 -mcmodel=medium -CB -Bstatic -o gen_vqs_2masters.exe ../UtilLib/schism_geometry.f90 gen_vqs_2masters.f90

      use schism_geometry_mod
      implicit real*8(a-h,o-z)
      integer, allocatable :: elnode(:,:),elside(:,:),kbp0(:),kbp(:),ic3(:,:),isdel(:,:),isidenode(:,:),m0(:)
      allocatable :: xnd(:),ynd(:),dp(:),xcj(:),ycj(:),eta2(:),znd(:,:),z1tmp(:),z2tmp(:)
      allocatable :: hsm(:,:),nv_vqs(:,:),z_mas(:,:,:) !,a_vqs(:)
      allocatable :: hsm0(:),nv_vqs0(:),z_mas0(:,:)
      allocatable :: xybp(:,:),dpbp(:),imap(:),transect_len(:),sigma_vqs(:,:),i34(:),iflag(:),ireg(:)

      print*, 'Want to output along a transect? (0: no; 1:yes)'
      read*, itran

      !# of master sets used in different regions (share same m_vqs)
      max_m=2

      !m_vqs: max # of grids/depths in each master grid
      m_vqs=10
      dz_bot_min=3 !min. bottom layer thickness [m]
      allocate(hsm(max_m,m_vqs),nv_vqs(max_m,m_vqs))

      !2nd master set: Black Sea
      hsm(2,:)=(/20.,60.,100.,140.,180.,220.,260.,2400.,3000.,4000./)
      nv_vqs(2,1:7)=(/15,19,23,26,29,31,33/) !# of levels for each master grid
      !grid #9, 10 are not used (exceed max depth of Black Sea)
      nv_vqs(2,8:m_vqs)=nv_vqs(2,7)+15

      !1st: rest
      hsm(1,:)=(/20.,60.,100.,140.,300.,500.,1.e3,2.e3,3.e3,6.e3/)
      nv_vqs(1,1:7)=(/15,19,23,26,31,33,39/)
      nv_vqs(1,8)=nv_vqs(1,7)+6
      nv_vqs(1,9)=nv_vqs(1,8)+5
      nv_vqs(1,10)=nv_vqs(1,9)+8

      if(m_vqs<2) then
        write(*,*)'Check vgrid.in:',m_vqs
        stop
      endif
      if(hsm(1,1)<0) stop 'hsm(1,1)<0'
      do m=2,m_vqs
        do i=1,max_m
          if(hsm(i,m)<=hsm(i,m-1)) then
            write(*,*)'Check hsm:',m,i,hsm(i,m),hsm(i,m-1)
            stop
          endif
        enddo !i
      enddo !m

!     Other consts.
!     Stretching const. for the 1st master grid (optional)
!     |a_vqs0|<=1 (1: skew toward bottom; -1: toward surface; 0: no bias)
      a_vqs0=-1

!     Generate a master vgrid (z_mas)
      etal=0 !used in master grid only; elev.
      if(etal<=-hsm(1,1)) then
        write(*,*)'elev<hsm:',etal
        stop
      endif

      nvrt_m=maxval(nv_vqs(:,m_vqs))
      print*, 'nvrt in master vgrid=',nvrt_m
      allocate(z_mas(max_m,nvrt_m,m_vqs))
      z_mas=-1.e5
      do m=1,7 !m_vqs
        do j=1,max_m
          do k=1,nv_vqs(j,m)
            sigma=(k-1.0)/(1.0-nv_vqs(j,m))

!           Alternative transformations below
!           Option 1: quadratic
!           a_vqs(m)=max(-1.d0,a_vqs0-(m-1)*0.03)
!           tmp=a_vqs(m)*sigma*sigma+(1+a_vqs(m))*sigma !transformed sigma
!           z_mas0(k,m)=tmp*(etal+hsm0(m))+etal

!           Option 2: S
            theta_b=0
            if(m<4) then
              theta_f=2.6
            else if(m==4) then
              theta_f=2.5
            else if(m==5) then
              theta_f=3.2
            else if(m==6) then
              theta_f=3.7
            else
              theta_f=4.3
            endif

            if(j==2) theta_f=2.6 !overwrite (make sure the first 3 grids match)

            cs=(1-theta_b)*sinh(theta_f*sigma)/sinh(theta_f)+ &
     &theta_b*(tanh(theta_f*(sigma+0.5))-tanh(theta_f*0.5))/2/tanh(theta_f*0.5)
            z_mas(j,k,m)=etal*(1+sigma)+hsm(j,1)*sigma+(hsm(j,m)-hsm(j,1))*cs
          enddo !k
        enddo !j
      enddo !m

      !Define last few master grids
      !1st set
      z_mas(1,1:nv_vqs(1,7),8)=z_mas(1,1:nv_vqs(1,7),7)
      z_mas(1,1+nv_vqs(1,7):nv_vqs(1,8),8)=(/-1100.,-1250.,-1400.,-1600.,-1800.,-2000./)
      z_mas(1,1:nv_vqs(1,8),9)=z_mas(1,1:nv_vqs(1,8),8)
      z_mas(1,1+nv_vqs(1,8):nv_vqs(1,9),9)=(/-2200.,-2400.,-2600.,-2800.,-3000./)
      z_mas(1,1:nv_vqs(1,9),10)=z_mas(1,1:nv_vqs(1,9),9)
      z_mas(1,1+nv_vqs(1,9):nv_vqs(1,10),10)=(/-3300.,-3600.,-4000.,-4400.,-4800.,-5200.,-5600.,-6000./)

      !2nd set
      z_mas(2,1:nv_vqs(2,7),8)=z_mas(2,1:nv_vqs(2,7),7)
      z_mas(2,nv_vqs(2,7)+1:nv_vqs(2,8),8)=-(/280,300,320,350, &
     &400,470,560,670,800,1000,1200,1400,1700,2000,2400/)
      !last 2 master grids are not used (exceed max depth)
      z_mas(2,1:nv_vqs(2,9)-1,9)=z_mas(2,1:nv_vqs(2,9)-1,8)
      z_mas(2,nv_vqs(2,9),9)=-3.e3
      z_mas(2,1:nv_vqs(2,10)-1,10)=z_mas(2,1:nv_vqs(2,10)-1,8)
      z_mas(2,nv_vqs(2,10),10)=-4.e3

!     Output master sets
      open(31,file='vgrid_master1.out',status='replace')
      open(32,file='vgrid_master2.out',status='replace')
      do m=1,m_vqs
        write(31,'(2(1x,i5),6000(1x,f12.4))')m,nv_vqs(1,m),hsm(1,m),z_mas(1,:,m)
        write(32,'(2(1x,i5),6000(1x,f12.4))')m,nv_vqs(2,m),hsm(2,m),z_mas(2,:,m)
      enddo !m
      close(31)

!      do k=1,nvrt_m
!        write(12,'(i5,6000(1x,f12.4))')k,z_mas(1,k,:)
!      enddo !k

      nvrt=nvrt_m
!     Read in hgrid
      open(14,file='hgrid.gr3',status='old')
      open(17,file='v_regions.gr3',status='old')
      read(14,*)
      read(14,*)ne,np
      read(17,*); read(17,*)
      allocate(xnd(np),ynd(np),dp(np),i34(ne),elnode(4,ne),kbp0(np),kbp(np),eta2(np), &
     &sigma_vqs(nvrt,np),znd(nvrt,np),m0(np),z1tmp(nvrt),z2tmp(nvrt),ireg(np),iflag(np))
      eta2=etal
      do i=1,np
        read(14,*)j,xnd(i),ynd(i),dp(i)
        read(17,*)j,tmp,tmp,tmp2
        ireg(i)=nint(tmp2)
      enddo !i
      do i=1,ne
        read(14,*)j,i34(i),elnode(1:i34(i),i)
      enddo !i
      close(14)

      !Check max depth
      dpmax=maxval(dp)
      if(dpmax>maxval(hsm(:,m_vqs))) then
        print*, 'Max depth exceeds master depth:',dpmax,maxval(hsm(:,m_vqs))
        stop
      endif

      call compute_nside(np,ne,i34,elnode,ns)
      print*, 'ns=',ns
      allocate(ic3(4,ne),elside(4,ne),isdel(2,ns),isidenode(2,ns),xcj(ns),ycj(ns))
      call schism_geometry_double(np,ne,ns,xnd,ynd,i34,elnode,ic3,elside,isdel,isidenode,xcj,ycj)
      !deallocate()

!     Compute zcoor
      znd=-1.e6
      iflag=0
      do i=1,np
        !Find master set
        if(ireg(i)==1) then !rest
          indx=1 !master set index
        else !Black Sea
          indx=2
        endif !ireg(i)

        if(dp(i)<=hsm(indx,1)) then !shallow; compute from 1st master
          if(etal+hsm(indx,1)<=0) then
            print*, 'Check etal+hsm(1):',etal+hsm(indx,1)
            stop
          endif

          iflag(i)=1
          kbp(i)=nv_vqs(indx,1)
          do k=1,nv_vqs(indx,1)
            !sigma=(k-1.0)/(1.0-nv_vqs(indx,1))
            sigma_vqs(k,i)=(z_mas(indx,k,1)-etal)/(etal+hsm(indx,1))
            znd(k,i)=sigma_vqs(k,i)*(eta2(i)+dp(i))+eta2(i)
          enddo !k
          cycle
        endif

        !Deep
        !Find a master vgrid
        m0(i)=0
        do m=2,m_vqs
          if(dp(i)>hsm(indx,m-1).and.dp(i)<=hsm(indx,m)) then
            m0(i)=m
            zrat=(dp(i)-hsm(indx,m-1))/(hsm(indx,m)-hsm(indx,m-1)) !(0,1]
            exit
          endif
        enddo !m
        if(m0(i)==0) then
          print*, 'Failed to find a master vgrid:',i,dp(i)
          stop
        endif

        kbp(i)=0
        do k=1,nv_vqs(indx,m0(i))
          z1=z_mas(indx,min(k,nv_vqs(indx,m0(i)-1)),m0(i)-1)
          z2=z_mas(indx,k,m0(i))
          z3=z1+(z2-z1)*zrat
          if(z3>=-dp(i)+dz_bot_min) then
            znd(k,i)=z3
          else
            kbp(i)=k; exit
          endif
        enddo !k
        if(kbp(i)==0) then
          print*, 'Failed to find a bottom:',i,dp(i),z3,z_mas(indx,1:nv_vqs(indx,m0(i)),m0(i))
          stop
        endif
        znd(kbp(i),i)=-dp(i)

        !Check order
        do k=2,kbp(i)
          if(znd(k-1,i)<=znd(k,i)) then
            print*, 'Inverted z:',i,dp(i),m0(i),k,znd(k-1,i),znd(k,i)
            stop
          endif
        enddo !k
        write(99,*)'Node:',i,real(dp(i)),real(znd(1:kbp(i),i))
      enddo !i=1,np

!     Extend beyond bottom for plotting
      do i=1,np
        znd(kbp(i)+1:nvrt,i)=-dp(i)
      enddo !i

!     Optional output along transect
      if(itran/=0) then
        open(9,file='transect.bp',status='old')
        read(9,*)
        read(9,*)npbp
        allocate(xybp(2,npbp),dpbp(npbp),imap(npbp),transect_len(npbp))
        transect_len(1)=0 !along transect distance
        do i=1,npbp
          read(9,*)j,xybp(1:2,i),dpbp(i)
          if(i>1) then
            rl=sqrt((xybp(1,i)-xybp(1,i-1))**2+(xybp(2,i)-xybp(2,i-1))**2)
            transect_len(i)=transect_len(i-1)+rl
          endif
        enddo !i
        close(9)
        !Find nearest node
        do i=1,npbp
          rlmin=huge(1.0d0)
          do j=1,np
            rl2=(xybp(1,i)-xnd(j))**2+(xybp(2,i)-ynd(j))**2
            if(rl2<rlmin) then
              imap(i)=j; rlmin=rl2
            endif
          enddo !j
        enddo !i

        open(13,file='transect1.out',status='replace')
        do i=1,npbp
          nd=imap(i)
          write(13,'(i6,1x,i4,2(1x,e16.6),10000(1x,f12.3))')i,kbp(nd),xybp(1:2,i),transect_len(i),dp(nd),dpbp(i),znd(:,nd)
        enddo !i
        close(13)
      endif !itran

      if(minval(kbp)<2) then
        print*, '# of levels <2:',minval(kbp)
        print*, 'Check parameters like dz_bot_min etc'
        stop
      endif

      nvrt=maxval(kbp)
      print*, 'Final nvrt=',nvrt
!     # of prisms
      nprism=0
      do i=1,ne
        kbpl=maxval(kbp(elnode(1:i34(i),i)))
        nprism=nprism+kbpl
      enddo !i

      print*, '# of prisms=',nprism
      print*, 'Average # of layers=',real(nprism)/ne

!     Output in SCHISM convention
      open(19,file='vgrid.in',status='replace')
      write(19,*)1 !ivcor
      write(19,*)nvrt
      do i=1,np
!        if(dp(i)<=hsm(1)) then
        if(iflag(i)/=0) then
          !sigma_vqs already assigned
        else
          sigma_vqs(1,i)=0
          sigma_vqs(kbp(i),i)=-1
          do k=2,kbp(i)-1
            sigma_vqs(k,i)=(znd(k,i)-eta2(i))/(eta2(i)+dp(i))
          enddo !k
        endif

        !Check order
        do k=2,kbp(i)
          if(sigma_vqs(k,i)>=sigma_vqs(k-1,i)) then
            print*, 'Inverted sigma:',i,k,dp(i),sigma_vqs(k,i),sigma_vqs(k-1,i)
            stop
          endif
        enddo !k

        write(19,'(2(1x,i10),10000(1x,f14.6))')i,nvrt+1-kbp(i),sigma_vqs(kbp(i):1:-1,i)
      enddo !i
      close(19)

!     Output horizontal map of # of levels for more adjustment
      open(13,file='nlev.gr3',status='replace')
      write(13,*)'# of levels at each node'
      write(13,*)ne,np
      do i=1,np
        write(13,*)i,xnd(i),ynd(i),kbp(i)
      enddo !i
      do i=1,ne
        write(13,*)i,i34(i),elnode(1:i34(i),i)
      enddo !i
      close(13)

      stop
      end

max_stars_repo_name: cseaton/schism

max_issues_repo_name: cseaton/schism

max_forks_repo_name: cseaton/schism


text:
	SUBROUTINE TI_CDTM  ( idate, itime, dattim, iret )
C************************************************************************
C* TI_CDTM								*
C*									*
C* This subroutine converts an integer date (YYMMDD) and time (HHMM) 	*
C* into a standard GEMPAK time.						*
C*									*
C* TI_CDTM  ( IDATE, ITIME, DATTIM, IRET )				*
C*									*
C* Input parameters:							*
C*	IDATE		INTEGER	 	Date (YYMMDD)			*
C*	ITIME		INTEGER	 	Time (HHMM)			*
C*									*
C* Output parameters:							*
C*	DATTIM		CHAR*	 	GEMPAK time			*
C*	IRET		INTEGER	 	Return code			*
C*					  0 = normal return		*
C*					 -1 = invalid time		*
C*					 -7 = invalid year		*
C*					 -8 = invalid month		*
C*					 -9 = invalid day		*
C*					-10 = invalid hour		*
C*					-11 = invalid minute		*
C**									*
C* Log:									*
C* I. Graffman/RDS	6/86						*
C* M. desJardins/GSFC	6/86	Added error check			*
C************************************************************************
	CHARACTER*(*)	dattim
	INTEGER		idtarr (5)
C------------------------------------------------------------------------
	idtarr (1) = idate / 10000
	idtarr (2) = ( idate - idtarr (1) * 10000 ) / 100
	idtarr (3) = MOD  ( idate, 100 )
	idtarr (4) = itime / 100
	idtarr (5) = MOD  ( itime, 100 )
	CALL TI_ITOC  ( idtarr, dattim, iret )
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: sgdecker/pygempak


        !COMPILER-GENERATED INTERFACE MODULE: Fri Mar  6 15:34:59 2020
        ! This source file is for reference only and may not completely
        ! represent the generated interface used by the compiler.
        MODULE ORB_PRINT__genmod
          INTERFACE
            SUBROUTINE ORB_PRINT(IYEAR_AD,ECCEN,OBLIQ,MVELP,MYPID,NROOT)
              INTEGER(KIND=4) :: IYEAR_AD
              REAL(KIND=8) :: ECCEN
              REAL(KIND=8) :: OBLIQ
              REAL(KIND=8) :: MVELP
              INTEGER(KIND=4) :: MYPID
              INTEGER(KIND=4) :: NROOT
            END SUBROUTINE ORB_PRINT
          END INTERFACE
        END MODULE ORB_PRINT__genmod


text:
!*******************************************************************************
!  File ipsl_mc.f
!  Contains module ipsl_mc

!  Module is part of the V3FIT stellarator equilibrium reconstruction code.
!  The postfix _mc in the name indicates Model Compute
!  The module contains the subroutine ipsl_mc_model_compute
!  Both an ipsl_desc and a model are needed for this computation.

!*******************************************************************************
!  MODULE ipsl_mc
!    (MDDC - Model Compute Calculations)
! SECTION I.     VARIABLE DECLARATIONS
! SECTION II.    INTERFACE BLOCKS
! SECTION III.   MAIN COMPUTATIONAL SUBROUTINE
!*******************************************************************************
      MODULE ipsl_mc

!*******************************************************************************
! SECTION I. VARIABLE DECLARATIONS
!*******************************************************************************

!-------------------------------------------------------------------------------
!  Type declarations, constants, utilities
!-------------------------------------------------------------------------------
      USE stel_kinds
      USE stel_constants
      USE v3_utilities

!-------------------------------------------------------------------------------
!  MDDC Derived Types
!-------------------------------------------------------------------------------
      USE ipsl_T

!-------------------------------------------------------------------------------
!  Model Derived Types
!-------------------------------------------------------------------------------
      USE model_T

!-------------------------------------------------------------------------------
!  IPSL specific modules
!-------------------------------------------------------------------------------
!      USE ip_beamline
      USE ipsl_T
      USE ip_global
      USE faraday_mod, ONLY: ipsl_integral

      IMPLICIT NONE

!*******************************************************************************
! SECTION II.  INTERFACE BLOCKS
!*******************************************************************************

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      CONTAINS

!*******************************************************************************
! SECTION III.  MAIN COMPUTATIONAL SUBROUTINE
!*******************************************************************************
!-------------------------------------------------------------------------------
!  Compute an MDDC signal
!
!    Information comes from the ipsl_desc and the model
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
!  Actual computation of the model signal is in this subroutine
!    s_type = diagnostic
!      d_type = ipsl
!        Magnetic Diagnostic-Dot Coil
!        signal_model_compute_ipsl
!           % data(1) - total flux
!           % data(2) - flux from the plasma currents
!           % data(3) - flux from all the external coils
!           % data(4) - data(3 + nextcur) - flux from individual external coils
!
!-------------------------------------------------------------------------------

      SUBROUTINE ipsl_mc_model_compute(a_ipsl,a_model,arr_data,                &
     &   arr_sigma)
!-------------------------------------------------------------------------
!
! FUNCTION: computes signals for interferometry/polarimetry straight
!           line" signals  (s_type = ipsl)
!
!           % data(1) -  Faraday Rotation angle      (ipsl_type = 'polarim')  OR
!                     -  Interferometry phase angle  (ipsl_type = 'interf')
!
!! created 6/19/07 by J. Shields
!-------------------------------------------------------------------------

!..........passed Arguments................................!
      TYPE (ipsl_desc), INTENT (inout), TARGET :: a_ipsl
      TYPE (model), INTENT (inout), TARGET :: a_model    ! equilbrium info
      REAL(rprec), POINTER, DIMENSION(:) :: arr_data, arr_sigma
!  Note - convention is that arr_data and arr_sigma
!    1) are deallocated in this subroutine, and then
!    2) are allocated in this subroutine

!.........Local Variables...............................!
      TYPE (eq_state), POINTER :: state => null()
      LOGICAL :: lfirst
      INTEGER :: n_data
      INTEGER :: i, k, mc, index
      REAL(rprec) :: line_integral_out
      INTEGER :: istat1, istat2
      CHARACTER(len=*), PARAMETER :: sub_name =                                &
     &  'ipsl_mc_model_compute: '

!      write(*,*) 'Now entering ', sub_name

!  Allocate data and sigma
      n_data = 1
      IF (ASSOCIATED(arr_data)) THEN
         DEALLOCATE(arr_data,arr_sigma, stat=istat1)
         CALL assert_eq(0,istat1,sub_name // 'arr_data, sigma dealloc')
      ENDIF
      ALLOCATE(arr_data(n_data),arr_sigma(n_data),stat=istat1)
      CALL assert_eq(0,istat1,sub_name // 'arr_data, sigma alloc')

      call IPSL_INTEGRAL(a_ipsl, line_integral_out)
      arr_data(1) = line_integral_out
      arr_sigma(1) = zero

      write(*,*) sub_name, 'output line integral = ', arr_data(1)

      RETURN
      END SUBROUTINE ipsl_mc_model_compute

!===================================================================================
      SUBROUTINE ipsl_model_init()
!  2009-03-03 JDH. Could this subroutine be eliminated?

!--------------------------------------------------------------------------------
!
! Function:  Tests VMEC/V3FIT interfacing routines for B transformations
!
!
!
!
!
!
!----------------------------------------------------------------------------------
      USE stel_kinds
      USE stel_constants
!      USE v3f_global   ! module containing V3FIT global data
      IMPLICIT NONE

!!..............dummy variables.............................................!
!      INTEGER, INTENT(in) :: i_beam     ! index of interferometer beam being integrated over

!............local variables.................................................!
      INTEGER      :: i, j, kp_junk
      REAL(rprec)               ::   blah
      CHARACTER(len=*), PARAMETER :: subname = 'V3F_GLOBAL_TEST: '

!      write(*,*) '===================================================='
!      write(*,*) ' '
      write(*,*) 'Hi.  Now in ', subname

!     call VMEC_B_INIT

      RETURN
      END SUBROUTINE ipsl_model_init

!*******************************************************************************
! SECTION X.  COMMENTS FOR DIFFERENT REVISIONS
!*******************************************************************************
!
!  JDH 2009-03-03. First version of module. Based on coding from previous
!    module signal_model
!
      END MODULE ipsl_mc


text:
!===============================================================================
!> $Id: 5d1a00ac7f333dc239ec9c7ca3826d6f0b4a7604 $
!===============================================================================
MODULE stagger_mod
  USE io_mod
  USE trace_mod
  USE stagger_16
  USE stagger_20
  USE stagger_24
  USE stagger_32
  USE stagger_36
  implicit none
  public
  integer, private:: verbose=0
  logical, save:: hardwire=.false.
  type, public:: stagger_t
    real(8):: ds(3)
    integer(kind=8):: flops=0
    integer(kind=8):: count=0
  contains
    procedure, nopass:: ddxdn1
    procedure, nopass:: ddydn1
    procedure, nopass:: ddzdn1
    procedure, nopass:: ddxup1
    procedure, nopass:: ddyup1
    procedure, nopass:: ddzup1
    procedure, nopass:: xdn1
    procedure, nopass:: ydn1
    procedure, nopass:: zdn1
    procedure, nopass:: xup1
    procedure, nopass:: yup1
    procedure, nopass:: zup1
    procedure, nopass:: ddx
    procedure, nopass:: ddy
    procedure, nopass:: ddz
    procedure, nopass:: ddxdn
    procedure, nopass:: ddydn
    procedure, nopass:: ddzdn
    procedure, nopass:: ddxup
    procedure, nopass:: ddyup
    procedure, nopass:: ddzup
    procedure, nopass:: xdn
    procedure, nopass:: ydn
    procedure, nopass:: zdn
    procedure, nopass:: xup
    procedure, nopass:: yup
    procedure, nopass:: zup
    procedure, nopass:: xsm
    procedure, nopass:: ysm
    procedure, nopass:: zsm
    procedure, nopass:: xyzsm
    procedure, nopass:: stagger_test
  end type
  type(stagger_t), public:: stagger
CONTAINS

!===============================================================================
FUNCTION ddxdn1 (ds, a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
  real(8):: ds(3)
    b= ddxdn (ds, a)
END FUNCTION ddxdn1

!*******************************************************************************
FUNCTION ddydn1 (ds, a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
  real(8):: ds(3)
    b= ddydn (ds, a)
END FUNCTION ddydn1

!*******************************************************************************
FUNCTION ddzdn1 (ds, a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
  real(8):: ds(3)
    b= ddzdn (ds, a)
END FUNCTION ddzdn1

!*******************************************************************************
FUNCTION ddxup1 (ds, a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
  real(8):: ds(3)
    b= ddxup (ds, a)
END FUNCTION ddxup1

!*******************************************************************************
FUNCTION ddyup1 (ds, a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
  real(8):: ds(3)
    b= ddyup (ds, a)
END FUNCTION ddyup1

!*******************************************************************************
FUNCTION ddzup1 (ds, a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
  real(8):: ds(3)
    b= ddzup (ds, a)
END FUNCTION ddzup1

!*******************************************************************************
FUNCTION xdn1 (a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
    b = xdn (a)
END FUNCTION xdn1

!*******************************************************************************
FUNCTION ydn1 (a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
    b = ydn (a)
END FUNCTION ydn1

!*******************************************************************************
FUNCTION zdn1 (a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
    b = zdn (a)
END FUNCTION zdn1

!*******************************************************************************
FUNCTION xup1 (a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
    b = xup (a)
END FUNCTION xup1

!*******************************************************************************
FUNCTION yup1 (a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
    b = yup (a)
END FUNCTION yup1

!*******************************************************************************
FUNCTION zup1 (a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
    b = zup (a)
END FUNCTION zup1

!*******************************************************************************
FUNCTION ddx (ds, a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
  real(8):: ds, c
  integer :: ix, iy, iz, n(3)
  integer, save:: itimer=0
!...............................................................................
  !call trace%begin ('stagger_mod::ddx', itimer=itimer)
  n = shape(a)
!  if (hardwire) then
!    if (all(n==[16,16,16])) then
!      call dxdn_16 (ds,a,b)
!    else if (all(n==[20,20,20])) then
!      call dxdn_20 (ds,a,b)
!    else if (all(n==[24,24,24])) then
!      call dxdn_24 (ds,a,b)
!    else if (all(n==[32,32,32])) then
!      call dxdn_32 (ds,a,b)
!    else if (all(n==[36,36,36])) then
!      call dxdn_36 (ds,a,b)
!    end if
!  else if (n(1) > 1) then
  if (n(1) > 1) then
    b(   1,:,:) = 0.0
    b(n(1),:,:) = 0.0
    c = 0.5/ds
    do iz=1,n(3)
    do iy=1,n(2)
    !$vector always assert
    do ix=2,n(1)-1
      b(ix,iy,iz) = c*(a(ix+1,iy,iz) - a(ix-1,iy,iz))
    end do
    end do
    end do
  else
    b = 0.0
  endif
  !call trace%end (itimer)
  stagger%flops = stagger%flops + 2*n(2)*n(3)*(n(1)-1)
  stagger%count = stagger%count + 1
END FUNCTION ddx

!*******************************************************************************
FUNCTION ddy (ds, a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
  real(8):: ds, c
  integer :: ix, iy, iz, n(3)
  integer, save:: itimer=0
!...............................................................................
  !call trace%begin ('stagger_mod::ddy', itimer=itimer)
  n = shape(a)
!  if (hardwire) then
!    if (all(n==[16,16,16])) then
!      call dxdn_16 (ds,a,b)
!    else if (all(n==[20,20,20])) then
!      call dxdn_20 (ds,a,b)
!    else if (all(n==[24,24,24])) then
!      call dxdn_24 (ds,a,b)
!    else if (all(n==[32,32,32])) then
!      call dxdn_32 (ds,a,b)
!    else if (all(n==[36,36,36])) then
!      call dxdn_36 (ds,a,b)
!    end if
!  else if (n(1) > 1) then
  if (n(2) > 1) then
    b(:,   1,:) = 0.0
    b(:,n(2),:) = 0.0
    c = 0.5/ds
    do iz=1,n(3)
    do iy=2,n(2)-1
    !$vector always assert
    do ix=1,n(1)
      b(ix,iy,iz) = c*(a(ix,iy+1,iz) - a(ix,iy-1,iz))
    end do
    end do
    end do
  else
    b = 0.0
  endif
  !call trace%end (itimer)
  stagger%flops = stagger%flops + 2*n(2)*n(3)*(n(1)-1)
  stagger%count = stagger%count + 1
END FUNCTION ddy

!*******************************************************************************
FUNCTION ddz (ds, a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
  real(8):: ds, c
  integer :: ix, iy, iz, n(3)
  integer, save:: itimer=0
!...............................................................................
  !call trace%begin ('stagger_mod::ddz', itimer=itimer)
  n = shape(a)
!  if (hardwire) then
!    if (all(n==[16,16,16])) then
!      call dxdn_16 (ds,a,b)
!    else if (all(n==[20,20,20])) then
!      call dxdn_20 (ds,a,b)
!    else if (all(n==[24,24,24])) then
!      call dxdn_24 (ds,a,b)
!    else if (all(n==[32,32,32])) then
!      call dxdn_32 (ds,a,b)
!    else if (all(n==[36,36,36])) then
!      call dxdn_36 (ds,a,b)
!    end if
!  else if (n(1) > 1) then
  if (n(3) > 1) then
    b(:,:,   1) = 0.0
    b(:,:,n(3)) = 0.0
    c = 0.5/ds
    do iz=2,n(3)-1
    do iy=1,n(2)
    !$vector always assert
    do ix=1,n(1)
      b(ix,iy,iz) = c*(a(ix,iy,iz+1) - a(ix,iy,iz-1))
    end do
    end do
    end do
  else
    b = 0.0
  endif
  !call trace%end (itimer)
  stagger%flops = stagger%flops + 2*n(2)*n(3)*(n(1)-1)
  stagger%count = stagger%count + 1
END FUNCTION ddz

!*******************************************************************************
FUNCTION ddxdn (ds, a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
  real(8):: ds(3), c
  integer :: ix, iy, iz, n(3)
  integer, save:: itimer=0
!...............................................................................
!  call trace%begin ('stagger_mod::ddxdn', itimer=itimer)
  n = shape(a)
!  if (hardwire) then
!    if (all(n==[16,16,16])) then
!      call dxdn_16 (ds,a,b)
!    else if (all(n==[20,20,20])) then
!      call dxdn_20 (ds,a,b)
!    else if (all(n==[24,24,24])) then
!      call dxdn_24 (ds,a,b)
!    else if (all(n==[32,32,32])) then
!      call dxdn_32 (ds,a,b)
!    else if (all(n==[36,36,36])) then
!      call dxdn_36 (ds,a,b)
!    end if
!  else if (n(1) > 1) then
  if (n(1) > 1) then
    b(1,:,:) = 0.0
    c = 1./ds(1)
    do iz=1,n(3)
    do iy=1,n(2)
    !$vector always assert
    do ix=2,n(1)
      b(ix,iy,iz) = c*(a(ix  ,iy,iz) - a(ix-1,iy,iz))
    end do
    end do
    end do
  else
    b = 0.0
  endif
!  call trace%end (itimer)
  stagger%flops = stagger%flops + 2*n(2)*n(3)*(n(1)-1)
  stagger%count = stagger%count + 1
END FUNCTION ddxdn

!*******************************************************************************
FUNCTION ddydn (ds, a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
  integer :: ix, iy, iz, n(3)
  real(8):: ds(3), c
  integer, save:: itimer=0
!...............................................................................
!  call trace%begin ('stagger_mod::ddydn', itimer=itimer)
  n = shape(a)
!  if (hardwire) then
!    if (all(n==[16,16,16])) then
!      call dydn_16 (ds,a,b)
!    else if (all(n==[20,20,20])) then
!      call dydn_20 (ds,a,b)
!    else if (all(n==[24,24,24])) then
!      call dydn_24 (ds,a,b)
!    else if (all(n==[32,32,32])) then
!      call dydn_32 (ds,a,b)
!    else if (all(n==[36,36,36])) then
!      call dydn_36 (ds,a,b)
!    end if
!  else if (n(2) > 1) then
  if (n(2) > 1) then
    b(:,1,:) = 0.0
    c = 1./ds(2)
    do iz=1,n(3)
    do iy=2,n(2)
    !$vector always assert
    do ix=1,n(1)
      b(ix,iy,iz) = c*(a(ix,iy  ,iz) - a(ix,iy-1,iz))
    end do
    end do
    end do
  else
    b = 0.0
  endif
!  call trace%end (itimer)
  stagger%flops = stagger%flops + 2*n(1)*n(3)*(n(2)-1)
  stagger%count = stagger%count + 1
END FUNCTION ddydn

!*******************************************************************************
FUNCTION ddzdn (ds, a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
  real(8):: ds(3), c
  integer :: ix, iy, iz, n(3)
  integer, save:: itimer=0
!...............................................................................
!  call trace%begin ('stagger_mod::ddzdn', itimer=itimer)
  n = shape(a)
!  if (hardwire) then
!    if (all(n==[16,16,16])) then
!      call dzdn_16 (ds,a,b)
!    else if (all(n==[20,20,20])) then
!      call dzdn_20 (ds,a,b)
!    else if (all(n==[24,24,24])) then
!      call dzdn_24 (ds,a,b)
!    else if (all(n==[32,32,32])) then
!      call dzdn_32 (ds,a,b)
!    else if (all(n==[36,36,36])) then
!      call dzdn_36 (ds,a,b)
!    end if
!  else if (n(3) > 1) then
  if (n(3) > 1) then
    b(:,:,1) = 0.0
    c = 1./ds(3)
    do iz=2,n(3)
    do iy=1,n(2)
    !$vector always assert
    do ix=1,n(1)
      b(ix,iy,iz) = c*(a(ix,iy,iz  ) - a(ix,iy,iz-1))
    end do
    end do
    end do
  else
    b = 0.0
  endif
!  call trace%end (itimer)
  stagger%flops = stagger%flops + 2*n(1)*n(2)*(n(3)-1)
  stagger%count = stagger%count + 1
END FUNCTION ddzdn

!*******************************************************************************
FUNCTION ddxup (ds, a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
  integer :: ix, iy, iz, n(3)
  real(8):: ds(3), c
  integer, save:: itimer=0
!...............................................................................
!  call trace%begin ('stagger_mod::ddxup', itimer=itimer)
  n = shape(a)
!  if (hardwire) then
!    if (all(n==[16,16,16])) then
!      call dxup_16 (ds,a,b)
!    else if (all(n==[20,20,20])) then
!      call dxup_20 (ds,a,b)
!    else if (all(n==[24,24,24])) then
!      call dxup_24 (ds,a,b)
!    else if (all(n==[32,32,32])) then
!      call dxup_32 (ds,a,b)
!    else if (all(n==[36,36,36])) then
!      call dxup_36 (ds,a,b)
!    end if
!  else if (n(1) > 1) then
  if (n(1) > 1) then
    b(n(1),:,:) = 0.0
    c = 1./ds(1)
    !print*,'mk'
    do iz=1,n(3)
    do iy=1,n(2)
    !$vector always assert
    do ix=1,n(1)-1
      b(ix,iy,iz) = c*(a(ix+1,iy,iz) - a(ix  ,iy,iz))
    end do
    end do
    end do
  else
    b = 0.0
  endif
!  call trace%end (itimer)
  stagger%flops = stagger%flops + 2*n(3)*n(2)*(n(1)-1)
  stagger%count = stagger%count + 1
END FUNCTION ddxup

!*******************************************************************************
FUNCTION ddyup (ds, a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
  integer :: ix, iy, iz, n(3)
  real(8):: ds(3), c
  integer, save:: itimer=0
!...............................................................................
!  call trace%begin ('stagger_mod::ddyup', itimer=itimer)
  n = shape(a)
!  if (hardwire) then
!    if (all(n==[16,16,16])) then
!      call dyup_16 (ds,a,b)
!    else if (all(n==[20,20,20])) then
!      call dyup_20 (ds,a,b)
!    else if (all(n==[24,24,24])) then
!      call dyup_24 (ds,a,b)
!    else if (all(n==[32,32,32])) then
!      call dyup_32 (ds,a,b)
!    else if (all(n==[36,36,36])) then
!      call dyup_36 (ds,a,b)
!    end if
!  else if (n(2) > 1) then ! x-derivative only non-zero iff x-dimension larger than 1
  if (n(2) > 1) then ! x-derivative only non-zero iff x-dimension larger than 1
    b(:,n(2),:) = 0.0
    c = 1./ds(2)
    do iz=1,n(3)
    do iy=1,n(2)-1
    !$vector always assert
    do ix=1,n(1)
      b(ix,iy,iz) = c*(a(ix,iy+1,iz) - a(ix,iy  ,iz))
    end do
    end do
    end do
  else
    b = 0.0
  endif
!  call trace%end (itimer)
  stagger%flops = stagger%flops + 2*n(1)*n(3)*(n(2)-1)
  stagger%count = stagger%count + 1
END FUNCTION ddyup

!*******************************************************************************
FUNCTION ddzup (ds, a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
  integer :: ix, iy, iz, n(3)
  real(8):: ds(3), c
  integer, save:: itimer=0
!...............................................................................
!  call trace%begin ('stagger_mod::ddzup', itimer=itimer)
  n = shape(a)
!  if (hardwire) then
!    if (all(n==[16,16,16])) then
!      call dzup_16 (ds,a,b)
!    else if (all(n==[20,20,20])) then
!      call dzup_20 (ds,a,b)
!    else if (all(n==[24,24,24])) then
!      call dzup_24 (ds,a,b)
!    else if (all(n==[32,32,32])) then
!      call dzup_32 (ds,a,b)
!    else if (all(n==[36,36,36])) then
!      call dzup_36 (ds,a,b)
!    end if
!  else if (n(3) > 1) then ! x-derivative only non-zero iff x-dimension larger than 1
  if (n(3) > 1) then ! x-derivative only non-zero iff x-dimension larger than 1
    b(:,:,n(3)) = 0.0
    c = 1./ds(3)
    do iz=1,n(3)-1
    do iy=1,n(2)
    !$vector always assert
    do ix=1,n(1)
      b(ix,iy,iz) = c*(a(ix,iy,iz+1) - a(ix,iy,iz  ))
    end do
    end do
    end do
  else
    b = 0.0
  endif
!  call trace%end (itimer)
  stagger%flops = stagger%flops + 2*n(1)*n(2)*(n(3)-1)
  stagger%count = stagger%count + 1
END FUNCTION ddzup

!*******************************************************************************
FUNCTION xdn (a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
  real    :: c
  integer :: ix, iy, iz, n(3)
  integer, save:: itimer=0
!...............................................................................
!  call trace%begin ('stagger_mod::xdn', itimer=itimer)
  n = shape(a)
  if (hardwire) then
    if (all(n==[16,16,16])) then
      call xdn_16 (a, b)
    else if (all(n==[20,20,20])) then
      call xdn_20 (a, b)
    else if (all(n==[24,24,24])) then
      call xdn_24 (a, b)
    else if (all(n==[32,32,32])) then
      call xdn_32 (a, b)
    else if (all(n==[36,36,36])) then
      call xdn_36 (a, b)
    end if
  else if (n(1) > 1) then ! x-derivative only non-zero iff x-dimension larger than 1
    c = 0.5
    b(1,:,:) = a(1,:,:)
    do iz=1,n(3)
    do iy=1,n(2)
    !$vector always assert
    do ix=2,n(1)
      b(ix,iy,iz) = c*(a(ix  ,iy,iz) + a(ix-1,iy,iz))
    end do
    end do
    end do
  else
    b = a
  end if
!  call trace%end (itimer)
  stagger%flops = stagger%flops + 2*n(2)*n(3)*(n(1)-1)
  stagger%count = stagger%count + 1
END FUNCTION xdn

!*******************************************************************************
FUNCTION ydn (a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
  real    :: c
  integer :: ix, iy, iz, n(3)
  integer, save:: itimer=0
!...............................................................................
!  call trace%begin ('stagger_mod::ydn', itimer=itimer)
  n = shape(a)
  if (hardwire) then
    if (all(n==[16,16,16])) then
      call ydn_16 (a, b)
    else if (all(n==[20,20,20])) then
      call ydn_20 (a, b)
    else if (all(n==[24,24,24])) then
      call ydn_24 (a, b)
    else if (all(n==[32,32,32])) then
      call ydn_32 (a, b)
    else if (all(n==[36,36,36])) then
      call ydn_36 (a, b)
    end if
  else if (n(2) > 1) then ! x-derivative only non-zero iff x-dimension larger than 1
    c = 0.5
    b(:,1,:) = a(:,1,:)
    do iz=1,n(3)
    do iy=2,n(2)
    !$vector always assert
    do ix=1,n(1)
      b(ix,iy,iz) = c*(a(ix,iy,iz) + a(ix,iy-1,iz))
    end do
    end do
    end do
  else
    b = a
  end if
!  call trace%end (itimer)
  stagger%flops = stagger%flops + 2*n(1)*n(3)*(n(2)-1)
  stagger%count = stagger%count + 1
END FUNCTION ydn

!*******************************************************************************
FUNCTION zdn (a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
  real    :: c
  integer :: ix, iy, iz, n(3)
  integer, save:: itimer=0
!...............................................................................
!  call trace%begin ('stagger_mod::zdn', itimer=itimer)
  n = shape(a)
  if (hardwire) then
    if (all(n==[16,16,16])) then
      call zdn_16 (a, b)
    else if (all(n==[20,20,20])) then
      call zdn_20 (a, b)
    else if (all(n==[24,24,24])) then
      call zdn_24 (a, b)
    else if (all(n==[32,32,32])) then
      call zdn_32 (a, b)
    else if (all(n==[36,36,36])) then
      call zdn_36 (a, b)
    end if
  else if (n(3) > 1) then ! x-derivative only non-zero iff x-dimension larger than 1
    c = 0.5
    b(:,:,1) = a(:,:,1)
    do iz=2,n(3)
    do iy=1,n(2)
    !$vector always assert
    do ix=1,n(1)
      b(ix,iy,iz) = c*(a(ix,iy,iz) + a(ix,iy,iz-1))
    end do
    end do
    end do
  else
    b = a
  end if
!  call trace%end (itimer)
  stagger%flops = stagger%flops + 2*n(1)*n(2)*(n(3)-1)
  stagger%count = stagger%count + 1
END FUNCTION zdn

!*******************************************************************************
FUNCTION xup (a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
  real    :: c
  integer :: ix, iy, iz, n(3)
  integer, save:: itimer=0
!...............................................................................
!  call trace%begin ('stagger_mod::xup', itimer=itimer)
  n = shape(a)
  if (hardwire) then
    if (all(n==[16,16,16])) then
      call xup_16 (a, b)
    else if (all(n==[20,20,20])) then
      call xup_20 (a, b)
    else if (all(n==[24,24,24])) then
      call xup_24 (a, b)
    else if (all(n==[32,32,32])) then
      call xup_32 (a, b)
    else if (all(n==[36,36,36])) then
      call xup_36 (a, b)
    end if
  else if (n(1) > 1) then ! x-derivative only non-zero iff x-dimension larger than 1
    c = 0.5
    b(n(1),:,:) = a(n(1),:,:)
    do iz=1,n(3)
    do iy=1,n(2)
    !$vector always assert
    do ix=1,n(1)-1
      b(ix,iy,iz) = c*(a(ix+1,iy,iz) + a(ix  ,iy,iz))
    end do
    end do
    end do
  else
    b = a
  endif
!  call trace%end (itimer)
  stagger%flops = stagger%flops + 2*n(2)*n(3)*(n(1)-1)
  stagger%count = stagger%count + 1
END FUNCTION

!*******************************************************************************
FUNCTION yup (a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
  real    :: c
  integer :: ix, iy, iz, n(3)
  integer, save:: itimer=0
!...............................................................................
!  call trace%begin ('stagger_mod::yup', itimer=itimer)
  n = shape(a)
  if (hardwire) then
    if (all(n==[16,16,16])) then
      call yup_16 (a, b)
    else if (all(n==[20,20,20])) then
      call yup_20 (a, b)
    else if (all(n==[24,24,24])) then
      call yup_24 (a, b)
    else if (all(n==[32,32,32])) then
      call yup_32 (a, b)
    else if (all(n==[36,36,36])) then
      call yup_36 (a, b)
    end if
  else if (n(2) > 1) then ! x-derivative only non-zero iff x-dimension larger than 1
    c = 0.5
    b(:,n(2),:) = a(:,n(2),:)
    do iz=1,n(3)
    do iy=1,n(2)-1
    !$vector always assert
    do ix=1,n(1)
      b(ix,iy,iz) = c*(a(ix,iy+1,iz) + a(ix,iy  ,iz))
    end do
    end do
    end do
  else
    b = a
  endif
!  call trace%end (itimer)
  stagger%flops = stagger%flops + 2*n(1)*n(3)*(n(2)-1)
  stagger%count = stagger%count + 1
END FUNCTION

!*******************************************************************************
FUNCTION zup (a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
  real    :: c
  integer :: ix, iy, iz, n(3)
  integer, save:: itimer=0
!...............................................................................
!  call trace%begin ('stagger_mod::zup', itimer=itimer)
  n = shape(a)
  if (hardwire) then
    if (all(n==[16,16,16])) then
      call zup_16 (a, b)
    else if (all(n==[20,20,20])) then
      call zup_20 (a, b)
    else if (all(n==[24,24,24])) then
      call zup_24 (a, b)
    else if (all(n==[32,32,32])) then
      call zup_32 (a, b)
    else if (all(n==[36,36,36])) then
      call zup_36 (a, b)
    end if
  else if (n(3) > 1) then ! x-derivative only non-zero iff x-dimension larger than 1
    c = 0.5
    b(:,:,n(3)) = a(:,:,n(3))
    do iz=1,n(3)-1
    do iy=1,n(2)
    !$vector always assert
    do ix=1,n(1)
      b(ix,iy,iz) = c*(a(ix,iy,iz+1) + a(ix,iy,iz  ))
    end do
    end do
    end do
  else
    b = a
  endif
!  call trace%end (itimer)
  stagger%flops = stagger%flops + 2*n(1)*n(2)*(n(3)-1)
  stagger%count = stagger%count + 1
END FUNCTION

!*******************************************************************************
FUNCTION xsm (a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
  real    :: c
  integer :: i, n
  integer, save:: itimer=0
!...............................................................................
!  call trace%begin ('stagger_mod::xsm', itimer=itimer)
  n = size(a,1)
  if (n > 1) then
    c = 0.25
    do i=2,n-1
      b(i,:,:) = (a(i-1,:,:) + a(i+1,:,:))*c +  a(i,:,:)*(1.-2.*c)
    end do
    do i=1,n,n-1
      b(i,:,:) = a(i,:,:)
    end do
  else
    b = a
  endif
!  call trace%end (itimer)
  stagger%flops = stagger%flops + 4*n*n*(n-2)
  stagger%count = stagger%count + 1
END FUNCTION

!*******************************************************************************
FUNCTION ysm (a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
  real    :: c
  integer :: i, n
  integer, save:: itimer=0
!...............................................................................
!  call trace%begin ('stagger_mod::ysm', itimer=itimer)
  n = size(a,2)
  if (n > 1) then
    c = 0.25
    do i=2,n-1
      b(:,i,:) = (a(:,i-1,:) + a(:,i+1,:))*c +  a(:,i,:)*(1.-2.*c)
    end do
    do i=1,n,n-1
      b(:,i,:)  = a(:,i,:)
    end do
  else
    b = a
  endif
!  call trace%end (itimer)
  stagger%flops = stagger%flops + 4*n*n*(n-2)
  stagger%count = stagger%count + 1
END FUNCTION

!*******************************************************************************
FUNCTION zsm (a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
  real    :: c
  integer :: i, n
  integer, save:: itimer=0
!...............................................................................
!  call trace%begin ('stagger_mod::zsm', itimer=itimer)
  n = size(a,3)
  if (n > 1) then
    c = 0.25
    do i=2,n-1
      b(:,:,i) = (a(:,:,i-1) + a(:,:,i+1))*c +  a(:,:,i)*(1.-2.*c)
    end do
    do i=1,n,n-1
      b(:,:,i) = a(:,:,i)
    end do
  else
    b = a
  endif
!  call trace%end (itimer)
  stagger%flops = stagger%flops + 4*n*n*(n-2)
  stagger%count = stagger%count + 1
END FUNCTION

!*******************************************************************************
FUNCTION xyzsm(a) RESULT (b)
  real, dimension(:,:,:), intent(in):: a
  real, dimension(size(a,1),size(a,2),size(a,3)):: b
  integer, save:: itimer=0
  call trace%begin ('stagger_mod::xyzsm', itimer=itimer)
  b = zsm(ysm(xsm(a)))
  call trace%end (itimer)
END FUNCTION

!===============================================================================
!> Unit test of 6th/5th order operators
!===============================================================================
SUBROUTINE stagger_test
  integer       :: m(3)
  real, pointer :: a(:,:,:)
  real, pointer :: b(:,:,:)
  real, pointer :: c(:,:,:)
  real(8)       :: ds(3), x, y, z, d
  integer       :: i
  logical       :: ok, allok
  real, parameter :: eps = 2e-6
  !-----------------------------------------------------------------------------
  m = [10,11,12]
  ds = 1d0/m
  allocate (a(m(1),m(2),m(3)))
  allocate (b(m(1),m(2),m(3)))
  allok = .true.
  !-----------------------------------------------------------------------------
  ! x-operators
  !-----------------------------------------------------------------------------
  do i=1,m(1)
    x = i*ds(1)
    a(i,:,:) = x
  end do
  b = xdn(a)
  ok = .true.
  do i=4,m(1)-2
    x = (i-0.5d0)*ds(1)
    d = x
    ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
  end do
  allok = allok .and. ok
  if (.not.ok) print *,'xdn:', ok
  b = xup(a)
  ok = .true.
  do i=3,m(1)-3
    x = (i+0.5d0)*ds(1)
    d = x
    ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
  end do
  allok = allok .and. ok
  if (.not.ok) print *,'xup:', ok
  !-----------------------------------------------------------------------------
  ! y-operators
  !-----------------------------------------------------------------------------
  do i=1,m(2)
    x = i*ds(2)
    a(:,i,:) = x
  end do
  b = ydn(a)
  ok = .true.
  do i=4,m(2)-2
    y = (i-0.5d0)*ds(2)
    d = y
    ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
  end do
  allok = allok .and. ok
  if (.not.ok) print *,'ydn:', ok
  b = yup(a)
  ok = .true.
  do i=3,m(2)-3
    y = (i+.5d0)*ds(2)
    d = y
    ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
  end do
  allok = allok .and. ok
  if (.not.ok) print *,'yup:', ok
  !-----------------------------------------------------------------------------
  ! z-operators
  !-----------------------------------------------------------------------------
  do i=1,m(3)
    x = i*ds(3)
    a(:,:,i) = x
  end do
  b = zdn(a)
  ok = .true.
  do i=4,m(3)-2
    z = (i-0.5d0)*ds(3)
    d = z
    ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
  end do
  allok = allok .and. ok
  if (.not.ok) print *,'zdn:', ok
  b = zup(a)
  ok = .true.
  do i=3,m(3)-3
    z = (i+0.5d0)*ds(3)
    d = z
    ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
  end do
  allok = allok .and. ok
  if (.not.ok) print *,'zup:', ok
  !-----------------------------------------------------------------------------
  ! x-operators
  !-----------------------------------------------------------------------------
  do i=1,m(1)
    x = i*ds(1)
    a(i,:,:) = 1d0 + x**2
  end do
  b = ddxdn(ds,a)
  ok = .true.
  do i=4,m(1)-2
    x = (i-0.5d0)*ds(1)
    d = 2d0*x
    ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
  end do
  allok = allok .and. ok
  if (.not.ok) print *,'ddxdn:', ok
  b = ddxup(ds,a)
  ok = .true.
  do i=3,m(1)-3
    x = (i+0.5d0)*ds(1)
    d = 2d0*x
    ok = ok .and. (maxval(abs(b(i,:,:)-d)) < eps)
  end do
  allok = allok .and. ok
  if (.not.ok) print *,'ddxup:', ok
  !-----------------------------------------------------------------------------
  ! y-operators
  !-----------------------------------------------------------------------------
  do i=1,m(2)
    x = i*ds(2)
    a(:,i,:) = 1d0 + x**2
  end do
  b = ddydn(ds,a)
  ok = .true.
  do i=4,m(2)-2
    y = (i-0.5d0)*ds(2)
    d = 2d0*y
    ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
  end do
  allok = allok .and. ok
  if (.not.ok) print *,'ddydn:', ok
  b = ddyup(ds,a)
  ok = .true.
  do i=3,m(2)-3
    y = (i+.5d0)*ds(2)
    d = 2d0*y
    ok = ok .and. (maxval(abs(b(:,i,:)-d)) < eps)
  end do
  allok = allok .and. ok
  if (.not.ok) print *,'ddyup:', ok
  !-----------------------------------------------------------------------------
  ! z-operators
  !-----------------------------------------------------------------------------
  do i=1,m(3)
    x = i*ds(3)
    a(:,:,i) = 1d0 + x**2
  end do
  b = ddzdn(ds,a)
  ok = .true.
  do i=4,m(3)-2
    z = (i-0.5d0)*ds(3)
    d = 2d0*z
    ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
  end do
  allok = allok .and. ok
  if (.not.ok) print *,'ddzdn:', ok
  b = ddzup(ds,a)
  ok = .true.
  do i=3,m(3)-3
    z = (i+0.5d0)*ds(3)
    d = 2d0*z
    ok = ok .and. (maxval(abs(b(:,:,i)-d)) < eps)
  end do
  allok = allok .and. ok
  if (.not.ok) print *,'ddzup:', ok
  print *,io%hl
  if (allok) then
    print *,'2nd/1st order stagger operator test passed'
    print *,io%hl
  else
    print *,'2nd/1st order stagger operator test failed!'
    print *,io%hl
    error stop 'stagger_test'
  end if
  !-----------------------------------------------------------------------------
  deallocate (a, b)
END SUBROUTINE stagger_test

END MODULE

max_stars_repo_name: applejwjcat/dispatch

max_issues_repo_name: applejwjcat/dispatch

max_forks_repo_name: applejwjcat/dispatch


text:
C     Test of .NE. handling in the false branch

      subroutine negneg(i,j)

      if(i.ne.j) then
         i = j
      else
         k = 3
      endif

      print *, i, j, k

      end


text:
! RUN: not %flang_fc1 -fdebug-unparse-with-symbols %s 2>&1 | FileCheck %s
! CHECK: '60' not a FORMAT
! CHECK: data transfer use of '60'

subroutine s(a)
  real a(10)
  write(*,60) "Hi there"
60 continue
70 format (i8)
end subroutine s

max_stars_repo_name: mkinsner/llvm

max_issues_repo_name: mkinsner/llvm

max_forks_repo_name: mkinsner/llvm


text:
	SUBROUTINE HINITA  ( dev, cdev, iunit, filnam, itype,
     +			     xsize, ysize, ncurwn, iret )
C************************************************************************
C* HINITA - RBK 				                        *
C*									*
C* This subroutine is called to initialize a new device driver.		*
C*									*
C* HINITA  ( DEV, CDEV, IUNIT, FILNAM, ITYPE, XSIZE, YSIZE, NCURWN,	*
C*	     IRET )							*
C*									*
C* Input parameters:							*
C*	DEV		CHAR*		Device name			*
C*	CDEV            CHAR*           Current device name		*
C*	IUNIT		INTEGER		Type of output device		*
C*					  (not used for UTF)	        *
C*	FILNAM		CHAR*		File name of output		*
C* 	ITYPE		INTEGER		Output file format		*
C*					    2 = NAFOS format		*
C*					  <>2 = AFOS format		*
C*	XSIZE		REAL		X size in pixels		*
C*	YSIZE		REAL		Y size in pixels		*
C*									*
C* Output parameters:							*
C*	NCURWN		INTEGER		Current window number		*
C*					  (not used for UTF)		*
C*	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* A. Hardy/GSC		9/98		Modified from utf's HINITA      *
C* A. Hardy/GSC		9/00		Modified base time YYMMDD/HHNN  *
C* B. Yin/SAIC          3/04   		Changed SS_GTIM to CSS_GTIM     *
C************************************************************************
	INCLUDE         'GEMPRM.PRM'
	INCLUDE         'ERROR.PRM'
	INCLUDE         'DVWNDW.CMN'
	INCLUDE         'DEVCHR.CMN'
	INCLUDE         'DEVACT.CMN'
C*
	CHARACTER*(*)	dev, cdev, filnam
C*
	CHARACTER	carr(2)*160, fname*160, ttime*20, ctime*20
	CHARACTER	hourst*4, datest*2
	INTEGER         ibase(5), itime
C------------------------------------------------------------------------
	iret   = NORMAL
	ncurwn = 0
C
C*	Save device common area.
C
	ddev   = dev
	curdev = cdev
	niunit = iunit
C
C*	Initialize device characteristics common area for specific
C*	device.
C
	CALL HINIT ( iret)
C
C*	Check to make sure bscalh has been set in device driver. If
C*	it has not, set it to a reasonable value.
C
	IF  ((bscalh .le. 0.) .or. (bscalh .gt. 1000.)) bscalh = 3.0
C
C*	Initialize the UTF output file.
C
	CALL ST_CLST ( filnam, ';', ' ', 2, carr, num, ierr )
	CALL ST_NULL ( carr(1), fname, lenf, ier )
C
C*	Retrieving full GEMPAK date/time base time string.
C*      If base time is given as DDHHNN convert to DD/HHNN.
C
        CALL ST_NOCC ( carr(2), '/', 1, ipntr, ier)
        IF ( ipntr .eq. 0 ) THEN
	    CALL ST_LSTR ( carr(2), leng, ier)
	    IF ( leng .eq. 6 ) THEN
              datest = carr(2)(1:2)
              hourst = carr(2)(3:6)
 	      ELSE IF (leng .eq. 5) THEN
              datest = '0' // carr(2)(1:1)
              hourst = carr(2)(2:5)
            END IF
            carr(2) = datest// '/' //hourst
        END IF
C
	itime = 1
        CALL CSS_GTIM ( itime, ctime, ier)
        CALL TI_STAN ( carr(2), ctime, ttime, ierr)
        IF ((ierr .eq. 0 ) .and. ( ttime(3:4) .le. '12') ) THEN
	    CALL TI_CTOI ( ttime, ibase, ier)
	  ELSE
	    DO i = 1,5
                ibase(i) = IMISSD
            END DO
            CALL ER_WMSG ('RBK', 1, ' ', ier)
        END IF
C
	CALL AINITA ( iunit, fname, lenf, itype, ibase, xsize, ysize,
     +		      ileft, ibot, iright, itop, nncolr, iret )
C
C*      isxoff = x offset of the 'S' coord relative to 'D' coord
C*      isyoff = y offset of the 'S' coord relative to 'D' coord
C*      iswdht = width of the 'S' coord system
C*      ishght = height of the  'S' coord system
C*      not used in PS driver
C
        isxoff  = 0
        isyoff  = 0
        iswdth  = ABS ( iright - ileft )
        ishght  = ABS ( itop - ibot )
C
C*	Do not change anything below this line.
C------------------------------------------------------------------------
	ispanx = ISIGN ( 1, ( iright - ileft ) )
	ispany = ISIGN ( 1, ( itop - ibot ) )
C*
	RETURN
	END

max_forks_repo_name: oxelson/gempak