text:
 MODULE ele18_db
   USE param_db,ONLY: mnam,midn,mlin
   USE mat_dba, ONLY: section,psecs,pmats,sect_search,mater,postv,pmater,snn
   USE c_input
   IMPLICIT NONE
     REAL(kind=8), PARAMETER :: sq3i =  0.577350269189626D+00

   ! Derived type for an ELE18 element
   ! Total Lagrangial (Log strain) 3-D brick element
   INTEGER(kind=4), PARAMETER :: nnode = 8    !number of nodes per element

   TYPE ele18
     INTEGER (kind=4) :: numel  ! label of element
     INTEGER (kind=4) :: matno  ! Material number
     INTEGER (kind=4) :: lnods(nnode)  ! Conectivities
     REAL (kind=8) :: angle(3)                  ! Euler angle between XYZ and orthotropic
     REAL (kind=8), POINTER :: dvol(:)          ! Initial volume
     REAL (kind=8), POINTER :: cartd(:,:,:)     !cartesian derivatives of SF
     REAL (kind=8), POINTER :: stint(:,:)       ! Actual Kirchhoff stresses
     REAL (kind=8), POINTER :: gausv(:,:)       !Gauss-point internal variables
     REAL (kind=8), POINTER :: nfdas(:,:),   &  ! Nodal Function Derivatives at the Assumed Strain points
                               jacin(:,:,:)     ! in-plane inverse jacobian
     TYPE (ele18), POINTER :: next              !pointer to next element
   END TYPE ele18

   TYPE ele18p
     TYPE(ele18), POINTER :: p
   END TYPE ele18p

   ! Derived type for a set of ELE18 elements
   TYPE ele18_set
     CHARACTER (len=mnam) :: sname ! set name
     INTEGER (kind=4) :: nelem, &  ! number of elements
                         nreqs, &  ! number of GP for hist. output
                         narch, &  ! number of output unit
                         ngaus     !number of integration points
     LOGICAL :: gauss           ! .FALSE. -> Initial constants not
                                !  defined or not updated
     LOGICAL :: small           ! .TRUE. -> use Green strain tensor
                                ! .FALSE. -> Use log strains (Default)
     INTEGER :: plstr           ! compute Plastic Strain Flag
         ! -1 from Cauchy stress  0 - do not   1 from 2nd P-K
     REAL (kind=8) :: angdf(3)    ! Default Euler angles between
                                  ! Global X-Y-Z and orthotropic system
     REAL(kind=8), POINTER :: gpc(:,:) !Gauss-points coordinates
     LOGICAL :: shell           ! .TRUE. -> modify transverse shear
                                ! .FALSE. -> use standard formulation
     LOGICAL :: bbar            ! .TRUE. -> use 1 GP for volume straian
                                ! .FALSE. -> use standard formulation
     LOGICAL :: check           ! .TRUE. -> modify transverse shear
     INTEGER :: locax           ! local x definition option
     REAL (kind=8) :: btscal      ! critical time scale factor

     TYPE (ele18), POINTER    :: head, tail !pointer to first and last elm.
     INTEGER (kind=4), POINTER :: ngrqs(:)  !gauss points for output
     TYPE (ele18_set), POINTER :: next      !pointer to next set
   END TYPE ele18_set
   TYPE (ele18_set), POINTER, SAVE :: head,tail !first and last elements sets

 CONTAINS
   SUBROUTINE ini_ele18 (head, tail)
     !initialize a list of ELE18 sets

     TYPE (ele18_set), POINTER :: head, tail

     NULLIFY (head, tail)

   END SUBROUTINE ini_ele18

   SUBROUTINE add_ele18 (new, head, tail)
     !This subroutine adds data to the end of the list
     !Dummy arguments
     TYPE (ele18_set), POINTER :: new, head, tail

     !Check if a list is empty
     IF (.NOT. ASSOCIATED (head)) THEN
       !list is empty, start it
       head => new
       tail => new
       NULLIFY (tail%next)

     ELSE
       !add a segment to the list
       tail%next => new
       NULLIFY (new%next)
       tail => new

     ENDIF
   END SUBROUTINE add_ele18

   SUBROUTINE srch_ele18 (head, anter, posic, name, found)
     !This subroutine searches for a set named "name"
     !Dummy arguments
     LOGICAL :: found
     CHARACTER (len=*) :: name ! set name
     TYPE (ele18_set), POINTER :: head, anter, posic

     found = .FALSE.                     !initializes flag
     NULLIFY (posic,anter)               !initializes pointers
     !Check if a list is empty
     IF (ASSOCIATED (head)) THEN         !if there are sets
       posic => head                     !point to first set
       DO
         IF (TRIM(posic%sname) == TRIM(name)) THEN    !check name
           found = .TRUE.                !set flag to .TRUE.
           EXIT                          !O.K. exit search
         END IF
         IF (ASSOCIATED(posic%next) ) THEN   !there is a next set
           anter => posic                    !point anter to present
           posic => posic%next               !point present to next
         ELSE
           EXIT                              !list exhausted, Exit search
         END IF
       END DO
     END IF
     IF (.NOT.found) NULLIFY (posic,anter)   !set not found, null pointers
   END SUBROUTINE srch_ele18

   SUBROUTINE del_ele18 (head, anter, posic)

     !This subroutine deletes a set pointed with posic

     TYPE (ele18_set), POINTER :: head, anter, posic

     TYPE (ele18), POINTER :: ea,ep
     INTEGER (kind=4) :: iel

     IF (.NOT.ASSOCIATED (anter)) THEN  !if anter pointer is null => head
       head => posic%next               !point first to next
     ELSE
       anter%next => posic%next         !skip posic
     END IF

     ! deallocation of the set memory is next done
     NULLIFY( ea )                  !nullify previous element in list
     ep => posic%head               !point present element to first
     DO iel = 1,posic%nelem         !for each element in the set
       CALL del_ele18e (posic%head,posic%tail, ea, ep )  !deletes element
     END DO

     NULLIFY (posic,anter)          !point to nothing
   END SUBROUTINE del_ele18

   ! ******* functions for a list of elements in a set ********

   SUBROUTINE ini_ele18e (head, tail)
     !initialize a list of ELE18 elements

     TYPE (ele18), POINTER :: head, tail

     NULLIFY (head, tail)       !initializes first and last pointer

   END SUBROUTINE ini_ele18e

   SUBROUTINE new_ele18e(elm)
   !Create a new element of ELE18 sets

     TYPE(ele18),POINTER:: elm

     ALLOCATE(elm)
     elm%numel = 0        !Initialize label of element
     elm%matno = 0        !     "     material number
     elm%lnods(1:8) = 0   !     "     conectivities
     elm%angle = 0d0      !     "     angle between dir 1 and orthotropic dir 1
     NULLIFY(elm%dvol,elm%cartd,elm%stint,elm%gausv,elm%nfdas,elm%jacin)
     NULLIFY(elm%next)

   RETURN
   END SUBROUTINE new_ele18e

   SUBROUTINE add_ele18e (new, head, tail)
     !This subroutine adds data to the end of the list
     !Dummy arguments
     TYPE (ele18), POINTER :: new, head, tail

     !Check if a list is empty
     IF (.NOT. ASSOCIATED (head)) THEN       !list is empty, start it
       head => new                           !first element
       tail => new                           !last element
       NULLIFY (tail%next)                   !last poit to nothing

     ELSE                                    !add a segment to the list
       tail%next => new                      !point to the new element
       NULLIFY (new%next)                    !nothing beyond the last
       tail => new                           !new last element

     END IF
   END SUBROUTINE add_ele18e

   SUBROUTINE srch_ele18e (head, anter, posic, kelem, found)
     !This subroutine searches for an element labeled "kelem"
     !Dummy arguments
     LOGICAL :: found
     INTEGER (kind=4) :: kelem
     TYPE (ele18), POINTER :: head, anter, posic

     found = .FALSE.             !initializes flag and pointers
     NULLIFY (posic,anter)

     IF (ASSOCIATED (head)) THEN          !Check if a list is empty
       posic => head                      !begin at top
       DO
         IF(posic%numel == kelem) THEN    !check if label found
           found = .TRUE.                 !Found
           EXIT                           !element found, EXIT
         END IF
         IF (ASSOCIATED(posic%next) ) THEN    !check if there are more els
           anter => posic                     !remember previous element
           posic => posic%next                !new present element
         ELSE
           EXIT                               !no more elements EXIT
         END IF
       END DO
     END IF
     IF (.NOT.found) NULLIFY (posic,anter)    !point to nothing
     RETURN
   END SUBROUTINE srch_ele18e

   SUBROUTINE del_ele18e (head, tail, anter, posic)

     !This subroutine deletes element pointed with posic

     TYPE (ele18), POINTER :: head, tail, anter, posic
     TYPE (ele18), POINTER :: e

     IF (.NOT.ASSOCIATED (anter)) THEN    !
       head => posic%next
     ELSE
       anter%next => posic%next
     END IF
     e => posic%next                       !keep pointer to next element
     IF( .NOT.ASSOCIATED(e) )tail => anter !last element in list
     DEALLOCATE (posic%dvol)               !deallocate vol arrays
     DEALLOCATE (posic%stint)              !deallocate stress array
     IF(ASSOCIATED(posic%gausv))DEALLOCATE (posic%gausv)  !deallocate variable arrays
     DEALLOCATE (posic%cartd)              !deallocate fixed space
     IF(ASSOCIATED(posic%nfdas))DEALLOCATE (posic%nfdas,posic%jacin)  !deallocate variable arrays
     DEALLOCATE (posic)                    !deallocate fixed space
     posic => e                            !point to next element
     ! NULLIFY (posic,anter)
     RETURN
   END SUBROUTINE del_ele18e

   SUBROUTINE cut_ele18e (head, anter, posic)
     !This subroutine deletes a set pointed with posic
     ! without nullifying anter    ???? what for ????
     TYPE (ele18), POINTER :: head, anter, posic

     IF (.NOT.ASSOCIATED (anter)) THEN
       head => posic%next
     ELSE
       anter%next => posic%next
     ENDIF
     NULLIFY (posic)
   END SUBROUTINE cut_ele18e

   include "acvd18.inc"
   include "bmat18.inc"
   include "bsma18.inc"
   include "comm18.inc"
   include "delt18.inc"
   include "dump18.inc"
   include "eige18.inc"
   include "elmd18.inc"
   include "expo18.inc"
   include "gaus18.inc"
   include "impo18.inc"
   !include "jacobi.inc"
   include "lgst18.inc"
   include "lcas18.inc"
   !include "luma18.inc"
   include "masm18.inc"
   include "mase18.inc"
   include "nods18.inc"
   include "outd18.inc"
   include "rest18.inc"
   include "resv18.inc"
   include "rubber3d.inc"
   include "rubber3dn.inc"
   include "slum18.inc"
   include "surf18.inc"
   include "secd18.inc"
   include "slno18.inc"
   include "updl18.inc"

 END MODULE ele18_db

max_stars_repo_name: jerebenitez/IFE-simpact-openfoam

max_issues_repo_name: jerebenitez/IFE-simpact-openfoam

max_forks_repo_name: jerebenitez/IFE-simpact-openfoam


text:
module slf_integrate

  use iso_fortran_env, only: r32 => real32, r64 => real64
  implicit none

  interface integrate
    module procedure integrate32
    module procedure integrate64
  end interface integrate

contains

  pure real(r32) function integrate32(y,x) result(intg)
    real(r32), dimension(:), intent(in) :: x, y
    integer :: n

    if (size(x) /= size(y)) error stop "size(x) /= size(y)"
    if (size(x) < 2) error stop 'size(x) < 2'

    n = size(x)
    intg = sum((y(2:n) + y(1:n-1)) * (x(2:n) - x(1:n-1))) / 2
  end function

  pure real(r64) function integrate64(y,x) result(intg)
    real(r64), dimension(:), intent(in) :: x, y
    integer :: n

    if (size(x) /= size(y)) error stop "size(x) /= size(y)"
    if (size(x) < 2) error stop 'size(x) < 2'

    n = size(x)
    intg = sum((y(2:n) + y(1:n-1)) * (x(2:n) - x(1:n-1))) / 2
  end function

end module slf_integrate

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

module slf_deriv

  use iso_fortran_env, only: r64 => real64
  implicit none

contains

  pure subroutine deriv(x,yin,yout)
    real(r64), intent(in) :: x(:), yin(:)
    real(r64), intent(out) :: yout(:)
    real(r64) :: dx2
    integer :: i,n

    n = size(yout)

    if ( n .lt. 3 ) then
      error stop 'n should be at least 3!'
    end if

    do i=2,n-1
      yout(i) = (yin(i+1)-yin(i-1))/(x(i+1)-x(i-1))
    end do
    dx2 = 2 * ( yin(3) + yin(1) - 2 * yin(2) ) / ( (x(3)-x(2))**2 + (x(2)-x(1))**2 )
    yout(1) = yout(2) + (x(1)-x(2))*dx2
    dx2 = 2 * ( yin(n) + yin(n-2) - 2 * yin(n-1) ) / ( (x(n)-x(n-1))**2 + (x(n-1)-x(n-2))**2 )
    yout(n) = yout(n-1) + (x(n)-x(n-1))*dx2
  end subroutine

  pure subroutine deriv2(x,yin,yout)
    real(r64), intent(in) :: x(:), yin(:)
    real(r64), intent(out) :: yout(:)
    integer :: i,n
    n = size(yout)

    do i=2,n-1
      yout(i) = 2 * ( yin(i+1) + yin(i-1) - 2 * yin(i) ) / ( (x(i+1)-x(i))**2 + (x(i-1)-x(i))**2 )
    end do
    yout(1) = yout(2) + ( x(1) - x(2) )*( yout(3)-yout(2) )/( x(3)-x(2) )
    yout(n) = yout(n-1) + ( x(n) - x(n-1) )*( yout(n-1)-yout(n-2) )/( x(n-1)-x(n-2) )
  end subroutine

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

  pure subroutine diffx(x, y)
    real(r64), intent(in) :: x(:)
    real(r64), intent(out) :: y(:)
    integer :: i, n

    if (size(x) /= size(y)) error stop 'diffx: size(x) /= size(y)'

    n = size(x)
    do concurrent (i = 2:n-1)
      y(i) = (x(i+1) - x(i-1)) / 2
    end do

    y(1) = x(2) - x(1)
    y(n) = x(n) - x(n-1)
  end subroutine

  pure function fdiffx(x) result(y)
    real(r64), intent(in) :: x(:)
    real(r64) :: y(size(x))

    call diffx(x, y)
  end function

   !----------------------------------------------------------------------------!
  ! computes log gradient of any function

  pure subroutine loggrad(x, y, d)
    real(r64), intent(in) :: x(:), y(:)
    real(r64), intent(out) :: d(:)
    integer :: i

    if (size(x) /= size(y) .or. size(d) /= size(y)) error stop 'loggrad: sizes dont match'

    do concurrent (i = 2:size(y)-1)
      d(i) = (y(i+1) - y(i-1)) * (x(i+1) + x(i-1)) &
      &   / ((y(i+1) + y(i-1)) * (x(i+1) - x(i-1)) )
    end do

    d(1) = d(2)
    d(size(d)) = d(size(d) - 1)
  end subroutine

end module

max_stars_repo_name: gronki/diskvert

max_issues_repo_name: gronki/diskvert

max_forks_repo_name: gronki/diskvert


text:
C----------------------------------------------------------------------
        subroutine ftcopy(iunit,ounit,moreky,status)

C       copies the CHDU from IUNIT to the CHDU of OUNIT.
C       This will also reserve space in the header for MOREKY keywords
C       if MOREKY > 0.

C       iunit   i  fortran unit number of the input file to be copied
C       ounit   i  fortran unit number of the output file to be copied to
C       moreky  i  create space in header for this many more keywords
C       status  i  output error status
C
C       written by Wm Pence, HEASARC/GSFC, Jan, 1992

        integer iunit,ounit,moreky,status

C       COMMON BLOCK DEFINITIONS:--------------------------------------------
        integer nb,ne
        parameter (nb = 20)
        parameter (ne = 200)
        integer bufnum,chdu,hdutyp,maxhdu,hdstrt,hdend,nxthdr,dtstrt
        integer nxtfld
        logical wrmode
        common/ft0001/bufnum(199),chdu(nb),hdutyp(nb),maxhdu(nb),
     &  wrmode(nb),hdstrt(nb,ne),hdend(nb),nxthdr(nb),dtstrt(nb),nxtfld
C       END OF COMMON BLOCK DEFINITIONS-----------------------------------

        integer ibuff,obuff,i,nkeys,nadd
        integer bitpix,naxis,naxes(99),pcount,gcount
        character hrec*80
        logical simple,extend

        if (status .gt. 0)return

        if (iunit .eq. ounit)then
                status=101
                return
        end if

        ibuff=bufnum(iunit)
        obuff=bufnum(ounit)

C       find out the number of keywords which exist in the input CHDU
        call ftghsp(iunit,nkeys,nadd,status)

C       copy the keywords one at a time to the output CHDU
        if ( (chdu(ibuff) .eq. 1 .and. chdu(obuff) .ne. 1) .or.
     &     (chdu(ibuff) .ne. 1 .and. chdu(obuff) .eq. 1) )then
C               copy primary array to image extension, or vise versa

C               copy the required keywords:
                simple=.true.
                extend=.true.
                call ftghpr(iunit,99,simple,bitpix,naxis,
     &          naxes,pcount,gcount,extend,status)
                if (status .gt. 0)return
                call ftphpr(ounit,simple,bitpix,naxis,
     &          naxes,pcount,gcount,extend,status)
                if (status .gt. 0)return

C               copy remaining keywords, excluding pcount, gcount and extend
                do 10 i=naxis+4,nkeys
                    call ftgrec(iunit,i,hrec,status)
                    if (hrec(1:8) .ne. 'PCOUNT  ' .and.
     &                  hrec(1:8) .ne. 'GCOUNT  ' .and.
     &                  hrec(1:8) .ne. 'EXTEND  ')then
                           call ftprec(ounit,hrec,status)
                    end if
10              continue
        else
C               just copy all the keys exactly from the input file to the output
                do 20 i=1,nkeys
                    call ftgrec(iunit,i,hrec,status)
                    call ftprec(ounit,hrec,status)
20              continue
        end if

C       reserve space for more keywords (if moreky > 0)
        call fthdef(ounit,moreky,status)

C       now ccopy the data from the input CHDU to the output CHDU
        call ftcpdt(iunit,ounit,status)

        end

ext: f

max_stars_repo_name: iraf-community/stsdas

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


text:
C
C	$Id: wmsect.f,v 1.4 2008-07-27 00:17:37 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE WMSECT(RECT1,RECT2,ISECT,RECT3)
C
C  Given two rectangles RECT1 and RECT2 (stored as XL,YB,XR,YT),
C  calculate the intersection of the two rectangles.  If the
C  intersection is not null, set ISECT to 1 and store the intersection
C  in RECT3.  If the intersection is null, set ISECT to 0 and return
C  all 0's in RECT3.
C
      DIMENSION RECT1(4),RECT2(4),RECT3(4)
C
      S2LLX = RECT2(1)
      S2LLY = RECT2(2)
      S2URX = RECT2(3)
      S2URY = RECT2(4)
      S1LLX = RECT1(1)
      S1LLY = RECT1(2)
      S1URX = RECT1(3)
      S1URY = RECT1(4)
C
      IF (S2LLX.LT.S1URX .AND. S1LLX.LT.S2URX .AND.
     +    S2LLY.LT.S1URY .AND. S1LLY.LT.S2URY)  THEN
        ISECT = 1
        RECT3(1) = MAX(S1LLX, S2LLX)
        RECT3(2) = MAX(S1LLY, S2LLY)
        RECT3(3) = MIN(S1URX, S2URX)
        RECT3(4) = MIN(S1URY, S2URY)
      ELSE
        ISECT = 0
        RECT3(1) = 0.
        RECT3(2) = 0.
        RECT3(3) = 0.
        RECT3(4) = 0.
      ENDIF
C
      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
      SUBROUTINE FRPRMN(P,N,FTOL,ITER,FRET)
      PARAMETER (NMAX=50,ITMAX=200,EPS=1.E-10)
      DIMENSION P(N),G(NMAX),H(NMAX),XI(NMAX)
      FP=FUNC(P)
      CALL DFUNC(P,XI)
      DO 11 J=1,N
        G(J)=-XI(J)
        H(J)=G(J)
        XI(J)=H(J)
11    CONTINUE
      DO 14 ITS=1,ITMAX
        ITER=ITS
        CALL LINMIN(P,XI,N,FRET)
        IF(2.*ABS(FRET-FP).LE.FTOL*(ABS(FRET)+ABS(FP)+EPS))RETURN
        FP=FUNC(P)
        CALL DFUNC(P,XI)
        GG=0.
        DGG=0.
        DO 12 J=1,N
          GG=GG+G(J)**2
C         DGG=DGG+XI(J)**2
          DGG=DGG+(XI(J)+G(J))*XI(J)
12      CONTINUE
        IF(GG.EQ.0.)RETURN
        GAM=DGG/GG
        DO 13 J=1,N
          G(J)=-XI(J)
          H(J)=G(J)+GAM*H(J)
          XI(J)=H(J)
13      CONTINUE
14    CONTINUE
      PAUSE 'FRPR maximum iterations exceeded'
      RETURN
      END

ext: for


text:
	PROGRAM GPLTLN
C************************************************************************
C* GPLTLN								*
C* This program sets up the graphics area and optionally draws a map,	*
C* lat/lon lines, and a title.						*
C**									*
C* Log:									*
C* I. Graffman/RDS	12/84						*
C* I. Graffman/RDS	 6/85   GEMPLT Version 3.1			*
C* I. Graffman/RDS	 5/88   DEVICE declaration			*
C* M. desJardins/GSFC	 7/88						*
C* G. Huffman/GSC	10/88	Error messages				*
C* S. Schotz/GSC	 8/90	Added lat/lon line, title, panel, 	*
C*				garea, proj features			*
C* K. Brill/NMC		10/91	PANEL*24  --> *48			*
C* S. Jacobs/EAI	11/92   Added call to GMESG and 'shrttl'	*
C* S. Jacobs/NMC	 3/94	Added satellite display routines	*
C* L. Williams/EAI	 3/94   Clean up declarations of variables	*
C* S. Jacobs/NMC	 6/94   DEVICE*12 --> *72			*
C* L. Williams/EAI	 7/94	Removed call to GPMUPD; added shrttl	*
C* P. Bruehl/Unidata	 8/94	Added animation				*
C* M. desJardins/NMC	 8/94	Added ST_FLST				*
C* S. Jacobs/NMC	 9/94	Added text; Added GG_STTL;		*
C*				Reorganized title plotting		*
C* J. Cowie/COMET	 1/95	Added SATFIL, RADFIL, image subsetting	*
C* S. Jacobs/NMC	 2/95	Moved IN_TEXT to before setting proj	*
C* J. Cowie/COMET	 8/95	Change GSATIM to IM_DROP, add IM_LUTF,	*
C*				use drop flag		 		*
C* S. Jacobs/NMC	10/95	Changed check for when to use GG_STTL	*
C* D. Plummer/NCEP	11/95	Added LUTFIL as a parameter             *
C* D. Keiser/GSC	12/95	Added STNPLT as a parameter		*
C* J. Lewis/AWC		 1/02	Added attribute tables to VG file parse	*
C* M. Li/SAIC		11/03	Added color bar for images		*
C* T. Piper/SAIC	08/04	Added gg_scal and mscale		*
C* T. Piper/SAIC	07/05	Added cplus variable			*
C* T. Piper/SAIC        01/08   Added GD_INIT; removed from IN_BDTA     *
C* L. Hinson/AWC        11/11   Added CTB_PFREAD  and CES_GTRTBL to load*
C*                              the prefs table and group table         *
C*                              Done to enable placement, if requested  *
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	LOGICAL		clear
	CHARACTER	map*72, device*72, title*72, panel*72, text*72,
     +			lat*72, lon*72, garea*72, proj*72, shrttl*72,
     +			satfil*72, radfil*72, lutfil*72, stnplt*72,
     +			marker*72, vgfile*72, plus*72, varr(3)*72,
     +			imcbar*72, mscale*(LLMXLN)
C*
	CHARACTER	ttlstr*72, ttlinp*72, label*20 , vgfl*72,
     +                  vgf2*72, vgf3*72
	REAL		rarr(2)
	LOGICAL		respnd, done
C*
	CHARACTER	imgfls(MXLOOP)*132, cplus*4, ucproj*72
	LOGICAL		first, proces
	DATA		cplus / '+' /
C-----------------------------------------------------------------------
C*  Initialize user interface and graphics.
C
	CALL IP_INIT  ( respnd, iperr )
	IF  ( iperr .ne. 0 )  THEN
	    CALL ER_WMSG  ( 'GPLTLN', -1, ' ', ier )
	    CALL SS_EXIT
	END IF
C
C*	Initialize preferences table.
C
	CALL CTB_PFREAD  ( iret )
	IF  ( iret .ne. 0 )  THEN
	    CALL ER_WMSG  ( 'GPLTLN', -1, ' ', ier )
	    CALL SS_EXIT
	END IF
C
C*	Initialize group id table
C
	CALL CES_GTRTBL  ( iret )
	IF  ( iret .ne. 0 )  THEN
	    CALL ER_WMSG  ( 'GPLTLN', -1, ' ', ier )
	    CALL SS_EXIT
	END IF
C
C*  Initialize graphics.
C
	CALL GG_INIT  ( 1, iret )
	IF  ( iret .eq. 0 )  THEN
C
C*  Initialize grid library common area grdcmn.cmn
C
	    CALL GD_INIT  ( iret )
	    done = .false.
	  ELSE
	    CALL ER_WMSG  ( 'GPLTLN', -3, ' ', ier )
	    done = .true.
	END IF
	CALL IP_IDNT ( 'GPLTLN', ier )
C
	DO WHILE  ( .not. done )
C
	    proces = .true.
	    first = .true.
	    numimg = 1
C
C*	    Get input parameters.
C
	    CALL GPTINP  ( device, map, mscale, garea, proj, satfil,
     +			   radfil, panel, title, text, imcbar, lat,
     +			   lon, clear, marker, lutfil, stnplt, vgfile,
     +			   plus, iperr )
	    IF  ( iperr .eq. 0 )  THEN
C
C*		Set device and projection.
C
		CALL GG_SDEV  ( device, ier )
		IF  ( ier .ne. 0 )  proces = .false.
		IF  ( proces )  THEN
		  CALL IN_TEXT ( text, ier )
C
C*		  If projection=SAT or RAD, see if multiple image
C*		  files have been specified.
C
		  CALL ST_LCUC ( proj, ucproj, ier )
		  IF  ( ucproj (1:3) .eq. 'SAT' )  THEN
		      CALL ST_FLST  ( satfil, ';', ' ', MXLOOP, imgfls,
     +				      numimg, ier )
		  ELSE IF  ( ucproj (1:3) .eq. 'RAD' )  THEN
		      CALL ST_FLST  ( radfil, ';', ' ', MXLOOP, imgfls,
     +				      numimg, ier )
		  END IF
C
C*		  Set map projection
C
		  CALL GG_MAPS ( proj, garea, imgfls (1), idrpfl, ier )
		  IF  ( ier .ne. 0 )  proces = .false.
C
C*		  Get the size and width of the '+' sign.
C
		  CALL ST_RLST ( plus, '/', RMISSD, 2, rarr, num, ier )
		  IF ( rarr(1) .eq. RMISSD ) THEN
		    splus = 1
		  ELSE
		    splus = rarr(1)
		  END IF
C*
		  IF ( rarr(2) .eq. RMISSD ) THEN
		    iwplus = 1
		  ELSE
		    iwplus = NINT ( rarr(2) )
		  END IF
C
C*		  Start loop over input image files.
C
		  IF  ( proces )  THEN
	      	     DO  ifile = 1, numimg
C
C*			Reset the projection for each image.
C
			IF ( ucproj (1:3) .eq. 'SAT' .or.
     +			     ucproj (1:3) .eq. 'RAD' )
     +                       CALL GG_MAPS ( proj, garea, imgfls (ifile),
     +          			    idrpfl, ier )

			IF ( ifile .eq. 1 ) THEN
C
C*			    Go to the first frame.
C
		  	    first = .true.
			    CALL GSTANM ( ier )
			  ELSE
C
C*			    Advance to the next frame.
C
			    first = .false.
			    CALL GSPLOT ( ier )
			END IF
C
C*			Display user options, allow program exit.
C
C			IF  ( first )
C     +			    CALL GPMOPT ( device, proj, garea, map,
C     +					  title, panel, latlon, clear,
C     +					  ier )
C
			IF  ( ier .ne. 0 )  proces = .false.
C
			IF  ( proces )  THEN
C
C*			    Clear the screen, if requested, and set
C*			    the panel.
C
			    IF  ( clear ) CALL GCLEAR ( iret )
			    CALL GG_PANL ( panel, iret )
C
C*			    Display satellite image, if desired.
C
			    IF  ( ( idrpfl .eq. 1 ) .or.
     +			          ( idrpfl .eq. 0 .and. clear ) ) THEN
     				CALL IM_DROP ( iret )
				CALL IM_CBAR ( imcbar, iret )
			    END IF
C
C*			    Apply LUT file
C
			    IF  ( ifile .eq. 1 )
     +				CALL IM_LUTF ( lutfil, ier )
C
C*			    Draw map and station ID/marker.
C
			    CALL GG_MAP  ( map, iret )
			    CALL IN_MARK  ( marker, mkcolr, iret )
			    CALL GSCOLR ( mkcolr, ier )
			    CALL GG_SCAL ( mscale, ier )
C
C*			    Draw the latitude lines.
C
			    CALL GQBND ( 'M', blatmn, blonmn, blatmx,
     +					 blonmx, ier )
C*
			    IF ( ( blatmn .lt. -80. ) .and.
     +				 ( blatmx .lt. -80. ) ) THEN
				rstrt = -999.
			    ELSE IF ( ( blatmn .lt. -80. ) .and.
     +				      ( blatmx .gt. -80. ) ) THEN
				rstrt = -80.
			    ELSE IF ( ( blatmn .gt. -80. ) .and.
     +				      ( blatmn .lt. 80. ) ) THEN
				rstrt = INT ((blatmn - 5) / 10) * 10.
			    ELSE IF ( blatmn .gt. 80. ) THEN
				rstrt = -999.
			    END IF
C*
			    IF ( blatmx .lt. -80. ) THEN
				rend = -999.
			    ELSE IF ( ( blatmx .gt. 80. ) .and.
     +				      ( blatmn .gt. 80. ) ) THEN
				rend = -999.
			    ELSE IF ( ( blatmx .gt. -80. ) .and.
     +				      ( blatmx .lt. 80. ) ) THEN
				rend = INT ((blatmx + 10) / 10) * 10.
			    ELSE IF ( ( blatmn .lt. 80. ) .and.
     +				      ( blatmx .gt. 80. ) ) THEN
				rend = 80.
			    END IF

			    blonstr = INT ( blonmn / 30. ) * 30. - 30.
			    blonend = INT ( blonmx / 30. ) * 30. + 30.
			    IF ( blonstr .lt. -180. ) blonstr = -180.
			    IF ( blonend .gt. 180. ) blonend = 180.
C
C Draw markers for the length of each latitude line, only up till -60/60
C latitude, past that we use a different marker spacing
C
                            rmstrt = rstrt
                            rmend  = rend
                            IF ( rmstrt .lt. -60. ) rmstrt = -60.
                            IF ( rmend  .gt.  60. ) rmend  =  60.
			    DO rii = rmstrt, rmend, 10
			      DO rjj = blonmn, blonmx, 1
			          CALL GMARK ( 'M', 1, rii, rjj, ier )
			      END DO
			    END DO
C*
C Draw latitude lines at every 5 past 65, but only put markers every 5
C
			    IF ( blatmn .lt. -65 ) THEN
			      DO rii = -70., -85., -5
			        DO rjj = blonmn, blonmx, 5
			          CALL GMARK ( 'M', 1, rii, rjj, ier )
			        END DO
			      END DO
			    END IF
C*
C Draw latitude lines at every 5 past 65, but only put markers every 5
C
			    IF ( blatmx .gt. 65 ) THEN
			      DO rii = 70., 85., 5
			        DO rjj = blonmn, blonmx, 5
			          CALL GMARK ( 'M', 1, rii, rjj, ier )
			        END DO
			      END DO
			    END IF
C
C Draw markers in the center of each square
C
C			    DO rii = rstrt, rend, 10
C			      DO rjj = blonstr, blonend, 10
C			        CALL GMARK ( 'M', 1, rii+5, rjj+5, ier )
C			      END DO
C			    END DO

C
C Draw latitude lables foreach ring
C
			    DO rii = rstrt, rend, 10
                              ires = MOD(NINT(rii),30)
                              IF ( ires .eq. 0) THEN
			        DO rjj = blonstr, blonend, 30
                                   CALL ST_INCH ( ABS(NINT(rii)), label,
     +                                             ier )
			            CALL ST_LSTR ( label, ilen, ier )
C			            IF ( rii .gt. 0. ) THEN
C			              label = label(:ilen) // 'N'
C			            ELSE IF ( rii .lt. 0. ) THEN
C			              label = label(:ilen) // 'S'
C			            END IF
C			            CALL ST_LSTR ( label, ilen, ier )
                                    rot = 0.
			            CALL GTEXT ( 'M',rii,rjj+5.,label,
     +			          		rot, 0, 0, ier )
			        END DO
                              END IF
			    END DO
C
C*			    Draw the longitude lines.
C
			    rstrt = INT ( ( blonmn - 5 ) / 10 ) * 10
			    IF ( rstrt .lt. -180 ) rstrt = -180.
			    rend =  INT ( ( blonmx + 10 ) / 10 ) * 10
			    IF ( rend .ge. 180. ) rend = 170.
			    IF ( ( blatmn .lt. -85 ) .and.
     +				 ( blatmx .lt. -85 ) ) THEN
				rlstrt = -999.
			    ELSE IF ( ( blatmn .lt. -85 ) .and.
     +				      ( blatmx .gt. -85 ) ) THEN
				rlstrt = -85.
			    ELSE IF ( ( blatmn .gt. -85 ) .and.
     +				      ( blatmn .lt. 85 ) ) THEN
				rlstrt = INT ((blatmn - 5) / 10) * 10.
			    ELSE IF ( blatmn .gt. 85 ) THEN
				rlstrt = -999.
			    END IF
C*
			    IF ( blatmx .lt. -85 ) THEN
				rlend = -999.
			    ELSE IF ( ( blatmx .gt. 85 ) .and.
     +				      ( blatmn .gt. 85 ) ) THEN
				rlend = -999.
			    ELSE IF ( ( blatmx .gt. -85 ) .and.
     +				      ( blatmx .lt. 85 ) ) THEN
				rlend = INT ((blatmx + 10) / 10) * 10.
			    ELSE IF ( ( blatmn .lt. 85 ) .and.
     +				      ( blatmx .gt. 85 ) ) THEN
				rlend = 85.
			    END IF
C*
C Draw markers for longitude lines and at the same time draw
C longitude labels, at 25 & 55 north and south
C
			    DO rii = rstrt, rend, 10
			      DO rjj = rlstrt, rlend, 1
			        CALL GMARK ( 'M', 1, rjj, rii, ier )
                                IF (( rjj .eq. -55. ) .or.
     +                              ( rjj .eq. -25. ) .or.
     +                              ( rjj .eq.  25. ) .or.
     +                              ( rjj .eq.  55. )) THEN
                                  ires = MOD(NINT(rii),30)
                                  IF (ires .eq. 0) THEN
				    CALL ST_INCH ( ABS(NINT(rii)), label,
     +				              	 ier )
				    CALL ST_LSTR ( label, ilen, ier )
C				    IF ( rii .gt. 0. ) THEN
C				      label = label(:ilen) // 'E'
C				    ELSE IF ( rii .lt. 0. ) THEN
C				      label = label(:ilen) // 'W'
C				    END IF
C				    CALL ST_LSTR ( label, ilen, ier )
                                    rot = 0.
			            CALL GTEXT ( 'M', rjj, rii, label,
     +						  rot, 0, 0, ier )
				  END IF
				END IF
			      END DO
			    END DO

C
C Draw crosshairs at the pole if it is visible on the map
C
C			    IF ( blatmn .lt. -65 ) THEN
C			      DO rii = -180, 170, 90
C			        DO rjj = -90., -80, 1
C			          CALL GMARK ( 'M', 1, rjj, rii, ier )
C			        END DO
C			      END DO
C			    END IF
C			    IF ( blatmx .gt. 75 ) THEN
C			      DO rii = -180, 170, 90
C			        DO rjj = 80., 90, 1
C			          CALL GMARK ( 'M', 1, rjj, rii, ier )
C			        END DO
C			      END DO
C			    END IF
C
C*			    Draw the pluses.
C
			    CALL GQTEXT ( itxfn, itxhw, sztext, itxwid,
     +                                    ibrdr, irrotn, ijust, iret )
			    CALL GSTEXT ( 0, 0, splus, iwplus,
     +                                    111, 2, 2, iret )
			    rstrt = INT ( ( blonmn - 5 ) / 10 ) * 10
			    IF ( rstrt .lt. -180. ) rstrt = -180.
			    rend =  INT ( ( blonmx + 10 ) / 10 ) * 10
			    IF ( rend .ge. 180. ) rend = 170.
			    rlstrt = INT ( ( blatmn - 5 ) / 10 ) * 10
			    IF ( rlstrt .lt. -60. ) rlstrt = -60.
			    rlend =  INT ( ( blatmx + 10 ) / 10 ) * 10
			    IF ( rlend .ge. 60. ) rlend = 60.
			    DO rii = rstrt, rend, 10
			     DO rjj = rlstrt, rlend, 10
			       CALL GTEXT ( 'M', rjj, rii, cplus, 0.,
     +					    0, 0, ier )
			     END DO
			    END DO
			    CALL GSTEXT ( itxfn, itxhw, sztext, itxwid,
     +					  ibrdr, irrotn, ijust, iret )
C
C
C*			    Plot station markers.
C
			    CALL GG_SPLT ( stnplt, iret )
C
C*                          Plot VGF file.
C
                            IF  ( vgfile .ne. ' ' )  THEN
                                CALL ST_CLST ( vgfile, '|', ' ', 3,
     +                                          varr, num, ier )
				CALL ST_NULL ( varr(1), vgfl, lenf, ier )
				CALL ST_NULL ( varr(2), vgf2, lenf, ier )
				CALL ST_NULL ( varr(3), vgf3, lenf, ier )
				icol = 0
				CALL CDS_INIT ( iret )
C
C*                              Set scaling factors.
C
                                IF  ( varr(2) .ne. ' ' )  THEN
                                    CALL CDS_SCAL ( vgf2, ier )
                                END IF
C
C*                              Set user attributes.
C
                                IF  ( varr(3) .ne. ' ' )  THEN
                                    CALL CDS_RTBL ( vgf3, ier )
                                END IF
C
C*                              Plot the contents of the VG file.
C
                                CALL GG_DVGF ( vgfl, icol, iret )
                                IF  ( iret .ne. 0 )  THEN
                                    CALL ER_WMSG  ( 'GG', iret, vgfile,
     +                                              ier )
                                END IF
C
C*                              Reset scaling factors.
C
                                IF  ( varr(2) .ne. ' ' )  THEN
                                    CALL CDS_RESS ( ier )
                                END IF
                            END IF
C
C*                          Reset the text attributes.
C
                            CALL IN_TEXT ( text, ier )
C
C*			    Decode title input and draw the title.
C
			    ipbar = INDEX ( title, '|' )
			    IF  ( ipbar .ne. 0 )  THEN
				shrttl = title (ipbar+1:)
				IF  ( ipbar .eq. 1 )  THEN
				    ttlinp = ' '
				  ELSE
				    ttlinp = title (:ipbar-1)
				END IF
			      ELSE
				CALL ST_LSTR ( garea, len1, ier )
				shrttl = 'MAP OF AREA: ' // garea(:len1)
				ttlinp = title
			    END IF
			    CALL IN_TITL ( ttlinp, -3, ititl, linttl,
     +					   ttlstr, iret )
			    IF  ( ( ucproj (1:3) .eq. 'SAT' .or.
     +				    ucproj (1:3) .eq. 'RAD' ) .and.
     +				  ( ttlstr .eq. ' ' ) )  THEN
				CALL GG_STTL ( ttlstr, iret )
			    END IF
			    IF  ( clear ) CALL GMESG ( shrttl, iret )
			    IF  ( ititl .gt. 0 )  THEN
				CALL GSCOLR  ( ititl, iret )
				CALL GG_WSTR ( ttlstr, linttl, iret )
			    END IF
C
C*			    Flush the graphics buffer.
C
			    CALL GEPLOT  ( iret)
			END IF
		      END DO
C
C*		      Mark the end of the animation sequence.
C
		      CALL GENANM ( iret)
		  END IF
		END IF
	      END IF
C
C*	    Call the dynamic tutor.
C
	    CALL IP_DYNM ( done, iret )
	END DO
C*
	IF  ( iperr .ne. 0 )  CALL ER_WMSG  ( 'GPLTLN', iperr, ' ', ier )
	CALL GENDP   ( 0, iret )
	CALL IP_EXIT ( iret )
C*
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
subroutine fgmax_frompatch(mx,my,meqn,mbc,maux,q,aux,dx,dy, &
           xlower,ylower,level,time,time_afterstep)

    ! Do the new fixed grid stuff on all fgrids, updating
    ! based on the patch passed in.

    use fgmax_module
    use geoclaw_module, only: sea_level, dry_tolerance

    implicit none
    integer, intent(in) :: mx,my,meqn,mbc,maux,level
    real(kind=8), intent(in) :: q(meqn, 1-mbc:mx+mbc, 1-mbc:my+mbc)
    real(kind=8), intent(in) :: aux(maux, 1-mbc:mx+mbc, 1-mbc:my+mbc)
    real(kind=8), intent(in) :: dx,dy,xlower,ylower,time,time_afterstep

    real(kind=8), allocatable, dimension(:,:) :: fg_values
    logical, allocatable, dimension(:) :: mask_fgrid
    type(fgrid), pointer :: fg
    integer :: ifg,k,mv
    real(kind=8) :: h,B,eta

    !print *, '+++ FG_DEBUG = ', FG_DEBUG

    if (FG_num_fgrids == 0) then
        return
        endif

    !write(61,*) '++++ In frompatch, level,mx,my,xlower,ylower'
    !write(61,*) level,mx,my,xlower,ylower
    do ifg=1,FG_num_fgrids
        fg => FG_fgrids(ifg)
        !write(6,*) '++++ ifg = ',ifg
        !write(6,*) '++++ fg%xNbb = ',fg%x1bb,fg%x2bb
        !write(61,*) '++++ frompatch xNbb', fg%x1bb, fg%x2bb

        if (FG_DEBUG) then
            write(61,61) ifg,level,time
 61         format('---------- In fgmax_frompatch ----------',/, &
               'ifg = ',i2,' level = ',i1,' time = ',d16.6)
            if (time_afterstep <= minval(fg%t_last_updated)+fg%dt_check) then
                write(61,68) time, minval(fg%t_last_updated)
     68         format('++++ Skipping update at t = ', e20.11,' min t_last = ',&
                       e20.11)
            else
                write(61,67) time, minval(fg%t_last_updated)
     67         format('++++ Doing update at t = ', e20.11,' min t_last = ', &
                       e20.11)
                endif
            endif

        if ((time >= fg%tstart_max) .and. (time <= fg%tend_max) .and. &
                (level >= fg%min_level_check) .and. &
                (level >= minval(fg%levelmax)) .and. &
                (time_afterstep > minval(fg%t_last_updated)+fg%dt_check)) &
                then
            ! Otherwise this level won't update any fg%valuemax elements.
            !write(61,*) '+++ Setting fg_values'

            allocate(mask_fgrid(1:fg%npts))
            allocate(fg_values(FG_NUM_VAL, 1:fg%npts))

            ! Interpolate from q on patch to desired values fg_values on fgrid.
            call fgmax_interpolate(mx,my,meqn,mbc,maux,q,aux, &
                 dx,dy,xlower,ylower,ifg,level,fg_values,mask_fgrid,fg%npts)

            if (FG_DEBUG) then
                write(61,*) '+++ updating level: ',level
                write(61,*) '+++ in fgmax_frompatch, count = ', count(mask_fgrid)
                endif

            do k=1,fg%npts
                ! fg_values is set only at points k where the fgrid intersects the
                ! patch, indicated by mask_fgrid(k) == .true.
                ! Also: only update valuemax if the current patch is at least
                ! as high a level as the patch last used to update:
                if (mask_fgrid(k) .and. (level >= fg%levelmax(k))) then
                    fg%t_last_updated(k) = time
                    do mv=1,FG_NUM_VAL
                        !print *,'+++ updating fg%valuemax(mv,k),fg_values(mv,k):',&
                        !      fg%valuemax(mv,k),fg_values(mv,k)
                        if ((level > fg%levelmax(k)) .or. &
                              (fg_values(mv,k) > fg%valuemax(mv,k))) then
                            fg%valuemax(mv,k) = fg_values(mv,k)
                            ! also keep track of time maximum happened:
                            fg%tmax(mv,k) = time
                            endif
                        enddo
                    fg%levelmax(k) = level
                    endif

                ! keep track of arrival time...
                ! This might not be correct when there is subsidence or uplift?
                h = fg_values(1,k)
                B = fg%aux(level,1,k)
                eta = h+B
                if (mask_fgrid(k) .and. &
                   (level >= fg%min_level_check) .and. &
                   (h>dry_tolerance) .and. &
                   (eta > sea_level + fg%arrival_tol) .and. &
                   (fg%arrival_time(k) == FG_NOTSET)) then
                        fg%arrival_time(k) = time
                    endif

                enddo
                deallocate(mask_fgrid, fg_values)
            endif
        enddo

end subroutine fgmax_frompatch

max_stars_repo_name: delgadom/geoclaw

max_issues_repo_name: delgadom/geoclaw

max_forks_repo_name: BrisaDavis/geoclaw


text:
      subroutine sonfath (ncorp, vnumnpl, ifath, nsons, tfath,
     &                    numpe, numnp,   maxnode)

c     This program writes to geom.dat.<1..numpe> the number sons and
c     fathers required for averaging in the z-direction (Turbulent
c     Boundary Layer), (x,y)-space (Turbulent Channel Flow),
c     or (x,y,z)-space (Isotropic Turbulence)

      integer inform,n1m,n2m,n3m,numpe,numpe1,numnp,maxnode, tfath,
     &     counter,nc8,nc4,node,ifil,itmp
      integer gfath(0:numnp),ifath(numpe,maxnode),ncorp(numpe,maxnode),
     &     nsons(tfath),dummy(numnp),iper(numnp), vnumnpl(numpe)
      integer, allocatable :: imap(:),invmap(:)
      allocate (imap(numnp))
      allocate (invmap(numnp))

      write(*,*)'Enter number of nodes in each direction'
      read(*,*)n1m,n2m,n3m

      write(*,*)'For Z averaging enter 1 For X-Y averaging enter 2'
      write(*,*)'For X-Y-Z averaging enter 3  W=0 SPEBC enter 4'
      write(*,*)'fully local model 5'
      read(*,*)inform

c
c.... read and generate map
c
      kmap = 15

      open(kmap,file='geom.kmap',status='old',err=183)

      do i=1,numnp
         read(kmap,*,err=183) iorig
         imap(i)=iorig+1  ! imap takes arguement of original and returns new
         invmap(iorig+1)=i ! invmap takes arguemnt of new and returns original
      enddo
      goto 184
 183  write(*,*) 'geom.kmap not read properly...assuming inform=5'
      inform=5
 184  continue

      nsoncor=1
      if(inform.eq.4) then
         inform=1
         nsoncor=0
      endif

      gfath(:) = 0

      nc4=n3m*(n2m-1)+1
      nc8=n3m*n2m*n1m+1-n3m

      if (numpe.eq.1)then
         numnp   = n1m*n2m*n3m
         maxnode = numnp
      endif

      if (inform.eq.1) then ! For averaging over Z - lines

         do i = 1, tfath
            nsons(i) = n3m-nsoncor
         enddo

         counter = 1
         do i = 1, nc8, n3m
            do j = i, i+n3m-1
               gfath(j) = counter
            enddo
            counter = counter + 1
         enddo

         if(numpe.gt.1)then
            do j = 1, numpe
               do i = 1, maxnode
                  ifath(j,i) = gfath(imap(ncorp(j,i)))
               enddo
            enddo
         else
            do i = 1, maxnode
               ifath(1,i) = gfath(imap(i))
            enddo
         endif

      endif

      if (inform.eq.2) then ! For averaging over X-Z planes

         do i = 1, tfath
            nsons(i) = n1m*n3m - (n1m+n3m-1)
         enddo

         counter = 1
         do j = 1, nc4, n3m
            do i = 1, n1m
               do k = j, j+n3m-1
                  node = (i-1)*n2m*n3m + k
                  gfath(node) = counter
               enddo
            enddo
            counter = counter + 1
         enddo

         if(numpe.gt.1)then
            do j = 1, numpe
               do i = 1, maxnode
                  ifath(j,i) = gfath(imap(ncorp(j,i)))
               enddo
            enddo
         else
            do i = 1, maxnode
               ifath(1,i) = gfath(imap(i))
            enddo
         endif

      endif

      if( inform.eq.3 )then !For averaging over X-Y-Z

         open (unit=96,FILE='geom.iper',form='unformatted',
     &        status='old')
         read(96) (iper(l),l=1,numnp)
         close(96)

         dummy(:) = 1

         do i = 1, numnp
            if(iper(i).eq.0)then
               iper(i)=i
            endif
         enddo
         do j = 1, numnp
            i = iper(j)
            if (i .ne. j) then
               dummy(j) = 0
            endif
         enddo
         itmp=sum(dummy)
         nsons(tfath) = itmp
         write(*,*)'nsons(1)=',nsons(1)
         ifath=1
      endif

      if(inform.eq.5) then      ! local clipping

         do i = 1, numpe
            ifil = 99

            nsons=1
            do j=1,vnumnpl(i)
               ifath(i,j)=j
            enddo

         enddo
      endif

      return
      end


text:
$ Declare
   SUBROUTINE SIMConditions (a, c, B, W, err)
!
! Stress intensity magnification factor Mm
! Mm according to BS 7910:2005 M.3 Flat Plates
!
! Parameters
   USE PublicPar
   IMPLICIT NONE
!
   DOUBLE PRECISION, INTENT(IN)    :: a         ! FT = 2: Crack depth
   DOUBLE PRECISION, INTENT(IN)    :: c         ! FT = 2: Crack width
   DOUBLE PRECISION, INTENT(IN)    :: B         ! FT = 2: Plate thickness
   DOUBLE PRECISION, INTENT(IN)    :: W         ! FT = 2: Plate width
   INTEGER, INTENT(OUT)            :: err       ! Error message, 1 = error
!
! Declarations
   DOUBLE PRECISION                :: M1, M2, M3
   DOUBLE PRECISION                :: as, cs
   DOUBLE PRECISION                :: g, radius, adivD, g1, g2
   DOUBLE PRECISION                :: ftheta
   DOUBLE PRECISION                :: fi
!
! Constants
   err = 0
!
!  cope
   IF (FlawType .EQ. 9) THEN
!
!  through-thickness deck plate crack
   ELSE IF (FlawType .EQ. 8) THEN
!
!  edge flaws in plates BS 7910:2005 M.3.5
   ELSE IF (FlawType .EQ. 7) THEN
!
!  surface flaws in bolts BS 7910:2005 M.6.3.2
   ELSE IF (FlawType .EQ. 6) THEN
      radius = B / 2.D0
	  adivD = a / (2.D0 * radius)
	  IF (adivD .GT. 0.65D0) THEN
		 err = 1
	  ELSE IF (adivD .LT. 0.004D0) THEN
		 err = 1
	  END IF
!
!  corner flaws in plates BS 7910:2005 M.3.6
   ELSE IF (FlawType .EQ. 5) THEN
	  IF ((a/c) .LT. 0.2D0) THEN
		 err = 1
	  ELSE IF ((a/c) .GT. 2.D0) THEN
		 err = 1
	  ELSE IF ((a/B) .GT. 1.D0) THEN
         err = 1
	  END IF
!
!  semi-circular surface flaw in round bar BS 7910:2005 M.3.x
   ELSE IF (FlawType .EQ. 4) THEN
      IF (a/radius .GT. 1.2D0) THEN
		 err = 1
	  END IF
!
!  through-width flaws in plates BS 7910:2005 M.3.3
   ELSE IF (FlawType .EQ. 3) THEN
      IF (a/B .GT. 0.6D0) THEN
	     err = 1
	  END IF
!
!  through-thickness flaws in plates BS 7910:2005 M.3.1
   ELSE IF (FlawType .EQ. 1) THEN
!
!  surface flaws in plates BS 7910:2005 M.3.2.2
   ELSE IF (FlawType .EQ. 2) THEN
      IF (a .LT. 0.0D0) THEN
	     err = 1
      ELSEIF (c .LT. 0.0D0) THEN
	     err = 1
      ELSEIF (a/B .GE. 1.D0) THEN
         err = 1
	  ELSEIF ((a/c) .GT. 2.D0) THEN
         err = 1
	  ELSEIF (c/W .GT. 0.4D0) THEN
	     err = 1
	  END IF
!
! Unknown flaw type
   ELSE
      STOP 'Incorrect input in Subroutine SIMConditions: Unknown flawtype'
   END IF
!
! End subroutine
   RETURN
   END

max_stars_repo_name: monicanicoreac/ProRail-ravelingen

max_issues_repo_name: monicanicoreac/ProRail-ravelingen

max_forks_repo_name: monicanicoreac/ProRail-ravelingen


text:
      SUBROUTINE DT_CHKINE(Pp1i,Ifp1,Pp2i,Ifp2,Pt1i,Ift1,Pt2i,Ift2,
     &                     Amch1,Amch1n,Amch2,Irej)

C***********************************************************************
C This subroutine replaces CORMOM.                                     *
C This version dated 05.01.95 is written by S. Roesler                 *
C***********************************************************************

      IMPLICIT NONE
      DOUBLE PRECISION Amch1 , Amch1n , Amch2 , amch22 , ech , p1 , p2 ,
     &                 pch , pp1 , Pp1i , pp2 , Pp2i , pt1 , Pt1i ,
     &                 pt2 , Pt2i , PYMASS , scale , TINY10 , xm1
      DOUBLE PRECISION xm2
      INTEGER i , Ifp1 , Ifp2 , Ift1 , Ift2 , Irej , irej1 , jmshl
      SAVE

      INCLUDE 'inc/dtflka'

      PARAMETER (TINY10=1.0D-10)

C flags for input different options
      INCLUDE 'inc/dtflg1'
C rejection counter
      INCLUDE 'inc/dtrejc'

      DIMENSION pp1(4) , pp2(4) , pt1(4) , pt2(4) , p1(4) , p2(4) ,
     &          Pp1i(4) , Pp2i(4) , Pt1i(4) , Pt2i(4)

      Irej = 0
      jmshl = IMShl

      scale = Amch1n/MAX(Amch1,TINY10)
      DO i = 1 , 4
         pp1(i) = Pp1i(i)
         pp2(i) = Pp2i(i)
         pt1(i) = Pt1i(i)
         pt2(i) = Pt2i(i)
         pp2(i) = pp2(i) + (1.0D0-scale)*pp1(i)
         pt2(i) = pt2(i) + (1.0D0-scale)*pt1(i)
         pp1(i) = scale*pp1(i)
         pt1(i) = scale*pt1(i)
      END DO
      IF ( (pp1(4).GE.0.0D0) .AND. (pp2(4).GE.0.0D0) .AND.
     &     (pt1(4).GE.0.0D0) .AND. (pt2(4).GE.0.0D0) ) THEN

         ech = pp2(4) + pt2(4)
         pch = SQRT((pp2(1)+pt2(1))**2+(pp2(2)+pt2(2))
     &         **2+(pp2(3)+pt2(3))**2)
         amch22 = (ech-pch)*(ech+pch)
         IF ( amch22.LT.0.0D0 ) THEN

            IF ( IOUlev(1).GT.0 .AND. LPRi.GT.4 ) WRITE (LOUt,'(1X,A)')
     &            'CHKINE: inconsistent treatment!'
            GOTO 200
         END IF

         Amch1 = Amch1n
         Amch2 = SQRT(amch22)

C put partons again on mass shell
 50      xm1 = 0.0D0
         xm2 = 0.0D0
         IF ( jmshl.EQ.1 ) THEN

            xm1 = PYMASS(Ifp1)
            xm2 = PYMASS(Ift1)

         END IF
         CALL DT_MASHEL(pp1,pt1,xm1,xm2,p1,p2,irej1)
         IF ( irej1.NE.0 ) THEN
            IF ( jmshl.EQ.0 ) GOTO 300
            jmshl = 0
            GOTO 50
         END IF
         jmshl = IMShl
         DO i = 1 , 4
            pp1(i) = p1(i)
            pt1(i) = p2(i)
         END DO
 100     xm1 = 0.0D0
         xm2 = 0.0D0
         IF ( jmshl.EQ.1 ) THEN

            xm1 = PYMASS(Ifp2)
            xm2 = PYMASS(Ift2)

         END IF
         CALL DT_MASHEL(pp2,pt2,xm1,xm2,p1,p2,irej1)
         IF ( irej1.NE.0 ) THEN
            IF ( jmshl.EQ.0 ) GOTO 300
            jmshl = 0
            GOTO 100
         END IF
         DO i = 1 , 4
            pp2(i) = p1(i)
            pt2(i) = p2(i)
         END DO
         DO i = 1 , 4
            Pp1i(i) = pp1(i)
            Pp2i(i) = pp2(i)
            Pt1i(i) = pt1(i)
            Pt2i(i) = pt2(i)
         END DO
         RETURN
      END IF

 200  IRChki(1) = IRChki(1) + 1
C*sr
C     GOTO 9999
      Irej = -1
      RETURN
C*
 300  IRChki(2) = IRChki(2) + 1
C*af
C9999 CONTINUE

      IF ( LPRi.GT.4 .AND. IOUlev(1).GT.0 ) WRITE (LOUt,*)
     &      'rejected 1 in CHKINE'
      Irej = 1
      END SUBROUTINE

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


text:
! Various subroutines for calculating quantities in Roche geometry
module roche
   use real_kind

contains

   ! Calculate the scaled (a = 1) x coordinate of the binary component Jstar.
   ! Internally, all calculations assume m1 > m2, q = m2/m1 < 1 and x1 > 0.
   ! This means that when the mass ratio inverts, the coordinate system changes.
   ! This can be handled entirely automatically, as long as the arguments passed
   ! in are consistent and all functions perform the required coordinate
   ! transformation.
   elemental function calc_scaled_x(m1, m2, Jstar)
   implicit none
   real(double), intent(in) :: m1, m2
   integer, intent(in) :: Jstar
   real(double) :: calc_scaled_x
   real(double) :: q
   integer :: n

   q = m2 / m1
   n = Jstar
   if (q > 1.0d0) then
      q = m1 / m2
      n = 3 - Jstar
   end if

   if (n == 1) then
      calc_scaled_x = q / (1.0d0 + q)
   else
      calc_scaled_x = -1.0d0 / (1.0d0 + q)
   end if
   end function calc_scaled_x

   elemental function calc_scaled_potential(m1, m2, x, y, z)
   implicit none
   real(double), intent(in) :: m1, m2, x, y, z
   real(double) :: calc_scaled_potential
   real(double) :: q, xx, qrat, yz2

   q = m2 / m1
   xx = x
   if (q > 1.0d0) then
      q = 1.0d0 / q
      xx = 1.0d0 - x
   end if

   qrat = 1.0d0 / (1.0d0 + q)
   yz2 = y**2 + z**2
   calc_scaled_potential = qrat * (1.0d0/sqrt((x - q*qrat)**2+yz2) + q / sqrt((x+qrat)**2+yz2)) + 0.5d0*(x**2+y**2)
   end function calc_scaled_potential

   elemental function calc_scaled_xl1(m1, m2)
   implicit none
   real(double), intent(in) :: m1, m2
   real(double) :: calc_scaled_xl1
   real(double) :: q,q1,q2,q3,q4,q5,q6,xl1
   integer :: ij

   q = m1 / m2
   if ( q < 1.0d0 ) q = 1.0d0/q
   q5 = 1.0d0 + q
   xl1 = (q/(3.0d0*q5))**(1.0d0 / 3.0d0)
   q3 = 3.0d0 + q
   q2 = q + q
   q4 = 3.0d0 + q2
   q1 = 2.0d0 + q
   q6 = 1.0d0/q5
  ! Newton-Raphson iteration for L1
   do ij = 1, 4
     xl1 = xl1 + &
          (q - xl1*(q2 - xl1*(q - xl1*(q3 - xl1*(q4 - xl1*q5)))))&
          /(q2 - xl1*(q2 - xl1*(3.0d0*q3 - xl1*(4.0d0*q4 - xl1*5.0d0*q5))))
   end do
   calc_scaled_xl1 = q6 - xl1
   end function calc_scaled_xl1

   elemental function calc_scaled_xl2(m1, m2)
   implicit none
   real(double), intent(in) :: m1, m2
   real(double) :: calc_scaled_xl2
   real(double) :: q,q1,q2,q3,q4,q5,q6,xl1,xl2
   integer :: ij

   q = m1 / m2
   if ( q < 1.0d0 ) q = 1.0d0/q
   q5 = 1.0d0 + q
   xl1 = (q/(3.0d0*q5))**(1.0d0 / 3.0d0)
   xl2 = 2.0d0 - q*xl1/q5
   q3 = 3.0d0 + q
   q2 = q + q
   q4 = 3.0d0 + q2
   q1 = 2.0d0 + q
   q6 = 1.0d0/q5
  ! Newton-Raphson iteration for L2
   do ij = 1, 4
     xl2 = xl2 + &
          (q - xl2*(q2 - xl2*(q1 - xl2*(q3 - xl2*(q4 - xl2*q5)))))&
          /(q2 - xl2*(2.0d0*q1 - xl2*(3.0d0*q3 - xl2*(4.0d0*q4 - xl2*5.0d0*q5))))
   end do
   calc_scaled_xl2 = q6 - xl2
   end function calc_scaled_xl2

   ! Calculate the correction factor for the cross section of the mass flow at level surface xx (dimensionless, L1=0, L2=1)
   elemental function stream_cross_section_primary(q, xx)
      use real_kind
      implicit none
      real(double), intent(in) :: q, xx
      real(double) :: stream_cross_section_primary
      real(double) :: a, b, c, d

      a = 4.50d-4 + (3.68d-4 + 3.15d-4 * q) * q
      if (q < 0.3d0) then
         b = 2.32d0 + (-1.31d0 + 2.87d0 * q) * q
      else
         b = 2.21d0 + 0.238d0 * q
      end if
      c = q * (1.22d0 + q*(-1.22d0 + q*0.564d0))
      d = 0.323d0 * q**0.230d0 - 7.39d-2

      stream_cross_section_primary = (c * xx + a)**b + d * xx
   end function stream_cross_section_primary

   ! Calculate the correction factor for the cross section of the mass flow at level surface xx (dimensionless, L1=0, L2=1)
   elemental function stream_cross_section_secondary(q, xx)
      use real_kind
      implicit none
      real(double), intent(in) :: q, xx
      real(double) :: stream_cross_section_secondary
      real(double) :: a, b, c, d

      a = 9.53d-4 * q**0.501d0 + 1.65d-4
      b = 3.25d-1 * q**0.326d0 + 2.03d0
      c = 6.33d-1 * q**0.271d0 - 7.18d-2
      d = 2.30d-1 * q**0.361d0 + 1.89d-2

      stream_cross_section_secondary = (c * xx + a)**b + d * xx
   end function stream_cross_section_secondary

   !> \brief Compute the Roche-lobe radius Rl/a
   !!
   !! \param x  q^(1/3)
   elemental function rlobe(x)
      use real_kind
      use constants

      implicit none
      real(double), intent(in) :: x
      real(double) :: rlobe,x2

      x2 = x**2
      rlobe = 0.49_dbl*x2 / (0.6_dbl*x2 + log(1.0_dbl + x))
   end function rlobe

end module

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse