text:
!----- LGPL --------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This library is free software; you can redistribute it and/or
!  modify it under the terms of the GNU Lesser General Public
!  License as published by the Free Software Foundation version 2.1.
!
!  This library is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!  Lesser General Public License for more details.
!
!  You should have received a copy of the GNU Lesser General Public
!  License along with this library; if not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: tree_struct.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/utils_lgpl/deltares_common/packages/deltares_common/src/tree_struct.f90 $
! tree_struct.f90 --
!    Module that implements a general tree structure in Fortran 90
!
! General information:
!    The tree is stored via a derived type TREE_DATA. A tree can
!    have an arbitrary number of nodes, each of which can again
!    have an arbitrary number of subnodes and so on.
!    The data type for a node is the same as for a tree - there
!    is no difference, except that the routine tree_create()
!    returns the initialised root of a new tree.
!    A node (and the root as well) can have the following
!    properties:
!    - A name (at most 80 characters)
!    - Arbitrary data
!    - A string indicating the type of data
!    - Zero, one or more subnodes
!    The storage is fairly efficient (via the TRANSFER() function
!    all data and strings are converted to arrays of default
!    integers). A node does not "know" its parent though and
!    there is no check on circularity.
!
module TREE_DATA_TYPES
   type TREE_DATA
      character(len=1), dimension(:), pointer         :: node_name
      character(len=1), dimension(:), pointer         :: node_data
      character(len=1), dimension(:), pointer         :: node_data_type
      integer                                         :: node_visit     !< Zeroed upon construction, incremented upon node_data request (properties.f90: prop_get_string)
      type(TREE_DATA_PTR), dimension(:), pointer :: child_nodes
   end type

   type TREE_DATA_PTR
      type(TREE_DATA), pointer               :: node_ptr
   end type
end module

module TREE_STRUCTURES
   use TREE_DATA_TYPES
   implicit none

   private
   !
   ! A variable that indicates the type of all data
   !
   character(len=1), dimension(1:1), public,save :: node_value

   !
   ! Auxiliary variable
   !
   integer, parameter                            :: maxlen = 300
   integer,save                                  :: traverse_level = 0

   !
   ! Public routines, types and parameters
   !
   public  :: TREE_DATA
   public  :: tree_create, tree_create_node, tree_add_node, tree_get_node_by_name, tree_num_nodes, &
              tree_get_data_ptr, tree_put_data, tree_get_name, tree_get_data,                      &
              tree_get_datatype, tree_get_data_string,                                             &
              tree_traverse, tree_traverse_level, print_tree,                                      &
              tree_fold, tree_destroy
   ! nested function has to be public for gfortran
   public ::  dealloc_tree_data

contains

! tree_create --
!    Create a new tree
!
! Arguments:
!    name         Name of the new tree
!    tree         Pointer to the new tree
! Result:
!    The argument tree points to a new, empty tree structure or is
!    not associated
!
subroutine tree_create( name, tree )
   character(len=*), intent(in)    :: name
   type(TREE_DATA), pointer        :: tree

   integer                         :: error
   integer                         :: newsize

!   GD: memory leak here
!   if(associated(tree)) then
!     deallocate(tree)
!   end if

   allocate( tree, stat = error )

   if ( error .ne. 0 ) then
      nullify( tree )
   else
      newsize = size( transfer( name, node_value ) )
      !GD: memory leak here
      !if(associated(tree%node_name)) deallocate(tree%node_name)
      allocate( tree%node_name(1:newsize), stat = error )
      if ( error .ne. 0 ) then
         deallocate( tree )
         return
      else
         tree%node_name(1:newsize) = transfer( name, node_value )
         tree%node_visit = 0
         nullify( tree%node_data )
         nullify( tree%node_data_type )
         nullify( tree%child_nodes )
      endif
   endif
end subroutine tree_create

! tree_create_node --
!    Create a new node to the given tree or node
!
! Arguments:
!    tree         The tree or node to which to append the new node
!    name         Name of the new node
!    node         Pointer to the new node
! Result:
!    The argument node points to a new, empty node or is
!    not associated
!
subroutine tree_create_node( tree, name, node )
   character(len=*), intent(in)    :: name
   type(TREE_DATA), pointer        :: tree
   type(TREE_DATA), pointer        :: node

   integer :: ierror

   !
   ! Check for uniqueness
   !
!   call tree_get_node_by_name( tree, name, node )
!   if ( associated( node ) ) then
!      return
!   endif

   !
   ! Create a new node, store it in the array of child nodes
   ! for this (sub)tree
   !
   call tree_create( name, node )

   if ( associated( node ) ) then
      call tree_add_node(tree, node, ierror)
      if (ierror /= 0) then
         deallocate(node)
      end if
   endif
end subroutine tree_create_node

!> Adds an existing tree node to the children array of a tree.
!! Both the tree and the new node are pointers, use this to efficiently
!! create or extend a tree with already existing subtrees.
subroutine tree_add_node(tree, node, ierror)
   type(TREE_DATA), pointer        :: tree   !< Pointer to the root of an existing tree, to which the node should be added.
   type(TREE_DATA), pointer        :: node   !< Pointer to an existing ('sub')tree, which will be added to the root tree.
   integer,            intent(out) :: ierror !< Error status, 0 if succesful.

   type(TREE_DATA_PTR), dimension(:), pointer :: children

   integer                         :: newsize

   ierror = 0

   if (.not. associated(tree)) then
      ierror = 1
      return
   end if

   if ( associated( node ) ) then
      newsize = 1
      if ( associated( tree%child_nodes ) ) then
         newsize = 1 + size( tree%child_nodes )
      endif

      allocate( children(1:newsize), stat = ierror )
      if ( ierror .ne. 0 ) then
         return
      else
         if ( newsize .gt. 1 ) then
            children(1:newsize-1) = tree%child_nodes
            deallocate( tree%child_nodes )
         endif

         tree%child_nodes => children
         tree%child_nodes(newsize)%node_ptr => node
      endif
   else
      ierror = 2
   endif
end subroutine tree_add_node

!> Returns the number of nodes in a tree.
!! This node count is NOT recursive: it represents the number of nodes directly under the root level.
!! When the tree pointer itself is not associated, 0 is returned.
function tree_num_nodes(tree) result(num_nodes)
   type(TREE_DATA), pointer :: tree !< Tree pointer for which to determine the number of child nodes.

   integer                  :: num_nodes !< Number of child nodes in the specified tree, or 0 if unassociated.

   num_nodes = 0

   if (associated(tree)) then
      if (associated(tree%child_nodes)) then
         num_nodes = size(tree%child_nodes)
      end if
   end if
end function tree_num_nodes

! tree_get_name --
!    Return the name of the tree or node
!
! Arguments:
!    tree         The tree or node
!
function tree_get_name( tree ) result( node_name )
   type(TREE_DATA), pointer :: tree
   character(len=80)        :: node_name

   integer                  :: length
   integer                  :: i

   length    = min(80, size( tree%node_name ))
   node_name = ' '
   do i=1,length
      node_name(i:i) =  tree%node_name(i)
   end do
end function tree_get_name

! tree_get_data --
!    Return the data of the tree or node
!
! Arguments:
!    tree         The tree or node
!
function tree_get_data( tree ) result( node_data )
   type(TREE_DATA), pointer :: tree
   character(len=80)        :: node_data

   integer                  :: length
   integer                  :: i

   length    = min(80, size( tree%node_data ))
   node_data = ' '
   do i=1,length
      node_data(i:i) =  tree%node_data(i)
   end do
end function tree_get_data

! tree_get_datatype --
!    Return the data type for the data stored in the tree or node
!
! Arguments:
!    tree         The tree or node
!
function tree_get_datatype( tree ) result( data_type )
   type(TREE_DATA), pointer        :: tree
   character(len=40)               :: data_type

   integer                         :: length
   integer                         :: i

   data_type = '?'
   if ( associated( tree%node_data_type ) ) then
      length    = size( tree%node_data_type )
      do i=1,length
         data_type(i:i) = tree%node_data_type(i)
      end do
   endif
end function tree_get_datatype

! tree_get_node_by_name --
!    Return the child node by name
!
! Arguments:
!    tree         The tree or node to which to append the new node
!    name         Name of the node to find
!    node         Pointer to the node or "null"
! Result:
!    The argument node points to a new, empty node or is
!    not associated
!
subroutine tree_get_node_by_name( tree, name, node, i_return )
   character(len=*), intent(in)    :: name
   type(TREE_DATA), pointer        :: tree
   type(TREE_DATA), pointer        :: node

   character(len=80)               :: node_name
   character(len=80)               :: low_name
   integer, optional               :: i_return

   integer                         :: i

   if (present(i_return)) i_return = 0
   nullify( node )
   low_name = name
   call lowercase(low_name,999)

   node_name = tree_get_name( tree )

   if ( node_name .eq. low_name ) then
      node => tree
   elseif ( associated(tree%child_nodes) ) then
      do i = 1,size(tree%child_nodes)
         node_name = tree_get_name( tree%child_nodes(i)%node_ptr )
         call lowercase(node_name,999)

         if ( node_name .eq. low_name ) then
            node => tree%child_nodes(i)%node_ptr
            if (present(i_return)) i_return = i
            exit
         endif
      enddo
   endif

end subroutine tree_get_node_by_name

! tree_get_data_ptr --
!    Return a pointer to the tree/node's data
!
! Arguments:
!    tree        The tree or node from which to get the data
!    data_ptr    Pointer to the node/tree data
!    data_type   String indicating the type
! Result:
!    The argument data_ptr points to the stored data or is
!    not associated
!
subroutine tree_get_data_ptr( tree, data_ptr, data_type )
    type(TREE_DATA), pointer        :: tree
    character(len=1), dimension(:), pointer  :: data_ptr
    character(len=*)                         :: data_type

    nullify( data_ptr )

    data_type = '?'
    if ( associated( tree%node_data) ) then
       data_ptr  => tree%node_data
       data_type =  tree_get_datatype( tree )
    endif

end subroutine tree_get_data_ptr

! tree_put_data --
!    Put (a copy of) the data in the tree/node
!
! Arguments:
!    tree        The tree or node with which to attach the data
!    data        Array of integers
!    data_type   Optional string indicating the type
!    success     True if all went well, false otherwise
! Result:
!    The tree structure points to a copy of the data
! Note:
!    A direct call to this routine will look something like:
!
!       call tree_put_data( tree, transfer( some_data, node_value ) )
!
!    where node_value acts as the mold for transferring the data
!
subroutine tree_put_data( tree, data, data_type, success )
    type(TREE_DATA), pointer        :: tree
    character(len=1), dimension(:)  :: data
    character(len=*), optional      :: data_type
    logical, intent(out), optional  :: success

    integer                         :: error

    if ( associated(tree%node_data) ) then
       deallocate( tree%node_data )
    endif

    if ( associated(tree%node_data_type) ) then
       deallocate( tree%node_data_type )
    endif

!    GD: memory leak
!    if(associated(tree%node_data)) deallocate(tree%node_data)
    allocate( tree%node_data(1:size(data)), stat = error )
    if ( error .eq. 0 ) then
       tree%node_data = data
       allocate( tree%node_data_type(1:len_trim(data_type)), &
          stat = error )
       if ( error .eq. 0 ) then
          tree%node_data_type = transfer( data_type, tree%node_data_type )
       endif
    endif

    if ( present( success ) ) then
       success = error .eq. 0
    endif

end subroutine tree_put_data

! tree_traverse_level --
!    Convenience function: level of the node during traversal
!
! Arguments:
!    None
!
integer function tree_traverse_level( )
   tree_traverse_level = traverse_level
end function tree_traverse_level

! tree_traverse --
!    Traverse a tree and handle the nodes by a depth-first method
!
! Arguments:
!    tree        The tree or node to traverse
!    handler     Routine to handle each node
!    data        Arbitrary data to be passed to the handler
!    stop        Whether to continue or stop (if set true)
! Result:
!    Each tree node is visited (unless the traversal is
!    prematurely ended by setting "stop" to true)
!
recursive subroutine tree_traverse( tree, handler, data, stop )
    type(TREE_DATA), pointer        :: tree
    character(len=1), dimension(:)  :: data
    logical, intent(out)            :: stop

    interface
       subroutine handler( node, data, stop )
          use TREE_DATA_TYPES
          type(TREE_DATA), pointer                    :: node
          character(len=1), dimension(:), intent(in)  :: data
          logical, intent(inout)                      :: stop
       end subroutine handler
    end interface

    integer                         :: i

    stop = .false.
    if ( .not. associated( tree ) ) then
       return
    endif

    !
    ! First call the handler for the current node/tree
    !
    call handler( tree, data, stop )
    if ( stop ) then
       return
    endif

    !
    ! Then recurse through the child nodes (if any)
    !
    if ( associated( tree%child_nodes) ) then
       do i = 1,size(tree%child_nodes)
          traverse_level = traverse_level + 1
          call tree_traverse( tree%child_nodes(i)%node_ptr, &
                              handler, data, stop )
          traverse_level = traverse_level - 1
          if ( stop ) then
             exit
          endif
       enddo
    endif

end subroutine tree_traverse

! tree_traverse_bottomup --
!    Traverse a tree and handle the nodes by a depth-first method.
!    The callback subroutine is first performed on the children of each
!    node and then on the node itself.
!
! Arguments:
!    tree        The tree or node to traverse
!    handler     Routine to handle each node
!    data        Arbitrary data to be passed to the handler
!    stop        Whether to continue or stop (if set true)
! Result:
!    Each tree node is visited (unless the traversal is
!    prematurely ended by setting "stop" to true)
!
recursive subroutine tree_traverse_bottomup( tree, handler, data, stop )
    type(TREE_DATA), pointer        :: tree
    character(len=1), dimension(:)  :: data
    logical, intent(out)            :: stop

    interface
       subroutine handler( node, data, stop )
          use TREE_DATA_TYPES
          type(TREE_DATA), pointer                    :: node
          character(len=1), dimension(:), intent(in)  :: data
          logical, intent(inout)                      :: stop
       end subroutine handler
    end interface

    integer                         :: i

    stop = .false.
    if ( .not. associated( tree ) ) then
       return
    endif

    !
    ! First recurse through the child nodes (if any)
    !
    if ( associated( tree%child_nodes) ) then
       do i = 1,size(tree%child_nodes)
          traverse_level = traverse_level + 1
          call tree_traverse_bottomup( tree%child_nodes(i)%node_ptr, &
                              handler, data, stop )
          traverse_level = traverse_level - 1
          if ( stop ) then
             exit
          endif
       enddo
    endif
    if ( stop ) then
       return
    endif

    !
    ! Then call the handler for the current node/tree
    !
    call handler( tree, data, stop )

end subroutine tree_traverse_bottomup

!> Destroys a tree freeing up all its memory. (don't use a nested subroutine)
subroutine tree_destroy(tree)
    type(TREE_DATA), pointer                   :: tree    !< Tree that should be destroyed.
    logical :: dummylog

    call tree_traverse_bottomup(tree, dealloc_tree_data, node_value, dummylog)
    nullify(tree)
end subroutine tree_destroy

!> Deallocates all node data for a tree root.
!! Assumes that all child nodes's data is already deallocated.
!! This subroutine is intended for use in tree_traverse_bottomup.
subroutine dealloc_tree_data(tree, data, stop)
    type(TREE_DATA), pointer                      :: tree !< Tree whose root should be printed.
    character(len=1), dimension(:), intent(in)    :: data !< Help data, not used, may be empty
    logical,                        intent(inout) :: stop !< Whether to continue or stop.

    integer :: error

    if (associated(tree)) then
        if (associated(tree%node_name)) then
            deallocate(tree%node_name, stat = error )
        endif
        if (associated(tree%node_data)) then
            deallocate(tree%node_data, stat = error )
        endif
        if (associated(tree%node_data_type)) then
            deallocate(tree%node_data_type, stat = error )
        endif
        if (associated(tree%child_nodes)) then
            deallocate(tree%child_nodes, stat = error )
        endif
        deallocate(tree, stat = error )
   end if
end subroutine dealloc_tree_data

!> 'Fold' a tree together, using operations on child data, in a bottomup fashion.
!! Two callback routines are required:
!! * one that operates on leaves and puts the result in the data(:) variable.
!! * one that operates on a tree node with the children already folded, using
!!   results in childdata(:,:).
!! The data and childdata can contain any type by  'transfer(.., node_value)'.
recursive subroutine tree_fold( tree, tree_handler, leaf_handler, data, stop )
    type(TREE_DATA), pointer        :: tree
    character(len=1), dimension(:)  :: data
    logical, intent(out)            :: stop

    interface
       subroutine tree_handler( node, childdata, data, stop )
          use TREE_DATA_TYPES
          type(TREE_DATA), pointer                        :: node
          character(len=1), dimension(:,:), intent(in)    :: childdata
          character(len=1), dimension(:),   intent(out)   :: data
          logical,                          intent(inout) :: stop
       end subroutine tree_handler
    end interface
    interface
       subroutine leaf_handler( node, data, stop )
          use TREE_DATA_TYPES
          type(TREE_DATA), pointer                        :: node
          character(len=1), dimension(:),   intent(out)   :: data
          logical,                          intent(inout) :: stop
       end subroutine leaf_handler
    end interface

    character(len=1), allocatable   :: childdata(:,:)

    integer                         :: i

    stop = .false.
    if ( .not. associated( tree ) ) then
       return
    endif

    !
    ! First recurse through the child nodes (if any)
    !
    if ( associated( tree%child_nodes) ) then
       allocate(childdata(size(data), size(tree%child_nodes)))
       do i = 1,size(tree%child_nodes)
          traverse_level = traverse_level + 1
          call tree_fold( tree%child_nodes(i)%node_ptr, &
                          tree_handler, leaf_handler, childdata(:,i), stop )
          traverse_level = traverse_level - 1
          if ( stop ) then
             exit
          endif
       enddo

       !
       ! Then call the handler for the childdata+current node/tree
       !
       call tree_handler( tree, childdata, data, stop )

       deallocate(childdata)
    else
       !
       ! Otherwise compute the data in a leaf node
       !
       call leaf_handler(tree, data, stop)
    endif

    if ( stop ) then
       return
    endif

end subroutine tree_fold

! tree_get_data_string --
!    Return data as a simple string
!
! Arguments:
!    tree        The tree or node from which to get the data
!    string      String to be filled
!    success     Whether successful or not
! Result:
!    The string is filled with the data stored in the node
!    not associated. The routine is successful if:
!    - there is data associated with the node/tree
!    - the data type is "STRING"
!    If the routine is not successful, the string is
!    not changed.
!
subroutine tree_get_data_string( tree, string, success )
   type(TREE_DATA), pointer                 :: tree
   character(len=*), intent(out)            :: string
   logical, intent(out)                     :: success

   character(len=1), dimension(:), pointer  :: data_ptr
   character(len=40)                        :: data_type
   integer                                  :: length
   integer                                  :: i

   success = .false.
   if ( associated(tree) ) then
      call tree_get_data_ptr( tree, data_ptr, data_type )

      if ( .not. associated(data_ptr) ) then
         return
      endif
      if ( data_type .ne. 'STRING' ) then
         return
      endif

      success          = .true.
      length           = size(data_ptr)
      string           = ' '
      if (length <= maxlen) then
         length           = min(length,len(string))
         do i=1,length
            string(i:i) = data_ptr(i)
         end do
      endif
   endif

end subroutine tree_get_data_string

subroutine print_tree( tree, data, stop )
   type(TREE_DATA), pointer               :: tree
   character(len=1), dimension(:), intent(in) :: data
   logical, intent(inout)                 :: stop

   character(len=1), dimension(:),pointer :: data_ptr
   character(len=60)                      :: string
   character(len=40)                      :: type_string

   integer,dimension(:),pointer    :: intarr
   integer       :: level
   integer       :: i
   real   ,dimension(:),pointer    :: realarr
   logical       :: success

   success = .true.
   level   = tree_traverse_level()
   write(*,*) ('   ', i=1,level), 'Node: ', trim(tree_get_name(tree))
   call tree_get_data_ptr( tree, data_ptr, type_string )

   select case (type_string)
   case ('STRING')
      string = '(no data)'
      call tree_get_data_string( tree, string, success )
      write(*,*) ('   ', i=1,level+1), trim(string), ' -- ', &
                 trim(type_string), ' -- ', success
   case ('INTEGER ARRAY')
      write(*,*) ('   ', i=1,level+1), transfer(tree%node_data,intarr), ' -- ', &
                 trim(type_string), ' -- ', success
   case ('REAL ARRAY')
      write(*,*) ('   ', i=1,level+1), transfer(tree%node_data,realarr), ' -- ', &
                 trim(type_string), ' -- ', success
   case default
      string = '(unknown data type)'
      write(*,*) ('   ', i=1,level+1), trim(string), ' -- ', &
                 trim(type_string), ' -- ', success
   end select
end subroutine print_tree
!
!
! --------------------------------------------------------------------
!   Subroutine: lowercase
!   Author:     Cor van der Schelde
!   Purpose:    Convert upper case characters to lower case
!   Context:    This is a copy of subroutine small in Delft3D-FLOW
!               Used inside properties module
!   Summary:
!               Scan string for upper case characters and
!               convert them.
!   Arguments:
!   string      String to be converted
!   lenstr      Length of string to be converted
! --------------------------------------------------------------------
!
subroutine lowercase(string    ,lenstr    )
    implicit none
    !
    ! Global variables
    !
    integer     , intent(in) :: lenstr
    character(*)             :: string
    !
    ! Local variables
    !
    integer :: i
    integer :: j
    integer :: newlen
    !
    !! executable statements -------------------------------------------------------
    !
    newlen = min(lenstr, len(string))
    do i = 1, newlen
       j = ichar(string(i:i))
       if ((j>64) .and. (j<91)) then
          j = j + 32
          string(i:i) = char(j)
       endif
    enddo
end subroutine lowercase

end module TREE_STRUCTURES

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
! MIT License
! Copyright (c) 2019 Sam Miller
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.

module mod_boundary_conditions
  !< Define the based boundary condition class

  use, intrinsic :: iso_fortran_env, only: ik => int32, rk => real64
  use mod_functional, only: operator(.reverse.)
  use mod_field, only: field_2d_t
  use mod_input, only: input_t
  use mod_grid_block_2d, only: grid_block_2d_t

  implicit none

  private
  public :: boundary_condition_t

  type, abstract :: boundary_condition_t
    character(len=32) :: name = ''
    character(len=2) :: location = '' !< Location (+x, -x, +y, or -y)
    real(rk) :: time = 0.0_rk ! Solution time (for time dependent bc's)
    real(rk) :: max_time = 0.0_rk !< Max time in source (e.g. stop after this)
    integer(ik) :: priority = 0 !< Certain b.c.'s take priority over others in terms of when they are applied (highest goes first)
    integer(ik) :: io_unit = 0

    integer(ik) :: n_ghost_layers = 0
    integer(ik), dimension(:), allocatable :: ilo_ghost !< (innermost:outermost ghost_layer); ghost layer indices; note: innermost (closest to real domain) layer is always first
    integer(ik), dimension(:), allocatable :: ihi_ghost !< (innermost:outermost ghost_layer); ghost layer indices; note: innermost (closest to real domain) layer is always first
    integer(ik), dimension(:), allocatable :: jlo_ghost !< (innermost:outermost ghost_layer); ghost layer indices; note: innermost (closest to real domain) layer is always first
    integer(ik), dimension(:), allocatable :: jhi_ghost !< (innermost:outermost ghost_layer); ghost layer indices; note: innermost (closest to real domain) layer is always first
    integer(ik) :: ilo = 0 !< ilo real domain cell index
    integer(ik) :: ihi = 0 !< ihi real domain cell index
    integer(ik) :: jlo = 0 !< jlo real domain cell index
    integer(ik) :: jhi = 0 !< jhi real domain cell index
  contains
    procedure, public :: set_time
    procedure, public :: get_time
    procedure, public :: set_indices
    procedure(apply), public, deferred :: apply

  endtype boundary_condition_t

  abstract interface
    subroutine apply(self, rho, u, v, p)
      import :: boundary_condition_t, field_2d_t
      import :: ik, rk
      class(boundary_condition_t), intent(inout) :: self
      class(field_2d_t), intent(inout) :: rho
      class(field_2d_t), intent(inout) :: u
      class(field_2d_t), intent(inout) :: v
      class(field_2d_t), intent(inout) :: p
    endsubroutine apply
  endinterface
contains
  pure subroutine set_time(self, time)
    class(boundary_condition_t), intent(inout) :: self
    real(rk), intent(in) :: time
    self%time = time
  endsubroutine

  pure real(rk) function get_time(self) result(time)
    class(boundary_condition_t), intent(in) :: self
    time = self%time
  endfunction

  subroutine set_indices(self, grid)
    !< Save the indices for the cells that are tagged as ghost cells. These indices will be used
    !< by the other procedures that apply the boundary conditions

    class(boundary_condition_t), intent(inout) :: self
    class(grid_block_2d_t), intent(in) :: grid
    integer(ik) :: i
    ! integer(ik), dimension(:, :), intent(in) :: ghost_layers
    !< (ilo_layers(n), ihi_layers(n), jlo_layers(n), jhi_layers(n)); indices to the ghost layers.
    !< The ilo_layers type var can be scalar or a vector, e.g. ilo_layers = [-1,0] or ilo_layers = 0

    ! Note: The indexing order starts with the inner-most index, e.g the index closest to the real domain

    ! Example:
    ! the cells are numbered as follows: [-1 0 | 1 2 3 4 | 5 6]
    ! the real cells are [1,2,3,4] and the ghost cells are [-1 0] and [5 6].
    ! Then ilo_ghost = [0, -1] and ihi_ghost = [5, 6]

    self%n_ghost_layers = grid%n_halo_cells

    allocate(self%ilo_ghost(self%n_ghost_layers))
    allocate(self%ihi_ghost(self%n_ghost_layers))
    allocate(self%jlo_ghost(self%n_ghost_layers))
    allocate(self%jhi_ghost(self%n_ghost_layers))

    if(self%n_ghost_layers == 1) then
      self%ilo_ghost = grid%lbounds_halo(1)
      self%ihi_ghost = grid%ubounds_halo(1)
      self%jlo_ghost = grid%lbounds_halo(2)
      self%jhi_ghost = grid%ubounds_halo(2)
      self%ilo = grid%lbounds(1)
      self%ihi = grid%ubounds(1)
      self%jlo = grid%lbounds(2)
      self%jhi = grid%ubounds(2)
    else if(self%n_ghost_layers > 1) then
      do i = 1, self%n_ghost_layers
        self%ilo_ghost(i) = grid%lbounds_halo(1) + i - 1
        self%ihi_ghost(i) = grid%ubounds_halo(1) - i + 1

        self%jlo_ghost(i) = grid%lbounds_halo(2) + i - 1
        self%jhi_ghost(i) = grid%ubounds_halo(2) - i + 1
      enddo
      self%ilo_ghost = .reverse.self%ilo_ghost
      self%jlo_ghost = .reverse.self%jlo_ghost
      self%ihi_ghost = .reverse.self%ihi_ghost
      self%jhi_ghost = .reverse.self%jhi_ghost

      self%ilo = grid%lbounds(1)
      self%ihi = grid%ubounds(1)
      self%jlo = grid%lbounds(2)
      self%jhi = grid%ubounds(2)
    else
      error stop "Error in boundary_condition_t%set_indicies(), n_ghost_layers <= 0"
    endif

  endsubroutine set_indices
endmodule mod_boundary_conditions

max_issues_repo_name: smillerc/cato

max_forks_repo_name: smillerc/cato


text:
Aggie Bucks are little notes you get for going to and supporting UCD Athletics UCD home games.  They may be redeemed for discounts on UC Davis apparel and merchandise sold at the Aggie Pride Outlet (APO) during select games, or possibly also a free item.  You need to swipe your reg card, or an Aggie Fanatics membership card at the entrance, and youre usually given $10 Aggie Bucks.  Theyre smaller than a dollar bill, and feature Gunrock in the middle, UCDs Mustang Mascot. More information could be found at http://www.ucdavisaggies.com UCDavisAggies.com.

Unfortunately, you can only use one coupon per item.  The free item is from the cookie jar.  The APO is located near the main entrance to most sports games, letting you easily redeem your Aggie Bucks on UCD Apparel and Merchandise. Go team, go.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
MODULE dictionary

  USE kinds

  IMPLICIT NONE

  SAVE

  !!!!!!!!!!!!!!!!!! read data variables !!!!!!!!!!!!!!!!!
  CHARACTER(len=100) :: filename
  CHARACTER(len=100) :: filename_bins
  CHARACTER(len=100) :: data_type
  REAL(REAL64), DIMENSION(:,:,:), ALLOCATABLE :: pos
  REAL(REAL64), DIMENSION(:,:,:), ALLOCATABLE :: pos_1
  REAL(REAL64), DIMENSION(:,:,:), ALLOCATABLE :: pos_2
  REAL(REAL64), DIMENSION(:,:,:), ALLOCATABLE :: pos_3
  INTEGER :: n_header ! number of lines within the header
  INTEGER :: n_particle ! number of particles
  INTEGER(INT64) :: n_lines ! number of lines the data file have
  INTEGER :: time_marks ! number of iterations
  INTEGER :: n_type_1, n_type_2, n_type_3
  REAL(REAL32), DIMENSION(2) :: XD, YD, ZD  ! size of simulation box
  INTEGER :: output_gap  ! time gap between outputs
  INTEGER :: time_marks_dens

  !!!!!!!!!!!!!!!!!!! analysis variables !!!!!!!!!!!!!!!!
  REAL(REAL64), DIMENSION(:,:), ALLOCATABLE :: dens, dens_1, dens_1_all
  REAL(REAL64) :: radius_tol ! radius cutoff for density
  INTEGER :: nbin_x, nbin_y, nbin_z
  REAL(REAL32), DIMENSION(:,:,:,:), ALLOCATABLE :: bins_1, bins_2, bins
  INTEGER, DIMENSION(:,:), ALLOCATABLE :: contact_line_index_1
  REAL(REAL32), DIMENSION(2) :: tol_1, tol_1_all
  REAL(REAL32), DIMENSION(2) :: XB, YB, ZB ! size of bins

END MODULE

max_stars_repo_name: JingBang-Liu/MovingContactLine

max_issues_repo_name: JingBang-Liu/MovingContactLine

max_forks_repo_name: JingBang-Liu/MovingContactLine


text:
!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
!> \brief Coupling Interface

module hecmw_couple_f

  use hecmw_util
  use hecmw_couple_struct_f
  use hecmw_couple_copy_c2f_f
  use hecmw_couple_copy_f2c_f

  implicit none
  private
  public :: hecmw_couple

contains

subroutine hecmw_couple(boundary_id, couple_value)

  character(len=HECMW_NAME_LEN), intent(in)    :: boundary_id
  type(hecmw_couple_value),      intent(inout) :: couple_value
  integer(kind=kint)                           :: ierr

  call hecmw_couple_exec_init_if(ierr)
  if(ierr /= 0) call hecmw_abort(hecmw_comm_get_comm())

  call hecmw_couple_copy_f2c(couple_value, ierr)
  if(ierr /= 0) call hecmw_abort(hecmw_comm_get_comm())

  call hecmw_couple_if(boundary_id, ierr)
  if(ierr /= 0) call hecmw_abort(hecmw_comm_get_comm())

  call hecmw_couple_copy_c2f(couple_value, ierr)
  if(ierr /= 0) call hecmw_abort(hecmw_comm_get_comm())

  call hecmw_couple_exec_finalize_if(ierr)
  if(ierr /= 0) call hecmw_abort(hecmw_comm_get_comm())

end subroutine hecmw_couple

end module hecmw_couple_f

max_stars_repo_name: masae-hayashi/FrontISTR

max_issues_repo_name: masae-hayashi/FrontISTR

max_forks_repo_name: masae-hayashi/FrontISTR


text:
!
!===============================================================================
! /****m*/ test_mpi
MODULE test_mpi
!
! PURPOSE:    Module contains tests for mpi_mod module
!
! CONTAINS:   Subroutine            Purpose
!             t_start_mpi()         - start MPI execution environment & get details
!             t_end_mpi()           - end MPI execution environment & stop ParaGEMS
!             t_parallel_setup()    - sets up partitioning, read mesh and nodal
!                                   locations, and sets connectivity
!             t_partition()         - partitions work across available processes
!             t_get_connectivity()  - gets process connectivity
!
! UPDATES:  created (PDB) :: 2019/08/21
!
!===============================================================================

  !-----------------------------------------------------------------------------
  ! use statements and implicit none
  !-----------------------------------------------------------------------------
  USE common_mod
  USE mpi_mod
  IMPLICIT NONE

!===============================================================================
CONTAINS
!===============================================================================
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! /****s*/ /src/mpi/modules/mpi_test|t_start_mpi
  SUBROUTINE t_start_mpi()
!
! PURPOSE:  Start MPI execution environment & get details
!
! TESTS:    Number  Purpose
!
! UPDATES:  created (PDB) :: 2019/08/21
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    IMPLICIT NONE

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! MAIN EXECUTION
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    !---------------------------------------------------------------------
    !  1. start MPI
    !---------------------------------------------------------------------
    CALL MPI_INIT(ier)
    IF (ier /= MPI_SUCCESS) THEN
      WRITE(*,'(A,A,I5)')'Error in MPI_INIT: errcode = ',ier
      STOP "ParaGEMS: shutdown"
    END IF

    !---------------------------------------------------------------------
    !  2. get local rank
    !---------------------------------------------------------------------
    CALL MPI_COMM_RANK(MPI_COMM_WORLD,rank,ier)
    IF (ier /= MPI_SUCCESS) THEN
      WRITE(*,'(A,A,I5)')'Error in MPI_COMM_RANK: errcode = ',ier
      CALL end_mpi()
    END IF

    !---------------------------------------------------------------------
    !  3. get the global number of processes
    !---------------------------------------------------------------------
    CALL MPI_COMM_SIZE(MPI_COMM_WORLD,num_procs,ier)
    IF (ier /= MPI_SUCCESS) THEN
      WRITE(*,'(A,A,I5)')'Error in MPI_COMM_SIZE: errcode = ',ier
      CALL end_mpi()
    END IF

    RETURN

  END SUBROUTINE
! mpi_test|t_start_mpi
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! /****s*/ /src/modules/mpi/mpi_test|t_end_mpi
  SUBROUTINE t_end_mpi()
!
! PURPOSE:  End MPI execution environment & stop ParaGEMS
!
! TESTS:
!
! UPDATES:  created (PDB) :: 2019/08/21
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    IMPLICIT NONE

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! MAIN EXECUTION
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    !---------------------------------------------------------------------
    !  1. end MPI execution environment and stop ParaGEMS
    !---------------------------------------------------------------------
    CALL MPI_FINALIZE(ier)
    STOP "ParaGEMS: shutdown: the program terminated successfully"

    RETURN

  END SUBROUTINE
! mpi_test|t_end_mpi
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! /****s*/ /src/modules/common/common_test|t_parallel_setup
  SUBROUTINE t_parallel_setup()
!
! PURPOSE:  sets up partitioning, read mesh and nodal locations, and sets
!           connectivity
!
! TESTS:
!
! UPDATES:  created (PDB) :: 2019/08/21
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    IMPLICIT NONE

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! MAIN EXECUTION
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    RETURN

  END SUBROUTINE
! common_test|t_parallel_setup
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! /****s*/ /src/modules/common/common_test|t_partition
  SUBROUTINE t_partition()
!
! PURPOSE:
!
! TESTS:
!
! UPDATES:  created (PDB) :: 2019/08/21
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    IMPLICIT NONE

    !---------------------------------------------------------------------
    ! local variables
    !---------------------------------------------------------------------
    CHARACTER(LEN=slen) :: fname          ! file name
    INTEGER             :: extras         ! modulus of the number of volumes
                                          ! and available processes
    INTEGER             :: int_buffer(2)  ! integer buffer for MPI comms
    INTEGER             :: junk           ! junk IO variable

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! MAIN EXECUTION
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    RETURN

  END SUBROUTINE
! common_test|t_partition
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!

END MODULE
! mpi_test
!===============================================================================
!

max_stars_repo_name: ParaFEM/ParaGEMS

max_issues_repo_name: ParaFEM/ParaGEMS

max_forks_repo_name: ParaFEM/ParaGEMS


text:
      SUBROUTINE FILLCPsuperH2(ISKIP_EDM
     .,NSMIN,NSSIN,SMPARA_H,SSPARA_H,NFLAG,IFLAG_H
     .,MCH,HMASS_H,OMIX_H
     .,STMASS_H,STMIX_H,SBMASS_H,SBMIX_H,STAUMASS_H,STAUMIX_H,SNU3MASS_H
     .,MC_H,UL_H,UR_H,MN_H,N_H,NCMAX,NHC_H,SHC_H,CHC_H
     .,NMNH,GAMBRN,NMCH,GAMBRC)
************************************************************************
*
************************************************************************
      IMPLICIT REAL*8(A-H,M,O-Z)
*
*-----------------------------------------------------------------------
*+CDE HC_ COMMON BLOCKS:
      COMMON /HC_SMPARA/ AEM_H,ASMZ_H,MZ_H,SW_H,ME_H,MMU_H,MTAU_H,MDMT_H
     .                  ,MSMT_H,MBMT_H,MUMT_H,MCMT_H,MTPOLE_H,GAMW_H
     .                  ,GAMZ_H,EEM_H,ASMT_H,CW_H,TW_H,MW_H,GW_H,GP_H
     .                  ,V_H,GF_H,MTMT_H
*
      COMMON /HC_RSUSYPARA/ TB_H,CB_H,SB_H,MQ3_H,MU3_H,MD3_H,ML3_H,ME3_H
*
      COMPLEX*16 MU_H,M1_H,M2_H,M3_H,AT_H,AB_H,ATAU_H
      COMMON /HC_CSUSYPARA/ MU_H,M1_H,M2_H,M3_H,AT_H,AB_H,ATAU_H
*
*NEW COMMON BLOCKS for V2
*
      REAL*8     RAUX_H(999)
      COMPLEX*16 CAUX_H(999)
      COMMON /HC_RAUX/ RAUX_H
      COMMON /HC_CAUX/ CAUX_H
      DATA NAUX/999/
*-----------------------------------------------------------------------
*ARRAYS:
      REAL*8 SMPARA_H(NSMIN),SSPARA_H(NSSIN)
*
      INTEGER*8 IFLAG_H(NFLAG)
*
      REAL*8     HMASS_H(3),OMIX_H(3,3)
      REAL*8     STMASS_H(2),SBMASS_H(2),STAUMASS_H(2),SNU3MASS_H
      REAL*8     MC_H(2),MN_H(4)
      COMPLEX*16 STMIX_H(2,2),SBMIX_H(2,2),STAUMIX_H(2,2)
      COMPLEX*16 UL_H(2,2),UR_H(2,2),N_H(4,4)
*
      COMPLEX*16 NHC_H(NCMAX,3)  ! 100 = NCMAX
      REAL*8     SHC_H(NCMAX)
      COMPLEX*16 CHC_H(NCMAX)
*
      REAL*8 GAMBRN(NMNH,3,3)
      REAL*8 GAMBRC(NMCH,3)
*-----------------------------------------------------------------------
*Local assays:
*
      REAL*8     HZP_WIDTH(3),HZP_WIDTH_SM(3)
*=======================================================================
      CALL FILLPARA2(NSMIN,NSSIN,SMPARA_H,SSPARA_H,NFLAG,IFLAG_H)
      MCH=SSPARA_H(2)
*
      IF(IFLAG_H(57).GT.0) GOTO 99
*=======================================================================
      CALL FILLHIGGS2(NSMIN,NSSIN,SMPARA_H,SSPARA_H,NFLAG,IFLAG_H,MCH
     .               ,HMASS_H,OMIX_H)
*      print*,'>0 ',(hmass_h(i),i=1,3)
*
      CALL INCL_STAU(NFLAG,IFLAG_H,HMASS_H,OMIX_H)
*      print*,'>1 ',(hmass_h(i),i=1,3)
*
      MCH=RAUX_H(10) ! Charged Higgs pole mass or effetive-pot. mass
*
      IERR1=IFLAG_H(50)+IFLAG_H(51)+IFLAG_H(52)+IFLAG_H(54)
     .     +IFLAG_H(56)
     .     +IFLAG_H(55)+IFLAG_H(60)
      IF(IERR1.GT.0) GOTO 99
*
      IF(IFLAG_H(53).EQ.1) THEN
       print*,'WARNING! IFLAG_H(53) = ',IFLAG_H(53)
       IFLAG_H(53)=0
      ENDIF
*=======================================================================
      CALL FILLCOUPL2(NSMIN,NSSIN,SMPARA_H,SSPARA_H,NFLAG,IFLAG_H,MCH
     . ,HMASS_H,OMIX_H,STMASS_H,STMIX_H,SBMASS_H,SBMIX_H,STAUMASS_H
     . ,STAUMIX_H,SNU3MASS_H,MC_H,UL_H,UR_H,MN_H,N_H,NCMAX,NHC_H,SHC_H
     . ,CHC_H)
*
      IERR2=IFLAG_H(54)+IFLAG_H(56)
      IF(IERR2.GT.0) GOTO 99
*=======================================================================
      CALL FILLGAMBR2(NSMIN,NSSIN,SMPARA_H,SSPARA_H,NFLAG,IFLAG_H
     . ,MCH,HMASS_H,NCMAX,NHC_H,SHC_H,CHC_H,STMASS_H
     . ,STMIX_H,SBMASS_H,SBMIX_H,STAUMASS_H,STAUMIX_H,SNU3MASS_H
     . ,MC_H,MN_H,NMNH,GAMBRN,NMCH,GAMBRC)
*
      CALL H2ZGAMMA_V0(ISKIP_EDM
     .,NSMIN,NSSIN,SMPARA_H,SSPARA_H,NFLAG,IFLAG_H
     .,MCH,HMASS_H,OMIX_H
     .,STMASS_H,STMIX_H,SBMASS_H,SBMIX_H,STAUMASS_H,STAUMIX_H,SNU3MASS_H
     .,MC_H,UL_H,UR_H,MN_H,N_H,NCMAX,NHC_H,SHC_H,CHC_H
     .,NMNH,GAMBRN,NMCH,GAMBRC,HZP_WIDTH,HZP_WIDTH_SM)
*
      CALL UPDATEGAMBR_0(ISKIP_EDM
     .,NSMIN,NSSIN,SMPARA_H,SSPARA_H,NFLAG,IFLAG_H
     .,MCH,HMASS_H,OMIX_H
     .,STMASS_H,STMIX_H,SBMASS_H,SBMIX_H,STAUMASS_H,STAUMIX_H,SNU3MASS_H
     .,MC_H,UL_H,UR_H,MN_H,N_H,NCMAX,NHC_H,SHC_H,CHC_H
     .,NMNH,GAMBRN,NMCH,GAMBRC,HZP_WIDTH,HZP_WIDTH_SM)
*
      MHSM=SMPARA_H(20)
      CALL FILLSMBRS(MHSM,GAMBRN,NMNH)
*=======================================================================
*For the subroutine FILLDHPG, we need \sqrt{s} value as an input.
      SQRTS=RAUX_H(101)
*
      CALL FILLDHPG(SQRTS,NFLAG,IFLAG_H
     . ,MCH,HMASS_H,OMIX_H,NCMAX,NHC_H,SHC_H,CHC_H
     . ,STMASS_H,STMIX_H,SBMASS_H,SBMIX_H,STAUMASS_H,STAUMIX_H
     . ,SNU3MASS_H,MC_H,UL_H,UR_H,MN_H,N_H,NMNH,GAMBRN,NMCH,GAMBRC)
*=======================================================================
      CALL FILLBOBS(NSMIN,NSSIN,SMPARA_H,SSPARA_H,NFLAG,IFLAG_H
     .             ,HMASS_H,OMIX_H,MCH,MC_H,UL_H,UR_H,STMASS_H,STMIX_H)
*=======================================================================
      IF(ISKIP_EDM.EQ.0) THEN
      CALL FILLEDMS(NSMIN,NSSIN,SMPARA_H,SSPARA_H,NFLAG,IFLAG_H
     . ,MCH,HMASS_H,OMIX_H
     . ,STMASS_H,STMIX_H,SBMASS_H,SBMIX_H,STAUMASS_H
     . ,STAUMIX_H,SNU3MASS_H
     . ,MC_H,UL_H,UR_H,MN_H,N_H,NCMAX,NHC_H,SHC_H,CHC_H)

      CALL HGRA_EDMS(NSMIN,NSSIN,SMPARA_H,SSPARA_H,NFLAG,IFLAG_H
     . ,MCH,HMASS_H,OMIX_H
     . ,STMASS_H,STMIX_H,SBMASS_H,SBMIX_H,STAUMASS_H
     . ,STAUMIX_H,SNU3MASS_H
     . ,MC_H,UL_H,UR_H,MN_H,N_H,NCMAX,NHC_H,SHC_H,CHC_H)
      ENDIF ! ISKIP_EDM
*=======================================================================
      IF(ISKIP_EDM.EQ.0) THEN
      CALL FILLMUON(NSMIN,NSSIN,SMPARA_H,SSPARA_H,NFLAG,IFLAG_H
     . ,MCH,HMASS_H,OMIX_H
     . ,STMASS_H,STMIX_H,SBMASS_H,SBMIX_H,STAUMASS_H
     . ,STAUMIX_H,SNU3MASS_H
     . ,MC_H,UL_H,UR_H,MN_H,N_H,NCMAX,NHC_H,SHC_H,CHC_H)
      ENDIF ! ISKIP_EDM
*=======================================================================
      CALL FILLCOLL(NSMIN,NSSIN,SMPARA_H,SSPARA_H,NFLAG,IFLAG_H
     .,MCH,HMASS_H,OMIX_H
     .,STMASS_H,STMIX_H,SBMASS_H,SBMIX_H,STAUMASS_H,STAUMIX_H,SNU3MASS_H
     .,MC_H,UL_H,UR_H,MN_H,N_H,NCMAX,NHC_H,SHC_H,CHC_H
     .,NMNH,GAMBRN,NMCH,GAMBRC)
*=======================================================================
      IF(IFLAG_H(30).EQ.1) THEN
      CALL FILLSLHA2(
     . NSMIN,NSSIN,SMPARA_H,SSPARA_H,NFLAG,IFLAG_H
     .,MCH,HMASS_H,OMIX_H
     .,STMASS_H,STMIX_H,SBMASS_H,SBMIX_H,STAUMASS_H,STAUMIX_H,SNU3MASS_H
     .,MC_H,UL_H,UR_H,MN_H,N_H,NCMAX,NHC_H,SHC_H,CHC_H
     .,NMNH,GAMBRN,NMCH,GAMBRC)
      ENDIF
*=======================================================================
 99   CONTINUE
*
      RETURN
      END


subroutine solver()
!determines the flow profile around the object. In order to do this we have to solve 2
!coupled equations:
! - calculate the streamfunction (psi) by solving Poisson's equation
! - update the vorticity by solving the vorticity transport equation (from the NS equations)
!
!First of all a potential flow is calculated around the object (solving Laplace's equation)
!then (optionally) the vorticity is turned on and alternately Poisson's equation and the
!vorticity transport equation are iterated to determine the fluid flow.

    use vars
    use parallel

    implicit none

    real :: time
    double precision :: tstart, tstop

    !$OMP PARALLEL
    !solve Laplace's equation for the irrotational flow profile (vorticity=0)
    call poisson(5000)
    !$OMP END PARALLEL

    !get the vorticity on the surface of the object

    call getvort()

    !write out the potential flow to file.

    call writetofile("potential.dat")

    !if the user has chosen to allow vorticity:
    if (vorticity) then

        ! set CFL conditions (dt << cell crossing time or cell diffusion time for stability)
        cfl_r0 = 0.25*dx*dy*R0
        cfl_v = 0.25*dx/1. !assume velocty is 1

        !dt is set according to the most restrictive CFL condition
        dt = minval((/ cfl_r0, cfl_v /))

        if (irank .eq. 0) then
            print*,''
            print*,"CFL (Reynolds, velocity)", cfl_r0, cfl_v
            print *, "Timestep=",dt
            print *, ""
        endif

        time=0.

        if (irank .eq. 0) tstart=MPI_Wtime()

        !$OMP PARALLEL

        do while (time .lt. real(nx)*crossing_times)

            !if (irank .eq. 0) print*, "t=",time,"of",nx*crossing_times

            call getv() !get the velocity

            call navier_stokes() !solve dw/dt=f(w,psi) for a timestep

            call poisson(2) !2 poisson 'relaxation' steps

            !$OMP SINGLE
            time=time+dt
            !$OMP END SINGLE
            !$OMP BARRIER

        enddo

        !$OMP END PARALLEL

        if (irank .eq. 0) then
            tstop=MPI_Wtime()
            print*,''
            print*, "Time to complete =",tstop-tstart
            print*,''
        endif

        call getvort()

    endif

end subroutine


text:
#
# CDDL HEADER START
#
# The contents of this file are subject to the terms of the
# Common Development and Distribution License (the "License").
# You may not use this file except in compliance with the License.
#
# You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
# or http://www.opensolaris.org/os/licensing.
# See the License for the specific language governing permissions
# and limitations under the License.
#
# When distributing Covered Code, include this CDDL HEADER in each
# file and include the License file at usr/src/OPENSOLARIS.LICENSE.
# If applicable, add the following below this CDDL HEADER, with the
# fields enclosed by brackets "[]" replaced with your own identifying
# information: Portions Copyright [yyyy] [name of copyright owner]
#
# CDDL HEADER END
#
#
# Copyright 2007 Sun Microsystems, Inc.  All rights reserved.
# Use is subject to license terms.
#

set $dir=/mlfs
#set $nfiles=5
#set $nfiles=40
#set $nfiles=200
#set $nfiles=500
#set $nfiles=1000
set $nfiles=1000
#set $meandirwidth=10000
set $meandirwidth=0
#set $filesize=cvar(type=cvar-gamma,parameters=mean:16384;gamma:1.5)
set $filesize=4k
set $nthreads=1
set $iosize=1m
set $meanappendsize=4k
set $iters=1
set $eventrate=500

define fileset name=bigfileset,path=$dir,size=$filesize,entries=$nfiles,dirwidth=0,prealloc=80

define process name=multiclient1,devid=4,instances=1
{
  thread name=filereaderthread,memsize=10m,instances=$nthreads
  {
    flowop bwlimit name=bwlimit4
    flowop deletefile name=deletefile41,filesetname=bigfileset,fd=2

    flowop createfile name=createfile41,filesetname=bigfileset,fd=1
    flowop appendfile name=appendfile41,iosize=$meanappendsize,iters=$iters,fd=1
    flowop closefile name=closefile41,filesetname=bigfileset,fd=1
    flowop deletefile name=deletefile42,filesetname=bigfileset,fd=1

    flowop createfile name=createfile42,filesetname=bigfileset,fd=2
    flowop appendfile name=appendfile42,iosize=$meanappendsize,iters=$iters,fd=2

    flowop fsync name=rsync4,filesetname=bigfileset,fd=2

    flowop closefile name=closefile42,filesetname=bigfileset,fd=2

  }
}

#debug 5

eventgen rate = $eventrate

#enable lathist

echo  "Single-process filemicro personality successfully loaded"

run 30

max_stars_repo_name: hayley-leblanc/SplitFS

max_issues_repo_name: hayley-leblanc/SplitFS

max_forks_repo_name: hayley-leblanc/SplitFS


text:
Subroutine finiteStrainDelta(AE, BE, MaxLRows, XLL, NDim, iDofT, NodElt, Sol0, Sol1, CommonPar, Param, JParam, DelT, DTm, Time)
!     ------------------------------------------------------------------
	use funcAux
	use loadingLib
	use finiteStrainLib

	implicit none

	!	PARAMETERS
	integer :: iSimplex, NGPMAX, iBu
	parameter (iBu = 0, NGPMAX = 11, iSimplex = 0) ! No Bubbles and Linear, ! iSimplex = 0 -> The element is a Simplex
!~ 	parameter (iBu = 0, pOrder = 2,NGP = 11, iSimplex = 0) ! No Bubbles and Linear, ! iSimplex = 0 -> The element is a Simplex

	!   ===== SUBROUTINE ARGUMENTS  =======
	integer :: MaxLRows,Ndim,iDofT, NodElt! all integers
	Real*8 :: DelT, DTm,Time ! all reals
	Real*8 :: AE(MaxLRows,MaxLRows), BE(MaxLRows), XLL(*), Sol0(*), Sol1(*), Param(*), JParam(*), CommonPar(*)

	!   =====   END ARGUMENTS  =======

	Integer :: nG, i,j,k,l,m,n,p,q, A, B, ApRow, BpCol , App,Bqp, ipDim, ipDimE, ipDof,ipDofmod ! counters
	Integer :: NparamElem, constLaw
	logical :: ifAnalytical
	integer , parameter ::  NdimE = 2 , NodG = 3 ! this code only makes sense in 2D and for triangles
	Real*8 :: PhiA,PhiB, dv , Det , bm(NdimE), MatPar(9)
	Real*8 :: Psi(NdimE,NGPMAX), W(NGPMAX), Phi(NodG), dPhiA(NdimE), dPhiB(NdimE),&
			  dPhi_G(NdimE,NodG), dPhi_L(NdimE,NodG), Jac(NDimE,NDimE)
	Real*8 :: GradU(NdimE,NdimE) , GradDeltaU(NdimE,NdimE), F(NdimE,NdimE)
	Real*8 :: D(NdimE,NdimE,NdimE,NdimE), DB(NdimE,NdimE,NdimE,NdimE), Ppk(NdimE,NdimE)
	Real*8 :: SolU(NodG*NdimE) , SolDeltaU(NodG*NdimE), XLLMod(NodG*NdimE), DBGradDeltaU(NdimE,NdimE)
	integer :: pOrder , NGP
	Real*8 :: LoadPar(6) , energy
	integer :: iShiftU, iShiftDeltaU , LoadType , LoadTypeProg

	iShiftU = nint(CommonPar(1))
	iShiftDeltaU = nint(CommonPar(2))
	constLaw = nint(CommonPar(3))
	matPar(1:3) = commonPar(4:6)
	pOrder = nint(commonPar(7))
	LoadType = nint(commonPar(8))
	LoadTypeProg = nint(commonPar(9))
	LoadPar(:) = CommonPar(10:15)

	call getSlice(SolU,Sol1,1,NodG,iShiftU + 1 ,iShiftU + NdimE, iDofT)
	call getSlice(SolDeltaU,Sol1,1,NodG,iShiftDeltaU + 1 ,iShiftDeltaU + NdimE, iDofT)
	call getSlice(XLLmod,XLL,1,NodG,1 ,NdimE, Ndim)

	select case(pOrder)
		case(1)
			NGP = 1
		case(2)
			NGP = 3
		case default
			write(0,*) "pOrder must be 1 or 2"
	end select

	W = 0.0d0
	Call GaussRule (Psi,W,NdimE,NGP,iSimplex)

	Do nG = 1, NGP ! LoopGauss

		Call LocalShapeDer(NodG,NdimE,dPhi_L,PSI(1,nG),pOrder,iBu) ! PSI(1,nG) is a pointer to the nG^th column
		if ( nG .eq. 1 ) then
			Call Jacobian(Jac,Det,NdimE,NodG,XLLmod,dPhi_L)
		endif
		Call GlobalShapeDer(NodG,NdimE,dPhi_L,dPhi_G,Jac)
		Call ShapeF (NodG,NdimE,Phi,PSI(1,nG),pOrder,iBu)

		call calcGradU(GradU,SolU,dPhi_G,NdimE,NodG,NdimE)
		call calcGradU(GradDeltaU,SolDeltaU,dPhi_G,NdimE,NodG,NdimE)

		call setF(F,LoadPar,Time,DelT,LoadTypeProg,LoadType)

		F = F + deltaKron(1:NdimE,1:NdimE) + GradU

		call calcPpk(Ppk,F,NdimE,matPar,constLaw)
		call calcD(D,F,NdimE,matPar,constLaw)

		energy = 0.0d0

!~ 		DB = 0.0d0

		DV=Det*W(nG)

!~ 		call T4xT2(DBGradDeltaU,DB,GradDeltaU)

		Do A=1, NodG !LoopRow
			ApRow  = (A-1)*iDofT + iShiftDeltaU
			PhiA = Phi(A)
			dPhiA=dPhi_G(:,A)

			Do p = 1 , NdimE
				App = ApRow+p

				BE(App) = BE(App) - dot_product(Ppk(p,:),dPhiA)*DV
!~ 				BE(App) = BE(App) + dot_product(DBGradDeltaU(p,:),dphiA)*DV

				do B= 1 , NodG ! LoopCol ! not considering the simmetry
					BpCol  = (B-1)*iDofT + iShiftDeltaU
					PhiB   = Phi(B)
					dPhiB=dPhi_G(:,B)

					do q=1,NdimE
						Bqp = BpCol+q

						AE(App,Bqp) = AE(App,Bqp) + DGradPhiBdotGradPhiA(D,dPhiB,dPhiA,p,q,NdimE)*DV

					end do ! loop Bq
				Enddo !LoopCol
			end do ! loop Ap
		Enddo !LoopRow
	EndDo !LoopGauss

 	Do A=1, NodG !LoopRow
		ApRow  = (A-1)*iDofT + iShiftU

		Do p = 1 , NdimE
			App = ApRow+p

			BE(App) = SolU((A-1)*NdimE + p)

			AE(App,App) = 1.0d0
			AE(App,App + iShiftDeltaU - iShiftU) = -1.0d0

		end do ! loop Ap
	Enddo !LoopRow

end subroutine

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Subroutine finiteStrainDeltaS(Coupling,CommonPar,iDofT,Ndim,MaxLRows,iAdd)
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    use funcAux
    IMPLICIT NONE
    Integer :: MaxLRows,Ndim,iDofT, iAdd
    Integer   Coupling(MaxLRows,MaxLRows)
    Real*8 CommonPar(*)

    Integer :: A,B,q,p, ApRow, BpCol,iShiftU , iShiftDeltaU
    integer , parameter :: NodG = 3 , NdimE = 2

    Coupling = 0

    iShiftU = nint(commonPar(1))
    iShiftDeltaU = nint(commonPar(2))

    do A = 1,NodG
        ApRow = (A-1) * iDofT + iShiftDeltaU

		do B = 1,NodG
			BpCol = (B-1) * iDofT + iShiftDeltaU

			do p=1,nDimE
			do q=1,nDimE

				Coupling (ApRow + p, BpCol + q) = 1

			enddo
            enddo
		enddo
	enddo

    do A = 1,NodG
        ApRow = (A-1) * iDofT + iShiftU
		do p=1,nDimE
			Coupling (ApRow + p, ApRow + p) = 1
			Coupling (ApRow + p, ApRow + p + iShiftDeltaU - iShiftU) = 1
		enddo
	enddo

	call numPrint(Coupling)

end Subroutine

max_stars_repo_name: felipefr/Piola

max_issues_repo_name: felipefr/Piola

max_forks_repo_name: felipefr/Piola


text:
C    Copyright(C) 2008-2017 National Technology & Engineering Solutions of
C    Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C    NTESS, the U.S. Government retains certain rights in this software.
C
C    Redistribution and use in source and binary forms, with or without
C    modification, are permitted provided that the following conditions are
C    met:
C
C    * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C    * Redistributions in binary form must reproduce the above
C      copyright notice, this list of conditions and the following
C      disclaimer in the documentation and/or other materials provided
C      with the distribution.
C
C    * Neither the name of NTESS nor the names of its
C      contributors may be used to endorse or promote products derived
C      from this software without specific prior written permission.
C
C    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C    OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
C
C=======================================================================
      SUBROUTINE COMAND (A, INLINE, INTYP, CFIELD, IFIELD, RFIELD,
     &   NAMECO, BLKTYP, NAMES, TIMES, IPTIMS,
     &   IDELB, VISELB, SELELB,
     &   QAREC, INFREC, RETVRB, MERR)
C=======================================================================
C   --*** COMAND *** (ALGEBRA) Process command
C   --   Written by Amy Gilkey - revised 05/18/88
C   --
C   --COMAND processes an ALGEBRA command.  The commands are:
C   --
C   --   TITLE                     Change database title
C   --
C   --   SAVE     {N,E,G,A,var,...} Save input database variable
C   --   DELETE   {var,var,...}    Delete assigned variable
C   --   ALIAS    {name,var,n}     Assign alias to vector/tensor variables
C   --
C   --   TMIN     {TMIN}           Set minimum selected time
C   --   TMAX     {TMAX}           Set maximum selected time
C   --   DELTIME  {DELT}           Set selected time interval
C   --   NINTV    {NINTV}          Change selected time interval
C   --   ZINTV    {NINTV}          Change selected time interval(zero interval)
C   --   ALLTIMES                  Select all times
C   --   TIMES    {t1,t2...}       Select specified times
C   --   STEPS    {n1,n2...}       Select specified steps
C   --
C   --   ZOOM     {x1,x2,y1,y2,z1,z2} Set zoomed mesh limits
C   --   VISIBLE  {id1,...}        Set element blocks to be written
C   --   BLOCKS   {id1,...}        Set selected element blocks
C   --   MATERIAL {id1,...}        Set selected element blocks
C   --
C   --   LIST     {option}         Display database information
C   --   SHOW     {option}         Display parameter setting
C   --   HELP     {option}         Help on ALGEBRA
C   --
C   --   END                       Exit to evaluate equations (also EXIT)
C   --   QUIT                      Exit but do not write database
C   --
C   --Parameters:
C   --   A         - IN - the dynamic memory base array
C   --   INLINE    - IN/OUT - the parsed input lines for the log file
C   --   INTYP     - IN - the input types from the free field reader
C   --   CFIELD    - IN - the character fields
C   --   IFIELD    - IN - the integer fields
C   --   RFIELD    - IN - the real fields
C   --   NAMECO    - IN - the coordinate names
C   --   BLKTYP    - IN - the element block names
C   --   NAMES     - IN - the global, nodal, and element variable names
C   --   TIMES     - IN - the database time steps
C   --   IPTIMS    - IN/OUT - the selected times steps
C   --   IDELB     - IN - the element block IDs
C   --   VISELB(i) - IN/OUT - true iff element block i is to be written
C   --   SELELB(i) - IN/OUT - true iff element block i is selected
C   --   QAREC     - IN - the QA records containing:
C   --              (1) - the analysis code name
C   --              (2) - the analysis code QA descriptor
C   --              (3) - the analysis date
C   --              (4) - the analysis time
C   --   INFREC    - IN - the information records
C   --   RETVRB    - OUT - returned action verb (END, QUIT, PRINT, LOG)
C   --
C   --Common Variables:
C   --   Uses NDIM, NSTEPS of /DBNUMS/
C   --   Sets TITLEO of /DBTITL/
C   --   Sets NPTIMS, TMIN, TMAX, DELT, NINTV, WHONLY of /TIMES/
C   --   Sets ISZOOM, ZMLIM of /ZOOM/
C   --   Uses FNCNAM of /FNCTBL/

      include 'exodusII.inc'
      include 'ag_namlen.blk'
      include 'ag_dbnums.blk'
      include 'ag_dbtitl.blk'
      include 'ag_times.blk'
      include 'ag_zoom.blk'
      include 'ag_filter.blk'
      include 'ag_remove.blk'
      include 'ag_fnctbc.blk'

      common /debugc/ cdebug
      common /debugn/ idebug
      character*(mxstln) cdebug

      LOGICAL FFEXST, FFNUMB, MATSTR

      DIMENSION A(*)
      CHARACTER*(*) INLINE(*)
      INTEGER       INTYP(*)
      CHARACTER*(*) CFIELD(*)
      INTEGER       IFIELD(*)
      REAL          RFIELD(*)
      CHARACTER*(namlen) NAMECO(*)
      CHARACTER*(MXSTLN) BLKTYP(*)
      CHARACTER*(namlen) NAMES(*)
      REAL TIMES(*)
      INTEGER IPTIMS(*)
      INTEGER IDELB(*)
      LOGICAL VISELB(NELBLK)
      LOGICAL SELELB(NELBLK)
      CHARACTER*(MXSTLN) QAREC(4,*)
      CHARACTER*(MXLNLN) INFREC(*)
      CHARACTER*(*) RETVRB
      CHARACTER*1 OPTION

      CHARACTER*(maxnam) VERB, WORD
      CHARACTER*(mxstln) HLPTYP
      LOGICAL ISON
      LOGICAL HELP
      LOGICAL XLIM, YLIM, ZLIM
      INTEGER MERR

      CHARACTER*(mxstln) CMDTBL(26)
      SAVE CMDTBL
C      --CMDTBL - the commands table

      CHARACTER*(mxstln) HLPTBL(4)
      SAVE HLPTBL
C      --HLPTBL - the HELP type table

      MERR = 0

C   --Command table follows.  Remember to change the dimensioned size when
C   --changing the table.
      DATA CMDTBL /
     1  'TITLE                           ',
     2  'SAVE                            ',
     3  'DELETE                          ',
     4  'ALIAS                           ',
     5  'TMIN                            ',
     6  'TMAX                            ',
     7  'DELTIME                         ',
     8  'NINTV                           ',
     9  'ZINTV                           ',
     *  'ALLTIMES                        ',
     1  'TIMES                           ',
     2  'STEPS                           ',
     3  'ZOOM                            ',
     4  'FILTER                          ',
     5  'VISIBLE                         ',
     6  'BLOCKS                          ',
     7  'MATERIAL                        ',
     8  'LOG                             ',
     9  'LIST                            ',
     *  'SHOW                            ',
     1  'HELP                            ',
     2  'END                             ',
     3  'EXIT                            ',
     4  'QUIT                            ',
     4  'REMOVE                          ',
     *  '                                ' /

C   --HELP type table follows.  Remember to change the dimensioned size when
C   --changing the table.
      DATA HLPTBL /
     &  'RULES                           ',
     *  'COMMANDS                        ',
     *  'FUNCTION                        ',
     &  '                                ' /

      RETVRB = ' '

      IFLD = 1
      CALL FFCHAR (IFLD, INTYP, CFIELD, ' ', WORD)
      CALL ABRSTR (VERB, WORD, CMDTBL)
      IF (VERB .EQ. ' ') VERB = WORD

      IF (VERB .EQ. 'TITLE') THEN
         CALL FFADDC (VERB, INLINE(1))
         CALL GETINP (0, 0, 'TITLE> ', TITLEO, IOSTAT)
         INLINE(2) = TITLEO

      ELSE IF (VERB .EQ. 'ALIAS') THEN
         CALL FFADDC (VERB, INLINE(1))
         CALL ALICMD (INLINE, INTYP(IFLD), CFIELD(IFLD), IFIELD(IFLD),
     &      NAMES, *150)

      ELSE IF (VERB .EQ. 'SAVE') THEN
         CALL FFADDC (VERB, INLINE(1))
         CALL SAVCMD (INLINE, INTYP(IFLD), CFIELD(IFLD), NAMES, *150)

      ELSE IF (VERB .EQ. 'DELETE') THEN
         CALL FFADDC (VERB, INLINE(1))
         CALL DELCMD (INLINE, INTYP(IFLD), CFIELD(IFLD), *150)

      ELSE IF ((VERB .EQ. 'TMIN') .OR. (VERB .EQ. 'TMAX')
     &   .OR. (VERB .EQ. 'DELTIME')
     &   .OR. (VERB .EQ. 'NINTV') .OR. (VERB .EQ. 'ZINTV')
     &   .OR. (VERB .EQ. 'ALLTIMES')
     &   .OR. (VERB .EQ. 'TIMES') .OR. (VERB .EQ. 'STEPS')) THEN
         CALL CMDTIM (INLINE, VERB, IFLD, INTYP, CFIELD, IFIELD,
     &                RFIELD, NSTEPS, TIMES, TMIN, TMAX, DELT,
     &                NINTV, NPTIMS, IPTIMS)
      ELSE IF (VERB .EQ. 'FILTER') THEN
         CALL FFADDC (VERB, INLINE(1))
C ... Get name of element variable to use for filtering
         CALL FFCHAR (IFLD, INTYP, CFIELD, ' ', WORD)
         if (.not. matstr(word, 'ELEMENTS', 1)) THEN
           CALL PRTERR ('CMDERR',
     &       'Expected "ELEMENTS".  Syntax: FILTER ELEMENTS {elem_var}')
           VERB = ' '
           GOTO 150
         ENDIF
         CALL FFADDC (word, INLINE(1))
         CALL FFCHAR (IFLD, INTYP, CFIELD, ' ', WORD)
         CALL DBVIX('E', 1, IEV)
         IDXFLT = LOCSTR (WORD, NVAREL, NAMES(IEV))
         if (idxflt .eq. 0) then
           CALL PRTERR ('CMDERR',
     &       'Variable name specified does not exist on database.')
           VERB = ' '
           GOTO 150
         end if
         CALL FFADDC (NAMES(idxflt+iev-1), INLINE(1))

C ... Check for comparison -- uses fortran LT, LE, EQ, NE, GT, GE
         CALL FFCHAR (IFLD, INTYP, CFIELD, ' ', WORD)
         if (.not. matstr(word, 'LT', 2) .and.
     *       .not. matstr(word, 'LE', 2) .and.
     *       .not. matstr(word, 'EQ', 2) .and.
     *       .not. matstr(word, 'NE', 2) .and.
     *       .not. matstr(word, 'GT', 2) .and.
     *       .not. matstr(word, 'GE', 2)) then
           CALL PRTERR ('CMDERR',
     &       'Expected "LT", "LE", "EQ", "NE", "GT", or "GE".')
           VERB = ' '
           GOTO 150
         end if
         CALL FFADDC (word, INLINE(1))
         if (matstr(word, 'LT', 2)) cmpflt = 1
         if (matstr(word, 'LE', 2)) cmpflt = 2
         if (matstr(word, 'EQ', 2)) cmpflt = 3
         if (matstr(word, 'NE', 2)) cmpflt = 4
         if (matstr(word, 'GT', 2)) cmpflt = 5
         if (matstr(word, 'GE', 2)) cmpflt = 6

         CALL FFREAL (IFLD, INTYP, RFIELD,
     &     'filter value', 0.0, VALFLT, *150)
         CALL FFADDR (VALFLT, INLINE(1))

         CALL FFCHAR (IFLD, INTYP, CFIELD, ' ', WORD)
         if (.not. matstr(word, 'TIME', 1)) THEN
           CALL PRTERR ('CMDERR',
     &       'Expected "TIME". Syntax:')
           CALL PRTERR ('CMDSPEC',
     *       'FILTER ELEMENTS {elem_var} {CMP} {value} TIME {time}')
           VERB = ' '
           GOTO 150
         ENDIF
         CALL FFADDC (word, INLINE(1))

         CALL FFREAL (IFLD, INTYP, RFIELD,
     &     'database time', 0.0, TIMFLT, *150)
C ... Check that specified time is in range of times on database. Assume TIMES() is sorted.
         if (timflt .lt. times(1)) then
           call prterr('WARNING',
     *       'Specified time is less than minimum database time.')
           call prterr('CMDSPEC',
     *       'Setting filter time to minimum database time.')
           timflt = times(1)
         else if (timflt .gt. times(nsteps)) then
           call prterr('WARNING',
     *       'Specified time is greater than maximum database time.')
           call prterr('CMDSPEC',
     *       'Setting filter time to maximum database time.')
           timflt = times(nsteps)
         end if
         CALL FFADDR (TIMFLT, INLINE(1))
         ISFILTER = .TRUE.

      ELSE IF (VERB .EQ. 'REMOVE') THEN
         CALL FFADDC (VERB, INLINE(1))
         CALL FFCHAR (IFLD, INTYP, CFIELD, ' ', WORD)
         if (.not. matstr(word, 'ELEMENTS', 1)) THEN
           CALL PRTERR ('CMDERR',
     &       'Expected "ELEMENTS".  Syntax: FILTER ELEMENTS {elem_var}')
           VERB = ' '
           GOTO 150
         ENDIF
         CALL FFADDC (word, INLINE(1))
C ... See if LOCAL or GLOBAL or nothing specified for id space...
         IF (INTYP(IFLD) .EQ. 0) THEN
           call ffchar(IFLD, INTYP, CFIELD, ' ', WORD)
           if (matstr(word, 'GLOBAL', 1)) then
             idsglobal = .true.
           else if (matstr(word, 'GLOBAL', 1)) then
             idsglobal = .false.
           end if
         else
C ... Default is local ids or last set value.
           if (irmcnt .eq. 0) then
             idsglobal = .false.
           end if
         END IF
C ... Gather ids of elements to delete...
 99      continue
         IF (FFEXST (IFLD, INTYP)) THEN
           CALL FFINTG (IFLD, INTYP, IFIELD,
     &         'element ID', 0, ID, *129)
           IRMCNT = IRMCNT + 1
           if (IRMCNT .gt. 1024) then
             CALL PRTERR ('CMDERR',
     &         'Can only specify a maximum of 1024 elements to delete')
             goto 129
           end if
           idsrem(irmcnt) = id
           goto 99
         END IF
 129     continue
         ISREMOVE = .TRUE.

      ELSE IF (VERB .EQ. 'ZOOM') THEN
         CALL FFADDC (VERB, INLINE(1))
         CALL FFNEED (IFLD, INTYP, 'R', MIN(3,NDIM)*2,
     &      'zoom mesh limits for all dimensions', *150)
         DO 100 I = 1, MIN(3,NDIM)*2
            CALL FFREAL (IFLD, INTYP, RFIELD,
     &         'zoom mesh limit', 0.0, ZMLIM(I), *150)
            CALL FFADDR (ZMLIM(I), INLINE(1))
  100    CONTINUE
         CALL FFCHAR (IFLD, INTYP, CFIELD, 'INSIDE', WORD)
         IF (MATSTR (WORD, 'INSIDE', 1)) THEN
           CALL FFADDC ('INSIDE', INLINE(1))
           ZOOMIN = .TRUE.
         ELSE IF (MATSTR (WORD, 'OUTSIDE', 1)) THEN
           CALL FFADDC ('OUTSIDE', INLINE(1))
           ZOOMIN = .FALSE.
         END IF
C        Set logical variables which determine if zoom mesh limits
C        are of the correct format (minlimit < maxlimit)
         XLIM = (ZMLIM(1) .LE. ZMLIM(2))
         YLIM = (ZMLIM(3) .LE. ZMLIM(4))
         IF (NDIM .EQ. 3) THEN
            ZLIM = (ZMLIM(5) .LE. ZMLIM(6))
         ELSE
            ZLIM = .TRUE.
         ENDIF
         IF ((.NOT. XLIM) .OR. (.NOT. YLIM) .OR.
     &       (.NOT. ZLIM)) THEN
C           Zoom limits are incorrect
            CALL PRTERR ('CMDERR',
     &      'Zoom limits are incorrect. (minlimit > maxlimit)')
            INLINE(1) = ' '
         ELSE
            ISZOOM = .TRUE.
         END IF
      ELSE IF (VERB .EQ. 'VISIBLE') THEN
         CALL FFADDC (VERB, INLINE(1))

         IF (.NOT. FFEXST (IFLD, INTYP)) THEN

C         --Select all element blocks if no parameters

            CALL INILOG (NELBLK, .TRUE., VISELB)
            OPTION = '+'

         ELSE IF (.NOT. FFNUMB (IFLD, INTYP)) THEN

C         --Strip off ADD or DELETE option
            CALL FFCHAR (IFLD, INTYP, CFIELD, ' ', WORD)
            IF (MATSTR (WORD, 'ADD', 1)) THEN
               CALL FFADDC ('ADD', INLINE(1))
               OPTION = '+'
            ELSE IF (MATSTR (WORD, 'DELETE', 1)) THEN
               CALL FFADDC ('DELETE', INLINE(1))
               OPTION = '-'
            ELSE
               CALL PRTERR ('CMDERR',
     &            'Expected "ADD" or "DELETE" or element block ID')
               GOTO 160
            END IF

         ELSE

C         --De-select all element blocks so only listed blocks are selected
            CALL INILOG (NELBLK, .FALSE., VISELB)
            OPTION = '+'
         END IF

  110    CONTINUE
         IF (FFEXST (IFLD, INTYP)) THEN
            CALL FFINTG (IFLD, INTYP, IFIELD,
     &         'element block ID', 0, ID, *120)
            IELB = LOCINT (ID, NELBLK, IDELB)
            IF (IELB .GT. 0) THEN
               CALL FFADDI (ID, INLINE(1))
               VISELB(IELB) = (OPTION .EQ. '+')
            ELSE
               CALL INTSTR (1, 0, ID, WORD, LSTR)
               CALL PRTERR ('CMDERR', 'Element block ID ' //
     &            WORD(:LSTR) // ' does not exist, ignored')
            END IF
  120       CONTINUE
            GOTO 110
         END IF

      ELSE IF ((VERB .EQ. 'BLOCKS') .OR. (VERB .EQ. 'MATERIAL')) THEN
         CALL FFADDC (VERB, INLINE(1))
         IF (VERB .EQ. 'MATERIAL') VERB = 'BLOCKS'

         IF (.NOT. FFEXST (IFLD, INTYP)) THEN

C         --Select all element blocks if no parameters
            CALL INILOG (NELBLK, .TRUE., SELELB)
            OPTION = '+'

         ELSE IF (.NOT. FFNUMB (IFLD, INTYP)) THEN

C         --Strip off ADD or DELETE option
            CALL FFCHAR (IFLD, INTYP, CFIELD, ' ', WORD)
            IF (MATSTR (WORD, 'ADD', 1)) THEN
               CALL FFADDC ('ADD', INLINE(1))
               OPTION = '+'
            ELSE IF (MATSTR (WORD, 'DELETE', 1)) THEN
               CALL FFADDC ('DELETE', INLINE(1))
               OPTION = '-'
            ELSE
               CALL PRTERR ('CMDERR',
     &            'Expected "ADD" or "DELETE" or element block ID')
               GOTO 160
            END IF

         ELSE

C         --De-select all element blocks so only listed blocks are selected
            CALL INILOG (NELBLK, .FALSE., SELELB)
            OPTION = '+'
         END IF

  130    CONTINUE
         IF (FFEXST (IFLD, INTYP)) THEN
            CALL FFINTG (IFLD, INTYP, IFIELD,
     &         'element block ID', 0, ID, *140)
            IELB = LOCINT (ID, NELBLK, IDELB)
            IF (IELB .GT. 0) THEN
               CALL FFADDI (ID, INLINE(1))
               SELELB(IELB) = (OPTION .EQ. '+')
            ELSE
               CALL INTSTR (1, 0, ID, WORD, LSTR)
               CALL PRTERR ('CMDERR', 'Element block ID ' //
     &            WORD(:LSTR) // ' does not exist, ignored')
            END IF
  140       CONTINUE
            GOTO 130
         END IF

         RETVRB = 'BLOCKS'

      ELSE IF ((VERB .EQ. 'LOG') .or. (verb .eq. 'SAVELOG')) THEN
         if (verb .eq. 'SAVELOG') then
            call prterr ('CMDSPEC', 'Please use the LOG command')
            verb = 'LOG'
         end if
         VERB = ' '
         RETVRB = 'LOG'

      ELSE IF (VERB .EQ. 'LIST') THEN
         VERB = ' '
         CALL FFCHAR (IFLD, INTYP, CFIELD, ' ', WORD)
         CALL DBLIST (WORD, A, NAMECO, BLKTYP, NAMES,
     &                TIMES, IDELB, QAREC, INFREC, MERR)
         IF (MERR .EQ. 1) RETURN

      ELSE IF (VERB .EQ. 'SHOW') THEN
         VERB = ' '
         CALL FFCHAR (IFLD, INTYP, CFIELD, ' ', WORD)
         CALL SHOW (WORD, CMDTBL, NAMECO, BLKTYP, NAMES,
     &        TIMES, IPTIMS, IDELB, VISELB, SELELB)

         IF (WORD .EQ. ' ') RETVRB = 'PRINT'

      ELSE IF (VERB .EQ. 'HELP') THEN
         VERB = ' '
         CALL ABRSTR (HLPTYP, WORD, HLPTBL)
         ISON = HELP ('ALGEBRA', HLPTYP, ' ')
         IF (.NOT. ISON) THEN
            IF (HLPTYP .EQ. 'RULES') THEN
               CONTINUE
            ELSE IF (HLPTYP .EQ. 'COMMANDS') THEN
               CALL SHOCMD ('COMMANDS', CMDTBL)
            ELSE IF (HLPTYP .EQ. 'FUNCTION') THEN
               CALL SHOCMD ('Available Functions:', FNCNAM)
            ELSE
               CALL SHOCMD ('HELP Options:', HLPTBL)
            END IF
         END IF

         RETVRB = 'PRINT'

      ELSE IF ((VERB .EQ. 'END') .OR. (VERB .EQ. 'EXIT')) THEN
         VERB = ' '
         CALL FFADDC ('END', INLINE(1))
         RETVRB = 'END'

      ELSE IF (VERB .EQ. 'QUIT') THEN
         VERB = ' '
         CALL FFADDC ('QUIT', INLINE(1))
         RETVRB = 'QUIT'

C   --This command allows selectable debugging
      else if (verb .eq. 'DEBUG') then
         verb = ' '
         call ffchar (ifld, intyp, cfield, ' ', cdebug)
         if (matstr (cdebug, 'EQUATION', 3)) then
            call prtdeb ('EQUATION', 0)
         else if (matstr (cdebug, 'VARIABLE', 3)) then
            call prtdeb ('VARIABLE', 0)
         else if (matstr (cdebug, 'ALL', 3)) then
            call prtdeb ('VARIABLE', 0)
            call prtdeb ('EQUATION', -1)
         end if

      ELSE
         CALL PRTERR ('CMDERR', '"' // VERB(:LENSTR(VERB))
     &      // '" is an invalid command')
         VERB = ' '
         GOTO 150
      END IF

      GOTO 160

  150 CONTINUE
      INLINE(1) = ' '

  160 CONTINUE
      IF (VERB .NE. ' ') THEN
         CALL SHOW (VERB, CMDTBL, NAMECO, BLKTYP, NAMES,
     &      TIMES, IPTIMS, IDELB, VISELB, SELELB)
      END IF

      RETURN
      END

max_stars_repo_name: mathstuf/seacas

max_stars_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_issues_repo_name: mathstuf/seacas

max_issues_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_forks_repo_name: mathstuf/seacas

max_forks_repo_licenses: BSD-3-Clause NetCDF Zlib MIT


!*****************************************************************************************
!>
!  Test of Fortran Astrodynamics Toolkit - integrate a trajectory near the Earth,
!  using the GGM03C geopotential model.

    program gravity_test

    use fortran_astrodynamics_toolkit

    implicit none

    !spacecraft propagation type:
    ! extend the rk class to include data used in the deriv routine
    type,extends(rk8_10_class) :: spacecraft
        type(geopotential_model_mueller) :: grav   !central body geopotential model
        logical :: first = .true.  !first point is being exported
    end type spacecraft

    integer,parameter  :: n=6    !number of state variables
    real(wp),parameter :: mu_earth  = 398600.4418_wp   !grav. param of earth [km3/s2]
    real(wp),parameter :: orbit_inc = 45.0_wp*deg2rad  !orbit inclination [rad]
    real(wp),parameter :: orbit_sma = 6778.0_wp        !orbit semimajor axis [km]

    !the coefficient file:
    character(len=*),parameter :: gravfile = '../grav/GGM03C.GEO'
    integer,parameter :: grav_n = 10    !max degree
    integer,parameter :: grav_m = 10    !max order

    type(spacecraft) :: s
    real(wp) :: t0,tf,x0(n),dt,xf(n),x02(n)
    logical :: status_ok
    integer :: iunit,istat

    !open output file:
    open(newunit=iunit,file='gravity_test.txt',status='REPLACE',iostat=istat)

    !constructor (main body is Earth):
    call s%initialize(n=n,f=deriv,report=twobody_report)
    call s%grav%initialize(gravfile,grav_n,grav_m,status_ok)
    if (.not. status_ok) stop 'Error'

    !get an initial cartesian state vector:
    ! [circular, 45 deg inclination, radius of 6778 km]
    call els3pv(mu_earth,[mu_earth/orbit_sma,orbit_sma,orbit_inc,zero,zero,zero],x0)

    write(iunit,'(A/,*(E30.16/))') 'Initial state:',x0

    !initial conditions:
    t0 = zero         !initial time (sec)
    dt = 100.0_wp     !time step (sec)
    tf = two*day2sec  !final time (sec)

    !integrate:
    s%first = .true.
    call s%integrate(t0,x0,dt,tf,xf)
    write(iunit,*) ''
    write(iunit,'(A/,*(E30.16/))') 'Final state:',xf

    !cleanup:
    call s%grav%destroy()
    close(iunit,iostat=istat)

!*****************************************************************************************
    contains

    !*********************************************************
        subroutine deriv(me,et,x,xdot)

        !! derivative routine

        implicit none

        class(rk_class),intent(inout)         :: me
        real(wp),intent(in)                   :: et    !! this is ephemeris time
        real(wp),dimension(me%n),intent(in)   :: x
        real(wp),dimension(me%n),intent(out)  :: xdot

        real(wp),dimension(3) :: r,rb,v,a_grav
        real(wp),dimension(3,3) :: rotmat

        select type (me)
        class is (spacecraft)

            !get state:
            r = x(1:3)
            v = x(4:6)

            !input state is inertial frame, have to convert
            ! to body-fixed Earth frame:
            rotmat = icrf_to_iau_earth(et)  !rotation matrix
            rb = matmul(rotmat,r)    !r in body-fixed frame

            !get the acceleration due to the geopotential:
            call me%grav%get_acc(rb,grav_n,grav_m,a_grav)

            !convert acceleration back to inertial frame:
            a_grav = matmul(transpose(rotmat),a_grav)

            !derivative vector:
            xdot(1:3) = v
            xdot(4:6) = a_grav

        end select

        end subroutine deriv
    !*********************************************************

    !*********************************************************
        subroutine twobody_report(me,t,x)

        !! report function - write time,state to console

        implicit none

        class(rk_class),intent(inout)        :: me
        real(wp),intent(in)                  :: t
        real(wp),dimension(me%n),intent(in)  :: x

        select type (me)
        class is (spacecraft)
            if (me%first) then  !print header
                write(iunit,*) ''
                write(iunit,'(*(A15,1X))')  't (sec)','x (km)','y (km)','z (km)',&
                                        'vx (km/s)','vy (km/s)','vz (km/s)'
                me%first = .false.
            end if
        end select

        write(iunit,'(*(F15.6,1X))') t,x

        end subroutine twobody_report
    !*********************************************************

    end program gravity_test
!*****************************************************************************************


text:
C    Copyright(C) 1999-2020 National Technology & Engineering Solutions
C    of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C    NTESS, the U.S. Government retains certain rights in this software.
C
C    See packages/seacas/LICENSE for details

      SUBROUTINE PDATA (MP, ML, MR, MSC, IPOINT, COOR, IPBOUN, ILINE,
     &   LTYPE, NINT, LCON, FACTOR, ILBOUN, ISBOUN, IREGN, IMAT, LINKP,
     &   LINKL, LINKR, LINKSC, RSIZE, SCHEME, DEFSCH, DEFSIZ, REXTRM,
     &   N, LABP, LABL, LABR, FULL, LABMD, LABI, LABF, LABPB, LABLB,
     &   LABSBD, LABSC, LABSZ, AXISD, TITLE, XMIN, XMAX, YMIN, YMAX,
     &   XX1, YY1, XX2, YY2, DEV1, VERSN)
C***********************************************************************

C  SUBROUTINE PDATA = PLOTS FLAGGED POINTS,  LINES,  AND REGIONS

C***********************************************************************

      DIMENSION IPOINT (MP), COOR (2, MP), IPBOUN (MP)
      DIMENSION ILINE (ML), LTYPE (ML), NINT (ML), LCON (3, ML)
      DIMENSION FACTOR (ML)
      DIMENSION ILBOUN (ML), ISBOUN (ML)
      DIMENSION IREGN (MR), IMAT (MR), REXTRM (4, MR), RSIZE (MR)
      DIMENSION SCHEME (MSC)
      DIMENSION LINKP (2, MP), LINKL (2, ML), LINKR (2, MR)
      DIMENSION LINKSC (2, MR)
      DIMENSION N (29), XDUM (2), YDUM (2)

      CHARACTER*72 DUMMY, SCHEME, DEFSCH, TITLE, DEV1*3
      CHARACTER*8 DATE, TIME, VERSN*10

      LOGICAL LABP, LABL, LABR, AXISD, LABMD, LABI, LABF
      LOGICAL LABPB, LABLB, LABSBD
      LOGICAL ADDLNK, CPUIFC, TEST, FULL, LABSC
      LOGICAL GETMAX, ADD, LABSZ

C  INITIALIZE THE PLOTTING SURFACE

      TEST = .FALSE.
      GETMAX = .FALSE.
      IF (TEST)OPEN (UNIT = 12, FILE = 'HP7580.DAT', STATUS = 'NEW')
      ADDLNK = .FALSE.
      CALL PLTBGN
      CALL PLTSTV (2, 160.)
      XDIMR = ABS (XMAX - XMIN)
      YDIMR = ABS (YMAX - YMIN)
      XDIMD = 1.
      YDIMD = .75
      CALL MPVIEW (0., XDIMD, 0., YDIMD)
      XRAT = XDIMR/XDIMD
      YRAT = YDIMR/YDIMD
      IF (XRAT.LT.YRAT) THEN
         XDIMR = XDIMD*YRAT
         XX1 =  (XMIN + XMAX - XDIMR)*.5
         XX2 =  (XMIN + XMAX + XDIMR)*.5
         XDIMR = XX2 - XX1
         YY1 = YMIN
         YY2 = YMAX
      ELSE
         YDIMR = YDIMD*XRAT
         YY1 =  (YMIN + YMAX - YDIMR)*.5
         YY2 =  (YMIN + YMAX + YDIMR)*.5
         YDIMR = YY2 - YY1
         XX1 = XMIN
         XX2 = XMAX
      ENDIF

C  SET UP SCALING EXTREMES FOR AXIS

      IF (TEST) THEN
         WRITE (12, 10000)'IN;SP6;;IP - 5710, -10060, 15710, 10060;'
         WRITE (12, 10010)
     &      'SC', INT (XX1*1000), ', ', INT (YY1*1000), ', ',
     &      INT (XX2*1000), ', ', INT (YY2*1000), ';'
      ENDIF
      IF (AXISD) THEN
         XDUM (1) = XX1 -  (XDIMR*.05)
         XDUM (2) = XX2 +  (XDIMR*.05)
         YDUM (1) = YY1 -  (YDIMR*.05)
         YDUM (2) = YY2 +  (YDIMR*.05)
         SHRINK = .2
      ELSE
         SHRINK = .1
      ENDIF

C  SHRINK TO FIT A BORDER ON THE PLOT

      XX1 = XX1 -  (XDIMR*SHRINK)
      XX2 = XX2 +  (XDIMR*SHRINK)
      YY1 = YY1 -  (YDIMR*SHRINK)
      YY2 = YY2 +  (YDIMR*SHRINK)
      CALL MPORT2 (XX1, XX2, YY1, YY2)
      CALL PLTFRM (0)

C  PLOT THE TITLE AND THE TRACE

      CALL STRLNG (TITLE, LEN)
      IF ( (LEN.GT.1) .OR. (TITLE (1:1).NE.' ')) THEN
         CALL PLTXHL (TITLE (1:LEN), XLEN)
         XBEGIN = AMAX1 (0.,  (XDIMD*.5 - XLEN*.5))
         CALL PLTXTH (XBEGIN, YDIMD*.95, TITLE (1:LEN))
      ENDIF
      DUMMY(1:10) = ' DRAWN BY '
      DUMMY(11:20) = VERSN
      DUMMY(21:22) = '  '
      CALL EXDATE (DATE)
      DUMMY(23:30) = DATE
      DUMMY(31:32) = '  '
      CALL EXTIME (TIME)
      DUMMY(33:40) = TIME
      CALL PLTXTH (0., 0., DUMMY(1:40))

C  DRAW THE AXIS IF REQUIRED,  AND SET CLIPPING WITHIN AXIS

      IF (AXISD)CALL SETAXS (XDUM, YDUM)
      IF (CPUIFC (.TRUE.))GOTO 130

C  PLOT THE POINTS FLAGGED

      IF ( (LABP) .OR. (LABPB)) THEN
         DO 100 I = 1, N (18)
            IF (CPUIFC (.TRUE.))GOTO 130
            CALL LTSORT (MP, LINKP, I, II, ADDLNK)
            IF (II.GT.0) THEN
               IF (IPOINT (II).LT.0) THEN
                  INUM =  - IPOINT (II)
                  CALL MP2PT (1, COOR (1, II), COOR (2, II),
     &               X1, Y1, MASK)
                  IF (MOD (MASK, 2).NE.0) THEN

C  PLOT THE POINT LABELS

                     IF (LABP) THEN
                        CALL PLTSTD (1, 1.)
                        CALL GETDUM (INUM, DUMMY, LEN)
                        CALL PLTXTH (X1, Y1, DUMMY (1:LEN))
                        CALL PLTXHE (X1, Y1)
                     ENDIF

C  PLOT THE POINBC FLAGS

                     IF ( ( (LABPB) .OR. ( (FULL) .AND. (LABP))) .AND.
     &                  (IPBOUN (II).GT.0)) THEN
                        CALL PLTSTD (1, 5.)
                        IF (LABP) THEN
                           CALL PLTXTH (X1, Y1, '/')
                           CALL PLTXHE (X1, Y1)
                        ENDIF
                        CALL GETDUM (IPBOUN (II), DUMMY, LEN)
                        CALL PLTXTH (X1, Y1, DUMMY (1:LEN))
                     ENDIF
                  ENDIF
               ENDIF
            ENDIF
  100    CONTINUE
      ENDIF

C  PLOT ALL LINES THAT HAVE BEEN FLAGGED

      DO 110 I = 1, N (19)
         IF (CPUIFC (.TRUE.))GOTO 130
         CALL LTSORT (ML, LINKL, I, II, ADDLNK)
         IF (II.GT.0) THEN
            IF (LABL) THEN
               ADD = .TRUE.
            ELSE
               ADD = .FALSE.
            ENDIF
            CALL PLTSTD (1, 7.)
            IF (ILINE (II).LT.0) THEN
               IF ( (LABL) .OR. (LABLB) .OR. (LABSBD) .OR. (LABF)
     &            .OR. (LABI)) THEN
                  KNUM =  - ILINE (II)
               ELSE
                  KNUM = 0
               ENDIF
               LT = LTYPE (II)
               IP1 = LCON (1, II)
               IP2 = LCON (2, II)
               IP3 = LCON (3, II)
               CALL LTSORT (MP, LINKP, IP1, IPNTR1, ADDLNK)
               CALL LTSORT (MP, LINKP, IP2, IPNTR2, ADDLNK)
               IF (IP3.NE.0) THEN
                  CALL LTSORT (MP, LINKP, IABS (IP3), IPNTR3, ADDLNK)
               ELSE
                  IPNTR3 = 0
               ENDIF
               IF ((IPNTR1.GT.0) .AND. (IPNTR2.GT.0) .AND.
     &            ((LT.EQ.1) .OR. (IPNTR3.GT.0)) ) THEN
                  CALL DLINE (MP, ML, COOR, LINKP, KNUM, LT, IP1, IP2,
     &               IP3, LABL, X1, Y1, TEST, GETMAX, DUM1, DUM2, DUM3,
     &               DUM4)

C  PLOT INTERVAL NUMBERS

                  IF ( ( (FULL) .AND. (LABL)) .OR. (LABI)) THEN
                     CALL PLTSTD (1, 5.)
                     IF (ADD) THEN
                        CALL PLTXHE (X1, Y1)
                        CALL PLTXTH (X1, Y1, '/')
                        CALL PLTXHE (X1, Y1)
                     ENDIF
                     CALL GETDUM (NINT (II), DUMMY, LEN)
                     IF (TEST) THEN
                        CALL PLTD2G (X1, Y1, XR, YR)
                        CALL PLTG2D (X1, Y1, XR, YR)
                        WRITE (12, 10020)'PU;PA', INT (XR*1000.),
     &                     ', ', INT (YR*1000.), ';LB',
     &                     DUMMY (1:LEN), CHAR (3)
                     ENDIF
                     CALL PLTXTH (X1, Y1, DUMMY (1:LEN))
                     ADD = .TRUE.
                  ENDIF

C  PLOT THE LINE FACTOR

                  IF ( ( (FULL) .AND. (LABL)) .OR. (LABF)) THEN
                     IF (ADD) THEN
                        CALL PLTSTD (1, 1.)
                        CALL PLTXHE (X1, Y1)
                        CALL PLTXTH (X1, Y1, '/')
                        CALL PLTXHE (X1, Y1)
                     ENDIF
                     CALL GTXDUM (FACTOR (II), DUMMY, LEN)
                     CALL PLTXTH (X1, Y1, DUMMY (1:LEN))
                     ADD = .TRUE.
                  ENDIF

C  PLOT THE LINEBC FLAGS

                  IF ( ( ( (FULL) .AND. (LABL)) .OR. (LABLB)) .AND.
     &               (ILBOUN (II).GT.0)) THEN
                     CALL PLTSTD (1, 2.)
                     IF (ADD) THEN
                        CALL PLTXHE (X1, Y1)
                        CALL PLTXTH (X1, Y1, '/')
                        CALL PLTXHE (X1, Y1)
                     ENDIF
                     CALL GETDUM (ILBOUN (II), DUMMY, LEN)
                     CALL PLTXTH (X1, Y1, DUMMY (1:LEN))
                     ADD = .TRUE.
                  ENDIF

C  PLOT THE SIDEBC FLAGS

                  IF ( ( ( (FULL) .AND. (LABL)) .OR. (LABSBD)) .AND.
     &               (ISBOUN (II).GT.0)) THEN
                     CALL PLTSTD (1, 3.)
                     IF (ADD) THEN
                        CALL PLTXHE (X1, Y1)
                        CALL PLTXTH (X1, Y1, '/')
                        CALL PLTXHE (X1, Y1)
                     ENDIF
                     CALL GETDUM (ISBOUN (II), DUMMY, LEN)
                     CALL PLTXTH (X1, Y1, DUMMY (1:LEN))
                  ENDIF
               ENDIF
            ENDIF
         ENDIF
  110 CONTINUE

C  PLOT ALL REGIONS FLAGGED

      IF ( (LABR) .OR. (LABMD) .OR. (LABSC) .OR. (LABSZ)) THEN
         IF (CPUIFC (.TRUE.))GOTO 130
         DO 120 I = 1, N (22)
            CALL LTSORT (MR, LINKR, I, II, ADDLNK)
            IF (II.GT.0) THEN
               IF (IREGN (II).LT.0) THEN
                  ADD = .FALSE.
                  INUM =  - IREGN (II)
                  XMID =  (REXTRM (1, II) + REXTRM (2, II))/2.
                  YMID =  (REXTRM (3, II) + REXTRM (4, II))/2.
                  CALL MP2PT (1, XMID, YMID, X1, Y1, MASK)
                  IF ( (MOD (MASK, 2).NE.0)) THEN

C  PLOT THE REGION NUMBER

                     IF (LABR) THEN
                        CALL PLTSTD (1, 2.)
                        CALL GETDUM (INUM, DUMMY, LEN)
                        CALL PLTXTH (X1, Y1, DUMMY (1:LEN))
                        ADD = .TRUE.
                     ENDIF

C  PLOT OUT THE MATERIAL NUMBER

                     IF (((FULL) .AND. (LABR)) .OR. (LABMD)) THEN
                        CALL PLTSTD (1, 1.)
                        IF (ADD) THEN
                           CALL PLTXHE (X1, Y1)
                           CALL PLTXTH (X1, Y1, '/')
                           CALL PLTXHE (X1, Y1)
                        ENDIF
                        ADD = .TRUE.
                        CALL GETDUM (IMAT (II), DUMMY, LEN)
                        CALL PLTXTH (X1, Y1, DUMMY (1:LEN))
                     ENDIF

C  PLOT OUT THE SIZE NUMBER FOR THE REGION

                     IF (((FULL) .AND. (LABR)) .OR. (LABSZ)) THEN
                        CALL PLTSTD (1, 1.)
                        IF (ADD) THEN
                           CALL PLTXHE (X1, Y1)
                           CALL PLTXTH (X1, Y1, '/')
                           CALL PLTXHE (X1, Y1)
                        ENDIF
                        ADD = .TRUE.
                        CALL GTXDUM (RSIZE (II), DUMMY, LEN)
                        CALL PLTXTH (X1, Y1, DUMMY (1:LEN))
                     ENDIF

C  PLOT OUT THE SCHEME

                     IF (((FULL) .AND. (LABR)) .OR. (LABSC)) THEN
                        CALL PLTSTD (1, 7.)
                        IF (ADD) THEN
                           CALL PLTXHE (X1, Y1)
                           CALL PLTXTH (X1, Y1, '/')
                           CALL PLTXHE (X1, Y1)
                        ENDIF
                        CALL LTSORT (MR, LINKSC, INUM, IPNTR, ADDLNK)
                        IF ( (INUM.LE.N (24)) .AND. (IPNTR.GT.0)) THEN
                           CALL STRLNG (SCHEME (IPNTR), LEN)
                           IF (TEST) THEN
                              CALL PLTD2G (X1, Y1, XR, YR)
                              WRITE (12, 10020)'PU;PA', INT (XR*1000.),
     &                           ', ', INT (YR*1000.), ';LB',
     &                           SCHEME (IPNTR) (1:LEN), CHAR (3)
                           ENDIF
                           CALL PLTXTH (X1, Y1, SCHEME (IPNTR) (1:LEN))
                        ELSE
                           CALL STRLNG (DEFSCH, LEN)
                           IF (TEST) THEN
                              CALL PLTD2G (X1, Y1, XR, YR)
                              WRITE (12, 10020)'PU;PA', INT (XR*1000.),
     &                           ', ', INT (YR*1000.), ';LB',
     &                           DEFSCH (1:LEN), CHAR (3)
                           ENDIF
                           CALL PLTXTH (X1, Y1, DEFSCH (1:LEN))
                        ENDIF
                     ENDIF
                  ENDIF
               ENDIF
            ENDIF
  120    CONTINUE
      ENDIF
  130 CONTINUE
      CALL PLTSTD (1, 7.)
      CALL PLTBEL
      CALL PLTFLU

      RETURN

10000 FORMAT (A)
10010 FORMAT (A2, I10, A1, I10, A1, I10, A1, I10, A1)
10020 FORMAT (A5, I10, A1, I10, A3, A, A1)

      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas