! { dg-do compile }
! { dg-additional-options "-O -fno-tree-loop-ivcanon -fno-tree-ccp -fno-tree-ch -finline-small-functions" }

module foo
  INTEGER, PARAMETER :: ONE = 1
end module foo
program test
  use foo
  integer :: a(ONE), b(ONE), c(ONE), d(ONE)
  interface
    function h_ext()
    end function h_ext
  end interface
  c = j()
  if (any (c .ne. check)) call myabort (7)
contains
  function j()
     integer :: j(ONE), cc(ONE)
     j = cc - j
  end function j
  function get_d()
  end function get_d
end program test


text:
	SUBROUTINE GBOUND ( imn, imx, jmn, jmx,
     +			    grid, hist, kx, ky, iret )
C************************************************************************
C* GBOUND                                                               *
C*                                                                      *
C* This subroutine performs the GEMPAK bounds search for the 		*
C* graph-to-grid algorithm. The bounds polygons are saved in the 	*
C* graph-to-grid arrays as closed contour lines with smoothing IMISSD	*
C* and contour value RMISSD. Grid points within the bounds areas are	*
C* set to RMISSD and the corresponding history value is set to BOUNDED.	*
C*                                                                      *
C* Contour lines may be used to simulate boundaries. Closed lines with	*
C* the label "-9999" will mask all points inside their area; lines with	*
C* "9999" will mask all points outside their area.			*
C*                                                                      *
C* GBOUND ( IMN, IMX, JMN, JMX, GRID, HIST, KX, KY, IRET )		*
C*                                                                      *
C* Input parameters:                                                    *
C*      IMN		INTEGER		Lower I range of GRID		*
C*      IMX		INTEGER		Upper I range of GRID		*
C*      JMN		INTEGER		Lower J range of GRID		*
C*      JMX		INTEGER		Upper J range of GRID		*
C*      KX		INTEGER		Number of grid points in X	*
C*      KY		INTEGER		Number of grid points in Y	*
C*                                                                      *
C* Input and Output parameters:                                         *
C*      HIST (KX,KY)	REAL		Grid to keep history		*
C*                                                                      *
C* Output parameters:                                                   *
C*      GRID (KX,KY)	REAL		Grid to calculate point values	*
C*      IRET            INTEGER         Return code                     *
C**                                                                     *
C* Log:                                                                 *
C* W.D.Plummer/NCEP     12/02                                           *
C* W.D.Plummer/NCEP     08/03	Treat missing-value lines as bounds	*
C* W.D.Plummer/NCEP     12/03	Bug fix for BOUNDS processing		*
C* D.W.Plummer/NCEP     02/04	Bug fix for BOUNDS w/ > MAXPPL points	*
C************************************************************************
C
        INCLUDE         'GEMPRM.PRM'
	INCLUDE		'grphgd.cmn'
        INCLUDE         'ERROR.PRM'
C
	REAL		grid(kx,ky), hist(kx,ky)
C
	INTEGER		inout(1200), qinout
	REAL		fiq(1200), fjq(1200)
	REAL            xlat (LLMXPT), ylon (LLMXPT)
	CHARACTER       sys_G*2, sys_M*2, sys_N*2
	CHARACTER*32    btyp, btag, carr(3), barr(10)
C
	LOGICAL		proces, done
C-----------------------------------------------------------------------
	iret = NORMAL
	qinout = GG_IN
	sys_G = 'G'
	CALL ST_NULL ( sys_G, sys_G, lens, ier )
	sys_M = 'M'
	CALL ST_NULL ( sys_M, sys_M, lens, ier )
	sys_N = 'N'
	CALL ST_NULL ( sys_N, sys_N, lens, ier )
C
	CALL CLO_INIT ( iret )
C
C*      Set the geographical limits for searching.
C
        CALL GQBND ( sys_M, rlatmn, rlonmn, dlatmx, dlonmx, iret )
        CALL CLO_BSAREA(rlatmn, rlonmn, dlatmx, dlonmx, iret)
C
	nprow = imx - imn + 1
C
C*      Loop over bounds specifications, setting the bounds type
C*      (and tag) for searching.
C
        CALL ST_CLST ( bnds, '+', ' ', 10, barr, numb, ier )
C
	DO  nb = 1, numb
C
	  proces = .true.
C
          CALL ST_CLST ( barr(nb), '|', ' ', 3, carr, num, ier )
          CALL ST_LSTR ( carr(1), lens, ier )
          btyp = carr(1)(1:lens) // CHNULL
          CALL ST_LCUC ( btyp, btyp, ier )
C
	  CALL CLO_BSTYPE( btyp, iret )
          IF  ( iret .ne. NORMAL )  proces = .false.
C
	  IF ( proces )  THEN
C
            CALL ST_LSTR ( carr(2), lentag, ier )
            IF ( lentag .gt. 0 )  THEN
              btag = carr(2)(1:lentag) // CHNULL
              CALL CLO_BSTAG( btag, iret )
              IF  ( iret .ne. NORMAL )  proces = .false.
            END IF
C
          END IF
C
	  IF ( proces )  THEN
C
C*	    Check third parm for TRUE or otherwise.
C*	    This controls whether to assign RMISSD to the grid (and
C*	    BOUNDED to the history grid) to gridpoints inside (TRUE, default)
C*	    or outside (otherwise) the bound area.
C*	    Logic: Assign BOUNDED to the history grid if the grid
C*	    point is inside a bounds part. After all bounds parts have been
C*	    checked, apply the TRUE/FALSE logic. Note that this procedure
C*	    will only work when bounds checking is performed before any other
C*	    grid point value assignment.
C
            CALL ST_LSTR ( carr(3), lenlog, ier )
	    IF ( lenlog .gt. 0 )  THEN
              CALL ST_LCUC ( carr(3), carr(3), ier )
	      IF ( carr(3)(1:1) .ne. 'T' )  qinout = GG_OUT
	    END IF
C
	  END IF
C
	  IF ( proces )  THEN
C
C*	   Loop over gridpoint rows.
C
	   minp   = 0
	   mxpts  = LLMXPT
	   filter = 0.0
	   DO  i = imn, imx
	     fiq(i) = i
	   END DO
C
	   ier = 0
	   DO WHILE ( ier .eq. 0 )
C
	    CALL CLO_BGNEXT ( minp, mxpts, filter, npoly,
     +			      xlat, ylon, ier)
C
            IF ( ier .eq. 0 ) THEN
C
	      IF ( ( xlat(1) .ne. xlat(npoly) ) .or.
     +             ( ylon(1) .ne. ylon(npoly) ) )  THEN
	        npoly = npoly + 1
	        xlat(npoly) = xlat(1)
	        ylon(npoly) = ylon(1)
	      END IF
C
	      DO  i = imn, imx
	        fiq(i) = i
	      END DO
C
	      DO  j = jmn, jmx
C
	        DO  i = imn, imx
	          fjq(i) = j
	        END DO
C
	        CALL CGR_INPOLY ( sys_G, nprow, fiq, fjq,
     +			        sys_M, npoly, xlat, ylon,
     +			        inout, ier )
C
C*		Flag all gridpoints inside as BOUNDED
C
	        DO  i = imn, imx
	          IF ( inout(i) .eq. GG_IN )  THEN
		    hist(i,j) = BOUNDED
		    grid(i,j) = RMISSD
	          END IF
	        END DO
C
	      END DO
C
	      nps = 1
	      npe = MIN ( (nps+MAXPPL-1), npoly )
	      done = .false.
	      DO WHILE ( .not. done )
		IF ( nlines .eq. MAXLIN )  THEN
		    done = .true.
		    iret = +4
		ELSE
	            nlines = nlines + 1
	            npts(nlines) = npe - nps + 1
	            value(nlines) = RMISSD
	            ismth(nlines) = IMISSD
	            closed(nlines) = CLSD
	            n = 0
	            DO  np = nps, npe
		        n = n + 1
	                flat(n, nlines) = xlat(np)
	                flon(n, nlines) = ylon(np)
	            END DO
		    if ( npe .lt. npoly )  THEN
		        nps = npe
		        npe = MIN ( (nps+MAXPPL-1), npoly )
		    ELSE
		        done = .true.
		    END IF
		END IF
	      END DO
C
	    END IF
C
	  END DO
C
C*	  If user wants all points outside of bounds areas to be missing,
C*	  then reverse the above assigments.
C*        Note that all history gridpoints should be either BOUNDED or INIT.
C
	  IF ( qinout .eq. GG_IN )  THEN
	      DO  i = imn, imx
	          DO  j = jmn, jmx
	              IF ( hist(i,j) .eq. BOUNDED )  THEN
			  hist(i,j) = INIT
		      ELSE IF ( hist(i,j) .eq. INIT )  THEN
			  hist(i,j) = BOUNDED
		      END IF
	          END DO
	      END DO
          END IF
C
	 END IF
C
	END DO
C
	DO  nl = 1, nlines
C
	  IF ( ( value(nl).eq.RMISSD .or. value(nl).eq.-RMISSD ) .and.
     &		ismth(nl).ne.IMISSD )  THEN
C
	    DO  nn = 1, npts(nl)
	      xlat(nn) = flat(nn,nl)
	      ylon(nn) = flon(nn,nl)
	    END DO
C
C*	    Loop over gridpoint rows.
C
	    DO  i = imn, imx
	      fiq(i) = i
	    END DO
C
	    DO  j = jmn, jmx
C
	      DO  i = imn, imx
	        fjq(i) = j
	      END DO
C
	      CALL CGR_INPOLY ( sys_G, nprow, fiq, fjq,
     +			        sys_M, npts(nl), xlat, ylon,
     +			        inout, ier )
C
	      DO  i = imn, imx
	        IF ( (value(nl).eq. RMISSD .and. inout(i).eq.1) .or.
     &		     (value(nl).eq.-RMISSD .and. inout(i).eq.0) )  THEN
		  hist(i,j) = BOUNDED
		  grid(i,j) = RMISSD
	        END IF
	      END DO
C
	    END DO
C
	    IF ( value(nl).eq.-RMISSD )  value(nl) = RMISSD
C
	  END IF
C
	END DO
C
	RETURN
C
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


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

C=======================================================================
      LOGICAL FUNCTION LXGTBP(STR,NS,CH)
      IMPLICIT INTEGER (A-Z)
      CHARACTER*504 ILINE
      COMMON /LXCOM1/ILINE
      COMMON /LXCOM2/JLINE,LXINIT
      CHARACTER*(*) STR
      CHARACTER CH
      LOGICAL QFLAG
      CHARACTER QCH

      NS = 0
      CH = ILINE(JLINE:JLINE)
      LXGTBP = (CH.EQ.'(') .OR. (CH.EQ.'[') .OR. (CH.EQ.'{')
      IF (.NOT.LXGTBP) THEN
         RETURN

      END IF

      J = JLINE
      LP = 1
      QFLAG = .FALSE.
 2440 CONTINUE
      J = J + 1
      CH = ILINE(J:J)
      IF (.NOT.QFLAG) THEN
         IF (CH.EQ.'(' .OR. CH.EQ.'[' .OR. CH.EQ.'{') THEN
            LP = LP + 1

         ELSE IF (CH.EQ.')' .OR. CH.EQ.']' .OR. CH.EQ.'}') THEN
            LP = LP - 1

         ELSE IF (CH.EQ.'''' .OR. CH.EQ.'"' .OR. CH.EQ.CHAR(96)) THEN
            QFLAG = .TRUE.
            QCH = CH

         ELSE IF (CH.EQ.CHAR(0)) THEN
            LXGTBP = .FALSE.
            RETURN

         END IF

         IF (LP.EQ.0) THEN
            GO TO 2460

         END IF

         NS = NS + 1
         STR(NS:NS) = CH

      ELSE
         NS = NS + 1
         STR(NS:NS) = CH
         IF (CH.EQ.QCH) THEN
            IF (CH.EQ.ILINE(J+1:J+1)) THEN
               NS = NS + 1
               STR(NS:NS) = CH
               J = J + 1
               GO TO 2450

            ELSE
               QFLAG = .FALSE.
            END IF

         ELSE IF (CH.EQ.CHAR(0)) THEN
            LXGTBP = .FALSE.
            RETURN

         END IF

      END IF

 2450 GO TO 2440

 2460 CONTINUE
      JLINE = J + 1
      LXGTBP = .TRUE.
      RETURN

      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


!!  N-dimensional complex system of arrays encapsulation
!
! This file is part of LIBPFASST.
!

!> !> Module to define and encapsulation for a system of N-ndimensional complex arrays
!!
!! When a new solution is created by a PFASST level, this encapsulation
!! uses the levels 'lev_shape' attribute to create a new multi-component array with that
!! shape.  Thus, the 'lev_shape' attributes of the PFASST levels should be
!! set appropriately.  The last component of lev_shape is the number of components in the system
!!
!! For example, before calling pf_pfasst_run we can
!! set the arr_shape of the coarsest level by doing:
!!
!!   allocate(pf%levels(1)%lev_shape(3))
!!   pf%levels(1)%lev_shape = [ nx, ny, 3 ]
!!
!! Which would imply that a 3 component system of two-dimensional solutions.
!!
!! The helper routines get_array1d, get_array2d, get_array3d, etc can be used to
!! extract pointers to a component of  encapsulated system without
!! performing any copies, or return a pointer to the whole array.
!!
module pf_mod_zndsysarray
  use iso_c_binding
  use pf_mod_dtype
  use pf_mod_stop
  implicit none

  !>  Type to create and destroy the arrays
  type, extends(pf_factory_t) :: pf_zndsysarray_factory_t
   contains
     procedure :: create_single  => zndsysarray_create_single
     procedure :: create_array  => zndsysarray_create_array
     procedure :: destroy_single => zndsysarray_destroy_single
     procedure :: destroy_array => zndsysarray_destroy_array
  end type pf_zndsysarray_factory_t

  !>  Type to extend the abstract encap and set procedure pointers
  type, extends(pf_encap_t) :: pf_zndsysarray_t
     integer             :: ndim    !  The spatial dimension of each component in system
     integer             :: ncomp  !  The number of components in the system
     integer             :: ndof   !  The number of variables in each component
     integer,    allocatable :: arr_shape(:)
     complex(pfdp), allocatable :: flatarray(:)
   contains
     procedure :: setval => zndsysarray_setval
     procedure :: copy => zndsysarray_copy
     procedure :: norm => zndsysarray_norm
     procedure :: pack => zndsysarray_pack
     procedure :: unpack => zndsysarray_unpack
     procedure :: axpy => zndsysarray_axpy
     procedure :: eprint => zndsysarray_eprint
     procedure, private  :: get_sysarray_1d, get_sysarray_2d,get_sysarray_3d,get_sysarray_4d
     generic :: get_sysarray => get_sysarray_1d ,get_sysarray_2d,get_sysarray_3d,get_sysarray_4d
     procedure, private  :: get_array_1d  ,get_array_2d,get_array_3d,get_array_4d
     generic :: get_array => get_array_1d ,get_array_2d,get_array_3d,get_array_4d
  end type pf_zndsysarray_t

contains
  !>  Subroutine to allocate the array and set the size parameters
  subroutine zndsysarray_build(q, arr_shape)
    class(pf_encap_t), intent(inout) :: q
    integer,           intent(in   ) :: arr_shape(:)

    integer :: ierr

    select type (q)
    class is (pf_zndsysarray_t)
       allocate(q%arr_shape(SIZE(arr_shape)),stat=ierr)
       if (ierr /=0) call pf_stop(__FILE__,__LINE__,'allocate fail, error=',ierr)

       q%ndim   = SIZE(arr_shape)-1
       q%ncomp = arr_shape(q%ndim+1)
       q%ndof = product(arr_shape(1:q%ndim))
       q%arr_shape = arr_shape

       allocate(q%flatarray(product(arr_shape)),stat=ierr)
       if (ierr /=0) call pf_stop(__FILE__,__LINE__,'allocate fail, error=',ierr)

    end select
  end subroutine zndsysarray_build

  !> Subroutine to  create a single array
  subroutine zndsysarray_create_single(this, x, level_index, lev_shape)
    class(pf_zndsysarray_factory_t), intent(inout)     :: this
    class(pf_encap_t),      intent(inout), allocatable :: x
    integer,                intent(in   )              :: level_index
    integer,                intent(in   )              :: lev_shape(:)
    integer :: i,ierr
    allocate(pf_zndsysarray_t::x,stat=ierr)
    if (ierr /=0) call pf_stop(__FILE__,__LINE__,'allocate fail, error=',ierr)
    call zndsysarray_build(x, lev_shape)
  end subroutine zndsysarray_create_single

  !> Subroutine to create an array of arrays
  subroutine zndsysarray_create_array(this, x, n, level_index,  lev_shape)
    class(pf_zndsysarray_factory_t), intent(inout)     :: this
    class(pf_encap_t),      intent(inout), allocatable :: x(:)
    integer,                intent(in   )              :: n
    integer,                intent(in   )              ::  level_index
    integer,                intent(in   )              ::  lev_shape(:)
    integer :: i,ierr
    allocate(pf_zndsysarray_t::x(n),stat=ierr)
    if (ierr /=0) call pf_stop(__FILE__,__LINE__,'allocate fail, error=',ierr)
    do i = 1, n
       call zndsysarray_build(x(i), lev_shape)
    end do
  end subroutine zndsysarray_create_array

  !>  Subroutine to destroy array
  subroutine zndsysarray_destroy(encap)
    class(pf_encap_t), intent(inout) :: encap
    type(pf_zndsysarray_t), pointer :: zndsysarray_obj

    zndsysarray_obj => cast_as_zndsysarray(encap)

    deallocate(zndsysarray_obj%arr_shape)
    deallocate(zndsysarray_obj%flatarray)

    nullify(zndsysarray_obj)

  end subroutine zndsysarray_destroy

  !> Subroutine to destroy an single array
  subroutine zndsysarray_destroy_single(this, x)
    class(pf_zndsysarray_factory_t), intent(inout)     :: this
    class(pf_encap_t),      intent(inout), allocatable :: x

    select type (x)
    class is (pf_zndsysarray_t)
       deallocate(x%arr_shape)
       deallocate(x%flatarray)
    end select
    deallocate(x)
  end subroutine zndsysarray_destroy_single

  !> Subroutine to destroy an array of arrays
  subroutine zndsysarray_destroy_array(this, x)
    class(pf_zndsysarray_factory_t), intent(inout)     :: this
    class(pf_encap_t),      intent(inout), allocatable :: x(:)
    integer                                            :: i

    select type(x)
    class is (pf_zndsysarray_t)
       do i = 1,SIZE(x)
          deallocate(x(i)%arr_shape)
          deallocate(x(i)%flatarray)
       end do
    end select
    deallocate(x)
  end subroutine zndsysarray_destroy_array

  !>  The following are the base subroutines that all encapsulations must provide
  !!

  !> Subroutine to set array to a scalare  value.
  subroutine zndsysarray_setval(this, val, flags)
    class(pf_zndsysarray_t), intent(inout)           :: this
    real(pfdp),     intent(in   )           :: val
    integer,        intent(in   ), optional :: flags
    this%flatarray = val
  end subroutine zndsysarray_setval

  !> Subroutine to copy an array
  subroutine zndsysarray_copy(this, src, flags)
    class(pf_zndsysarray_t),    intent(inout)           :: this
    class(pf_encap_t), intent(in   )           :: src
    integer,           intent(in   ), optional :: flags
    select type(src)
    type is (pf_zndsysarray_t)
       this%flatarray = src%flatarray
    class default
       call pf_stop(__FILE__,__LINE__,'unsupported type ')
    end select
  end subroutine zndsysarray_copy

  !> Subroutine to pack an array into a flat array for sending
  subroutine zndsysarray_pack(this, z, flags)
    class(pf_zndsysarray_t), intent(in   ) :: this
    real(pfdp),     intent(  out) :: z(:)
    integer,     intent(in   ), optional :: flags
    integer :: ntot
    ntot = this%ndof*this%ncomp

    z(1:ntot) = REAL(this%flatarray,pfdp)
    z(ntot+1:2*ntot) = AIMAG(this%flatarray)
  end subroutine zndsysarray_pack

  !> Subroutine to unpack a flatarray after receiving
  subroutine zndsysarray_unpack(this, z, flags)
    class(pf_zndsysarray_t), intent(inout) :: this
    real(pfdp),     intent(in   ) :: z(:)
    integer,     intent(in   ), optional :: flags

    integer :: ntot
    ntot = this%ndof*this%ncomp
    this%flatarray = z(1:ntot)
!    this%flatarray = this%flatarray + cmplx(0.0,1.0,pfdp)*z(ntot+1:2*ntot)
    this%flatarray =  cmplx(z(1:ntot),z(ntot+1:2*ntot),pfdp)
  end subroutine zndsysarray_unpack

  !> Subroutine to define the norm of the array (here the max norm)
  function zndsysarray_norm(this, flags) result (norm)
    class(pf_zndsysarray_t), intent(in   ) :: this
    integer,     intent(in   ), optional :: flags
    real(pfdp) :: norm
    norm = maxval(abs(this%flatarray))
  end function zndsysarray_norm

  !> Subroutine to compute y = a x + y where a is a scalar and x and y are arrays
  subroutine zndsysarray_axpy(this, a, x, flags)
    class(pf_zndsysarray_t),    intent(inout)           :: this
    class(pf_encap_t), intent(in   )           :: x
    real(pfdp),        intent(in   )           :: a
    integer,           intent(in   ), optional :: flags
    if (a .eq. 0.0_pfdp) return
    select type(x)
    type is (pf_zndsysarray_t)
       if (a .eq. 1.0_pfdp) then
          this%flatarray = x%flatarray + this%flatarray
       elseif (a .eq. -1.0_pfdp) then
          this%flatarray = -x%flatarray + this%flatarray
       else
          this%flatarray = a*x%flatarray + this%flatarray
       end if
    class default
       call pf_stop(__FILE__,__LINE__,'unsupported type ')
    end select
  end subroutine zndsysarray_axpy

  !>  Subroutine to print the array to the screen (mainly for debugging purposes)
  subroutine zndsysarray_eprint(this,flags)
    class(pf_zndsysarray_t), intent(inout) :: this
    integer,           intent(in   ), optional :: flags
    !  Just print the first few values
        print *, this%flatarray(1:10)
    !print *, this%flatarray
  end subroutine zndsysarray_eprint

  function cast_as_zndsysarray(encap_polymorph) result(zndsysarray_obj)
    class(pf_encap_t), intent(in), target :: encap_polymorph
    type(pf_zndsysarray_t), pointer :: zndsysarray_obj

    select type(encap_polymorph)
    type is (pf_zndsysarray_t)
       zndsysarray_obj => encap_polymorph
    end select
  end function cast_as_zndsysarray

  !>  Helper function to return the array part, these are called with get_array
  subroutine get_sysarray_1d(this,r,flags)
    class(pf_zndsysarray_t), target, intent(in) :: this
    complex(pfdp), pointer, intent(inout) :: r(:)
    integer,    intent(in   ), optional :: flags
    r => this%flatarray
  end subroutine get_sysarray_1d
  subroutine get_sysarray_2d(this,r,flags)
    class(pf_zndsysarray_t), target, intent(in) :: this
    complex(pfdp), pointer, intent(inout) :: r(:,:)
    integer,    intent(in   ), optional :: flags
    r(1:this%arr_shape(1),1:this%arr_shape(2)) => this%flatarray
  end subroutine get_sysarray_2d
  subroutine get_sysarray_3d(this,r,flags)
    class(pf_zndsysarray_t), target, intent(in) :: this
    complex(pfdp), pointer, intent(inout) :: r(:,:,:)
    integer,    intent(in   ), optional :: flags
    r(1:this%arr_shape(1),1:this%arr_shape(2),1:this%arr_shape(3)) => this%flatarray
  end subroutine get_sysarray_3d
  subroutine get_sysarray_4d(this,r,flags)
    class(pf_zndsysarray_t), target, intent(in) :: this
    complex(pfdp), pointer, intent(inout) :: r(:,:,:,:)
    integer,    intent(in   ), optional :: flags
    r(1:this%arr_shape(1),1:this%arr_shape(2),1:this%arr_shape(3),1:this%arr_shape(4)) => this%flatarray
  end subroutine get_sysarray_4d

  !>  Helper function to return the array part or whole
  function get_array1d(x,n,flags) result(r)
    class(pf_encap_t), target,intent(in) :: x
    integer, intent(in) :: n
    integer,           intent(in   ), optional :: flags
    complex(pfdp), pointer :: r(:)
    select type (x)
    type is (pf_zndsysarray_t)
       if (n .eq. 0) then  !  Return the whole array
          if (x%ncomp .eq. 1 .and. x%ndim .eq. 1) then
             r => x%flatarray
          else
             call pf_stop(__FILE__,__LINE__,'bad dimension, must be 1. ndim=',x%ndim)
          end if
       else                !  Return the nth component
          if (x%ndim .eq. 1) then
             r => x%flatarray(x%ndof*(n-1)+1:x%ndof*n)
          else
             call pf_stop(__FILE__,__LINE__,'bad dimension, must be 1. ndim=',x%ndim)
          end if
       end if
    end select
  end function get_array1d

  function get_array2d(x,n,flags) result(r)
    class(pf_encap_t), target,intent(in) :: x
    integer, intent(in) :: n
    integer,           intent(in   ), optional :: flags
    complex(pfdp), pointer :: r(:,:)

    select type (x)
    type is (pf_zndsysarray_t)
       if (n .eq. 0) then
          if (x%ndim .eq. 1) then
             r(1:x%arr_shape(1),1:x%arr_shape(2)) => x%flatarray
          else
             call pf_stop(__FILE__,__LINE__,'bad dimension, must be 1. ndim=',x%ndim)
          end if
       else                   !  Return the nth component
          if (x%ndim .eq. 2) then
             r(1:x%arr_shape(1),1:x%arr_shape(2)) => x%flatarray(x%ndof*(n-1)+1:x%ndof*n)
          else
             call pf_stop(__FILE__,__LINE__,'bad dimension, must be 2. ndim=',x%ndim)
          end if

       endif
    end select
  end function get_array2d

  function get_array3d(x,n,flags) result(r)
    class(pf_encap_t), target,intent(in) :: x
    integer, intent(in) :: n
    integer,           intent(in   ), optional :: flags
    complex(pfdp), pointer :: r(:,:,:)

    select type (x)
    type is (pf_zndsysarray_t)
       if (n .eq. 0) then        !  Return the whole array
          if (x%ndim .eq. 2) then
             r(1:x%arr_shape(1),1:x%arr_shape(2),1:x%arr_shape(3)) => x%flatarray
          else
             call pf_stop(__FILE__,__LINE__,'bad dimension, must be 2. ndim=',x%ndim)
          end if
       else                     !  Return the nth component
          if (x%ndim .eq. 3) then
             r(1:x%arr_shape(1),1:x%arr_shape(2),1:x%arr_shape(3)) => x%flatarray(x%ndof*(n-1)+1:x%ndof*n)
          else
             call pf_stop(__FILE__,__LINE__,'bad dimension, must be 3. ndim=',x%ndim)
          end if
       end if
    end select
  end function get_array3d

  function get_array4d(x,n,flags) result(r)
    class(pf_encap_t), target,intent(in) :: x
    integer, intent(in) :: n
    integer,           intent(in   ), optional :: flags
    complex(pfdp), pointer :: r(:,:,:,:)

    select type (x)
    type is (pf_zndsysarray_t)
       if (n .eq. 0) then        !  Return the whole array
          if (x%ndim .eq. 3) then
             r(1:x%arr_shape(1),1:x%arr_shape(2),1:x%arr_shape(3),1:x%arr_shape(4)) => x%flatarray
          else
             call pf_stop(__FILE__,__LINE__,'bad dimension, must be 3. ndim=',x%ndim)
          end if
       else                       !  Return the nth component
          if (x%ndim .eq. 4) then
             r(1:x%arr_shape(1),1:x%arr_shape(2),1:x%arr_shape(3),1:x%arr_shape(4)) => x%flatarray(x%ndof*(n-1)+1:x%ndof*n)
          else
             call pf_stop(__FILE__,__LINE__,'bad dimension, must be 4. ndim=',x%ndim)
          end if
       end if
    end select
  end function get_array4d

  !>  Helper function to return the array part or whole
  subroutine get_array_1d(x,r,n,flags)
    class(pf_zndsysarray_t), target,intent(in) :: x
    integer, intent(in) :: n
    integer,           intent(in   ), optional :: flags
    complex(pfdp), pointer :: r(:)
    if (n .eq. 0) then  !  Return the whole array
       if (x%ncomp .eq. 1 .and. x%ndim .eq. 1) then
          r => x%flatarray
          else
             call pf_stop(__FILE__,__LINE__,'bad dimension, must be 1. ndim=',x%ndim)
          end if
       else                !  Return the nth component
          if (x%ndim .eq. 1) then
             r => x%flatarray(x%ndof*(n-1)+1:x%ndof*n)
          else
             call pf_stop(__FILE__,__LINE__,'bad dimension, must be 1. ndim=',x%ndim)
          end if
       end if

  end subroutine get_array_1d

  subroutine get_array_2d(x,r,n,flags)
    class(pf_zndsysarray_t), target,intent(in) :: x
    integer, intent(in) :: n
    integer,           intent(in   ), optional :: flags
    complex(pfdp), pointer :: r(:,:)

    if (n .eq. 0) then
       if (x%ndim .eq. 1) then
          r(1:x%arr_shape(1),1:x%arr_shape(2)) => x%flatarray
       else
          call pf_stop(__FILE__,__LINE__,'bad dimension, must be 1. ndim=',x%ndim)
       end if
    else                   !  Return the nth component
       if (x%ndim .eq. 2) then
          r(1:x%arr_shape(1),1:x%arr_shape(2)) => x%flatarray(x%ndof*(n-1)+1:x%ndof*n)
       else
          call pf_stop(__FILE__,__LINE__,'bad dimension, must be 2. ndim=',x%ndim)
       end if

    endif
   end subroutine get_array_2d

  subroutine get_array_3d(x,r,n,flags)
    class(pf_zndsysarray_t), target,intent(in) :: x
    integer, intent(in) :: n
    integer,           intent(in   ), optional :: flags
    complex(pfdp), pointer :: r(:,:,:)

    if (n .eq. 0) then        !  Return the whole array
       if (x%ndim .eq. 2) then
          r(1:x%arr_shape(1),1:x%arr_shape(2),1:x%arr_shape(3)) => x%flatarray
       else
          call pf_stop(__FILE__,__LINE__,'bad dimension, must be 2. ndim=',x%ndim)
       end if
    else                     !  Return the nth component
       if (x%ndim .eq. 3) then
          r(1:x%arr_shape(1),1:x%arr_shape(2),1:x%arr_shape(3)) => x%flatarray(x%ndof*(n-1)+1:x%ndof*n)
       else
          call pf_stop(__FILE__,__LINE__,'bad dimension, must be 3. ndim=',x%ndim)
       end if
    end if

  end subroutine get_array_3d

  subroutine get_array_4d(x,r,n,flags)
    class(pf_zndsysarray_t), target,intent(in) :: x
    integer, intent(in) :: n
    integer,           intent(in   ), optional :: flags
    complex(pfdp), pointer :: r(:,:,:,:)

    select type (x)
    type is (pf_zndsysarray_t)
       if (n .eq. 0) then        !  Return the whole array
          if (x%ndim .eq. 3) then
             r(1:x%arr_shape(1),1:x%arr_shape(2),1:x%arr_shape(3),1:x%arr_shape(4)) => x%flatarray
          else
             call pf_stop(__FILE__,__LINE__,'bad dimension, must be 3. ndim=',x%ndim)
          end if
       else                       !  Return the nth component
          if (x%ndim .eq. 4) then
             r(1:x%arr_shape(1),1:x%arr_shape(2),1:x%arr_shape(3),1:x%arr_shape(4)) => x%flatarray(x%ndof*(n-1)+1:x%ndof*n)
          else
             call pf_stop(__FILE__,__LINE__,'bad dimension, must be 4. ndim=',x%ndim)
          end if
       end if
    end select
  end subroutine get_array_4d

end module pf_mod_zndsysarray


text:
module modA

use mod_kinds
use mod_global

IMPLICIT NONE

type, public :: type1
    integer(kind=irg)   :: a

    contains
        private
        procedure, pass(self) :: writeHDFnamelist_
        procedure, pass(self) :: p_

        generic, public :: writeHDFnamelist => writeHDFnamelist_
        generic, public :: p => p_

end type type1

contains

subroutine writeHDFnamelist_(self, n)

class(type1), INTENT(INOUT) :: self
integer(kind=irg),INTENT(IN):: n

write (*,*) '================'
write (*,*) 'type1:writeHDFnamelist: n = ', n
write (*,*) 'type1:writeHDFnamelist: a = ', self%a
write (*,*) '================'

end subroutine writeHDFnamelist_

subroutine p_(self, n)

class(type1), INTENT(INOUT) :: self
integer(kind=irg),INTENT(IN):: n

write (*,"(' ====>', I1,'<==== ')") n

end subroutine p_

end module modA

! ===================================================
module modB

use mod_kinds
use mod_global
use modA

IMPLICIT NONE

type, public, extends(type1) :: type2
    integer(kind=irg)   :: b

    contains
        procedure, pass(self) :: writeHDFnamelist_ => writelist

end type type2

contains

subroutine writelist(self, n)

class(type2), INTENT(INOUT) :: self
integer(kind=irg),INTENT(IN):: n

write (*,*) '================'
write (*,*) 'type2:writeHDFnamelist: n = ', n
write (*,*) 'type2:writeHDFnamelist: a = ', self%a
write (*,*) 'type2:writeHDFnamelist: b = ', self%b
write (*,*) '================'

end subroutine writelist

end module modB

program TestProgram

use mod_kinds
use mod_global
use modA
use modB
use mod_memory

! use mod_EMsoft
! use mod_io
! use mod_timing
! use mod_quaternions
! use mod_rng
use mod_HDFsupport
use ISO_C_BINDING
IMPLICIT NONE

type(type1)  :: t1
type(type2)  :: t2

character(fnlen)        :: progname = 'this is the program name'
character(kind=c_char)  :: Cprogname(fnlen)
integer(kind=irg)       :: i, slen
! character(fnlen)        :: progdesc = 'and this is the descriptor'
! character(fnlen)        :: m

! type(EMsoft_T)     :: EMsoft
! type(IO_T)         :: Message
! type(Timing_T)     :: Timing

! integer(kind=irg)       :: io_int(2), status
! type(rng_t)             :: seed
! type(QuaternionArray_T) :: qra , qrb

integer(kind=ish), allocatable :: ar(:)
real(kind=sgl), allocatable    :: ar3(:,:,:)
complex(kind=dbl), allocatable :: ar2(:,:)
character(fnlen), allocatable  :: str(:)
type(memory_T)                 :: mem, memth

! memth = memory_T(nt = 2)
mem = memory_T()
call mem%toggle_verbose()

call mem%alloc( ar, (/250/), 'ar', 15_ish, startdims = (/ -250 /) )
!call memth%alloc( ar3, (/50, 60, 70/), 'ar3', 10.0_sgl, TID=1 )
call mem%alloc( ar2, (/20, 20/), 'ar2', (0.D0, -1.D0), startdims = (/-20, -20/) )
call mem%alloc( str, (/ 15 /), 'str', 'test')

write (*,*) ar(-5:5)
write (*,*) ar2(1,6)
! write (*,*) ar3(1,2,3)
write (*,*) trim(str(5))

call mem%allocated_memory_use()
! call memth%thread_memory_use()

call mem%dealloc( ar, 'ar' )
! call memth%dealloc3( ar3, 'ar3', TID=1)
! call mem%thread_memory_use()
call mem%dealloc( ar2, 'ar2')
call mem%dealloc(str, 'str')

call mem%allocated_memory_use()
! call memth%allocated_memory_use()

! write (*,*) '->'//trim(progname)//'<-', len_trim(progname)

! Cprogname = carstringify(progname)

! do i=1,len_trim(progname)
!     write (*,"(A1$)") Cprogname(i:i)
! end do
! write (*,*) ''

! progname = ''
! i=1
! do while(Cprogname(i).ne.C_NULL_CHAR)
!   progname(i:i) = Cprogname(i)
!   i = i+1
! end do

! write (*,*) '->'//trim(progname)//'<-', len_trim(progname)

! stop

! t1%a = 10
! call t1%writeHDFnamelist(40)
! call t1%p(1)

! t2%a = 20
! t2%b = 30
! call t2%writeHDFnamelist(60)
! call t2%p(2)

! Timing = Timing_T( showDateTime = .TRUE. )

! call Timing % Time_tick()

! EMsoft = EMsoft_T(progname, progdesc, showconfig=.TRUE., makeconfig=.FALSE.)

! call EMsoft%setConfigParameter('WyckoffPositionsfilename','no idea where this file is located ')

! write (*,*) 'modified string = ',trim(EMsoft%getConfigParameter('WyckoffPositionsfilename'))

! Message = IO_T()

! m = 'Revision = '//trim( EMsoft % getConfigParameter('EMsoftRevision') )
! call Message % printMessage(m)

! m = 'Completed path = '//trim( EMsoft % generateFilePath('EMsofttestpath','test.h5'))
! call Message % printMessage(m)

! m = 'Completed path = '//trim( EMsoft % generateFilePath('Templatecodefilename',''))
! call Message % printMessage(m)

! io_int = (/ 10, 20 /)
! call Message % WriteValue('test message ', io_int, 2)
! call Message % WriteValue('test message ', io_int, 2, "(I10,'f----f',I10)")

! call Message % printWarning('hmmm',(/'why not ... '/))

!     Message = IO_T()

!          call Message % printWarning('EMdatapathname was not defined in the json file', &
!                      (/ 'EMDATAPATHNAME environment variable was NOT defined as a backup.', &
!                         '----> using absolute path convention                            '/) )

!     status = 999001
!     call Message % printError('EMsoftpathname was not defined in the json file', &
!                    status, (/ 'EMSOFTPATHNAME environment variable was NOT defined as a backup.' /) )

! call Message % printError('blah','more blah')

! call Message % printMessage( 'Time : '//Timing % getTimeString() )
! call Message % printMessage( 'Date : '//Timing % getDateString() )

! ! call Timing % printTimeStamp(redirect=10)

! call sleep(2)

! call Timing % Time_tick(2)

! call sleep(1)
! call Timing % Time_tock(2)
! call Timing % Time_tock(1)

! write (*,*) 'Interval 1 : ', Timing % getInterval(1)
! write (*,*) 'Interval 2 : ', Timing % getInterval(2)

! call Timing % Time_reset(2)

! write (*,*) 'Interval 1 : ', Timing % getInterval(1)
! write (*,*) 'Interval 2 : ', Timing % getInterval(2)

! write (*,*) 'calling Marsaglia generator for 10 random quaternions '

! qra = quat_randomArray(10, 'd', seed, northern=.TRUE.)
! call qra%quat_print()

! qrb = quat_randomArray(10, 's', seed, northern=.TRUE.)
! call qrb%quat_print()

end program TestProgram

max_stars_repo_name: EMsoft-org/EMsoftOO

max_issues_repo_name: EMsoft-org/EMsoftOO

max_forks_repo_name: EMsoft-org/EMsoftOO


text:
c------------------------------------------------------------
      subroutine up(txt,len)
c
c convert character string to upper case
c
      character*(*) txt
      character*80 save
      character*26 ualpha,lalpha
      data ualpha /'ABCDEFGHIJKLMNOPQRSTUVWXYZ'/
      data lalpha /'abcdefghijklmnopqrstuvwxyz'/

      do 9000 i=1,len
        if((txt(i:i).ge.'a').and.(txt(i:i).le.'z')) then
          match = index(lalpha,txt(i:i))
          save(i:i) = ualpha(match:match)
        else
          save(i:i) = txt(i:i)
        end if
c      end do
9000	continue

      txt = save
      return
      end

max_stars_repo_name: caixiuhong/Stable-MCCE

max_issues_repo_name: caixiuhong/Stable-MCCE

max_forks_repo_name: caixiuhong/Stable-MCCE


text:
C Copyright (C) 2000, International Business Machines
C Corporation and others.  All Rights Reserved.
       SUBROUTINE EVALNP(VAL, POINTS, J, POLY, I, N, LPNT, LPOLY)

C
C  *****************************************************************
C  THIS SUBROUTINE COMPUTES THE VALUE OF THE I-TH NEWTON FUNDAMENTAL
C  (STORED IN ARRAY 'POLY') AT THE J-TH POINT IN ARRAY 'POINTS'
C
C  (IF THE POLYNOMIAL IS QUADRATIC WE COMPUTE THE VALUE BY THE HORNER'S
C   RULE)
C
C  PARAMETERS
C
C  VAL    (OUTPUT) VALUE OF THE I-TH POLYNOMIAL AT THE J-TH POINT
C
C  POINTS (INPUT)  THE ARRAY OF POINTS
C
C  J      (INPUT)  THE INDEX OF THE EVALUATION POINT
C
C  N      (INPUT)  PROBLEM DIMENSION
C
C  POLY   (INPUT)  ARRAY OF NEWTON FUNDAMENTAL POLYNOMIALS
C
C  I      (INPUT)  INDEX OF THE EVALUATED POLYNOMIAL
C
C  LPNT   (INPUT)  LENGTH OF THE ARRAY 'POINTS'
C
C  LPOLY  (INPUT)  LENGTH OF THE ARRAY 'POLY'
C  *****************************************************************
C

       DOUBLE PRECISION POINTS(LPNT), POLY(LPOLY), VAL

       INTEGER          I, J, LPNT, LPOLY, N

C
C  LOCAL VARIABLES
C

       INTEGER          II, IPOINT, JJ , K,  L,  NP1

       DOUBLE PRECISION V

       NP1    = N+1

C
C  THE POINTER TO WHERE ELEMENTS OF J-TH POINT BEGIN IN ARRAY 'POINTS'
C
       IPOINT = (J-1)*N

       IF (I.EQ.1) THEN
C
C  IF THE I-TH POLYNOMIAL IS A  CONSTANT
C
         VAL = POLY(1)
       ELSE IF (I.LE.NP1) THEN
C
C  IF THE I-TH POLYNOMIAL IS LINEAR
C
         K    = ( I - 2 )*( NP1 ) + 2
         VAL  = POLY(K)
C
C  ADD UP THE TERMS
C

         DO 10 L=1,N
           VAL = VAL + POLY( K + L )*POINTS( IPOINT + L )
 10      CONTINUE
       ELSE
C
C  IF THE I-TH POLYNOMIAL IS QUADRATIC
C

         K   = ( I - NP1 - 1 )*(N+2)*NP1/2 + N*NP1 + 2
         VAL = POLY( K )
         JJ  = K + NP1
         DO 30 L=1,N
           V = POLY( K + L )
           DO 20 II=L, N
             V  = V + POLY( JJ )*POINTS( IPOINT + II )
             JJ = JJ + 1
 20        CONTINUE
           VAL=VAL+V*POINTS( IPOINT + L )
 30      CONTINUE
       ENDIF

       RETURN
       END

************************************

*         NEXT SUBROUTINE          *

************************************

       SUBROUTINE EVALX(VAL, X, POLY, I, N, LPOLY)

C
C  *****************************************************************
C  THIS SUBROUTINE COMPUTES THE VALUE OF THE I-TH NEWTON FUNDAMENTAL
C  (STORED IN ARRAY 'POLY') AT A GIVEN POINT X
C
C  IF THE POLYNOMIAL IS QUADRATIC WE COMPUTE THE VALUE BY THE HORNER'S
C  RULE
C
C  PARAMETERS
C
C  VAL    (OUTPUT) VALUE OF THE I-TH POLYNOMIAL AT THE J-TH POINT
C
C  X      (INPUT)  THE EVALUATION POINT
C
C  N      (INPUT)  PROBLEM DIMENSION
C
C  POLY   (INPUT)  ARRAY OF NEWTON FUNDAMENTAL POLYNOMIALS
C
C  I      (INPUT)  INDEX OF THE EVALUATED POLYNOMIAL
C
C  LPNT   (INPUT)  LENGTH OF THE ARRAY 'POINTS'
C
C  LPOLY  (INPUT)  LENGTH OF THE ARRAY 'POLY'
C  *****************************************************************
C

       DOUBLE PRECISION X(N), POLY(LPOLY), VAL

       INTEGER          I, LPOLY, N

C
C  LOCAL VARIABLES
C

       INTEGER          II, J , K,  L,  NP1

       DOUBLE PRECISION V
       NP1    = N+1

       IF (I.EQ.1) THEN
C
C  IF THE I-TH POLYNOMIAL IS A  CONSTANT
C
         VAL = POLY(1)
       ELSE IF (I.LE.NP1) THEN
C
C  IF THE I-TH POLYNOMIAL IS LINEAR
C

         K   = ( I - 2 )*( NP1 ) + 2
         VAL = POLY(K)
         DO 10 L=1,N
           VAL = VAL + POLY( K + L )*X( L )
 10     CONTINUE
       ELSE
C
C  IF THE I-TH POLYNOMIAL IS QUADRATIC
C
         K   = ( I - NP1 - 1 )*(N+2)*NP1/2 + N*NP1 + 2
         J   = K + NP1
         VAL = POLY( K )
         DO 30 L=1,N
           V = POLY( K + L )
           DO 20 II=L,N
             V = V + POLY( J )*X( II )
             J = J + 1
 20        CONTINUE
           VAL=VAL+V*X( L )
 30      CONTINUE
       ENDIF
       RETURN
       END

max_stars_repo_name: kant/CMMPPT

max_issues_repo_name: kant/CMMPPT

max_forks_repo_name: kant/CMMPPT


text:
!==========================================================================
pure subroutine gsw_turner_rsubrho (sa, ct, p, tu, rsubrho, p_mid)
!==========================================================================
!
!  Calculates the Turner angle and the Rsubrho as a function of pressure
!  down a vertical water column.  These quantities express the relative
!  contributions of the vertical gradients of Conservative Temperature
!  and Absolute Salinity to the vertical stability (the square of the
!  Brunt-Vaisala Frequency squared, N^2).  Tu and Rsubrho are evaluated at
!  the mid pressure between the individual data points in the vertical.
!  density in terms of SA, CT and p (IOC et al., 2010).
!
!  Note that in the double-diffusive literature, papers concerned with
!  the "diffusive" form of double-diffusive convection often define the
!  stability ratio as the reciprocal of what is defined here as the
!  stability ratio.
!
! sa      : Absolute Salinity         (a profile (length nz))     [g/kg]
! ct      : Conservative Temperature  (a profile (length nz))     [deg C]
! p       : sea pressure              (a profile (length nz))     [dbar]
! tu      : Turner angle, on the same (nz-1) grid as p_mid.
!           Turner angle has units of:           [ degrees of rotation ]
! rsubrho : Stability Ratio, on the same (nz-1) grid as p_mid.
!           rsubrho is dimensionless.                       [ unitless ]
! p_mid   : Mid pressure between p grid  (length nz-1)           [dbar]
!--------------------------------------------------------------------------

use gsw_mod_toolbox, only : gsw_specvol_alpha_beta

use gsw_mod_teos10_constants, only : rad2deg

use gsw_mod_error_functions, only : gsw_error_code

use gsw_mod_kinds

implicit none

real (r8), intent(in) :: sa(:), ct(:), p(:)
real (r8), intent(out) :: tu(:), rsubrho(:), p_mid(:)

integer :: nz, k
real (r8), dimension(:), allocatable :: dsa, sa_mid, dct, ct_mid, dp
real (r8), dimension(:), allocatable :: alpha_mid, beta_mid

character (*), parameter :: func_name = "gsw_turner_rsubrho"

nz = size(sa)
if (size(tu).lt.nz-1 .or. size(rsubrho).lt.nz-1 .or. size(p_mid).lt.nz-1) then
    tu = gsw_error_code(1,func_name)
    rsubrho = tu(1)
    p_mid = tu(1)
    return
end if

allocate (dsa(nz-1), sa_mid(nz-1), dct(nz-1), ct_mid(nz-1), dp(nz-1))
allocate (alpha_mid(nz-1), beta_mid(nz-1))

forall (k = 1: nz-1)
    dsa(k) = (sa(k) - sa(k+1))
    sa_mid(k) = 0.5_r8*(sa(k) + sa(k+1))
    dct(k) = (ct(k) - ct(k+1))
    ct_mid(k) = 0.5_r8*(ct(k) + ct(k+1))
    dp(k) = (p(k) - p(k+1))
    p_mid(k) = 0.5_r8*(p(k) + p(k+1))
end forall

call gsw_specvol_alpha_beta(sa_mid,ct_mid,p_mid(1:nz-1), &
                            alpha=alpha_mid,beta=beta_mid)

tu(1:nz-1) = rad2deg*atan2((alpha_mid*dct + beta_mid*dsa), &
                              (alpha_mid*dct - beta_mid*dsa))

where (dsa .ne. 0.0_r8)
   rsubrho = (alpha_mid*dct)/(beta_mid*dsa)
elsewhere
   rsubrho = gsw_error_code(2,func_name)
end where

deallocate (dsa, sa_mid, dct, ct_mid, dp)
deallocate (alpha_mid, beta_mid)

return
end subroutine

!--------------------------------------------------------------------------

max_stars_repo_name: minsukji/ci-debug

max_issues_repo_name: minsukji/ci-debug

max_forks_repo_name: minsukji/ci-debug


subroutine initial
!
! Subroutine reads in parameter file
! and sets up analysis
!

use tachedata
implicit none

integer :: ifile
character(100) :: suffix,eigensuffix

! Get input parameters and set up header:
  print*, " "
  print*, "-----------------------------------------------"
  print*, " TACHE: TensoriAl Classification of Hydrodynamic Elements"
  print*, "     dh4gan (current version: 1st Nov 2017)    "
  print*, "	                                       	  "
  print*, "-----------------------------------------------"
  print*, " "

  print*, " "
  print*, " "
  print*, " input parameters in ./",trim(paramfile)
  print*, "-----------------------------------------------"
  print*, " "

  open(10,file=paramfile, status='old')
  read(10,*) listfile ! File containing list of dumps to analyse
  read(10,*) filetype ! SPH, grid or mesh
  read(10,*) fileformat ! sphNG_wkmr, sphNG_iab
  read(10,*) tensorchoice ! tidal, velocity
  read(10,*) splitdumpchoice ! (y/n)
  read(10,*) threshold

  close(10)

  call checksetup

  ! Setup kernel interpolation tables if this is an SPH run
  if(filetype=='SPH') call ktable

  ! Read listfile and generate array of filenames

  open(20, file=listfile, form='formatted')

  read(20,*) nfiles

  allocate(filename(nfiles))
  do ifile=1,nfiles
     read(20,*) filename(ifile)
  enddo

  ! Generate relevant filenames:

  if(tensorchoice=='tidal') then
     tensorchar = 'T'
  else if(tensorchoice=='velocity') then
     tensorchar='V'
  endif

  allocate(gravfile(nfiles))
  allocate(potfile(nfiles))
  allocate(eigenfile(nfiles))
  allocate(vectorfile(nfiles))
  allocate(memberfile(nfiles))

  do ifile=1,nfiles

     write(suffix, '("_",A)'),trim(filename(ifile))
     write(eigensuffix,'(A1,"_",A)'),tensorchar,trim(filename(ifile))

     write(gravfile(ifile),'("gravforce",A)') trim(suffix)
     write(potfile(ifile), '("gravpotential",A)') trim(suffix)

     write(eigenfile(ifile), '("eigenvalues",A)') trim(eigensuffix)
     write(vectorfile(ifile), '("eigenvectors",A)')trim(eigensuffix)
     write(memberfile(ifile), '("tensorclass",A)') trim(eigensuffix)

  enddo

end subroutine initial


text:
      SUBROUTINE CMFIRR2(I,NNB,KLIST,XI,XID,FIRR,FD)
*
*
*       Irregular force correction from J > N (singles or unpert KS).
*       -------------------------------------------------------------
*
      INCLUDE 'common6.h'
      COMMON/PREDICT/ TPRED(NMAX)
      REAL*8  XI(3),XID(3),FIRR(3),FD(3),DX(3),DV(3),
     &        FP(6),FPD(6),FCM(3),FCMD(3),CMX(3),CMV(3),XK(6),VK(6)
      INTEGER KLIST(LMAX)
*
*
*       Set neighbour number & list index of the last single particle.
      NNB1 = NNB + 1
      NNB2 = NNB1
   25 IF (KLIST(NNB2).LE.N) GO TO 30
      NNB2 = NNB2 - 1
      IF (NNB2.GT.1) GO TO 25
*       Include special case of only c.m. neighbours.
      GO TO 40
*
*       Exit if no binaries on neighbour list.
   30 IF (NNB2.EQ.NNB1) GO TO 60

*       Perform differential correction for regularized c.m. neighbours.
   40 DO 50 LL = NNB2+1,NNB1
          J = KLIST(LL)
*       See whether to sum over binary components.
          JPAIR = J - N
          J1 = 2*JPAIR - 1
*       Skip unperturbed binary (treated as single particle).
          IF (LIST(1,J1).EQ.0) THEN
              GO TO 50
          END IF
*
*       Obtain coordinates & velocity by prediction or copying.
          IPRED = 1
          IF (TPRED(J).EQ.TIME) IPRED = 0
          CALL KSRES3(JPAIR,J1,J2,IPRED,CMX,CMV,XK,VK)
*
*       Obtain individual c.m. force with single particle approximation.
          A1 = CMX(1) - X(1,I)
          A2 = CMX(2) - X(2,I)
          A3 = CMX(3) - X(3,I)
          RIJ2 = A1*A1 + A2*A2 + A3*A3
*
          DV(1) = CMV(1) - XDOT(1,I)
          DV(2) = CMV(2) - XDOT(2,I)
          DV(3) = CMV(3) - XDOT(3,I)
          DR2I = 1.0/RIJ2
          DR3I = BODY(J)*DR2I*SQRT(DR2I)
          DRDV = 3.0*(A1*DV(1) + A2*DV(2) + A3*DV(3))*DR2I
*
          FCM(1) = A1*DR3I
          FCM(2) = A2*DR3I
          FCM(3) = A3*DR3I
          FCMD(1) = (DV(1) - A1*DRDV)*DR3I
          FCMD(2) = (DV(2) - A2*DRDV)*DR3I
          FCMD(3) = (DV(3) - A3*DRDV)*DR3I
*
*       Evaluate perturbation on first component of body #J.
          dr2 = 0.0
          drdv = 0.0
          DO 42 L = 1,3
              dx(L) = XK(L) - X(L,I)
              dv(L) = VK(L) - XDOT(L,I)
              dr2 = dr2 + dx(L)**2
              drdv = drdv + dx(L)*dv(L)
   42     CONTINUE
*
          dr2i = 1.0/dr2
          dr3i = BODY(J1)*dr2i*SQRT(dr2i)
          drdv = 3.0*drdv*dr2i
*
          DO 45 L = 1,3
              FP(L) = dx(L)*dr3i
              FPD(L) = (dv(L) - dx(L)*drdv)*dr3i
   45     CONTINUE
*
*       Evaluate perturbation on second component of body #J.
          dr2 = 0.0
          drdv = 0.0
          DO 46 L = 1,3
              dx(L) = XK(L+3) - X(L,I)
              dv(L) = VK(L+3) - XDOT(L,I)
              dr2 = dr2 + dx(L)**2
              drdv = drdv + dx(L)*dv(L)
   46     CONTINUE
*
          dr2i = 1.0/dr2
          dr3i = BODY(J2)*dr2i*SQRT(dr2i)
          drdv = 3.0*drdv*dr2i
*
          DO 48 L = 1,3
              FP(L+3) = dx(L)*dr3i
              FPD(L+3) = (dv(L) - dx(L)*drdv)*dr3i
   48     CONTINUE
*
*       Accumulate individual correction terms together.
          DO 49 L = 1,3
              FIRR(L) = FIRR(L) + (FP(L) + FP(L+3) - FCM(L))
              FD(L) = FD(L) + (FPD(L) + FPD(L+3) - FCMD(L))
   49     CONTINUE
   50 CONTINUE
*
   60 RETURN
*
      END


text:
      SUBROUTINE alltoallv_real( a, ssize, sdisp, b, rsize, rdisp )
      IMPLICIT NONE
      include 'mpif.h'
      REAL*8  a( * ), b( * )
      INTEGER ssize( * ), sdisp( * ), rsize( * ), rdisp( * )
      INTEGER ierr
      CALL MPI_ALLTOALLV( a, ssize, sdisp, MPI_DOUBLE_PRECISION,
     &        b, rsize, rdisp, MPI_DOUBLE_PRECISION,
     &        MPI_COMM_WORLD, ierr )

      RETURN
      END

max_stars_repo_name: sscerr/eTF

max_issues_repo_name: sscerr/eTF

max_forks_repo_name: sscerr/eTF


text:
      SUBROUTINE SNPIC
!     THIS SUBROUTINE WILL ANALYZE THE INPUT DATA AND PRINT A PICTURE
      character*1 ILINE(133),ICHR(50),IBLNK
      DIMENSION INUM(50)
      COMMON ISET
      DATA IBLNK/' '/
      DO 4 I=1,121
    4 ILINE(I)=IBLNK
      K=1
   10 READ (2,1000) (INUM(I),ICHR(I),I=1,ISET)
      i=1
!!      DO 40 I = 1, ISET
      do
      IF (INUM(I)+1)100,11,100
!     HERE WE WRITE A LINE TO THE PRINTER AND GO BUILD ANOTHER
   11 DO 15 L = K, 121
   15 ILINE(L)=ICHR(I)
      WRITE (3,2000) (ILINE(K),K=1,121)
      ILINE(1)=IBLNK
      DO 20 K = 2, 121
   20 ILINE(K)=ICHR(I)
      K=1
      I=I+1
  100 ITST=INUM(I)
      IF (ITST + 2) 101,200,101
  101 IF (ITST) 102,40,102
  102 DO 30 J = 1, ITST
      ILINE(K)=ICHR(I)
      K=K+1
   30 CONTINUE
   40 CONTINUE
      i=i+1
      if(i .gt. iset)goto 10
      enddo
      GOTO 10
!     HERE WE EXIT THE PICTURE AND RETURN TO THE CALLING PROGRAM
  200 RETURN
!     FORMAT STATEMENTS
 1000 FORMAT (25(I2,A1))
 2000 FORMAT (121A1)
      END
!     ****************************************************************
!     * IBM 1130 EMULATOR VERSION (KYM FARNIK)                       *
!     *                                                              *
!     * COPYRIGHT - THIS IS CONSIDERED PUBLIC DOMAIN FOR THE         *
!     * FOLLOWING REASON                                             *
!     *    DURING THE 60'S AND 70'S THIS CODE (AND VARIATIONS)       *
!     *    WAS WIDELY SHARED BETWEEN GROUPS SUCH AS SHARE/GUIDE(IBM) *
!     *    AND DECUS(DIGITAL) AND THE ORGINAL AUTHOR IS UNKNOWN AND  *
!     *    HAS NEVER ASSERTED ANY RIGHTS OVER THIS CODE.             *
!     *                                                              *
!     * THIS VERSION WAS MODIFIED FROM CODE SOURCED FROM DECUS       *
!     * SEE: http://www.ibiblio.org/pub/academic/computer-science    *
!     *          /history/pdp-11/rsts/decus/sig87/087018/            *
!     * WHICH IN TURN WAS SOURCED FROM AN IBM FORTRAN+BAL VERSION    *
!     *                                                              *
!     * IBM 1130 MODS: (FORTRAN IV SUBSET)                           *
!     *   ADDED IBM 1130 JCL                                         *
!     *   REMOVE OPEN/CLOSE                                          *
!     *   CHANGED LOGIGAL IF TO ARITHMETIC IF                        *
!     *   CHANGED TO 5 CHARACTER VARIABLE NAMES                      *
!     *   CHANGED TO 120 COLUMNS IBM 1132 PRINTER                    *
!     *   EXTENED PRECISION REAL (3*16 BIT WORDS)                    *
!     *                                                              *
!     ****************************************************************
!     * PRINTS CALENDAR, ONE MONTH PER PAGE WITH PICTURES OPTIONAL.  *
!     *                                                              *
!     * BEGINNING MONTH AND YEAR, ENDING MONTH AND YEAR MUST BE PRO- *
!     * VIDED IN 4(I6)  FORMAT ON A CARD IMMEDIATELY FOLLOWING       *
!     * CARD 98 OF DECK.                                             *
!     *                                                              *
!     * IF GRID LINES ARE DESIRED, A 1 MUST APPEAR IN COLUMN 30 OF   *
!     * ABOVE CARD.  A BLANK OR ZERO WILL SUPPRESS GRID LINES.       *
!     *                                                              *
!     * ALL PICTURE DATA DECKS MUST BE TERMINATED WITH CODE -2.      *
!     * CONSECUTIVE -2*S WILL RESULT IN NO PICTURE BEING PRINTED     *
!     * FOR THAT MONTH.                                              *
!     *                                                              *
!     * PICTURE FORMAT CODES --                                      *
!     *    -1    END OF LINE                                         *
!     *    -2    END OF PICTURE                                      *
!     *    -3    LIST CARDS, ONE PER LINE, FORMAT 13A6               *
!     *    -4    LIST CARDS, TWO PER LINE, FORMAT 11A6/11A6          *
!     *    -5    LIST CARDS, TWO PER LINE, FORMAT 12A6/10A6          *
!     ****************************************************************
      character*6 AMNTH (12,7,13), ANAM(22), ANUM(2,10,5)
      character*6 BLANK,ONE,ALIN1,ALIN2,ALIN3,ALIN4
      integer NODS(12)
      character*6 CAL(60,22)
      COMMON ISET
      READ (2,1) (((AMNTH(I,J,K),K=1,13),J=1,7),I=1,12)
      READ (2,2) (ANAM(I),I=1,22)
      READ (2,3) (((ANUM(I,J,K),J=1,10),K=1,5),I=1,2)
      READ (2,4) (NODS(I),I=1,12)
      READ (2,1) BLANK,ONE,ALIN1,ALIN2,ALIN3,ALIN4
      READ (2,4) MF,IYR,MTLST,IYLST,LNSW
1     FORMAT (13A6)
2     FORMAT (11A6)
3     FORMAT (10A6)
4     FORMAT (12I6)
5     FORMAT (A1,20A6)

      ISET=25
      DO 10 I=1,60
      DO 10 J=1,22
10    CAL(I,J)= BLANK
      CAL(1,1)= ONE
      DO 20 J=1,22
20    CAL(11,J)=ANAM(J)
      IF (LNSW) 122,142,122
122   DO 125 I=20,60,8
      DO 125 J=1,22
125   CAL(I,J)=ALIN2
      DO 140 J=4,19,3
      I=13
127   DO 130 L=1,7
      CAL(I,J)=ALIN1
130   I=I+1
      IF (I-55) 135,135,140
135   CAL(I,J)=ALIN3
      I=I+1
      GO TO 127
140   CONTINUE
      DO 141 I=20,60,8
141   CAL(I,1)=ALIN4
142   IDOW=(IYR-1751)+(IYR-1753)/4-(IYR-1701)/100+(IYR-1601)/400
      IDOW=IDOW-7*((IDOW-1)/7)
55    IF (IYR-IYLST) 60,65,100
60    ML=12
      GO TO 70
65    ML=MTLST
70    IY1=IYR/1000
      NUMB=IYR-1000*IY1
      IY2=NUMB/100
      NUMB=NUMB-100*IY2
      IY3=NUMB/10
      NUMB=NUMB-10*IY3
      IY4=NUMB
      DO 72 J=1,5
      CAL(J+3,2)=ANUM(2,IY1+1,J)
      CAL(J+1,3)=ANUM(2,IY2+1,J)
      CAL(J+1,19)=ANUM(2,IY3+1,J)
72    CAL(J+3,20)=ANUM(2,IY4+1,J)
      LPYSW=0
      IF (IYR-4*(IYR/4)) 90,75,90
75    IF (IYR-100*(IYR/100)) 85,80,85
80    IF (IYR-400*(IYR/400)) 90,85,90
85    LPYSW=1
90    NODS(2)=NODS(2)+LPYSW
      IF (MF-1) 100,110,95
95    MF=MF-1
      DO 105 MONTH=1,MF
105   IDOW=IDOW+NODS(MONTH)
      IDOW=IDOW-7*((IDOW-1)/7)
      MF=MF+1
110   DO 51 MONTH=MF,ML
      LSTDY=NODS(MONTH)
      DO 115 I=1,7
      DO 115 JM=1,13
      J=JM+4
115   CAL(I,J)=AMNTH(MONTH,I,JM)
      IF (IDOW-1) 160,160,120
120   ID=IDOW-1
      J=2
      DO 155 K=1,ID
      DO 150 I=14,18
      CAL (I,J)= BLANK
150   CAL(I,J+1)= BLANK
      J=J+3
155   CONTINUE
160   IDAY=1
      II=14
25    J=3*IDOW-1
      N=IDAY/10+1
      I=II
      DO 30 K=1,5
      CAL(I,J)=ANUM(1,N,K)
30    I=I+1
      N=IDAY-10*N+11
      J=J+1
      I=II
      DO 35 K=1,5
      CAL(I,J)=ANUM(2,N,K)
35    I=I+1
      IDOW=IDOW+1
      IF (IDOW-7) 45,45,40
40    IDOW=1
      II=II+8
45    IDAY=IDAY+1
      IF (IDAY-LSTDY) 25,25,50
50    ID=IDOW
205   I=II
      J=3*ID-1
      DO 210 K=1,5
      CAL(I,J)= BLANK
      CAL(I,J+1)= BLANK
210   I=I+1
      IF (ID-7) 215,220,220
215   ID=ID+1
      GO TO 205
220   IF (II-54) 225,230,230
225   II=54
      ID=1
      GO TO 205
230   CALL SNPIC
      WRITE (3,5) ((CAL(I,J),J=1,21),I=1,60)
51    CONTINUE
      IF (IYR-IYLST) 235,100,100
235   NODS(2)=NODS(2)-LPYSW
      IYR=IYR+1
      MF=1
      GO TO 55
100   CALL EXIT
      END

max_stars_repo_name: timkphd/examples

max_issues_repo_name: timkphd/examples

max_forks_repo_name: timkphd/examples


text:
c
c    This program finds the irreducible fragment that can be used
c    to generate a molecule given the symmetry point group.
c
      implicit real*8 (a-h,o-z)
      parameter(maxops=192,max_cent=2000)
      logical oprint
      dimension cdist(3),cang(3)
      double precision sym_ops(maxops*3,4)
      double precision coords(3,max_cent)
      double precision coords_new(3,max_cent)
      character*10 sym_molgnames(46)
      character*10 group_name
      data sym_molgnames/
     & 'C1','CS','CI','C2','C3',
     & 'C4','C5','C6','C7','C8',
     & 'D2','D3','D4','D5','D6',
     & 'C2V','C3V','C4V','C5V','C6V',
     & 'C2H','C3H','C4H','C5H','C6H',
     & 'D2H','D3H','D4H','D5H','D6H',
     & 'D8H','D2D','D3D','D4D','D5D',
     & 'D6D','S4','S6','S8','T',
     & 'TH','TD','O','OH','I',
     & 'IH'/
c
c  read in information from standard input.
c
      read(5,*) nata
      if (nata.gt.max_cent) call exit(2)
      read(5,*) threquiv
      read(5,100) group_name
  100 format(a10)
      do i = 1, nata
        read(5,*) coords(1,i),coords(2,i),coords(3,i)
      end do
c
c  parameters used in gensym and dosymops calls. itype = 0 for
c  molecules, itype = 1,2,3 for 1d,2d,3d periodic systems (not
c  currently implemented). For molecular systems, numset should
c  probably be 0 or 1 (probably 1).
c
      itype  = 0
      numset = 1
      do i = 1, 46
       if (group_name.eq.sym_molgnames(i)) numgrp = i
      end do
      oprint = .false.
c
c   lattice vectors & angles (cell constants)
c   values should have been read in and initialized properly for
c   each system dimension (not currently implemented)
c
      do 1000 i=1,3
        cdist(i)= 0.0d00
        cang(i) = 0.0d00
 1000 continue
c
c    generate symmetry operators for whatever group (crystals, surfaces
c    polymers or molecules) have been requested.
c
      call gensym(itype,numgrp,numset,sym_ops,nops,oprint,group_name)
c
c    find irreducible fragment
c
      call sym_avg(sym_ops,nops,coords,nata,
     &     itype,coords_new,cdist,cang,threquiv,ierr)
c
c
      write(6,300) nata
  300 format(i8)
      do i = 1, nata
        write(6,400) (coords_new(j,i),j=1,3)
  400   format(3e24.14)
      end do
c
c   We are enforcing error checking by having the
c   code exit with a non-zero status using the exit()
c   subroutine call. This is a non-standard call and may
c   conflict with the F90 EXIT statement if this code is
c   ever compiled under Fortran 90.
c
      if (ierr.ne.0) call exit(ierr)
      stop
      end

max_stars_repo_name: mattbernst/ECCE

max_issues_repo_name: mattbernst/ECCE

max_forks_repo_name: mattbernst/ECCE


text:
module modis28

	use CKD_common

	implicit none
	private

	integer, parameter :: ni = 5
	integer, parameter :: nch4 = 5

	public :: ch28, init_ch28

contains
	subroutine init_ch28(modis_channel)

		type(ckd_type), intent(inout) :: modis_channel
		integer :: all_weights

		all_weights = ni*nch4

		modis_channel%num_weights = all_weights
		allocate(modis_channel%weights(all_weights))
		allocate(modis_channel%taus(all_weights, max_layers))

	end subroutine init_ch28

      subroutine ch28(u0, p0,t0, ux, fac,modis_channel)
!
!      This routine was developed for MODIS Channel 28: 1335--1395 cm^{-1}
!
!-----------------------------------------------------------------
!     Channel 28   7.175 - 7.475 microns
!-----------------------------------------------------------------
!  INPUTS:
!     u0 --> water vapor amount [g/cm2/km]
!     p0 --> pressure in atmospheres
!     t0 --> temp in K
!	  ux --> ozone (g/cm^2/km)
!     fac --> layer thickness in km
!-----------------------------------------------------------------

		real, dimension(:), intent(in) :: u0, p0, t0, ux, fac
		type(ckd_type), intent(inout) :: modis_channel

		integer, parameter :: nlev = max_levels
		integer, parameter :: nlay = max_layers
		integer, parameter :: mlv = nlay

		real :: p(nlay), t(nlay),u(nlay), xlog, &
				  uch4(nlay)

		real :: tau(ni,nlay),tauch4(nch4, nlay), &
				 f(ni),  fch4(nch4), twin

		integer :: i,m,ia, index

      do  m=1,use_layers
        t(m)=(t0(m)+t0(m+1))/2.0
        p(m)=(p0(m)+p0(m+1))/2.0
!       u(m)=fac(m)*(u0(m)+u0(m+1))/2.0*6.023e+23/2.687e+19/18.01534
        u(m)=fac(m)*(u0(m)+u0(m+1))/2.0
        uch4(m)=p(m)*1.01325e+06/(1.3805e-16*t(m))*1.75e-06*fac(m)* &
     &       1.0e+05/6.023e+23*16.04303
      enddo

      call ck_28(u,f,p,t,tau)
      call ckch4_28(uch4,fch4,p,t,tauch4)

      do  m=1,use_layers
        call window_28(u(m),twin,p(m),t(m),fac(m))
          do ia=1,nch4 ! 5 k's for ch4
        do i=1,ni  ! 5 k's for h2o
			index = (ia-1)*ni + i
            modis_channel%taus(index,m)=tau(i,m)+twin+tauch4(ia,m)
 			if (m == 1) modis_channel%weights(index) = f(i)*fch4(ia)
          end do
        end do
      enddo

	end subroutine ch28

! *********************************************************************
!     h2o 1335--1395 cm^{-1}
      subroutine ck_28(u,f,p,t,tau)
		real, dimension(:), intent(in) :: p, t, u

		real, dimension(:), intent(inout) :: f
		real, dimension(:,:), intent(inout) :: tau

		real :: k(ni)
		real :: coefk(ni,3,num_pressures)
		integer :: i, jp, jt
      f(1)=0.526064
      f(2)=0.309565
      f(3)=0.117418
      f(4)=0.041242
      f(5)=0.005711
      k(1)=1.7591042
      do i=2,ni
        k(i)=8.0*k(i-1)
      end do
      data ( (coefk(1,jt,jp), jp = 1, 19), jt = 1, 3)/ &
     & 0.00062,0.00098,0.00154,0.00242,0.00381,0.00599, &
     & 0.00939,0.01471,0.02305,0.03606,0.05636,0.08828, &
     & 0.13929,0.22014,0.34626,0.53718,0.81621,1.22661, &
     & 1.80673, &
     & 6.500E-06,1.025E-05,1.613E-05,2.538E-05,3.963E-05,6.213E-05, &
     & 9.675E-05,1.514E-04,2.360E-04,3.688E-04,5.795E-04,9.074E-04, &
     & 1.431E-03,2.270E-03,3.581E-03,5.594E-03,8.534E-03,1.287E-02, &
     & 1.904E-02, &
     & 1.250E-08,1.875E-08,3.438E-08,5.938E-08,8.438E-08,1.281E-07, &
     & 2.000E-07,3.219E-07,4.875E-07,7.437E-07,1.212E-06,1.916E-06, &
     & 3.000E-06,4.819E-06,7.669E-06,1.271E-05,2.050E-05,3.285E-05, &
     & 4.776E-05/
      data ( (coefk(2,jt,jp), jp = 1, 19), jt = 1, 3)/ &
     & 0.00096,0.00147,0.00226,0.00347,0.00534,0.00825, &
     & 0.01276,0.01977,0.03076,0.04771,0.07370,0.11322, &
     & 0.17291,0.25854,0.37848,0.55473,0.82280,1.21548, &
     & 1.74685, &
     & 1.088E-05,1.637E-05,2.487E-05,3.762E-05,5.725E-05,8.800E-05, &
     & 1.355E-04,2.096E-04,3.231E-04,4.994E-04,7.725E-04,1.182E-03, &
     & 1.806E-03,2.692E-03,3.979E-03,5.992E-03,9.123E-03,1.369E-02, &
     & 1.967E-02, &
     & 3.438E-08,4.687E-08,5.937E-08,8.437E-08,1.312E-07,2.000E-07, &
     & 3.062E-07,4.781E-07,6.719E-07,9.906E-07,1.562E-06,2.372E-06, &
     & 3.675E-06,5.878E-06,9.709E-06,1.389E-05,1.807E-05,2.477E-05, &
     & 3.312E-05/
      data ( (coefk(3,jt,jp), jp = 1, 19), jt = 1, 3)/ &
     & 0.00196,0.00281,0.00407,0.00594,0.00870,0.01285, &
     & 0.01912,0.02843,0.04214,0.06186,0.08957,0.12728, &
     & 0.18341,0.26963,0.39774,0.58307,0.84273,1.16653, &
     & 1.47341, &
     & 2.925E-05,3.950E-05,5.437E-05,7.587E-05,1.071E-04,1.531E-04, &
     & 2.199E-04,3.173E-04,4.546E-04,6.528E-04,9.160E-04,1.298E-03, &
     & 1.931E-03,2.918E-03,4.393E-03,6.534E-03,9.599E-03,1.356E-02, &
     & 1.745E-02, &
     & 1.562E-07,1.812E-07,2.156E-07,2.531E-07,3.281E-07,4.469E-07, &
     & 5.531E-07,7.375E-07,9.406E-07,1.344E-06,1.763E-06,2.972E-06, &
     & 4.247E-06,5.406E-06,7.066E-06,1.065E-05,1.552E-05,2.320E-05, &
     & 3.507E-05/
      data ( (coefk(4,jt,jp), jp = 1, 19), jt = 1, 3)/ &
     & 0.02474,0.02716,0.03081,0.03609,0.04349,0.05438, &
     & 0.06979,0.09164,0.12349,0.17027,0.23742,0.33297, &
     & 0.45910,0.60364,0.75125,0.88178,0.97072,1.02249, &
     & 1.00497, &
     & 4.100E-04,4.280E-04,4.606E-04,5.182E-04,6.022E-04,7.234E-04, &
     & 9.085E-04,1.172E-03,1.547E-03,2.079E-03,2.894E-03,4.029E-03, &
     & 5.545E-03,7.384E-03,9.206E-03,1.095E-02,1.171E-02,1.222E-02, &
     & 1.193E-02, &
     & 3.275E-06,3.306E-06,3.422E-06,3.612E-06,3.956E-06,4.078E-06, &
     & 4.463E-06,5.184E-06,5.691E-06,5.903E-06,7.572E-06,9.675E-06, &
     & 1.257E-05,1.655E-05,2.163E-05,2.821E-05,2.849E-05,2.838E-05, &
     & 2.990E-05/
      data ( (coefk(5,jt,jp), jp = 1, 19), jt = 1, 3)/ &
     & 2.47130,2.46377,2.45640,2.44544,2.42681,2.40072, &
     & 2.36233,2.31177,2.24881,2.16870,2.07667,1.96404, &
     & 1.82806,1.66439,1.48253,1.29109,1.09475,0.90080, &
     & 0.71857, &
     & 2.732E-02,2.725E-02,2.721E-02,2.712E-02,2.697E-02,2.681E-02, &
     & 2.645E-02,2.605E-02,2.544E-02,2.467E-02,2.370E-02,2.244E-02, &
     & 2.090E-02,1.900E-02,1.687E-02,1.450E-02,1.260E-02,1.029E-02, &
     & 7.818E-03, &
     & 4.491E-05,4.539E-05,4.539E-05,4.364E-05,4.401E-05,4.425E-05, &
     & 4.409E-05,4.453E-05,4.300E-05,4.284E-05,4.264E-05,4.075E-05, &
     & 3.800E-05,3.592E-05,3.070E-05,2.342E-05,2.162E-05,1.951E-05, &
     & 1.032E-05/

		call apply_ckd(ni, k, u, f, p, t, tau, coefk, .false.)

      end subroutine ck_28
! *********************************************************************
!     ch4 1335--1395 cm^{-1}
      subroutine ckch4_28(u,f,p,t,tau)
		real, dimension(:), intent(in) :: p, t, u

		real, dimension(:), intent(inout) :: f
		real, dimension(:,:), intent(inout) :: tau

		real :: k(nch4)
		real :: coefk(nch4,3,num_pressures)
		integer :: i, jp, jt
      f(1)=0.426433
      f(2)=0.340613
      f(3)=0.166590
      f(4)=0.065309
      f(5)=0.001055
      k(1)=15.6509586
      do i=2,nch4
        k(i)=8.0*k(i-1)
      end do
      data ( (coefk(1,jt,jp), jp = 1, 19), jt = 1, 3)/ &
     & 0.00000,0.00000,0.00000,0.00000,0.00008,0.00029, &
     & 0.00082,0.00193,0.00425,0.00898,0.01875,0.04152, &
     & 0.08782,0.16846,0.31236,0.52235,0.80696,1.21861, &
     & 1.73129, &
     &  0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 1.250E-06, 4.750E-06, &
     &  1.363E-05, 3.175E-05, 6.413E-05, 1.294E-04, 2.410E-04, 4.560E-04, &
     &  8.580E-04, 1.494E-03, 2.569E-03, 4.527E-03, 7.177E-03, 1.153E-02, &
     &  1.783E-02, &
     &  0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 1.137E-15, 1.875E-08, &
     &  5.938E-08, 1.688E-07, 3.281E-07, 6.781E-07, 1.225E-06, 1.594E-06, &
     &  2.006E-06, 4.303E-06, 1.306E-06,-1.872E-06,-3.906E-06,-1.012E-05, &
     & -9.453E-06/
      data ( (coefk(2,jt,jp), jp = 1, 19), jt = 1, 3)/ &
     & 0.00027,0.00029,0.00033,0.00059,0.00108,0.00212, &
     & 0.00423,0.00838,0.01613,0.02932,0.05124,0.08652, &
     & 0.14456,0.23771,0.37767,0.57391,0.83532,1.19534, &
     & 1.66934, &
     &  1.062E-05, 1.138E-05, 1.200E-05, 1.562E-05, 2.362E-05, 3.738E-05, &
     &  6.275E-05, 1.058E-04, 1.706E-04, 2.501E-04, 3.873E-04, 5.951E-04, &
     &  9.106E-04, 1.358E-03, 1.935E-03, 2.680E-03, 3.623E-03, 4.858E-03, &
     &  6.160E-03, &
     &  1.281E-07, 1.406E-07, 1.500E-07, 1.469E-07, 1.969E-07, 2.594E-07, &
     &  3.813E-07, 5.125E-07, 6.094E-07, 8.094E-07, 9.125E-07, 1.116E-06, &
     &  1.478E-06, 2.119E-06, 2.588E-06, 1.216E-06, 5.719E-07,-1.653E-06, &
     & -6.644E-06/
      data ( (coefk(3,jt,jp), jp = 1, 19), jt = 1, 3)/ &
     & 0.00245,0.00255,0.00279,0.00408,0.00655,0.01082, &
     & 0.01690,0.02530,0.03707,0.05464,0.08205,0.12517, &
     & 0.19330,0.29851,0.44980,0.64644,0.87803,1.11719, &
     & 1.33794, &
     &  6.575E-05, 6.825E-05, 6.750E-05, 7.600E-05, 9.987E-05, 1.271E-04, &
     &  1.564E-04, 1.872E-04, 2.242E-04, 2.740E-04, 3.414E-04, 4.508E-04, &
     &  6.183E-04, 8.915E-04, 1.301E-03, 1.913E-03, 2.686E-03, 3.347E-03, &
     &  3.775E-03, &
     &  5.750E-07, 6.125E-07, 6.500E-07, 5.375E-07, 6.094E-07, 6.094E-07, &
     &  6.469E-07, 7.313E-07, 8.750E-07, 9.188E-07, 7.344E-07, 3.000E-07, &
     & -6.187E-07,-1.894E-06,-4.491E-06,-6.503E-06,-8.113E-06,-7.216E-06, &
     & -9.734E-06/
      data ( (coefk(4,jt,jp), jp = 1, 19), jt = 1, 3)/ &
     & 0.47004,0.47577,0.48301,0.49214,0.50136,0.51167, &
     & 0.52363,0.54155,0.57086,0.61648,0.68084,0.76979, &
     & 0.87775,0.98090,1.04569,1.05977,1.02796,0.96499, &
     & 0.87486, &
     &  3.241E-03, 3.261E-03, 3.252E-03, 3.223E-03, 3.161E-03, 3.051E-03, &
     &  2.861E-03, 2.616E-03, 2.324E-03, 2.063E-03, 1.916E-03, 1.867E-03, &
     &  1.926E-03, 2.047E-03, 2.058E-03, 1.941E-03, 1.730E-03, 1.466E-03, &
     &  1.213E-03, &
     &  3.438E-08,-1.562E-08, 9.313E-12,-5.969E-07,-5.719E-07,-6.656E-07, &
     & -8.813E-07,-6.219E-07,-6.031E-07,-1.522E-06,-2.325E-06,-3.216E-06, &
     & -4.716E-06,-6.494E-06,-6.875E-06,-6.991E-06,-7.169E-06,-8.038E-06, &
     & -6.550E-06/
      data ( (coefk(5,jt,jp), jp = 1, 19), jt = 1, 3)/ &
     & 8.90708,8.95982,9.00443,9.01234,8.93422,8.74311, &
     & 8.39020,7.87647,7.20785,6.40016,5.57208,4.67567, &
     & 3.78880,2.98249,2.26846,1.65120,1.18555,0.83498, &
     & 0.58795, &
     & -4.827E-03,-4.787E-03,-5.012E-03,-5.177E-03,-5.327E-03,-5.256E-03, &
     & -5.075E-03,-4.295E-03,-2.663E-03, 1.191E-04, 2.788E-03, 5.592E-03, &
     &  6.374E-03, 5.920E-03, 5.053E-03, 4.087E-03, 2.991E-03, 2.124E-03, &
     &  1.317E-03, &
     & -7.298E-05,-7.411E-05,-7.494E-05,-8.078E-05,-7.758E-05,-7.528E-05, &
     & -7.372E-05,-6.625E-05,-6.545E-05,-6.135E-05,-6.455E-05,-6.573E-05, &
     & -4.441E-05,-2.553E-05,-2.307E-05,-1.565E-05,-1.322E-05,-9.872E-06, &
     & -5.906E-06/
		call apply_ckd(nch4, k, u, f, p, t, tau, coefk, .false.)

      end subroutine ckch4_28
! *********************************************************************
      subroutine window_28 (amnt,twin,patm,temp,fac)
!     Parameterization of CKD_2.4 continuum over 1335 to 1395 cm-1 band
!     Fit inaccurate for amnt < 0.01 g/cm*2;
!          however optical depths < 0.001.
!     INPUT:
!     amnt = h2O layer amount (g/cm**2)
!     patm= pressure (atm)
!     temp = temperature (k)
!     OUTPUT:
!     twin = parameterized CKD_2.4 optical depth over band
!
	  real, intent(in) :: amnt, patm, temp, fac
	  real, intent(inout) :: twin

      integer, parameter :: ncoef=7

      real ::  aa(ncoef), ph2o, tau_log
!
!
!     data aa/0.8065,1.019,-0.003620,1.0180,15.940,-8.057e-03,-0.03242/
      data aa/0.0000,1.200,-0.003620,1.0180,15.940,-8.057e-03,-0.03242/
!
      ph2o = amnt*(8.314e+07*temp)/(fac*1.0e+05*18.01534*1.01325e+06)
!
          tau_log   = aa(1)              +  &
     &                 aa(2) * log(amnt)  + &
     &                 aa(3) * temp       + &
     &                 aa(4) * patm       + &
     &                 aa(5) * ph2o       + &
     &                 aa(6) * amnt       + &
     &                 aa(7) * log(ph2o)

        twin = exp ( tau_log )
!
        end subroutine window_28

end module modis28

max_stars_repo_name: GEOS-ESM/AeroApps

max_issues_repo_name: GEOS-ESM/AeroApps

max_forks_repo_name: GEOS-ESM/AeroApps


text:
Address(Ottowa Avenue) is a residential street in West Davis.

 Intersecting Streets

  Garrison Street
  Rubicon Avenue
  Seine Avenue

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
program ferecovery

  use boundary
  use cubature
  use element
  use fsystem
  use genoutput
  use hadaptaux
  use hadaptivity
  use io
  use linearsystemblock
  use linearsystemscalar
  use pprocgradients
  use spatialdiscretisation
  use storage
  use triangulation
  use vectorio

  implicit none

  type(t_blockDiscretisation) :: rblockDiscrSolution, rblockDiscrGradient
  type(t_boundary) :: rbnd
  type(t_triangulation) :: rtria
  type(t_vectorBlock) :: rsolution, rgradient
  real(DP), dimension(:), pointer :: p_Ddata
  character(len=256) :: sarg,scubature,selement,sgrad,smesh,sname,ssol,stoken
  integer ::  i,iarg,iel,ilinelen,ios,itoken,iunit,ntoken,ndim
  integer(I32), dimension(:), allocatable :: Celement
  logical :: bbnd

  ! Initialise system-wide settings:
  call sys_init()

  ! Initialise the output system.
  call output_init()

  ! Initialise the FEAT 2.0 storage management:
  call storage_init(999, 100)

  ! Print help
  if(sys_ncommandLineArgs() .lt. 1) then
    call output_lbrk()
    call output_line("USAGE: ferecovery <options>")
    call output_lbrk()
    call output_line("Valid options:")
    call output_line("-read1d <mesh>              Read 1D mesh from <mesh>.tri")
    call output_line("-read2d <mesh>              Read 2D mesh from <mesh>.tri/prm")
    call output_line("-read3d <mesh>              Read 3D mesh from <mesh>.tri")
    call output_line("-sol <file>                 Read nodal solution values from <file>")
    call output_line("-grad <file>                Write nodal gradient values to <file>")
    call output_line("-elem <elementID>           Specify type(s) of element by <elementID>")
    call output_lbrk()
    call sys_halt()
  end if

  ! Initialise standard parameters
  ndim = 0
  bbnd = .false.
  smesh = ''
  sgrad = ''
  ssol = ''

  ! Loop over all arguments
  iarg = 1
  do while(iarg .lt. sys_ncommandLineArgs())
    ! fetch argument string
    call sys_getcommandLineArg(iarg, sarg)
    iarg = iarg + 1

    ! check argument
    if(sarg .eq. '-read1d') then
      call sys_getcommandLineArg(iarg, smesh)
      iarg = iarg + 1
      ndim = 1
    else if(sarg .eq. '-read2d') then
      call sys_getcommandLineArg(iarg, smesh)
      iarg = iarg + 1
      ndim = 2
    else if(sarg .eq. '-read3d') then
      call sys_getcommandLineArg(iarg, smesh)
      iarg = iarg + 1
      ndim = 3
    else if(sarg .eq. '-sol') then
      call sys_getcommandLineArg(iarg, ssol)
      iarg = iarg + 1
    else if(sarg .eq. '-grad') then
      call sys_getcommandLineArg(iarg, sgrad)
      iarg = iarg + 1
    else if(sarg .eq. '-elem') then
      call sys_getcommandLineArg(iarg, selement)
      call sys_countTokens(selement, ntoken, ',')
      allocate(Celement(ntoken))
      i = 1
      do itoken = 1,ntoken
        call sys_getNextToken(selement, stoken, i, ',')
        Celement(itoken) = elem_igetID(trim(stoken))
      end do
      iarg = iarg + 1
    else
      ! unknown parameter
      call output_line("ERROR: unknown parameter '" //trim(sarg)//"'")
      call sys_halt()
    end if
  end do

  ! Read boundary and mesh
  select case(ndim)
  case (1)
    call output_line("Reading mesh from '"//"./"//trim(smesh)//".tri'...")
    call tria_readTriFile1D (rtria, "/."//trim(smesh)//'.tri')
    bbnd = .false.
  case (2)
    inquire(file="./"//trim(smesh)//'.prm', exist=bbnd)
    if (bbnd) then
      call output_line("Reading mesh from '"//"./"//trim(smesh)//".tri/prm'...")
      call boundary_read_prm(rbnd, "./"//trim(smesh)//'.prm')
      call tria_readTriFile2D (rtria, "./"//trim(smesh)//'.tri', rbnd)
    else
      call output_line("Reading mesh from '"//"./"//trim(smesh)//".tri'...")
      call tria_readTriFile2D (rtria, "./"//trim(smesh)//'.tri')
    end if
  case (3)
    call output_line("Reading mesh from '"//"./"//trim(smesh)//".tri'...")
    call tria_readTriFile3D (rtria, "./"//trim(smesh)//'.tri')
    bbnd = .false.
  case default
    call output_line("ERROR: no input mesh specified")
    call sys_halt()
  end select

  ! Initialise block and spatial discretisations
  if (bbnd) then
    call spdiscr_initBlockDiscr(rblockDiscrSolution, 1, rtria, rbnd)
    call spdiscr_initBlockDiscr(rblockDiscrGradient, ndim, rtria, rbnd)
    select case(ndim)
    case (1)
      call spdiscr_initDiscr_simple(rblockDiscrSolution%RspatialDiscr(1),&
          Celement(1), rtria, rbnd)
      call spdiscr_initDiscr_simple(rblockDiscrGradient%RspatialDiscr(1),&
          Celement(1), rtria, rbnd)
    case (2)
      if (size(Celement,1) .eq. 1) then
        call spdiscr_initDiscr_simple(rblockDiscrSolution%RspatialDiscr(1),&
            Celement(1), rtria, rbnd)
        call spdiscr_initDiscr_simple(rblockDiscrGradient%RspatialDiscr(1),&
            Celement(1), rtria, rbnd)
        call spdiscr_initDiscr_simple(rblockDiscrGradient%RspatialDiscr(2),&
            Celement(1), rtria, rbnd)
      else
        call spdiscr_initDiscr_simple(rblockDiscrSolution%RspatialDiscr(1),&
            Celement(1), Celement(2), rtria, rbnd)
        call spdiscr_initDiscr_simple(rblockDiscrGradient%RspatialDiscr(1),&
            Celement(1), Celement(2), rtria, rbnd)
        call spdiscr_initDiscr_simple(rblockDiscrGradient%RspatialDiscr(2),&
            Celement(1), Celement(2), rtria, rbnd)
      end if
    case (3)
      call spdiscr_initDiscr_simple(rblockDiscrSolution%RspatialDiscr(1),&
          Celement(1), rtria, rbnd)
      call spdiscr_initDiscr_simple(rblockDiscrGradient%RspatialDiscr(1),&
          Celement(1), rtria, rbnd)
      call spdiscr_initDiscr_simple(rblockDiscrGradient%RspatialDiscr(2),&
          Celement(1), rtria, rbnd)
      call spdiscr_initDiscr_simple(rblockDiscrGradient%RspatialDiscr(3),&
          Celement(1), rtria, rbnd)
    end select
  else
    call spdiscr_initBlockDiscr(rblockDiscrSolution, 1, rtria)
    call spdiscr_initBlockDiscr(rblockDiscrGradient, ndim, rtria)
    select case(ndim)
    case (1)
      call spdiscr_initDiscr_simple(rblockDiscrSolution%RspatialDiscr(1),&
          Celement(1), rtria)
      call spdiscr_initDiscr_simple(rblockDiscrGradient%RspatialDiscr(1),&
          Celement(1), rtria)
    case (2)
      if (size(Celement,1) .eq. 1) then
        call spdiscr_initDiscr_simple(rblockDiscrSolution%RspatialDiscr(1),&
            Celement(1), rtria)
        call spdiscr_initDiscr_simple(rblockDiscrGradient%RspatialDiscr(1),&
            Celement(1), rtria)
        call spdiscr_initDiscr_simple(rblockDiscrGradient%RspatialDiscr(2),&
            Celement(1), rtria)
      else
        call spdiscr_initDiscr_simple(rblockDiscrSolution%RspatialDiscr(1),&
            Celement(1), Celement(2), rtria)
        call spdiscr_initDiscr_simple(rblockDiscrGradient%RspatialDiscr(1),&
            Celement(1), Celement(2), rtria)
        call spdiscr_initDiscr_simple(rblockDiscrGradient%RspatialDiscr(2),&
            Celement(1), Celement(2), rtria)
      end if
    case (3)
      call spdiscr_initDiscr_simple(rblockDiscrSolution%RspatialDiscr(1),&
          Celement(1), rtria)
      call spdiscr_initDiscr_simple(rblockDiscrGradient%RspatialDiscr(1),&
          Celement(1), rtria)
      call spdiscr_initDiscr_simple(rblockDiscrGradient%RspatialDiscr(2),&
          Celement(1), rtria)
      call spdiscr_initDiscr_simple(rblockDiscrGradient%RspatialDiscr(3),&
          Celement(1), rtria)
    end select
  end if

  ! Read solution from file
  call lsysbl_createVector(rblockDiscrSolution, rsolution, .false.)
  call output_line("Reading solution from '"//"./"//trim(ssol)//"'...")
  call vecio_readBlockVectorHR(rsolution, sname, .true., 0, "./"//trim(ssol), .true.)

  ! Recover gradient by superconvergent patch recovery
  call lsysbl_createVector(rblockDiscrGradient, rgradient, .false.)
  call ppgrd_calcGradient(rsolution%RvectorBlock(1), rgradient,&
      PPGRD_ZZTECHNIQUE, PPGRD_NODEPATCH)

  ! Write gradient to file
  call output_line("Writing gradient to '"//"./"//trim(sgrad)//"'...")
  call vecio_writeBlockVectorHR(rgradient, 'grad', .true., 0, "./"//trim(sgrad), '(E20.10)')

  ! Clean up
  deallocate(Celement)
  call lsysbl_releaseVector(rsolution)
  call lsysbl_releaseVector(rgradient)
  call spdiscr_releaseBlockDiscr(rblockDiscrSolution)
  call spdiscr_releaseBlockDiscr(rblockDiscrGradient)
  call tria_done(rtria)
  if(bbnd) call boundary_release(rbnd)

  ! Clean up the storage management, finish
  call storage_done()

end program ferecovery

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
!> \brief  This module provides function to calcualte to do updates
module m_fstr_Update
  use m_fstr
  implicit none

  private :: Update_abort

contains

  !=====================================================================*
  !  UPDATE_C3
  !>  \brief 変位／応力・ひずみ／内力のアップデート
  !>  \author     K. Sato(Advancesoft), X. YUAN(AdavanceSoft)
  !>  \version    0.00
  !> \if AS
  !>    \par  サブルーチン構成
  !>    -# 変位の更新              \f$ u_{n+1}^{(k)} = u_{n+1}^{(k-1)} + \delta u^{(k)} \f$
  !>    -# ひずみ・応力の更新       \f$ \varepsilon_{n+1}^{(k)} = \varepsilon_{n+1}^{(k-1)} + \delta \varepsilon^{(k)} \f$, \f$ \sigma_{n+1}^{(k)} = \sigma_{n+1}^{(k-1)} + \delta \sigma^{(k)} \f$
  !>    -# 内力（等価節点力）の計算 \f$ Q_{n+1}^{(k-1)} ( u_{n+1}^{(k-1)} ) \f$
  !> \endif
  subroutine fstr_UpdateNewton ( hecMESH, hecMAT, fstrSOLID, time, tincr,iter, strainEnergy)
    !=====================================================================*
    use m_static_lib

    type (hecmwST_matrix)       :: hecMAT    !< linear equation, its right side modified here
    type (hecmwST_local_mesh)   :: hecMESH   !< mesh information
    type (fstr_solid)           :: fstrSOLID !< we need boundary conditions of curr step
    real(kind=kreal),intent(in) :: time      !< current time
    real(kind=kreal),intent(in) :: tincr     !< time increment
    integer, intent(in)         :: iter      !< NR iterations

    integer(kind=kint) :: nodLOCAL(20)
    real(kind=kreal)   :: ecoord(3, 20), stiff(60, 60)
    real(kind=kreal)   :: thick
    integer(kind=kint) :: ndof, itype, is, iE, ic_type, nn, icel, iiS, i, j

    real(kind=kreal)   :: total_disp(6, 20), du(6, 20), ddu(6, 20)
    real(kind=kreal)   :: tt(20), tt0(20), ttn(20), qf(20*6), coords(3, 3)
    integer            :: ig0, ig, ik, in, ierror, isect, ihead, cdsys_ID
    integer            :: ndim, initt

    real(kind=kreal), optional :: strainEnergy
    real(kind=kreal) :: tmp
    real(kind=kreal)   :: ddaux(3,3)

    ndof = hecMAT%NDOF
    fstrSOLID%QFORCE=0.0d0

    tt0 = 0.d0
    ttn = 0.d0
    tt = 0.d0

    ! if initial temperature exists
    initt = 0
    if( associated(g_InitialCnd) ) then
        do j=1,size(g_InitialCnd)
          if( g_InitialCnd(j)%cond_name=="temperature" ) then
            initt=j
            exit
          endif
        end do
    endif

    ! --------------------------------------------------------------------
    !      updated
    !   1. stress and strain  : ep^(k) = ep^(k-1)+dep^(k)
    !                           sgm^(k) = sgm^(k-1)+dsgm^(k)
    !   2. Internal Force     : Q^(k-1) ( u^(k-1) )
    ! --------------------------------------------------------------------
    ! ----------------------------------------------------------------------------------
    !      calculate the Strain and Stress and Internal Force ( Equivalent Nodal Force )
    ! ----------------------------------------------------------------------------------

    do itype = 1, hecMESH%n_elem_type
      is = hecMESH%elem_type_index(itype-1)+1
      iE = hecMESH%elem_type_index(itype  )
      ic_type= hecMESH%elem_type_item(itype)
      if (hecmw_is_etype_link(ic_type)) cycle
      if (hecmw_is_etype_patch(ic_type)) cycle
      nn = hecmw_get_max_node(ic_type)
      if( nn>20 ) stop "Elemental nodes>20!"

      !element loop
      !$omp parallel default(none), &
        !$omp&  private(icel,iiS,j,nodLOCAL,i,ecoord,ddu,du,total_disp, &
        !$omp&  cdsys_ID,coords,thick,qf,isect,ihead,tmp,ndim,ddaux), &
        !$omp&  shared(iS,iE,hecMESH,nn,fstrSOLID,ndof,hecMAT,ic_type,fstrPR, &
        !$omp&         strainEnergy,iter,time,tincr,initt,g_InitialCnd), &
        !$omp&  firstprivate(tt0,ttn,tt)
      !$omp do
      do icel = is, iE

        ! ----- nodal coordinate
        iiS = hecMESH%elem_node_index(icel-1)

        do j = 1, nn
          nodLOCAL(j) = hecMESH%elem_node_item (iiS+j)
          do i = 1, 3
            ecoord(i,j) = hecMESH%node(3*nodLOCAL(j)+i-3)
          enddo
          if( fstrSOLID%TEMP_ngrp_tot > 0 .or. fstrSOLID%TEMP_irres > 0 ) then
            if( isElastoplastic(fstrSOLID%elements(icel)%gausses(1)%pMaterial%mtype) .or. &
                fstrSOLID%elements(icel)%gausses(1)%pMaterial%mtype == NORTON ) then
              tt0(j)=fstrSOLID%last_temp( nodLOCAL(j) )
            else
              tt0(j) = 0.d0
              if( hecMESH%hecmw_flag_initcon == 1 ) tt0(j) = hecMESH%node_init_val_item(nodLOCAL(j))
              if( initt>0 ) tt0(j) = g_InitialCnd(initt)%realval(nodLOCAL(j))
            endif
            ttn(j) = fstrSOLID%last_temp( nodLOCAL(j) )
            tt(j)  = fstrSOLID%temperature( nodLOCAL(j) )
          endif
        enddo

        ! nodal displacement
        do j = 1, nn
          nodLOCAL(j) = hecMESH%elem_node_item (iiS+j)
          do i = 1, ndof
            ddu(i,j) = hecMAT%X(ndof*nodLOCAL(j)+i-ndof)
            du(i,j)  = fstrSOLID%dunode(ndof*nodLOCAL(j)+i-ndof)
            total_disp(i,j) = fstrSOLID%unode(ndof*nodLOCAL(j)+i-ndof)
          enddo
        enddo

        isect = hecMESH%section_ID(icel)
        cdsys_ID = hecMESH%section%sect_orien_ID(isect)
        if( cdsys_ID > 0 ) call get_coordsys(cdsys_ID, hecMESH, fstrSOLID, coords)

        ! ===== calculate the Internal Force
        if( getSpaceDimension( ic_type ) == 2 ) thick = 1.0d0

        if( ic_type == 241 .or. ic_type == 242 .or. ic_type == 231 .or. ic_type == 232 .or. ic_type == 2322 ) then

          if( fstrSOLID%TEMP_ngrp_tot > 0 .or. fstrSOLID%TEMP_irres > 0 ) then
            call UPDATE_C2( ic_type,nn,ecoord(1:3,1:nn),fstrSOLID%elements(icel)%gausses(:), &
              thick,fstrSOLID%elements(icel)%iset,                             &
              total_disp(1:2,1:nn), ddu(1:2,1:nn), qf(1:nn*ndof),              &
              tt(1:nn), tt0(1:nn), ttn(1:nn)  )
          else
            call UPDATE_C2( ic_type,nn,ecoord(1:3,1:nn),fstrSOLID%elements(icel)%gausses(:), &
              thick,fstrSOLID%elements(icel)%iset,           &
              total_disp(1:2,1:nn), ddu(1:2,1:nn), qf(1:nn*ndof) )
          endif

        else if( ic_type == 301 ) then
          isect= hecMESH%section_ID(icel)
          ihead = hecMESH%section%sect_R_index(isect-1)
          thick = hecMESH%section%sect_R_item(ihead+1)
          call UPDATE_C1( ic_type,nn,ecoord(:,1:nn), thick, total_disp(1:3,1:nn), du(1:3,1:nn), &
            qf(1:nn*ndof),fstrSOLID%elements(icel)%gausses(:) )

        else if( ic_type == 361 ) then

          if( fstrSOLID%sections(isect)%elemopt361 == kel361FI ) then ! full integration element
            if( fstrSOLID%TEMP_ngrp_tot > 0 .or. fstrSOLID%TEMP_irres > 0 ) then
              call UPDATE_C3( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_ID, coords, &
                qf(1:nn*ndof), fstrSOLID%elements(icel)%gausses(:), iter, time, tincr, tt(1:nn), tt0(1:nn), ttn(1:nn)  )
            else
              call UPDATE_C3( ic_type,nn,ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_ID, coords, &
                qf(1:nn*ndof), fstrSOLID%elements(icel)%gausses(:), iter, time, tincr )
            endif
          else if( fstrSOLID%sections(isect)%elemopt361 == kel361BBAR ) then ! B-bar element
            if( fstrSOLID%TEMP_ngrp_tot > 0 .or. fstrSOLID%TEMP_irres > 0 ) then
              call UPDATE_C3D8Bbar( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_ID, coords,    &
                qf(1:nn*ndof), fstrSOLID%elements(icel)%gausses(:), iter, time, tincr, tt(1:nn), tt0(1:nn), ttn(1:nn)  )
            else
              call Update_C3D8Bbar( ic_type,nn,ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_ID, coords, &
                qf(1:nn*ndof), fstrSOLID%elements(icel)%gausses(:), iter, time, tincr )
            endif
          else if( fstrSOLID%sections(isect)%elemopt361 == kel361IC ) then ! incompatible element
            if( fstrSOLID%TEMP_ngrp_tot > 0 .or. fstrSOLID%TEMP_irres > 0 ) then
              call UPDATE_C3D8IC( ic_type,nn,ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), ddu(1:3,1:nn), cdsys_ID, coords,&
                qf(1:nn*ndof), fstrSOLID%elements(icel)%gausses(:), iter, time, tincr, &
                fstrSOLID%elements(icel)%aux, ddaux(1:3,1:3), tt(1:nn), tt0(1:nn), ttn(1:nn) )
            else
              call UPDATE_C3D8IC( ic_type,nn,ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), ddu(1:3,1:nn), cdsys_ID, coords,&
                qf(1:nn*ndof), fstrSOLID%elements(icel)%gausses(:), iter, time, tincr, &
                fstrSOLID%elements(icel)%aux, ddaux(1:3,1:3) )
            endif
            fstrSOLID%elements(icel)%aux(1:3,1:3) = fstrSOLID%elements(icel)%aux(1:3,1:3) + ddaux(1:3,1:3)
          else if( fstrSOLID%sections(isect)%elemopt361 == kel361FBAR ) then ! F-bar element
            if( fstrSOLID%TEMP_ngrp_tot > 0 .or. fstrSOLID%TEMP_irres > 0 ) then
              call UPDATE_C3D8Fbar( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_ID, coords,    &
                qf(1:nn*ndof), fstrSOLID%elements(icel)%gausses(:), iter, time, tincr, tt(1:nn), tt0(1:nn), ttn(1:nn)  )
            else
              call Update_C3D8Fbar( ic_type,nn,ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_ID, coords, &
                qf(1:nn*ndof), fstrSOLID%elements(icel)%gausses(:), iter, time, tincr )
            endif
          endif

        else if (ic_type == 341 .or. ic_type == 351 .or. ic_type == 342 .or. ic_type == 352 .or. ic_type == 362 ) then
          if( fstrSOLID%TEMP_ngrp_tot > 0 .or. fstrSOLID%TEMP_irres > 0 ) then
            call UPDATE_C3( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_ID, coords, &
              qf(1:nn*ndof), fstrSOLID%elements(icel)%gausses(:), iter, time, tincr, tt(1:nn), tt0(1:nn), ttn(1:nn)  )
          else
            call UPDATE_C3( ic_type,nn,ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_ID, coords, &
              qf(1:nn*ndof), fstrSOLID%elements(icel)%gausses(:), iter, time, tincr )
          endif

        else if( ic_type == 611) then
          if( fstrPR%nlgeom ) call Update_abort( ic_type, 2 )
          isect = hecMESH%section_ID(icel)
          ihead = hecMESH%section%sect_R_index(isect-1)
          !CALL STF_Beam(ic_type, nn, ecoord, hecMESH%section%sect_R_item(ihead+1:), &
            !     &   material%variables(M_YOUNGS), material%variables(M_POISSON), stiffness(1:nn*ndof,1:nn*ndof))

        else if( ic_type == 641 ) then
          if( fstrPR%nlgeom ) call Update_abort( ic_type, 2 )
          isect = hecMESH%section_ID(icel)
          ihead = hecMESH%section%sect_R_index(isect-1)
          call UpdateST_Beam_641(ic_type, nn, ecoord, total_disp(1:ndof,1:nn), du(1:ndof,1:nn), &
            &    fstrSOLID%elements(icel)%gausses(:), hecMESH%section%sect_R_item(ihead+1:), qf(1:nn*ndof))

        else if( ( ic_type == 741 ) .or. ( ic_type == 743 ) .or. ( ic_type == 731 ) ) then
          if( fstrPR%nlgeom ) call Update_abort( ic_type, 2 )
          isect = hecMESH%section_ID(icel)
          ihead = hecMESH%section%sect_R_index(isect-1)
          thick = hecMESH%section%sect_R_item(ihead+1)
          call UpdateST_Shell_MITC(ic_type, nn, ndof, ecoord(1:3, 1:nn), total_disp(1:ndof,1:nn), du(1:ndof,1:nn), &
            &              fstrSOLID%elements(icel)%gausses(:), qf(1:nn*ndof), thick, 0)

        else if( ic_type == 761 ) then   !for shell-solid mixed analysis
          if( fstrPR%nlgeom ) call Update_abort( ic_type, 2 )
          isect = hecMESH%section_ID(icel)
          ihead = hecMESH%section%sect_R_index(isect-1)
          thick = hecMESH%section%sect_R_item(ihead+1)
          call UpdateST_Shell_MITC33(731, 3, 6, ecoord(1:3, 1:3), total_disp(1:ndof,1:nn), du(1:ndof,1:nn), &
            &              fstrSOLID%elements(icel)%gausses(:), qf(1:nn*ndof), thick, 2)

        else if( ic_type == 781 ) then   !for shell-solid mixed analysis
          if( fstrPR%nlgeom ) call Update_abort( ic_type, 2 )
          isect = hecMESH%section_ID(icel)
          ihead = hecMESH%section%sect_R_index(isect-1)
          thick = hecMESH%section%sect_R_item(ihead+1)
          call UpdateST_Shell_MITC33(741, 4, 6, ecoord(1:3, 1:4), total_disp(1:ndof,1:nn), du(1:ndof,1:nn), &
            &              fstrSOLID%elements(icel)%gausses(:), qf(1:nn*ndof), thick, 1)

        else if ( ic_type == 3414 ) then
          if(fstrSOLID%elements(icel)%gausses(1)%pMaterial%mtype /= INCOMP_NEWTONIAN) then
            write(*, *) '###ERROR### : This element is not supported for this material'
            write(*, *) 'ic_type = ', ic_type, ', mtype = ', fstrSOLID%elements(icel)%gausses(1)%pMaterial%mtype
            call hecmw_abort(hecmw_comm_get_comm())
          endif
          call UPDATE_C3_vp                                                       &
            ( ic_type, nn, ecoord(:,1:nn), total_disp(1:4,1:nn), du(1:4,1:nn), &
            fstrSOLID%elements(icel)%gausses(:) )
          qf = 0.0d0

          !      else if ( ic_type==731) then
          !        call UPDATE_S3(xx,yy,zz,ee,pp,thick,local_stf)
          !        call fstr_local_stf_restore_temp(local_stf, nn*ndof, stiffness)
          !      else if ( ic_type==741) then
          !        call UPDATE_S4(xx,yy,zz,ee,pp,thick,local_stf)
          !        call fstr_local_stf_restore_temp(local_stf, nn*ndof, stiffness)

        else
          write(*, *) '###ERROR### : Element type not supported for nonlinear static analysis'
          write(*, *) ' ic_type = ', ic_type
          call hecmw_abort(hecmw_comm_get_comm())

        endif

        ! ----- calculate the global internal force ( Q(u_{n+1}^{k-1}) )
        do j = 1, nn
          do i = 1, ndof
            !$omp atomic
            fstrSOLID%QFORCE(ndof*(nodLOCAL(j)-1)+i) = fstrSOLID%QFORCE(ndof*(nodLOCAL(j)-1)+i)+qf(ndof*(j-1)+i)
          enddo
        enddo

        ! ----- calculate strain energy
        if(present(strainEnergy))then
          ndim = getSpaceDimension( fstrSOLID%elements(icel)%etype )
          do j = 1, nn
            do i = 1, ndim
              tmp = 0.5d0*( fstrSOLID%elements(icel)%equiForces(ndim*(j-1)+i)+qf(ndim*(j-1)+i) )*ddu(i,j)
              !$omp atomic
              strainEnergy = strainEnergy+tmp
              fstrSOLID%elements(icel)%equiForces(ndim*(j-1)+i) = qf(ndim*(j-1)+i)
            enddo
          enddo
        endif

      enddo ! icel
      !$omp end do
      !$omp end parallel
    enddo   ! itype

    !C
    !C Update for fstrSOLID%QFORCE
    !C
    if( ndof==3 ) then
      call hecmw_update_3_R(hecMESH,fstrSOLID%QFORCE,hecMESH%n_node)
    else if( ndof==2 ) then
      call hecmw_update_2_R(hecMESH,fstrSOLID%QFORCE,hecMESH%n_node)
    else if( ndof==4 ) then
      call hecmw_update_4_R(hecMESH,fstrSOLID%QFORCE,hecMESH%n_node)
    else if( ndof==6 ) then
      call hecmw_update_m_R(hecMESH,fstrSOLID%QFORCE,hecMESH%n_node,6)
    endif

  end subroutine fstr_UpdateNewton

  !> Update elastiplastic status
  subroutine fstr_UpdateState( hecMESH, fstrSOLID, tincr)
    use m_fstr
    use m_static_lib
    use m_ElastoPlastic
    use mCreep
    use mViscoElastic
    type(hecmwST_local_mesh) :: hecMESH     !< hecmw mesh
    type(fstr_solid) :: fstrSOLID   !< fstr_solid
    real(kind=kreal) :: tincr
    integer(kind=kint) :: itype, is, iE, ic_type, icel, ngauss, i

    if( associated( fstrSOLID%temperature ) ) then
      do i = 1, hecMESH%n_node
        fstrSOLID%last_temp(i) = fstrSOLID%temperature(i)
      end do
    endif

    do itype = 1, hecMESH%n_elem_type
      is = hecMESH%elem_type_index(itype-1) + 1
      iE = hecMESH%elem_type_index(itype  )
      ic_type= hecMESH%elem_type_item(itype)
      if( ic_type == 301 ) ic_type = 111
      if( hecmw_is_etype_link(ic_type) ) cycle
      if( hecmw_is_etype_patch(ic_type) ) cycle

      ngauss = NumOfQuadPoints( ic_type )
      do icel = is, iE
        if( isElastoplastic( fstrSOLID%elements(icel)%gausses(1)%pMaterial%mtype ) ) then
          do i = 1, ngauss
            call updateEPState( fstrSOLID%elements(icel)%gausses(i) )
          enddo
        elseif( fstrSOLID%elements(icel)%gausses(1)%pMaterial%mtype == NORTON ) then
          if( tincr>0.d0 ) then
            do i = 1, ngauss
              call updateViscoState( fstrSOLID%elements(icel)%gausses(i) )
            enddo
          endif
        elseif( isViscoelastic( fstrSOLID%elements(icel)%gausses(1)%pMaterial%mtype ) ) then
          if( tincr > 0.d0 ) then
            do i = 1, ngauss
              call updateViscoElasticState( fstrSOLID%elements(icel)%gausses(i) )
            enddo
          endif
        endif

        do i = 1, ngauss
          fstrSOLID%elements(icel)%gausses(i)%strain_bak = fstrSOLID%elements(icel)%gausses(i)%strain
          fstrSOLID%elements(icel)%gausses(i)%stress_bak = fstrSOLID%elements(icel)%gausses(i)%stress
        enddo
      enddo
    enddo
  end subroutine fstr_UpdateState

  subroutine Update_abort( ic_type, flag )
    integer(kind=kint), intent(in) :: ic_type
    integer(kind=kint), intent(in) :: flag

    if( flag == 1 ) then
      write(*,*) '###ERROR### : Element type not supported for static analysis'
    else if( flag == 2 ) then
      write(*,*) '###ERROR### : Element type not supported for nonlinear static analysis'
    endif
    write(*,*) ' ic_type = ', ic_type
    call hecmw_abort(hecmw_comm_get_comm())
  end subroutine

end module m_fstr_Update

max_stars_repo_name: matsuo-san/FrontISTR

max_issues_repo_name: matsuo-san/FrontISTR

max_forks_repo_name: matsuo-san/FrontISTR


text:
      function atdy(z,x,y)
      atdy=y
      return
      end

max_stars_repo_name: sergev/grafor

max_issues_repo_name: sergev/grafor

max_forks_repo_name: sergev/grafor


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

20090816 05:37:29 nbsp Hi and welcome to the Wiki, and welcome back to Davis. Your Aggie Liquor comment about their previous location is a helpful addition to the history of the place. Were they one of the liquor stores that had to be outside of the city limits due to liquor laws (unsure of the year that changed)? There were a couple of liquor stores just east of what is now Cindys, on Chiles Road. Users/robinlaughlin