! This file is a part of JuliaFEM.
! License is MIT: see https://github.com/JuliaFEM/UMAT.jl/blob/master/LICENSE
!***********************************************************************
! File drucker_prager_plasticity.f90
! Including:
!
!   subroutine UMAT
!   subroutine UMAT_Overstress
!   subroutine UMAT_Consistency
!

!***********************************************************************
!=======================================================================
! Subroutine UMAT
!
! ABAQUS User Material Subroutine, Calls real material subroutines
! depending on the material type
!-----------------------------------------------------------------------
!
!     Variables to be defined in all situtions
!     - - - - - - - - - - - - - - - - - - - - -
!      DDSDDE(NTENS,NTENS) Jacobian matrix of the constitutive model
!      STRESS(NTENS) This array is passed in as the stress tensor at the
!           beginning of the
!           increment and must be updated in this routine to be the stress
!           tensor at the end of the increment. If initial stresses are
!           specified using the *INITIAL CONDITIONS, TYPE=STRESS option,
!           this array will contain the initial stresses at the start of
!           the analysis. The size of this array depends on the value of
!           NTENS as defined below. In finite-strain problems the stress
!           tensor has already been rotated to account or rigid body motion
!           in the increment before UMAT is called, so that only the
!           corotational part of the stress integration should be done
!           in UMAT. The measure of stress used is "true" (Cauchy) stress.
!      STATEV(NSTATV) An array containing the solution-dependent state
!           variables. These are passed in as the values at the beginning of
!           the increment unless they are updated in user subroutines
!           USDFLD (``USDFLD,'' Section 24.2.37) or UEXPAN (``UEXPAN,''
!           Section 24.2.20), in which case the updated values are passed
!           in. In all cases STATEV must be returned as the values at
!           the end of the increment. The size of the array is that
!           defined on the *DEPVAR option associated with this material.
!           In finite-strain problems any vector-valued or tensor-valued
!           state variables must be rotated to account for rigid body
!           motion of the material, in addition to any update in the
!           values associated with constitutive behavior. The rotation
!           increment matrix, drott, is provided for this purpose.
!      SSE, SPD, SCD
!           Specific elastic strain energy, plastic dissipation, and
!           "creep" dissipation, respectively. These are passed in as
!           the values at the start of the increment and should be
!           updated to the corresponding specific energy values at the
!           end of the increment. They have no effect on the solution,
!           except that they are used for the energy output options.
!
!  The following only in a fully coupled temperature-displacement analysis:
!  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!      RPL  Volumetric heat generation per unit time at the end of the
!           increment caused by mechanical working of the material.
!      DDSDDT(NTENS) Variation of the stress increments with respect to
!           the temperature.
!      DRPLDE(NTENS) Variation of RPL with respect to the strain increments.
!      DRPLDT Variation of RPL with respect to the temperature.
!
!      Variable that can be updated
!      - - - - - - - - - - - - - - -
!       PNEWDT  Ratio of suggested new time increment to the time increment
!           being used (DTIME, see discussion later in this section). This
!           variable allows the user to provide input to the automatic
!           time incrementation algorithms in ABAQUS (if automatic time
!           incrementation is chosen)....
!
!  Variables passed in for information
!  - - - - - - - - - - - - - - - - - -
!        STRAN(NTENS)
!           An array containing the total strains at the beginning of
!           the increment. If the *EXPANSION option is used in the same
!           material definition, the strains passed into UMAT are the
!           mechanical strains only (that is, the thermal expansion strains
!           defined by the *EXPANSION option have been subtracted from the
!           total strains). These strains are available for output as the
!           "elastic" strains. In finite-strain problems the strain
!           components have been rotated to account for rigid body motion
!           in the increment before UMAT is called and are approximations
!           to logarithmic strain.
!        DSTRAN(NTENS)
!           Array of strain increments. If the *EXPANSION option is used
!           in the same material definition, these are the mechanical
!           strain increments (the total strain increments minus the thermal
!           strain increments).
!        TIME(1)
!           Value of step time at the beginning of the current increment.
!        TIME(2)
!           Value of total time at the beginning of the current increment.
!        DTIME
!           Time increment.
!        TEMP
!           Temperature at the start of the increment.
!        DTEMP
!           Increment of temperature.
!        PREDEF
!           Array of interpolated values of predefined field variables at
!           this point at the start of the increment, based on the values
!           read in at the nodes in the *FIELD option.
!        DPRED
!           Array of increments of predefined field variables.
!        CMNAME
!           Name given on the *MATERIAL option, left justified.
!        NDI
!           Number of direct stress components at this point.
!        NSHR
!           Number of engineering shear stress components at this point.
!        NTENS
!           Size of the stress or strain component array (NDI + NSHR).
!        NSTATV
!           Number of solution-dependent state variables that are associated
!           with this material type (as defined in the *DEPVAR option).
!        PROPS(NPROPS)
!           Array of material constants entered in the *USER MATERIAL option
!           for this material.
!        NPROPS
!           Number of material constants (the value given to the CONSTANTS
!           parameter on the *USER MATERIAL option).
!        COORDS
!           An array containing the coordinates of this point. These are
!           the current coordinates if geometric nonlinearity is accounted
!           for during the step (see ``Procedures: overview,'' Section 6.1.1);
!           otherwise, the array contains
!           the original coordinates of the point.
!        drott(3,3)
!           Rotation increment matrix. This matrix represents the increment
!           of rigid body rotation of the basis system in which the components
!           of stress (STRESS) and strain (STRAN) are stored. It is provided so
!           that vector- or tensor-valued state variables can be rotated
!           appropriately in this subroutine: stress and strain components are
!           already rotated by this amount before UMAT is called. This matrix
!           is passed in as a unit matrix for small-displacement analysis and
!           for large-displacement analysis if the basis system for the material
!           point rotates with the material (as in a shell element or when
!           the *ORIENTATION option is used).
!        CELENT
!           Characteristic element length, which is a typical length of a
!           line across an element for a first-order element; it is half of
!           the same typical length for a second-order element. For beams and
!           trusses it is a characteristic length along the element axis.
!           For membranes and shells it is a characteristic length in the
!           reference surface. For axisymmetric elements it is a characteristic
!           length in the  plane only.
!        DFGRD0(3,3)
!           Array containing the deformation gradient at the beginning of
!           the increment. See the discussion regarding the availability of
!           the deformation gradient
!           for various element types.
!        DFGRD1(3,3)
!           Array containing the deformation gradient at the end of the
!           increment. The components of this array are set to zero if
!           the NLGEOM parameter is not used on the *STEP option associated
!           with this increment. See the discussion regarding the availability
!           of the deformation gradient for various element types.
!        NOEL
!           Element number.
!        NPT
!           Integration point number.
!        LAYER
!           Layer number (for composite shells and layered solids).
!        KSPT
!           Section point number within the current layer.
!        KSTEP
!           Step number.
!        KINC
!           Increment number.
!
!-----------------------------------------------------------------------
! Originally coded 19/03/2003 KKO
! last modified
!=======================================================================
  subroutine umat(stress, statev, ddsdde, sse, spd, scd, rpl,           &
          ddsddt, drplde, drpldt, stran, dstran, time, dtime, temp,     &
          dtemp, predef, dpred, cmname, ndi, nshr, ntens, nstatv, props,&
          nprops, coords, drott, pnewdt, celent, dfgrd0, dfgrd1, noel,  &
          npt, layer, kspt, kstep, kinc)
!
	implicit none
!  include 'aba_param.inc'
! implicit real*8(a-h,o-z) !aba_param.inc
  INTEGER nprecd           !aba_param.inc
  parameter (nprecd=2)
!
  character*8 cmname
  real*8 :: stress, statev, ddsdde, ddsddt, drplde, stran, dstran,     &
      predef, dpred, props, coords, drott, dfgrd0, dfgrd1,sse,spd,scd, &
      rpl,drpldt,time,dtime,temp,dtemp,pnewdt,celent
	Integer ntens,nstatv,nprops,noel,npt,layer,kspt,kstep,kinc,ndi,      &
          nshr
  dimension stress(ntens), statev(nstatv), ddsdde(ntens, ntens),       &
      ddsddt(ntens), drplde(ntens), stran(ntens), dstran(ntens),       &
      predef(1), dpred(1), props(nprops), coords(3), drott(3, 3),      &
      dfgrd0(3, 3), dfgrd1(3, 3)
!- local variables -----------------------------------------------------
  Integer, Parameter :: DATF=6,MSGF=7,TESTFIL1=105,TESTFIL2=106
  Integer :: II,LENGTH,ISCAN,ELN,IPN,MTYP
  character*80 :: OPEN_File_Name1, OPEN_File_Name2
  Logical :: LOPEN,LNAME,From_Back
  real*8 :: stress_guess(ntens)

! Open files for printing results

  IF (KINC <= 1 .and. KSTEP == 1 .and. npt == 1) Then
    INQUIRE (UNIT=TESTFIL1, OPENED = LOPEN, NAMED=LNAME,NAME=OPEN_File_Name1)
    IF(.NOT. LOPEN)   Then
      INQUIRE (UNIT=6, OPENED = LOPEN, NAMED=LNAME,NAME=OPEN_File_Name1)
      From_Back = .True.
      ISCAN = SCAN(OPEN_File_Name1,'.',From_Back)
      If(ISCAN /= 0) Then
        OPEN_File_Name1(ISCAN:ISCAN+5) ='-1.txt'
        OPEN_File_Name2= OPEN_File_Name1
        OPEN_File_Name1(ISCAN:ISCAN+5) ='-2.txt'
      Else
        LENGTH = LEN_TRIM(OPEN_File_Name1)
        OPEN_File_Name1(LENGTH+1:LENGTH+6) = '-1.txt'
        OPEN_File_Name2= OPEN_File_Name1
        OPEN_File_Name2(LENGTH+1: LENGTH+6) = '-1.txt'
      Endif
      Open(unit=TESTFIL1,status='UNKNOWN',file=OPEN_File_Name1)
      Open(unit=TESTFIL2,status='UNKNOWN',file=OPEN_File_Name2)
    Endif

   ! DO II=1,100
   !   INQUIRE (UNIT=II, OPENED = LOPEN, NAMED=LNAME,NAME=OPEN_File_Name1)
   !   IF(LOPEN) THEN
   !     Write(MSGF,*) 'II,LOPEN',II,LOPEN
   !     Write(MSGF,*) 'LNAME',LNAME
   !     Write(MSGF,*) 'FNAME',OPEN_File_Name1
   !   ENDIF
   ! ENDDO

  Endif

  ELN  =  INT(PROPS(15))
  IPN  =  INT(PROPS(16))
  MTYP =  INT(PROPS(17))

! Update stresses for printing:
  IF (KINC >= 1 .and. NOEL == ELN .and. NPT == IPN ) &
    stress_guess = stress

  If(MTYP == 1) Then
    CALL UMAT_Overstress(stress, statev, ddsdde, sse, spd, scd, stran,  &
          dstran, time, dtime, cmname, ndi, nshr, ntens, nstatv, props, &
          nprops, coords, drott, pnewdt, celent, noel,                  &
          npt, layer, kspt, kstep, kinc)

  ElseIf(MTYP == 2 .or. MTYP == 4)Then
    CALL UMAT_Consistency(stress, statev, ddsdde, sse, spd, scd, stran, &
          dstran, time, dtime, cmname, ndi, nshr, ntens, nstatv, props, &
          nprops, coords, drott, pnewdt, celent, noel,                  &
          npt, layer, kspt, kstep, kinc)
  Endif

! Updete stresses for printing:
  IF (KINC >= 1 .and. NOEL == ELN .and. NPT == IPN ) &
    stress_guess = stress_guess + matmul(ddsdde,dstran)

  IF (KINC >= 1 .and. NOEL == ELN .and. NPT == IPN ) &
    WRITE(TESTFIL1,'(1P,21E13.4)') (stran(II)+dstran(II),II=1,6),&
                (stress(II),II=1,6),statev(1),(stress_guess(II),II=1,6)

  RETURN
	end subroutine umat
!
!***********************************************************************
!=======================================================================
! Subroutine UMAT_Overstress
!
! ABAQUS User Material Subroutine
!-----------------------------------------------------------------------
!
! input parameters
! ----------------
!  PROPS(1)  - Sigma_t
!  PROPS(2)  - Sigma_c
!  PROPS(3)  - K_Zero
!  PROPS(4)  - C_const
!  PROPS(5)  - N_Power
!  PROPS(6)  - Eta
!  PROPS(7)  - A_Const
!  PROPS(8)  - B_Const
!  PROPS(9)  - E
!  PROPS(10) - NU
!  PROPS(11) - DT (If real DTIME = 0, gien as negative)
!  PROPS(12) - 1.0       ! Print All to DAT
!  PROPS(13) - FTOLER (tolerance for local iteration
!  PROPS(14) - MAXITER
!  PROPS(15) - Element number for printing
!  PROPS(16) - Integration point number for printing
!  PROPS(17) - Plasticity type
!         1  - Perzyna Drucker-Prager (overstress)
!         2  - Drucker-Prager Consistency model for Perzyna's overstress
!              modell
!  PROPS(18) - Hardening type
!         1  - Linear hardening
!         2  - Exponential hardening
!
! input/output parameters
! ----------------
!
!  STATEV(1) - Equivalent plastic strain
!  stress    - Stress vector, to be updated
!  ddsdde    - Consistent tangent, to be updated
!
!-----------------------------------------------------------------------
! Originally coded 19/03/2003 KKO
! last modified
!=======================================================================
!-----------------------------------------------------------------------
subroutine UMAT_Overstress(stress, statev, ddsdde, sse, spd, scd, stran,    &
          dstran, time, dtime, cmname, ndi, nshr, ntens, nstatv, props, &
          nprops, coords, drott, pnewdt, celent, noel,                  &
          npt, layer, kspt, kstep, kinc)
  implicit none
!  include 'aba_param.inc'
! implicit real*8(a-h,o-z) !aba_param.inc
  integer nprecd           !aba_param.inc
  parameter (nprecd=2)
  Integer, Parameter :: dp = Selected_real_kind(12,50)
  character*8 cmname
  Real(dp) :: stress, statev, ddsdde, ddsddt, drplde, stran, dstran,   &
      predef, dpred, props, coords, drott, dfgrd0, dfgrd1,sse,spd,scd, &
      rpl,drpldt,time,dtime,temp,dtemp,pnewdt,celent
	Integer ntens,nstatv,nprops,noel,npt,layer,kspt,kstep,kinc,ndi,      &
	    nshr
  dimension stress(ntens), statev(nstatv), ddsdde(ntens, ntens),       &
      ddsddt(ntens), drplde(ntens), stran(ntens), dstran(ntens),       &
      predef(1), dpred(1), props(nprops), coords(3), drott(3, 3),      &
      dfgrd0(3, 3), dfgrd1(3, 3)
!- local variables -----------------------------------------------------
  Real(dp) :: EMOD,ENU,EBULK3,EG2,EG,ELAM,ONE,TWO,THREE
  Real(dp) :: Stress_Old(ntens),YIELDF,F_Tilde,Dev_Stress(ntens),SM, &
              ALPHA,A_const,B_const, DF_DSigma(ntens),               &
              DDF_DDSigma(ntens,ntens),Lambda,Eta,H(ntens,ntens),    &
              N_Power,CINV(ntens,ntens),Resid,HN(ntens)
  Real(dp) :: TOLER,DT,DLambda_DF,DPhi_DF,DKappa_DLambda,DF_DKappa,  &
              LambdaDot,DF_DLambda,DDF_DSigmaDLambda(ntens),         &
              DF_DLambdaDot,DDF_DSigmaDLambdaDot(ntens),             &
              Statev_Old(nstatv),DPHIDF,AA
  Character*80 :: string,form,CFLAG,DumCHR
	Integer K1,K2,II,JJ,ITER,MAXITER,ICONVE, Hardening_type,PRNTYP,ELN,IPN
  Integer, Parameter :: DATF=6,LOGF=5,MSGF=7,TESTFIL1=105,TESTFIL2=106
  Parameter(ONE=1.0_dp,TWO=2.0_dp,THREE=3.0_dp)

  Real(dp) :: SLEN,BEETA,SIGMA_T,SIGMA_C,J2INVARIANT,HARDENING_K,DUM

  If (NDI.NE.3) then
    Write(6,1)
 1  Format(//,30X,'***ERROR - THIS UMAT MAY ONLY BE USED FOR ',&
               'ELEMENTS WITH THREE DIRECT STRESS COMPONENTS')
  Endif

  N_Power = PROPS(5)
  Eta     = PROPS(6)
  A_Const = PROPS(7)
  B_Const = PROPS(8)
  PRNTYP  = INT(PROPS(12))
  TOLER   = PROPS(13)
  ELN     = INT(PROPS(15))
  IPN     = INT(PROPS(16))
  !MTYP    = INT(PROPS(17))
  MAXITER = INT(PROPS(14))
  IF(MAXITER < 2) MAXITER = 2
  DT      = DTIME
  IF(PROPS(11) < 0.0_dp ) DT = -PROPS(11)
  CFLAG = 'Drucker-Prager_Overstress_A'

! Check Hardening type
  If(INT(PROPS(18)) == 1) Then
    Hardening_type = 1
  ElseIf(INT(PROPS(18)) == 2) Then
    Hardening_type = 2
  Else
    Write(DumCHR,*) PROPS(18)
    string = ' Unknown Hardening type: No = '//Trim(Adjustl(DumCHR))
    form = '(/,'//''' *** ERROR in UMAT_Vers01 *** '''//',/,A)'
    Call WriteOut(string,form,4,KSTEP,KINC,NOEL,NPT)
    Call xit
  Endif

!     ELASTIC PROPERTIES

  EMOD=PROPS(9)
  ENU=PROPS(10)
  IF(ENU.GT.0.4999.AND.ENU.LT.0.5001) ENU=0.499
  EBULK3=EMOD/(ONE-TWO*ENU)
  EG2=EMOD/(ONE+ENU)
  EG=EG2/TWO
  ELAM=(EBULK3-EG2)/THREE

!     ELASTIC STIFFNESS

  Do K1=1,NTENS
    Do K2=1,NTENS
       DDSDDE(K2,K1)=0.0
    Enddo
  Enddo

  Do K1=1,NDI
    Do K2=1,NDI
       DDSDDE(K2,K1)=ELAM
    Enddo
    DDSDDE(K1,K1)=EG2+ELAM
  Enddo

  Do K1=NDI+1,NTENS
    DDSDDE(K1,K1)=EG
  Enddo

  Stress_Old = stress
  stress = stress + MATMUL(DDSDDE,dstran)

! Check if Yielding or not?
  F_Tilde = YIELDF(NTENS, STRESS, PROPS, STATEV, CFLAG,Hardening_type)

  If(F_Tilde <= 0.0_dp) Then
!   tee elastic temput
    RETURN
  Endif

! Print if PRNTYP == 1

  IF(PRNTYP == 1 .and. noel == ELN) Then
    string = ' Subroutine UMAT_Overstress, material type: '//Trim(cmname)
    form = '(A)'
    Call WriteOut(string,form,4,KSTEP,KINC,NOEL,NPT)
    Write(DumCHR,'(1P,E13.4)') F_Tilde
    string = ' '
    form = '('''//' F_Tilde = '//Trim(Adjustl(DumCHR))//''','//' A)'
    Call WriteOut(string,form,4,0,KINC,NOEL,NPT)
  Endif
  IF(PRNTYP /= 0  .and. NPT == IPN .and. noel == ELN) &
     Write(MSGF,'(A,1P,4E12.3)') '** DT, DTIME',DT, DTIME
  IF(PRNTYP /=0  .and. NPT == IPN .and. noel == ELN) &
    WRITE(MSGF,'(2A)') &
      'ITER,F           Lambda      Kappa,      K,',&
             '          DLambda_DF, slen        R_norm'

! Local Iteration ( Wang et al )

  Lambda  = 0.0_dp
  ICONVE  = 0
  Statev_Old = Statev
  Call invert( DDSDDE, CINV, ntens )

! > > > > > > > > > > > > > > > > > > > >

  Do ITER = 1, MAXITER

    F_Tilde = YIELDF(NTENS,STRESS,PROPS,STATEV,CFLAG,Hardening_type)
    Resid = (F_Tilde / A_Const+B_Const)**N_Power -Lambda/Eta/DT
    If(Abs(Resid) < TOLER) Then
      ICONVE = 1
      EXIT
    Endif

    Call DEVSTRESS(NTENS, STRESS,SM,DEV_STRESS)
    Call DF_Derivatives(NTENS, CFLAG,Hardening_type,STRESS,DEV_STRESS,&
                          PROPS,                                      &
                          STATEV,LambdaDot,DF_DLambda,DF_DLambdaDot,  &
                          DDF_DDSigma,DDF_DSigmaDLambda,              &
                          DDF_DSigmaDLambdaDot,Alpha)
    H = CINV + Lambda*DDF_DDSigma
    Call invert( H, H, ntens )

    Call DFDSIGMA(NTENS,STRESS,DEV_STRESS,PROPS,STATEV,CFLAG,         &
                Hardening_type,DF_DSIGMA)
    HN = MATMUL(H, Lambda*DDF_DSigmaDLambda + DF_DSigma)
    DPhi_DF= DPHIDF(F_Tilde,PROPS,STRESS,NTENS, STATEV, CFLAG,        &
                    Hardening_type)
    AA = DPhi_DF*Dot_Product(DF_DSigma,HN) +1.0_dp/Eta/DT             &
        -DPhi_DF*DF_DLambda
    Lambda = Lambda + Resid/AA
    Stress = Stress_Old + Matmul(DDSDDE,dstran-Lambda*DF_DSigma)
    Statev(1) = Statev_Old(1) + (1+Alpha)*Lambda

!   Printing:

    IF(PRNTYP /= 0  .and. NPT == IPN .and. noel == ELN ) Then
      SLEN    = SQRT(2.0_dp * J2INVARIANT(NTENS, DEV_STRESS) )
      WRITE(MSGF,'(1P,I2,1x,9E12.4)') ITER,F_Tilde,Lambda,Statev(1),&
             HARDENING_K(PROPS, STATEV, Hardening_type),1.0/DF_DLambda,&
             SLEN,Resid
    EndIf

  Enddo
! < < < < < < < < < < < < < < < < < < < <

!   Printing:

    IF(PRNTYP /= 0  .and. NPT == IPN .and. noel == ELN) Then
      SLEN    = SQRT(2.0_dp * J2INVARIANT(NTENS, DEV_STRESS) )
      WRITE(MSGF,'(1P,I2,1x,9E12.4)') ITER,F_Tilde,Lambda,Statev(1),&
             HARDENING_K(PROPS, STATEV, Hardening_type),1.0/DF_DLambda,&
             SLEN,Resid
    EndIf
    IF(PRNTYP == 2  .and. NPT == IPN .and. noel == ELN) Then
      WRITE(MSGF,'(/,A)') 'ITER,F,Lambda,Kappa,K,DLambda_DF,0.0, &
                        Dot_Product(Resid,Resid),DPhi_DF,DKappa_DLambda'
      WRITE(MSGF,'(1P,I2,1x,13E12.4)') ITER,F_Tilde,Lambda,Statev(1), &
                           HARDENING_K(PROPS, STATEV, Hardening_type),&
                           1.0/DF_DLambda,slen,Resid,DPhi_DF,1.0+alpha
      WRITE(MSGF,'(/,A,2I3)') ' STRESS, KINC,ITER',KINC,ITER
      WRITE(MSGF,'(1P,6E13.4)') (STRESS(ii),ii=1,6)
      WRITE(MSGF,'(/,A,2I3)') ' Dev_STRESS, KINC,ITER',KINC,ITER
      WRITE(MSGF,'(1P,6E13.4)') (Dev_STRESS(ii),ii=1,6)
      WRITE(MSGF,'(/,A,2I3)') ' DF_DSigma, KINC,ITER',KINC,ITER
      WRITE(MSGF,'(1P,6E13.4)') (DF_DSigma(ii),ii=1,6)
      WRITE(MSGF,'(/,A,2I3)') ' DDF_DDSigma, KINC,ITER',KINC,ITER
      WRITE(MSGF,'(1P,6E13.4)') ((DDF_DDSigma(ii,jj),ii=1,6),jj=1,6)
    ENDIF

    IF(PRNTYP /= 0.and. NPT==IPN .and. noel==ELN .and. ICONVE == 1) Then
      SIGMA_T = PROPS(1)
      SIGMA_C = PROPS(2)
      Call DEVSTRESS(NTENS, STRESS,SM,DEV_STRESS)
      ALPHA   = (SIGMA_C - SIGMA_T) / (SIGMA_C + SIGMA_T)
      BEETA   = (1.0_dp + ALPHA) * SIGMA_T
      SLEN    = SQRT(2.0_dp * J2INVARIANT(NTENS, DEV_STRESS) )

      WRITE(TESTFIL2,'(1P,6E13.4)') &
                    3.0*SM,SQRT(1.5_dp)*SLEN,BEETA-ALPHA*3*SM,&
                    HARDENING_K(PROPS, STATEV, Hardening_type),Statev(1)
    ENDIF

  IF (ICONVE == 1) Then

!   "Algorithmic Tangent"

    If(ITER > 1 ) Then
      Stress_Old = MATMUL(H,DF_DSigma) ! Stress_Old used as dummy vector
      AA = DPhi_DF/AA                  ! Now: HN = Stress_Old
      DO II=1,ntens
        DO JJ=1,ntens
          DDSDDE(JJ,II) = H(JJ,II) - AA*HN(JJ)*Stress_Old(II)
        ENDDO
      ENDDO
    Endif

!   Suggestion for new increment length
    DUM = ITER*1.0/(1.0*MAXITER)
    IF(DUM>0.8) PNEWDT = 0.8D0

  Else
    write(DumCHR,'(A,1P,E13.4,A,E13.4)') ' ABS(F)',ABS(F_Tilde),' > ',TOLER
    string =' Material type: '//Trim(cmname)//' '//Trim(Adjustl(DumCHR))
    form = '(/,'//''' *** ERROR in UMAT_Vers01 *** '''//',/,A)'
    Call WriteOut(string,form,4,KSTEP,KINC,NOEL,NPT)
    !Call xit
    PNEWDT = 0.75D0
  Endif

End subroutine UMAT_Overstress
!***********************************************************************
!=======================================================================
! Subroutine UMAT_Consistency
!
! ABAQUS User Material Subroutine For Drucker Prager (Perzyna Consistency
! model)
!-----------------------------------------------------------------------
!
! input parameters
! ----------------
!  PROPS(1)  - Sigma_t
!  PROPS(2)  - Sigma_c
!  PROPS(3)  - K_Zero
!  PROPS(4)  - C_const
!  PROPS(5)  - N_Power
!  PROPS(6)  - Eta
!  PROPS(7)  - A_Const
!  PROPS(8)  - B_Const
!  PROPS(9)  - E
!  PROPS(10) - NU
!  PROPS(11) - DT (If real DTIME = 0, gien as negative)
!  PROPS(12) - 1.0       ! Print All to DAT
!  PROPS(13) - FTOLER
!  PROPS(14) - MAXITER
!  PROPS(15) - Element number for printing
!  PROPS(16) - Integration point number for printing
!  PROPS(17) - Plasticity type
!         1  - Perzyna Drucker-Prager (overstress)
!         2  - Drucker-Prager Consistency model for Perzyna's overstress
!              modell
!  PROPS(18) - Hardening type
!         1  - Linear hardening
!         2  - Exponential hardening
!
! input/output parameters
! ----------------
!
!  STATEV(1) - Equivalent plastic strain
!  stress    - Stress vector, to be updated
!  ddsdde    - Consistent tangent, to be updated
!
!-----------------------------------------------------------------------
! Originally coded 19/03/2003 KKO
! last modified
!=======================================================================
!-----------------------------------------------------------------------
subroutine UMAT_Consistency(stress, statev, ddsdde, sse, spd, scd, stran,&
          dstran, time, dtime, cmname, ndi, nshr, ntens, nstatv, props,  &
          nprops, coords, drott, pnewdt, celent, noel,                   &
          npt, layer, kspt, kstep, kinc)
  implicit none
!  include 'aba_param.inc'
! implicit real*8(a-h,o-z) !aba_param.inc
  integer nprecd           !aba_param.inc
  parameter (nprecd=2)
  Integer, Parameter :: dp = Selected_real_kind(12,50)
  character*8 cmname
  Real(dp) :: stress, statev, ddsdde, ddsddt, drplde, stran, dstran,   &
      predef, dpred, props, coords, drott, dfgrd0, dfgrd1,sse,spd,scd, &
      rpl,drpldt,time,dtime,temp,dtemp,pnewdt,celent
	Integer ntens,nstatv,nprops,noel,npt,layer,kspt,kstep,kinc,ndi,      &
	    nshr
  dimension stress(ntens), statev(nstatv), ddsdde(ntens, ntens),       &
      ddsddt(ntens), drplde(ntens), stran(ntens), dstran(ntens),       &
      predef(1), dpred(1), props(nprops), coords(3), drott(3, 3),      &
      dfgrd0(3, 3), dfgrd1(3, 3)
!- local variables -----------------------------------------------------
  Real(dp) :: EMOD,ENU,EBULK3,EG2,EG,ELAM,ONE,TWO,THREE
  Real(dp) :: Stress_Old(ntens),YIELDF,F_Tilde,Dev_Stress(ntens),SM, &
              ALPHA,A_const,B_const, DF_DSigma(ntens),               &
              DDF_DDSigma(ntens,ntens),Lambda,Eta,H(ntens,ntens),    &
              N_Power,CINV(ntens,ntens),Resid,HN(ntens)
  Real(dp) :: TOLER,DT,DLambda_DF,DPhi_DF,DKappa_DLambda,DF_DKappa,  &
              LambdaDot,DF_DLambda,DDF_DSigmaDLambda(ntens),         &
              DF_DLambdaDot,DDF_DSigmaDLambdaDot(ntens),             &
              Statev_Old(nstatv),DPHIDF,AA,Lambda_Dot
  Character*80 :: string,form,CFLAG,DumCHR
	Integer K1,K2,II,JJ,ITER,MAXITER,ICONVE, Hardening_type,PRNTYP,ELN,IPN
  Integer, Parameter :: DATF=6,LOGF=5,MSGF=7,TESTFIL1=105,TESTFIL2=106
  Parameter(ONE=1.0_dp,TWO=2.0_dp,THREE=3.0_dp)

  Real(dp) :: SLEN,BEETA,SIGMA_T,SIGMA_C,J2INVARIANT,HARDENING_K,DUM

  If (NDI.NE.3) then
    Write(6,1)
 1  Format(//,30X,'***ERROR - THIS UMAT MAY ONLY BE USED FOR ',&
               'ELEMENTS WITH THREE DIRECT STRESS COMPONENTS')
  Endif

  N_Power = PROPS(5)
  Eta     = PROPS(6)
  A_Const = PROPS(7)
  B_Const = PROPS(8)
  PRNTYP  = INT(PROPS(12))
  TOLER   = PROPS(13)
  MAXITER = INT(PROPS(14))
  ELN     = INT(PROPS(15))
  IPN     = INT(PROPS(16))
  IF(MAXITER < 2) MAXITER = 2
  DT      = DTIME
  IF(PROPS(11) < 0.0_dp ) DT = -PROPS(11)
  CFLAG = 'Drucker-Prager_Consistency_A'

! Check Hardening type

  If(INT(PROPS(18)) == 1) Then
    Hardening_type = 1
  ElseIf(INT(PROPS(18)) == 2) Then
    Hardening_type = 2
  Else
    Write(DumCHR,*) PROPS(18)
    string = ' Unknown Hardening type: No = '//Trim(Adjustl(DumCHR))
    form = '(/,'//''' *** ERROR in UMAT_Vers01 *** '''//',/,A)'
    Call WriteOut(string,form,4,KSTEP,KINC,NOEL,NPT)
    Call xit
  Endif

!     ELASTIC PROPERTIES

  EMOD=PROPS(9)
  ENU=PROPS(10)
  IF(ENU.GT.0.4999.AND.ENU.LT.0.5001) ENU=0.499
  EBULK3=EMOD/(ONE-TWO*ENU)
  EG2=EMOD/(ONE+ENU)
  EG=EG2/TWO
  ELAM=(EBULK3-EG2)/THREE

!     ELASTIC STIFFNESS

  Do K1=1,NTENS
    Do K2=1,NTENS
       DDSDDE(K2,K1)=0.0
    Enddo
  Enddo

  Do K1=1,NDI
    Do K2=1,NDI
       DDSDDE(K2,K1)=ELAM
    Enddo
    DDSDDE(K1,K1)=EG2+ELAM
  Enddo

  Do K1=NDI+1,NTENS
    DDSDDE(K1,K1)=EG
  Enddo

  Stress_Old = stress
  stress = stress + MATMUL(DDSDDE,dstran)

! Check if Yielding or not?
  F_Tilde = YIELDF(NTENS, STRESS, PROPS, STATEV, CFLAG,Hardening_type)

  If(F_Tilde <= 0.0_dp) Then
!   tee elastic temput
    RETURN
  Endif

  Lambda    = 0.0_dp
  Lambda_Dot= 0.0_dp
  Statev_Old = Statev
  ICONVE = 0
  Call invert( DDSDDE, CINV, ntens )

! Print if PRNTYP == 1

  IF(PRNTYP == 1 .and. noel == ELN) Then
    string = ' Subroutine UMAT_Consistency, material type: '//Trim(cmname)
    form = '(A)'
    Call WriteOut(string,form,4,KSTEP,KINC,NOEL,NPT)
    Write(DumCHR,'(1P,E13.4)') F_Tilde
    string = ' '
    form = '('''//' F_Tilde = '//Trim(Adjustl(DumCHR))//''','//' A)'
    Call WriteOut(string,form,4,0,KINC,NOEL,NPT)
  Endif
  IF(PRNTYP /= 0  .and. NPT == IPN .and. noel == ELN) &
     Write(MSGF,'(A,1P,4E12.3)') '** DT, DTIME',DT, DTIME
  IF(PRNTYP /=0  .and. NPT == IPN .and. noel == ELN) &
    WRITE(MSGF,'(2A)') &
      'ITER,F           Lambda      Kappa,      K,',&
             '          DLambda_DF, slen        R_norm'

! > > > > > > > > > > > > > > > > > > > >

! Local Iteration ( Wang et al )

  Do ITER = 1, MAXITER

    Lambda_Dot = Lambda/DT
    F_Tilde = YIELDF(NTENS,STRESS,PROPS,STATEV,CFLAG,Hardening_type) &
              -  A_Const*(B_Const+(Lambda_Dot/Eta)**(1.0_dp/N_Power))
    Resid = F_Tilde
    If(Abs(Resid) < TOLER) Then
      ICONVE = 1
      EXIT
    Endif

    Call DEVSTRESS(NTENS, STRESS,SM,DEV_STRESS)
    Call DF_Derivatives(NTENS, CFLAG,Hardening_type,STRESS,DEV_STRESS,&
                          PROPS,                                      &
                          STATEV,LambdaDot,DF_DLambda,DF_DLambdaDot,  &
                          DDF_DDSigma,DDF_DSigmaDLambda,              &
                          DDF_DSigmaDLambdaDot,Alpha)
    H = CINV + Lambda*DDF_DDSigma
    Call invert( H, H, ntens )

    Call DFDSIGMA(NTENS,STRESS,DEV_STRESS,PROPS,STATEV,CFLAG,         &
                Hardening_type,DF_DSIGMA)
    HN = MATMUL(H, Lambda*DDF_DSigmaDLambda + DF_DSigma               &
                + Lambda/DT*DDF_DSigmaDLambdaDot)
    AA = Dot_Product(DF_DSigma,HN) - DF_DLambda - DF_DLambdaDot/DT
    Lambda = Lambda + Resid/AA
    Stress = Stress_Old + Matmul(DDSDDE,dstran-Lambda*DF_DSigma)
    Statev(1) = Statev_Old(1) + (1+Alpha)*Lambda

!   Printing:

    IF(PRNTYP /= 0  .and. NPT == IPN .and. noel == ELN ) Then
      SLEN    = SQRT(2.0_dp * J2INVARIANT(NTENS, DEV_STRESS) )
      WRITE(MSGF,'(1P,I2,1x,9E12.4)') ITER,F_Tilde,Lambda,Statev(1),&
             HARDENING_K(PROPS, STATEV, Hardening_type),1.0/DF_DLambda,&
             SLEN,Resid
    EndIf

  Enddo

! < < < < < < < < < < < < < < < < < < < <

!   Printing:

    IF(PRNTYP /= 0  .and. NPT == IPN .and. noel == ELN) Then
      SLEN    = SQRT(2.0_dp * J2INVARIANT(NTENS, DEV_STRESS) )
      WRITE(MSGF,'(1P,I2,1x,9E12.4)') ITER,F_Tilde,Lambda,Statev(1),&
             HARDENING_K(PROPS, STATEV, Hardening_type),1.0/DF_DLambda,&
             SLEN,Resid
    EndIf
    IF(PRNTYP == 2  .and. NPT == IPN .and. noel == ELN) Then
      WRITE(MSGF,'(/,A)') 'ITER,F,Lambda,Kappa,K,DLambda_DF,0.0, &
                        Dot_Product(Resid,Resid),DF_DLambdaDot,DKappa_DLambda'
      WRITE(MSGF,'(1P,I2,1x,13E12.4)') ITER,F_Tilde,Lambda,Statev(1), &
                           HARDENING_K(PROPS, STATEV, Hardening_type),&
                           1.0/DF_DLambda,slen,Resid,DF_DLambdaDot,1.0+alpha
      WRITE(MSGF,'(/,A,2I3)') ' STRESS, KINC,ITER',KINC,ITER
      WRITE(MSGF,'(1P,6E13.4)') (STRESS(ii),ii=1,6)
      WRITE(MSGF,'(/,A,2I3)') ' Dev_STRESS, KINC,ITER',KINC,ITER
      WRITE(MSGF,'(1P,6E13.4)') (Dev_STRESS(ii),ii=1,6)
      WRITE(MSGF,'(/,A,2I3)') ' DF_DSigma, KINC,ITER',KINC,ITER
      WRITE(MSGF,'(1P,6E13.4)') (DF_DSigma(ii),ii=1,6)
      WRITE(MSGF,'(/,A,2I3)') ' DDF_DDSigma, KINC,ITER',KINC,ITER
      WRITE(MSGF,'(1P,6E13.4)') ((DDF_DDSigma(ii,jj),ii=1,6),jj=1,6)
    ENDIF

    IF(PRNTYP /= 0.and. NPT==IPN .and. noel==ELN .and. ICONVE == 1) Then
      SIGMA_T = PROPS(1)
      SIGMA_C = PROPS(2)
      Call DEVSTRESS(NTENS, STRESS,SM,DEV_STRESS)
      ALPHA   = (SIGMA_C - SIGMA_T) / (SIGMA_C + SIGMA_T)
      BEETA   = (1.0_dp + ALPHA) * SIGMA_T
      SLEN    = SQRT(2.0_dp * J2INVARIANT(NTENS, DEV_STRESS) )

      WRITE(TESTFIL2,'(1P,6E13.4)') &
                    3.0*SM,SQRT(1.5_dp)*SLEN,BEETA-ALPHA*3*SM,&
                    HARDENING_K(PROPS, STATEV, Hardening_type),Statev(1)
    ENDIF

  IF (ICONVE == 1) Then

!   "Algorithmic Tangent"

    If(ITER > 1 ) Then
      Stress_Old = MATMUL(H,DF_DSigma) ! Stress_Old used as dummy vector
      AA = 1.0_dp/AA                   ! Now: HN = Stress_Old
      DO II=1,ntens
        DO JJ=1,ntens
          DDSDDE(JJ,II) = H(JJ,II) - AA*HN(JJ)*Stress_Old(II)
        ENDDO
      ENDDO
    Endif

!   Suggestion for new increment length
    DUM = ITER*1.0/(1.0*MAXITER)
    IF(DUM>0.8) PNEWDT = 0.8D0

  Else
    write(DumCHR,'(A,1P,E13.4,A,E13.4)') ' ABS(F)',ABS(F_Tilde),' > ',TOLER
    string =' Material type: '//Trim(cmname)//' '//Trim(Adjustl(DumCHR))
    form = '(/,'//''' *** ERROR in UMAT_Vers01 *** '''//',/,A)'
    Call WriteOut(string,form,4,KSTEP,KINC,NOEL,NPT)
    !Call xit
    PNEWDT = 0.75D0
  Endif

End subroutine UMAT_Consistency


text:
program string
implicit none
character(50) :: str
    str = "Hi, I'm using string"
    str = 'Use "Double" "quote"'
    str = "'random'' ab''cd""fg''h'ij""klmno'pqr's'''"
    str = 'random'' a"b''cd""fg""hij"kl"mno"p"qr"s'''
    str = """equal''"
    str = '"Double" ''Quote'''
    str = "'Single' ""Quote"""

    print *, "Hi, I'm using string"
    print *, 'This is a "sting"'
    print *, "'random'' ab''cd""fg''h'ij""klmno'pqr's'''"
end program

max_stars_repo_name: Thirumalai-Shaktivel/lfortran

max_issues_repo_name: Thirumalai-Shaktivel/lfortran

max_forks_repo_name: Thirumalai-Shaktivel/lfortran


text:
       PROGRAM  QSAR
C#######################################################################
C  NAME OF ROUTINE
C      QSAR ( Quadrilateral Segment Averaging Routine )
C             Quadrilateral means rectangle here. SAR part is a misnomer.
C
C  PURPOSE
C      THIS IS THE STANDARD MAIN PROGRAM USED FOR TAE/VICAR PROGRAMS.
C      THIS MODULE CALLS SUBROUTINE MAIN44 TO ENTER INTO THE BODY OF THE
C      PROGRAM.
C      Program QSAR is used to add or subtract values from the data numbers in
C      rectangular sections of an image.
C  PREPARED FOR USE ON MIPL SYSTEM BY
C      STEVE POHORSKY   INFORMATICS GENERAL CORPORATION    OCTOBER 1983
C  FOR
C      EARTH RESOURCES APPLICATIONS
C
C  ORIGINAL QSAR PROGRAM BY
C      HOWARD FRIEDEN with modifications by JOEL MOSHER
C
C  REVISION HISTORY
C
C  CALLING SEQUENCE (TAE COMMAND LINE)
C      The following command line formats show the major allowable forms:
C
C      QSAR INP=a OUT=b SIZE=(sl,ss,nl,ns) optional parameters
C      QSAR INP=a OUT=b SL=sl SS=ss NL=nl NS=ns optional parameters
C      QSAR a b (sl,ss,nl,ns) optional parameters
C      QSAR a b optional parameters
C
C       Here 'a' represents the input image file name,
C       'b' represents the output image file name.
C
C  INPUT PARAMETERS (listed by keyword)
C      INP    - Input file name.
C      OUT    - Output file name.
C      SIZE   - Standard Vicar size field:  (SL,SS,NL,NS)
C               SL = Starting line number.
C               SS = Starting sample number.
C               NL = Number of lines.
C               NS = Number of samples.
C      AREA   - Specifies the rectangles and the values to be added to data
C               numbers.
C  OUTPUT PARAMETERS
C      The output image produced is written to the output file.
C  PROGRAM LIMITATIONS
C      1. The input and output images must be byte or halfword data.
C      2. The maximum number of pixels per line is 64000.
C      3. The maximum number of rectangles that can be specified is 120,
C         as of this writing.
C  SUBROUTINES CALLED
C      MAIN44
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

        INCLUDE 'VICMAIN_FOR'
      SUBROUTINE MAIN44
C#######################################################################
C  NAME OF ROUTINE
C     MAIN44 (name for top level subroutine by VICAR convention)
C
C  PURPOSE
C      MAIN44 is used to add or subtract values from the data numbers in
C      rectangular sections of an image.
C
C  CONVERTED FOR USE ON MIPL SYSTEM BY
C      STEVE POHORSKY   INFORMATICS GENERAL CORPORATION     OCTOBER 1983
C  FOR
C      EARTH RESOURCES APPLICATIONS
C  ENVIRONMENT
C      VAX 11/780    VMS  with TAE/VICAR2 EXECUTIVE       FORTRAN-77
C
C  REVISION HISTORY
C    04-91  REA  ADAPTED FOR UNIX/VICAR
C    07-85  JHR  CONVERTED TO VICAR2
C    10-83  SP   CONVERTED FROM IBM VICAR VERSION: MISCELLANEOUS CLEANUP.
C    10-83  SP   INCREASED MAXIMUM LINE SIZE TO 64000 PIXELS.
C    10-83  SP   CHANGED LABEL PROCESSING OF INPUT FILE BY USING LABELC
C                TO FIND THE DATA FORMAT AND USING FORMAT IN LABEL AS DEFAULT.
C    10-83  SP   CHANGED HALF PARAMETER TO FORMAT=HALF/BYTE SO USER CAN
C                OVERRIDE LABEL WHEN LABEL SAYS HALF.
C    10-83  SP   CHANGED LABEL PROCESSING OF OUTPUT FILE TO ALWAYS INCLUDE
C                CORRECT DATA FORMAT IN LABEL.
C    10-83  SP   MODIFIED PROGRAM SO THAT STARTING PIXEL VALUE IN THE  AREA
C                PARAMETER REFERS TO POSITION IN INPUT IMAGE EVEN WHEN THE
C                STARTING SAMPLE IN THE SIZE FIELD IS NOT 1. (THE STARTING
C                LINE VALUE IN THE  AREA PARAMETER ALREADY WAS SET TO REFER
C                TO THE POSITION IN INPUT IMAGE EVEN WHEN THE
C                STARTING LINE IN THE SIZE FIELD IS NOT 1. )
C    10-83  SP   CORRECTED PROGRAM TO NOT STORE VARIABLE INT INTO IBUF ARRAY
C                WHEN THE VALUE OF INT HAD NOT BEEN SET PROPERLY ( I.E. WHEN
C                LOOP INDEX K WAS NOT IN RANGE.)
C  CALLING SEQUENCE
C      Standard subroutine call and return.
C  INPUT AND OUTPUT PARAMETERS
C      SEE UNDER PROGRAM QSAR.
C
C  CALLED BY
C      QSAR
C  SUBROUTINES CALLED
C      The library routines  ABEND, XVOPEN, XVCLOSE, XVPARM, QPRINT,
C                            XVREAD, XVWRIT, XVCHECK, XVUNIT, XVGET
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C * 27 JUNE 1975...DAH...CHANGES FOR CONVERSION TO 360/OS
C *   25 JULY 79   ...JBS...   MAX LINE LENGTH TO 4000 BYTES
C *   5 JUNE 80   ...JAM...   ADD HALFWORD OPTION
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C   THIS PROGRAM FOLLOWS THE STANDARD FORTRAN NAMING CONVENTION FOR VARIABLES:
C   VARIABLES STARTING WITH I-N ARE INTEGERS UNLESS EXPLICITLY DECLARED.

      PARAMETER ( LSIZE_PAR = 64000 )      ! MAX NUMBER OF PIXELS PER LINE.
      INTEGER*4 OUNIT,STAT,SS,SLO,SSO,IPARM(1800),ES,ADJ
      INTEGER*2 IBUF(LSIZE_PAR)
      CHARACTER*8 FORMAT

C          OPEN INPUT DATA SET
      CALL XVUNIT(IUNIT,'INP',1,STAT,' ')
      CALL XVOPEN(IUNIT,STAT,'U_FORMAT','HALF','OPEN_ACT','SA',
     +            'IO_ACT','SA',' ')

C        GET DATA FORMAT AND CHECK
      CALL XVGET(IUNIT,STAT,'FORMAT',FORMAT,' ')
      IF(FORMAT.NE.'BYTE'.AND.FORMAT.NE.'HALF') THEN
	 CALL XVMESSAGE(' QSAR ACCEPTS BYTE AND HALFWORD DATA ONLY',' ')
	 CALL ABEND
      ENDIF

C        GET SIZE INFORMATION AND CHECK
      CALL XVSIZE(SLO,SSO,NLO,NSO,NLI,NSI)
      IF(SLO+NLO-1 .GT. NLI) THEN
	 CALL XVMESSAGE(
     +		' NUMBER OF LINES REQUESTED EXCEEDS INPUT SIZE',' ')
	 CALL ABEND
      ENDIF
      IF(SSO+NSO-1 .GT. NSI) THEN
	 CALL XVMESSAGE(
     +		' NUMBER OF SAMPLES REQUESTED EXCEEDS INPUT SIZE',' ')
	 CALL ABEND
      ENDIF

C        OPEN OUTPUT DATA SET
      CALL XVUNIT(OUNIT,'OUT',1,STAT,' ')
      CALL XVOPEN(OUNIT,STAT,'OP','WRITE','U_FORMAT','HALF',
     &          'U_NL',NLO,'U_NS',NSO,'OPEN_ACT','SA','IO_ACT','SA',' ')

C        PROCESS PARAMETERS
      CALL XVPARM('AREA',IPARM,ICOUNT,IDEF,1800)
      IF(MOD(ICOUNT,5).NE.0) THEN
	 CALL XVMESSAGE(
     +		'NUMBER OF PARAMETERS NOT A MULTIPLE OF FIVE',' ')
	 CALL ABEND
      ENDIF

      LLO=SLO+NLO-1                 ! LAST LINE.

C  COPY THE INPUT IMAGE TO OUTPUT MAKING THE REQUESTED MODIFICATIONS.

      DO 100 LINE=SLO,LLO
         CALL XVREAD(IUNIT,IBUF,STAT,'LINE',LINE,' ')

         IF(ICOUNT.EQ.0) GO TO 55      ! IF NO MODIFICATIONS, JUST COPY.

         DO 50 J = 1,ICOUNT,5
C              CHECK FOR MODS TO THIS LINE
            IF(LINE.LT.IPARM(J).OR.LINE.GE.IPARM(J)+IPARM(J+2)) GO TO 50
            SS=IPARM(J+1)                ! STARTING PIXEL.
            IF(SS.LT.1) THEN
		CALL XVMESSAGE(' INVALID SS PARAMETER SPECIFIED',' ')
		CALL ABEND
	    ENDIF
            ES=IPARM(J+3)+SS-1           ! ENDING PIXEL.
            IF(ES.GT.NSI) THEN
		CALL XVMESSAGE(' INPUT SAMPLE SIZE EXCEEDED',' ')
		CALL ABEND
	    ENDIF
            ADJ=IPARM(J+4)               ! ADJUSTMENT TO BE ADDED.
C
            DO K=SS,ES
               INT=IBUF(K)+ADJ
               IF(FORMAT.EQ.'HALF')  THEN
                  IF(INT.LT.-32768) INT=-32768
                  IF(INT.GT.32767)  INT= 32767
               ELSE
                  IF (INT.LT.0) INT=0
                  IF (INT.GT.255) INT=255
               END IF
               IBUF(K)=INT
            END DO

50       CONTINUE

55       CALL XVWRIT(OUNIT,IBUF(SSO),STAT,'NSAMPS',NSO,' ')
100   CONTINUE

C        CLOSE DATA SETS
      CALL XVCLOSE(IUNIT,STAT,' ')
      CALL XVCLOSE(OUNIT,STAT,' ')
C
      RETURN
      END

ext: f

max_stars_repo_name: NASA-AMMOS/VICAR

max_issues_repo_name: NASA-AMMOS/VICAR

max_forks_repo_name: NASA-AMMOS/VICAR


text:
      SUBROUTINE SB_SF_004(P,ANS_SUMMED)
C
C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
C     By the MadGraph5_aMC@NLO Development Team
C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
C
C
C     Return the sum of the split orders which are required in
C      orders.inc (BORN_ORDERS)
C     Also the values needed for the counterterms are stored in the
C      C_BORN_CNT common block
C
C
C     Process: u u~ > t t~ [ real = QCD QED ] QCD^2<=4 QED^2<=2
C     Process: c c~ > t t~ [ real = QCD QED ] QCD^2<=4 QED^2<=2
C     spectators: 2 3

C
C
C     CONSTANTS
C
      IMPLICIT NONE
      INCLUDE 'nexternal.inc'
      INTEGER NSQAMPSO
      PARAMETER (NSQAMPSO=3)
C
C     ARGUMENTS
C
      REAL*8 P(0:3,NEXTERNAL), ANS_SUMMED
C
C     VARIABLES
C
      INTEGER I,J
      INCLUDE 'orders.inc'
      REAL*8 ANS(0:NSQAMPSO)
      LOGICAL KEEP_ORDER_CNT(NSPLITORDERS, NSQAMPSO)
      COMMON /C_KEEP_ORDER_CNT/ KEEP_ORDER_CNT
      INTEGER AMP_ORDERS(NSPLITORDERS)
      DOUBLE PRECISION TINY
      PARAMETER (TINY = 1D-12)
      DOUBLE PRECISION MAX_VAL
C
C     FUNCTIONS
C
      INTEGER GETORDPOWFROMINDEX_B
      INTEGER ORDERS_TO_AMP_SPLIT_POS
C
C     BEGIN CODE
C
      CALL SB_SF_004_SPLITORDERS(P,ANS)
C     color-linked borns are called for QCD-type emissions
      ANS_SUMMED = 0D0
      MAX_VAL = 0D0

C     reset the amp_split_cnt array
      AMP_SPLIT_CNT(1:AMP_SPLIT_SIZE,1:2,1:NSPLITORDERS) = DCMPLX(0D0
     $ ,0D0)

      DO I = 1, NSQAMPSO
        MAX_VAL = MAX(MAX_VAL, ABS(ANS(I)))
      ENDDO

      DO I = 1, NSQAMPSO
        IF (KEEP_ORDER_CNT(QCD_POS, I)) THEN
          ANS_SUMMED = ANS_SUMMED + ANS(I)
          DO J = 1, NSPLITORDERS
            AMP_ORDERS(J) = GETORDPOWFROMINDEX_B(J, I)
C           take into account the fact that this is for QCD
            IF (J.EQ.QCD_POS) AMP_ORDERS(J) = AMP_ORDERS(J) + 2
          ENDDO
            !amp_split_cnt(orders_to_amp_split_pos(amp_orders),1,qcd_pos) = ans(I)
          IF(ABS(ANS(I)).GT.MAX_VAL*TINY)
     $      AMP_SPLIT_CNT(ORDERS_TO_AMP_SPLIT_POS(AMP_ORDERS),1
     $     ,QCD_POS) = ANS(I)
        ENDIF
      ENDDO

C     this is to avoid fake non-zero contributions
      IF (ABS(ANS_SUMMED).LT.MAX_VAL*TINY) ANS_SUMMED=0D0

      RETURN
      END

      SUBROUTINE SB_SF_004_SPLITORDERS(P1,ANS)
C
C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
C     By the MadGraph5_aMC@NLO Development Team
C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
C
C     RETURNS AMPLITUDE SQUARED SUMMED/AVG OVER COLORS
C     AND HELICITIES
C     FOR THE POINT IN PHASE SPACE P(0:3,NEXTERNAL-1)
C
C     Process: u u~ > t t~ [ real = QCD QED ] QCD^2<=4 QED^2<=2
C     Process: c c~ > t t~ [ real = QCD QED ] QCD^2<=4 QED^2<=2
C     spectators: 2 3

C
      IMPLICIT NONE
C
C     CONSTANTS
C
      INCLUDE 'nexternal.inc'
      INTEGER     NCOMB
      PARAMETER ( NCOMB=  16 )
      INTEGER NSQAMPSO
      PARAMETER (NSQAMPSO=3)
      INTEGER    THEL
      PARAMETER (THEL=NCOMB*12)
      INTEGER NGRAPHS
      PARAMETER (NGRAPHS=   3)
C
C     ARGUMENTS
C
      REAL*8 P1(0:3,NEXTERNAL-1),ANS(0:NSQAMPSO)
C
C     LOCAL VARIABLES
C
      INTEGER IHEL,IDEN,I,J
      DOUBLE PRECISION T(NSQAMPSO)
      INTEGER IDEN_VALUES(12)
      DATA IDEN_VALUES /36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36/
C
C     GLOBAL VARIABLES
C
      LOGICAL GOODHEL(NCOMB,12)
      COMMON /C_GOODHEL/ GOODHEL
      DOUBLE PRECISION SAVEMOM(NEXTERNAL-1,2)
      COMMON/TO_SAVEMOM/SAVEMOM
      LOGICAL CALCULATEDBORN
      COMMON/CCALCULATEDBORN/CALCULATEDBORN
      INTEGER NFKSPROCESS
      COMMON/C_NFKSPROCESS/NFKSPROCESS
C     ----------
C     BEGIN CODE
C     ----------
      IDEN=IDEN_VALUES(NFKSPROCESS)
      IF (CALCULATEDBORN) THEN
        DO J=1,NEXTERNAL-1
          IF (SAVEMOM(J,1).NE.P1(0,J) .OR. SAVEMOM(J,2).NE.P1(3,J))
     $      THEN
            CALCULATEDBORN=.FALSE.
            WRITE(*,*) 'Error in sb_sf: momenta not the same in the'
     $       //' born'
            STOP
          ENDIF
        ENDDO
      ELSE
        WRITE(*,*) 'Error in sb_sf: color_linked borns should be'
     $   //' called only with calculatedborn = true'
        STOP
      ENDIF
      DO I=0,NSQAMPSO
        ANS(I) = 0D0
      ENDDO
      DO IHEL=1,NCOMB
        IF (GOODHEL(IHEL,NFKSPROCESS)) THEN
          CALL B_SF_004(IHEL,T)
          DO I=1,NSQAMPSO
            ANS(I)=ANS(I)+T(I)
          ENDDO
        ENDIF
      ENDDO
      DO I=1,NSQAMPSO
        ANS(I)=ANS(I)/DBLE(IDEN)
        ANS(0)=ANS(0)+ANS(I)
      ENDDO
      END

      SUBROUTINE B_SF_004(HELL,ANS)
C
C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
C     By the MadGraph5_aMC@NLO Development Team
C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
C     RETURNS AMPLITUDE SQUARED SUMMED/AVG OVER COLORS
C     FOR THE POINT WITH EXTERNAL LINES W(0:6,NEXTERNAL-1)

C     Process: u u~ > t t~ [ real = QCD QED ] QCD^2<=4 QED^2<=2
C     Process: c c~ > t t~ [ real = QCD QED ] QCD^2<=4 QED^2<=2
C     spectators: 2 3

C
      IMPLICIT NONE
C
C     CONSTANTS
C
      INTEGER NAMPSO, NSQAMPSO
      PARAMETER (NAMPSO=2, NSQAMPSO=3)
      INTEGER     NGRAPHS
      PARAMETER ( NGRAPHS = 3 )
      INTEGER NCOLOR1, NCOLOR2
      PARAMETER (NCOLOR1=2, NCOLOR2=2)
      REAL*8     ZERO
      PARAMETER (ZERO=0D0)
      COMPLEX*16 IMAG1
      PARAMETER (IMAG1 = (0D0,1D0))
      INCLUDE 'nexternal.inc'
      INCLUDE 'born_nhel.inc'
C
C     ARGUMENTS
C
      INTEGER HELL
      REAL*8 ANS(NSQAMPSO)
C
C     LOCAL VARIABLES
C
      INTEGER I,J,M,N
      REAL*8 CF(NCOLOR2,NCOLOR1)
      COMPLEX*16 ZTEMP, AMP(NGRAPHS), JAMP1(NCOLOR1,NAMPSO),
     $  JAMP2(NCOLOR2,NAMPSO)
      COMPLEX*16 TMP_JAMP(1)
C
C     GLOBAL VARIABLES
C
      DOUBLE COMPLEX SAVEAMP(NGRAPHS,MAX_BHEL)
      COMMON/TO_SAVEAMP/SAVEAMP
      LOGICAL CALCULATEDBORN
      COMMON/CCALCULATEDBORN/CALCULATEDBORN
C
C     FUNCTION
C
      INTEGER SQSOINDEXB
C
C     COLOR DATA
C
      DATA (CF(I,  1),I=  1,  2) /9.000000000000000D+00
     $ ,3.000000000000000D+00/
      DATA (CF(I,  2),I=  1,  2) /3.000000000000000D+00
     $ ,9.000000000000000D+00/
C     ----------
C     BEGIN CODE
C     ----------
      IF (.NOT. CALCULATEDBORN) THEN
        WRITE(*,*) 'Error in b_sf: color_linked borns should be called'
     $   //' only with calculatedborn = true'
        STOP
      ELSEIF (CALCULATEDBORN) THEN
        DO I=1,NGRAPHS
          AMP(I)=SAVEAMP(I,HELL)
        ENDDO
      ENDIF
C     JAMPs contributing to orders QCD=2 QED=0
      JAMP1(1,1) = (1.666666666666667D-01)*AMP(2)
      JAMP1(2,1) = (-5.000000000000000D-01)*AMP(2)
C     JAMPs contributing to orders QCD=0 QED=2
      JAMP1(1,2) = (-1.000000000000000D+00)*AMP(1)+(
     $ -1.000000000000000D+00)*AMP(3)
C     JAMPs contributing to orders QCD=2 QED=0
      JAMP2(1,1) = (-2.777777777777778D-01)*AMP(2)
      JAMP2(2,1) = (1.666666666666667D-01)*AMP(2)
C     JAMPs contributing to orders QCD=0 QED=2
      TMP_JAMP(1) = AMP(1) +  AMP(3)  ! used 2 times
      JAMP2(1,2) = (1.666666666666667D-01)*TMP_JAMP(1)
      JAMP2(2,2) = (-5.000000000000000D-01)*TMP_JAMP(1)
      DO I = 1, NSQAMPSO
        ANS(I) = 0D0
      ENDDO
      DO M = 1, NAMPSO
        DO I = 1, NCOLOR1
          ZTEMP = (0.D0,0.D0)
          DO J = 1, NCOLOR2
            ZTEMP = ZTEMP + CF(J,I)*JAMP2(J,M)
          ENDDO
          DO N = 1, NAMPSO
            ANS(SQSOINDEXB(M,N))=ANS(SQSOINDEXB(M,N))+ZTEMP
     $       *DCONJG(JAMP1(I,N))
          ENDDO
        ENDDO
      ENDDO
      END


!        Generated by TAPENADE     (INRIA, Ecuador team)
!  tapenade 3.x
!
!  Differentiation of lz in forward (tangent) mode:
!   variations   of useful results: lz
!   with respect to varying inputs: *temp *propunit
!   Plus diff mem management of: temp:in propunit:in
! MIT License
!
! Copyright (c) 2020 SHEMAT-Suite
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.
!>    @brief calculates effective thermal conductivity of the two phase
!>    @param[in] i cell index, direction I0
!>    @param[in] j cell index, direction J0
!>    @param[in] k cell index, direction K0
!>    @param[in] ismpl local sample index
!>    @return  thermal conductivity                lz[W/(m*K)]
!>    @details
!>    calculates effective thermal conductivity of the two phase\n
!>    system matrix-porosity .\n
!>    input:\n
!>      porosity                            porlocal [-]\n
!>      pressure                            plocal [Mpa]\n
!>      temperature                         tlocal in [C]\n
DOUBLE PRECISION FUNCTION g_LZ(i, j, k, ismpl, lz)
  USE ARRAYS

  USE g_ARRAYS

  USE ICE
  USE MOD_TEMP
  IMPLICIT NONE
  INTEGER :: i, j, k, ui, ismpl
  DOUBLE PRECISION :: tlocal, lsolid, lfluid, porlocal, fm, fi, ff
  DOUBLE PRECISION :: g_tlocal, g_lsolid, g_lfluid, g_porlocal, &
& g_fm, g_fi, g_ff
  DOUBLE PRECISION :: t0, theta, dtheta, w0, lice
  DOUBLE PRECISION :: g_theta, g_dtheta, g_lice
  EXTERNAL LAMF, LAMM, LAMI, POR
  EXTERNAL g_LAMF, g_LAMM, g_LAMI, g_POR
  DOUBLE PRECISION :: LAMM, LAMF, LAMI, POR
  DOUBLE PRECISION :: g_LAMM, g_LAMF, g_LAMI, g_POR
  INTRINSIC ABS
  DOUBLE PRECISION :: abs0
  DOUBLE PRECISION :: pwr1
  DOUBLE PRECISION :: g_pwr1
  DOUBLE PRECISION :: pwr2
  DOUBLE PRECISION :: g_pwr2
  DOUBLE PRECISION :: pwr3
  DOUBLE PRECISION :: g_pwr3
  DOUBLE PRECISION :: tmpresult
  DOUBLE PRECISION :: temp0
  DOUBLE PRECISION :: lz
  ui = uindex(i, j, k)
  g_tlocal = g_temp(i, j, k, ismpl)
  tlocal = temp(i, j, k, ismpl)
  t0 = liq(i, j, k)
  IF (liq(i, j, k) - sol(i, j, k) .GE. 0.) THEN
    abs0 = liq(i, j, k) - sol(i, j, k)
  ELSE
    abs0 = -(liq(i, j, k)-sol(i, j, k))
  END IF
  w0 = abs0/2.d0
  CALL g_FTHETA(tlocal, g_tlocal, t0, w0, theta, g_theta, dtheta, &
&           g_dtheta, ismpl)
  g_lfluid = g_LAMF(i, j, k, ismpl, lfluid)
  g_lice = g_LAMI(i, j, k, ismpl, lice)
  g_lsolid = g_propunit(ui, idx_lz, ismpl)
  lsolid = propunit(ui, idx_lz, ismpl)
  g_lsolid = g_LAMM(lsolid, g_lsolid, tlocal, g_tlocal, tref, &
&   ismpl, tmpresult)
  lsolid = tmpresult
  g_porlocal = g_POR(i, j, k, ismpl, porlocal)
  g_fm = -g_porlocal
  fm = 1.d0 - porlocal
  g_ff = theta*g_porlocal + porlocal*g_theta
  ff = porlocal*theta
  g_fi = g_porlocal - g_ff
  fi = porlocal - ff
  temp0 = lsolid**fm
  IF (lsolid .LE. 0.0 .AND. (fm .EQ. 0.0 .OR. fm .NE. INT(fm))) THEN
    g_pwr1 = 0.D0
  ELSE IF (lsolid .LE. 0.0) THEN
    g_pwr1 = fm*lsolid**(fm-1)*g_lsolid
  ELSE
    g_pwr1 = fm*lsolid**(fm-1)*g_lsolid + temp0*LOG(lsolid)*g_fm
  END IF
  pwr1 = temp0
  temp0 = lfluid**ff
  IF (lfluid .LE. 0.0 .AND. (ff .EQ. 0.0 .OR. ff .NE. INT(ff))) THEN
    g_pwr2 = 0.D0
  ELSE IF (lfluid .LE. 0.0) THEN
    g_pwr2 = ff*lfluid**(ff-1)*g_lfluid
  ELSE
    g_pwr2 = ff*lfluid**(ff-1)*g_lfluid + temp0*LOG(lfluid)*g_ff
  END IF
  pwr2 = temp0
  temp0 = lice**fi
  IF (lice .LE. 0.0 .AND. (fi .EQ. 0.0 .OR. fi .NE. INT(fi))) THEN
    g_pwr3 = 0.D0
  ELSE IF (lice .LE. 0.0) THEN
    g_pwr3 = fi*lice**(fi-1)*g_lice
  ELSE
    g_pwr3 = fi*lice**(fi-1)*g_lice + temp0*LOG(lice)*g_fi
  END IF
  pwr3 = temp0
  g_lz = pwr3*(pwr2*g_pwr1+pwr1*g_pwr2) + pwr1*pwr2*g_pwr3
  lz = pwr1*pwr2*pwr3
  RETURN
END FUNCTION g_LZ


text:
      SUBROUTINE AG07BD( JOBE, N, M, A, LDA, E, LDE, B, LDB, C, LDC,
     $                   D, LDD, AI, LDAI, EI, LDEI, BI, LDBI, CI, LDCI,
     $                   DI, LDDI, INFO )
C
C     SLICOT RELEASE 5.5.
C
C     Copyright (c) 2002-2012 NICONET e.V.
C
C     PURPOSE
C
C     To compute the inverse (Ai-lambda*Ei,Bi,Ci,Di) of a given
C     descriptor system (A-lambda*E,B,C,D).
C
C     ARGUMENTS
C
C     Mode Parameters
C
C     JOBE    CHARACTER*1
C             Specifies whether E is a general square or an identity
C             matrix as follows:
C             = 'G':  E is a general square matrix;
C             = 'I':  E is the identity matrix.
C
C     Input/Output Parameters
C
C     N       (input) INTEGER
C             The order of the square matrices A and E;
C             also the number of rows of matrix B and the number of
C             columns of matrix C.  N >= 0.
C
C     M       (input) INTEGER
C             The number of system inputs and outputs, i.e., the number
C             of columns of matrices B and D and the number of rows of
C             matrices C and D.  M >= 0.
C
C     A       (input) DOUBLE PRECISION array, dimension (LDA,N)
C             The leading N-by-N part of this array must contain the
C             state matrix A of the original system.
C
C     LDA     INTEGER
C             The leading dimension of the array A.  LDA >= MAX(1,N).
C
C     E       (input) DOUBLE PRECISION array, dimension (LDE,N)
C             If JOBE = 'G', the leading N-by-N part of this array must
C             contain the descriptor matrix E of the original system.
C             If JOBE = 'I', then E is assumed to be the identity
C             matrix and is not referenced.
C
C     LDE     INTEGER
C             The leading dimension of the array E.
C             LDE >= MAX(1,N), if JOBE = 'G';
C             LDE >= 1,        if JOBE = 'I'.
C
C     B       (input) DOUBLE PRECISION array, dimension (LDB,M)
C             The leading N-by-M part of this array must contain the
C             input matrix B of the original system.
C
C     LDB     INTEGER
C             The leading dimension of the array B.  LDB >= MAX(1,N).
C
C     C       (input) DOUBLE PRECISION array, dimension (LDC,N)
C             The leading M-by-N part of this array must contain the
C             output matrix C of the original system.
C
C     LDC     INTEGER
C             The leading dimension of the array C.  LDC >= MAX(1,M).
C
C     D       (input) DOUBLE PRECISION array, dimension (LDD,M)
C             The leading M-by-M part of this array must contain the
C             feedthrough matrix D of the original system.
C
C     LDD     INTEGER
C             The leading dimension of the array D.  LDD >= MAX(1,M).
C
C     AI      (output) DOUBLE PRECISION array, dimension (LDAI,N+M)
C             The leading (N+M)-by-(N+M) part of this array contains
C             the state matrix Ai of the inverse system.
C             If LDAI = LDA >= N+M, then AI and A can share the same
C             storage locations.
C
C     LDAI    INTEGER
C             The leading dimension of the array AI.
C             LDAI >= MAX(1,N+M).
C
C     EI      (output) DOUBLE PRECISION array, dimension (LDEI,N+M)
C             The leading (N+M)-by-(N+M) part of this array contains
C             the descriptor matrix Ei of the inverse system.
C             If LDEI = LDE >= N+M, then EI and E can share the same
C             storage locations.
C
C     LDEI    INTEGER
C             The leading dimension of the array EI.
C             LDEI >= MAX(1,N+M).
C
C     BI      (output) DOUBLE PRECISION array, dimension (LDBI,M)
C             The leading (N+M)-by-M part of this array contains
C             the input matrix Bi of the inverse system.
C             If LDBI = LDB >= N+M, then BI and B can share the same
C             storage locations.
C
C     LDBI    INTEGER
C             The leading dimension of the array BI.
C             LDBI >= MAX(1,N+M).
C
C     CI      (output) DOUBLE PRECISION array, dimension (LDCI,N+M)
C             The leading M-by-(N+M) part of this array contains
C             the output matrix Ci of the inverse system.
C             If LDCI = LDC, CI and C can share the same storage
C             locations.
C
C     LDCI    INTEGER
C             The leading dimension of the array CI.  LDCI >= MAX(1,M).
C
C     DI      (output) DOUBLE PRECISION array, dimension (LDDI,M)
C             The leading M-by-M part of this array contains
C             the feedthrough matrix Di = 0 of the inverse system.
C             DI and D can share the same storage locations.
C
C     LDDI    INTEGER
C             The leading dimension of the array DI.  LDDI >= MAX(1,M).
C
C     Error Indicator
C
C     INFO    INTEGER
C             = 0:  successful exit;
C             < 0:  if INFO = -i, the i-th argument had an illegal
C                   value.
C
C     METHOD
C
C     The matrices of the inverse system are computed with the formulas
C
C                ( E  0 )        ( A  B )         (  0 )
C           Ei = (      ) , Ai = (      ) ,  Bi = (    ),
C                ( 0  0 )        ( C  D )         ( -I )
C
C           Ci = ( 0  I ),  Di = 0.
C
C     FURTHER COMMENTS
C
C     The routine does not perform an invertibility test. This check can
C     be performed by using the SLICOT routines AB08NX or AG08BY.
C
C     CONTRIBUTORS
C
C     A. Varga, German Aerospace Center, Oberpfaffenhofen, July 2000.
C
C     REVISIONS
C
C     V. Sima, Research Institute for Informatics, Bucharest, Mar. 2001.
C
C     KEYWORDS
C
C     Descriptor system, inverse system, state-space representation.
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
C     .. Scalar Arguments ..
      CHARACTER          JOBE
      INTEGER            INFO, LDA, LDAI, LDB, LDBI, LDC, LDCI,
     $                   LDD, LDDI, LDE, LDEI, M, N
C     .. Array Arguments ..
      DOUBLE PRECISION   A(LDA,*), AI(LDAI,*), B(LDB,*), BI(LDBI,*),
     $                   C(LDC,*), CI(LDCI,*), D(LDD,*), DI(LDDI,*),
     $                   E(LDE,*), EI(LDEI,*)
C     .. Local Scalars ..
      LOGICAL            UNITE
C     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
C     .. External Subroutines ..
      EXTERNAL           DLACPY, DLASET, XERBLA
C     .. Intrinsic Functions ..
      INTRINSIC          MAX
C     .. Executable Statements ..
C
      INFO = 0
C
C     Test the input scalar arguments.
C
      UNITE = LSAME( JOBE, 'I' )
      IF( .NOT. ( LSAME( JOBE, 'G' ) .OR. UNITE ) ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( M.LT.0 ) THEN
         INFO = -3
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -5
      ELSE IF( LDE.LT.1 .OR. ( .NOT.UNITE .AND. LDE.LT.N ) ) THEN
         INFO = -7
      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
         INFO = -9
      ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
         INFO = -11
      ELSE IF( LDD.LT.MAX( 1, M ) ) THEN
         INFO = -13
      ELSE IF( LDAI.LT.MAX( 1, N+M ) ) THEN
         INFO = -15
      ELSE IF( LDEI.LT.MAX( 1, N+M ) ) THEN
         INFO = -17
      ELSE IF( LDBI.LT.MAX( 1, N+M ) ) THEN
         INFO = -19
      ELSE IF( LDCI.LT.MAX( 1, M ) ) THEN
         INFO = -21
      ELSE IF( LDDI.LT.MAX( 1, M ) ) THEN
         INFO = -23
      END IF
C
      IF ( INFO.NE.0 ) THEN
C
C        Error return.
C
         CALL XERBLA( 'AG07BD', -INFO )
         RETURN
      END IF
C
C     Quick return if possible.
C
      IF ( M.EQ.0 )
     $   RETURN
C
C     Form Ai.
C
      CALL DLACPY( 'Full', N, N, A, LDA, AI, LDAI )
      CALL DLACPY( 'Full', M, N, C, LDC, AI(N+1,1), LDAI )
      CALL DLACPY( 'Full', N, M, B, LDB, AI(1,N+1), LDAI )
      CALL DLACPY( 'Full', M, M, D, LDD, AI(N+1,N+1), LDAI )
C
C     Form Ei.
C
      IF( UNITE ) THEN
         CALL DLASET( 'Full', N+M, N, ZERO, ONE, EI, LDEI )
      ELSE
         CALL DLACPY( 'Full', N, N, E, LDE, EI, LDEI )
         CALL DLASET( 'Full', M, N, ZERO, ZERO, EI(N+1,1), LDEI )
      END IF
      CALL DLASET( 'Full', N+M, M, ZERO, ZERO, EI(1,N+1), LDEI )
C
C     Form Bi.
C
      CALL DLASET( 'Full', N, M, ZERO, ZERO, BI, LDBI )
      CALL DLASET( 'Full', M, M, ZERO, -ONE, BI(N+1,1), LDBI )
C
C     Form Ci.
C
      CALL DLASET( 'Full', M, N, ZERO, ZERO, CI, LDCI )
      CALL DLASET( 'Full', M, M, ZERO, ONE, CI(1,N+1), LDCI )
C
C     Set Di.
C
      CALL DLASET( 'Full', M, M, ZERO, ZERO, DI, LDDI )
C
      RETURN
C *** Last line of AG07BD ***
      END

max_stars_repo_name: bgin/MissileSimulation

max_issues_repo_name: bgin/MissileSimulation

max_forks_repo_name: bgin/MissileSimulation


text:
C Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
C See https://llvm.org/LICENSE.txt for license information.
C SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception

	double precision a0(100), a1(100)
	double precision b0(100), b1(100)

	do i = 1,100
	 a0(i) = i
	 a1(i) = 3
	 b0(i) = 0
	 b1(i) = 0
	 b1(i) = mod(a0(i),a1(i))
	 b0(i) = mod(a0(i),a1(i))
	enddo

	call checkd( b0, b1, 100 )
	end

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang


! { dg-do run }
! PR61499
program read_internal

  integer :: x(9),i,iostat
  character(len=512) :: iomsg
  character(kind=1,len=30), dimension(:), allocatable, save :: source
  allocate(source(3))
  source=["  1   1  -1","  1  -1   1"," -1   1   1"]      !This fails
  read(source,*) (x(i), i=1,6)
end program read_internal


        !COMPILER-GENERATED INTERFACE MODULE: Fri Mar  6 15:34:34 2020
        ! This source file is for reference only and may not completely
        ! represent the generated interface used by the compiler.
        MODULE TWODEE_EMBM__genmod
          INTERFACE
            SUBROUTINE TWODEE_EMBM(TEMPER,TEMP1,IMAX,JMAX,IX,IY,SCALE)
              INTEGER(KIND=4) :: IY
              INTEGER(KIND=4) :: IX
              INTEGER(KIND=4) :: JMAX
              INTEGER(KIND=4) :: IMAX
              REAL(KIND=8) :: TEMPER(IMAX,JMAX)
              REAL(KIND=8) :: TEMP1(IX,IY)
              REAL(KIND=8) :: SCALE
            END SUBROUTINE TWODEE_EMBM
          END INTERFACE
        END MODULE TWODEE_EMBM__genmod


text:
! Issue # 92: https://github.com/OpenFortranProject/open-fortran-parser/issues/92
!
! When a function is declared with an initial type of 'character()', the parser gives the following error message:
! no viable alternative at input 'f'
! Such a declaration is legal -- see Fortran 90 Specification (N692) 5.1.1.5 rule R509, list item (3).
!
! If the alternative syntax 'character(len=*)' is used, the function parses correctly.
!
character(*) function name()
   name = "kindness"
end function

max_stars_repo_name: OpenFortranProject/ofp-sdf

max_issues_repo_name: OpenFortranProject/ofp-sdf

max_forks_repo_name: OpenFortranProject/ofp-sdf


text:
c     # We tag for coarsening if this coarsened patch isn't tagged for refinement
      subroutine clawpack5_tag4coarsening(mx,my,mbc,meqn,
     &      xlower,ylower,dx,dy, blockno, q0, q1, q2, q3,
     &      coarsen_threshold, init_flag, tag_patch)
      implicit none

      integer mx,my, mbc, meqn, tag_patch
      integer blockno, init_flag
      double precision xlower(0:3), ylower(0:3), dx, dy
      double precision coarsen_threshold
      double precision q0(meqn,1-mbc:mx+mbc,1-mbc:my+mbc)
      double precision q1(meqn,1-mbc:mx+mbc,1-mbc:my+mbc)
      double precision q2(meqn,1-mbc:mx+mbc,1-mbc:my+mbc)
      double precision q3(meqn,1-mbc:mx+mbc,1-mbc:my+mbc)

      integer mq
      double precision qmin, qmax

      if (init_flag .ne. 0) then
            tag_patch = 0
            return
      end if

c     # Assume that we will coarsen a family unless we find a grid
c     # that doesn't pass the coarsening test.
      tag_patch = 1
      mq = 1
      qmin = q0(mq,1,1)
      qmax = q0(mq,1,1)

c     # If we find that (qmax-qmin > coarsen_threshold) on any
c     # grid, we return immediately, since the family will then
c     # not be coarsened.

      call user5_get_minmax(blockno, mx,my,mbc,meqn,mq,q0,qmin,qmax,
     &      dx,dy,xlower(0),ylower(0),coarsen_threshold,tag_patch)
      if (tag_patch == 0) return

      call user5_get_minmax(blockno, mx,my,mbc,meqn,mq,q1,qmin,qmax,
     &      dx,dy,xlower(1),ylower(1),coarsen_threshold,tag_patch)
      if (tag_patch == 0) return

      call user5_get_minmax(blockno, mx,my,mbc,meqn,mq,q2,qmin,qmax,
     &      dx,dy,xlower(2),ylower(2),coarsen_threshold,tag_patch)
      if (tag_patch == 0) return

      call user5_get_minmax(blockno, mx,my,mbc,meqn,mq,q3,qmin,qmax,
     &      dx,dy,xlower(3),ylower(3),coarsen_threshold,tag_patch)

      end

      subroutine user5_get_minmax(blockno, mx,my,mbc,meqn,mq,q,
     &      qmin,qmax,dx,dy,xlower,ylower,coarsen_threshold,tag_patch)

      implicit none
      integer blockno, mx,my,mbc,meqn,mq,tag_patch
      double precision dx,dy,xlower, ylower
      double precision coarsen_threshold
      double precision qmin,qmax
      double precision q(meqn,1-mbc:mx+mbc,1-mbc:my+mbc)
      integer i,j

      double precision xc,yc,quad(-1:1,-1:1)
      integer ii,jj

      logical exceeds_th, minmax_exceeds_th

      do i = 1,mx
          do j = 1,my
              xc = xlower + (i-0.5)*dx
              yc = ylower + (j-0.5)*dy
              qmin = min(qmin,q(mq,i,j))
              qmax = max(qmax,q(mq,i,j))
              do ii = -1,1
                  do jj = -1,1
                      quad(ii,jj) = q(mq,i+ii,j+jj)
                  end do
              end do
              exceeds_th = minmax_exceeds_th(blockno,
     &                     q(mq,i,j),qmin,qmax,quad, dx,dy,xc,yc,
     &                     coarsen_threshold)

              if (exceeds_th) then
c                 # We won't coarsen this family because at least one
c                 # grid fails the coarsening test.
                  tag_patch = 0
                  return
              endif
         enddo
      enddo

      end

max_stars_repo_name: jbsnively/forestclaw

max_issues_repo_name: jbsnively/forestclaw

max_forks_repo_name: jbsnively/forestclaw


text:
      subroutine writepolar(ngroups,nch0,nat0,lab0,c0,sym,numb)

      implicit double precision(a-h,o-z)

      dimension nch0(ngroups),nat0(ngroups)

      dimension c0(ngroups,60,3)

      character*2 lab0(ngroups,60),sym(110)

      character*20 ca,ca1

      dimension num(60)

      dimension numb(110)

      write(24,*)' The number of polarisation calcs is'
      write(24,*)ngroups
      write(24,*)' The number of atoms and electrons in each calc is'
      do n=1,ngroups

c calculate the number of electrons
      call Atnum(60,nat0(n),lab0(n,:),num,sym,numb,nflag)
        if(nflag.eq.1)then
         write(6,*)' n ',n
         write(6,*)' nat0(n) ',nat0(n)
         do i=1,nat0(n)
          write(6,*)i,lab0(n,i)
         enddo
         stop
        endif
      ne=nch0(n)
      do i=1,nat0(n)
       ne=ne+num(i)
      enddo

      ne2=ne/2
      ne3=ne-2*ne2
      if(ne3.eq.0)then
       multiplicity=1
      else
       multiplicity=2
      endif
c make the filename
       call filelabel(n,ca)
       n1=index(ca,' ')-1
       ca1='nb.'//ca(1:n1)//'.0-polar.com'
#ifdef __GFORTRAN__
       open(unit=1,file=ca1,status='unknown')
#else
       open(unit=1,file=ca1,status='unknown',buffered='YES')
#endif
       write(1,*)nch0(n),multiplicity
       do i=1,nat0(n)
        write(1,100)lab0(n,i),(c0(n,i,k),k=1,3)
       enddo
100    format(a2,3f13.6)
c add a blank line
       write(1,*)
c end the loop over groups
       close(unit=1)
c write out the number of electrons for information
c     write(6,*)'writing 24'
      write(24,*)nat0(n),ne,ca1

      enddo

      return
      end

max_stars_repo_name: mickcollins/SMFA

max_issues_repo_name: mickcollins/SMFA

max_forks_repo_name: mickcollins/SMFAPAC


text:
subroutine ex_elastic(C,lC,nC,nT,col_temps,j_lev,e_levs,e_num,C0)
!subroutine to obtain the elastic collision rates
implicit none
integer nC,nT,nP,e_num
integer dj,d_j,max_dj,min_dj
integer i,j,i1,i2,j1,j2

integer, dimension(2,nC) :: lC
integer, dimension (e_num) :: j_lev
integer, dimension (nC) :: dj_nums

double precision dE,E1,E2,g1,g2,x,y,Temp,C0

double precision, dimension(nT) :: col_temps
double precision, dimension(nT,nC) :: C
double precision, dimension(e_num) :: e_levs
double precision, dimension(2,2) :: B
double precision, dimension(2) :: a
double precision, allocatable, dimension(:,:) :: const

!based on Keto & Rybicki 2005
!and de Jong et al. 1970

min_dj = 100
max_dj = 0
do i = 1,nC
  dj = lC(2,i) - lC(1,i)
  if (dj.gt.max_dj) max_dj = dj
  if (dj.lt.min_dj) min_dj = dj
end do

allocate(const(2,min_dj:max_dj))
!now were gonna find the a(\Delta J) and b(\Delta J) parameters
do dj = min_dj,max_dj
  B(:,:) = 0.d0
  a(:)   = 0.d0

!  write(*,*)dj
  do i = 1,nC
    i1 = lC(1,i)
    i2 = lC(2,i)

    j1 = j_lev(i1)
    j2 = j_lev(i2)

    g1 = 1.d0 * (2.d0 * j1 + 1.d0)
    g2 = 1.d0 * (2.d0 * j2 + 1.d0)

    E1 = e_levs(i1)
    E2 = e_levs(i2)

    d_j = j2 - j1
    if (dj.eq.d_j) then
      dE = E2 - E1
      do j = 1,nT
        Temp = col_temps(j)
        x = 4.779E-11 * dE/Temp   !h/kB * \Delta v / T
        y = dlog(C(j,i) * ((g2/g1)*(1.d0+x))**(-1.d0))

        B(1,1) = B(1,1) + 1.d0
        B(1,2) = B(1,2) - dsqrt(x)
        B(2,2) = B(2,2) + x

        a(1) = a(1) + y
        a(2) = a(2) - dsqrt(x) * y

!        write(*,*)dsqrt(x),y

      end do
    endif
  end do

  B(2,1) = B(1,2)

  !now solve the 2x2 inversion problem
  call invert_22(B)

  const(1,dj) = B(1,1)*a(1) + B(1,2)*a(2)       !ln(a(\Delta J))
  const(2,dj) = B(2,1)*a(1) + B(2,2)*a(2)       !b(\Delta J)
!  write(*,*)
  write(*,*)dj,dexp(const(1,dj)),const(2,dj)
!  write(*,*)
!  write(*,*)
end do

!now let's fit the a-constant

B(:,:) = 0.d0
do dj = min_dj,max_dj
  B(1,1) = B(1,1) + 1.d0
  B(1,2) = B(1,2) + 1.d0 * dj
  B(2,2) = B(2,2) + 1.d0 * dj*dj

  a(1) = a(1) + 1.d0 * const(1,dj)
  a(2) = a(2) + dj * const(1,dj)
end do

call invert_22(B)

C0 = B(1,1) * a(1) + B(1,2) * a(2)
C0 = dexp(C0)

!write(*,*)C0
end subroutine

subroutine invert_22(A)
!subroutine to invert a 2x2 matrix
implicit none

double precision det,a11

double precision, dimension(2,2) :: A

det = (A(1,1) * A(2,2) - A(2,1) * A(1,2))**(-1.d0)

a11 = A(1,1)

A(1,1) = det * A(2,2)
A(1,2) = -det * A(1,2)
A(2,1) = -det * A(2,1)
A(2,2) = det * a11

end subroutine

subroutine ex_elastic_singleJ(J_state,C,lC,nC,nT,col_temps,j_lev,e_levs,e_num,C0)
!subroutine to obtain the elastic collision rates for a single level J by
!only considering the temperature variation
implicit none
integer nC,nC_red,nT,nP,e_num
integer dj,d_j,max_dj,min_dj
integer i,j,i1,i2,j1,j2,J_state

integer, dimension(2,nC) :: lC
integer, dimension(2,nC) :: lC_red
integer, dimension(nC) :: red_to_real
integer, dimension (e_num) :: j_lev
integer, dimension (nC) :: dj_nums

double precision dE,E1,E2,g1,g2,x,y,Temp,C0

double precision, dimension(nT) :: col_temps
double precision, dimension(nT,nC) :: C
double precision, dimension(e_num) :: e_levs
double precision, dimension(2,2) :: B
double precision, dimension(2) :: a
double precision, allocatable, dimension(:,:) :: const

!based on Keto & Rybicki 2005
!and de Jong et al. 1970

min_dj = 100
max_dj = 0

nC_red = 0
do i = 1,nC
  if (lC(1,i).eq.J_state) then
    nC_red = nC_red + 1
    dj = lC(2,i) - lC(1,i)

    red_to_real(nC_red) = i

    lC_red(1,nC_red) = lC(1,i)
    lC_red(2,nC_red) = lC(2,i)

    if (dj.gt.max_dj) max_dj = dj
    if (dj.lt.min_dj) min_dj = dj
  endif
end do

allocate(const(2,min_dj:max_dj))
!now were gonna find the a(\Delta J) and b(\Delta J) parameters
do dj = min_dj,max_dj
  B(:,:) = 0.d0
  a(:)   = 0.d0

  write(*,*)dj
  do i = 1,nC_red
    i1 = lC_red(1,i)
    i2 = lC_red(2,i)

    j1 = j_lev(i1)
    j2 = j_lev(i2)

    g1 = 1.d0 * (2.d0 * j1 + 1.d0)
    g2 = 1.d0 * (2.d0 * j2 + 1.d0)

    E1 = e_levs(i1)
    E2 = e_levs(i2)

    d_j = j2 - j1
    if (dj.eq.d_j) then
      dE = E2 - E1
      do j = 1,nT
        Temp = col_temps(j)
        x = 4.779E-11 * dE/Temp   !h/kB * \Delta v / T
        y = dlog(C(j,red_to_real(i)) * (g2/g1) * ((1.d0+x))**(-1.d0))

        B(1,1) = B(1,1) + 1.d0
        B(1,2) = B(1,2) - dsqrt(x)
        B(2,2) = B(2,2) + x

        a(1) = a(1) + y
        a(2) = a(2) - dsqrt(x) * y

        write(*,*)dsqrt(x),y

      end do
    endif
  end do

  B(2,1) = B(1,2)

  !now solve the 2x2 inversion problem
  call invert_22(B)

  const(1,dj) = B(1,1)*a(1) + B(1,2)*a(2)       !ln(a(\Delta J))
  const(2,dj) = B(2,1)*a(1) + B(2,2)*a(2)       !b(\Delta J)
!  write(*,*)
!  write(*,*)const(1,dj),const(2,dj)
!  write(*,*)
!  write(*,*)
end do

!now let's fit the a-constant

B(:,:) = 0.d0
do dj = min_dj,max_dj
  B(1,1) = B(1,1) + 1.d0
  B(1,2) = B(1,2) + 1.d0 * dj
  B(2,2) = B(2,2) + 1.d0 * dj*dj

  a(1) = a(1) + 1.d0 * dlog(const(1,dj))
  a(2) = a(2) + dj * dlog(const(1,dj))
end do

call invert_22(B)

C0 = B(1,1) * a(1) + B(1,2) * a(2)

end subroutine

subroutine de_polarization(C,lC,nC,nT,e_num,kmax,D0)
!subroutine to obtain the de-polarization rates in the sudden
!approximation
use fwigxjpf
implicit none
integer nC,nT,e_num
integer i,j,k,l,m,kmax

integer, dimension(2,nC) :: lC
integer, dimension (e_num) :: j_lev
integer, dimension (nC) :: dj_nums

double precision lj2,fac,w1,w2,w3,fac2

double precision, dimension(nT,nC) :: C
double precision, dimension(nT,e_num,2:kmax) :: D0
double precision, dimension(nT,e_num) :: C_j0

!input:
!C      -- collision rates
!lC     -- collisional transitions
!nC     -- number of collisional transitions
!nT     -- number of collision temperatures tabulated
!e_num  -- number of levels
!kmax   -- maximum number of irreducible tensor-rank

!output:
!D0     -- all depolarization coefficients, temperature dependent

!based on Morris 1985 -- Eq. (B11)
!same definition of D^(K)(\alpha J) as was
!used in L&L06, Eq. (7.102)

do i = 1,nC
  if (lC(1,i).eq.1) C_j0(:,lC(2,i)) = C(:,i)
end do
!if lC(1,i) is the zero-level, then use these as the sudden-approximation
!paramaters
!C_j0(:,1) = 0 -- we don't need it

!write(*,*)C_j0(1,2)
!write(*,*)C_j0(1,3)

!only works if j_lev(i) = i-1
D0(:,:,:) = 0.d0
do i = 2,e_num
  j = i-1
  lj2 = (2.d0*j+1.d0)**2.d0
  do k = 2,min(2*j,kmax),2
    do l = 1,nT
      D0(l,i,k) = 0.d0
      !Wigner 3j = 0 if J+J < j
      do m = 1,2*j

        fac = (-1.d0)**(m) * (2.d0*m+1.d0)

        w1 = fwig3jj (2* j , 2* j    , 2* m ,&
        &              2* 0 , 2* 0 , 0    )

        w2 = fwig6jj (2* j , 2* j , 2* m ,&
        &              2* j , 2* j  , 2* 0    )

        w3 = fwig6jj (2* j , 2* j , 2* m ,&
        &              2* j , 2* j  , 2* k    )

        fac2 = w2 - w3

        D0(l,i,k) = D0(l,i,k) + lj2 * fac * w1 * w1 * fac2 * C_j0(l,m)
!        write(*,*)m,k
!        write(*,*)lj2 * fac * w1 * w1 * fac2,C_j0(l,m+1)
      end do
    end do
  end do
end do

end subroutine

subroutine make_depolarization_temp(D0_all,nlev,nT,kmax,col_temps,Temp,nH2,D0)
!interpolate the temperature-dependent de-polarization coefficients
implicit none
integer nlev,nT,kmax,i

double precision Temp,ddT,nH2

double precision, dimension(nT,nlev,2:kmax) :: D0_all
double precision, dimension(nlev,2:kmax) :: D0
double precision, dimension(nlev,2:kmax) :: f1,f2
double precision, dimension(nT) :: col_temps

!input:
!D0_all         -- all depolarization coefficients, temperature dependent
!nlev           -- number of levels
!nT             -- number of temperatures
!kmax           -- maximum irreducible tensor rank
!col_temps      -- tabulated collision temperatures
!Temp           -- temperature we want coefficients for
!nH2            -- H2 number density

!output
!D0             -- depolarization coefficients for specific temperature

if (Temp.le.col_temps(1)) D0(:,:) = D0_all(1,:,:)

do i = 2,nT
  if (col_temps(i).gt.Temp) then
    ddT = (Temp - col_temps(i-1)) / (col_temps(i) - col_temps(i-1))

    f1(:,:) = D0_all(i-1,:,:)
    f2(:,:) = D0_all(i,:,:)

    D0(:,:) = f1(:,:) + ddT*(f2(:,:) - f1(:,:))
  endif
end do

D0(:,:) = D0(:,:) * nH2 * 1d-6

end subroutine

max_stars_repo_name: blankhaar/PORTAL

max_issues_repo_name: blankhaar/PORTAL

max_forks_repo_name: blankhaar/PORTAL


text:
c This file has been generated automatically by FeynRules

c Version: 1.4   Date: 25. 10. 2010,    11:39

      subroutine coup(readlha)

      implicit none
      logical readlha

      CALL coup1(readlha)
      CALL coup2(readlha)

      return
      end


text:
C................................................................       DGS   10
      SUBROUTINE DGSUB(I, Z, DG)                                        DGS   20
      REAL*8 Z(4), DG(4)
      DO 10 J=1,4
        DG(J) = 0.
   10 CONTINUE
      GO TO (20, 30, 20, 30), I
   20 DG(1) = 1.
      RETURN
   30 DG(3) = 1.
      RETURN
      END

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg


text:
! { dg-do run }

! PR fortran/35681
! Check that dependencies of MVBITS arguments are resolved correctly by using
! temporaries if both arguments refer to the same variable.

  integer, dimension(10) :: ila1 = (/1,2,3,4,5,6,7,8,9,10/)
  integer, dimension(20) :: ila2
  integer, dimension(10), target :: ila3
  integer, pointer :: ila3_ptr(:)
  integer, parameter :: SHOULD_BE(10) = (/17,18,11,4,13,22,7,16,9,18/)
  integer, parameter :: INDEX_VECTOR(10) = (/9,9,6,2,4,9,2,9,6,10/)

  ila2(2:20:2) = ila1
  ila3 = ila1

  ! Argument is already packed.
  call mvbits (ila1(INDEX_VECTOR), 2, 4, ila1, 3)
  write (*,'(10(I3))') ila1
  if (any (ila1 /= SHOULD_BE)) call abort ()

  ! Argument is not packed.
  call mvbits (ila2(2*INDEX_VECTOR), 2, 4, ila2(2:20:2), 3)
  write (*,'(10(I3))') ila2(2:20:2)
  if (any (ila2(2:20:2) /= SHOULD_BE)) call abort ()

  ! Pointer and target
  ila3_ptr => ila3
  call mvbits (ila3(INDEX_VECTOR), 2, 4, ila3_ptr, 3)
  write (*,'(10(I3))') ila3
  if (any (ila3 /= SHOULD_BE)) call abort ()

  end

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
!                          Futility Development Group                          !
!                             All rights reserved.                             !
!                                                                              !
! Futility is a jointly-maintained, open-source project between the University !
! of Michigan and Oak Ridge National Laboratory.  The copyright and license    !
! can be found in LICENSE.txt in the head directory of this repository.        !
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
PROGRAM testGeom_Points
#include "UnitTest.h"
USE ISO_FORTRAN_ENV
USE UnitTest
USE IntrType
USE Constants_Conversion
USE ParameterLists
USE Geom_Points
USE Geom

IMPLICIT NONE

TYPE(PointType) :: point,point2,point3
TYPE(PointType) :: points(2),points2(2),points3(2)
TYPE(LinkedListPointType),POINTER :: firstPoint,thisPoint
INTEGER(SIK) :: i
REAL(SRK) :: d,s(2)
LOGICAL(SBK) :: bool

CREATE_TEST('Test Geom')
CALL eParams%setQuietMode(.TRUE.)
CALL eParams%setStopOnError(.FALSE.)

REGISTER_SUBTEST('Test Points',TestPoints)

FINALIZE_TEST()
!
!===============================================================================
CONTAINS
!
!-------------------------------------------------------------------------------
SUBROUTINE TestPoints

  !Initialize by hand
  COMPONENT_TEST('%clear()')
  !or point=PointType(2,(/0.5_SRK,0.3_SRK/))
  point%dim=1
  ALLOCATE(point%coord(1))
  point%coord(1)=0.5_SRK
!
!Test clear routine
  CALL point%clear()
  bool = .NOT.(point%dim /= 0 .OR. ALLOCATED(point%coord))
  ASSERT(bool, 'point%clear()')
  WRITE(*,*) '  Passed: CALL point%clear()'
  CALL point%clear() !Test redundant call doesn't have an error
!
!Test initialization
  COMPONENT_TEST('%init()')
  CALL point%init(DIM=1,X=0.5_SRK) !test 1-D
  bool = .NOT.(point%dim /= 1 .OR. point%coord(1) /= 0.5_SRK)
  ASSERT(bool, 'point%init(DIM=1,X=0.5)')

  !Test redundant call
  CALL point%init(DIM=2,X=0.6_SRK,Y=0.7_SRK)
  bool = .NOT.(point%dim /= 1 .OR. point%coord(1) /= 0.5_SRK)
  ASSERT(bool, 'point%init(DIM=1,X=0.5)')
  CALL point%clear()

  CALL point%init(DIM=2,X=0.2_SRK,Y=0.3_SRK) !test 2-D
  bool = .NOT.(point%dim /= 2 .OR. point%coord(1) /= 0.2_SRK .OR. &
      point%coord(2) /= 0.3_SRK)
  ASSERT(bool, 'point%init(DIM=2,X=0.2_SRK,Y=0.3_SRK)')
  CALL point%clear()

  CALL point%init(DIM=3,X=0.7_SRK,Y=0.8_SRK,Z=0.9_SRK) !test 3-D
  bool = .NOT.(point%dim /= 3 .OR. point%coord(1) /= 0.7_SRK .OR. &
      point%coord(2) /= 0.8_SRK .OR. point%coord(3) /= 0.9_SRK)
  ASSERT(bool, 'point%init(DIM=3,X=0.7_SRK,Y=0.8_SRK,Z=0.9_SRK)')
  CALL point%clear()

  CALL point%init(COORD=(/1.7_SRK,1.8_SRK,1.9_SRK,1.0_SRK/)) !test N-D
  bool = .NOT.(point%dim /= 4 .OR. point%coord(1) /= 1.7_SRK .OR. &
      point%coord(2) /= 1.8_SRK .OR. point%coord(3) /= 1.9_SRK .OR. &
      point%coord(4) /= 1.0_SRK)
  ASSERT(bool, 'point%init(COORD=(/.../))')
!
!Test getCoordString
  COMPONENT_TEST('%getCoordString()')
  bool = .NOT.(TRIM(point%getCoordString()) /= '( 1.70000000000000E+00, '// &
      '1.80000000000000E+00, 1.90000000000000E+00, '// &
      '1.00000000000000E+00)' .OR. TRIM(point%getCoordString( &
      'f6.3')) /= '( 1.700, 1.800, 1.900, 1.000)')
  ASSERT(bool, 'point%getCoordString(...))')
#ifdef DEBUG
  bool = .NOT.(LEN_TRIM(point%getCoordString('i6')) /= 0 .OR. &
      LEN_TRIM(point2%getCoordString()) /= 0)
#else
  bool = .NOT.(LEN_TRIM(point2%getCoordString()) /= 0)
#endif
  ASSERT(bool, 'point%getCoordString(...))')
!
!Test Distance
  COMPONENT_TEST('Distance()')
  CALL point2%init(DIM=1,X=0.5_SRK)
  CALL point3%init(DIM=1,X=0.1_SRK)
  d=Distance(point2,point3)
  ASSERT(d == 0.4_SRK, '1-D Distance(...)') !Check 1-D
  CALL point2%clear()
  CALL point3%clear()

  CALL point2%init(DIM=2,X=0.5_SRK,Y=0.6_SRK)
  CALL point3%init(DIM=2,X=0.3_SRK,Y=0.2_SRK)
  d=Distance(point2,point3)
  bool = d .APPROXEQ. 0.447213595499958_SRK
  ASSERT(bool, '2-D Distance(...)')
  CALL point2%clear()
  CALL point3%clear()

  CALL point2%init(DIM=3,X=0.5_SRK,Y=0.6_SRK,Z=0.7_SRK)
  CALL point3%init(DIM=3,X=0.3_SRK,Y=0.2_SRK,Z=0.1_SRK)
  d=Distance(point2,point3)
  bool = d .APPROXEQ. 0.748331477354788_SRK
  ASSERT(bool, '3-D Distance(...)')

  !Redundant call to test error check
  ASSERT(Distance(point2,point) == 0.0_SRK, 'Distance(point2,point)')
  CALL point2%clear()
  CALL point3%clear()

  CALL point2%init(DIM=4,COORD=(/0.5_SRK,0.6_SRK,0.7_SRK,0.8_SRK/))
  CALL point3%init(DIM=4,COORD=(/0.3_SRK,0.2_SRK,0.1_SRK,0.0_SRK/))
  d=Distance(point2,point3)
  bool = d .APPROXEQ. 1.09544511501033_SRK
  ASSERT(bool, 'N-D Distance(...)')
  CALL point2%clear()
  CALL point3%clear()
  d=Distance(point2,point3) !Test for empty points
!
!Test midpoint
  COMPONENT_TEST('midPoint()')
  CALL point2%init(DIM=1,X=0.5_SRK)
  CALL point3%init(DIM=1,X=0.1_SRK)
  point=midPoint(point2,point3)
  bool = point%coord(1) .APPROXEQ. 0.3_SRK
  ASSERT(bool, '1-D midPoint(...)')
  CALL point2%clear()
  CALL point3%clear()

  !Test 2-D
  CALL point2%init(DIM=2,X=0.5_SRK,Y=0.6_SRK)
  CALL point3%init(DIM=2,X=0.3_SRK,Y=0.2_SRK)
  point=midPoint(point2,point3)
  bool = .NOT.(.NOT.(point%coord(1) .APPROXEQ. 0.4_SRK) .OR. &
      .NOT.(point%coord(2) .APPROXEQ. 0.4_SRK))
  ASSERT(bool, '2-D midPoint(...)')
  CALL point2%clear()
  CALL point3%clear()

  !Test 3-D
  CALL point2%init(DIM=3,X=0.5_SRK,Y=0.6_SRK,Z=0.7_SRK)
  CALL point3%init(DIM=3,X=0.3_SRK,Y=0.2_SRK,Z=0.1_SRK)
  point=midPoint(point2,point3)
  bool = .NOT.(.NOT.(point%coord(1) .APPROXEQ. 0.4_SRK) .OR. &
      .NOT.(point%coord(2) .APPROXEQ. 0.4_SRK) .OR. &
      .NOT.(point%coord(3) .APPROXEQ. 0.4_SRK))
  ASSERT(bool, '3-D midPoint(...)')
  CALL point2%clear()
  CALL point3%clear()

  !Test N-D
  CALL point2%init(DIM=4,COORD=(/0.5_SRK,0.6_SRK,0.7_SRK,0.8_SRK/))
  CALL point3%init(DIM=4,COORD=(/0.3_SRK,0.2_SRK,0.1_SRK,0.0_SRK/))
  point=midPoint(point2,point3) !Test for incorrect dimensions
  CALL point3%init(DIM=4,COORD=(/0.3_SRK,0.2_SRK,0.1_SRK,0.0_SRK/))
  point=midPoint(point2,point3)
  bool = .NOT.(.NOT.(point%coord(1) .APPROXEQ. 0.4_SRK) .OR. &
      .NOT.(point%coord(2) .APPROXEQ. 0.4_SRK) .OR. &
      .NOT.(point%coord(3) .APPROXEQ. 0.4_SRK) .OR. &
      .NOT.(point%coord(4) .APPROXEQ. 0.4_SRK))
  ASSERT(bool, 'N-D midPoint(...)')
  CALL point2%clear()
  CALL point3%clear()
  point=midPoint(point2,point3) !Test for empty points
!
!Test innerAngle
  COMPONENT_TEST('innerAngle()')
  CALL point%clear()
  CALL point%init(DIM=2,X=0.0_SRK,Y=2.0_SRK)
  CALL point2%init(DIM=2,X=0.0_SRK,Y=0.0_SRK)
  CALL point3%init(DIM=2,X=5.0_SRK,Y=0.0_SRK)
  ASSERT(innerAngle(point,point2,point3) .APPROXEQA. HALFPI,'90 right angle')
  CALL point%clear()
  CALL point2%clear()
  CALL point3%clear()
  CALL point%init(DIM=2,X=1.75_SRK,Y=2.134189828107240_SRK)
  CALL point2%init(DIM=2,X=0.25_SRK,Y=1.5_SRK)
  CALL point3%init(DIM=2,X=1.75_SRK,Y=0.865810171892757_SRK)
  ASSERT(innerAngle(point,point2,point3) .APPROXEQA. 0.8_SRK,'acute angle')
  CALL point%clear()
  CALL point2%clear()
  CALL point3%clear()
  CALL point%init(DIM=2,X=1.5_SRK+2.0_SRK*COS(QTRPI*0.25_SRK), &
      Y=0.5_SRK+2.0_SRK*SIN(QTRPI*0.25_SRK))
  CALL point2%init(DIM=2,X=1.5_SRK,Y=0.5_SRK)
  CALL point3%init(DIM=2,X=1.5_SRK+4.0_SRK*COS(QTRPI*0.25_SRK), &
      Y=0.5_SRK+4.0_SRK*SIN(-QTRPI*0.25_SRK))
  ASSERT(innerAngle(point,point2,point3) .APPROXEQA. QTRPI*0.5_SRK,'45 acute angle')
  CALL point%clear()
  CALL point2%clear()
  CALL point3%clear()
  CALL point%init(DIM=2,X=-40.0_SRK,Y=61.5_SRK)
  CALL point2%init(DIM=2,X=20.0_SRK,Y=1.5_SRK)
  CALL point3%init(DIM=2,X=60.0_SRK,Y=1.5_SRK)
  ASSERT(innerAngle(point,point2,point3) .APPROXEQA. PI-QTRPI,'135 obtuse angle')
  CALL point%clear()
  CALL point2%clear()
  CALL point3%clear()
  CALL point%init(DIM=2,X=-40.0_SRK,Y=0.5_SRK)
  CALL point2%init(DIM=2,X=20.0_SRK,Y=0.5_SRK)
  CALL point3%init(DIM=2,X=60.0_SRK,Y=0.5_SRK)
  ASSERT(innerAngle(point,point2,point3) .APPROXEQA. PI,'180 straight angle')
!
!Test outerAngle
  COMPONENT_TEST('outerAngle()')
  CALL point%clear()
  CALL point2%clear()
  CALL point3%clear()
  CALL point%init(DIM=2,X=0.0_SRK,Y=2.0_SRK)
  CALL point2%init(DIM=2,X=0.0_SRK,Y=0.0_SRK)
  CALL point3%init(DIM=2,X=5.0_SRK,Y=0.0_SRK)
  ASSERT(outerAngle(point,point2,point3) .APPROXEQA. TWOPI-HALFPI,'270 right angle')
  CALL point%clear()
  CALL point2%clear()
  CALL point3%clear()
  CALL point%init(DIM=2,X=1.75_SRK,Y=2.134189828107240_SRK)
  CALL point2%init(DIM=2,X=0.25_SRK,Y=1.5_SRK)
  CALL point3%init(DIM=2,X=1.75_SRK,Y=0.865810171892757_SRK)
  ASSERT(outerAngle(point,point2,point3) .APPROXEQA. TWOPI-0.8_SRK,'acute angle')
  CALL point%clear()
  CALL point2%clear()
  CALL point3%clear()
  CALL point%init(DIM=2,X=1.5_SRK+2.0_SRK*COS(QTRPI*0.25_SRK), &
      Y=0.5_SRK+2.0_SRK*SIN(QTRPI*0.25_SRK))
  CALL point2%init(DIM=2,X=1.5_SRK,Y=0.5_SRK)
  CALL point3%init(DIM=2,X=1.5_SRK+4.0_SRK*COS(QTRPI*0.25_SRK), &
      Y=0.5_SRK+4.0_SRK*SIN(-QTRPI*0.25_SRK))
  ASSERT(outerAngle(point,point2,point3) .APPROXEQA. TWOPI-QTRPI*0.5_SRK,'315 acute angle')
  CALL point%clear()
  CALL point2%clear()
  CALL point3%clear()
  CALL point%init(DIM=2,X=-40.0_SRK,Y=61.5_SRK)
  CALL point2%init(DIM=2,X=20.0_SRK,Y=1.5_SRK)
  CALL point3%init(DIM=2,X=60.0_SRK,Y=1.5_SRK)
  ASSERT(outerAngle(point,point2,point3) .APPROXEQA. PI+QTRPI,'225 obtuse angle')
  CALL point%clear()
  CALL point2%clear()
  CALL point3%clear()
  CALL point%init(DIM=2,X=-40.0_SRK,Y=0.5_SRK)
  CALL point2%init(DIM=2,X=20.0_SRK,Y=0.5_SRK)
  CALL point3%init(DIM=2,X=60.0_SRK,Y=0.5_SRK)
  ASSERT(outerAngle(point,point2,point3) .APPROXEQA. PI,'180 straight angle')

!
!Test Operators
  COMPONENT_TEST('OPERATOR(+)')
  CALL point%clear()
  CALL point2%clear()
  CALL point3%clear()
  CALL point2%init(DIM=1,X=0.5_SRK)
  CALL point3%init(DIM=1,X=0.1_SRK)
  point=point2+point3
  bool = .NOT.(point%dim /= 1 .OR. point%coord(1) /= 0.6_SRK)
  ASSERT(bool, '1-D PointType')
  CALL point2%clear()
  CALL point3%clear()

  CALL point2%init(DIM=2,X=0.5_SRK,Y=0.6_SRK)
  CALL point3%init(DIM=2,X=0.3_SRK,Y=0.2_SRK)
  point=point2+point3
  bool = .NOT.(point%dim /= 2 .OR. ANY(point%coord /= 0.8_SRK))
  ASSERT(bool, '2-D PointType OPERATOR(+)')
  CALL point2%clear()
  CALL point3%clear()

  CALL point2%init(DIM=3,X=0.5_SRK,Y=0.6_SRK,Z=0.7_SRK)
  CALL point3%init(DIM=3,X=0.3_SRK,Y=0.2_SRK,Z=0.1_SRK)
  point=point2+point3
  bool = .NOT.(point%dim /= 3 .OR. ANY(.NOT.(point%coord .APPROXEQ. 0.8_SRK)))
  ASSERT(bool, '3-D PointType OPERATOR(+)')
  CALL point2%clear()
  CALL point3%clear()

  CALL point2%init(DIM=4,COORD=(/0.5_SRK,0.6_SRK,0.7_SRK,0.8_SRK/))
  CALL point3%init(DIM=4,COORD=(/0.3_SRK,0.2_SRK,0.1_SRK,0.0_SRK/))
  point=point2+point3
  bool = .NOT.(point%dim /= 4 .OR. ANY(.NOT.(point%coord .APPROXEQ. 0.8_SRK)))
  ASSERT(bool, 'N-D PointType OPERATOR(+)')
  CALL point2%clear()
  CALL point3%clear()

  !Redundant calls for error checking
  point=point2+point3 !Empty case
  bool = .NOT.(point%dim /= 0 .OR. ALLOCATED(point%coord))
  ASSERT(bool, 'Empty PointType OPERATOR(+)')
  CALL point2%init(DIM=4,COORD=(/0.5_SRK,0.6_SRK,0.7_SRK,0.8_SRK/))
  point2=point2+point3 !mismatched dimensions case
  bool = .NOT.(point2%dim /= 0 .OR. ALLOCATED(point2%coord))
  ASSERT(bool, 'Mismatched PointType OPERATOR(+)')

  !Test subtraction
  COMPONENT_TEST('OPERATOR(-)')
  CALL point2%init(DIM=1,X=0.5_SRK)
  CALL point3%init(DIM=1,X=0.1_SRK)
  point=point2-point3
  bool = .NOT.(point%dim /= 1 .OR. point%coord(1) /= 0.4_SRK)
  ASSERT(bool, '1-D PointType OPERATOR(-)')
!      ENDIF
  CALL point2%clear()
  CALL point3%clear()

  CALL point2%init(DIM=2,X=0.5_SRK,Y=0.6_SRK)
  CALL point3%init(DIM=2,X=0.3_SRK,Y=0.2_SRK)
  point=point2-point3
  bool = .NOT.(point%dim /= 2 .OR. .NOT.(point%coord(1) .APPROXEQ. 0.2_SRK) .OR. &
      .NOT.(point%coord(2) .APPROXEQ. 0.4_SRK))
  ASSERT(bool, '2-D PointType OPERATOR(-)')
  CALL point2%clear()
  CALL point3%clear()

  CALL point2%init(DIM=3,X=0.5_SRK,Y=0.6_SRK,Z=0.7_SRK)
  CALL point3%init(DIM=3,X=0.3_SRK,Y=0.2_SRK,Z=0.1_SRK)
  point=point2-point3
  bool = .NOT.(point%dim /= 3 .OR. .NOT.(point%coord(1) .APPROXEQ. 0.2_SRK) .OR. &
      .NOT.(point%coord(2) .APPROXEQ. 0.4_SRK) .OR. &
      .NOT.(point%coord(3) .APPROXEQ. 0.6_SRK))
  ASSERT(bool, '3-D PointType OPERATOR(-)')
  CALL point2%clear()
  CALL point3%clear()

  CALL point2%init(DIM=4,COORD=(/0.5_SRK,0.6_SRK,0.7_SRK,0.8_SRK/))
  CALL point3%init(DIM=4,COORD=(/0.3_SRK,0.2_SRK,0.1_SRK,0.0_SRK/))
  point=point2-point3
  bool = .NOT.(point%dim /= 4 .OR. .NOT.(point%coord(1) .APPROXEQ. 0.2_SRK) .OR. &
      .NOT.(point%coord(2) .APPROXEQ. 0.4_SRK) .OR. &
      .NOT.(point%coord(3) .APPROXEQ. 0.6_SRK) .OR. &
      .NOT.(point%coord(4) .APPROXEQ. 0.8_SRK))
  ASSERT(bool, 'N-D PointType OPERATOR(-)')
  CALL point2%clear()
  CALL point3%clear()

  !Redundant calls for error checking
  point=point2-point3 !Empty case
  bool = .NOT.(point%dim /= 0 .OR. ALLOCATED(point%coord))
  ASSERT(bool, 'Empty PointType OPERATOR(-)')
  CALL point2%init(DIM=4,COORD=(/0.5_SRK,0.6_SRK,0.7_SRK,0.8_SRK/))
  point2=point2-point3 !mismatched dimensions case
  bool = .NOT.(point2%dim /= 0 .OR. ALLOCATED(point2%coord))
  ASSERT(bool, 'Mismatched PointType OPERATOR(-)')

  COMPONENT_TEST('OPERATOR(==)')
  CALL point2%init(DIM=4,COORD=(/0.5_SRK,0.6_SRK,0.7_SRK,0.8_SRK/))
  CALL point3%init(DIM=4,COORD=(/0.5_SRK,0.6_SRK,0.7_SRK,0.80000000000002_SRK/))
  point=point2
  bool = .NOT.(.NOT.(point == point2) .OR. point == point3)
  ASSERT(bool, 'PointType OPERATOR(==)')

  COMPONENT_TEST('OPERATOR(/=)')
  bool = .NOT.((point /= point2) .OR. .NOT.(point /= point3))
  ASSERT(bool, 'PointType OPERATOR(/=)')

  COMPONENT_TEST('OPERATOR(.APPROXEQ.)')
  point2%coord(4)=0.80000000000001_SRK
  bool = .NOT.(.NOT.(point .APPROXEQA. point2) .OR. (point .APPROXEQA. point3))
  ASSERT(bool, 'PointType OPERATOR(.APPROXEQ.)')
  CALL point2%clear()
  CALL point3%clear()

  CALL point%init(DIM=1,X=0.5_SRK)
  points=point
!
!Test clear routine
  COMPONENT_TEST('Elemental %clear()')
  CALL points%clear()
  bool = .NOT.(ANY(points%dim /= 0) .OR. ALLOCATED(points(1)%coord) .OR. &
      ALLOCATED(points(2)%coord))
  ASSERT(bool, 'point%clear()')
!
!Test distance routine
  COMPONENT_TEST('Elemental Distance()')
  CALL points(1)%init(DIM=3,X=0.5_SRK,Y=0.6_SRK,Z=0.7_SRK)
  CALL points2(1)%init(DIM=3,X=0.3_SRK,Y=0.2_SRK,Z=0.1_SRK)
  CALL points(2)%init(DIM=2,X=0.5_SRK,Y=0.6_SRK)
  CALL points2(2)%init(DIM=2,X=0.3_SRK,Y=0.2_SRK)
  s=Distance(points,points2)
  bool = .NOT.(.NOT.(s(1) .APPROXEQ. 0.748331477354788_SRK) .OR. &
      .NOT.(s(2) .APPROXEQ. 0.447213595499958_SRK))
  ASSERT(bool, 'Distance(...)')
!
!Test midPoint routine
  COMPONENT_TEST('Elemental Midpoint()')
  CALL points(1)%init(DIM=3,X=0.5_SRK,Y=0.6_SRK,Z=0.7_SRK)
  CALL points2(1)%init(DIM=3,X=0.3_SRK,Y=0.2_SRK,Z=0.1_SRK)
  CALL points(2)%init(DIM=2,X=0.5_SRK,Y=0.6_SRK)
  CALL points2(2)%init(DIM=2,X=0.3_SRK,Y=0.2_SRK)
  points3=midPoint(points,points2)
  bool = .NOT.(ANY(.NOT.(points3(1)%coord .APPROXEQ. 0.4_SRK)) .OR. &
      ANY(.NOT.(points3(2)%coord .APPROXEQ. 0.4_SRK)))
  ASSERT(bool, 'midPoint(...)')
!
!Test Operators
  !Addition
  COMPONENT_TEST('Elemental OPERATOR(+)')
  CALL points(1)%init(DIM=3,X=0.5_SRK,Y=0.6_SRK,Z=0.7_SRK)
  CALL points2(1)%init(DIM=3,X=0.3_SRK,Y=0.2_SRK,Z=0.1_SRK)
  CALL points(2)%init(DIM=2,X=0.5_SRK,Y=0.6_SRK)
  CALL points2(2)%init(DIM=2,X=0.3_SRK,Y=0.2_SRK)
  points3=points+points2
  bool = .NOT.(ANY(.NOT.(points3(1)%coord .APPROXEQ. 0.8_SRK)) .OR. &
      ANY(.NOT.(points3(2)%coord .APPROXEQ. 0.8_SRK)))
  ASSERT(bool, 'PointType Array OPERATOR(+)')

  !Subtraction
  COMPONENT_TEST('Elemental OPERATOR(-)')
  points3=points-points2
  bool = .NOT.(.NOT.(points3(1)%coord(1) .APPROXEQ. 0.2_SRK) .OR. &
      .NOT.(points3(1)%coord(2) .APPROXEQ. 0.4_SRK) .OR. &
      .NOT.(points3(1)%coord(3) .APPROXEQ. 0.6_SRK) .OR. &
      .NOT.(points3(2)%coord(1) .APPROXEQ. 0.2_SRK) .OR. &
      .NOT.(points3(2)%coord(2) .APPROXEQ. 0.4_SRK))
  ASSERT(bool, 'PointType Array OPERATOR(-)')

  !Equal to
  COMPONENT_TEST('Elemental OPERATOR(==)')
  points3=points
  bool = .NOT.(ANY(.NOT.(points == points3)) .OR. ANY(points2 == points))
  ASSERT(bool, 'PointType Array OPERATOR(==)')

  !Not equal to
  COMPONENT_TEST('Elemental OPERATOR(/=)')
  bool = .NOT.(ANY(.NOT.(points2 /= points)) .OR. ANY(points /= points3))
  ASSERT(bool, 'PointType Array OPERATOR(/=)')

  !Approximately equal to
  COMPONENT_TEST('Elemental OPERATOR(.APPROXEQ.)')
  points3=points
  points2=points
  points2(1)%coord(1)=0.50000000000002_SRK
  points2(2)%coord(1)=0.50000000000002_SRK
  points3(1)%coord(1)=0.50000000000001_SRK
  points3(2)%coord(1)=0.50000000000001_SRK
  bool = .NOT.(ANY(.NOT.(points .APPROXEQA. points3)) .OR. ANY(points2 .APPROXEQA. points))
  ASSERT(bool, 'PointType Array OPERATOR(.APPROXEQ.)')
!
!Test linked list of points
  COMPONENT_TEST('LinkedListPointType %insert()')
  ALLOCATE(firstPoint)
  thisPoint => firstPoint
  DO i=1,5
    CALL thisPoint%p%init(DIM=3,X=REAL(i,SRK),Y=0.1_SRK,Z=0.2_SRK)
    thisPoint%sortval=REAL(i,SRK)
    ALLOCATE(thisPoint%next)
    thisPoint => thisPoint%next
  ENDDO
  NULLIFY(thisPoint)
  CALL firstPoint%insert(thisPoint)
  ALLOCATE(thisPoint)
  thisPoint%sortval=-1._SRK
  CALL firstPoint%insert(thisPoint)
  thisPoint%sortval=101._SRK
  CALL firstPoint%insert(thisPoint)
  thisPoint%sortval=1._SRK
  CALL firstPoint%insert(thisPoint)
  CALL thisPoint%p%init(DIM=3,X=firstPoint%p%coord(1)+EPSREAL, &
      Y=firstPoint%p%coord(2)+EPSREAL,Z=firstPoint%p%coord(3)-EPSREAL)
  thisPoint%sortval=1._SRK+Distance(firstPoint%p,thisPoint%p)
  CALL firstPoint%insert(thisPoint)
  bool = ASSOCIATED(thisPoint)
  ASSERT(bool, 'firstPoint%insert(thisPoint) for nearby point FAILED')
  thisPoint%sortval=1._SRK-99._SRK*EPSREAL
  CALL firstPoint%insert(thisPoint)
  bool = ASSOCIATED(thisPoint)
  ASSERT(bool, 'firstPoint%insert(thisPoint) for fuzzy points < 1e-12')
  thisPoint%sortval=1._SRK+101._SRK*EPSREAL
  CALL firstPoint%insert(thisPoint)
  bool = .NOT.(ASSOCIATED(thisPoint))
  ASSERT(bool, 'firstPoint%insert(thisPoint) for fuzzy points > 1e-12')
  ALLOCATE(thisPoint)
  thisPoint%sortval=1.5_SRK
  CALL firstPoint%insert(thisPoint)
  bool = .NOT.(ASSOCIATED(thisPoint) .OR. firstPoint%next%next%sortval /= 1.5_SRK)
  ASSERT(bool, 'firstPoint%insert(thisPoint)')

  COMPONENT_TEST('LinkedListPointType %clear()')
  CALL ClearLinkedListPointType(firstPoint)
  ASSERT(.NOT.ASSOCIATED(firstPoint), 'ClearLinkedListPointType(...)')
ENDSUBROUTINE TestPoints

ENDPROGRAM testGeom_Points

max_stars_repo_name: wgurecky/Futility

max_issues_repo_name: wgurecky/Futility

max_forks_repo_name: wgurecky/Futility


text:
!
!
! =========================================================
      subroutine tfluct(maxmx,meqn,mwaves,mbc,mx,ql,qr,auxl,auxr, &
                 s,adq)
! =========================================================
!
!
      implicit double precision (a-h,o-z)
!
!
      write(*,*) "Error: method(8)=1, but you have not defined"
      write(*,*) "a function tfluct1()."
      stop
!
      return
      end

max_stars_repo_name: ketch/pyclaw

max_issues_repo_name: ketch/pyclaw

max_forks_repo_name: ketch/pyclaw