text:
      PROGRAM TSTRML
C
C Define the error file, the Fortran unit number, the workstation type,
C and the workstation ID to be used in calls to GKS routines.
C
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=1,  IWKID=1)   ! NCGM
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=8,  IWKID=1)   ! X Windows
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=11, IWKID=1)   ! PDF
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=20, IWKID=1)   ! PostScript
C
      PARAMETER (IERRF=6, LUNIT=2, IWTYPE=1,  IWKID=1)
C
C OPEN GKS, OPEN WORKSTATION OF TYPE 1, ACTIVATE WORKSTATION
C
      CALL GOPKS (IERRF, ISZDM)
      CALL GOPWK (IWKID, LUNIT, IWTYPE)
      CALL GACWK (IWKID)
C
C INVOKE DEMO DRIVER
C
      CALL STRML(IERR)
C
C DEACTIVATE AND CLOSE WORKSTATION, CLOSE GKS.
C
      CALL GDAWK (IWKID)
      CALL GCLWK (IWKID)
      CALL GCLKS
C
      STOP
      END
C
      SUBROUTINE STRML (IERROR)
C
C PURPOSE                To provide a simple demonstration of STRMLN.
C
C USAGE                  CALL STRML (IERROR)
C
C ARGUMENTS
C
C ON OUTPUT              IERROR
C                          An integer variable
C                          = 0, if the test was successful,
C                          = 1, the test was not successful.
C
C I/O                    If the test is successful, the message
C
C               STRMLN TEST EXECUTED--SEE PLOT TO CERTIFY
C
C                        is printed on unit 6.  In addition, 1
C                        frame is produced on the machine graphics
C                        device.  In order to determine if the test
C                        was successful, it is necessary to examine
C                        the plot.
C
C PRECISION              Single
C
C LANGUAGE               FORTRAN 77
C
C REQUIRED ROUTINES      STRMLN
C
C REQUIRED GKS LEVEL     0A
C
C ALGORITHM              Routine TSTRML calls routine STRMLN to
C                        produce a plot which depicts the flow and
C                        magnitude of a vector field.
C
      REAL            U(21,25)   ,V(21,25)   ,WRK(1050)
C
C Specify coordinates for plot titles.  The values TX and TY
C define the center of the title string in a 0. to 1. range.
C
      DATA TX/.5/,TY/.9765/
C
C Set the grid dimensions.
C
      DATA NH,NV/21,25/
C
C Initialize the error parameter.
C
      IERROR = 1
C
C Specify horizontal and vertical vector components U and V on
C the rectangular grid.
C
      TPIMX = 2.*3.14/REAL(NH)
      TPJMX = 2.*3.14/REAL(NV)
      DO  20 J=1,NV
         DO  10 I=1,NH
         U(I,J) = SIN(TPIMX*(REAL(I)-1.))
         V(I,J) = SIN(TPJMX*(REAL(J)-1.))
   10    CONTINUE
   20 CONTINUE
C
C Select normalization transformation 0.
C
      CALL GSELNT (0)
C
C Call PLCHLQ to write the plot title.
C
      CALL PLCHLQ (TX,TY,'DEMONSTRATION PLOT FOR ROUTINE STRMLN',16.,
     1            0.,0.)
C
C Define normalization transformation 1, and set up log scaling.
C
      CALL SET(0.1, 0.9, 0.1, 0.9,1.0, 21., 1.0, 25.,1)
C
C Draw the plot perimeter.
C
      CALL PERIM(1,0,1,0)
C
C Call STRMLN for vector field streamlines plot.
C
      CALL STRMLN (U,V,WRK,NH,NH,NV,1,IER)
C
      CALL FRAME
C
      IERROR = 0
      WRITE (6,1001)
      RETURN
C
 1001 FORMAT (' STRMLN TEST EXECUTED--SEE PLOT TO CERTIFY')
C
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


! Copyright (c) 2018-2019, NVIDIA CORPORATION.  All rights reserved.
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.

!ERROR: No explicit type declared for 'f'
function f()
  implicit none
end

!ERROR: No explicit type declared for 'y'
subroutine s(x, y)
  implicit none
  integer :: x
end

subroutine s2
  implicit none
  block
    !ERROR: No explicit type declared for 'i'
    i = 1
  end block
contains
  subroutine s3
    !ERROR: No explicit type declared for 'j'
    j = 2
  end subroutine
end subroutine

module m1
  implicit none
contains
  subroutine s1
    implicit real (a-h)
    a1 = 1.
    h1 = 1.
    !ERROR: No explicit type declared for 'i1'
    i1 = 1
    !ERROR: No explicit type declared for 'z1'
    z1 = 2.
  contains
    subroutine ss1
      implicit integer(f-j) ! overlap with host scope import is OK
      a2 = 1.
      h2 = 1
      i2 = 1
      !ERROR: No explicit type declared for 'z2'
      z2 = 2.
    contains
      subroutine sss1
        implicit none
        !ERROR: No explicit type declared for 'a3'
        a3 = 1.
      end subroutine
    end subroutine
  end subroutine
  subroutine s2
    !ERROR: No explicit type declared for 'b1'
    b1 = 1.
  end subroutine
end module


text:
! module file describing the C++-to-Fortran interface found in f_qd.cpp.

module qdext
  implicit none

  interface
    pure subroutine f_qd_add(a, b, c)
      real*8, intent(in) :: a(4), b(4)
      real*8, intent(out) :: c(4)
    end subroutine

    pure subroutine f_qd_add_qd_d(a, b, c)
      real*8, intent(in) :: a(4), b
      real*8, intent(out) :: c(4)
    end subroutine

    pure subroutine f_qd_add_d_qd(a, b, c)
      real*8, intent(in) :: a, b(4)
      real*8, intent(out) :: c(4)
    end subroutine

    pure subroutine f_qd_sub(a, b, c)
      real*8, intent(in) :: a(4), b(4)
      real*8, intent(out) :: c(4)
    end subroutine

    pure subroutine f_qd_sub_qd_d(a, b, c)
      real*8, intent(in) :: a(4), b
      real*8, intent(out) :: c(4)
    end subroutine

    pure subroutine f_qd_sub_d_qd(a, b, c)
      real*8, intent(in) :: a, b(4)
      real*8, intent(out) :: c(4)
    end subroutine

    pure subroutine f_qd_mul(a, b, c)
      real*8, intent(in) :: a(4), b(4)
      real*8, intent(out) :: c(4)
    end subroutine

    pure subroutine f_qd_mul_qd_d(a, b, c)
      real*8, intent(in) :: a(4), b
      real*8, intent(out) :: c(4)
    end subroutine

    pure subroutine f_qd_mul_d_qd(a, b, c)
      real*8, intent(in) :: a, b(4)
      real*8, intent(out) :: c(4)
    end subroutine

    pure subroutine f_qd_div(a, b, c)
      real*8, intent(in) :: a(4), b(4)
      real*8, intent(out) :: c(4)
    end subroutine

    pure subroutine f_qd_div_qd_d(a, b, c)
      real*8, intent(in) :: a(4), b
      real*8, intent(out) :: c(4)
    end subroutine

    pure subroutine f_qd_div_d_qd(a, b, c)
      real*8, intent(in) :: a, b(4)
      real*8, intent(out) :: c(4)
    end subroutine

    pure subroutine f_qd_sqrt(a, b)
      real*8, intent(in) :: a(4)
      real*8, intent(out) :: b(4)
    end subroutine

    pure subroutine f_qd_sqr(a, b)
      real*8, intent(in) :: a(4)
      real*8, intent(out) :: b(4)
    end subroutine

    pure subroutine f_qd_abs(a, b)
      real*8, intent(in) :: a(4)
      real*8, intent(out) :: b(4)
    end subroutine

    pure subroutine f_qd_npwr(a, n, b)
      real*8, intent(in) :: a(4)
      integer, intent(in) :: n
      real*8, intent(out) :: b(4)
    end subroutine

    pure subroutine f_qd_nroot(a, n, b)
      real*8, intent(in) :: a(4)
      integer, intent(in) :: n
      real*8, intent(out) :: b(4)
    end subroutine

    pure subroutine f_qd_nint(a, b)
      real*8, intent(in) :: a(4)
      real*8, intent(out) :: b(4)
    end subroutine

    pure subroutine f_qd_aint(a, b)
      real*8, intent(in) :: a(4)
      real*8, intent(out) :: b(4)
    end subroutine

    pure subroutine f_qd_floor(a, b)
      real*8, intent(in) :: a(4)
      real*8, intent(out) :: b(4)
    end subroutine

    pure subroutine f_qd_ceil(a, b)
      real*8, intent(in) :: a(4)
      real*8, intent(out) :: b(4)
    end subroutine

    pure subroutine f_qd_log(a, b)
      real*8, intent(in) :: a(4)
      real*8, intent(out) :: b(4)
    end subroutine

    pure subroutine f_qd_log10(a, b)
      real*8, intent(in) :: a(4)
      real*8, intent(out) :: b(4)
    end subroutine

    pure subroutine f_qd_exp(a, b)
      real*8, intent(in) :: a(4)
      real*8, intent(out) :: b(4)
    end subroutine

    pure subroutine f_qd_sin(a, b)
      real*8, intent(in) :: a(4)
      real*8, intent(out) :: b(4)
    end subroutine

    pure subroutine f_qd_cos(a, b)
      real*8, intent(in) :: a(4)
      real*8, intent(out) :: b(4)
    end subroutine

    pure subroutine f_qd_tan(a, b)
      real*8, intent(in) :: a(4)
      real*8, intent(out) :: b(4)
    end subroutine

    pure subroutine f_qd_asin(a, b)
      real*8, intent(in) :: a(4)
      real*8, intent(out) :: b(4)
    end subroutine

    pure subroutine f_qd_acos(a, b)
      real*8, intent(in) :: a(4)
      real*8, intent(out) :: b(4)
    end subroutine

    pure subroutine f_qd_atan(a, b)
      real*8, intent(in) :: a(4)
      real*8, intent(out) :: b(4)
    end subroutine

    pure subroutine f_qd_atan2(a, b, c)
      real*8, intent(in) :: a(4), b(4)
      real*8, intent(out) :: c(4)
    end subroutine

    pure subroutine f_qd_sinh(a, b)
      real*8, intent(in) :: a(4)
      real*8, intent(out) :: b(4)
    end subroutine

    pure subroutine f_qd_cosh(a, b)
      real*8, intent(in) :: a(4)
      real*8, intent(out) :: b(4)
    end subroutine

    pure subroutine f_qd_tanh(a, b)
      real*8, intent(in) :: a(4)
      real*8, intent(out) :: b(4)
    end subroutine

    pure subroutine f_qd_asinh(a, b)
      real*8, intent(in) :: a(4)
      real*8, intent(out) :: b(4)
    end subroutine

    pure subroutine f_qd_acosh(a, b)
      real*8, intent(in) :: a(4)
      real*8, intent(out) :: b(4)
    end subroutine

    pure subroutine f_qd_atanh(a, b)
      real*8, intent(in) :: a(4)
      real*8, intent(out) :: b(4)
    end subroutine

    pure subroutine f_qd_sincos(a, s, c)
      real*8, intent(in) :: a(4)
      real*8, intent(out) :: s(4), c(4)
    end subroutine

    pure subroutine f_qd_sincosh(a, s, c)
      real*8, intent(in) :: a(4)
      real*8, intent(out) :: s(4), c(4)
    end subroutine

    subroutine f_qd_swrite(a, prec, str, maxlen)
      real*8, intent(in) :: a(4)
      integer, intent(in) :: prec, maxlen
      character, intent(out) :: str(*)
    end subroutine

    subroutine f_qd_rand(a)
      real*8, intent(out) :: a(4)
    end subroutine

    pure subroutine f_qd_comp(a, b, r)
      real*8, intent(in) :: a(4), b(4)
      integer, intent(out) :: r
    end subroutine

    pure subroutine f_qd_comp_qd_d(a, b, r)
      real*8, intent(in) :: a(4), b
      integer, intent(out) :: r
    end subroutine

    pure subroutine f_qd_comp_d_qd(a, b, r)
      real*8, intent(in) :: a, b(4)
      integer, intent(out) :: r
    end subroutine

    pure subroutine f_qd_pi(a)
      real*8, intent(out) :: a(4)
    end subroutine

    pure subroutine f_qd_nan(a)
      real*8, intent(out) :: a(4)
    end subroutine

  end interface
end

ext: f

max_stars_repo_name: t-hishinuma/QD

max_issues_repo_name: t-hishinuma/QD

max_forks_repo_name: t-hishinuma/QD


text:
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!   Copyright by The HDF Group.                                               *
!   Copyright by the Board of Trustees of the University of Illinois.         *
!   All rights reserved.                                                      *
!                                                                             *
!   This file is part of HDF5.  The full HDF5 copyright notice, including     *
!   terms governing use, modification, and redistribution, is contained in    *
!   the files COPYING and Copyright.html.  COPYING can be found at the root   *
!   of the source code distribution tree; Copyright.html can be found at the  *
!   root level of an installed copy of the electronic HDF5 document set and   *
!   is linked from the top-level documents page.  It can also be found at     *
!   http://hdfgroup.org/HDF5/doc/Copyright.html.  If you do not have          *
!   access to either file, you may request a copy from help@hdfgroup.org.     *
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!
!
! This file is automatically generated and contains HDF5 Fortran90 type definitions.
!
       MODULE H5FORTRAN_TYPES
         !
         !  HDF5 integers
         !
        INTEGER, PARAMETER :: HADDR_T = 8
        INTEGER, PARAMETER :: HSIZE_T = 8
        INTEGER, PARAMETER :: HSSIZE_T = 8
        INTEGER, PARAMETER :: SIZE_T = 8
        INTEGER, PARAMETER :: Fortran_INTEGER = 4
        INTEGER, PARAMETER :: Fortran_INTEGER_1 = 1
        INTEGER, PARAMETER :: Fortran_INTEGER_2 = 2
        INTEGER, PARAMETER :: Fortran_INTEGER_4 = 4
        INTEGER, PARAMETER :: Fortran_INTEGER_8 = 8
        INTEGER, PARAMETER :: Fortran_REAL_4 = 4
        INTEGER, PARAMETER :: Fortran_REAL_8 = 8
        INTEGER, PARAMETER :: Fortran_REAL_16 = 8
        INTEGER, PARAMETER :: HID_T = 4
        INTEGER, PARAMETER :: Fortran_REAL = 4
        INTEGER, PARAMETER :: Fortran_DOUBLE = 8
        INTEGER, PARAMETER :: H5R_DSET_REG_REF_BUF_SIZE_F = 12

        INTEGER(SIZE_T), PARAMETER :: OBJECT_NAMELEN_DEFAULT_F = -1

        END MODULE H5FORTRAN_TYPES

max_stars_repo_name: CodingNowNow/jiaoyi

max_issues_repo_name: nvsnvyu/hikyuu

max_forks_repo_name: nvsnvyu/hikyuu


text:
        subroutine intrp (itab, xtab, ytab, ntab, x, y, ierr)
c
c Interpolator using CODIM1 algorithm which is admittedly
c obscure  but works well.
c
c       itab - a label between 1 and 20 to identify the table and its
c              most recent search index
c       xtab - array of length ntab containing the x-values
c       ytab -                                     y-values
c       ntab - number of x,y pairs in the table
c          x - independent for which a y-value is desired
c          y - returned interpolated (or extrapolated) value
c       ierr - =0 for ok, -1 for extrapolation
c
        double precision xtab(ntab), ytab(ntab), x, y
        integer itab, ierr, index
        double precision t(4), u(4)
c       integer savind
c       data savind/-1/
c
c----- Only 1 pt in table
        if (ntab .eq. 1) then
                y = ytab(1)
                ierr = 0
                return
        endif
c
c-----
c Locate search index
        call srch (itab, x, xtab, ntab, index, ierr)
c       if (index .eq. savind) go to 2000
c       savind = index
c
c-----
c Set interpolator index flags
        i1 = 2
        i2 = 3
        iload = max0 (index-2, 1)
c
        if (ntab .gt. 2) then
                if (index.eq.   2) i2 = 4
c
                if (index.eq.ntab) i1 = 1
        endif
c
        if (index.gt.2 .and. index.lt.ntab) then
                i1 = 1
                i2 = 4
        endif
c-----
c Load interpolation arrays
        do 1000 i = i1, i2
                j = iload + (i-i1)
                t(i) = xtab (j)
                u(i) = ytab (j)
1000    continue
c
c-----
c Get interpolated value
2000    call codim1 (x, t, u, i1, i2, y)
        return
        end
c
c--------------------------------------------------------------
c
        subroutine srch (itab, x, xtab, ntab, index, ierr)
c
c Search table of x-values to bracket the desired interpolant, x
c
c The returned search index will be:
c       2 - if extrapolation below the table is required
c    ntab -                  above
c   index - points to value just above x in the table if bounded.
c
c The index is saved as a starting point for subsequent entries
c in an array indexed through 'itab' which serves to label the
c set of saved search indices. Itab may be between 1 and 20.
c
c       itab - The table identifier (1-20)
c          x - The value for which an index is desired
c       xtab - The table containing the x-values (array of length ntab)
c       ntab - number of elements in the table
c      index - returned index into the table (points just above x)
c       ierr - 0 for ok, -1 for extrapolation
c
c       Modified to remove entry points.  3/20/86 Valdes
c
        integer ntab, index, init
        double precision xtab(ntab), x
c
c common for subroutines intrp0 and intrpi
c
        common /insvcm/ insave(20)
c
c initialize
        data init/0/
c
c-----
c Initialize
        if (init.eq.0) then
            do 1110 i = 1, 20
1110            insave(i) = 0
            init = 1
        endif
c
c Determine direction of table, ascending or descending
        idir = sign (1.0d0, xtab(ntab) - xtab(1))
c
c-----
c Reset error flag
        ierr = 0
c
c-----
c Check for previous insaved index
        last = insave(itab)
        if (last .eq. 0 .or. last .gt. ntab) then
c
c-----
c no previous entry
                isrch = 1
c check for extrapolation
                if ((x-xtab(   1)) * idir .lt. 0.0d0) go to 2000
                if ((x-xtab(ntab)) * idir .gt. 0.0d0) go to 2100
        else
c
c-----
c previous entry left a valid index
                isrch = last
c
c check for still wihin bounds - difference from above should be opposite
c                                sign of difference from below
c
                if ((xtab(last)-x) * (xtab(last-1)-x) .lt. 0.0d0) then
                        index = last
                        return
                endif
        endif
c
c -----
c Begin searching - first determine direction
c
c This change made because x = xtab(1) was considered extrapolation.
c       if ((x - xtab(isrch)) * idir .gt. 0.0d0) then
        if ((x - xtab(isrch)) * idir .ge. 0.0d0) then
c forward
                do 1100 i = isrch+1, ntab
                        if ((x-xtab(i)) * idir .gt. 0.0d0) go to 1100
                        go to 1500
1100            continue
c fall thru implies extrapolation required at high end
                go to 2100
        else
c
c-----
c negative direction search
                do 1200 i = isrch-1,1,-1
                        if ((x-xtab(i)) * idir .lt. 0.0d0) go to 1200
                        go to 1400
1200            continue
c fall through implies extrapolation at low end
                go to 2000
        endif
c
c-----
c point has been bounded
1400    index = i + 1
        go to 3000
1500    index = i
        go to 3000
c
c-----
c extrapolations
2000    index = 2
        ierr = -1
        go to 3000
2100    index = ntab
        ierr = -1
        go to 3000
c
c-----
c insave index
3000    insave(itab) = index
        end
c
c------
c Subroutine to reset saved index
        subroutine intrp0 (itab)
        integer itab
        common /insvcm/ insave(20)
c
        insave(itab) = 0
        end
c
c-----
c Subroutine to return current index
        subroutine intrpi (itab, ind)
        integer itab, ind
        common /insvcm/ insave(20)
c
        ind = insave(itab)
        end
c
c-------------------------------------------------------------------
c
        subroutine codim1 (x, t, u, i1, i2, y)
c
c this subroutine performs an interposlation in a fashion
c not really understandable, but it works well.
c
c       x - input independent variable
c       t - array of 4 table independents surrounding x if possible
c       u - array of 4 table dependents corresponding to the t array
c
c  i1, i2 - indicators as follows:
c
c           i1 = 1, i2 = 4  :  4 pts available in t and u arrays
c           i1 = 1, i2 = 3  :  3 pts available (x near right edge of table)
c           i1 = 2, i2 = 4  :                  (x near left  edge of table)
c           i1 = 2, i2 = 3  :  2 pts available
c           i1 = 3, i3 = 3  :  1 pt  available
c
c       y - output interpolated (or extrapolated) dependent value
c
        double precision t(4), u(4), x, y
        integer i1, i2
        double precision s, v, z, a1, a2, a3, c1, c2, c3, a4, c4, c5, c6
	double precision e1, e2, p1, p2, slope1, slope2, al, bt, xe

c
c variable xk affects the extrapolation procedure. a value of -1.0
c appears to be a reliable value.
c
        data xk/-1.0d0/
c
        v = x
c the following code is extracted from an original source
c
      a2=v-t(2)
      al=a2/(t(3)-t(2))
      s=al*u(3)+(1.-al)*u(2)
      if(i1.gt.1.and.i2.lt.4)goto1530
      a3=v-t(3)
      if(i1.gt.1)goto1185
1180  a1=v-t(1)
      c1=a2/(t(1)-t(2))*a3/(t(1)-t(3))
      c2=a1/(t(2)-t(1))*a3/(t(2)-t(3))
      c3=a1/(t(3)-t(1))*a2/(t(3)-t(2))
      p1=c1*u(1)+c2*u(2)+c3*u(3)
      if(i2.lt.4)goto1400
1185  a4=v-t(4)
      c4=a3/(t(2)-t(3))*a4/(t(2)-t(4))
      c5=a2/(t(3)-t(2))*a4/(t(3)-t(4))
      c6=a2/(t(4)-t(2))*a3/(t(4)-t(3))
      p2=c4*u(2)+c5*u(3)+c6*u(4)
      if(i1.eq.1)goto1500
1200  if(xk.lt.0.)goto1230
      xe=xk
      goto1260
1230  slope1=abs((u(4)-u(3))/(t(4)-t(3)))
      slope2=abs((u(3)-u(2))/(t(3)-t(2)))
      xe=1.0d0
      if(slope1+slope2.ne.0.)xe=1.-abs(slope1-slope2)/(slope1+slope2)
1260  p1=s+xe*(p2-s)
      goto1500
1400  if(xk.lt.0.)goto1430
      xe=xk
      goto1460
1430  slope1=abs((u(2)-u(1))/(t(2)-t(1)))
      slope2=abs((u(3)-u(2))/(t(3)-t(2)))
      xe=1.0d0
      if(slope1+slope2.ne.0.)xe=1.-abs(slope1-slope2)/(slope1+slope2)
1460  p2=s+xe*(p1-s)
1500  e1=abs(p1-s)
      e2=abs(p2-s)
      if(e1+e2.gt.0.)goto1560
1530  z=s
      goto1700
1560  bt=(e1*al)/(e1*al+(1.-al)*e2)
      z=bt*p2+(1.-bt)*p1
c
1700    y = z
        return
        end
c
c----------------------------------------------------------------------
c
        subroutine lintrp (itab, xtab, ytab, ntab, x, y, ierr)
c
c Linear interpolator with last index save
c
c Arguments are identical to INTRP, and uses the same index search
c scheme so that values for ITAB should not clash with calls
c to INTRP and LINTRP.
c
        double precision xtab(ntab), ytab(ntab), x , y
        integer itab, ierr
c
c----- Only 1 pt in table
        if (ntab .eq. 1) then
                y = ytab (1)
                ierr = 0
                return
        endif
c
c-----locate search index
        call srch (itab, x, xtab, ntab, index, ierr)
c
c----- index points just above x
        y = ytab(index-1) + (x - xtab(index-1)) *
     1   (ytab(index) - ytab(index-1)) / (xtab(index) - xtab(index-1))
c
        return
        end

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


module StrainClass
    use, intrinsic :: iso_fortran_env
    use MathClass
    use ShapeFunctionClass
    implicit none

    type ::  Strain_
        ! Finite strain theory
        real(real64),allocatable ::   F(:,:)
        real(real64),allocatable :: F_n(:,:)
        real(real64),allocatable ::   C(:,:)
        real(real64),allocatable :: C_n(:,:)
        real(real64),allocatable ::   b(:,:)
        real(real64),allocatable ::  Cp(:,:)
        real(real64),allocatable ::Cp_n(:,:)
        real(real64)             :: detF

        ! Hypo-elasto-plasticity
        real(real64),allocatable :: d(:,:)
        real(real64),allocatable ::de(:,:)
        real(real64),allocatable ::dp(:,:)
        real(real64),allocatable :: l(:,:)
        real(real64),allocatable :: w(:,:)

        ! small-strain
        real(real64),allocatable :: eps(:,:)
        real(real64),allocatable :: eps_n(:,:)
        real(real64),allocatable :: eps_p(:,:)
        real(real64),allocatable :: eps_p_n(:,:)
        real(real64),allocatable :: eps_e(:,:)
        real(real64),allocatable :: eps_e_n(:,:)

        integer(int32) :: TheoryID

        character*40 :: StrainTheory

        ! Please input one of following keyphrases

        ! Finite_Elasticity
        ! Finite_ElastoPlasticity
        ! Infinitesimal_Elasticity
        ! Infinitesimal_ElastoPlasticity
        ! Small_strain

    contains
        procedure,public :: init => InitStrain
        procedure,public :: import => importStrain
        procedure,public :: get => getStrain
        procedure,public :: getAll => getAllStrain
        procedure,public :: delete => deleteStrain
    end type

contains

! ###############################
subroutine InitStrain(obj,StrainTheory)
    class(Strain_),intent(inout) :: obj
    character(*),intent(in) :: StrainTheory
    real(real64) :: delta(3,3)

    delta(:,:)=0.0d0
    delta(1,1)=1.0d0
    delta(2,2)=1.0d0
    delta(3,3)=1.0d0

    if(allocated(obj%F) )then
        ! delete old obj
        call obj%delete()
    endif

    obj%StrainTheory=StrainTheory

    if(    trim(obj%StrainTheory)=="Finite_Elasticity")then
        obj%theoryID=1

        ! Finite strain theory
        allocate(obj%    F(3,3) )
        allocate(obj%  F_n(3,3) )
        allocate(obj%    C(3,3) )
        allocate(obj%    b(3,3) )
        allocate(obj%  C_n(3,3) )
        allocate(obj%   Cp(3,3) )
        allocate(obj% Cp_n(3,3) )

        ! Hypo-elasto-plasticity
        allocate(obj%  d(0,0) )
        allocate(obj% de(0,0) )
        allocate(obj% dp(0,0) )
        allocate(obj%  l(0,0) )
        allocate(obj%  w(0,0) )

        ! small-strain
        allocate(obj%  eps(0,0) )
        allocate(obj%  eps_n(0,0) )
        allocate(obj%  eps_e(0,0) )
        allocate(obj%  eps_e_n(0,0) )
        allocate(obj%  eps_p(0,0) )
        allocate(obj%  eps_p_n(0,0) )

        ! Initialize
        obj%    F(:,:) = delta(:,:)
        obj%  F_n(:,:) = delta(:,:)
        obj%    C(:,:) = delta(:,:)
        obj%  C_n(:,:) = delta(:,:)
        obj%   Cp(:,:) = delta(:,:)
        obj% Cp_n(:,:) = delta(:,:)

    elseif(trim(obj%StrainTheory)=="Finite_ElastoPlasticity")then
        obj%theoryID=2

        ! Finite strain theory
        allocate(obj%    F(3,3) )
        allocate(obj%  F_n(3,3) )
        allocate(obj%    C(3,3) )
        allocate(obj%  C_n(3,3) )
        allocate(obj%    b(3,3) )
        allocate(obj%   Cp(3,3) )
        allocate(obj% Cp_n(3,3) )

        ! Hypo-elasto-plasticity
        allocate(obj%  d(0,0) )
        allocate(obj% de(0,0) )
        allocate(obj% dp(0,0) )
        allocate(obj%  l(0,0) )
        allocate(obj%  w(0,0) )

        ! small-strain
        allocate(obj%  eps(0,0) )
        allocate(obj%  eps_n(0,0) )
        allocate(obj%  eps_e(0,0) )
        allocate(obj%  eps_e_n(0,0) )
        allocate(obj%  eps_p(0,0) )
        allocate(obj%  eps_p_n(0,0) )

        obj%    F(:,:) = delta(:,:)
        obj%  F_n(:,:) = delta(:,:)
        obj%    C(:,:) = delta(:,:)
        obj%  C_n(:,:) = delta(:,:)
        obj%   Cp(:,:) = delta(:,:)
        obj% Cp_n(:,:) = delta(:,:)

    elseif(trim(obj%StrainTheory)=="Infinitesimal_Elasticity")then
        obj%theoryID=3

        ! Finite strain theory
        allocate(obj%    F(0,0) )
        allocate(obj%  F_n(0,0) )
        allocate(obj%    C(0,0) )
        allocate(obj%  C_n(0,0) )
        allocate(obj%    b(0,0) )
        allocate(obj%   Cp(0,0) )
        allocate(obj% Cp_n(0,0) )

        ! Hypo-elasto-plasticity
        allocate(obj%  d(3,3) )
        allocate(obj% de(0,0) )
        allocate(obj% dp(0,0) )
        allocate(obj%  l(3,3) )
        allocate(obj%  w(3,3) )

        ! small-strain
        allocate(obj%  eps(3,3) )
        allocate(obj%  eps_n(3,3) )
        allocate(obj%  eps_e(3,3) )
        allocate(obj%  eps_e_n(3,3) )
        allocate(obj%  eps_p(3,3) )
        allocate(obj%  eps_p_n(3,3) )

        ! initialize
        obj%  d(:,:) =delta(:,:)
        obj%  l(:,:) =delta(:,:)
        obj%  w(:,:) =0.0d0
        obj%  eps(:,:) =delta(:,:)
        obj%  eps_n(:,:) =delta(:,:)
        obj%  eps_e(:,:) =delta(:,:)
        obj%  eps_e_n(:,:) =delta(:,:)
        obj%  eps_p(:,:) =0.0d0
        obj%  eps_p_n(:,:) =0.0d0

    elseif(trim(obj%StrainTheory)=="Infinitesimal_ElastoPlasticity")then
        obj%theoryID=4

        ! Finite strain theory
        allocate(obj%    F(0,0) )
        allocate(obj%  F_n(0,0) )
        allocate(obj%    C(0,0) )
        allocate(obj%  C_n(0,0) )
        allocate(obj%    b(0,0) )
        allocate(obj%   Cp(0,0) )
        allocate(obj% Cp_n(0,0) )

        ! Hypo-elasto-plasticity
        allocate(obj%  d(3,3) )
        allocate(obj% de(3,3) )
        allocate(obj% dp(3,3) )
        allocate(obj%  l(3,3) )
        allocate(obj%  w(3,3) )

        ! small-strain
        allocate(obj%  eps(3,3) )
        allocate(obj%  eps_n(3,3) )
        allocate(obj%  eps_e(3,3) )
        allocate(obj%  eps_e_n(3,3) )
        allocate(obj%  eps_p(3,3) )
        allocate(obj%  eps_p_n(3,3) )

        ! initialize
        obj%  d(:,:) =delta(:,:)
        obj% de(:,:) =delta(:,:)
        obj% dp(:,:) =delta(:,:)
        obj%  l(:,:) =delta(:,:)
        obj%  w(:,:) =0.0d0
        obj%  eps(:,:) =delta(:,:)
        obj%  eps_n(:,:) =delta(:,:)
        obj%  eps_e(:,:) =delta(:,:)
        obj%  eps_e_n(:,:) =delta(:,:)
        obj%  eps_p(:,:) =0.0d0
        obj%  eps_p_n(:,:) =0.0d0

    elseif(trim(obj%StrainTheory)=="Small_strain")then
        obj%theoryID=5

        ! Finite strain theory
        allocate(obj%    F(0,0) )
        allocate(obj%  F_n(0,0) )
        allocate(obj%    C(0,0) )
        allocate(obj%  C_n(0,0) )
        allocate(obj%    b(0,0) )
        allocate(obj%   Cp(0,0) )
        allocate(obj% Cp_n(0,0) )

        ! Hypo-elasto-plasticity
        allocate(obj%  d(0,0) )
        allocate(obj% de(0,0) )
        allocate(obj% dp(0,0) )
        allocate(obj%  l(0,0) )
        allocate(obj%  w(0,0) )

        ! small-strain
        allocate(obj%  eps(3,3) )
        allocate(obj%  eps_n(3,3) )
        allocate(obj%  eps_e(3,3) )
        allocate(obj%  eps_e_n(3,3) )
        allocate(obj%  eps_p(3,3) )
        allocate(obj%  eps_p_n(3,3) )

        obj%  eps(:,:) =delta(:,:)
        obj%  eps_n(:,:) =delta(:,:)
        obj%  eps_e(:,:) =delta(:,:)
        obj%  eps_e_n(:,:) =delta(:,:)
        obj%  eps_p(:,:) =0.0d0
        obj%  eps_p_n(:,:) =0.0d0

    else
        print *, trim(StrainTheory)
        print *, "Please input one of following keyphrases"
        print *, " ---->"
        print *, "Finite_Elasticity"
        print *, "Finite_ElastoPlasticity"
        print *, "Infinitesimal_Elasticity"
        print *, "Infinitesimal_ElastoPlasticity"
        print *, "Small_strain"
        print *, " <----"
    endif

end subroutine
! ###############################

! ###############################
subroutine deleteStrain(obj)
    class(Strain_),intent(inout) :: obj

        ! Finite strain theory
    deallocate(obj%   F )
    deallocate(obj% F_n )
    deallocate(obj%   C )
    deallocate(obj% C_n )
    deallocate(obj%   b )
    deallocate(obj%  Cp )
    deallocate(obj%Cp_n )
    obj%detF = 0.0d0

    ! Hypo-elasto-plasticity
    deallocate(obj% d )
    deallocate(obj%de )
    deallocate(obj%dp )
    deallocate(obj% l )
    deallocate(obj% w )

    ! small-strain
    deallocate(obj% eps )
    deallocate(obj% eps_n )
    deallocate(obj%  eps_e )
    deallocate(obj%  eps_e_n )
    deallocate(obj%  eps_p )
    deallocate(obj%  eps_p_n )

    obj%TheoryID = 0

    obj%StrainTheory = " "
end subroutine
! ###############################

! ###############################
subroutine importStrain(obj,F,F_n,C,C_n,b,Cp,Cp_n,d,de,dp,l,w,eps,eps_n)
    class(Strain_),intent(inout) :: obj
    real(real64),optional,intent(in) :: F(:,:),F_n(:,:),C(:,:),C_n(:,:),b(:,:)
    real(real64),optional,intent(in) :: Cp(:,:),Cp_n(:,:),d(:,:),de(:,:)
    real(real64),optional,intent(in) :: dp(:,:),l(:,:),w(:,:),eps(:,:),eps_n(:,:)

    if(present(F))then
        if(size(obj%F,1)/= size(F,1) .or. size(obj%F,2) /= size(F,2) )then
            print *, "ERROR :: importStrain :: invalid size "
            return
        else
            obj%F(:,:) = F(:,:)
        endif
    endif
    if(present(F_n))then
        if(size(obj%F_n,1)/= size(F_n,1) .or. size(obj%F_n,2) /= size(F_n,2) )then
            print *, "ERROR :: importStrain :: invalid size "
            return
        else
            obj%F_n(:,:) = F_n(:,:)
        endif
    endif
    if(present(C))then
        if(size(obj%C,1)/= size(C,1) .or. size(obj%C,2) /= size(C,2) )then
            print *, "ERROR :: importStrain :: invalid size "
            return
        else
            obj%C(:,:) = C(:,:)
        endif
    endif
    if(present(C_n))then
        if(size(obj%C_n,1)/= size(C_n,1) .or. size(obj%C_n,2) /= size(C_n,2) )then
            print *, "ERROR :: importStrain :: invalid size "
            return
        else
            obj%C_n(:,:) = C_n(:,:)
        endif
    endif
    if(present(b))then
        if(size(obj%b,1)/= size(b,1) .or. size(obj%b,2) /= size(b,2) )then
            print *, "ERROR :: importStrain :: invalid size "
            return
        else
            obj%b(:,:) = b(:,:)
        endif
    endif
    if(present(Cp))then
        if(size(obj%Cp,1)/= size(Cp,1) .or. size(obj%Cp,2) /= size(Cp,2) )then
            print *, "ERROR :: importStrain :: invalid size "
            return
        else
            obj%Cp(:,:) = Cp(:,:)
        endif
    endif
    if(present(Cp_n))then
        if(size(obj%Cp_n,1)/= size(Cp_n,1) .or. size(obj%Cp_n,2) /= size(Cp_n,2) )then
            print *, "ERROR :: importStrain :: invalid size "
            return
        else
            obj%Cp_n(:,:) = Cp_n(:,:)
        endif
    endif
    if(present(d))then
        if(size(obj%d,1)/= size(d,1) .or. size(obj%d,2) /= size(d,2) )then
            print *, "ERROR :: importStrain :: invalid size "
            return
        else
            obj%d(:,:) = d(:,:)
        endif
    endif
    if(present(de))then
        if(size(obj%de,1)/= size(de,1) .or. size(obj%de,2) /= size(de,2) )then
            print *, "ERROR :: importStrain :: invalid size "
            return
        else
            obj%de(:,:) = de(:,:)
        endif
    endif
    if(present(dp))then
        if(size(obj%dp,1)/= size(dp,1) .or. size(obj%dp,2) /= size(dp,2) )then
            print *, "ERROR :: importStrain :: invalid size "
            return
        else
            obj%dp(:,:) = dp(:,:)
        endif
    endif
    if(present(l))then
        if(size(obj%l,1)/= size(l,1) .or. size(obj%l,2) /= size(l,2) )then
            print *, "ERROR :: importStrain :: invalid size "
            return
        else
            obj%l(:,:) = l(:,:)
        endif
    endif
    if(present(w))then
        if(size(obj%w,1)/= size(w,1) .or. size(obj%w,2) /= size(w,2) )then
            print *, "ERROR :: importStrain :: invalid size "
            return
        else
            obj%w(:,:) = w(:,:)
        endif
    endif
    if(present(eps))then
        if(size(obj%eps,1)/= size(eps,1) .or. size(obj%eps,2) /= size(eps,2) )then
            print *, "ERROR :: importStrain :: invalid size "
            return
        else
            obj%eps(:,:) = eps(:,:)
        endif
    endif
    if(present(eps_n))then
        if(size(obj%eps_n,1)/= size(eps_n,1) .or. size(obj%eps_n,2) /= size(eps_n,2) )then
            print *, "ERROR :: importStrain :: invalid size "
            return
        else
            obj%eps_n(:,:) = eps_n(:,:)
        endif
    endif

end subroutine
! ###############################

! ###############################
subroutine getallStrain(obj,ShapeFunction)
    class(Strain_),intent(inout) :: obj
    class(ShapeFunction_),intent(in)::ShapeFunction
    real(real64),allocatable :: fmat(:,:) ,Jmat_n(:,:),ddudgzi(:,:),ddudx_n(:,:),ddudx(:,:),&
        dudx(:,:),JmatInv_n(:,:),dudgzi(:,:)

    ! for Finite strain theory
    if(size(obj%F,1) ==3 )then
        allocate(fmat(3,3) ,Jmat_n(3,3),ddudgzi(3,3),ddudx_n(3,3),JmatInv_n(3,3))
        fmat(:,:) = 0.0d0
        fmat(1,1) = 1.0d0
        fmat(2,2) = 1.0d0
        fmat(3,3) = 1.0d0
        Jmat_n(:,:) = matmul(ShapeFunction%dNdgzi,ShapeFunction%ElemCoord_n)
        ddudgzi(:,:) = matmul(ShapeFunction%dNdgzi,ShapeFunction%du)
        call inverse_rank_2(Jmat_n,JmatInv_n)
        ddudx_n(:,:) = matmul(ddudgzi, JmatInv_n  )
        fmat(:,:) =fmat(:,:) + ddudx_n(:,:)
        obj%F = matmul(fmat,obj%F_n)
        call obj%get(C=.true.,b=.true.,detF=.true.)
    endif

    ! for infinitesimal strain theory
    if(size(obj%d,1) == 3 )then
        allocate(ddudx(3,3),ddudgzi(3,3) )
        ddudgzi(:,:) = matmul(ShapeFunction%dNdgzi,ShapeFunction%du)
        ddudx(:,:) = matmul(ddudgzi, ShapeFunction%JmatInv  )
        obj%l(:,:) = ddudx(:,:)
        ! from velocity gradient tensor l, spin tensor w and stretch tensor d is computed.
        call obj%get(d=.true.)
        call obj%get(w=.true.)
        if(size(obj%de,1)==3 .and. size(obj%dp,1)==3 )then
            call obj%get(de=.true.)
        endif
    endif

    ! for small strain theory
    if(size(obj%eps,1 )==3 )then
        if( size(obj%d,1) ==3)then
            ! forward Euler
            obj%eps(:,:) =obj%eps_n(:,:)+obj%d(:,:)
            ! other integral scheme will be implemented.
        else
            ! small strain
            ! Caution :: du here is seen as u
            allocate(dudgzi(3,3),dudx(3,3) )
            dudgzi(:,:) = matmul(ShapeFunction%dNdgzi,ShapeFunction%du)
            dudx(:,:) = matmul(dudgzi, ShapeFunction%JmatInv  )
            ! Here is the small strain tensor.
            obj%eps(:,:) = 0.50d0*(dudx(:,:) + transpose(dudx))
        endif
    endif

end subroutine
! ###############################

! ###############################
subroutine getStrain(obj,C,b,d,w,de,detF)
    class(Strain_),intent(inout)::obj
    logical,optional,intent(in) :: C,b,d,w,detF,de

    if(present(C) )then
        if(C .eqv. .true.)then
            obj%C(:,:) = matmul(transpose(obj%F),obj%F)
        endif
    endif

    if(present(b) )then
        if(b .eqv. .true.)then
            obj%b(:,:) = matmul(obj%F,transpose(obj%F))
        endif
    endif

    if(present(d) )then
        if(d .eqv. .true.)then
            obj%d(:,:) = 0.50d0*(obj%l(:,:) + transpose(obj%l) )
        endif
    endif

    if(present(de) )then
        if(de .eqv. .true.)then
            obj%de(:,:) = obj%d(:,:) - obj%dp(:,:)
        endif
    endif

    if(present(w) )then
        if(w .eqv. .true.)then
            obj%w(:,:) = 0.50d0*(obj%l(:,:) - transpose(obj%l) )
        endif
    endif

    if(present(detF) )then
        if(detF .eqv. .true.)then
            obj%detF = det_mat(obj%F,size(obj%F,1) )
        endif
    endif

end subroutine
! ###############################

end module


text:
module init_phi_module

  use multifab_module

  implicit none

  private

  public :: init_phi

contains

  subroutine init_phi(phi,dx,prob_lo)

    type(multifab) , intent(inout) :: phi
    real(kind=dp_t), intent(in   ) :: dx
    real(kind=dp_t), intent(in   ) :: prob_lo(phi%dim)

    ! local variables
    integer :: lo(phi%dim), hi(phi%dim)
    integer :: dm, ng, i

    real(kind=dp_t), pointer :: dp(:,:,:,:)

    ! set these here so we don't have to pass them into the subroutine
    dm = phi%dim
    ng = phi%ng

    do i=1,nfabs(phi)
       dp => dataptr(phi,i)
       lo = lwb(get_box(phi,i))
       hi = upb(get_box(phi,i))
       select case(dm)
       case (2)
          call init_phi_2d(dp(:,:,1,1), ng, lo, hi, prob_lo, dx)
       case (3)
          call init_phi_3d(dp(:,:,:,1), ng, lo, hi, prob_lo, dx)
       end select
    end do

    ! fill ghost cells
    ! this only fills periodic ghost cells and ghost cells for neighboring
    ! grids at the same level.  Physical boundary ghost cells are filled
    ! using multifab_physbc.  But this problem is periodic, so this
    ! call is sufficient.
    call multifab_fill_boundary(phi)

  end subroutine init_phi

  subroutine init_phi_2d(phi, ng, lo, hi, prob_lo, dx)

    integer          :: lo(2), hi(2), ng
    double precision :: phi(lo(1)-ng:,lo(2)-ng:)
    double precision :: prob_lo(2)
    double precision :: dx

    ! local varables
    integer          :: i,j
    double precision :: x,y,r2

    do j=lo(2),hi(2)
         y = prob_lo(2) + (dble(j)+0.5d0) * dx
         do i=lo(1),hi(1)
            x = prob_lo(1) + (dble(i)+0.5d0) * dx

            r2 = ((x-0.25d0)**2 + (y-0.25d0)**2) / 0.01d0
            phi(i,j) = 1.d0 + exp(-r2)

         end do
      end do

    end subroutine init_phi_2d

    subroutine init_phi_3d(phi, ng, lo, hi, prob_lo, dx)

    integer          :: lo(3), hi(3), ng
    double precision :: phi(lo(1)-ng:,lo(2)-ng:,lo(3)-ng:)
    double precision :: prob_lo(3)
    double precision :: dx

    ! local varables
    integer          :: i,j,k
    double precision :: x,y,z,r2

    !$omp parallel do private(i,j,k,x,y,z,r2)
    do k=lo(3),hi(3)
       z = prob_lo(3) + (dble(k)+0.5d0) * dx
       do j=lo(2),hi(2)
          y = prob_lo(2) + (dble(j)+0.5d0) * dx
          do i=lo(1),hi(1)
             x = prob_lo(1) + (dble(i)+0.5d0) * dx

             r2 = ((x-0.25d0)**2 + (y-0.25d0)**2 + (z-0.25d0)**2) / 0.01d0
             phi(i,j,k) = 1.d0 + exp(-r2)

          end do
       end do
    end do
    !$omp end parallel do

  end subroutine init_phi_3d

end module init_phi_module

max_stars_repo_name: memmett/BoxLib

max_issues_repo_name: memmett/BoxLib

max_forks_repo_name: memmett/BoxLib


program main
    implicit none
    character(len=100) :: string
    character(len=100) :: deleteSpace

    write (*,'(A)') 'Please input a string contains some spaces:'
    read  (*,'(A)') string
    write (*,'(A,/,A)') 'The new string after deleting spaces is:', trim(deleteSpace(string))
end program main

function deleteSpace(string)
    implicit none
    character(len=*) :: string
    character(len=*) :: deleteSpace
    integer :: counter, counter2 = 0

    deleteSpace = ''
    do counter = 0, len(string)
        if (string(counter:counter) /= ' ') then
            deleteSpace(counter2:counter2) = string(counter:counter)
            counter2 = counter2 + 1
        end if
    end do
    return
end function


! { dg-do compile }
module op

   implicit none

   type a
      integer i
   end type a

   type b
      real i
   end type b

   interface operator(==)
      module procedure f1
   end interface operator(.eq.)
   interface operator(.eq.)
      module procedure f2
   end interface operator(==)

   interface operator(/=)
      module procedure f1
   end interface operator(.ne.)
   interface operator(.ne.)
      module procedure f2
   end interface operator(/=)

   interface operator(<=)
      module procedure f1
   end interface operator(.le.)
   interface operator(.le.)
      module procedure f2
   end interface operator(<=)

   interface operator(<)
      module procedure f1
   end interface operator(.lt.)
   interface operator(.lt.)
      module procedure f2
   end interface operator(<)

   interface operator(>=)
      module procedure f1
   end interface operator(.ge.)
   interface operator(.ge.)
      module procedure f2
   end interface operator(>=)

   interface operator(>)
      module procedure f1
   end interface operator(.gt.)
   interface operator(.gt.)
      module procedure f2
   end interface operator(>)

   contains

      function f2(x,y)
         logical f2
         type(a), intent(in) :: x, y
      end function f2

      function f1(x,y)
         logical f1
         type(b), intent(in) :: x, y
      end function f1

end module op


        !COMPILER-GENERATED INTERFACE MODULE: Fri Mar  6 15:34:27 2020
        ! This source file is for reference only and may not completely
        ! represent the generated interface used by the compiler.
        MODULE GWTCNR__genmod
          INTERFACE
            SUBROUTINE GWTCNR(ALATO,JG)
              REAL(KIND=8) :: ALATO(*)
              INTEGER(KIND=4) :: JG
            END SUBROUTINE GWTCNR
          END INTERFACE
        END MODULE GWTCNR__genmod


!==============================================================================!
  subroutine Grid_Mod_Find_Nodes_Cells(grid)
!------------------------------------------------------------------------------!
!   Cells around each node are needed for Lagrangian particle tracking.        !
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  type(Grid_Type) :: grid
!-----------------------------------[Locals]-----------------------------------!
  integer              :: c, ln, n, run, s1, s2, n1, n2, ln1, ln2, nc1, nc2, nu
  integer              :: max_n_cells
  integer, allocatable :: cell_list(:)
  real                 :: x1, x2, y1, y2, z1, z2

! Just for checking, erase this later
! integer :: lc
! real    :: xn, xc, yn, yc, zn, zc, max_del, min_del, del
!==============================================================================!

  ! Allocate memory for node coordinates
  n = grid % n_nodes
  allocate(grid % nodes_n_cells(1:n))
  grid % nodes_n_cells(:) = 0

  !--------------------------------------------------------!
  !   Find maximum number of cells surrounding each node   !
  !   and allocate memory in run 1, store them in run 2    !
  !--------------------------------------------------------!

  do run = 1, 2
    grid % nodes_n_cells(:) = 0  ! re-initialize the cell count

    do c = -grid % n_bnd_cells, grid % n_cells
      do ln = 1, grid % cells_n_nodes(c)  ! local node number
        n = grid % cells_n(ln, c)         ! global node number

        ! Increase number of cells surrounding the this node by one
        grid % nodes_n_cells(n) = grid % nodes_n_cells(n) + 1

        ! Store the current cell in the second run
        if(run .eq. 2) then
          grid % nodes_c(grid % nodes_n_cells(n), n) = c
        end if
      end do
    end do

    ! Allocate memory for cells surrounding each node
    if(run .eq. 1) then
      max_n_cells = maxval(grid % nodes_n_cells)
      allocate(grid % nodes_c(1:max_n_cells, 1:grid % n_nodes))
      allocate(cell_list(max_n_cells*2));  cell_list(:) = 0
    end if

  end do

  !-------------------------------------------------------------------!
  !   Connect nodes' cell lists through periodic boundaries as well   !
  !-------------------------------------------------------------------!

  ! Run through three possible periodic directions.
  ! (If you do it less, some cells might get missed
  !  during the browsing and gathering procedure)
  do run = 1, 3
    do s1 = 1, grid % n_faces
      s2 = grid % faces_s(s1)

      ! Take only faces with shadows - ony those have shadow nodes
      if(s2 > 0) then
        do ln1 = 1, grid % faces_n_nodes(s1)
          n1 = grid % faces_n(ln1, s1)
          x1 = grid % xn(n1)
          y1 = grid % yn(n1)
          z1 = grid % zn(n1)
          do ln2 = 1, grid % faces_n_nodes(s2)
            n2 = grid % faces_n(ln2, s2)
            x2 = grid % xn(n2) + grid % dx(s2)
            y2 = grid % yn(n2) + grid % dy(s2)
            z2 = grid % zn(n2) + grid % dz(s2)

            ! Nodes are matching, add cells to each other's list
            if(Math_Mod_Distance_Squared(x1, y1, z1, x2, y2, z2) < PICO) then

              ! Gather list of cells for both nodes and make a unique sore
              nc1 = grid % nodes_n_cells(n1)  ! number of cells arund node 1
              nc2 = grid % nodes_n_cells(n2)  ! number of cells arund node 2
              cell_list(    1:nc1    ) = grid % nodes_c(1:nc1, n1)
              cell_list(nc1+1:nc1+nc2) = grid % nodes_c(1:nc2, n2)
              call Sort_Mod_Unique_Int(cell_list(1:nc1+nc2), nu)

              ! Copy unique list of cells to both nodes' lists
              grid % nodes_n_cells(n1) = nu
              grid % nodes_n_cells(n2) = nu
              grid % nodes_c(1:nu, n1) = cell_list(1:nu)
              grid % nodes_c(1:nu, n2) = cell_list(1:nu)
            end if
          end do
        end do
      end if
    end do
  end do

! ! Check #1, save those nodes' cells
! do c = 1, grid % n_cells - grid % comm % n_buff_cells
!   do ln = 1, grid % cells_n_nodes(c)  ! local node number
!     n = grid % cells_n(ln, c)         ! global node number
!     write(600+this_proc, '(a,i2,a,36i6)')  &
!           'nc=', grid % nodes_n_cells(n),  &
!           ' cell list=', grid % nodes_c(1:max_n_cells, n)
!     write(600+this_proc, '(a,36i6)')  &
!           '           proc=',  &
!           grid % comm % cell_proc(grid % nodes_c(1:max_n_cells, n))
!   end do
! end do

! ! Check #2, erase this later
! max_del = -HUGE
! min_del = +HUGE
!
! do n = 1, grid % n_nodes
!   do lc = 1, grid % nodes_n_cells(n)  ! local cell number
!     c = grid % nodes_c(lc, n)         ! global cell number
!
!     xn = grid % xn(n)
!     yn = grid % yn(n)
!     zn = grid % zn(n)
!
!     xc = grid % xc(c)
!     yc = grid % yc(c)
!     zc = grid % zc(c)
!
!     del = sqrt( (xn-xc)**2 + (yn-yc)**2 + (zn-zc)**2 )
!
!     min_del = min(del, min_del)
!     max_del = max(del, max_del)
!   end do
! end do
!
! print *, '# Checking: min and max del: ', min_del, max_del
! stop

  end subroutine


text:
      SUBROUTINE CINVPR (EED,METHOD,NFOUND)
C
C     GIVEN REAL OR COMPLEX MATRICIES, CINVPR WILL SOLVE FOR ALL OF
C     THE EIGENVALUES AND EIGENVECTORS WITHIN A SPECIFIED REGION
C
C     DEFINITION OF INPUT AND OUTPUT PARAMETERS
C
C     FILEK(7) = MATRIX CONTROL BLOCK FOR THE INPUT STIFFNESS MATRIX K
C     FILEM(7) = MATRIX CONTROL BLOCK FOR THE INPUT MASS MATRIX M
C     FILEB(7) = MATRIX CONTROL BLOCK FOR THE INPUT DAMPING MATRIX B
C     FILELM(7)= MATRIX CONTROL BLOCK FOR THE OUTPUT EIGENVALUES
C     FILEVC(7)= MATRIX CONTROL BLOCK FOR THE OUTPUT EIGENVECTORS
C     DMPFIL   = FILE CONTAINING THE EIGENVALUE SUMMARY
C     SR1FIL-  = SCRATCH FILES USED INTERNALLY
C     SR0FIL
C     EPS      = CONVERGENCE CRITERIA
C     NOREG    = NUMBER OF REGIONS INPUT
C     REG(1,I) = X1 FOR REGION I
C     REG(2,I) = Y1 FOR REGION I
C     REG(3,I) = X2 FOR REGION I
C     REG(4,I) = Y2 FOR REGION I
C     REG(5,I) = L1 FOR REGION I
C     REG(6,I) = NO. OF DESIRED  ROOTS FOR REGION I
C     REG(7,I) = NO. OF ESTIMATED ROOTS IN REGION I
C
C
      LOGICAL          NOLEFT
      INTEGER          METHOD    ,EED      ,EIGC(2)  ,PHIDLI  ,
     1                 SWITCH    ,SCRFIL   ,IHEAD(10),IREG(7,1)
      INTEGER          NAME(2)   ,FILELM   ,FILEVC   ,
     1                 REAL      ,RDP      ,TYPEK    ,TYPEM    ,
     2                 TYPEB     ,COMFLG   ,IZ(1)    ,DMPFIL   ,
     3                 TIMED     ,FILEK    ,T1       ,T2       ,
     4                 FILEB     ,FILEM
      REAL             L         ,L1       ,MAXMOD
      DOUBLE PRECISION LAM1      ,DZ(1)    ,MINDIA
      DOUBLE PRECISION LAMBDA    ,LMBDA    ,DTEMP(2)
      CHARACTER        UFM*23
      COMMON /XMSSG /  UFM
      COMMON /CDCMPX/  IDUM(30)  ,MINDIA
      COMMON /CINVPX/  FILEK(7)  ,FILEM(7) ,FILEB(7) ,FILELM(7),
     1                 FILEVC(7) ,DMPFIL   ,SCRFIL(11),NOREG   ,
     2                 EPS       ,REG(7,10),PHIDLI
      COMMON /NAMES /  RD        ,RDREW    ,WRT       ,WRTREW  ,
     1                 REW       ,NOREW    ,EOFNRW    ,RSP     ,
     1                 RDP
      COMMON /SYSTEM/  KSYSTM(65)
      COMMON /ZZZZZZ/  Z(1)
      COMMON /OUTPUT/  HEAD(1)
      COMMON /CINVXX/  LAMBDA(2) ,SWITCH   ,COMFLG    ,LMBDA(2),
     1                 ITER      ,TIMED    ,NOCHNG    ,RZERO   ,
     2                 IND       ,IVECT    ,KREG      ,REAL    ,
     3                 LEFT      ,NORTHO   ,NOROOT    ,NZERO   ,
     4                 LAM1(2)   ,MAXMOD   ,NODES     ,NOEST   ,
     5                 ISTART    ,IND1     ,ITERX     ,ISYM
      EQUIVALENCE      (KSYSTM(1),ISYS )   ,(IREG(1,1),REG(1,1))
      EQUIVALENCE      (FILEK(5) ,TYPEK)   ,(FILEM(5),TYPEM)   ,
     1                 (FILEB(5) ,TYPEB)   ,(IZ(1),Z(1))
      EQUIVALENCE      (ANODES   ,NODES)   ,(ANOEST,NOEST)     ,
     1                 (Z(1)     ,DZ(1))   ,(KSYSTM(2),NOUT )
      DATA     IHEAD/  0,1009,2,7*0 /
      DATA     EIGC /  207,2        /
      DATA     NAME /  4HCINV,4HPR  /
      DATA     SIGN /  1.0          /
C
C     DEFINITION OF INTERNAL PARAMETERS
C
C     REAL     = 0 - ALL MATRICIES ARE REAL
C                1 - AT LEAST ONE MATRIX IS COMPLEX
C     NSHIFT   = NO. OF SHIFT POINTS IN A REGION
C     NODES    = NO. OF DESIRED ROOTS IN A REGION
C     NOEST    = NO. OF ESTIMATED ROOTS IN A REGION
C     SHIFT    = INDEX OF THE CURRENT SHIFT POINT
C     ISHIFT   = INDEX OF THE CURRENT SHIFT POINT
C     IMIN     = LOWEST INDEX OF THE COMPLETED SHIFT POINTS
C     IMAX     = HIGHEST INDEX OF COMPLETED SHIFT POINTS
C
C     FILE ALLOCATION
C
C     SR1FIL CONTAINS (LAMBDA**2*M + LAMBDA*B + K)
C     SR2FIL CONTAINS -(B+LAMBDA*M)
C     SR3FIL CONTAINS THE LOWER TRIANGLE OF THE DECOMPOSED DYNAMIC MTRX
C     SR4FIL CONTAINS THE UPPER TRIANGLE OF THE DECOMPOSED DYNAMIC MTRX
C     SR5FIL IS USED AS A SCRATCH FOR CDCOMP
C     SR6FIL IS USED AS A SCRATCH FOR CDCOMP
C     SR7FIL IS USED AS A SCRATCH FOR CDCOMP
C     SR8FIL CONTAINS THE LOWER TRIANGLE L
C     SR9FIL CONTAINS THE UPPER TRIANGLE U
C     SR0FIL CONTAINS THE LEFT EIGENVECTORS
C     S11FIL CONTAINS  -(B + LAMBDA*M)
C
C     DEFINITION OF INTERNAL PARAMETERS
C
C     IND      = AN INDEX FOR GENERATING VARIOUS STARTING VECTORS
C     ITER     = TOTAL NUMBER OF ITERATIONS
C     NODCMP   = TOTAL NUMBER OF DECOMPOSITIONS
C     NOSTRT   = NUMBER OF STARTING POINTS USED
C     NOMOVS   = NUMBER OF TIMES A STARTING POINT HAD TO BE MOVED
C     RZERO    = DISTANCE FROM THE STARTING POINT TO THE CORNER OF THE
C                PARALELOGRAM
C     NOCHNG   = COUNT OF THE NUMBER OF MOVES WHILE LOOKING FOR ONE ROO
C     COMFLG   = 0 -
C              = 1 -
C              = 2 -
C              = 3 -
C              = 4 -
C              = 5 -
C              = 6 -
C     SWITCH   =
C     IVECT    =
C     KREG     = 0-NO VECTORS FOUND IN SEARCH AREA YET
C                1- A VECTOR HAS BEEN FOUND IN THE SEARCH AREA
C     ISING    = SINGULARITY FLAG
C     ITERM    = REASON FOR TERMINATING
C              = 1 - 2SINGULARITIES IN A ROW
C              = 2 - 4 MOVES WHILE TRACKING ONE ROOT
C              = 3 - ALL REGIONS COMPLETED
C              = 4 - 3*NOEST FOUND
C              = 5 - ALL ROOTS FOUND
C              = 8 - 200 ITERATIONS WITH ONE MOVE WITHOUR CONVERGING
C     TIMED    = TIME TOO FORM AND DECOMPOSE THE DYNAMIC MATRIX
C     LEFT     = 1 - DECOMPOSE MATRIX FOR THE COMPUTATION OF THE LEFT
C                EIGENVECTORS
C
      CALL SSWTCH (12,IDIAG)
      IND1 = 0
      NZ   = KORSZ(Z)
      CALL KLOCK (ISTART)
      IBUF = NZ - ISYS - 2
      IFILE= FILELM(1)
      CALL OPEN  (*500,FILELM,Z(IBUF),WRTREW)
      CALL CLOSE (FILELM,REW)
      IFILE = FILEVC (1)
      CALL OPEN  (*500,FILEVC,Z(IBUF),WRTREW)
      CALL CLOSE (FILEVC,REW)
      CALL GOPEN (DMPFIL,Z(IBUF),WRTREW)
      CALL CLOSE (DMPFIL,EOFNRW)
      IFILE = SCRFIL(10)
      CALL OPEN (*500,IFILE,Z(IBUF),WRTREW)
      CALL CLOSE (IFILE,REW)
      NOLEFT = .FALSE.
      IZ(1)  = 204
      CALL RDTRL (IZ)
      IF (IZ(1) .LT. 0) NOLEFT = .TRUE.
      NORTHO = 0
      NROW   = 2*FILEK(3)
      NROW2  = 2*NROW
      ISYM   = 1
      IF (FILEK(1).NE.0 .AND. FILEK(4).NE.6) GO TO 2
      IF (FILEM(1).NE.0 .AND. FILEM(4).NE.6) GO TO 2
      IF (FILEB(1).NE.0 .AND. FILEB(4).NE.6) GO TO 2
      ISYM = 0
    2 CONTINUE
C
C     PICK UP REGION PARAMETERS
C
      CALL PRELOC (*500,Z(IBUF),EED)
      CALL LOCATE (*500,Z(IBUF),EIGC(1),FLAG)
    6 CALL FREAD (EED,IREG,10,0)
      IF (METHOD.EQ.IREG(1,1) .OR. METHOD.EQ.-1) GO TO 8
    7 CALL FREAD (EED,IREG,7,0)
      IF (IREG(6,1) .NE. -1) GO TO 7
      GO TO 6
    8 JREG = 1
      EPS  = .0001
      IF (REG(1,2) .NE. 0.) EPS = REG(1,2)
   11 CALL FREAD (EED,IREG(1,JREG),7,0)
      IF (IREG(6,JREG) .EQ. -1) GO TO 9
      JREG = JREG + 1
      IF (JREG .GT. 10) GO TO 9
      GO TO 11
    9 CALL CLOSE (EED,REW)
      NOREG = JREG - 1
      JREG  = 0
C
C     PICK UP PARAMETERS FOR REGION I
C
    5 JREG   = JREG + 1
      ITER   = 0
      NODCMP = 0
      NOSTRT = 0
      NOMOVS = 0
      X1     = REG(1,JREG)
      Y1     = REG(2,JREG)
      X2     = REG(3,JREG)
      Y2     = REG(4,JREG)
      L      = REG(5,JREG)
      ANOEST = REG(6,JREG)
      ANODES = REG(7,JREG)
      IF (NODES.EQ. 0) NODES = 3*NOEST
      NSHIFT = SQRT((X1-X2)**2+(Y1-Y2)**2)/L + 1.
      L1     = L*.5
      NOROOT = 0
C
C
C     FIND SHIFT POINT CLOSEST TO THE ORIGIN
C
      R = SQRT((X1-X2)**2 + (Y1-Y2)**2)
      IF (R) 10,10,15
   10 WRITE  (NOUT,12) UFM
   12 FORMAT (A23,' 2366, REGION IMPROPERLY DEFINED ON EIGC CARD.')
      CALL MESAGE (-61,0,0)
   15 CONTINUE
      D = (FLOAT(NSHIFT)*L-R)/2.0
      XX = X1 + D*(X1-X2)/R
      X2 = X2 + D*(X2-X1)/R
      X1 = XX
      YY = Y1 + D*(Y1-Y2)/R
      Y2 = Y2 + D*(Y2-Y1)/R
      Y1 = YY
      IF (IDIAG .EQ. 0) GO TO 7000
      WRITE  (NOUT,1000) X1,Y1,X2,Y2,L1,NODES,NOEST,NSHIFT
 1000 FORMAT (1H1,5F10.2,3I5)
 7000 CONTINUE
      DELTX = (X1-X2)/FLOAT(NSHIFT)
      DELTY = (Y1-Y2)/FLOAT(NSHIFT)
      XX    = X2 + DELTX/2.
      YY    = Y2 + DELTY/2.
      RANGE = XX**2  + YY**2
      N     = NSHIFT - 1
      SHIFT = 1.
      IF (DELTX .NE. 0.) GO TO 20
      ANUM1 = L1
      ANUM2 = 0.
      GO TO 25
   20 SLOPE = DELTY/DELTX
      ARG   = SQRT(1.+SLOPE**2)
      ANUM1 = SLOPE*L1/ARG
      ANUM2 = L1/ARG
   25 CONTINUE
      IF (N .EQ. 0) GO TO 40
      DO 30 I = 1,N
      XX    = XX + DELTX
      YY    = YY + DELTY
      RANG  = XX**2 + YY**2
      IF (RANG .GE. RANGE) GO TO 40
      RANGE = RANG
   30 SHIFT = I + 1
C
C     COMPUTE COORDINATES OF CORNERS OF THE REGION
C
   40 XL2   = X2 + ANUM1
      YL2   = Y2 - ANUM2
      IMIN  = SHIFT
      IMAX  = SHIFT
C
C     FIND THE MAXIMUM MODULUS OF THE SEARCH REGION
C
      MAXMOD = XL2**2 + YL2**2
      XX     = X2 - ANUM1
      YY     = Y2 + ANUM2
      MAXMOD = AMAX1(MAXMOD,XX**2+YY**2)
      XX     = X1 + ANUM1
      YY     = Y1 - ANUM2
      MAXMOD = AMAX1(MAXMOD,XX**2+YY**2)
      XX     = X1 - ANUM1
      YY     = Y1 + ANUM2
      MAXMOD = AMAX1(MAXMOD,XX**2+YY**2)
C
C     INITIALIZE
C
      IND    = 0
      LEFT   = 0
   45 ISHIFT = SHIFT
C
C     EVALUATE THE VALUE OF LAMBDA IN THE CENTER OF THE CURRENT SEARCH
C     REGION
C
      LAMBDA(1) = X2 + (SHIFT-.5)*DELTX
      LAMBDA(2) = Y2 + (SHIFT-.5)*DELTY
      IF (LAMBDA(2) .EQ. 0.0D0) LAMBDA(2) = .01*DELTY
C
C     COMPUTE DISTANCE TO FARTHEST CORNER OF THE SQUARE SEARCH REGION
C
      XX = XL2 + SHIFT*DELTX
      YY = YL2 + SHIFT*DELTY
      RZERO = (LAMBDA(1)-XX)**2 + (LAMBDA(2)-YY)**2
      RZERO = SQRT(RZERO)*1.05
      IF (IDIAG .EQ. 0) GO TO 7001
      WRITE  (NOUT,1216)RZERO
 1216 FORMAT (//,10H RZERO =     ,F10.4)
 7001 CONTINUE
      NOSTRT = NOSTRT + 1
      COMFLG = 0
   61 LMBDA(1) = LAMBDA(1)
      LMBDA(2) = LAMBDA(2)
      NOCHNG = 0
      SWITCH = 0
      IVECT  = 0
      KREG   = 0
      IND    = IND + 1
      IF (IABS (IND) .EQ. 13) IND = 1
      ISING  = 0
      GO TO 100
   80 ISING  = 0
      SWITCH = 1
  100 IF (NOCHNG .GE. 4) GO TO 220
      NOCHNG = NOCHNG + 1
      CALL KLOCK (T1)
C
C     CALL IN ADD LINK TO FORM (LAMBDA**2*M + LAMBDA*B + K)
C
      CALL CINVP1
C
C     CALL IN CD COMP TO DECOMPOSE THE MATRIX
C
      IF (IDIAG .EQ. 0) GO TO 7002
      WRITE  (NOUT,1001) LAMBDA
 1001 FORMAT (10H1LAMBDA =   ,2D15.5)
 7002 CONTINUE
      NODCMP = NODCMP + 1
      CALL CINVP2 (*110)
      CALL KLOCK (T2)
      GO TO 120
  110 IF (ISING .EQ. 1) GO TO 210
C
C     SINGULAR MATRIX. INCREMENT LAMBDA AND TRY ONCE MORE
C
      ISING = 1
      LAMBDA(1) = LAMBDA(1) + .02*RZERO
      LAMBDA(2) = LAMBDA(2) + .02*RZERO
      GO TO 100
C
C     DETERMINE THE TIME REQUIRED TO FORM AND DECOMPOSE THE DYNAMIC
C     MATRIX
C
  120 TIMED = T2 - T1
      IF (TIMED .EQ. 0) TIMED = 1
C
C     CALL IN MAIN LINK TO ITERATE FOR EIGENVALUES
C
  121 CALL CINVP3
      IF (LEFT   .EQ. 1) GO TO 130
      IF (COMFLG .EQ. 2) GO TO 125
      IF (COMFLG .EQ. 1) GO TO 80
      GO TO 140
  125 NOMOVS = NOMOVS + 1
      GO TO 61
C
C     CALL IN LINK TO COMPUTE THE LEFT EIGENVECTOR
C
  130 DTEMP(1)  = LAMBDA(1)
      DTEMP(2)  = LAMBDA(2)
      LAMBDA(1) = LAM1(1)
      LAMBDA(2) = LAM1(2)
  131 SWITCH    = -1
      CALL CINVP1
C
C     DECOMPOSE THE DYNAMIC MATRIX AT THE EIGENVALUE TO OBTAIN THE LEFT
C     EIGENVECTOR BY THE DETERMINATE METHOD
C
      IF (IDIAG .EQ. 0) GO TO 132
      WRITE (NOUT,1001) LAMBDA
  132 CALL CINVP2 (*138)
C
C     BUILD LOAD FOR FBS
C
      D1 = NROW/2
      D2 = NORTHO
      DO 133 I = 1,NROW,2
      K  = (I+1)/2
      DZ(I) = SIGN*MINDIA/(1.D0+(1.D0-FLOAT(K)/D1)*D2)
  133 DZ(I+1) = 0.0D0
      SIGN = -SIGN
      CALL CDIFBS (DZ(1),Z(IBUF))
      LAMBDA(1) = DTEMP(1)
      LAMBDA(2) = DTEMP(2)
      SWITCH = 0
C
C     NORMALIZE AND STORE THE LEFT EIGENVECTOR
C
      CALL CNORM1 (DZ(1),FILEK(2))
      IF (IDIAG .EQ. 0) GO TO 135
      WRITE  (NOUT,134) (DZ(I),I=1,NROW)
  134 FORMAT (///,15H LEFT VECTOR   ,//,(10D12.4))
  135 CONTINUE
      IF (NOLEFT .OR. ISYM.EQ.0) GO TO 136
      IFILE = PHIDLI
      CALL OPEN  (*500,IFILE,Z(IBUF),WRT)
      CALL WRITE (IFILE,DZ(1),NROW2,1)
      CALL CLOSE (IFILE,NOREW)
  136 IFILE = SCRFIL(10)
      CALL GOPEN (IFILE,Z(IBUF),RD)
      CALL BCKREC (IFILE)
      CALL FREAD (IFILE,DZ(NROW+2),NROW2,1)
      CALL BCKREC (IFILE)
      CALL CLOSE (IFILE,NOREW)
C
C     COMPUTE REAL LEFT VECTOR SCALED
C
      CALL CX TRN Y (DZ(1),DZ(NROW+2),DTEMP)
      CALL CDIVID (DZ(1),DZ(1),DTEMP,NROW)
      CALL OPEN (*500,IFILE,Z(IBUF),WRT)
      CALL WRITE (IFILE,DZ(1),NROW2,1)
      CALL CLOSE (IFILE,REW)
      GO TO 121
  138 LAMBDA(1) = 1.01*LAMBDA(1)
      LAMBDA(2) = 1.01*LAMBDA(2)
      GO TO 131
  140 IF (COMFLG .GE. 3) GO TO 200
      IF (COMFLG .EQ. 0) GO TO 160
      IF (IDIAG  .EQ. 0) GO TO 150
      WRITE  (NOUT,145) NOREG,JREG
  145 FORMAT (2I10)
  150 IF (NOREG .EQ. JREG) RETURN
      GO TO 5
C
C     FIND NEXT SHIFT POINT WHICH IS CLOSEST TO THE ORIGIN
C
  160 IF (IMIN .NE. 1) GO TO 170
      IF (IMAX .EQ. NSHIFT) GO TO 250
  165 SHIFT = SHIFT + 1.
      IMAX  = IMAX  + 1
      LAMBDA(1) = LMBDA(1) + DELTX
      LAMBDA(2) = LMBDA(2) + DELTY
      GO TO 45
  170 IF (IMAX .NE. NSHIFT) GO TO 180
  175 SHIFT = SHIFT - 1.
      IMIN  = IMIN  - 1
      LAMBDA(1) = LMBDA(1) - DELTX
      LAMBDA(2) = LMBDA(2) - DELTY
      GO TO 45
  180 XX   = LMBDA(1) - DELTX
      YY   = LMBDA(2) - DELTY
      RANG = XX**2 + YY**2
      XX   = LMBDA(1) + DELTX
      YY   = LMBDA(2) + DELTY
      RANGE= XX**2 + YY**2
      IF (RANGE-RANG) 175,175,165
  200 ITERM = COMFLG
      GO TO 260
C
C     SINGULARITY ENCOUNTERED TWICE IN A ROW
C
  210 ITERM = 1
      GO TO 260
C
C     4 MOVES WHILE TRACKING ONE ROOT
C
  220 ITERM = 2
      GO TO 260
C
C     REGIONS COMPLETED
C
  250 ITERM = 3
C
C     SET UP THE SUMMARY FILE
C
  260 IFILE = DMPFIL
      CALL OPEN  (*500,DMPFIL,Z(IBUF),WRT)
      CALL WRITE (DMPFIL,IHEAD(1),10,0)
      I       = 0
      IZ(I+2) = NORTHO
      IZ(I+3) = NOSTRT
      IZ(I+4) = NOMOVS
      IZ(I+5) = NODCMP
      IZ(I+6) = ITER
      IZ(I+7) = ITERM
      DO 270 I = 8,12
 270  IZ(I) = 0
      I = 2
      CALL WRITE (DMPFIL,IZ(I),40,0)
      CALL WRITE (DMPFIL,HEAD(1),96,1)
      CALL WRITE (DMPFIL,IZ(1),0,1)
      CALL CLOSE (DMPFIL,EOFNRW)
C
C     WRITE DUMMY TRAILER
      IXX = FILEK(1)
      FILEK(1) = DMPFIL
      CALL WRTTRL (FILEK(1))
      FILEK(1) = IXX
      NFOUND = NORTHO
      IF (IDIAG .EQ. 0) GO TO 350
      J = 12
      WRITE  (NOUT,300)(IZ(I),I=1,J)
  300 FORMAT (///,12I10)
  350 CONTINUE
      IF (ITERM .EQ. 5) RETURN
      GO TO 150
C
  500 CALL MESAGE (-1,IFILE,NAME)
      RETURN
      END


!-------------------------------------------------------------------------------
! Copyright (c) 2016 The University of Tokyo
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
!>   This module provide common functions of Plane deformation elements
module m_static_LIB_2d
  use hecmw, only : kint, kreal
  use elementInfo
  implicit none

contains
  !***********************************************************************
  !  2D Element:
  !  STF_C2   :Calculate stiff matrix of 2d elements
  !  DL_C2    :Deal with DLOAD conditions
  !  TLOAD_C2 :Deal with thermal expansion force
  !  UpdateST_C2  : Update Strain and stress
  !***********************************************************************
  !----------------------------------------------------------------------*
  subroutine STF_C2( ETYPE,NN,ecoord,gausses,PARAM1,stiff         &
      ,ISET, u)
    !----------------------------------------------------------------------*
    use mMechGauss
    use m_MatMatrix
    integer(kind=kint), intent(in) :: ETYPE
    integer(kind=kint), intent(in) :: NN
    type(tGaussStatus), intent(in) :: gausses(:)
    real(kind=kreal), intent(in)   :: ecoord(2,NN),PARAM1
    real(kind=kreal), intent(out)  :: stiff(:,:)
    integer(kind=kint),intent(in)  :: ISET
    real(kind=kreal), intent(in), optional :: u(:,:)

    ! LOCAL VARIABLES
    integer(kind=kint) :: flag
    integer(kind=kint) NDOF
    parameter(NDOF=2)
    real(kind=kreal) D(4,4),B(4,NDOF*NN),DB(4,NDOF*NN)
    real(kind=kreal) H(NN),stress(4)
    real(kind=kreal) THICK,PAI
    real(kind=kreal) DET,RR,WG
    real(kind=kreal) localcoord(2)
    real(kind=kreal) gderiv(NN,2), cdsys(3,3)
    integer(kind=kint) J,LX
    real(kind=kreal) gdispderiv(2,2)
    real(kind=kreal) B1(4,nn*ndof)
    real(kind=kreal) Smat(4,4)
    real(kind=kreal) BN(4,nn*ndof), SBN(4,nn*ndof)

    stiff =0.d0
    flag = gausses(1)%pMaterial%nlgeom_flag
    ! THICKNESS
    THICK=PARAM1
    !  FOR AX-SYM. ANALYSIS
    if(ISET==2) then
      THICK=1.d0
      PAI=4.d0*atan(1.d0)
    endif
    !* LOOP OVER ALL INTEGRATION POINTS
    do LX=1,NumOfQuadPoints( ETYPE )
      call MatlMatrix( gausses(LX), ISET, D, 1.d0, 1.d0,cdsys )
      if( .not. present(u) ) flag=INFINITE    ! enforce to infinite deformation analysis

      if( flag==1 .and. ISET == 2 ) then
        write(*,'(a)') '    PROGRAM STOP : non-TL element for axixsymmetric element'
        stop
      endif

      call getQuadPoint( ETYPE, LX, localcoord )
      call getGlobalDeriv( etype, nn, localcoord, ecoord, det, gderiv )
      !
      if(ISET==2) then
        call getShapeFunc( ETYPE, localcoord, H(:) )
        RR=dot_product( H(1:NN), ecoord(1,1:NN) )
        WG=getWeight( ETYPE, LX )*DET*RR*2.d0*PAI
      else
        RR=THICK
        H(:)=0.d0
        WG=getWeight( ETYPE, LX )*DET*RR
      end if
      do J=1,NN
        B(1,2*J-1)=gderiv(J,1)
        B(2,2*J-1)=0.d0
        B(3,2*J-1)=gderiv(J,2)
        B(1,2*J  )=0.d0
        B(2,2*J  )=gderiv(J,2)
        B(3,2*J  )=gderiv(J,1)
        B(4,2*J-1)=H(J)/RR
        B(4,2*J  )=0.d0
      enddo
      !
      ! ----- calculate the BL1 matrix ( TOTAL LAGRANGE METHOD )
      if( flag==1 ) then
        !       ----- dudx(i,j) ==> gdispderiv(i,j)
        gdispderiv(1:ndof, 1:ndof) = matmul( u(1:ndof,1:nn), gderiv(1:nn,1:ndof) )
        B1(1:4,1:nn*ndof) = 0.d0
        do j=1,nn
          B1(1,2*j-1) = gdispderiv(1,1)*gderiv(j,1)
          B1(2,2*j-1) = gdispderiv(1,2)*gderiv(j,2)
          B1(3,2*j-1) = gdispderiv(1,2)*gderiv(j,1)+gdispderiv(1,1)*gderiv(j,2)
          B1(1,2*j  ) = gdispderiv(2,1)*gderiv(j,1)
          B1(2,2*j  ) = gdispderiv(2,2)*gderiv(j,2)
          B1(3,2*j  ) = gdispderiv(2,2)*gderiv(j,1)+gdispderiv(2,1)*gderiv(j,2)
          B1(4,2*j-1) = 0.d0
          B1(4,2*j  ) = 0.d0
        enddo
        do j=1,nn*ndof
          B(:,j) = B(:,j)+B1(:,j)
        enddo
      endif
      !
      DB(1:4,1:nn*ndof) = 0.d0
      DB(1:4,1:nn*ndof) = matmul( D(1:4,1:4), B(1:4,1:nn*ndof) )
      stiff(1:nn*ndof,1:nn*ndof) = stiff(1:nn*ndof,1:nn*ndof) +             &
        matmul( transpose( B(1:4,1:nn*ndof) ), DB(1:4,1:nn*ndof) )*WG
      !
      !    ----- calculate the stress matrix ( TOTAL LAGRANGE METHOD )
      if( flag==1 ) then
        stress(1:4)=gausses(LX)%stress(1:4)
        BN(1:4,1:nn*ndof) = 0.d0
        do j=1,nn
          BN(1,2*j-1) = gderiv(j,1)
          BN(2,2*j  ) = gderiv(j,1)
          BN(3,2*j-1) = gderiv(j,2)
          BN(4,2*j  ) = gderiv(j,2)
        enddo
        Smat(:,:) = 0.d0
        do j=1,2
          Smat(j  ,j  ) = stress(1)
          Smat(j  ,j+2) = stress(3)
          Smat(j+2,j  ) = stress(3)
          Smat(j+2,j+2) = stress(2)
        enddo
        SBN(1:4,1:nn*ndof) = matmul( Smat(1:4,1:4), BN(1:4,1:nn*ndof) )
        stiff(1:nn*ndof,1:nn*ndof) = stiff(1:nn*ndof,1:nn*ndof) +           &
          matmul( transpose( BN(1:4,1:nn*ndof) ), SBN(1:4,1:nn*ndof) )*WG
      endif
      !
    enddo

  end subroutine STF_C2

  !----------------------------------------------------------------------*
  subroutine DL_C2(ETYPE,NN,XX,YY,RHO,PARAM1,LTYPE,PARAMS,VECT,NSIZE,ISET)
    !----------------------------------------------------------------------*
    !**
    !**  Deal with DLOAD conditions
    !**
    !   BX   LTYPE=1  :BODY FORCE IN X-DIRECTION
    !   BY   LTYPE=2  :BODY FORCE IN Y-DIRECTION
    !   GRAV LTYPE=4  :GRAVITY FORCE
    !   CENT LTYPE=5  :CENTRIFUGAL LOAD
    !   P1   LTYPE=10 :TRACTION IN NORMAL-DIRECTION FOR FACE-1
    !   P2   LTYPE=20 :TRACTION IN NORMAL-DIRECTION FOR FACE-2
    !   P3   LTYPE=30 :TRACTION IN NORMAL-DIRECTION FOR FACE-3
    !   P4   LTYPE=40 :TRACTION IN NORMAL-DIRECTION FOR FACE-4
    ! I/F VARIABLES
    integer(kind=kint), intent(in) :: ETYPE,NN
    real(kind=kreal), intent(in)   :: XX(NN),YY(NN),PARAMS(0:6)
    real(kind=kreal), intent(out)  :: VECT(NN*2)
    real(kind=kreal) RHO,PARAM1
    integer(kind=kint) LTYPE,NSIZE,ISET
    ! LOCAL VARIABLES
    integer(kind=kint), parameter :: NDOF=2
    real(kind=kreal) H(NN)
    real(kind=kreal) PLX(NN),PLY(NN)
    real(kind=kreal) XJ(2,2),DET,RR,WG
    real(kind=kreal) PAI,val,THICK
    integer(kind=kint) IVOL,ISUF
    real(kind=kreal) AX,AY,RX,RY
    real(kind=kreal) normal(2)
    real(kind=kreal) COEFX,COEFY,XCOD,YCOD,HX,HY,PHX,PHY
    integer(kind=kint) NOD(NN)
    integer(kind=kint) LX,I,SURTYPE,NSUR
    real(kind=kreal) VX,VY,localcoord(2),deriv(NN,2),elecoord(2,NN)
    !  CONSTANT
    PAI=4.d0*atan(1.d0)
    ! SET VALUE
    val=PARAMS(0)
    ! THICKNESS
    THICK=PARAM1
    ! CLEAR VECT
    NSIZE=NN*NDOF
    VECT(1:NSIZE)=0.d0
    !
    ! SELCTION OF LOAD TYPE
    !
    IVOL=0
    ISUF=0
    if( LTYPE.LT.10 ) then
      IVOL=1
      if( LTYPE.EQ.5 ) then
        AX=PARAMS(1)
        AY=PARAMS(2)
        RX=PARAMS(4)
        RY=PARAMS(5)
      endif
    else if( LTYPE.GE.10 ) then
      ISUF=1
      call getSubFace( ETYPE, LTYPE/10, SURTYPE, NOD )
      NSUR = getNumberOfNodes( SURTYPE )
    endif
    !** SURFACE LOAD
    if( ISUF==1 ) then
      do I=1,NSUR
        elecoord(1,i)=XX(NOD(I))
        elecoord(2,i)=YY(NOD(i))
      enddo
      !** INTEGRATION OVER SURFACE
      do LX=1,NumOfQuadPoints( SURTYPE )
        call getQuadPoint( SURTYPE, LX, localcoord(1:1) )
        call getShapeFunc( SURTYPE, localcoord(1:1), H(1:NSUR) )
        normal=EdgeNormal( SURTYPE, NSUR, localcoord(1:1), elecoord(:,1:NSUR) )
        WG = getWeight( SURTYPE, LX )
        if( ISET==2 ) then
          RR=0.d0
          do I=1,NSUR
            RR=RR+H(I)*XX(NOD(I))
          enddo
          WG=WG*RR*2.d0*PAI
        else
          WG=WG*THICK
        endif
        do I=1,NSUR
          VECT(2*NOD(I)-1)=VECT(2*NOD(I)-1)+val*WG*H(I)*normal(1)
          VECT(2*NOD(I)  )=VECT(2*NOD(I)  )+val*WG*H(I)*normal(2)
        enddo
      enddo
    endif
    !** VOLUME LOAD
    if( IVOL==1 ) then
      PLX(:)=0.d0
      PLY(:)=0.d0
      do LX=1,NumOfQuadPoints( ETYPE )
        call getQuadPoint( ETYPE, LX, localcoord(1:2) )
        call getShapeDeriv( ETYPE, localcoord(1:2), deriv )
        call getShapeFunc( ETYPE, localcoord(1:2), H(1:NN) )
        XJ(1,1:2)=matmul( XX(1:NN), deriv(1:NN,1:2) )
        XJ(2,1:2)=matmul( YY(1:NN), deriv(1:NN,1:2) )

        DET=XJ(1,1)*XJ(2,2)-XJ(2,1)*XJ(1,2)

        WG = getWeight( ETYPE, LX )
        if(ISET==2) then
          RR=dot_product( H(1:NN),XX(1:NN) )
          WG=WG*DET*RR*2.d0*PAI
        else
          RR=THICK
          WG=WG*DET*RR
        end if
        COEFX=1.d0
        COEFY=1.d0
        ! CENTRIFUGAL LOAD
        if( LTYPE==5 ) then
          XCOD=dot_product( H(1:NN),XX(1:NN) )
          YCOD=dot_product( H(1:NN),YY(1:NN) )
          HX=AX + ( (XCOD-AX)*RX+(YCOD-AY)*RY )/(RX**2+RY**2)*RX
          HY=AY + ( (XCOD-AX)*RX+(YCOD-AY)*RY )/(RX**2+RY**2)*RY
          PHX=XCOD-HX
          PHY=YCOD-HY
          COEFX=RHO*val*val*PHX
          COEFY=RHO*val*val*PHY
        end if
        do I=1,NN
          PLX(I)=PLX(I)+H(I)*WG*COEFX
          PLY(I)=PLY(I)+H(I)*WG*COEFY
        enddo
      enddo
      if( LTYPE.EQ.1) then
        do I=1,NN
          VECT(2*I-1)=val*PLX(I)
          VECT(2*I  )=0.d0
        enddo
      else if( LTYPE.EQ.2 ) then
        do I=1,NN
          VECT(2*I-1)=0.d0
          VECT(2*I  )=val*PLY(I)
        enddo
      else if( LTYPE.EQ.4 ) then
        VX=PARAMS(1)
        VY=PARAMS(2)
        VX=VX/sqrt( PARAMS(1)**2 + PARAMS(2)**2 )
        VY=VY/sqrt( PARAMS(1)**2 + PARAMS(2)**2 )
        do I=1,NN
          VECT(2*I-1)=val*PLX(I)*RHO*VX
          VECT(2*I  )=val*PLY(I)*RHO*VY
        enddo
      else if( LTYPE==5 ) then
        do I=1,NN
          VECT(2*I-1)=PLX(I)
          VECT(2*I  )=PLY(I)
        enddo
      end if
    endif
  end subroutine DL_C2

  !----------------------------------------------------------------------*
  subroutine TLOAD_C2(ETYPE,NN,XX,YY,TT,T0,gausses,PARAM1,ISET,VECT)
    !----------------------------------------------------------------------*
    !
    ! THERMAL LOAD OF PLANE ELEMENT
    !
    use mMaterial
    use m_MatMatrix
    use m_ElasticLinear
    ! I/F VARIABLES
    integer(kind=kint), intent(in) :: ETYPE,NN
    real(kind=kreal),intent(in)    :: XX(NN),YY(NN),TT(NN),T0(NN),PARAM1
    type(tGaussStatus), intent(in) :: gausses(:)          !< status of qudrature points
    real(kind=kreal),intent(out)   :: VECT(NN*2)
    integer(kind=kint) ISET
    ! LOCAL VARIABLES
    integer(kind=kint) NDOF
    parameter(NDOF=2)
    real(kind=kreal) ALP,PP,D(4,4),B(4,NDOF*NN)
    real(kind=kreal) H(NN)
    real(kind=kreal) EPS(4),SGM(4),localcoord(2)
    real(kind=kreal) deriv(NN,2),DNDE(2,NN)
    real(kind=kreal) XJ(2,2),DET,RR,DUM
    real(kind=kreal) XJI(2,2)
    real(kind=kreal) PAI,THICK,WG
    integer(kind=kint) J,LX
    real(kind=kreal) TEMPC,TEMP0,THERMAL_EPS
    !*************************
    !  CONSTANT
    !*************************
    PAI=4.d0*atan(1.d0)
    ! CLEAR VECT
    VECT(1:NN*NDOF)=0.d0
    ! THICKNESS
    THICK=PARAM1
    !  FOR AX-SYM. ANALYSIS
    if(ISET==2) THICK=1.d0
    ! We suppose material properties doesn't varies inside element
    call calElasticMatrix( gausses(1)%pMaterial, ISET, D  )
    ALP = gausses(1)%pMaterial%variables(M_EXAPNSION)
    pp = gausses(1)%pMaterial%variables(M_POISSON)
    !* LOOP OVER ALL INTEGRATION POINTS
    do LX=1,NumOfQuadPoints( ETYPE )
      call getQuadPoint( ETYPE, LX, localcoord )
      call getShapeFunc( ETYPE, localcoord, H(1:NN) )
      call getShapeDeriv( ETYPE, localcoord, deriv(:,:) )
      XJ(1,1:2)=matmul( XX(1:NN), deriv(1:NN,1:2) )
      XJ(2,1:2)=matmul( YY(1:NN), deriv(1:NN,1:2) )

      DET=XJ(1,1)*XJ(2,2)-XJ(2,1)*XJ(1,2)

      TEMPC=dot_product(H(1:NN),TT(1:NN))
      TEMP0=dot_product(H(1:NN),T0(1:NN))
      if(ISET==2) then
        RR=dot_product(H(1:NN),XX(1:NN))
        WG=getWeight( ETYPE, LX )*DET*RR*2.d0*PAI
      else
        RR=THICK
        WG=getWeight( ETYPE, LX )*DET*RR
      end if
      DUM=1.d0/DET
      XJI(1,1)= XJ(2,2)*DUM
      XJI(1,2)=-XJ(2,1)*DUM
      XJI(2,1)=-XJ(1,2)*DUM
      XJI(2,2)= XJ(1,1)*DUM

      DNDE=matmul( XJI, transpose(deriv) )
      do J=1,NN
        B(1,2*J-1)=DNDE(1,J)
        B(2,2*J-1)=0.d0
        B(3,2*J-1)=DNDE(2,J)
        B(1,2*J  )=0.d0
        B(2,2*J  )=DNDE(2,J)
        B(3,2*J  )=DNDE(1,J)
        B(4,2*J-1)=0.d0
        B(4,2*J  )=0.d0
        if( ISET==2 ) then
          B(4,2*J-1)=H(J)/RR
        endif
      enddo
      !**
      !** THERMAL EPS
      !**
      THERMAL_EPS=ALP*(TEMPC-TEMP0)
      if( ISET .EQ. 2 ) then
        EPS(1)=THERMAL_EPS
        EPS(2)=THERMAL_EPS
        EPS(3)=0.d0
        EPS(4)=THERMAL_EPS
      else if( ISET.EQ.0 ) then
        EPS(1)=THERMAL_EPS*(1.d0+PP)
        EPS(2)=THERMAL_EPS*(1.d0+PP)
        EPS(3)=0.d0
        EPS(4)=0.d0
      else if( ISET.EQ.1 ) then
        EPS(1)=THERMAL_EPS
        EPS(2)=THERMAL_EPS
        EPS(3)=0.d0
        EPS(4)=0.d0
      endif
      !**
      !** SET SGM  {s}=[D]{e}
      !**
      SGM = matmul( EPS(1:4), D(1:4,1:4) )
      !**
      !** CALCULATE LOAD {F}=[B]T{e}
      !**
      VECT(1:NN*NDOF)=VECT(1:NN*NDOF)+matmul( SGM(1:4), B(1:4,1:NN*NDOF) )*WG
    enddo

  end subroutine TLOAD_C2
  !
  !
  !> Update strain and stress inside element
  !---------------------------------------------------------------------*
  subroutine UPDATE_C2( etype,nn,ecoord,gausses,PARAM1,iset,         &
      u, ddu, qf, TT, T0, TN )
    !---------------------------------------------------------------------*
    use m_fstr
    use mMechGauss
    use m_MatMatrix
    ! I/F VARIAVLES
    integer(kind=kint), intent(in)     :: etype, nn
    real(kind=kreal),   intent(in)     :: ecoord(3,nn),PARAM1
    integer(kind=kint), intent(in)     :: ISET
    type(tGaussStatus), intent(inout)  :: gausses(:)
    real(kind=kreal),   intent(in)     :: u(2,nn)
    real(kind=kreal),   intent(in)     :: ddu(2,nn)
    real(kind=kreal),   intent(out)    :: qf(:)
    real(kind=kreal), intent(in), optional :: TT(nn)   !< current temperature
    real(kind=kreal), intent(in), optional :: T0(nn)   !< reference temperature
    real(kind=kreal), intent(in), optional :: TN(nn)   !< reference temperature

    integer(kind=kint), parameter :: ndof=2
    real(kind=kreal)   :: D(4,4), B(4,ndof*nn), B1(4,ndof*nn)
    real(kind=kreal)   :: H(nn)
    real(kind=kreal)   :: thick,pai,rr
    real(kind=kreal)   :: det, WG
    real(kind=kreal)   :: localCoord(2)
    real(kind=kreal)   :: gderiv(nn,2), dstrain(4), ttc, tt0, ttn
    real(kind=kreal)   :: gdispderiv(2,2), cdsys(3,3)
    integer(kind=kint) :: j, LX
    real(kind=kreal)   :: totaldisp(2*nn)
    real(kind=kreal)   :: EPSTH(4), alp, alp0, ina(1), outa(1)
    logical            :: ierr
    !
    qf(:)              = 0.d0
    do j=1,nn
      totaldisp((j-1)*2+1) = u(1,j) + ddu(1,j)
      totaldisp(j*2)       = u(2,j) + ddu(2,j)
    enddo
    !
    ! THICKNESS
    THICK=PARAM1
    !  FOR AX-SYM. ANALYSIS
    if(ISET==2) then
      THICK=1.d0
      PAI=4.d0*atan(1.d0)
    endif
    !* LOOP OVER ALL INTEGRATION POINTS
    do LX=1, NumOfQuadPoints(etype)
      call MatlMatrix( gausses(LX), ISET, D, 1.d0, 1.d0, cdsys  )

      call getQuadPoint( etype, LX, localCoord(:) )
      call getGlobalDeriv( etype, nn, localcoord, ecoord, det, gderiv )
      !
      EPSTH = 0.d0
      if( present(TT) .AND. present(T0) ) then
        call getShapeFunc( ETYPE, localcoord, H(:) )
        ttc = dot_product(TT(:), H(:))
        tt0 = dot_product(T0(:), H(:))
        ttn = dot_product(TN(:), H(:))
        ina(1) = ttc
        call fetch_TableData( MC_THEMOEXP, gausses(LX)%pMaterial%dict, outa(:), ierr, ina )
        if( ierr ) outa(1) = gausses(LX)%pMaterial%variables(M_EXAPNSION)
        alp = outa(1)
        ina(1) = tt0
        call fetch_TableData( MC_THEMOEXP, gausses(LX)%pMaterial%dict, outa(:), ierr, ina )
        if( ierr ) outa(1) = gausses(LX)%pMaterial%variables(M_EXAPNSION)
        alp0 = outa(1)
        EPSTH(1:2)=alp*(ttc-ref_temp)-alp0*(tt0-ref_temp)
      end if
      !
      WG=getWeight( etype, LX )*DET
      if(ISET==2) then
        call getShapeFunc( ETYPE, localcoord, H(:) )
        RR=dot_product( H(1:NN), ecoord(1,1:NN) )
      else
        RR=THICK
        H(:)=0.d0
      end if
      do J=1,NN
        B(1,2*J-1)=gderiv(J,1)
        B(2,2*J-1)=0.d0
        B(3,2*J-1)=gderiv(J,2)
        B(1,2*J  )=0.d0
        B(2,2*J  )=gderiv(J,2)
        B(3,2*J  )=gderiv(J,1)
        B(4,2*J-1)=H(J)/RR
        B(4,2*J  )=0.d0
      enddo

      gausses(LX)%strain(1:4) = matmul( B(:,:), totaldisp )
      gausses(LX)%stress(1:4) = matmul( D(1:4, 1:4), gausses(LX)%strain(1:4)-EPSTH(1:4) )

      !set stress and strain for output
      gausses(LX)%strain_out(1:4) = gausses(LX)%strain(1:4)
      gausses(LX)%stress_out(1:4) = gausses(LX)%stress(1:4)

      !
      !    ----- calculate the Internal Force
      qf(1:nn*ndof) = qf(1:nn*ndof) +                                     &
        matmul( gausses(LX)%stress(1:4), B(1:4,1:nn*ndof) )*WG
      !
    enddo
    !
  end subroutine UPDATE_C2
  !
  !----------------------------------------------------------------------*
  subroutine NodalStress_C2(ETYPE,NN,gausses,ndstrain,ndstress)
    !----------------------------------------------------------------------*
    !
    ! Calculate Strain and Stress increment of solid elements
    !
    use mMechGauss

    integer(kind=kint), intent(in) :: ETYPE,NN
    type(tGaussStatus), intent(in) :: gausses(:)
    real(kind=kreal), intent(out)  :: ndstrain(NN,4)
    real(kind=kreal), intent(out)  :: ndstress(NN,4)

    integer          :: i,ic
    real(kind=kreal) :: TEMP(8)

    TEMP(:)=0.d0
    IC = NumOfQuadPoints(etype)
    do i=1,IC
      TEMP(1:4) = TEMP(1:4) + gausses(i)%strain_out(1:4)
      TEMP(5:8) = TEMP(5:8) + gausses(i)%stress_out(1:4)
    enddo
    TEMP(1:8) = TEMP(1:8)/IC
    forall( i=1:NN )
      ndstrain(i,1:4) = TEMP(1:4)
      ndstress(i,1:4) = TEMP(5:8)
    end forall

  end subroutine

  !----------------------------------------------------------------------*
  subroutine ElementStress_C2(ETYPE,gausses,strain,stress)
    !----------------------------------------------------------------------*
    !
    ! Calculate Strain and Stress increment of solid elements
    !
    use mMechGauss
    integer(kind=kint), intent(in) :: ETYPE
    type(tGaussStatus), intent(in) :: gausses(:)
    real(kind=kreal), intent(out)  :: strain(4)
    real(kind=kreal), intent(out)  :: stress(4)

    integer          :: i,ic

    strain(:)=0.d0; stress(:)=0.d0
    IC = NumOfQuadPoints(etype)
    do i=1,IC
      strain(:) = strain(:) + gausses(i)%strain_out(1:4)
      stress(:) = stress(:) + gausses(i)%stress_out(1:4)
    enddo
    strain(:) = strain(:)/IC
    stress(:) = stress(:)/IC

  end subroutine

end module m_static_LIB_2d


module io_manager_base
  use netcdf_layer_base
  implicit none

  type :: IOManager_
     logical :: parallel = .false.
     class(NetCDF_layer_),allocatable :: netcdf_layer
   contains
  end type IOManager_

  interface IOManager_
     module procedure IOManager_init
  end interface IOManager_

contains

  type(IOManager_) function IOManager_init(parallel)
    implicit none

    logical, optional, intent(in) :: parallel

    if(.not.present(parallel)) then
       IOManager_init%parallel = .false.
    else
       IOManager_init%parallel = parallel
    end if

    if( IOManager_init%parallel .eqv. .false.) then
       allocate(NetCDF_serial_ :: IOManager_init%netcdf_layer)
    else
       allocate(NetCDF_parallel_ :: IOManager_init%netcdf_layer)
    end if

  end function IOManager_init

end module io_manager_base


text:
program hello_images_ordered
  ! Synchronizing parallel images to print in order.
  implicit none
  integer :: n
  do n = 1, num_images()
    if (this_image() == n) &
      print *, 'Hello from image', this_image(), 'of', num_images()
    sync all
  end do
end program hello_images_ordered

max_stars_repo_name: wilsonify/c-consumer

max_issues_repo_name: wilsonify/c-consumer

max_forks_repo_name: wilsonify/c-consumer


text:
      subroutine tdrwdt(kv,y,jp,ymin,ymax,ls,le,frs,fre,ale,np,
     $     only,monly,patt,
     1     latt,twiss,idp,pos,imon,emon,nmon)
      use tfstk
      use ffs
      use tffitcode
      use ffs_pointer, only:idvalc,idtypec,pnamec,lpnamec,gammab
      use mathfun
      implicit real*8 (a-h,o-z)
      logical*4 sol,over,fin,monel,only,monly,tmatch
      character*(*) patt
      real*8 y(2,np),twiss(nlat,-ndim:ndim,ntwissfun),pos(nlat),
     $     emon(nmona,4),twisss(ntwissfun),ymin,ymax,
     $     frs,fre,hbs
      integer*8 latt(nlat),ip
      integer*4 kv,ls,le,imon(nmona,4)
      real*8 beam(42),sv(5),vsave(100),ftwiss(ntwissfun)
      include 'inc/common.inc'
      theta=0d0
      altotal=pos(nlat)-pos(1)
      alen=pos(le)-pos(ls)
      if(alen .le. 0.d0)then
        alen=alen+altotal
      endif
      jp=0
      sol=.false.
      istrt=ls
      if(ls. gt. le)then
        istop=nlat-1
      else
        istop=le
      endif
      pa=0.d0
      fin=.false.
 111  LOOP_I: do i=istrt,istop-1
        if(only .and. .not. tmatch(pname(ilist(2,latt(i))),patt))then
        elseif(monly .and. .not. monel(imon,nmon,i,m1))then
        elseif(kv .gt. 32 .and. kv .le. 34) then
          if(monel(imon,nmon,i,m1)) then
            jp=jp+1
            y(1,jp)=pos(i)+pa
            ll=imon(m1,2)
            nq=imon(ll,4)
            y(2,jp)=twiss(i,idp,2*kv-51)-twiss(i,ndim,2*kv-51)
     1           -rlist(latt(nq)+kv-28)-emon(ll,kv-32)
          endif
        else
          id=idtypec(i)
          if(id .ge. icMARK)then
            cycle LOOP_I
          endif
          if(i .eq. istrt .and. frs .ne. 0.d0)then
            f=1.d0-frs
            dpa=(pos(i+1)-pos(i))*f
          elseif(i .eq. istop-1 .and. fre .ne. 0.d0)then
            f=fre
            dpa=0.d0
          else
            f=1.d0
            dpa=0.d0
          endif
          if(kv .le. 27 .or. kv .eq. 35)then
            jp=jp+1
            y(1,jp)=pos(i)+pa+dpa
            y(2,jp)=drwfun(kv,i,twiss,gammab,idp,amass,nlat,ndim)
          elseif(kv .le. 32) then
            jp=jp+1
            y(1,jp)=pos(i)+pa+dpa
            if(trsize)then
              call tftrb(latt,i,sv)
              beam(1)=sv(3)
              beam(4)=sv(5)
              beam(6)=sv(4)
            else
              call tfbeam(i,theta,beam)
            endif
            if(kv.ne.32) then
              ia=(kv-27)*(kv-26)/2
              y(2,jp)=sqrt(beam(ia))
            else
              y(2,jp)=atan2(-2d0*beam(4),beam(1)-beam(6)) /2d0
            endif
          endif
          ip=latt(i)
          if(id .le. icdodeca .or. id .eq. icmult
     $         .or. id .eq. iccavi .or. id .eq. icTCAV)then
            al=rlist(ip+kytbl(kwL,id))
            ndiv=int(al*f/ale)
            if(ndiv .gt. 0)then
              nkey=kytbl(kwmax,id)-1
              call tmov(rlist(ip+1),vsave,nkey)
              gbs=gammab(i+1)
              hs=p2h(gbs)
c              hs=gbs*sqrt(1.d0+1.d0/gbs**2)
              ha=p2h(gammab(i))
c              ha=gammab(i)*sqrt(1.d0+1.d0/gammab(i)**2)
              dh=hs-ha
              do k=1,ntwissfun
                twisss(k)=twiss(i+1,idp,k)
              enddo
              if(kv .gt. 27 .and. kv .le. 32 .and. trsize)then
                call tracke(latt,i,sv,np,'STANDBY',' ',0)
                call tracke(latt,i,sv,np,'TRACK',' ',0)
                call tracke(latt,i,sv,np,'SAVE',' ',0)
              endif
              do n=1,ndiv
                rn=dble(n)/(ndiv+1)
                call qtwissfrac(ftwiss,i,rn*f,over)
                do ii=1,ntwissfun
                  twiss(i+1,idp,ii)=ftwiss(ii)
                enddo
                jp=jp+1
                y(1,jp)=pos(i)+pa+dpa+al*rn*f
                if(kv .le. 27 .or. kv .eq. 35) then
                  gammab(i+1)=h2p(rn*f*dh+ha)
c                  gammab(i+1)=sqrt((rn*f*dh+ha)**2-1.d0)
                  y(2,jp)=drwfun(kv,i+1,twiss,gammab,
     $                 idp,amass,nlat,ndim)
                  gammab(i+1)=gbs
                elseif(kv.le.32)then
                  if(trsize)then
c                    call tpara(ilist(2,latt(i)))
c                    call tracke(latt,i+1,sv,np,'TRACK',' ',0)
c                    beam(1)=sv(3)
c                    beam(4)=sv(5)
c                    beam(6)=sv(4)
c                    if(ilist(2,latt(i)) .ne. 0)then
c                      call tfree(ilist(2,latt(i)))
c                      ilist(2,latt(i))=0
c                    endif
                  else
                    hbs=rn*f*dh+ha
                    gammab(i+1)=h2p(hbs)
c                    gammab(i+1)=hbs*sqrt(1.d0-1.d0/hbs**2)
                    call tfbeam(i+1,theta,beam)
                    gammab(i+1)=gbs
                  endif
                  if(kv.ne.32) then
                    ia=(kv-27)*(kv-26)/2
                    y(2,jp)=sqrt(beam(ia))
                  else
                    y(2,jp)=atan2(-2d0*beam(4),beam(1)-beam(6)) /2d0
                  endif
                endif
              enddo
              if(kv .gt. 27 .and. kv .le. 32 .and. trsize)then
                call tracke(latt,i+1,sv,np,'RESET',' ',0)
              endif
              call tmov(vsave,rlist(ip+1),nkey)
              gammab(i+1)=gbs
              do k=1,ntwissfun
                twiss(i+1,idp,k)=twisss(k)
              enddo
            endif
          endif
        endif
c        write(*,*)'tdrwdt ',i,id,jp
      enddo LOOP_I
      i=istop
      if(fre .ne. 0.d0)then
        if(i .eq. nlat-1)then
          dpa=0.d0
        else
          dpa=(pos(i)-pos(i-1))*(fre-1.d0)
        endif
      else
        dpa=0.d0
      endif
      if(only .and. .not. tmatch(pnamec(i),patt))then
      elseif(monly .and. .not. monel(imon,nmon,i,m1))then
      else
        if(kv .le. 27 .or. kv .eq. 35)then
          jp=jp+1
          y(1,jp)=pos(i)+pa+dpa
          y(2,jp)=drwfun(kv,i,twiss,gammab,idp,amass,nlat,ndim)
        elseif(kv .le. 32) then
          jp=jp+1
          y(1,jp)=pos(i)+pa+dpa
          if(trsize)then
            call tftrb(latt,i,sv)
            beam(1)=sv(3)
            beam(4)=sv(5)
            beam(6)=sv(4)
          else
            call tfbeam(i,theta,beam)
          endif
          if(kv.ne.32) then
            ia=(kv-27)*(kv-26)/2
            y(2,jp)=sqrt(beam(ia))
          else
            y(2,jp)=atan2(-2d0*beam(4),beam(1)-beam(6)) /2d0
          endif
        elseif(kv.le.34) then
          if(monel(imon,nmon,i,m1)) then
            jp=jp+1
            y(1,jp)=pos(i)+pa+dpa
            ll=imon(m1,2)
            nq=imon(ll,4)
            y(2,jp)=twiss(i,idp,2*kv-51)-twiss(i,ndim,2*kv-51)
     1           -rlist(latt(nq)+kv-28)-emon(ll,kv-32)
          endif
        endif
      endif
      if(ls .gt. le .and. .not. fin) then
        pa=pos(nlat)-pos(1)
        istrt=1
        istop=le
        fin=.true.
        goto 111
      endif
c      write(*,*)'tdrwdt ',jp,np
      ymin= 1d20
      ymax=-1d20
      do i=1,jp
        ymin=min(y(2,i),ymin)
        ymax=max(y(2,i),ymax)
      enddo
      if(kv .eq. mfitbx .or. kv .eq. mfitby
     $     .or. (kv .ge. 28 .and .kv .le. 31))then
        ymin=0d0
      endif
      return
      end

      real*8 function drwfun(kv,i,twiss,gammab,idp,
     $     amass,nlat,ndim)
      use tffitcode
      implicit none
      integer*4 kv,i,idp,nlat,ndim
      real*8 twiss(nlat,-ndim:ndim,ntwissfun),gammab(nlat),
     $     amass,u(4),v(4)
      if(kv .le. mfitdpy) then
        drwfun=twiss(i,idp,kv)
      elseif(kv .le. 22) then
c        u(1)=twiss(i,idp,mfitex)
c        u(2)=twiss(i,idp,mfitepx)
c        u(3)=twiss(i,idp,mfitey)
c        u(4)=twiss(i,idp,mfitepy)
c        call mc2to4(twiss,idp,i,u,v)
        call tgetphysdisp(i,v)
        drwfun=v(kv-18)
      elseif(kv .eq. 23) then
        drwfun=twiss(i,idp,mfitr1)*twiss(i,idp,mfitr4)
     z       -twiss(i,idp,mfitr2)*twiss(i,idp,mfitr3)
      elseif(kv .le. 27) then
        u(1)=twiss(i,idp,mfitdx)
        u(2)=twiss(i,idp,mfitdpx)
        u(3)=twiss(i,idp,mfitdy)
        u(4)=twiss(i,idp,mfitdpy)
        call mc4to2(twiss,idp,i,u,v)
        drwfun=v(kv-23)
      elseif(kv .eq. 35) then
        drwfun=amass*gammab(i)
      else
        drwfun=0.d0
      endif
      return
      end

max_stars_repo_name: noboruatkek/SAD

max_issues_repo_name: noboruatkek/SAD

max_forks_repo_name: noboruatkek/SAD


text:
      module comrlp
!***********************************************************************
! Copyright 2008 Los Alamos National Security, LLC  All rights reserved
! Unless otherwise indicated,  this information has been authored by an
! employee or employees of the Los Alamos National Security, LLC (LANS),
! operator of the  Los  Alamos National  Laboratory  under Contract  No.
! DE-AC52-06NA25396  with  the U. S. Department  of  Energy.  The  U. S.
! Government   has   rights  to  use,  reproduce,  and  distribute  this
! information.  The  public may copy  and  use this  information without
! charge, provided that this  Notice and any statement of authorship are
! reproduced on all copies.  Neither  the  Government nor LANS makes any
! warranty,   express   or   implied,   or   assumes  any  liability  or
! responsibility for the use of this information.
!***********************************************************************
!D1
!D1 PURPOSE
!D1
!D1 Global file for relative permeability and capillary pressure.
!D1
!***********************************************************************

      integer :: maxrp = 0, maxcp = 0, ntable =0, ntblines=0
      integer :: max_rp, max_rpf, max_cp
      integer :: ishisrlp, num_sat,phase_comp(30,3),nphases(10)
      integer, allocatable :: cap_coupling(:,:)
      integer, allocatable :: cap_pos(:,:)
      integer, allocatable :: cap_type(:), cap_type2(:,:)
      integer, allocatable :: rlp_group(:)
      integer, allocatable :: rlp_phase(:,:)
      integer, allocatable :: rlp_pos(:,:)
      integer, allocatable :: rlp_type(:), rlp_type2(:,:)
      integer, allocatable :: tblindx(:,:)

      real(8) :: delta_sat
      real(8), allocatable :: cap_param(:,:)
      real(8), allocatable :: cap_fparam(:,:)
      real(8), allocatable :: rlp_param(:,:)
      real(8), allocatable :: rlp_fparam(:,:)
      real(8), allocatable :: vg1(:,:), vg2(:,:), vg3(:,:), vg4(:,:)
      real(8), allocatable :: vg1f(:,:), vg2f(:,:), vg3f(:,:), vg4f(:,:)
      real(8), allocatable :: cp1(:,:), cp2(:,:), cp1f(:,:), cp2f(:,:)
      real(8), allocatable :: vg5(:,:), vg6(:,:), vg7(:,:), vg8(:,:)
      real(8), allocatable :: cp1s(:,:), cp2s(:,:), cp3s(:,:), cp4s(:,:)
      real(8), allocatable :: rlp_table(:,:)
      real(8), allocatable :: sat_out(:)

      logical :: rlpnew = .false.

      parameter(max_rp = 4, max_rpf = 8, max_cp = 6)

      end module comrlp

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM


text:
C     -*- fortran -*-
C     This file is autogenerated with f2py (version:2)
C     It contains Fortran 77 wrappers to fortran functions.

      subroutine f2pywrapsplint (splintf2pywrap, t, n, c, k, a, b,
     & wrk)
      external splint
      integer n
      integer k
      real*8 a
      real*8 b
      real*8 t(n)
      real*8 c(n)
      real*8 wrk(n)
      real*8 splintf2pywrap, splint
      splintf2pywrap = splint(t, n, c, k, a, b, wrk)
      end

      subroutine f2pywrapdblint (dblintf2pywrap, tx, nx, ty, ny, c
     &, kx, ky, xb, xe, yb, ye, wrk)
      external dblint
      integer nx
      integer ny
      integer kx
      integer ky
      real*8 xb
      real*8 xe
      real*8 yb
      real*8 ye
      real*8 tx(nx)
      real*8 ty(ny)
      real*8 c((nx-kx-1)*(ny-ky-1))
      real*8 wrk(nx+ny-kx-ky-2)
      real*8 dblintf2pywrap, dblint
      dblintf2pywrap = dblint(tx, nx, ty, ny, c, kx, ky, xb, xe, y
     &b, ye, wrk)
      end