text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2015 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine sensitivitys(inpc,textpart,nmethod,iperturb,isolver,
     &  istep,istat,n,tinc,tper,tmin,tmax,idrct,iline,ipol,inl,ipoinp,
     &  inp,ithermal,cs,ics,tieset,istartset,
     &  iendset,ialset,ipompc,nodempc,coefmpc,nmpc,nmpc_,ikmpc,
     &  ilmpc,mpcfree,mcs,set,nset,labmpc,ipoinpc,iexpl,cfd,ttime,
     &  iaxial,nelcon,nmat,tincf)
!
!     reading the input deck: *SENSITIVITY
!
!     isolver=0: SPOOLES
!             2: iterative solver with diagonal scaling
!             3: iterative solver with Cholesky preconditioning
!             4: sgi solver
!             5: TAUCS
!             7: pardiso
!
!      iexpl==0:  structure:implicit, fluid:incompressible
!
      implicit none
!
      logical timereset
!
      character*1 inpc(*)
      character*20 labmpc(*),solver
      character*81 set(*),tieset(3,*)
      character*132 textpart(16)
!
      integer nmethod,iperturb,isolver,istep,istat,n,key,i,idrct,
     &  iline,ipol,inl,ipoinp(2,*),inp(3,*),ithermal,ics(*),iexpl,
     &  istartset(*),iendset(*),ialset(*),ipompc(*),nodempc(3,*),
     &  nmpc,nmpc_,ikmpc(*),ilmpc(*),mpcfree,nset,mcs,ipoinpc(0:*),
     &  cfd,iaxial,nelcon(2,*),nmat
!
      real*8 tinc,tper,tmin,tmax,cs(17,*),coefmpc(*),ttime,tincf
!
      idrct=0
      tper=1.d0
      tmin=0.d0
      tmax=0.d0
      timereset=.false.
!
      if(istep.lt.1) then
         write(*,*) '*ERROR reading *SENSITIVITY: *SENSITIVITY can
     &only be used within a STEP'
         call exit(201)
      endif
!
!     default solver
!
      solver='                    '
      if(isolver.eq.0) then
         solver(1:7)='SPOOLES'
      elseif(isolver.eq.2) then
         solver(1:16)='ITERATIVESCALING'
      elseif(isolver.eq.3) then
         solver(1:17)='ITERATIVECHOLESKY'
      elseif(isolver.eq.4) then
         solver(1:3)='SGI'
      elseif(isolver.eq.5) then
         solver(1:5)='TAUCS'
      elseif(isolver.eq.7) then
         solver(1:7)='PARDISO'
      endif
!
      do i=2,n
         if(textpart(i)(1:7).eq.'SOLVER=') then
            read(textpart(i)(8:27),'(a20)') solver
         elseif((textpart(i)(1:6).eq.'DIRECT').and.
     &          (textpart(i)(1:9).ne.'DIRECT=NO')) then
            idrct=1
         elseif(textpart(i)(1:9).eq.'TIMERESET') then
            timereset=.true.
         elseif(textpart(i)(1:17).eq.'TOTALTIMEATSTART=') then
            read(textpart(i)(18:37),'(f20.0)',iostat=istat) ttime
         else
            write(*,*)
     &        '*WARNING reading *SENSITIVITY: parameter not
     &recognized:'
            write(*,*) '         ',
     &                 textpart(i)(1:index(textpart(i),' ')-1)
            call inputwarning(inpc,ipoinpc,iline,
     &"*STATIC%")
         endif
      enddo
!
      if(solver(1:7).eq.'SPOOLES') then
         isolver=0
      elseif(solver(1:16).eq.'ITERATIVESCALING') then
         isolver=2
      elseif(solver(1:17).eq.'ITERATIVECHOLESKY') then
         isolver=3
      elseif(solver(1:3).eq.'SGI') then
         isolver=4
      elseif(solver(1:5).eq.'TAUCS') then
         isolver=5
      elseif(solver(1:7).eq.'PARDISO') then
         isolver=7
      else
         write(*,*) '*WARNING reading *SENSITIVITY: unknown solver;'
         write(*,*) '         the default solver is used'
      endif
!
      nmethod=12
!
      call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
     &        ipoinp,inp,ipoinpc)
!
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
      MODULE GWFSPGMODULE
        INTEGER,SAVE,POINTER  ::NSEEP,MXSPG,NSPGVL,ISPGCB,IPRSPG
        INTEGER,SAVE,POINTER  ::NPSPG,ISPGPB,NNPSPG
        CHARACTER(LEN=16),SAVE, DIMENSION(:),     POINTER   ::SPGAUX
        REAL,             SAVE, DIMENSION(:,:),   POINTER   ::SEEP
        REAL,             SAVE, DIMENSION(:,:,:), POINTER   ::QSPG
      TYPE GWFSPGTYPE
        INTEGER,POINTER  ::NSEEP,MXSPG,NSPGVL,ISPGCB,IPRSPG
        INTEGER,POINTER  ::NPSPG,ISPGPB,NNPSPG
        CHARACTER(LEN=16), DIMENSION(:),     POINTER     ::SPGAUX
        REAL,              DIMENSION(:,:),   POINTER     ::SEEP
        REAL,              DIMENSION(:,:,:), POINTER     ::QSPG
      END TYPE
      TYPE(GWFSPGTYPE), SAVE:: GWFSPGDAT(10)
      END MODULE GWFSPGMODULE

      SUBROUTINE GWF2SPG7AR(IN,IGRID)
C     ******************************************************************
C     ALLOCATE ARRAY STORAGE FOR SEEPS AND READ PARAMETER DEFINITIONS
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,       ONLY:IOUT,NCOL,NROW,NLAY,IFREFM
      USE GWFSPGMODULE, ONLY:NSEEP,MXSPG,NSPGVL,ISPGCB,IPRSPG,NPSPG,
     1                       ISPGPB,NNPSPG,SPGAUX,SEEP,QSPG
      CHARACTER*200 LINE
C     ------------------------------------------------------------------
      ALLOCATE(NSEEP,MXSPG,NSPGVL,ISPGCB,IPRSPG)
      ALLOCATE(NPSPG,ISPGPB,NNPSPG)
C
C1------IDENTIFY PACKAGE AND INITIALIZE NSEEP.
      WRITE(IOUT,1)IN
    1 FORMAT(1X,/1X,'SPG -- SEEPAGE PACKAGE, VERSION 7, 7/4/2021',
     1' INPUT READ FROM UNIT ',I4)
      NSEEP=0
      NNPSPG=0
C
C2------READ MAXIMUM NUMBER OF SEEPS AND UNIT OR FLAG FOR
C2------CELL-BY-CELL FLOW TERMS.
      CALL URDCOM(IN,IOUT,LINE)
      CALL UPARLSTAL(IN,IOUT,LINE,NPSPG,MXPD)
      IF(IFREFM.EQ.0) THEN
         READ(LINE,'(2I10)') MXACTS,ISPGCB
         LLOC=21
      ELSE
         LLOC=1
         CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,MXACTS,R,IOUT,IN)
         CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISPGCB,R,IOUT,IN)
      END IF
      WRITE(IOUT,3) MXACTS
    3 FORMAT(1X,'MAXIMUM OF ',I6,' ACTIVE SEEPAGES AT ONE TIME')
      IF(ISPGCB.LT.0) WRITE(IOUT,7)
    7 FORMAT(1X,'CELL-BY-CELL FLOWS WILL BE PRINTED WHEN ICBCFL NOT 0')
         IF(ISPGCB.GT.0) WRITE(IOUT,8) ISPGCB
    8 FORMAT(1X,'CELL-BY-CELL FLOWS WILL BE SAVED ON UNIT ',I4)
C
C3------READ AUXILIARY VARIABLES AND CBC ALLOCATION OPTION.
      ALLOCATE (SPGAUX(20))
      NAUX=0
      IPRSPG=1
   10 CALL URWORD(LINE,LLOC,ISTART,ISTOP,1,N,R,IOUT,IN)
      IF(LINE(ISTART:ISTOP).EQ.'AUXILIARY' .OR.
     1        LINE(ISTART:ISTOP).EQ.'AUX') THEN
         CALL URWORD(LINE,LLOC,ISTART,ISTOP,1,N,R,IOUT,IN)
         IF(NAUX.LT.20) THEN
            NAUX=NAUX+1
            SPGAUX(NAUX)=LINE(ISTART:ISTOP)
            WRITE(IOUT,12) SPGAUX(NAUX)
   12       FORMAT(1X,'AUXILIARY SEEPAGE VARIABLE: ',A)
         END IF
         GO TO 10
      ELSE IF(LINE(ISTART:ISTOP).EQ.'NOPRINT') THEN
         WRITE(IOUT,13)
   13    FORMAT(1X,'LISTS OF SEEPAGE CELLS WILL NOT BE PRINTED')
         IPRSPG = 0
         GO TO 10
      END IF
C3A-----THERE ARE FIVE INPUT DATA VALUES PLUS ONE LOCATION FOR
C3A-----CELL-BY-CELL FLOW.
      NSPGVL=6+NAUX
C
C4------ALLOCATE SPACE FOR SEEPAGE ARRAYs.
      ISPGPB=MXACTS+1
      MXSPG=MXACTS+MXPD
      ALLOCATE (SEEP(NSPGVL,MXSPG))
      ALLOCATE (QSPG(NCOL,NROW,NLAY))
C
C5------READ NAMED PARAMETERS.
      WRITE(IOUT,1000) NPSPG
 1000 FORMAT(1X,//1X,I5,' Seepage parameters')
      IF(NPSPG.GT.0) THEN
        LSTSUM=ISPGPB
        DO 120 K=1,NPSPG
          LSTBEG=LSTSUM
          CALL UPARLSTRP(LSTSUM,MXSPG,IN,IOUT,IP,'SPG','SPG',1,
     &                   NUMINST)
          NLST=LSTSUM-LSTBEG
          IF(NUMINST.EQ.0) THEN
C5A-----READ PARAMETER WITHOUT INSTANCES
            CALL ULSTRD(NLST,SEEP,LSTBEG,NSPGVL,MXSPG,1,IN,IOUT,
     &      'SEEP. NO.  LAYER   ROW   COL     SEEP. EL.  STRESS FACTOR',
     &        SPGAUX,20,NAUX,IFREFM,NCOL,NROW,NLAY,5,5,IPRSPG)
          ELSE
C5B-----READ INSTANCES
            NINLST=NLST/NUMINST
            DO 110 I=1,NUMINST
            CALL UINSRP(I,IN,IOUT,IP,IPRSPG)
            CALL ULSTRD(NINLST,SEEP,LSTBEG,NSPGVL,MXSPG,1,IN,IOUT,
     &      'SEEP. NO.  LAYER   ROW   COL     SEEP. EL.  STRESS FACTOR',
     &        SPGAUX,20,NAUX,IFREFM,NCOL,NROW,NLAY,5,5,IPRSPG)
            LSTBEG=LSTBEG+NINLST
  110       CONTINUE
          END IF
  120   CONTINUE
      END IF
C
C6------RETURN
      CALL SGWF2SPG7PSV(IGRID)
      RETURN
      END
      SUBROUTINE GWF2SPG7RP(IN,IGRID)
C     ******************************************************************
C     READ SEEPAGE HEAD, CONDUCTANCE AND BOTTOM ELEVATION
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,       ONLY:IOUT,NCOL,NROW,NLAY,IFREFM
      USE GWFSPGMODULE, ONLY:NSEEP,MXSPG,NSPGVL,IPRSPG,NPSPG,
     1                       ISPGPB,NNPSPG,SPGAUX,SEEP
C     ------------------------------------------------------------------
      CALL SGWF2SPG7PNT(IGRID)
C
C1------READ ITMP (NUMBER OF SEEPAGES OR FLAG TO REUSE DATA) AND
C1------NUMBER OF PARAMETERS.
      IF(NPSPG.GT.0) THEN
         IF(IFREFM.EQ.0) THEN
            READ(IN,'(2I10)') ITMP,NP
         ELSE
            READ(IN,*) ITMP,NP
         END IF
      ELSE
         NP=0
         IF(IFREFM.EQ.0) THEN
            READ(IN,'(I10)') ITMP
         ELSE
            READ(IN,*) ITMP
         END IF
      END IF
C
C------CALCULATE SOME CONSTANTS
      NAUX=NSPGVL-6
      IOUTU = IOUT
      IF(IPRSPG.EQ.0) IOUTU=-IOUT
C
C2------DETERMINE THE NUMBER OF NON-PARAMETER SEEPAGES.
      IF(ITMP.LT.0) THEN
         WRITE(IOUT,7)
    7    FORMAT(1X,/1X,
     1        'REUSING NON-PARAMETER SEEPAGES FROM LAST STRESS PERIOD')
      ELSE
         NNPSPG=ITMP
      END IF
C
C3------IF THERE ARE NEW NON-PARAMETER SEEPAGES, READ THEM.
      MXACTS=ISPGPB-1
      IF(ITMP.GT.0) THEN
         IF(NNPSPG.GT.MXACTS) THEN
            WRITE(IOUT,99) NNPSPG,MXACTS
   99       FORMAT(1X,/1X,'THE NUMBER OF ACTIVE SEEPAGES (',I6,
     1                     ') IS GREATER THAN MXACTS(',I6,')')
            CALL USTOP(' ')
         END IF
         CALL ULSTRD(NNPSPG,SEEP,1,NSPGVL,MXSPG,1,IN,IOUT,
     1     'SEEP. NO.  LAYER   ROW   COL     SEEP. EL.  CONDUCTANCE',
     2     SPGAUX,20,NAUX,IFREFM,NCOL,NROW,NLAY,5,5,IPRSPG)
      END IF
      NSEEP=NNPSPG
C
C1C-----IF THERE ARE ACTIVE SPG PARAMETERS, READ THEM AND SUBSTITUTE
      CALL PRESET('SPG')
      IF(NP.GT.0) THEN
         NREAD=NSPGVL-1
         DO 30 N=1,NP
         CALL UPARLSTSUB(IN,'SPG',IOUTU,'SPG',SEEP,NSPGVL,MXSPG,NREAD,
     1                MXACTS,NSEEP,5,5,
     2     'SEEP. NO.  LAYER   ROW   COL     SEEP. EL.  CONDUCTANCE',
     3            SPGAUX,20,NAUX)
   30    CONTINUE
      END IF
C
C3------PRINT NUMBER OF SEEPAGES IN CURRENT STRESS PERIOD.
      WRITE (IOUT,101) NSEEP
  101 FORMAT(1X,/1X,I6,' SEEPAGES')
C
C8------RETURN.
  260 RETURN
      END
      SUBROUTINE GWF2SPG7FM(KPER,IGRID)
C     ******************************************************************
C     SETS GROUNDWATER LEVEL EQUAL TO SEEPAGE LEVEL IF CELL IS SEEPING
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,       ONLY:NCOL,NROW,NLAY,IBOUND,HNEW,CV,CC,CR,IUNIT
      USE GWFSPGMODULE, ONLY:NSEEP,SEEP,QSPG
C
      DOUBLE PRECISION HHNEW,EEL
C     ------------------------------------------------------------------
      CALL SGWF2SPG7PNT(IGRID)
C
      ZERO=0.0
C1------INITIALIZE QSPG.
      DO 50 IL=1,NLAY
      DO 50 IR=1,NROW
      DO 50 IC=1,NCOL
      QSPG(IC,IR,IL)=ZERO
 50   CONTINUE
C
C2------IF NSEEP<=0 THERE ARE NO SEEPAGES. RETURN.
      IF(NSEEP.LE.0) RETURN
C
C3------COLLECT FLUXES IN EACH CELL, FROM ALL ACTIVE STRESSES.
C
C3B-----ADD RECHARGE FLUX TO QSPG.
      IF(IUNIT(8).GT.0) CALL SGWF2SPG7H()
C
C3L-----ADD BCF FLUX TO QSPG.
      IF(IUNIT(1).GT.0) CALL SGWF2SPG7R(KPER)
C
C4------PROCESS EACH CELL IN THE SEEPAGE LIST TO CHECK FLOW.
C4------CELL DISCHARGES GROUNDWATER IF Q<0,
C4------IS INTERMEDIATE IF STRESS FLUXES>Q>0 AND
C4------IS INFILTRATING IF Q>STRESS FLUXES.
C4------THEREFOR SUBSTRACT EXTERNAL (NON SEEPAGE)
C4------FLUXES, QSPG, FROM Q.
C4------NOTE WE CHECK FLOW FOR ALL SEEPAGE CELLS BEFORE CHECKING
C4------HEAD (C5) TO AVOID USING NEW HEAD WHEN CHECKING FLOW.
      DO 101 L=1,NSEEP
C
C4A-----GET COLUMN, ROW AND LAYER OF CELL CONTAINING SEEPAGE.
      IL=SEEP(1,L)
      IR=SEEP(2,L)
      IC=SEEP(3,L)
C
C4B-----IF SEEPAGE IS NOT ACTIVE SKIP CELL.
      IF(IBOUND(IC,IR,IL).NE.-777) GOTO 101
C
C4C-----GET HEAD AT THE CELL.
      HHNEW=HNEW(IC,IR,IL)
C
C4D-----CALCULATE FIRST VERTICAL FLUX TO SEEPAGE CELL
C4D-----SEEPAGE CELL IS IN A ONE-LAYER SYSTEM.
      FLUXV=ZERO
      IF(NLAY.GT.1) THEN
C
C4D1----SEEPAGE CELL IS IN FIRST LAYER OF A MULTI-LAYER SYSTEM.
         IF(IL.EQ.1) THEN
            FLUXV=CV(IC,IR,IL)*HHNEW-CV(IC,IR,IL)*HNEW(IC,IR,IL+1)
C
C4D2----SEEPAGE CELL IS BETWEEN FIRST LAYER AND
C4D3----BOTTOM LAYER OF A MULTI-LAYER SYSTEM.
         ELSE IF(IL.LT.NLAY.AND.IL.GT.1) THEN
            FLUXV=CV(IC,IR,IL-1)*HHNEW-
     1           CV(IC,IR,IL-1)*HNEW(IC,IR,IL-1) +
     1           CV(IC,IR,IL)*HHNEW-CV(IC,IR,IL)*HNEW(IC,IR,IL+1)
C
C4D4----SEEPAGE CELL IS IN LOWEST LAYER OF A MULTI-LAYER SYSTEM.
         ELSE
            FLUXV=CV(IC,IR,IL-1)*HHNEW-
     1           CV(IC,IR,IL-1)*HNEW(IC,IR,IL-1)
         END IF
      ENDIF
C
C4E-----FORMULATE FLUX TO SEEPAGE CELL.
      Q=FLUXV
      IF(IR.LT.NROW) Q=Q+CC(IC,IR,IL)*HHNEW-
     1     CC(IC,IR,IL)*HNEW(IC,IR+1,IL)
      IF(IR.GT.1) Q=Q+CC(IC,IR-1,IL)*HHNEW-
     1     CC(IC,IR-1,IL)*HNEW(IC,IR-1,IL)
      IF(IC.LT.NCOL) Q=Q+CR(IC,IR,IL)*HHNEW-
     1     CR(IC,IR,IL)*HNEW(IC+1,IR,IL)
      IF(IC.GT.1) Q=Q+CR(IC-1,IR,IL)*HHNEW-
     1     CR(IC-1,IR,IL)*HNEW(IC-1,IR,IL)
      Q=Q+QSPG(IC,IR,IL)
C
C4F-----IF FLUX IS LARGER THAN ZERO CELL MUST BE TURNED BACK TO INTERNAL.
C4F-----HERE WE ONLY TAG IT TO AVOID DEALING WITH IT IN C5.
      IF(Q.GT.ZERO) IBOUND(IC,IR,IL)=-778
 101  CONTINUE

C5------PROCESS EACH CELL IN THE SEEPAGE LIST TO CHECK HEAD.
      DO 102 L=1,NSEEP
C
C5A-----GET COLUMN, ROW AND LAYER OF CELL CONTAINING SEEPAGE.
      IL=SEEP(1,L)
      IR=SEEP(2,L)
      IC=SEEP(3,L)
C
C5B-----GET SEEPAGE PARAMETERS FROM SEEPAGE LIST.
      EL=SEEP(4,L)
      EEL=EL
C
C5C-----GET HEAD AT THE CELL.
      HHNEW=HNEW(IC,IR,IL)
C
C5D-----IF CELL IS INTERNAL AND HEAD HIGHER THAN SEEPAGE LEVEL
C5D-----TURN IT TO CONSTANT HEAD.
      IF(IBOUND(IC,IR,IL).GT.0 .AND. HHNEW.GT.EEL) THEN
         IBOUND(IC,IR,IL)=-777
         HNEW(IC,IR,IL)=EEL
C5E-----TURN PREVIOUSLY TAGGED CELL (C4F) INTO INTERNAL CELL.
      ELSEIF(IBOUND(IC,IR,IL).EQ.-778) THEN
         IBOUND(IC,IR,IL)=1
      ENDIF
 102  CONTINUE
C
C16-----RETURN.
      RETURN
      END
      SUBROUTINE GWF2SPG7BD(KSTP,KPER,IGRID)
C     ******************************************************************
C     CALCULATE VOLUMETRIC BUDGET FOR SEEPAGES
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,      ONLY:IOUT,NCOL,NROW,NLAY,IBOUND,HNEW,BUFF,CV,CC,
     1                      CR,IUNIT
      USE GWFBASMODULE,ONLY:MSUM,ICBCFL,IAUXSV,DELT,PERTIM,TOTIM,
     1                      VBVL,VBNM
      USE GWFSPGMODULE,ONLY:NSEEP,ISPGCB,SEEP,NSPGVL,SPGAUX,QSPG
C
      CHARACTER*16 TEXT
      DOUBLE PRECISION HHNEW,EEL,CCSPG,CEL,RATOUT,QQ
C
      DATA TEXT /'        SEEPAGES'/
C     ------------------------------------------------------------------
      CALL SGWF2SPG7PNT(IGRID)
C
C1------INITIALIZE CELL-BY-CELL FLOW TERM FLAG (IBD) AND
C1------ACCUMULATOR (RATOUT).
      ZERO=0.
      RATOUT=ZERO
      IBD=0
      IF(ISPGCB.LT.0 .AND. ICBCFL.NE.0) IBD=-1
      IF(ISPGCB.GT.0) IBD=ICBCFL
      IBDLBL=0
C
C2------IF CELL-BY-CELL FLOWS WILL BE SAVED AS A LIST, WRITE HEADER.
      IF(IBD.EQ.2) THEN
         NAUX=NSPGVL-6
         IF(IAUXSV.EQ.0) NAUX=0
         CALL UBDSV4(KSTP,KPER,TEXT,NAUX,SPGAUX,ISPGCB,NCOL,NROW,NLAY,
     1          NSEEP,IOUT,DELT,PERTIM,TOTIM,IBOUND)
      END IF
C
C3------CLEAR THE BUFFER AND QSPG.
      DO 50 IL=1,NLAY
      DO 50 IR=1,NROW
      DO 50 IC=1,NCOL
      BUFF(IC,IR,IL)=ZERO
      QSPG(IC,IR,IL)=ZERO
50    CONTINUE
C
C4------IF THERE ARE NO SEEPAGES THEN DO NOT ACCUMULATE SEEPAGE FLOW.
      IF(NSEEP.LE.0) GO TO 200
C
C5------COLLECT FLUXES IN EACH CELL, FROM ALL ACTIVE STRESSES.
C
C5B-----ADD RECHARGE FLUX TO QSPG.
      IF(IUNIT(8).GT.0) CALL SGWF2SPG7H()
C
C5L-----ADD BCF FLUX TO QSPG.
      IF(IUNIT(1).GT.0) CALL SGWF2SPG7R(KPER)
C
C6------LOOP THROUGH EACH SEEPAGE CALCULATING FLOW.
      DO 100 L=1,NSEEP
      Q=ZERO
C
C6A-----GET LAYER, ROW & COLUMN OF CELL CONTAINING SEEPAGE.
      IL=SEEP(1,L)
      IR=SEEP(2,L)
      IC=SEEP(3,L)
C
C6B-----IF SEEPAGE IS NOT ACTIVE SKIP CALCULATION.
      IF(IBOUND(IC,IR,IL).NE.-777) GO TO 60
C
C6C-----SET TEMPORARY CONSTANT HEAD CELL BACK TO INTERNAL CELL.
      IBOUND(IC,IR,IL)=1
C
C6C-----GET HEAD AT THE CELL.
      HHNEW=HNEW(IC,IR,IL)
C
C7A----SEEPAGE CELL IS IN A ONE-LAYER SYSTEM.
      FLUXV=ZERO
      IF(NLAY.GT.1) THEN
C
C7B----SEEPAGE CELL IS IN FIRST LAYER OF A MULTI-LAYER SYSTEM.
         IF(IL.EQ.1) THEN
            FLUXV=CV(IC,IR,IL)*HHNEW-CV(IC,IR,IL)*HNEW(IC,IR,IL+1)
C
C7C----SEEPAGE CELL IS BETWEEN FIRST LAYER AND
C7C----BOTTOM LAYER OF A MULTI-LAYER SYSTEM.
         ELSE IF(IL.LT.NLAY.AND.IL.GT.1) THEN
            FLUXV=CV(IC,IR,IL-1)*HHNEW-
     1           CV(IC,IR,IL-1)*HNEW(IC,IR,IL-1) +
     1           CV(IC,IR,IL)*HHNEW-CV(IC,IR,IL)*HNEW(IC,IR,IL+1)
C
C7D----SEEPAGE CELL IS IN LOWEST LAYER OF A MULTI-LAYER SYSTEM.
         ELSE
            FLUXV=CV(IC,IR,IL-1)*HHNEW-
     1           CV(IC,IR,IL-1)*HNEW(IC,IR,IL-1)
         END IF
      ENDIF
C
C7E----FORMULATE FLUX TO SEEPAGE CELL.
      Q=FLUXV
      IF(IR.LT.NROW) Q=Q+CC(IC,IR,IL)*HHNEW-
     1     CC(IC,IR,IL)*HNEW(IC,IR+1,IL)
      IF(IR.GT.1) Q=Q+CC(IC,IR-1,IL)*HHNEW-
     1     CC(IC,IR-1,IL)*HNEW(IC,IR-1,IL)
      IF(IC.LT.NCOL) Q=Q+CR(IC,IR,IL)*HHNEW-
     1     CR(IC,IR,IL)*HNEW(IC+1,IR,IL)
      IF(IC.GT.1) Q=Q+CR(IC-1,IR,IL)*HHNEW-
     1     CR(IC-1,IR,IL)*HNEW(IC-1,IR,IL)
      Q=Q+QSPG(IC,IR,IL)
C
C8----ACCUMULATE FLOW OUT OF CELL IF THERE IS GROUNDWATER DISCHARGE.
      IF(Q.LT.ZERO) THEN
         RATOUT=RATOUT-Q
      ELSE
         WRITE(IOUT,59) Q,IL,IR,IC
 59      FORMAT(1X,'WARNING THERE IS AN UNACCOUNTED FLUX OF (',
     1        1PG15.6,') IN SEEPAGE CELL (',I3,',',I3,',',I3,').',
     2        'MAYBE YOUR CONVERGENCE CRITERIA IS TOO BIG (CHECK '
     3        'TOTAL WATER BALANCE) OR YOU USE AN UNSUPPORTED ',
     4        'PACKAGE.')
      ENDIF
C
C9-----PRINT THE INDIVIDUAL RATES IF REQUESTED(ISPGCB<0).
 60   IF(IBD.LT.0) THEN
         IF(IBDLBL.EQ.0) WRITE(IOUT,61) TEXT,KPER,KSTP
   61    FORMAT(1X,/1X,A,'   PERIOD ',I4,'   STEP ',I3)
         WRITE(IOUT,62) L,IL,IR,IC,Q
   62    FORMAT(1X,'SEEPAGE ',I6,'   LAYER ',I3,'   ROW ',I5,
     1       '   COL ',I5,'   RATE ',1PG15.6)
         IBDLBL=1
      END IF
C
C10-----ADD Q TO BUFFER.
      BUFF(IC,IR,IL)=BUFF(IC,IR,IL)+Q
C
C11-----IF SAVING CELL-BY-CELL FLOWS IN A LIST, WRITE FLOW.  ALSO
C11-----COPY FLOW TO SEEP.
   99 IF(IBD.EQ.2) CALL UBDSVB(ISPGCB,NCOL,NROW,IC,IR,IL,Q,
     1                  SEEP(:,L),NSPGVL,NAUX,6,IBOUND,NLAY)
      SEEP(NSPGVL,L)=Q
  100 CONTINUE
C
C12------IF CELL-BY-CELL FLOW WILL BE SAVED AS A 3-D ARRAY,
C12------CALL UBUDSV TO SAVE THEM.
      IF(IBD.EQ.1) CALL UBUDSV(KSTP,KPER,TEXT,ISPGCB,BUFF,NCOL,NROW,
     1                          NLAY,IOUT)
C
C13------MOVE RATES,VOLUMES & LABELS INTO ARRAYS FOR PRINTING.
  200 ROUT=RATOUT
      VBVL(3,MSUM)=ZERO
      VBVL(4,MSUM)=ROUT
      VBVL(2,MSUM)=VBVL(2,MSUM)+ROUT*DELT
      VBNM(MSUM)=TEXT
C
C14------INCREMENT BUDGET TERM COUNTER.
      MSUM=MSUM+1
C
C15------RETURN.
      RETURN
      END
      SUBROUTINE SGWF2SPG7H()
C     ******************************************************************
C     ADD RECHARGE TO QSPG
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,      ONLY:NCOL,NROW,NLAY,IBOUND
      USE GWFSPGMODULE,ONLY:QSPG
      USE GWFRCHMODULE,ONLY:NRCHOP,IRCH,RECH
C     ------------------------------------------------------------------
C
C1------IF NRCHOP IS 1 RECHARGE IS IN TOP LAYER.
      IF(NRCHOP.NE.1) GO TO 15
      DO 10 IR=1,NROW
      DO 10 IC=1,NCOL
C
C2A-----IF CELL IS EXTERNAL THERE IS NO RECHARGE INTO IT.
         IF(IBOUND(IC,IR,1).LE.0 .AND. IBOUND(IC,IR,1).NE.-777) GO TO 10
C
C2B-----ADD RECHARGE RATE TO QSPG.
         QSPG(IC,IR,1)=QSPG(IC,IR,1)-RECH(IC,IR)
 10   CONTINUE
      GO TO 100
C
C3------IF OPTION IS 2 THEN RECHARGE IS IN INDICATOR ARRAY.
   15 IF(NRCHOP.NE.2) GO TO 25
      DO 20 IR=1,NROW
      DO 20 IC=1,NCOL
C
C3A-----LAYER INDEX IS IN INDICATOR ARRAY.
         IL=IRCH(IC,IR)
C
C3B-----IF THE CELL IS EXTERNAL THERE IS NO RECHARGE INTO IT.
         IF(IBOUND(IC,IR,IL).LE.0 .AND. IBOUND(IC,IR,IL).NE.-777)
     1      GO TO 20
C
C3C-----ADD RECHARGE RATE TO QSPG.
         QSPG(IC,IR,IL)=QSPG(IC,IR,IL)-RECH(IC,IR)
 20   CONTINUE
      GO TO 100
C
C4------IF OPTION IS 3 RECHARGE IS INTO HIGHEST INTERNAL CELL.
   25 IF(NRCHOP.NE.3) GO TO 100
      DO 30 IR=1,NROW
      DO 30 IC=1,NCOL
      DO 28 IL=1,NLAY
C
C4A-----IF CELL IS CONSTANT HEAD MOVE ON TO NEXT HORIZONTAL LOCATION.
         IF(IBOUND(IC,IR,IL).LT.0 .AND. IBOUND(IC,IR,IL).NE.-777)
     1      GO TO 30
C
C4B-----IF CELL IS INACTIVE MOVE DOWN A LAYER.
         IF (IBOUND(IC,IR,IL).EQ.0) GO TO 28
C
C4C-----SET LAYER INDEX ARRAY EQUAL TO ACTIVE LAYER NUMBER.
         QSPG(IC,IR,IL)=QSPG(IC,IR,IL)-RECH(IC,IR)
         GO TO 30
 28   CONTINUE
 30   CONTINUE
 100  CONTINUE
C
C5------RETURN
      RETURN
      END
      SUBROUTINE SGWF2SPG7R(KPER)
C     ******************************************************************
C     ADD STORAGE AND LEAKAGE CORRECTION TO QSPG
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,      ONLY:NCOL,NROW,NLAY,IBOUND,HNEW,HOLD,CV,BOTM,
     1     LBOTM,ISSFLG
      USE GWFSPGMODULE,ONLY:QSPG
      USE GWFBASMODULE,ONLY:DELT
      USE GWFBCFMODULE,ONLY:LAYCON,SC1,SC2
C     ------------------------------------------------------------------
C
      ISS=ISSFLG(KPER)
C
      KT=0
      ONE=1.
      TLED=ONE/DELT
C
C1------IF THE SIMULATION IS TRANSIENT DETERMINE FLUX FROM STORAGE
      IF(ISS.NE.0) GO TO 201
      KT=0
      DO 200 K=1,NLAY
C
C2------SEE IF THIS LAYER IS CONVERTIBLE OR NON-CONVERTIBLE.
      IF(LAYCON(K).EQ.3 .OR. LAYCON(K).EQ.2) GO TO 150
C3------NON-CONVERTIBLE LAYER, SO USE PRIMARY STORAGE.
      DO 140 I=1,NROW
      DO 140 J=1,NCOL
      IF(IBOUND(J,I,K).LE.0 .AND. IBOUND(J,I,K).NE.-777) GO TO 140
      RHO=SC1(J,I,K)*TLED
      QSPG(J,I,K)=QSPG(J,I,K)-(RHO*HOLD(J,I,K)-RHO*HNEW(J,I,K))
  140 CONTINUE
      GO TO 200
C
C4------A CONVERTIBLE LAYER, SO CHECK OLD AND NEW HEADS TO DETERMINE
C4------WHEN TO USE PRIMARY AND SECONDARY STORAGE
  150 KT=KT+1
      DO 180 I=1,NROW
      DO 180 J=1,NCOL
C
C5A-----IF THE CELL IS EXTERNAL THEN SKIP IT.
      IF(IBOUND(J,I,K).LE.0 .AND. IBOUND(J,I,K).NE.-777) GO TO 180
      TP=BOTM(J,I,LBOTM(K)-1)
      RHO2=SC2(J,I,KT)*TLED
      RHO1=SC1(J,I,K)*TLED
C
C5B-----FIND STORAGE FACTOR AT START OF TIME STEP.
      SOLD=RHO2
      IF(HOLD(J,I,K).GT.TP) SOLD=RHO1
C
C5C-----FIND STORAGE FACTOR AT END OF TIME STEP.
      HTMP=HNEW(J,I,K)
      SNEW=RHO2
      IF(HTMP.GT.TP) SNEW=RHO1
C
C5D-----ADD FLUX TO QSPG.
      QSPG(J,I,K)=QSPG(J,I,K)-(SOLD*(HOLD(J,I,K)-TP)-SNEW*(HTMP-TP))
C
  180 CONTINUE
C
  200 CONTINUE
C
C6------FOR EACH LAYER DETERMINE IF CORRECTION TERMS ARE NEEDED FOR
C6------FLOW DOWN INTO PARTIALLY SATURATED LAYERS.
  201 DO 300 K=1,NLAY
C
C7------SEE IF CORRECTION IS NEEDED FOR LEAKAGE FROM ABOVE.
      IF(LAYCON(K).NE.3 .AND. LAYCON(K).NE.2) GO TO 250
      IF(K.EQ.1) GO TO 250
C
C7A-----FOR EACH CELL MAKE THE CORRECTION IF NEEDED.
      DO 220 I=1,NROW
      DO 220 J=1,NCOL
C
C7B-----IF THE CELL IS EXTERNAL(IBOUND<=0) THEN SKIP IT.
      IF(IBOUND(J,I,K).LE.0 .AND. IBOUND(J,I,K).NE.-777) GO TO 220
      HTMP=HNEW(J,I,K)
C
C7C-----IF HEAD IS ABOVE TOP THEN CORRECTION NOT NEEDED
      IF(HTMP.GE.BOTM(J,I,LBOTM(K)-1)) GO TO 220
C
C7D-----WITH HEAD BELOW TOP ADD CORRECTION TERMS TO QSPG.
      QSPG(J,I,K)=QSPG(J,I,K)+CV(J,I,K-1)*(BOTM(J,I,LBOTM(K)-1)-HTMP)
  220 CONTINUE
C
C8------SEE IF THIS LAYER MAY NEED CORRECTION FOR LEAKAGE TO BELOW.
  250 IF(K.EQ.NLAY) GO TO 300
      IF(LAYCON(K+1).NE.3 .AND. LAYCON(K+1).NE.2) GO TO 300
C
C8A-----FOR EACH CELL MAKE THE CORRECTION IF NEEDED.
      DO 280 I=1,NROW
      DO 280 J=1,NCOL
C
C8B-----IF CELL IS EXTERNAL (IBOUND<=0) THEN SKIP IT.
      IF(IBOUND(J,I,K).LE.0 .AND. IBOUND(J,I,K).NE.-777) GO TO 280
C
C8C-----IF HEAD IN THE LOWER CELL IS LESS THAN TOP ADD CORRECTION
C8C-----TERM TO QSPG.
      HTMP=HNEW(J,I,K+1)
      IF(HTMP.LT.BOTM(J,I,LBOTM(K+1)-1)) QSPG(J,I,K)=QSPG(J,I,K)
     1     -CV(J,I,K)*(BOTM(J,I,LBOTM(K+1)-1)-HTMP)
  280 CONTINUE
  300 CONTINUE
C
C9------RETURN
      RETURN
      END
      SUBROUTINE GWF2SPG7DA(IGRID)
C  Deallocate SPG MEMORY
      USE GWFSPGMODULE
C
        CALL SGWF2SPG7PNT(IGRID)
        DEALLOCATE(NSEEP)
        DEALLOCATE(MXSPG)
        DEALLOCATE(NSPGVL)
        DEALLOCATE(ISPGCB)
        DEALLOCATE(IPRSPG)
        DEALLOCATE(NPSPG)
        DEALLOCATE(ISPGPB)
        DEALLOCATE(NNPSPG)
        DEALLOCATE(SPGAUX)
        DEALLOCATE(SEEP)
        DEALLOCATE(QSPG)
C
      RETURN
      END
      SUBROUTINE SGWF2SPG7PNT(IGRID)
C  Change SPG data to a different grid.
      USE GWFSPGMODULE
C
        NSEEP=>GWFSPGDAT(IGRID)%NSEEP
        MXSPG=>GWFSPGDAT(IGRID)%MXSPG
        NSPGVL=>GWFSPGDAT(IGRID)%NSPGVL
        ISPGCB=>GWFSPGDAT(IGRID)%ISPGCB
        IPRSPG=>GWFSPGDAT(IGRID)%IPRSPG
        NPSPG=>GWFSPGDAT(IGRID)%NPSPG
        ISPGPB=>GWFSPGDAT(IGRID)%ISPGPB
        NNPSPG=>GWFSPGDAT(IGRID)%NNPSPG
        SPGAUX=>GWFSPGDAT(IGRID)%SPGAUX
        SEEP=>GWFSPGDAT(IGRID)%SEEP
        QSPG=>GWFSPGDAT(IGRID)%QSPG
C
      RETURN
      END
      SUBROUTINE SGWF2SPG7PSV(IGRID)
C  Save SPG data for a grid.
      USE GWFSPGMODULE
C
        GWFSPGDAT(IGRID)%NSEEP=>NSEEP
        GWFSPGDAT(IGRID)%MXSPG=>MXSPG
        GWFSPGDAT(IGRID)%NSPGVL=>NSPGVL
        GWFSPGDAT(IGRID)%ISPGCB=>ISPGCB
        GWFSPGDAT(IGRID)%IPRSPG=>IPRSPG
        GWFSPGDAT(IGRID)%NPSPG=>NPSPG
        GWFSPGDAT(IGRID)%ISPGPB=>ISPGPB
        GWFSPGDAT(IGRID)%NNPSPG=>NNPSPG
        GWFSPGDAT(IGRID)%SPGAUX=>SPGAUX
        GWFSPGDAT(IGRID)%SEEP=>SEEP
        GWFSPGDAT(IGRID)%QSPG=>QSPG
C
      RETURN
      END

ext: f

max_stars_repo_name: etiennebresciani/modflow-spg

max_issues_repo_name: etiennebresciani/modflow-spg

max_forks_repo_name: etiennebresciani/modflow-spg


text:
Jake Lippmans Page... BORANG more coming...

20050630 00:41:08 nbsp hi Users/JakeLippman

20050630 10:48:13 nbsp hi. Users/GeorgeLewis

20071229 14:11:17 nbsp is this the martinez jake lippman? Users/gager

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      MODULE MADLIB	!Messing with COMMON is less convenient.
       INTEGER MSG,KBD,INF		!I/O unit numbers.
       DATA MSG,KBD,INF/6,5,10/		!Output, input, some disc file.
       INTEGER LSTASH,NSTASH,MSTASH	!Prepare a common text stash.
       PARAMETER (LSTASH = 246810, MSTASH = 6666)	!LSTASH characters for MSTASH texts.
       CHARACTER*(LSTASH) STASH		!The pool.
       INTEGER ISTASH(MSTASH + 1)	!Fingers start positions, and thus end positions by extension.
       DATA NSTASH,ISTASH(1)/0,1/	!Empty pool: no entries, first available character is at 1.
       INTEGER MANYLINES,MANYTESTS	!I also want some lists of texts.
       PARAMETER (MANYLINES = 1234)	!This is to hold the story.
       INTEGER NSTORY,STORY(MANYLINES)	!Fingering texts in the stash.
       PARAMETER (MANYTESTS = 1234)	!Likewise, some target/replacement texts.
       INTEGER NTESTS,TARGET(MANYTESTS),REPLACEMENT(MANYTESTS)	!Thus.
       DATA NSTORY,NTESTS/0,0/		!No story lines, and no tests.
       INTEGER STACKLIMIT		!A recursion limit.
       PARAMETER (STACKLIMIT = 28)	!This should suffice.

       CONTAINS
        SUBROUTINE CROAK(GASP)	!A dying remark.
         CHARACTER*(*) GASP	!The last words.
         WRITE (MSG,*) "Oh dear."	!Shock.
         WRITE (MSG,*) GASP		!Aargh!
         STOP "How sad."		!Farewell, cruel world.
        END SUBROUTINE CROAK	!Farewell...

        SUBROUTINE SHOWSTASH(BLAH,I)	!One might be wondering.
         CHARACTER*(*) BLAH		!An annotation.
         INTEGER I			!The desired stashed text.
          WRITE (MSG,1) BLAH,I,STASH(ISTASH(I):ISTASH(I + 1) - 1)	!Whee!
    1     FORMAT (A,': Text(',I0,')="',A,'"')	!Hopefully, helpful.
        END SUBROUTINE SHOWSTASH	!Ah, debugging.

        INTEGER FUNCTION EATTEXT(IN)	!Add a text to STASH and finger it.
Co-opts the as-yet unused space in STASH as its scratchpad.
         INTEGER IN	!Input from this I/O unit number.
         INTEGER I,N,L	!Fingers.
          I = ISTASH(NSTASH + 1)!First available position in STASH.
          N = LSTASH - I + 1	!Number of characters yet unused.
          IF (N.LT.666) CALL CROAK("Insufficient STASH space remains!")
          READ (IN,1,END = 66) L,STASH(I:I + MIN(L,N) - 1)	!Calculated during the read.
    1     FORMAT (Q,A)		!Obviously, Q = character count incoming, A = accept all of them.
          L = I + MIN(L,N) - 1	!The last character read.
   10     IF (L.LT.I) GO TO 66	!A blank line! Deemed end-of-file.
          IF (ICHAR(STASH(L:L)).LE.ICHAR(" ")) THEN	!A trailing space?
            L = L - 1		!Yes. Pull back.
            GO TO 10		!And try again.
          END IF		!So much for trailing spaces.
          IF (NSTASH.GE.MSTASH) CALL CROAK("Too many texts!")
          NSTASH = NSTASH + 1	!Admit another text.
          ISTASH(NSTASH + 1) = L + 1	!The start point of the following text.
          EATTEXT = NSTASH	!STASH(ISTASH(n):ISTASH(n + 1) - 1) has text n.
         RETURN			!All well.
   66     EATTEXT = 0		!Sez: "No text".
        END FUNCTION EATTEXT	!Rather odd side effects.

        INTEGER FUNCTION ADDSTASH(TEXT)	!Appends an arbitrary text to the pool of stashed texts.
         CHARACTER*(*) TEXT	!The stuff.
         INTEGER I		!A finger.
          IF (NSTASH.GE.MSTASH) CALL CROAK("The text pool is crowded!")	!Alas.
          I = ISTASH(NSTASH + 1)	!First unused character.
          IF (I + LEN(TEXT).GT.LSTASH) CALL CROAK("Overtexted!")	!Alack.
          STASH(I:I + LEN(TEXT) - 1) = TEXT	!Place.
          NSTASH = NSTASH + 1			!Count in another entry.
          ISTASH(NSTASH + 1) = I + LEN(TEXT)	!The new "first available" position.
          ADDSTASH = NSTASH	!Pass a finger back to the caller.
        END FUNCTION ADDSTASH	!Just an integer.

        INTEGER FUNCTION ANOTHER(TEXT)	!Possibly add TEXT to the table of target texts.
Collects TARGET REPLACEMENT pairs (increasing NTESTS) as directed by INSPECT.
         CHARACTER*(*) TEXT	!The text of the target.
         INTEGER I,IT		!Steppers.
          ANOTHER = 0		!Possibly, the text is already in the table.
          DO I = 1,NTESTS	!So, step through the known target texts.
            IT = TARGET(I)		!Finger a target text.
            IF (TEXT.EQ.STASH(ISTASH(IT):ISTASH(IT + 1) - 1)) RETURN	!Already have this one.
          END DO		!Otherwise, try the next.
          IF (NTESTS.GE.MANYTESTS) CALL CROAK("Too many tests!")	!Oh dear.
          NTESTS = NTESTS + 1		!Count in another.
          TARGET(NTESTS) = ADDSTASH(TEXT)!Stash its text and get a finger to it.
          ANOTHER = NTESTS		!My caller will want to know which test.
          WRITE (MSG,1) TEXT		!Now request the replacement text.
    1     FORMAT ("Enter your text for ",A,": ",$)	!Obviously, the $ indicates "no new line".
          REPLACEMENT(NTESTS) = EATTEXT(KBD)	!Zero for "no text".
        END FUNCTION ANOTHER	!Produces entries for TARGET and REPLACEMENT.

        SUBROUTINE INSPECT(X)	!Examine text number X for the special <...> sequence.
Calls for inspection of REPLACEMENT texts as well, should ANOTHER report a new entry.
         INTEGER X	!Fingers the text in STASH via ISTASH(X).
         INTEGER MARK	!Recalls where the < was found.
         INTEGER IT,NEW	!Fingers to entries in STASH.
         INTEGER I	!A stepper.
         INTEGER SP,STACK(STACKLIMIT)	!Prepare for some recursion.
          SP = 1		!Start with the starter.
          STACK(1) = X		!Stack up.
          DO WHILE(SP.GT.0)	!While texts are yet uninspected,
            IT = STACK(SP)		!Finger one.
            SP = SP - 1			!Working down the stack.
            MARK = 0			!Uninitialised variables are bad.
            DO I = ISTASH(IT),ISTASH(IT + 1) - 1!Step through the stashed text.
              IF (STASH(I:I).EQ."<") THEN	!Is it the starter?
                MARK = I			!Yes. Remember where it is.
              ELSE IF (STASH(I:I).EQ.">") THEN	!The ender?
                IF (MARK.LE.0) CALL CROAK("A > with no preceeding <!")	!Bah.
                NEW = ANOTHER(STASH(MARK:I))	!Consider the spanned text.
                IF (NEW.GT.0) THEN		!If that became a new table entry,
                  IF (SP.GE.STACKLIMIT) CALL CROAK("Stack overflow!")	!Its replacement is to be inspected.
                  SP = SP + 1			!But I'm still busy with the current text.
                  STACK(SP) = REPLACEMENT(NEW)	!So, stack it for later.
                END IF			!So much for that <...> apparition.
                MARK = 0		!Be ready to check afresh for the next.
              END IF		!So much for that character.
            END DO		!On to the next.
          END DO	!So much for that stacked entry.
        END SUBROUTINE INSPECT	!WRITESTORY will rescan the story lines.

        SUBROUTINE READSTORY(IN)!Read and stash the lines.
         INTEGER IN		!Input from here.
         INTEGER LINE		!A finger to the story line.
   10    LINE = EATTEXT(IN)	!So, grab a line.
         IF (LINE.GT.0) THEN	!A live line?
           NSTORY = NSTORY + 1	!Yes.Count it in.
           STORY(NSTORY) = LINE	!Save it in the story list.
           CALL INSPECT(LINE)	!Look for trouble as well.
           GO TO 10		!And go for the next line.
         END IF			!Oh for while (Line:=EatText(in)) > 0 do SaveAndInspect(Line);
        END SUBROUTINE READSTORY!Simple enough, anyway.

        SUBROUTINE WRITESTORY(WIDTH)	!Applying the replacements, with replacement replacement too.
Co-opts the as-yet unused space in STASH as its output scratchpad.
Can't rely on changing the index and bounds of a DO-loop on the fly.
         INTEGER WIDTH
         INTEGER LINE,IT,I,J	!Steppers.
         INTEGER L,L0,N		!Fingers.
         INTEGER TAIL,MARK,LAST	!Scan choppers.
         INTEGER SP,STACKI(STACKLIMIT),STACKL(STACKLIMIT)	!Ah, recursion.
          L0 = ISTASH(NSTASH + 1)	!The first available place in the stash.
          L = L0 - 1			!Syncopation for my output finger.
       LL:DO LINE = 1,NSTORY		!Step through the lines of the story.
            SP = 0			!Start with the task in hand.
            IT = STORY(LINE)		!Finger the stashed line.
            LAST = ISTASH(IT + 1) - 1	!Find its last character in STASH.
            I = ISTASH(IT)		!Find its first character in STASH.
            TAIL = I - 1		!Syncopation. No text from this line yet.
            IF (STASH(I:I).LE." ") THEN	!The line starts with a space?
              CALL BURP			!Yes. Flush, so as to start a new paragraph.
            ELSE IF (LINE.GT.1) THEN	!Otherwise, the line is a continuation.
              L = L + 1			!So, squeeze in a space as a separator.
              STASH(L:L) = " "		!Since its text follows on.
            END IF			!Now for the content of the line.
  666    II:DO WHILE(I.LE.LAST)		!Step along its text.
              IF (STASH(I:I).EQ."<") THEN	!Trouble starter?
                MARK = I			!Yes. Remember where.
              ELSE IF (STASH(I:I).EQ.">") THEN	!The corresponding ender?
                CALL APPEND(TAIL + 1,MARK - 1)	!Waiting text up to the mark.
             JJ:DO J = 1,NTESTS		!Step through the target texts.
                  IT = TARGET(J)		!Finger one.
                  IF (STASH(ISTASH(IT):ISTASH(IT + 1) - 1)	!Its stashed text.
     1            .EQ.STASH(MARK:I)) THEN		!Matches the suspect text?
                    IT = REPLACEMENT(J)		!Yes! Finger the replacement text.
                    IF (IT.GT.0) THEN	!Null replacements can be ignored.
                      IF (SP.GE.STACKLIMIT) CALL CROAK("StackOverflow!")	!Always diff. messages.
                      SP = SP + 1		!Interrupt the current scan.
                      STACKI(SP) = I		!Remember where we're up to,
                      STACKL(SP) = LAST		!And the end of the text.
                      I = ISTASH(IT) - 1	!One will be added shortly, at JJ+1.
                      LAST = ISTASH(IT + 1) - 1	!Preempt the scan-in-progress.
                    END IF			!To work along the replacement text.
                    EXIT JJ		!Found the target, so the search is finished.
                  END IF		!Otherwise,
                END DO JJ		!Try the next target text.
                TAIL = I		!Normal text resumes at TAIL + 1.
              END IF			!Enough analysis of that character from the story line.
              I = I + 1			!The next to consider.
            END DO II		!Perhaps we've finished this text.
            IF (SP.GT.0) THEN	!Yes! But, were we interrupted in a previous scan?
              CALL APPEND(TAIL + 1,LAST)!Yes! Roll the tail of the just-finished scan.
              TAIL = STACKI(SP)		!The stacked value of I was fingering a >.
              LAST = STACKL(SP)		!And this was the end of the text.
              SP = SP - 1		!So we've recovered where the scan was.
              I = TAIL + 1		!And this is the next to look at.
              GO TO 666			!Proceed to do so.
            END IF		!But if all is unstacked,
            CALL APPEND(TAIL + 1,LAST)	!Don't forget the tail end.
          END DO LL			!On to the next story line.
          CALL BURP		!Any waiting text must be less than WIDTH.
         CONTAINS		!Some assistants, defined after usage...
          SUBROUTINE APPEND(IST,LST)	!Has access to L.
           INTEGER IST,LST		!To copy STASH(IST:LST) to the scratchpad.
           INTEGER N			!The number of characters to copy.
            N = LST - IST + 1		!So find out.
            IF (N.LE.0) RETURN		!Avoid relying on zero-length action.
            IF (L + N.GT.LSTASH) CALL CROAK("Out of stash!")	!Oh dear.
            STASH(L + 1:L + N) = STASH(IST:LST)	!There they go.
            L = L + N			!Advance my oputput finger.
            IF (L - L0 + 1.GE.WIDTH) CALL BURP	!Enough to be going on with?
          END SUBROUTINE APPEND		!Few invocations, if with tricky parameters.
          SUBROUTINE BURP		!Flushes forth up to WIDTH characters.
           INTEGER N,W,L1		!And slides any remnant back.
            N = L - L0 + 1		!So, how many characters are waiting?
            IF (N.LE.WIDTH) THEN	!Too many for one line of output?
              L1 = L			!Nope. Roll the lot.
             ELSE			!Otherwise, a partial flush.
              W = L0 + WIDTH - 1	!Last character that can be fitted into WIDTH.
              DO L1 = W,L0,-1		!Look for a good split.
                IF (STASH(L1:L1).LE." ") EXIT	!Like, at a space.
              END DO			!Keep winding back.
              IF (L1.LE.L0) L1 = W	!No pleasing split found. Just roll a full width.
            END IF			!Ready to roll.
            WRITE (MSG,"(A)") STASH(L0:L1)	!Thus!
            IF (N.LE.WIDTH) THEN	!If the whole text was written,
              L = L0 - 1		!Then there is no text in the scratchpad.
             ELSE			!If only L0:L1 were written of L0:L,
              W = L0 + L - L1 - 1	!How far will the remaining text extend?
              STASH(L0:W) = STASH(L1 + 1:L)	!Shift it.
              L = W			!Finger the last used character position.
            END IF			!One trim is enough, even if the scracchpad contains multiple widths' worth..
          END SUBROUTINE BURP		!Since I don't want to flush the lot.
        END SUBROUTINE WRITESTORY	!Just a sequence of lines.
      END MODULE MADLIB		!Enough of that.

      PROGRAM MADLIBBER	!See, for example, https://en.wikipedia.org/wiki/Mad_Libs
      USE MADLIB
      WRITE (MSG,1)	!It's polite to explain.
    10FORMAT ("Reads a story in template form, containing special ",
     1 "entries such as <dog's name> amongst the text.",/,
     2 "You will be invited to supply a replacement text for each "
     3 "such entry, as encountered,",/,
     4 "after which the story will be presented with your ",
     5 "substitutions made.",//,
     6 "Here goes... Reading file Madlib.txt",/)
      OPEN(INF,STATUS="OLD",ACTION="READ",FORM="FORMATTED",
     1 FILE = "Madlib.txt")
      CALL READSTORY(INF)
      CLOSE(INF)
      WRITE (MSG,*)
      WRITE (MSG,*) "  Righto!"
      WRITE (MSG,*)
      CALL WRITESTORY(66)
      END

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
      subroutine fcomp1(indic2,ibloc,indi,h,g,d,w,w1,n,nr,ncs,
     &dga,delta,prop,acc,scale)
c     Copyright INRIA
      implicit double precision (a-h,o-z)
      dimension ibloc(n),indi(n),h(*),g(n),d(n),
     &w(n),w1(n),scale(n)
c
      ncs=0
      if(nr.eq.n) return
      zm=0.d0
      if(indic2.eq.1) go to 900
      delta=0.d0
      nh=nr*(nr+1)/2
      nrr=n-nr
      call fmlag1(n,nr,h,d,w)
      do 500 i=1,n
      ibi=ibloc(i)
      if(ibi.eq.0) go to 500
      gi=g(i)
      inc=indi(i)
      inc1=inc-1
      inr=inc-nr
      winc=w(inc)
      dmu=winc+gi
      am=dmin1(dabs(gi),dabs(dmu))
      if(2.d0*dabs(winc).ge.am) go to 500
      if(ibi.eq.-1.and.dmu.ge.0.d0) go to 500
      if(ibi.eq.1.and.dmu.le.0.d0) go to 500
      dmu=dabs(dmu)
      if(dmu*scale(i).le.acc) go to 500
      dmu1=dmu*dmu
      k=inr
      nh1=(inc1)*(n+1)-(inc1)*inc/2+1
      z=h(nh1)
      if(nr.eq.0) go to 350
      do 200 j=1,nr
      w1(j)=h(nh+k)
  200 k=k+nrr
      call fmc11e(h,nr,w1,w1,nr)
      k=inr
      do 250 j=1,nr
      z=z-w1(j)*h(nh+k)
  250 k=k+nrr
  350 dmu1=dmu1/z
      if(dmu1.le.delta) go to 500
      delta=dmu1
      ncs=i
      zm=dmu
  500 continue
      if(ncs.eq.0) return
      if(delta.le.-prop*dga)ncs=0
      return
  900 do 910 i=1,n
      ibi=ibloc(i)
      if(ibi.eq.0) go to 910
      dmu=g(i)
      if(ibi.eq.-1.and.dmu.ge.0.d0) go to 910
      if(ibi.eq.1.and.dmu.le.0.d0) go to 910
      dmu=dabs(dmu)*scale(i)
      if(dmu.le.zm) go to 910
      zm=dmu
      ncs=i
  910 continue
      if(zm.le.acc) ncs=0
      return
      end

max_stars_repo_name: alpgodev/numx

max_issues_repo_name: alpgodev/numx

max_forks_repo_name: alpgodev/numx


text:
use LinearSolverClass
use TimeClass
implicit none

type(LinearSolver_) :: solver
type(Time_) :: time
real(real64) :: Amat(5,5),mm(5)
integer(int32) :: i, j

Amat(1,:) = [ 1.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0 ]
Amat(2,:) = [ 2.0d0, 2.0d0, 0.0d0, 0.0d0, 0.0d0 ]
Amat(3,:) = [ 0.0d0, 3.0d0, 3.0d0, 0.0d0, 0.0d0 ]
Amat(4,:) = [ 2.0d0, 2.0d0, 2.0d0, 1.0d0, 0.0d0 ]
Amat(5,:) = [ 0.0d0, 0.0d0, 0.0d0, 0.0d0, 1.0d0 ]

solver%val = [1.0d0,2.0d0,2.0d0,3.0d0,3.0d0,2.0d0,2.0d0,2.0d0,1.0d0,1.0d0]
solver%Index_I = [1, 2, 2, 3, 3, 4,4,4,4,5]
solver%Index_J = [1, 1, 2, 2, 3, 1,2,3,4,5]

solver%b = [1.0d0, 2.0d0, 3.0d0,4.0d0,5.0d0]
solver%b_Index_J = [1, 2, 3, 4, 5]
call solver%convertCOOtoCRS()

print *, solver%CRS_val
print *, solver%CRS_Index_I
print *, solver%CRS_Index_J

call time%start()
mm = matmul(Amat,solver%b)
call time%show()
print *, mm

mm(:) = 0.0d0
call time%start()
mm = solver%matmulCRS()
call time%show()
print *, mm

end

max_stars_repo_name: kazulagi/plantFEM


text:
less_toxic(ll1,aa1).
less_toxic(bb1,t1).
less_toxic(x1,t1).
less_toxic(q1,p1).
less_toxic(s1,v1).
less_toxic(kk1,d1).
less_toxic(g1,l1).
less_toxic(g1,u1).
less_toxic(n1,l1).
less_toxic(ii1,d1).
less_toxic(j1,i1).
less_toxic(bb1,dd1).
less_toxic(ff1,w1).
less_toxic(b1,i1).
less_toxic(kk1,h1).
less_toxic(l1,p1).
less_toxic(u1,p1).
less_toxic(n1,v1).
less_toxic(r1,i1).
less_toxic(x1,aa1).
less_toxic(n1,e1).
less_toxic(u1,f1).
less_toxic(ii1,v1).
less_toxic(kk1,f1).
less_toxic(j1,e1).
less_toxic(a1,d1).
less_toxic(ee1,w1).
less_toxic(x1,w1).
less_toxic(y1,a1).
less_toxic(kk1,a1).
less_toxic(h1,p1).
less_toxic(k1,a1).
less_toxic(bb1,i1).
less_toxic(n1,cc1).
less_toxic(s1,f1).
less_toxic(x1,dd1).
less_toxic(n1,f1).
less_toxic(q1,w1).
less_toxic(cc1,d1).
less_toxic(m1,e1).
less_toxic(o1,l1).
less_toxic(j1,f1).
less_toxic(ll1,u1).
less_toxic(k1,z1).
less_toxic(q1,aa1).
less_toxic(hh1,aa1).
less_toxic(g1,f1).
less_toxic(f1,c1).
less_toxic(r1,p1).
less_toxic(kk1,cc1).
less_toxic(k1,b1).
less_toxic(j1,l1).
less_toxic(k1,aa1).
less_toxic(x1,f1).
less_toxic(g1,v1).
less_toxic(hh1,z1).
less_toxic(ff1,h1).
less_toxic(w1,e1).
less_toxic(r1,f1).
less_toxic(w1,d1).
less_toxic(v1,d1).
less_toxic(ee1,u1).
less_toxic(j1,d1).
less_toxic(s1,cc1).
less_toxic(hh1,l1).
less_toxic(k1,w1).
less_toxic(w1,t1).
less_toxic(ll1,w1).
less_toxic(ii1,l1).
less_toxic(ee1,z1).
less_toxic(y1,b1).
less_toxic(r1,e1).
less_toxic(n1,aa1).
less_toxic(q1,d1).
less_toxic(j1,v1).
less_toxic(b1,e1).
less_toxic(y1,c1).
less_toxic(ll1,dd1).
less_toxic(w1,i1).
less_toxic(hh1,e1).
less_toxic(jj1,p1).
less_toxic(n1,i1).
less_toxic(cc1,a1).
less_toxic(r1,l1).
less_toxic(ee1,d1).
less_toxic(b1,a1).
less_toxic(q1,jj1).
less_toxic(ee1,aa1).
less_toxic(ff1,p1).
less_toxic(x1,jj1).
less_toxic(w1,aa1).
less_toxic(m1,d1).
less_toxic(dd1,c1).
less_toxic(s1,p1).
less_toxic(y1,i1).
less_toxic(g1,c1).
less_toxic(f1,d1).
less_toxic(ii1,p1).
less_toxic(o1,c1).
less_toxic(j1,cc1).
less_toxic(q1,e1).
less_toxic(q1,a1).
less_toxic(ll1,jj1).
less_toxic(b1,u1).
less_toxic(o1,b1).
less_toxic(y1,f1).
less_toxic(y1,w1).
less_toxic(o1,cc1).
less_toxic(b1,l1).
less_toxic(kk1,aa1).
less_toxic(ll1,h1).
less_toxic(s1,i1).
less_toxic(m1,cc1).
less_toxic(x1,u1).
less_toxic(q1,h1).
less_toxic(ff1,i1).
less_toxic(j1,h1).
less_toxic(o1,t1).
less_toxic(cc1,u1).
less_toxic(o1,u1).
less_toxic(hh1,f1).
less_toxic(m1,i1).
less_toxic(g1,aa1).
less_toxic(y1,aa1).
less_toxic(bb1,b1).
less_toxic(q1,f1).
less_toxic(ii1,u1).
less_toxic(i1,f1).
less_toxic(r1,h1).
less_toxic(j1,w1).
less_toxic(k1,t1).
less_toxic(r1,cc1).
less_toxic(m1,b1).
less_toxic(ii1,t1).
less_toxic(ee1,v1).
less_toxic(ff1,v1).
less_toxic(o1,aa1).
less_toxic(w1,v1).
less_toxic(g1,jj1).
less_toxic(b1,h1).
less_toxic(a1,f1).
less_toxic(ll1,i1).
less_toxic(q1,dd1).
less_toxic(y1,l1).
less_toxic(cc1,i1).
less_toxic(x1,l1).
less_toxic(x1,b1).
less_toxic(w1,a1).
less_toxic(kk1,l1).
less_toxic(p1,c1).
less_toxic(n1,h1).
less_toxic(kk1,e1).
less_toxic(s1,w1).
less_toxic(i1,d1).
less_toxic(bb1,e1).
less_toxic(a1,p1).
less_toxic(e1,p1).
less_toxic(ii1,a1).
less_toxic(ii1,dd1).
less_toxic(n1,jj1).
less_toxic(aa1,d1).
less_toxic(q1,z1).
less_toxic(g1,a1).
less_toxic(v1,f1).
less_toxic(r1,a1).
less_toxic(w1,f1).
less_toxic(ff1,cc1).
less_toxic(j1,u1).
less_toxic(q1,b1).
less_toxic(s1,l1).
less_toxic(m1,w1).
less_toxic(g1,z1).
less_toxic(p1,f1).
less_toxic(m1,v1).
less_toxic(r1,u1).
less_toxic(n1,d1).
less_toxic(e1,d1).
less_toxic(z1,d1).
less_toxic(ll1,l1).
less_toxic(bb1,v1).
less_toxic(ll1,b1).
less_toxic(y1,jj1).
less_toxic(b1,v1).
less_toxic(ll1,z1).
less_toxic(ff1,aa1).
less_toxic(m1,l1).
less_toxic(cc1,v1).
less_toxic(dd1,p1).
less_toxic(j1,dd1).
less_toxic(ii1,z1).
less_toxic(k1,v1).
less_toxic(o1,e1).
less_toxic(cc1,aa1).
less_toxic(q1,t1).
less_toxic(y1,dd1).
less_toxic(ii1,jj1).
less_toxic(hh1,b1).
less_toxic(k1,f1).
less_toxic(b1,aa1).
less_toxic(m1,a1).
less_toxic(hh1,cc1).
less_toxic(b1,d1).
less_toxic(j1,p1).
less_toxic(z1,p1).
less_toxic(q1,v1).
less_toxic(m1,f1).
less_toxic(m1,t1).
less_toxic(r1,t1).
less_toxic(kk1,u1).
less_toxic(g1,w1).
less_toxic(k1,h1).
less_toxic(x1,c1).
less_toxic(g1,b1).
less_toxic(k1,cc1).
less_toxic(x1,d1).
less_toxic(m1,u1).
less_toxic(t1,f1).
less_toxic(j1,jj1).
less_toxic(ff1,dd1).
less_toxic(t1,p1).
less_toxic(bb1,aa1).
less_toxic(kk1,dd1).
less_toxic(q1,c1).
less_toxic(k1,dd1).
less_toxic(cc1,h1).
less_toxic(x1,cc1).
less_toxic(q1,cc1).
less_toxic(x1,a1).
less_toxic(t1,c1).
less_toxic(g1,t1).
less_toxic(u1,d1).
less_toxic(kk1,t1).
less_toxic(cc1,e1).
less_toxic(ii1,e1).
less_toxic(ii1,aa1).
less_toxic(g1,dd1).
less_toxic(b1,jj1).
less_toxic(j1,b1).
less_toxic(hh1,v1).
less_toxic(ee1,f1).
less_toxic(s1,z1).
less_toxic(m1,z1).
less_toxic(bb1,u1).
less_toxic(m1,h1).
less_toxic(ii1,i1).
less_toxic(ii1,c1).
less_toxic(q1,i1).
less_toxic(s1,u1).
less_toxic(hh1,a1).
less_toxic(bb1,w1).
less_toxic(g1,h1).
less_toxic(hh1,d1).
less_toxic(y1,v1).
less_toxic(r1,jj1).
less_toxic(hh1,h1).
less_toxic(bb1,l1).
less_toxic(y1,z1).
less_toxic(ee1,p1).
less_toxic(o1,p1).
less_toxic(x1,p1).
less_toxic(bb1,c1).
less_toxic(ii1,w1).
less_toxic(kk1,b1).
less_toxic(v1,p1).
less_toxic(ll1,a1).
less_toxic(r1,dd1).
less_toxic(k1,d1).
less_toxic(n1,t1).
less_toxic(s1,d1).
less_toxic(jj1,d1).
less_toxic(ff1,l1).
less_toxic(ee1,a1).
less_toxic(n1,w1).
less_toxic(x1,z1).
less_toxic(j1,t1).
less_toxic(y1,cc1).
less_toxic(ff1,d1).
less_toxic(kk1,jj1).
less_toxic(b1,dd1).
less_toxic(x1,e1).
less_toxic(x1,h1).
less_toxic(h1,d1).
less_toxic(v1,c1).
less_toxic(k1,l1).
less_toxic(r1,v1).
less_toxic(ff1,t1).
less_toxic(g1,i1).
less_toxic(s1,a1).
less_toxic(o1,w1).
less_toxic(z1,c1).
less_toxic(h1,f1).
less_toxic(n1,c1).
less_toxic(y1,u1).
less_toxic(y1,e1).
less_toxic(hh1,t1).
less_toxic(x1,v1).
less_toxic(h1,c1).
less_toxic(ll1,v1).
less_toxic(ff1,c1).
less_toxic(ll1,cc1).
less_toxic(u1,c1).
less_toxic(kk1,i1).
less_toxic(ee1,c1).
less_toxic(g1,e1).
less_toxic(bb1,f1).
less_toxic(kk1,c1).
less_toxic(ee1,t1).
less_toxic(r1,d1).
less_toxic(ii1,b1).
less_toxic(k1,u1).
less_toxic(o1,jj1).
less_toxic(ee1,cc1).
less_toxic(q1,l1).
less_toxic(g1,p1).
less_toxic(s1,dd1).
less_toxic(l1,c1).
less_toxic(s1,c1).
less_toxic(y1,p1).
less_toxic(n1,a1).
less_toxic(l1,f1).
less_toxic(ll1,f1).
less_toxic(dd1,f1).
less_toxic(i1,c1).
less_toxic(j1,aa1).
less_toxic(b1,f1).
less_toxic(i1,y1).
less_toxic(jj1,r1).
less_toxic(dd1,s1).
less_toxic(b1,x1).
less_toxic(t1,y1).
less_toxic(i1,ll1).
less_toxic(i1,n1).
less_toxic(i1,o1).
less_toxic(i1,m1).
less_toxic(t1,ll1).
less_toxic(u1,x1).
less_toxic(p1,w1).
less_toxic(dd1,cc1).
less_toxic(a1,ii1).
less_toxic(c1,hh1).
less_toxic(p1,cc1).
less_toxic(i1,g1).
less_toxic(p1,ll1).
less_toxic(i1,kk1).
less_toxic(w1,bb1).
less_toxic(f1,ii1).
less_toxic(c1,l1).
less_toxic(f1,p1).
less_toxic(p1,k1).
less_toxic(aa1,bb1).
less_toxic(v1,r1).
less_toxic(z1,ii1).
less_toxic(f1,u1).
less_toxic(d1,f1).
less_toxic(w1,hh1).
less_toxic(aa1,n1).
less_toxic(b1,y1).
less_toxic(z1,s1).
less_toxic(cc1,ll1).
less_toxic(h1,k1).
less_toxic(a1,r1).
less_toxic(v1,k1).
less_toxic(p1,g1).
less_toxic(c1,x1).
less_toxic(f1,b1).
less_toxic(i1,s1).
less_toxic(f1,n1).
less_toxic(z1,g1).
less_toxic(cc1,y1).
less_toxic(w1,ll1).
less_toxic(f1,hh1).
less_toxic(dd1,ee1).
less_toxic(l1,q1).
less_toxic(l1,cc1).
less_toxic(dd1,j1).
less_toxic(w1,s1).
less_toxic(d1,ii1).
less_toxic(h1,kk1).
less_toxic(dd1,hh1).
less_toxic(f1,y1).
less_toxic(i1,ff1).
less_toxic(v1,n1).
less_toxic(v1,bb1).
less_toxic(f1,i1).
less_toxic(aa1,r1).
less_toxic(i1,x1).
less_toxic(f1,m1).
less_toxic(aa1,cc1).
less_toxic(dd1,m1).
less_toxic(p1,u1).
less_toxic(f1,h1).
less_toxic(e1,s1).
less_toxic(jj1,m1).
less_toxic(f1,o1).
less_toxic(c1,i1).
less_toxic(v1,x1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


text:
!-*- mode: compilation; default-directory: "/tmp/" -*-
!Compilation started at Thu Jun  5 01:52:03
!
!make f && for a in '' a bark book treat common squad confuse ; do echo $a | ./f ; done
!gfortran -std=f2008 -Wall -fopenmp -ffree-form -fall-intrinsics -fimplicit-none -g f.f08 -o f
! T
! T  A                    NA
! T  BARK                 BO NA RE XK
! F  BOOK                 OB BO -- --
! T  TREAT                GT RE ER NA TG
! F  COMMON               PC OB ZM -- -- --
! T  SQUAD                FS DQ HU NA QD
! T  CONFUSE              CP BO NA FS HU FS RE
!
!Compilation finished at Thu Jun  5 01:52:03

program abc
  implicit none
  integer, parameter :: nblocks = 20
  character(len=nblocks) :: goal
  integer, dimension(nblocks) :: solution
  character(len=2), dimension(0:nblocks) :: blocks_copy, blocks = &
       &(/'--','BO','XK','DQ','CP','NA','GT','RE','TG','QD','FS','JW','HU','VI','AN','OB','ER','FS','LY','PC','ZM'/)
  logical :: valid
  integer :: i, iostat
  read(5,*,iostat=iostat) goal
  if (iostat .ne. 0) goal = ''
  call ucase(goal)
  solution = 0
  blocks_copy = blocks
  valid = assign_block(goal(1:len_trim(goal)), blocks, solution, 1)
  write(6,*) valid, ' '//goal, (' '//blocks_copy(solution(i)), i=1,len_trim(goal))

contains

  recursive function assign_block(goal, blocks, solution, n) result(valid)
    implicit none
    logical :: valid
    character(len=*), intent(in) :: goal
    character(len=2), dimension(0:), intent(inout) :: blocks
    integer, dimension(:), intent(out) :: solution
    integer, intent(in) :: n
    integer :: i
    character(len=2) :: backing_store
    valid = .true.
    if (len(goal)+1 .eq. n) return
    do i=1, size(blocks)
       if (index(blocks(i),goal(n:n)) .ne. 0) then
          backing_store = blocks(i)
          blocks(i) = ''
          solution(n) = i
          if (assign_block(goal, blocks, solution, n+1)) return
          blocks(i) = backing_store
       end if
    end do
    valid = .false.
    return
  end function assign_block

  subroutine ucase(a)
    implicit none
    character(len=*), intent(inout) :: a
    integer :: i, j
    do i = 1, len_trim(a)
       j = index('abcdefghijklmnopqrstuvwxyz',a(i:i))
       if (j .ne. 0) a(i:i) = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'(j:j)
    end do
  end subroutine ucase

end program abc