! Ensure that a non-scalar dummy arguments which are implicitly used inside
! offloaded regions are properly mapped using present_or_copy.

! { dg-do run }

program main
  implicit none

  integer,parameter :: size = 100
  integer :: array(size), i, n

  n = size

  !$acc data copy(array)
  call kernels(array, n)

  !$acc update host(array)

  do i = 1, n
     if (array(i) .ne. i) call abort
  end do

  call parallel(array, n)
  !$acc end data

  do i = 1, n
     if (array(i) .ne. i+i) call abort
  end do
end program main

subroutine kernels (array, n)
  integer, dimension (n) :: array
  integer :: n, i

  !$acc kernels
  do i = 1, n
     array(i) = i
  end do
  !$acc end kernels
end subroutine kernels

subroutine parallel (array, n)
  integer, dimension (n) :: array
  integer :: n, i

  !$acc parallel
  do i = 1, n
     array(i) = i+i
  end do
  !$acc end parallel
end subroutine parallel


!*==dscal.f90  processed by SPAG 7.51RB at 20:37 on  3 Mar 2022
!> \brief \b DSCAL
!
!  =========== DOCUMENTATION ===========
!
! Online html documentation available at
!            http://www.netlib.org/lapack/explore-html/
!
!  Definition:
!  ===========
!
!       SUBROUTINE DSCAL(N,DA,DX,INCX)
!
!       .. Scalar Arguments ..
!       DOUBLE PRECISION DA
!       INTEGER INCX,N
!       ..
!       .. Array Arguments ..
!       DOUBLE PRECISION DX(*)
!       ..
!
!
!> \par Purpose:
!  =============
!>
!> \verbatim
!>
!>    DSCAL scales a vector by a constant.
!>    uses unrolled loops for increment equal to 1.
!> \endverbatim
!
!  Arguments:
!  ==========
!
!> \param[in] N
!> \verbatim
!>          N is INTEGER
!>         number of elements in input vector(s)
!> \endverbatim
!>
!> \param[in] DA
!> \verbatim
!>          DA is DOUBLE PRECISION
!>           On entry, DA specifies the scalar alpha.
!> \endverbatim
!>
!> \param[in,out] DX
!> \verbatim
!>          DX is DOUBLE PRECISION array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
!> \endverbatim
!>
!> \param[in] INCX
!> \verbatim
!>          INCX is INTEGER
!>         storage spacing between elements of DX
!> \endverbatim
!
!  Authors:
!  ========
!
!> \author Univ. of Tennessee
!> \author Univ. of California Berkeley
!> \author Univ. of Colorado Denver
!> \author NAG Ltd.
!
!> \date November 2017
!
!> \ingroup double_blas_level1
!
!> \par Further Details:
!  =====================
!>
!> \verbatim
!>
!>     jack dongarra, linpack, 3/11/78.
!>     modified 3/93 to return if incx .le. 0.
!>     modified 12/3/93, array(1) declarations changed to array(*)
!> \endverbatim
!>
!  =====================================================================
      SUBROUTINE DSCAL(N,Da,Dx,Incx)
      IMPLICIT NONE
!*--DSCAL83
!
!  -- Reference BLAS level1 routine (version 3.8.0) --
!  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
!  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
!     November 2017
!
!     .. Scalar Arguments ..
      DOUBLE PRECISION Da
      INTEGER Incx , N
!     ..
!     .. Array Arguments ..
      DOUBLE PRECISION Dx(*)
!     ..
!
!  =====================================================================
!
!     .. Local Scalars ..
      INTEGER i , m , mp1 , nincx
!     ..
!     .. Intrinsic Functions ..
      INTRINSIC MOD
!     ..
      IF ( N<=0 .OR. Incx<=0 ) RETURN
      IF ( Incx==1 ) THEN
!
!        code for increment equal to 1
!
!
!        clean-up loop
!
         m = MOD(N,5)
         IF ( m/=0 ) THEN
            DO i = 1 , m
               Dx(i) = Da*Dx(i)
            ENDDO
            IF ( N<5 ) RETURN
         ENDIF
         mp1 = m + 1
         DO i = mp1 , N , 5
            Dx(i) = Da*Dx(i)
            Dx(i+1) = Da*Dx(i+1)
            Dx(i+2) = Da*Dx(i+2)
            Dx(i+3) = Da*Dx(i+3)
            Dx(i+4) = Da*Dx(i+4)
         ENDDO
      ELSE
!
!        code for increment not equal to 1
!
         nincx = N*Incx
         DO i = 1 , nincx , Incx
            Dx(i) = Da*Dx(i)
         ENDDO
      ENDIF
      END SUBROUTINE DSCAL


text:
	SUBROUTINE TI_DAYM  ( iyear, imon, iday, iret )
C************************************************************************
C* TI_DAYM								*
C*									*
C* This subroutine returns the number of days in the given month.	*
C* The year must be a full four-digit year.				*
C*									*
C* TI_DAYM  ( IYEAR, IMON, IDAY, IRET )					*
C*									*
C* Input parameters:							*
C*	IYEAR 		INTEGER		Year (YYYY)			*
C*	IMON		INTEGER		Month (MM)			*
C*									*
C* Output parameters:							*
C*	IDAY		INTEGER		Number of days in the month	*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					 -1 = invalid month		*
C**									*
C* Log:									*
C* S. Jacobs/NCEP	 4/96						*
C* I. Durham/GSC	 9/98	Added comment line			*
C************************************************************************
	INTEGER 	month (12)
	LOGICAL		LEAP
C*
	LEAP (iyr) =	( ( MOD ( iyr, 4 ) .eq. 0 ) .and.
     +			  ( ( MOD ( iyr, 100 ) .ne. 0 ) .or.
     +			    ( MOD ( iyr, 400 ) .eq. 0 ) ) )
C*
	DATA		month / 31, 28, 31, 30, 31, 30, 31, 31, 30, 31,
     +				30, 31 /
C------------------------------------------------------------------------
	IF ( imon .gt. 0 .and. imon .lt. 13 ) THEN
C
C*	Pick the number of days for the given month.
C
	    iday = month ( imon )
	    IF  ( ( imon .eq. 2 ) .and. LEAP ( iyear ) )  iday = iday + 1
	    iret = 0
	ELSE
	    print *, "TI_DAYM:  WARNING:  Invalid month = ", imon
	    iret = -1
	END IF
C*
	RETURN
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: sgdecker/pygempak


text:
function leapy(yr)
  implicit none
  logical           :: leapy
  real, intent(in)  :: yr

  if (mod(yr, 4.0) .ne. 0) then
     leapy = .false.
  elseif(mod(yr, 25.0) .ne. 0) then
     leapy = .true.
  elseif(mod(yr, 16.0) .ne. 0) then
     leapy = .false.
  else
     leapy = .true.
  endif

end function leapy

max_stars_repo_name: asnowmang/Calendar

max_issues_repo_name: asnowmang/Calendar

max_forks_repo_name: asnowmang/Calendar


! sdl2_ttf.f90
!
! Fortran 2008 interfaces to SDL2_ttf.
!
! Author:  Philipp Engel
! GitHub:  https://github.com/interkosmos/fortran-sdl2/
! Licence: ISC
module sdl2_ttf
    use, intrinsic :: iso_c_binding
    implicit none
    private

    public :: ttf_init
    public :: ttf_open_font
    public :: ttf_render_text_shaded
    public :: ttf_render_text_solid
    public :: ttf_close_font
    public :: ttf_quit

    interface
        ! int TTF_Init(void)
        function ttf_init() bind(c, name='TTF_Init')
            import :: c_int
            integer(kind=c_int) :: ttf_init
        end function ttf_init

        ! TTF_Font *TTF_OpenFont(const char *file, int ptsize)
        function ttf_open_font(file, ptsize) bind(c, name='TTF_OpenFont')
            import :: c_char, c_int, c_ptr
            character(kind=c_char), intent(in)        :: file
            integer(kind=c_int),    intent(in), value :: ptsize
            type(c_ptr)                               :: ttf_open_font
        end function ttf_open_font

        ! SDL_Surface *TTF_RenderText_Shaded(TTF_Font *font, const char *text, SDL_Color fg, SDL_Color bg)
        function ttf_render_text_shaded_(font, text, fg, bg) bind(c, name='TTF_RenderText_Shaded')
            use :: sdl2, only: sdl_color
            import :: c_char, c_ptr
            type(c_ptr),            intent(in), value :: font
            character(kind=c_char), intent(in)        :: text
            type(sdl_color),        intent(in), value :: fg
            type(sdl_color),        intent(in), value :: bg
            type(c_ptr)                               :: ttf_render_text_shaded_
        end function ttf_render_text_shaded_

        ! SDL_Surface *TTF_RenderText_Solid(TTF_Font *font, const char *text, SDL_Color fg)
        function ttf_render_text_solid_(font, text, fg) bind(c, name='TTF_RenderText_Solid')
            use :: sdl2, only: sdl_color
            import :: c_char, c_ptr
            type(c_ptr),            intent(in), value :: font
            character(kind=c_char), intent(in)        :: text
            type(sdl_color),        intent(in), value :: fg
            type(c_ptr)                               :: ttf_render_text_solid_
        end function ttf_render_text_solid_

        ! void TTF_CloseFont(TTF_Font *font)
        subroutine ttf_close_font(font) bind(c, name='TTF_CloseFont')
            import :: c_ptr
            type(c_ptr), intent(in), value :: font
        end subroutine ttf_close_font

        ! void TTF_Quit(void)
        subroutine ttf_quit() bind(c, name='TTF_Quit')
        end subroutine ttf_quit
    end interface
contains
    ! SDL_Surface *TTF_RenderText_Shaded(TTF_Font *font, const char *text, SDL_Color fg, SDL_Color bg)
    function ttf_render_text_shaded(font, text, fg, bg)
        !! Calls `ttf_render_text_shaded_()` and converts the returned
        !! C pointer to derived type `sdl_surface`.
        use :: sdl2, only: sdl_color, sdl_surface
        type(c_ptr),            intent(in) :: font
        character(kind=c_char), intent(in) :: text
        type(sdl_color),        intent(in) :: fg
        type(sdl_color),        intent(in) :: bg
        type(sdl_surface),      pointer    :: ttf_render_text_shaded
        type(c_ptr)                        :: ptr

        ptr = ttf_render_text_shaded_(font, text, fg, bg)
        call c_f_pointer(ptr, ttf_render_text_shaded)
    end function ttf_render_text_shaded

    ! SDL_Surface *TTF_RenderText_Solid(TTF_Font *font, const char *text, SDL_Color fg)
    function ttf_render_text_solid(font, text, fg)
        !! Calls `ttf_render_text_solid_()` and converts the returned
        !! C pointer to derived type `sdl_surface`.
        use :: sdl2, only: sdl_color, sdl_surface
        type(c_ptr),            intent(in) :: font
        character(kind=c_char), intent(in) :: text
        type(sdl_color),        intent(in) :: fg
        type(sdl_surface),      pointer    :: ttf_render_text_solid
        type(c_ptr)                        :: ptr

        ptr = ttf_render_text_solid_(font, text, fg)
        call c_f_pointer(ptr, ttf_render_text_solid)
    end function ttf_render_text_solid
end module sdl2_ttf


text:
c----67--1---------2---------3---------4---------5---------6---------712
c
c    File: plate1.f   Revised April 23, 1998
c    plate1.f is the primary subroutine for the FEACode mesh generation
c
c    "write(*,...)" changed to "write(iws,...)", 7/29/97
c
c----67--1---------2---------3---------4---------5---------6---------712
c
ccc
c  Convert the main program "plate3d" to the main subroutine called
c  from the shell/DLL accessible subroutine.
c  (GVT 7/29/97, note: look for the "ccc" triple comment for new
c                      additions )
CJF      program plate3d
      subroutine plate1(errorFlag,filename,jobname)
c
c  filename = input file name ("file.in") for the input data file
c             giving the geometric parameters for the mesh generation,
c             use a long character string to allow for a lengthy file
c             name, filename should have no leading blanks
c
c  errorFlag = integer error flag to indicate if the program has had
c              an error and the user should check the *.sta status
c              file in the same directory as the input file
c              (4 byte integer)
c           = 0 = no error
c           = 1 = error, check the *.sta status file
c
c
c----67--1---------2---------3---------4---------5---------6---------712
c
c  This program generates a 3-D FE-model of a plate containing a
c  surface crack of elliptical shape. The crack front is embeddied
c  in a tube following the ellips. A focused mesh is used around
c  the crack front within the tube. In-put deck are generated for
c  the FE-programs WARP3D and ABAQUS.
c
c          VERSION July 11 1997   / Jonas Faleskog /
c
c          Modified  Nov. 22 1999 /JF/
c
c----67--1---------2---------3---------4---------5---------6---------712
c
      implicit none
c
c Array bounds   Zone S: nods(i,j,k) nsm = ns1*ns2*ns3
c                Zone A: noda(i,j,k) nam = na1*na2*na3
c                Zone B: nodb(i,j,k) nbm = nb1*nb2*nb3
c
      include 'plate_common_eln.f'
      include 'plate_common_nod.f'
ccc  Declare the input arguments
      integer errorFlag
      character filename*200,jobname*200
ccc
      integer    nsm,nam,nbm, nesm,nesm2
       parameter (nsm=60000,nam=900000,nbm=1700000)
ccc      parameter (nsm=20000,nam=300000,nbm=600000)
ccc      parameter (nsm=20000,nam=140000,nbm=140000)
ccc      parameter (nsm=20000,nam=80000,nbm=80000)
      parameter (nesm=200, nesm2=nesm*nesm)
c
      integer  iws, ns1,ns2,ns3, na1,na2,na3, nb1,nb2,nb3
c
      integer  nea2,neb2, estk_s(2,nesm), estk_a(2,nesm2),
     &         estk_b(2,nesm2),nstk_s(3,401),nstk_a(3,401),
     &         nstk_b(3,401), tcr,tge,tbc,keyhole,analysis,
     &         numLoadStep,eqsolver,ncq, i
c
      integer  nods(nsm),noda(nam),nodb(nbm),
     1         etyp,elnum,no_of_nodes
c
      integer  eset_size
      parameter (eset_size = 20000)
      integer  efront(2,eset_size),ecrsur(2,eset_size),
     &         ecrtip(2,eset_size),easym(2,eset_size),
     &         ecsym(2,eset_size),erem(2,eset_size),
     &         esid(2,eset_size),etop(2,eset_size)
c
      double precision dmat(2000), rgeo(20),rload(10),
     1                 yg(0:200),xag(0:200),xbg(0:200),zg(2,0:200),
     2                 rfm(50),bias(10),crfnvec(3),cqvec(3,200)
c Make the status filename "stfile" the same size as the input filename
      character    prog*20,NeutralFormat*20,stfile*200
      logical      demoFlag, statusFlag
ccc
c
      integer      mf,mr,mv,mh,m1,m2,ma,na,mb,nb,lt,lred
      common /nod/ mf,mr,mv,mh,m1,m2,ma,na,mb,nb,lt,lred
      integer      ims,jms,kms,ima,jma,kma,imb,jmb
      common /max/ ims,jms,kms,ima,jma,kma,imb,jmb
      integer         kstart,kstep
      common /nblock/ kstart,kstep
      integer         ksr1,kar1,kar2,rtype,sfred,sfred_type,sjred_type
      common /reduce/ ksr1,kar1,kar2,rtype,sfred,sfred_type,sjred_type
c
      double precision t,w,c,a,alfa
      common /geom/    t,w,c,a,alfa
ccc
      errorFlag = 0
c Check for leading blank spaces in the filename string and remove
cgt      filename = ADJUSTL(filename)
ccc
c
c
c..Input is read from file mesh3d_scp.in
c
      call read_input( prog,etyp,analysis,numLoadStep,eqsolver,
     &          tcr,tge,tbc,rload, dmat,keyhole,
     &          nods,ns1,ns2,ns3,nsm, noda,na1,na2,na3,nam,
     &          nodb,nb1,nb2,nb3,nbm, nesm,nesm2,nea2,neb2,
     &          estk_s,estk_a,estk_b, rgeo,yg,xag,xbg,zg,rfm,bias,
     &          NeutralFormat, demoFlag,statusFlag, iws,errorFlag,
     &          filename,jobname,stfile)
      if ( errorFlag .gt. 0 ) goto 999
c
c TEST
c
c Note!  "write(iws,...)" to write to the status *.sta file
c
      if (demoFlag) write(iws,'(/t1,a/)') '>>> Demo Version <<<'
      write(iws,'(t5,a,i2,a,i4,a,i2)')  'TCR =',tcr,
     &                    '   TGE =',tge,'   TBC =',tbc
      write(iws,'(t5,3(a,g12.4))') 'W0 =',w,'  L0 =',yg(lt),'  t0 =',t
      write(iws,'(t5,2(a,g12.4)/)') 'Rc =',rgeo(4),'  Rp =',rgeo(5)
      write(iws,*) ' dmat:'
      do i=1, (int(dmat(1)+0.5) + 1)
         write(iws,*) real(dmat(i))
      enddo
c
c      write(iws,*) ' After read_input'
c      write(iws,*) ' indicies'
c      write(iws,'(t1,a,8i4)') ' mf,mr: ',mf,mr
c      write(iws,'(t1,a,8i4)') ' mv,mh,m1,m2,ma,na: ',mv,mh,m1,m2,ma,na
c      write(iws,'(t1,a,8i4)') ' mb,nb: ',mb,nb
c      write(iws,'(t1,a,8i4)') ' lt,lred: ',lt,lred
c      write(iws,'(t1,a,8i4)') ' ims,jms,kms: ',ims,jms,kms
c      write(iws,'(t1,a,8i4)') ' ima,jma,kma: ',ima,jma,kma
c      write(iws,'(t1,a,8i4)') ' imb,jmb: ',imb,jmb
c      write(iws,'(t1,a,8i4)') ' ksr1,kar1,kar2: ',ksr1,kar1,kar2
c      write(iws,'(t1,a,8i4)') ' rtype,sfred,sfred_type,sjred_type: ',
c     1                        rtype,sfred,sfred_type,sjred_type
c   write(iws,'(t1,a,5g12.4)') 't,w:',t,w
c   write(iws,'(t1,a,5g12.4)') 'c,a,alfa:',c,a,alfa
c      write(iws,'(t3,a)') 'BIAS'
c      do i=1, 8
c         write(iws,'(g15.6)') bias(i)
c      enddo
c      write(iws,'(t3,a)') 'rfm'
c      do i=1, mr
c         write(iws,'(g15.6)') rfm(i)
c      enddo
c      write(iws,'(t3,a)') 'yg, zg'
c      do i=0, lt
c         write(iws,'(3g15.6)') yg(i),zg(1,i), zg(2,i)
c      enddo
c      write(iws,'(t3,a,i2)') 'xag;  m2 =',m2
c      do i=0, m2
c         if (i.le.1) then
c            write(iws,'(g15.6)') xag(i)
c         else
c            write(iws,'(2g15.6)') xag(i),
c     &                          (xag(i)-xag(i-1))/(xag(i-1)-xag(i-2))
c         endif
c      enddo
c      write(iws,'(t3,a,i2)') 'xbg;  mb =',mb
c      do i=0, mb
c         if (i.le.1) then
c            write(iws,'(g15.6)') xbg(i)
c         else
c            write(iws,'(2g15.6)') xbg(i),
c     &                          (xbg(i)-xbg(i-1))/(xbg(i-1)-xbg(i-2))
c         endif
c      enddo
c
c Generate node numbers, store in the arrays NODS,NODA & NODB
c
      call nodnumber(nods,ns1,ns2,ns3,noda,na1,na2,na3,
     &               nodb,nb1,nb2,nb3, iws,errorFlag)
      if ( errorFlag .gt. 0 ) then
         write(iws,*)'Stopped after error in nodnumber subroutine'
         goto 999
      endif
c
c Generate node coordinates
c
      if (statusFlag) then
c         close(unit=iws)
c         open(unit=iws, file=stfile, position='append')
      endif
      write(iws,'(/t1,a)') '* generating coordinates'
      write(*,'(/t1,a)') '* generating coordinates'
      call node_coordinates(iws,nods,ns1,ns2,ns3, noda,na1,na2,na3,
     &                      nodb,nb1,nb2,nb3, yg,xag,xbg,rfm, tcr,
     &                      keyhole,errorFlag)
      if ( errorFlag .gt. 0 ) then
        write(iws,*)'Stopped after error in node_coordinates subroutine'
         goto 999
      endif
c
c..Generate elements
c
      if (statusFlag) then
         close(unit=iws)
         open(unit=iws, file=stfile, position='append')
      endif
      write(iws,'(/t1,a)') '* generating elements'
      write(*,'(/t1,a)') '* generating elements'
      call make_element(nods,ns1,ns2,ns3, noda,na1,na2,na3,
     1                  nodb,nb1,nb2,nb3, elnum, analysis,keyhole,
     2                  nea2,neb2, estk_s,estk_a,estk_b, prog,
     3                  iws,errorFlag)
      if ( errorFlag .gt. 0 ) then
         write(iws,*)'Stopped after error in make_element subroutine'
         goto 999
      end if
      write(iws,'(t15,a,i5,a)') '=> # of elements = ',elnum
      write(*,'(t15,a,i5,a)') '=> # of elements = ',elnum
c
c..Sort out the nodes in the model
c
      if (statusFlag) then
         close(unit=iws)
         open(unit=iws, file=stfile, position='append')
      endif
      write(iws,'(/t1,a)') '* sort out the unique nodes'
      write(*,'(/t1,a)') '* sort out the unique nodes'
      call sort_out_nodes(elnum,no_of_nodes,etyp,prog,
     1                    nstk_s,nstk_a,nstk_b, nods,ns1,ns2,ns3,
     2                    noda,na1,na2,na3, nodb,nb1,nb2,nb3, iws)
c
c Coordinate transformation - type defined by parameter TGE
c
      if (statusFlag) then
         close(unit=iws)
         open(unit=iws, file=stfile, position='append')
      endif
      write(iws,'(/t1,a)') '* Coordinate Transformation'
      write(*,'(/t1,a)') '* Coordinate Transformation'
      call coord_transf(tge,elnum,etyp,rgeo,crfnvec,cqvec,ncq,
     &                  ns1,ns2,ns3,nods)
c
c Generate element sets containing element on all the surfacres
c of the model.
c
      if (statusFlag) then
         close(unit=iws)
         open(unit=iws, file=stfile, position='append')
      endif
      write(iws,'(/t1,a)') '* Retrieve element surfaces'
      write(*,'(/t1,a)')'* Retrieve element surfaces - may take a while'
      call gen_surface(efront,ecrsur,ecrtip,easym,ecsym,erem,esid,etop,
     1                 nods,ns1,ns2,ns3, noda,na1,na2,na3,
     2                 nodb,nb1,nb2,nb3, elnum,iws,errorFlag)
c
c Generate plot files containing wire frames of key surfaces
c
ccc      if (statusFlag) then
ccc         close(unit=iws)
ccc         open(unit=iws, file=stfile, position='append')
ccc      endif
ccc      write(iws,'(/t1,a)') '* Generate Graphics files'
ccc      write(*,'(/t1,a)') '* Generate Graphics files'
ccc      call gen_graph(efront,ecrsur,easym,ecsym,erem,esid,etop,
ccc     &              tge,no_of_nodes,elnum,iws,jobname)
c
c Generate FEA in-put deck and an info-file for Post-Processing
c
      if ( demoFlag ) goto 999
      if (statusFlag) then
         close(unit=iws)
         open(unit=iws, file=stfile, position='append')
      endif
      if (NeutralFormat.eq.'PATRAN') then
         call patran_neu(etyp,elnum,no_of_nodes,dmat,rgeo,rload,
     &                   analysis,tcr,tge,tbc,ecrsur,ecrtip,ecsym,etop,
     &                   erem,esid, nods,ns1,ns2,ns3, noda,na1,na2,na3,
     &                   nodb,nb1,nb2,nb3, jobname,iws,errorFlag)
      elseif (NeutralFormat.eq.'FEMAP') then
         call femap_neu(etyp,elnum,dmat,jobname)
      elseif (prog.eq.'ABAQUS') then
         write(iws,'(/t1,a,a)') '* Generate in-put deck for ',prog
         write(*,'(/t1,a)') '* Generate ABAQUS in-put deck'
         call abaqus_inp(etyp,elnum,analysis,numLoadStep,keyhole,
     &                   tge,tcr,tbc,rgeo,rload,
     &                   dmat,ecrsur,ecsym,erem,esid,etop,
     &                   crfnvec,cqvec,ncq,
     &                   nods,ns1,ns2,ns3, noda,na1,na2,na3,
     &                   nodb,nb1,nb2,nb3, jobname,iws,errorFlag)
         call abaqus_post(tge,tcr,tbc,dmat,nods,ns1,ns2,ns3,jobname)
c
      elseif (prog .eq. 'WARP3D') then
         write(*,'(/t1,a)') '* Generate WARP3D in-put deck'
         call warp3d_inp(etyp,elnum,no_of_nodes,dmat,rgeo,rload,
     &                   analysis,keyhole,numLoadStep,
     &                   tcr,tge,tbc,eqsolver,crfnvec,
     &                   ecrsur,ecrtip,ecsym,etop, erem,esid,
     &                   nods,ns1,ns2,ns3, noda,na1,na2,na3,
     &                   nodb,nb1,nb2,nb3,
     &                   jobname,iws,errorFlag)
         call warp3d_post(tge,tcr,tbc,dmat,nods,ns1,ns2,ns3,jobname)
c
      endif
c
      if (statusFlag) then
         close(unit=iws)
         open(unit=iws, file=stfile, position='append')
      endif
c
      write(iws,'(/t1,a)') '*  Done with mesh generation'
c
999   continue
      close(iws)
c
      return
      end
c
c----67--1---------2---------3---------4---------5---------6---------712
c
      subroutine read_input( prog,etyp,analysis,numLoadStep,eqsolver,
     &           tcr,tge,tbc,rload,dmat,keyhole,
     &           nods,ns1,ns2,ns3,nsm, noda,na1,na2,na3,nam,
     &           nodb,nb1,nb2,nb3,nbm, nesm,nesm2,nea2,neb2,
     &           estk_s,estk_a,estk_b, rgeo,yg,xag,xbg,zg,rfm,bias,
     &           NeutralFormat, demoFlag,statusFlag, iws,errorFlag,
     &           filename,jobname,stfile)
c
c The subroutine reads in-put parameters from the file "plate3d.in"
c mesh3d_nsct.in. Some of the parameters are checked so that
c they are located in a reasonable interval.
c
      implicit none
      include 'plate_common_nod.f'
c
      integer ns1,ns2,ns3,nsm,nods(nsm), na1,na2,na3,nam,noda(nam),
     1        nb1,nb2,nb3,nbm,nodb(nbm)
      integer analysis,numLoadStep,eqsolver,tcr,tge,tbc,keyhole
      integer iws,nesm,nesm2,nes2,nea2,neb2,estk_s(2,nesm),
     1        estk_a(2,nesm2),estk_b(2,nesm2)
      double precision  rload(*),dmat(*),rgeo(*),yg(0:200),xag(0:200),
     1                  xbg(0:200),zg(2,0:200),rfm(50),bias(10),b
      integer errorFlag
      character filename*200,jobname*200,stfile*200
c
      integer i,j,k,n,num,etyp,max_node_no,io,ivek(20),
     1        m1min,m2min,m12m,m12e,m12,mb1,nmin
      double precision  length,rvek(15),dx0,xmax,dxmax,
     3                  x(200),xi,beta,beta0,betam,c1,c2,psi,tmp,
     4                  lambda,lambda0,lambdam,lambda1
      double precision  bias_crack
ccc  Make the inpfile string the same size as filename
      character prog*20,NeutralFormat*20,row*200,inpfile*200
ccc  Old string sizes
ccc      character prog*20,row*200,inpfile*10
      logical   ok, demoFlag,statusFlag
      double precision  pi,one
      parameter (one=1.0)
c
      integer      mf,mr,mv,mh,m1,m2,ma,na,mb,nb,lt,lred
      common /nod/ mf,mr,mv,mh,m1,m2,ma,na,mb,nb,lt,lred
      integer      ims,jms,kms,ima,jma,kma,imb,jmb
      common /max/ ims,jms,kms,ima,jma,kma,imb,jmb
      integer         ksr1,kar1,kar2,rtype,sfred,sfred_type,sjred_type
      common /reduce/ ksr1,kar1,kar2,rtype,sfred,sfred_type,sjred_type
      double precision t,w,c,a,alfa
      common    /geom/ t,w,c,a,alfa
c
      pi = 4.*atan(one)
c
ccc
c  Use the filename given by the calling routine
c      inpfile = 'plate3d.in'
      inpfile = filename
      inquire( file=inpfile, exist=ok )
      iws = 44
c      open( unit=iws, file='tmp.sta',  status='scratch' )
      open( unit=iws, status='scratch' )
      if (.not.ok) then
         write(iws,'(t1,3a)')'>> The File: ',inpfile(1:20),
     &                       ' does not exist!'
         errorFlag = 1
         return
      endif
      write(iws,201) '* reading in-put-data from infile:',inpfile
      close(iws)
201   format(/t1,a,/,a/)
c
      io = 12
      open(unit=io,file=inpfile,status='old')
c
c DemoFlag & statusFlag:
c
c     If DemoFlag set this program version is only for demonstration
c     purposes. Graphics files will still be generated for all
c     geometries but FE in-put deck will only be generated for the
c     basic PLATE model.
c
c     If statusFlag is set a "job.sta" will be generated otherwise not.
c     The status file is opened on UNIT=10
c
      read(io,'(a)') row
      call char_to_int(row,ivek,num)
cc      if (num.ne.2) goto 99
c Remove Demoflag-option as in input option /JF Nov. 22 1999/
c
      if (num.ne.1) goto 99
      ivek(2) = ivek(1)
      ivek(1) = 0
      if (ivek(1).eq.0) then
         demoFlag = .false.
      else
         demoFlag = .true.
      endif
c
c Open a mesh-generating status file on UNIT=10
c
ccc  Append the given filename with the "*.sta" extension
c      i = index(jobname,' ',.false.)
c      stfile = jobname(1:i)//'.sta'
      call appendFile(jobname,'.sta',stfile)
      iws = 10
      if (ivek(2).eq.0) then
         statusFlag = .false.
         open(unit=iws, status='scratch')
c         open(unit=iws, file=stfile, status='scratch')
      else
         statusFlag = .true.
         open(unit=iws, file=stfile, status='replace')
      endif
      call writeDate(iws)
c
c FE-code: prog (read and left adjust)
c
      read(io,'(a)') row
      i = 0
 10   continue
         i = i + 1
      if (row(i:i).eq.' ') goto 10
         j = i
 15   continue
         j = j + 1
      if (row(j:j).ne.' ') goto 15
      j = j - 1
      do k=1, 20
         prog(k:k) = ' '
      enddo
      prog(1:(j-i+1)) = row(i:j)
c
c Make sure prog is all upper case
c
      call ToUpper( prog )
c
c If in-put should be generated on neutral format use
c WARP3D format as default until the specific neutral format
c is generated.
c
      if     ( prog .eq. 'PATRAN' ) then
         NeutralFormat = prog
         prog = 'WARP3D    '
      elseif ( prog .eq. 'FEMAP' ) then
         NeutralFormat = prog
         prog = 'WARP3D    '
      else
         NeutralFormat = '          '
      endif
c
c Type of analysis:
c
c      analysis = 1   static,  linear elastic
c      analysis = 2   static,  nonlinear elastic
c      analysis = 3   static,  incremental elastic-plastic
c      analysis = 4 static,  incremental finite strain el.-pl.
c      analysis = 5   dynamic, linear elastic
c      analysis = 6   dynamic, nonlinear elastic
c      analysis = 7   dynamic, incremental elastic-plastic
c      analysis = 8   ...
c      analysis = 9   ...
c
c   Equation Solver   Direct Sparse Solver      eqsolver=1
c                     Conjugate Gradient Solver eqsolver=2
c
      read(io,'(a)') row
      call char_to_int(row,ivek,num)
      if ((num.ne.2).and.(num.ne.3)) goto 99
      analysis = ivek(1)
      numLoadStep = ivek(2)
      if ( (num.eq.3) .and. (ivek(3).eq.1)) then
         eqsolver = 1
      else
         eqsolver = 2
      endif
c
c Material parameters; elastic: E, ny
c
c   Material parameters are stored in array dmat(*)
c
c      analysis = 1   static,  linear elastic
c        dmat(1) = no of material constants defined = 2
c        dmat(2) = Young's modulii
c        dmat(3) = Poisson's ratio
c
c      analysis = 2  static, nonlinear elastic (deformation Plasticity)
c        dmat(1) = no of material constants defined = 4 (5 if abaqus)
c        dmat(2) = Young's modulii
c        dmat(3) = Poisson's ratio
c        dmat(4) = sigY
c        dmat(5) = n = hardening
c        dmat(6) = alpha  = offset strain (Ramberg-Osgood model)
c
c      analysis = 3   static,  incremental elastic-plastic
c        dmat(1) = no of material constants defined >= 4
c        dmat(2) = Young's modulii
c        dmat(3) = Poisson's ratio
c        dmat(4) = sig1
c        dmat(5) = eps1
c        dmat(6) = sig2
c        dmat(7) = eps3
c        dmat(6) = ..
c        dmat(7) = ..
c
      read(io,'(a)') row
      call char_to_real(row,rvek,num)
      if (num.ne.1) goto 99
      dmat(1) = rvek(1)
      n = int( dmat(1) + 0.5 )
      j = 0
c
c Material data constants must be written in the in-put file
c using columns 1 to 200.
c
 20   continue
         read(io,'(a)') row
         call char_to_real(row,rvek,num)
         do i=1, num
            j = j + 1
            dmat(j+1) = rvek(i)
         enddo
      if (j.lt.n) then
         goto 20
      elseif (j.gt.n) then
         write(iws,*) 'ERROR: no. of material constants does not match'
         goto 99
      endif
c
c Geometry parameters: L, W, T, Rc, Rp, Phi0
c
      read(io,'(a)') row
      call char_to_real(row,rvek,num)
      call determ_geom(tge,w,length,t,rgeo,rvek,num, iws,errorFlag)
      if ( errorFlag .gt. 0 ) then
         write(iws,*)'Stop after error in determ_geom subroutine'
         return
      end if
c
c Crack geomtry: a, c
c
c  TCR = 1  embedded crack
c      = 2  semi elliptical crack (c - major axis = free surface)
c      = 3  semi elliptical crack (a - minor axis = free surface)
c      = 4  quarter (corner) crack (a,c - both free surfaces)
c
      read(io,'(a)') row
      call char_to_real(row,rvek,num)
      if (num.ne.4) goto 99
      tcr = int(rvek(1) + 0.5)
      a    = rvek(2)
      c    = rvek(3)
      keyhole = int( rvek(4) + 0.5 )
c
c Test overall dimensions of the body in relation to crack dimensions.
c
      if (a .gt. (0.99*t)) then
         write(iws,'(t2,a)') '>> a > 0.99*t => decrease a'
         errorFlag = 1
         return
      endif
      if (w .lt. (c + 0.01*t)) then
         write(iws,'(t2,a)') '>> W < c + 0.01*t => increase W'
         errorFlag = 1
         return
      endif
c
c Load definition
c
c  TBC = 1   prescribed end displacements
c  TBC = 2   prescribed end tractions and   internal pressure
c                           (cylinder, sphere, pipe)
c  TBC = 3   prescribed tractions on the crack face
c
      do i=1, 10
         rload(i) = 0
      enddo
      read(io,'(a)') row
      call char_to_real(row,rvek,num)
      tbc = int(rvek(1) + 0.5)
      do i=2, num
         rload(i-1) = rvek(i)
      enddo
c
c Mesh parameters: elemtype, mf
c
      read(io,'(a)') row
      call char_to_real(row,rvek,num)
      if (num.eq.1) then
         etyp = int(rvek(1)+0.5)
c . . . Default values
         mf = 10
         if (keyhole.ne.0) then
            mr = 7
         else
            mr = 5
         endif
      elseif (num.eq.2) then
         etyp = int(rvek(1)+0.5)
         mr   = int(rvek(2)+0.5)
      elseif (num.eq.3) then
         etyp = int(rvek(1)+0.5)
         mr   = int(rvek(2)+0.5)
         mf   = int(rvek(3)+0.5)
      else
         goto 99
      endif
c
      close(io)
c
c Check input data, so that parameters are within permitted limits
c
c    =>   To be implemented
c
c Determine model parameters
c
c  # of elements
c
c  Zone S
c
      sfred_type = 1
      sjred_type = 1
      na = 12
      nb = 8
c     na = 12*(c/a)**0.15
      na = nb + (na-nb)*(c/a)**0.5
      if (mod(na,2) .ne. 0) na = na+1
c
c if s = 0 is a free surface => increase  na
c
      if ((tcr.eq.3) .or. (tcr.eq.4)) na = na + 2
      sfred = mr - 2
      if (sfred.lt.2) sfred = mr - 1
      call determ_mv_mh(mf,mv,mh,sfred_type)
c
c  Zone A
c
      m1min = 2
      m2min = 3
      m12m  = 4
      psi = a/t
      if (psi .le. 0.5) then
      m12e  = 5
         m1  = (1.-2.*psi)*real(m1min) + 2.*psi*(real(m12m)+0.99)
         m12 = 2*int((1.-2.*psi)*real(m12e) + 2.*psi*(real(m12m)+0.99))
         m2  = m12 - m1
         if (mod((m1+m2),2).ne.0) m2 = m2
      else
      m12e  = 6
         psi = (1.-a/t)
         m2  = (1.-2.*psi)*real(m2min) + 2.*psi*(real(m12m)+0.99)
         m12 = 2*int((1.-2.*psi)*real(m12e) + 2.*psi*(real(m12m)+0.99))
         m1  = m12 - m2
         if (mod((m1+m2),2).ne.0) m1 = m1
      endif
      ma = m1 + m2 + 2*mh
c
c Geometry
c
c  Zone S
c
      alfa = sqrt(c*c-a*a)/2.
      beta0 = 0.35
      betam = 0.50
      psi  = abs(2.*(a/t)-1.)
      beta = (1.0-psi)*betam + psi*beta0
      c1 = 1.0 / ( 1.0 + beta*real(m1min) )
      c2 = 4.0 / ( 2.0 + beta*real(2*m12m) ) - 2.0*c1
      psi = a/t
      if (psi .gt. 0.5 ) psi = 1.0 - psi
      b = t*psi*(c1 + c2*psi)
c
 30   continue
c
c  Determine the bias in the focused mesh
      bias(1) =  bias_crack(b,b*beta,mr)
      xi = 1 - mr
      rfm(1) = beta*b*bias(1)**xi
      do i=2, mr
         xi = i - 1
         rfm(i) = rfm(i-1) + rfm(1)*bias(1)**xi
      enddo
c
      bias(2) = 1.2
      xi = mv
      yg(0) = 0.0
      yg(1) = b * ( bias(2)-1. ) / ( bias(2)**xi-1. )
      do i=2, mv
         xi = i-1
         yg(i) = yg(i-1) + yg(1)*bias(2)**xi
      enddo
c
c  Zone A
c
      bias(4) = 1.5
      dx0 = yg(mv)-yg(mv-1)
      dxmax = 7.5*t / real(ma/2)
      if (dxmax.gt.t) then
         if ( (tge.eq.201).or.(tge.eq.202).or.(tge.eq.301).or.
     &        (tge.eq.302).or.(tge.eq.305).or.(tge.eq.306) ) then
            dxmax = t
         endif
      endif
      xmax = length - yg(mv-1)
      nmin = 1
      call determ_remote_mesh(nmin,dx0,dxmax,xmax,bias(4),x,n)
c      lt = n + mv-1
      lt = n + mv
      do i=1, n
         j = i+mv
         yg(j) = yg(mv) + x(i)
      enddo
c
      do i=0, lt
         zg(1,i) = 0.0
         zg(2,i) = t
      enddo
      if ( length .lt. t) then
         rtype = 0
      elseif (length .lt. (1.6*t)) then
         rtype = 1
         lred = 0
         do i=mv, lt
            if ((yg(i).gt.(0.5*t)).and.(lred.eq.0)) lred = i
         enddo
      else
         rtype = 2
         lred = 0
         do i=mv, lt
            if ((yg(i).gt.(0.6*t)).and.(lred.eq.0)) lred = i
         enddo
      endif
      if (tge .gt. 200) then
         rtype = 2
         lred = 0
         do i=mv, lt
            if ((yg(i).gt.(0.5*t)).and.(lred.eq.0)) lred = i
         enddo
      endif
c
c If rtype = 0 no mesh coarsening
c    rtype = 1 mesh coarsening 2 to 1 in Zone A
c              (one extra layer is included at the red. site)
c    rtype = 2 mesh coarsening 4 to 1 in Zone A and 2 to 1 in Zone B
c              (currently not used)
c
      if (rtype.eq.0) then
         kma=2*lt+1
      elseif (rtype.ge.1) then
         lt=lt+1
         kma=2*lt+1
         do i=lt, lred+1, -1
            yg(i)   = yg(i-1)
            zg(1,i) = zg(1,i-1)
            zg(2,i) = zg(2,i-1)
         enddo
         yg(lred)   = 0.5*( yg(lred-1)   + yg(lred+1) )
         zg(1,lred) = 0.5*( zg(1,lred-1) + zg(1,lred+1) )
         zg(2,lred) = 0.5*( zg(2,lred-1) + zg(2,lred+1) )
      endif
      if (rtype.eq.2) then
         lt=lt+1
         kma=2*lt+1
         do i=lt,lred+3, -1
            yg(i)   = yg(i-1)
            zg(1,i) = zg(1,i-1)
            zg(2,i) = zg(2,i-1)
         enddo
         yg(lred+2)   = 0.5*( yg(lred+1)   + yg(lred+3) )
         zg(1,lred+2) = 0.5*( zg(1,lred+1) + zg(1,lred+3) )
         zg(2,lred+2) = 0.5*( zg(2,lred+1) + zg(2,lred+3) )
      endif
c
c  Mesh parameters in region X > c+b; Zone A and Zone B
c
      psi = (a+b)/t
      lambda0 = 1.8
      lambdam = 1.1
      lambda1 = 1.8
      if (psi .lt. 0.5) then
         lambda = (1.-2.*psi)*lambda0 + 2.*psi*lambdam
      else
         lambda = 2.*(1.-psi)*lambdam + (2.*psi-1.)*lambda1
      endif
      bias(5) = lambda
c
c  define element lengths in X-dir. in zone A & B
c
      if ( ((a+b)/t) .lt. 0.45 ) then
         nmin = m2
      else
         mb1 = 2
         mb  = int( 2.*(1.-psi) + (2.*psi-1.)*real(mb1) + 0.5 )
         nmin = m2 + mb
      endif
c   decrease b if W < tmp
      tmp = c + b + beta*b*real(nmin)*0.8
      if ( w .lt. tmp ) then
ccc         write(*,'(t2,a,g11.4)') 'old b =',b
         write(iws,'(t2,a,g11.4)') 'old b =',b
         b = 0.99*(w-c) / ( 1.0 + beta*real(nmin)*0.8 )
ccc         write(*,'(t2,a,g11.4)') 'new b =',b
         write(iws,'(t2,a,g11.4)') 'new b =',b
         goto 30
      endif
c
c  If TGE > 200  => Relate dxmax to the inner radius of the Cylinder
c
      dx0 = rfm(mr) - rfm(mr-1)
      if ( rtype .lt. 2) then
         dxmax = 7.5*t / real(nb)
      else
         dxmax = 7.5*t / real(nb/2)
      endif
      if (dxmax.gt.t) then
         if ( (tge.eq.203).or.(tge.eq.204).or.(tge.eq.303).or.
     &        (tge.eq.304).or.(tge.eq.307).or.(tge.eq.308) ) then
            dxmax = t
         endif
      endif
      xmax = w - (c + b)
      call determ_remote_mesh(nmin,dx0,dxmax,xmax,bias(5),x,n)
      xag(0) = c + b
      do i=1, n
         xag(i) = xag(0) + x(i)
      enddo
      mb = n - m2
      do i=m2, n
         xbg(i-m2) = xag(i)
      enddo
c
c  Indicies in node number arrays
c
c  Zone S:
c
      ims = 2*mf+1
      jms = 2*na*sjred_type+1
      kms = 2*mr + 1
      if ( (sfred_type.eq.1) .and. (sjred_type.gt.1) ) then
         kms = kms + 2
      elseif ( (sfred_type.gt.1) .and. (sjred_type.eq.1) ) then
         kms = kms + 2
      elseif ( (sfred_type.gt.1) .and. (sjred_type.gt.1) ) then
         kms = kms + 4
      endif
      ksr1 = 2*sfred+1
c
c  Zone A:
      ima = 2*ma+1
      jma = 2*na+1
      kar1 = 2*lred+1
      kar2 = 2*(lred+2)+1
c
c  Zone B:
      if (mb.gt.0) then
         imb=2*mb+1
      else
         imb = 1
      endif
      jmb=2*nb+1
c
c TEST
c
c      write(iws,'(t3,a)') 'Inside read_input'
c      write(iws,'(t3,a)') 'BIAS'
c      do i=1, 8
c         write(iws,'(g15.6)') bias(i)
c      enddo
c      write(iws,'(t3,a)') 'rfm'
c      do i=1, mr
c         write(iws,'(g15.6)') rfm(i)
c      enddo
c      write(iws,'(t3,a)') 'yg, zg'
c      do i=0, lt
c         write(iws,'(3g15.6)') yg(i),zg(1,i), zg(2,i)
c      enddo
c      write(iws,'(t3,a)') 'xag'
c      do i=1, m2
c         write(iws,'(g15.6)') xag(i)
c      enddo
c      write(iws,'(t3,a)') 'xbg'
c      do i=1, mb
c         write(iws,'(g15.6)') xbg(i)
c      enddo
c
c  Set size of matrices
c
      ns1 = ims
      ns2 = jms
      ns3 = kms
      if ((ns1*ns2*ns3).gt.nsm) then
         write(iws,'(t1,a)') '>> increase the size of array nods(i,j,k)'
         write(iws,'(t4,a,i8)') 'current  size = ', nsm
         write(iws,'(t4,a,i8)') 'required size = ', ns1*ns2*ns3
         errorFlag = 1
         return
      endif
      do i=1, nsm
         nods(i) = 0
      enddo
c
      na1 = ima
      na2 = jma
      na3 = kma
      if ((na1*na2*na3).gt.nam) then
         write(iws,'(t1,a)') '>> increase the size of array noda(i,j,k)'
         write(iws,'(t4,a,i8)') 'current  size = ', nam
         write(iws,'(t4,a,i8)') 'required size = ', na1*na2*na3
         errorFlag = 1
         return
      endif
      do i=1, nam
         noda(i) = 0
      enddo
c
      nb1 = imb
      nb2 = jmb
      nb3 = kma
      if ((nb1*nb2*nb3).gt.nbm) then
         write(iws,'(t1,a)') '>> increase the size of array nodb(i,j,k)'
         write(iws,'(t4,a,i8)') 'current  size = ', nbm
         write(iws,'(t4,a,i8)') 'required size = ', nb1*nb2*nb3
         errorFlag = 1
         return
      endif
      do i=1, nbm
         nodb(i) = 0
      enddo
c
c  Check if maximum size of npos (coordinate array) is sufficient
c
      max_node_no = jms*ims*kms + (jma*ima+jmb*imb)*kma
      if (max_node_no.gt.inm) then
        write(iws,'(t1,a)') '>> increase the size of array npos(inm,3)'
        write(iws,'(t4,a,i8)') 'current  size = ', inm
        write(iws,'(t4,a,i8)') 'required size = ', max_node_no
        errorFlag = 1
        return
      endif
c
      nes2 = (kms+1)/2
      if (nes2.gt.nesm) then
        write(iws,'(t1,a)') '>> increase the array size parameter nesm'
        write(iws,'(t4,a,i8)') 'current  size = ', nesm
        write(iws,'(t4,a,i8)') 'required size = ', nes2
        errorFlag = 1
        return
      endif
      nea2 = (kma+1)/2
      if ( (nea2*(jma+1)/2) .gt. (nesm*nesm) ) then
        write(iws,'(t1,a)') '>> increase the array size parameter nesm'
        write(iws,'(t4,a,i8)') 'current  size = ', nesm
        tmp = real(nea2*(jma+1)/2)
        write(iws,'(t4,a,i8)') 'required size = ', int( sqrt(tmp) )+1
        errorFlag = 1
        return
      endif
      neb2 = (kma+1)/2
      if ( (neb2*(imb+1)/2) .gt. (nesm*nesm) ) then
        write(iws,'(t1,a)') '>> increase the array size parameter nesm'
        write(iws,'(t4,a,i8)') 'current  size = ', nesm
        tmp = real(neb2*(imb+1)/2)
        write(iws,'(t4,a,i8)') 'required size = ', int( sqrt(tmp) )+1
        errorFlag = 1
        return
      endif
      if ( (neb2) .gt. (nesm) ) then
        write(iws,'(t1,a)') '>> increase the array size parameter nesm'
        write(iws,'(t4,a,i8)') 'current  size = ', nesm
        write(iws,'(t4,a,i8)') 'required size = ', neb2
        errorFlag = 1
        return
      endif
      do i=1, nesm
         do j=1,2
            estk_s(j,i) = 0
         enddo
      enddo
      do i=1, nesm2
         do j=1,2
            estk_a(j,i) = 0
            estk_b(j,i) = 0
         enddo
      enddo
c
      return
c
99    write(iws,'(t10,a/t10,a)')
     1  '>> the number of data given does not match the number of',
     2  '   data required at each line, => correct the in-put file!!'
      close(io)
      return
      end
c
c----67--1---------2---------3---------4---------5---------6---------712
c
      subroutine determ_geom(tge,w0,l0,t0,rgeo,rvek,num, iws,errorFlag)
c
c Determine the basic plate geometric parameters (W,L,t), given TGE
c
      implicit none
      integer          tge,num,iws,errorFlag
      double precision w0,l0,t0,rgeo(*),phi0,rvek(*), l,rc,rp,pi,one
      parameter (one = 1.0)
c
      tge = int(rvek(1) + 0.5)
      if ( (tge.lt.299) .and. (num.ne.4) ) then
         write(iws,'(t5,a)') '>> IN-PUT ERROR: geometric parameters'
         errorFlag = 1
         return
      elseif ( (tge.gt.299) .and. (num.ne.5) ) then
         write(iws,'(t5,a)') '>> IN-PUT ERROR: geometric parameters'
         errorFlag = 1
         return
      endif
c
      pi = 4.0*atan(one)
c
c  Basic Plate
c
      if (tge.eq.101) then
         w0 = rvek(2)
         l0 = rvek(3)
         t0 = rvek(4)
c
c  Cylinder
c
      elseif ((tge.gt.200).and.(tge.lt.300)) then
         rc = rvek(2)
         l  = rvek(3)
         t0 = rvek(4)
         if (tge.eq.201) then
            l0 = pi * rc
            w0 = l
         elseif (tge.eq.202) then
            l0 = pi * ( rc + t0 )
            w0 = l
         elseif (tge.eq.203) then
            l0 = l
            w0 = pi * rc
         elseif (tge.eq.204) then
            l0 = l
            w0 = pi * ( rc + t0 )
         else
            write(iws,'(t5,a)') '>> ERROR: cylider input'
            errorFlag = 1
            return
         endif
         rgeo(4) = rc
         rgeo(5) = l
c
c  Pipe - Elbow
c
      elseif ((tge.gt.300).and.(tge.lt.400)) then
         rc   = rvek(2)
         t0   = rvek(3)
         rp   = rvek(4)
         phi0 = rvek(5)
         if (tge.eq.301) then
            l0 = pi * rc
            w0 = phi0 * ( rp + rc )
         elseif (tge.eq.302) then
            l0 = pi * ( rc + t0 )
            w0 = phi0 * ( rp + rc + t0 )
         elseif (tge.eq.303) then
            l0 = phi0 * ( rp + rc )
            w0 = pi * rc
         elseif (tge.eq.304) then
            l0 = phi0 * ( rp + rc + t0 )
            w0 = pi * ( rc + t0 )
         elseif (tge.eq.305) then
            l0 = pi * rc
            w0 = phi0 * ( rp - rc )
         elseif (tge.eq.306) then
            l0 = pi * ( rc + t0 )
            w0 = phi0 * ( rp - rc - t0 )
         elseif (tge.eq.307) then
            l0 = phi0 * ( rp - rc )
            w0 = pi * rc
         elseif (tge.eq.308) then
            l0 = phi0 * ( rp - rc - t0 )
            w0 = pi * ( rc + t0 )
         else
            write(iws,'(t5,a)') '>> ERROR: pipe-elbow input'
            errorFlag = 1
            return
         endif
         rgeo(4) = rc
         rgeo(5) = rp
         rgeo(6) = phi0
      endif
c
      rgeo(1) = w0
      rgeo(2) = l0
      rgeo(3) = t0
c
      return
      end
c
c----67--1---------2---------3---------4---------5---------6---------712
c
      subroutine determ_mv_mh(mf,mv,mh,sfred_type)
      implicit none
      integer  mf,mv,mh,sfred_type, mfa
c
      if (sfred_type.eq.1) then
         if (mod(mf,2).ne.0) mf = mf + mod(mf,2)
         mfa = mf
      elseif (sfred_type.eq.2) then
         if (mod(mf,4).ne.0) mf = mf + mod(mf,4)
         mfa = mf/2
      elseif (sfred_type.eq.3) then
         if (mod(mf,6).ne.0) mf = mf + mod(mf,6)
         mfa = mf/3
      endif
      if (mfa.lt.8)  mfa=8
      if (mfa.gt.20) mfa=20
      if     (mfa .le. 8) then
         mv = 2
      elseif (mfa .le. 12) then
         mv = 3
      elseif (mfa .lt. 16) then
         mv = 4
      else
         mv = 5
      endif
      mh = (mfa-mv-mv)/2
      mf = sfred_type*mfa
      return
      end
c
c----67--1---------2---------3---------4---------5---------6---------712
c
      double precision function bias_crack(xtot,xn,n)
c
c Given: Xtot, xn, and n;  Xtot = x1+x2+...+xn
c Find : bias (=b below) in a geometric series  xn > x1
c
      implicit none
      integer  n,i,it
      double precision xtot,xn,rtot,b,db,tol,xi,r,rold
      parameter( tol = 1.e-6)
c
      rtot = xtot/xn
      b  = 0.8
      db = 0.2
      r  = 0.0
      it = 0
10    continue
         it = it + 1
         b = b + db
         rold = r
         r = 1.0
         do i=2, n
            xi = 1 - i
            r = r + b**xi
         enddo
         if (it.eq.1) goto 10
         if (abs((rtot-r)/rtot) .gt. tol) then
            if ( ((rtot-r)*(rtot-rold)) .lt. 0) db = - db/2.0
            goto 10
         endif
      bias_crack = b
      return
      end
c
c----67--1---------2---------3---------4---------5---------6---------712
c
      subroutine determ_remote_mesh(nmin,dx0,dxmax,xm,bias,x,n)
      implicit none
      integer  nmin,n, i,k
      double precision dx0,dxmax,xm,bias,x(*),xnmin,xstar,dx,dx1,xi
      double precision calc_lambda
c
c nmin = minimum elements allowed
c nmax = maximum elements allowed (currently not implemented)
c dx0 = initial minimum reference length
c dxmax = maximum element length allowed
c xm = total length to be divided into segments
c bias = initial ratio between adjacent elements
c
      if ((xm/real(nmin)) .lt. dx0) then
         n = nmin
         do i=1, n
            x(i) = xm*real(i)/real(n)
         enddo
         return
      endif
c
      dx1 = bias*dx0
      xi = nmin
      xnmin = dx1*(bias**xi-1.0)/(bias-1.0)
c
      if (xnmin .ge. xm) then
         n = nmin
         bias = exp( calc_lambda( dx0, xm+dx0, n+1) )
         dx1 = dx0*bias
         x(1) = dx1
         do i=2, n
            xi = i - 1
            x(i) = x(i-1) + dx1*bias**xi
         enddo
         return
      endif
c
      k = 1
      x(k) = dx0*bias
10    continue
         k = k + 1
         xi = k
         dx = dx0*bias**xi
         x(k) = x(k-1) + dx
      if (dx .lt. dxmax) goto 10
      xstar = x(k-1)
c
      if (xstar .gt. xm) then
         n = 0
20       continue
            n = n + 1
         if (x(n) .lt. xm) goto 20
         bias = exp( calc_lambda( dx0, xm+dx0, n+1) )
         dx1 = dx0*bias
         x(1) = dx1
         do i=2, n
            xi = i - 1
            x(i) = x(i-1) + dx1*bias**xi
         enddo
         return
      endif
c
      n = k
30    continue
         n = n + 1
         xstar = xstar + dxmax
      if (xstar .le. (xm+0.2*dxmax)) goto 30
      call calc_bias_int1(dx0,xm,n+1,k+1,bias)
      dx1 = dx0*bias
      x(1) = dx1
      do i=2, n
         xi = i-1
         if ( i .ge. k ) xi = k - 1
         x(i) = x(i-1) + dx1*bias**xi
      enddo
      return
      end
c
c----67--1---------2---------3---------4---------5---------6---------712
c
      subroutine Xdeterm_remote_mesh(bias,dx1,dxmax,xmax,x,n)
      implicit none
      integer  i,k,n
      double precision bias,dx1,dxmax,xmax,x(*),dx(200),xi
      double precision calc_lambda
      i = 1
      x(i) = dx1
10    continue
         i = i + 1
         xi = i - 1
         dx(i) = dx1*bias**xi
         x(i) = x(i-1) + dx(i)
      if (dx(i) .lt. dxmax) goto 10
      k = i - 1
c
c Case 1: Xgradation < Xmax
c
      if (x(k) .lt. xmax ) then
         i = k
20       continue
            i = i + 1
            x(i) = x(i-1) + dxmax
         if (x(i) .lt. xmax) goto 20
         n = i
30       continue
            xi = n - k
            dx(k+1) = (xmax-x(k)) / xi
         if (dx(k+1) .lt. dx(k)) then
            k = k - 1
            goto 30
         endif
         do i=k+1, n
            x(i) = x(i-1) + dx(k+1)
         enddo
      else
c
c Case 2: Xgradation > Xmax
c
         i = k
40       continue
            i = i - 1
         if (x(i) .gt. xmax) goto 40
         if ((xmax-x(i)) .le. (0.1+0.025*i*dx(i)) ) then
            n = i
         else
            n = i + 1
         endif
         bias = exp(calc_lambda(dx1,xmax,n))
         do i=2, n
            xi = i-1
            x(i) = x(i-1) + dx1*bias**xi
         enddo
      endif
      return
      end
c
c----67--1---------2---------3---------4---------5---------6---------712
c
      subroutine char_to_real(row,rvek,num)
c  The routine transform a character to a real number
      implicit   none
      integer    num,n,k,i,exponent,es, nr
      double precision  rvek(*),sign,xi,xe
      character  row*200
c
c  nr = to maximum position of characters in row
c
      nr = 200
      num=0
      n=1
 5    continue
         i=0
         k=0
         exponent=0
         es=1
         sign=1.0
 10      continue
c    Searcing for the 1:st digit !
         if ( ( (row(n:n).lt.'0').or.(row(n:n).gt.'9') ).and.
     &          (row(n:n).ne.'.') )  then
            n=n+1
            if (n.gt.nr) goto 99
            goto 10
         endif
c    Negative sign !
         if (row((n-1):(n-1)).eq.'-') sign=-1.0
 20      continue
c    Integer part
         if ( (row(n:n).ge.'0').and.(row(n:n).le.'9') ) then
            i=10*i+ichar(row(n:n))-48
            n=n+1
            if (n.gt.nr) goto 99
            goto 20
         endif
c
         if (row(n:n).ne.'.') goto 40
         n=n+1
         if (n.gt.nr) goto 99
 30      continue
c    Decimal part !
         if ( (row(n:n).ge.'0').and.(row(n:n).le.'9') ) then
            i=10*i+(ichar(row(n:n))-48)
            n=n+1
            if (n.gt.nr) goto 99
            k=k-1
            goto 30
         endif
 40      continue
         if ( (row(n:n).eq.'E').or.(row(n:n).eq.'e').or.
     &          (row(n:n).eq.'D').or.(row(n:n).eq.'d') ) then
c    Exponent
            n=n+1
            if (n.gt.nr) goto 99
            if (row(n:n).eq.'+') then
c    Positive or negative sign !
               es=1
               n=n+1
               if (n.gt.nr) goto 99
            elseif (row(n:n).eq.'-') then
               es=-1
               n=n+1
               if (n.gt.nr) goto 99
            endif
 50         if ( (row(n:n).ge.'0').and.(row(n:n).le.'9') ) then
               exponent=10*exponent+ichar(row(n:n))-48
               n=n+1
               if (n.gt.nr) goto 99
               goto 50
            endif
         endif
         xi = 1000*i
         k = k - 3
c         xi = i
c         k = k
         if (es.gt.0) then
            xe = k + exponent
         else
            xe = k - exponent
         endif
         num=num+1
      if (sign.gt.0) then
         rvek(num) =  xi*10.0**xe
      else
         rvek(num) = -xi*10.0**xe
      endif
      if (num.lt.nr) goto 5
 99   return
      end
c
c----67--1---------2---------3---------4---------5---------6---------712
c
      subroutine char_to_int(row,ivek,num)
c-- The routine transform a character to an integer number
      integer    ivek(*),num,n,i
      character  row*200
      logical    minus
c
c  nr = to maximum position of characters in row
c
      nr = 200
      num=0
      n=1
5     continue
         i=0
         minus = .false.
10        continue
c... Search for the 1:st digit
         if ( (row(n:n).lt.'0') .or. (row(n:n).gt.'9') ) then
            n=n+1
            if (n.gt.nr) goto 99
            goto 10
         endif
c... Negative sign !
         if (row((n-1):(n-1)).eq.'-') minus=.true.
20       continue
         if ( (row(n:n).ge.'0').and.(row(n:n).le.'9') ) then
            i=10*i+ichar(row(n:n))-48
            n=n+1
            if (n.gt.nr) goto 99
            goto 20
         endif
         num=num+1
         ivek(num)=i
         if (minus) ivek(num)= - ivek(num)
      if (num.lt.nr) goto 5
99    return
      end
c
c----67--1---------2---------3---------4---------5---------6---------712


text:
module inc

contains

integer(c_int) function f_inc(i) result(r) bind(c)
    use add
    use iso_c_binding, only: c_int
    integer(c_int), value :: i
    print *, 'f_inc()'
    r = f_add(i, 1)
end

end

max_stars_repo_name: okhlybov/cmakeworks

max_issues_repo_name: okhlybov/cmakeworks

max_forks_repo_name: okhlybov/cmakeworks


text:
      CHARACTER*35 FILNAM
      INTEGER PLTREC(750)
      LOGICAL LIBM, LSKIP
      DATA LSKIP/.FALSE./
C
      COMMON /NASTCOM/ ICMND,ICNTRL,IR,IS,IT,IU
      COMMON /SCOPE/ IBAUD
C
C     GET INPUT FILENAME & TYPE (VAX OR IBM)
10	  CALL NAMTYP(FILNAM,LIBM)
C     OPEN PLOTFILE
	  OPEN(UNIT=13,NAME=FILNAM,TYPE='OLD',
!    *         ACCESS='SEQUENTIAL',RECORDTYPE='FIXED',
!    *         RECORDSIZE=3000,BLOCKSIZE=3000,
!    *         FORM='FORMATTED',ORGANIZATION='SEQUENTIAL',
!    *	       READONLY,CARRIAGECONTROL='NONE', ERR=10)
     *	       READONLY,ERR=10)
C     INITIALIZE PLOT SCREEN
	  TYPE 1000, 'ENTER TERMINAL BAUD RATE:'
	  ACCEPT *, IBAUD
	  IBAUD = IBAUD/10  !! CHANGE TO CHAR/SEC
	  CALL INITT(IBAUD)
C
      DO 200 IREC=1,1000000
C	  GET NEXT PLOT RECORD
	      READ(13,2000,END=250) PLTREC
C	  CONVERT IBM RECORD
	      IF (LIBM) CALL CNVIBM(PLTREC)
	  DO 100 I=1,100
C	      GET NEXT COMMAND (TOTAL OF 100 COMMANDS/RECORD)
		  CALL GETCMD(PLTREC)
C	      EXECUTE COMMAND
		  IF (ICMND .EQ. 0)  CONTINUE     !! NULL COMMAND
		  IF (ICMND .EQ. 1)  CALL BGNPLT(LSKIP) !START NEW PLOT
		  IF (LSKIP)  GOTO 100
		  IF (ICMND .EQ. 2)  CONTINUE     !! SELECT CAMERA
		  IF (ICMND .EQ. 3)  CONTINUE     !! SKIP FRAME
		  IF (ICMND .EQ. 4)  CALL DRWCHR  !! TYPE A CHARACTER
		  IF (ICMND .EQ. 5)  CALL DRWLIN  !! DRAW A LINE
		  IF (ICMND .EQ. 6)  CALL DRWLIN  !! DRAW AN AXIS
		  IF (ICMND .EQ. 15) CALL DRWLIN  !! DRAW A LINE
		  IF (ICMND .EQ. 16) CALL DRWLIN  !! DRAW AN AXIS
100	      CONTINUE
200	  CONTINUE
C
C     END OF PLOT FILE
250	  CALL RESET
	  CALL SEETW(MINX,MAXX,MINY,MAXY)
	  CALL FINITT(MINX,MAXY)
	  STOP 'END OF PLOT FILE'
C
1000  FORMAT(1X,A25,1X,$)
2000  FORMAT(750A4)
      END
C------------------------------------------------------------------------------
      SUBROUTINE BGNPLT(LSKIP)
C
      CHARACTER*1 YESNO
      LOGICAL LSKIP
      COMMON /NASTCOM/ ICMND,ICNTRL,IR,IS,IT,IU
      COMMON /SCOPE/ IBAUD
C
C     DUMP BUFFER FOR PREVIOUS PLOT
	  CALL RESET
	  CALL SEETW(IMINX,IMAXX,IMINY,IMAXY)
	  CALL HOME
	  CALL ANMODE
	  PAUSE 'TYPE "CONTINUE" FOR NEXT PLOT:'
	  LSKIP = .FALSE.
	  TYPE 1000, 'SKIP NEXT FRAME', IR, '? (Y OR N):'
	  ACCEPT 2000, YESNO
	  IF (YESNO .EQ. 'Y')  LSKIP = .TRUE.
C
C     INITIALIZE FOR NEXT PLOT
	  CALL INITT(IBAUD)
	  CALL SEETW(IMINX,IMAXX,IMINY,IMAXY)
	  IMAXY = NINT(.95 * FLOAT(IMAXY))
	  XMAX = IS
	  YMAX = IT
	  RATIO = YMAX/XMAX
	  LENY = IMAXY
	  LENX = NINT( FLOAT(LENY) / RATIO )
	  IF (LENX .GT. IMAXX)  THEN
	      LENX = IMAXX
	      LENY = NINT( FLOAT(LENX) * RATIO )
	      END IF
C	  SET SCREEN WINDOW
	      MINX = (IMAXX-LENX)/2
	      MINY = (IMAXY-LENY)/2
	      CALL SWINDO(MINX,LENX,MINY,LENY)
C	  SET VIRTUAL (USER) WINDOW
	      CALL DWINDO(0.0,XMAX,0.0,YMAX)
      RETURN
C
1000  FORMAT(1X,A15,' (PLOT NO.',I4,')',A11,1X,$)
2000  FORMAT(A1)
      END
C------------------------------------------------------------------------------
      SUBROUTINE CNVIBM(RECORD)
      INTEGER RECORD(750)
      BYTE BYTE(4), BTEMP(4)
      EQUIVALENCE (IWORD,BYTE(1)) , (ITEMP,BTEMP(1))
C
C     CONVERT FROM IBM FORMAT TO VAX FORMAT:
C     REVERSE THE ORDER OF THE BYTES IN EACH WORD OF THE RECORD
	  DO 100 I=1,750
	      ITEMP = RECORD(I)
	      BYTE(1) = BTEMP(4)
              BYTE(2) = BTEMP(3)
              BYTE(3) = BTEMP(2)
	      BYTE(4) = BTEMP(1)
              RECORD(I) = IWORD
100           CONTINUE
      RETURN
      END
C------------------------------------------------------------------------------
      SUBROUTINE DRWCHR
C
C     THIS (CHARACTER-TYPING) CAPABILITY IS NOT YET IMPLEMENTED
C     FOR THIS PLOTTING PROGRAM.
C
      RETURN
      END
C------------------------------------------------------------------------------
      SUBROUTINE DRWLIN
C
      COMMON /NASTCOM/ ICMND,ICNTRL,IR,IS,IT,IU
C
C     MOVE BEAM TO STARTING POINT OF LINE
	  XSTART = IR
          YSTART = IS
	  CALL MOVEA(XSTART,YSTART)
C     DRAW LINE FROM STARTING POINT TO END POINT
	  XFINAL = IT
	  YFINAL = IU
	  CALL DRAWA(XFINAL,YFINAL)
      RETURN
      END
C------------------------------------------------------------------------------
      SUBROUTINE GETCMD(PLTREC)
      IMPLICIT INTEGER (A-Z)
      DIMENSION PLTREC(750),Q(30),MASK(4)
      COMMON /NASTCOM/ CMND,CNTRL,R,S,T,U
      COMMON /RECPNT/ KWORD,KBYTE
C     EQUATE COMMAND INTEGERS TO MNEMONIC VARIABLES
          EQUIVALENCE
     *     (Q( 1),PC), (Q( 2),CI),
     *     (Q( 3),R4), (Q( 4),R3), (Q( 5),R2), (Q( 6),R1), (Q( 7),R0),
     *     (Q( 8),S4), (Q( 9),S3), (Q(10),S2), (Q(11),S1), (Q(12),S0),
     *     (Q(13),T4), (Q(14),T3), (Q(15),T2), (Q(16),T1), (Q(17),T0),
     *     (Q(18),U4), (Q(19),U3), (Q(20),U2), (Q(21),U1), (Q(22),U0)
C     INITIALIZE BYTE MASKS, AND PLTREC WORD & BYTE POINTERS
          DATA MASK(1) /'000000FF'X/ ,
     *         MASK(2) /'0000FF00'X/ ,
     *         MASK(3) /'00FF0000'X/ ,
     *         MASK(4) /'FF000000'X/
	  DATA KWORD /1/, KBYTE /4/
C
C     BREAK OUT 30 INTEGERS OF NEXT COMMAND
          DO 100 I=1,30
              Q(I) = IAND( PLTREC(KWORD) , MASK(KBYTE) )
              Q(I) = ISHFT( Q(I) , -8*(KBYTE-1) )
C             SET POINTERS TO GET NEXT BYTE
                  KBYTE = KBYTE - 1
                  IF (KBYTE .EQ. 0) KWORD = KWORD + 1
                  IF (KBYTE .EQ. 0) KBYTE = 4
		  IF (KWORD .EQ. 751) KWORD = 1
100           CONTINUE
C     CALCULATE CMND,CNTRL,R,S,T,&U
          CMND = PC
          CNTRL = CI
          R = R0 + 10*R1 + 100*R2 + 1000*R3 + 10000*R4
          S = S0 + 10*S1 + 100*S2 + 1000*S3 + 10000*S4
          T = T0 + 10*T1 + 100*T2 + 1000*T3 + 10000*T4
          U = U0 + 10*U1 + 100*U2 + 1000*U3 + 10000*U4
      RETURN
      END
C------------------------------------------------------------------------------
      SUBROUTINE NAMTYP(FILNAM,LIBM)
      CHARACTER FILNAM*(*) , QFLAG*1
      LOGICAL LIBM
C
C     GET INPUT FILENAME
          TYPE 1000
	  ACCEPT 2000, FILNAM
C     IS THIS AN IBM FORMAT FILE ?
10	  TYPE 3000
	  ACCEPT 4000, QFLAG
	  IF (QFLAG .EQ. 'Y') THEN
	      LIBM = .TRUE.
          ELSE IF (QFLAG .EQ. 'N') THEN
              LIBM = .FALSE.
	  ELSE
C	      INPUT ERROR; TRY AGAIN
	      GOTO 10
	  ENDIF
      RETURN
C
1000  FORMAT(' ENTER NAME OF INPUT PLT2 FILE: ',$)
2000  FORMAT(A35)
3000  FORMAT(' CONVERT FILE FROM IBM FORMAT ? (Y OR N): ',$)
4000  FORMAT(A1)
      END

ext: for


text:
      subroutine setprob
      implicit double precision (a-h,o-z)
      character*12 fname
      common /comrp/ grav
      common /comsrc/ ndim
      common/cdisc/ x0,y0,alf,beta,r0,idisc
      common /comic/ hin,hout
c
c     # Set the material parameters for the acoustic equations
c     # Passed to the Riemann solver rp1.f in a common block
c
c
      iunit = 7
      fname = 'setprob.data'
c     # open the unit with new routine from Clawpack 4.4 to skip over
c     # comment lines starting with #:
      call opendatafile(iunit, fname)

c     # ndim = space dimensions (2 = cylindrical symmetry, 3 = spherical)
      read(7,*) ndim

c     # gravitational constant:
      read(7,*) grav

c     # data for radial dam-break problem:
      read(7,*) x0,y0,r0
      read(7,*) hin,hout
c
      return
      end

max_issues_repo_name: che-wenchao/D-Claw

max_forks_repo_name: che-wenchao/D-Claw


text:
program test
! Test the CIJ_global_VTI function

use anisotropy_ajn

implicit none

real(8) :: C(6,6),vp,vs,rho,xi,phi,eta

! At 2750 in AK135
vp = 13650.
vs = 7249.
rho = 5423.
xi = 1.3
phi = 0.8
eta = 1.3

C = CIJ_global_VTI(vp,vs,rho,xi,phi,eta)

write(*,*) C

end program test

max_stars_repo_name: JackWalpole/seismo-fortran-fork

max_issues_repo_name: JackWalpole/seismo-fortran-fork

max_forks_repo_name: JackWalpole/seismo-fortran-fork