text:
! Time of File Save by ERB: 3/15/2004 11:23AM
C     Last change:  ERB   3 Sep 2002    2:20 pm
      SUBROUTINE OBS1RIV6AL(IURVOB,NQ,NQC,NQT,IOUT,NQRV,NQTRV,IOBSUM,
     &                     LCOBRIV,ITMXP,LCSSRV,ISUM,IOBS)
C     VERSION 20000125
C     ******************************************************************
C     ALLOCATE ARRAY STORAGE FOR FLOW OBSERVATIONS AT RIVER CELLS
C     ******************************************************************
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      INTEGER IOUT, IURVOB, NQ, NQC, NQT, NQRV, NQCRV, NQTRV
      CHARACTER*200 LINE
C     ------------------------------------------------------------------
C     IDENTIFY PROCESS
      WRITE(IOUT,490) IURVOB
  490 FORMAT(/,' OBS1RIV6 -- OBSERVATION PROCESS (RIVER FLOW ',
     &    'OBSERVATIONS)',/,' VERSION 1.0, 10/15/98',/,
     &    ' INPUT READ FROM UNIT ',I4)
C
C  Turn off observation package if OBS is not active
      IF(IOBS.LE.0) THEN
        WRITE(IOUT,610)
610     FORMAT(/,1X,'WARNING: OBSERVATION (OBS) FILE IS NOT LISTED BUT',
     &      ' THE RIV OBSERVATION',/,' FILE (RVOB) IS',
     &     ' LISTED -- TURNING OFF RIV OBSERVATIONS (OBS1RIV6AL)')
        IURVOB = 0
        RETURN
      ENDIF
C
C  Read RVOB data
      CALL URDCOM(IURVOB,IOUT,LINE)
      LLOC = 1
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,NQRV,DUM,IOUT,IURVOB)
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,NQCRV,DUM,IOUT,IURVOB)
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,NQTRV,DUM,IOUT,IURVOB)
      WRITE (IOUT,500) NQRV, NQCRV, NQTRV
  500 FORMAT (/,
     &     ' NUMBER OF FLOW-OBSERVATION RIVER-CELL GROUPS.....: ',I6,/,
     &     '   NUMBER OF CELLS IN RIVER-CELL GROUPS...........: ',I6,/,
     &     '   NUMBER OF RIVER-CELL FLOWS.....................: ',I6)
C
      NQ = NQ + NQRV
      NQC = NQC + NQCRV
      NQT = NQT + NQTRV
      LCSSRV = ISUM
      ISUM = ISUM + ITMXP + 1
C     POINTER TO OBSERVATION ARRAYS
      LCOBRIV = IOBSUM
      IOBSUM = IOBSUM + NQTRV
C
      RETURN
      END
C=======================================================================
      SUBROUTINE OBS1RIV6RP(NCOL,NROW,NPER,IURVOB,IOUT,OBSNAM,NHT,JT,
     &                      IBT,NQOB,NQCL,IQOB,QCLS,IERR,HOBS,TOFF,
     &                      WTQ,IOWTQ,IPRN,NDMH,NSTP,PERLEN,TSMULT,
     &                      ISSA,ITRSS,NQAR,NQCAR,NQTAR,IQ1,NQT1,NDD,
     &                      IURV,NQRV,NQTRV,NT,NC,IPLOT,NAMES,ND,IPR,
     &                      MPR,IOWTQRV,OTIME)
C     VERSION 20010921 ERB
C     ******************************************************************
C     READ, CHECK AND STORE FLOW-OBSERVATION DATA FOR RIVER BOUNDARIES.
C     ******************************************************************
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      REAL BLANK, EVFRV, HOBS, PERLEN, QCLS, TOFF, TOFFSET, TOMULTRV,
     &     TSMULT, WTQ
      INTEGER I, I4, IBT, IERR, IOUT, IOWTQ, IPRN, IQ, IQOB, IURVOB,
     &        IWT, J, JT, L, N, NC, NC1, NC2, NCOL, NDMH, NHT, NPER,
     &        NQCL, NQOB, NROW, NSTP, NT, NT1, NT2, IURV, ISSA
      INTEGER IPLOT(ND+IPR+MPR)
      CHARACTER*12 OBSNAM(NDD), NAMES(ND+IPR+MPR)
      CHARACTER*20 FMTIN*20, ANAME*43
      DIMENSION IBT(2,NQAR), NQOB(NQAR), NQCL(NQAR), IQOB(NQTAR),
     &          QCLS(5,NQCAR),HOBS(ND), TOFF(ND), NSTP(NPER),
     &          PERLEN(NPER), TSMULT(NPER), ISSA(NPER)
      DIMENSION WTQ(NDMH,NDMH)
      CHARACTER*10 STATYP(0:2)
      REAL OTIME(ND)
      DATA (STATYP(I),I=0,2)/'VARIANCE','STD. DEV.','COEF. VAR.'/
      DATA ANAME/'COVARIANCE OF RIVER-CELL FLOW OBSERVATIONS '/
C     ------------------------------------------------------------------
  500 FORMAT (15X,2F5.0,F10.0)
  505 FORMAT (8F10.0)
  517 FORMAT (/,' RIVER-CELL FLOW OBSERVATION VARIANCES ARE',
     &        ' MULTIPLIED BY: ',G15.4)
  520 FORMAT (/,' OBSERVED RIVER-CELL FLOW DATA',/,' -- TIME OFFSETS',
     &        ' ARE MULTIPLIED BY: ',G12.5)
  525 FORMAT (/,'   GROUP NUMBER: ',I6,'   BOUNDARY TYPE: ',A,
     &'   NUMBER OF CELLS IN GROUP: ',I6,/,
     &'   NUMBER OF FLOW OBSERVATIONS: ',I6,//,
     &40X,'OBSERVED',/,
     &20X,'REFER.',13X,'RIVER FLOW',/,
     &7X,'OBSERVATION',2X,'STRESS',4X,'TIME',5X,'GAIN (-) OR',14X,
     &'STATISTIC   PLOT',/,
     &2X,'OBS#    NAME',6X,'PERIOD   OFFSET',5X,'LOSS (+)',
     &4X,'STATISTIC     TYPE      SYM.')
  535 FORMAT (I6,1X,A12,2X,I4,2X,G11.4,1X,G11.4,1X,G11.4,2X,A10,
     &1X,I5)
  540 FORMAT (/,'       LAYER  ROW  COLUMN    FACTOR')
  550 FORMAT (4X,F8.0,F6.0,F7.0,F9.2)
  555 FORMAT (4X,F8.0,F6.0,F9.2)
  560 FORMAT (/,' FOR OBS ',I6,' STATISTIC RELATED TO WEIGHT < OR =0 --'
     &        ,' STOP EXECUTION (OBS1RIV6RP)',/)
  565 FORMAT (/,' RIVER PACKAGE',
     &        ' IS NOT OPEN -- STOP EXECUTION (OBS1RIV6RP)')
  590 FORMAT (/,' ROW OR COLUMN NUMBER INVALID',
     &        ' -- STOP EXECUTION (OBS1RIV6RP)',/)
  605 FORMAT (/,' OBSERVATION ',I6,' EQUALS ZERO, THE STATISTIC ',
     &        'CAN NOT BE A',/,' COEFFICIENT OF VARIATION (ISTAT=2)',
     &        ' -- STOP EXECUTION(OBS1RIV6RP)')
  615 FORMAT (//,1X,A,/,1X,42('-'))
  620 FORMAT (/,1X,'ERROR: SEARCH ABOVE FOR ERROR MESSAGE(S)',/,
     &' -- STOP EXECUTION (OBS1RIV6RP)')
C
C-----READ TIME-OFFSET MULTIPLIER FOR FLOW-OBSERVATION TIMES AND INPUT
C     ERROR VARIANCE FOR FLOW OBSERVATIONS
      READ(IURVOB,*) TOMULTRV, EVFRV, IOWTQRV
      IF (IOWTQRV.GT.0) IOWTQ = 1
C
C-------INITIALIZE VARIABLES
      IQ2 = IQ1 + NQRV - 1
      DO 20 IQ = IQ1,IQ2
        IBT(2,IQ) = 0
   20 CONTINUE
C-------WRITE TITLE AND LOOP THROUGH CELL GROUPS (READ ITEMS 3, 4, AND 5)
      WRITE (IOUT,517) EVFRV
      WRITE (IOUT,520) TOMULTRV
      DO 120 IQ = IQ1,IQ2
C       READ ITEM 3
        READ (IURVOB,*) NQOB(IQ), NQCL(IQ)
        IBT(1,IQ) = 1
        WRITE (IOUT,525) IQ, 'RIV', NQCL(IQ), NQOB(IQ)
C----------SET FLAG FOR SETTING ALL FACTORS TO 1
        I4 = 0
        IF (NQCL(IQ).LT.0) THEN
          I4 = 1
          NQCL(IQ) = -NQCL(IQ)
        ENDIF
C----------READ TIME STEPS, MEASURED FLOWS, AND WEIGHTS.
        NT1 = NT + 1
        NT2 = NT + NQOB(IQ)
        DO 30 J = NT1, NT2
          N = NHT + J
C---------READ ITEM 4
          IF (IOWTQRV.GT.0) THEN
            READ (IURVOB,*) OBSNAM(N), IREFSP, TOFFSET, HOBS(N), STAT,
     &                      ISP, IPLOT(N)
            NAMES(N) = OBSNAM(N)
            WRITE (IOUT,535) N,OBSNAM(N),IREFSP,TOFFSET,HOBS(N),0.0,
     &                       STATYP(ISP),IPLOT(N)
          ELSE
            READ (IURVOB,*) OBSNAM(N), IREFSP, TOFFSET, HOBS(N),
     &                      WTQ(J,J), IWT, IPLOT(N)
            NAMES(N) = OBSNAM(N)
            WRITE (IOUT,535) N,OBSNAM(N),IREFSP,TOFFSET,HOBS(N),
     &                       WTQ(J,J),STATYP(IWT),IPLOT(N)
            IF (HOBS(N).EQ.0 .AND. IWT.EQ.2) THEN
              WRITE (IOUT,605) N
              IERR = 1
            ENDIF
            IF (IWT.EQ.2) WTQ(J,J) = WTQ(J,J)*HOBS(N)
            IF (IWT.GT.0) WTQ(J,J) = WTQ(J,J)*WTQ(J,J)
            WTQ(J,J) = WTQ(J,J)*EVFRV
            IF (WTQ(J,J).LE.0.) THEN
              WRITE (IOUT,560) N
              IERR = 1
            ENDIF
          ENDIF
          CALL UOBSTI(OBSNAM(N),IOUT,ISSA,ITRSS,NPER,NSTP,IREFSP,
     &                IQOB(J),PERLEN,TOFF(N),TOFFSET,TOMULTRV,TSMULT,1,
     &                OTIME(N))
C----------ERROR CHECKING
          IF (IQOB(J).GE.JT) THEN
            JT = IQOB(J)
            IF (TOFF(N).GT.0.) JT = JT+1
          ENDIF
          IF (IURV.EQ.0) THEN
            WRITE (IOUT,565)
            IERR = 1
          ENDIF
   30   CONTINUE
C----------READ LAYER, ROW, COLUMN, AND FACTOR (ITEM 5)
        NC1 = NC + 1
        NC2 = NC + NQCL(IQ)
        WRITE (IOUT,540)
        DO 40 L = NC1, NC2
          READ (IURVOB,*) (QCLS(I,L),I=1,4)
          IF (I4.EQ.1) QCLS(4,L) = 1.
          WRITE (IOUT,550) (QCLS(I,L),I=1,4)
          I = QCLS(2,L)
          J = QCLS(3,L)
          IF (J.LE.0 .OR. J.GT.NCOL .OR. I.LE.0 .OR. I.GT.NROW) THEN
            WRITE (IOUT,590)
            IERR = 1
          ENDIF
   40   CONTINUE
C-------UPDATE COUNTERS
        NC = NC2
        NT = NT2
  120 CONTINUE
      IQ1 = IQ2 + 1
C
C-------READ FULL COVARIANCE MATRIX ON RIVER-CELL FLOW-OBSERVATION DATA
      IPRN = 0
      NQT2 = NQT1 + NQTRV - 1
      IF (IOWTQRV.GT.0 .AND. NQTRV.GT.0) THEN
C       READ ITEM 6
        READ (IURVOB,*) FMTIN, IPRN
C       READ ITEM 7
        DO 140 I = NQT1,NQT2
          READ (IURVOB,FMTIN) (BLANK,J=NQT1,I-1), (WTQ(I,J),J=I,NQT2)
          DO 130 J = I, NQT2
            WTQ(I,J) = WTQ(I,J)*EVFRV
            IF (I.EQ.J) THEN
              IF (WTQ(I,J).LT.0.) WTQ(I,J) = -WTQ(I,J)
            ELSE
              WTQ(J,I) = WTQ(I,J)
            ENDIF
  130     CONTINUE
  140   CONTINUE
        IF (IPRN.GE.0) THEN
          WRITE (IOUT,615) ANAME
          CALL UARRSUBPRW(WTQ,NDMH,NDMH,NQT1,NQT2,NQT1,NQT2,IPRN,IOUT,
     &                    OBSNAM(NHT+1),NDMH)
        ENDIF
      ENDIF
      NQT1 = NQT2 + 1
C
      IF (IERR.GT.0) THEN
        WRITE(IOUT,620)
        CALL USTOP(' ')
      ENDIF
C
      RETURN
      END
C=======================================================================
      SUBROUTINE OBS1RIV6FM(NQ,NQOB,NQCL,IQOB,QCLS,IBT,MXRIVR,NRIVER,
     &                      RIVR,HNEW,NCOL,NROW,NLAY,IOUT,IBOUND,
     &                      NHT,OBSNAM,H,TOFF,WTQ,NDMH,ITS,NQAR,NQCAR,
     &                      NQTAR,NRIVVL,ND)
C     VERSION 19981020 ERB
C     ******************************************************************
C     CALCULATE SIMULATED EQUIVALENTS TO OBSERVED FLOWS FOR THE RIVER
C     PACKAGE
C     ******************************************************************
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      REAL C, FACT, H, HB, QCLS, RBOT, RIVR, TOFF, WTQ, ZERO
      INTEGER I, IBOUND, IBT, IBT1, IFLAG, II, IOUT, IQ,
     &        IQOB, IRBOT, ITS, J, JJ, JRBOT, K, KK, KRBOT,
     &        MXRIVR, N, NB, NBN, NC, NC1, NC2, NCOL, NDMH, NHT,
     &        NLAY, NQ, NQCL, NQOB, NRIVER, NROW, NT,
     &        NT1, NT2
      CHARACTER*12 OBSNAM(ND)
      DOUBLE PRECISION HNEW(NCOL,NROW,NLAY)
      DOUBLE PRECISION HH, HHNEW
      DIMENSION RIVR(NRIVVL,MXRIVR), IBOUND(NCOL,NROW,NLAY),
     &          IBT(2,NQAR), NQOB(NQAR), NQCL(NQAR), IQOB(NQTAR),
     &          QCLS(5,NQCAR), H(ND), TOFF(ND), WTQ(NDMH,NDMH)
      INCLUDE 'param.inc'
C     ------------------------------------------------------------------
  500 FORMAT (/,
     &' HEADS AT RIVER CELLS ARE BELOW THE',
     &' BOTTOM OF THE RIVER BED AT THE CELLS LISTED',/,
     &' BELOW.  THESE CONDITIONS DIMINISH THE IMPACT',
     &' OF THE OBSERVATION ON ESTIMATES OF',/,
     &' ALL PARAMETERS EXCEPT THOSE THAT CONTROL THE HYDRAULIC',
     &' CONDUCTIVITY OF THE',/,
     &' RIVER BED.  (SEE TEXT FOR MORE INFORMATION).')
  505 FORMAT (/,' OBS# ',I6,', ID ',A,', TIME STEP ',I5)
  510 FORMAT ('    LAYER   ROW  COLUMN')
  520 FORMAT (3I7)
  530 FORMAT (I7,' OF THE',I7,' CELLS USED TO SIMULATE THE',
     &        ' GAIN OR LOSS ARE',/,22X,'AFFECTED.')
  535 FORMAT (' ALL CELLS INCLUDED IN THIS OBSERVATION ARE INACTIVE.  ',
     &        'THE OBSERVATION WILL',/
     &        ,' BE OMITTED FROM THIS PARAMETER-ESTIMATION ITERATION')
  540 FORMAT (' CELL OR REACH # ',I6,
     &        ' OF HEAD-DEP. BOUNDARY GAIN OR LOSS OBS#',I5,' ID=',A,/,
     &        ' NOT FOUND IN CELLS LISTED FOR RIVER PACKAGE',/,
     &        ' -- STOP EXECUTION (OBS1RIV6FM)')
C
C-------INITIALIZE VARIABLES
      ZERO = 0.0
      NC = 0
      NT1 = 1
      JRBOT = 0
C-------LOOP THROUGH BOUNDARY FLOWS
      DO 60 IQ = 1, NQ
        IBT1 = IBT(1,IQ)
        NT2 = NT1 + NQOB(IQ) - 1
        IF (IBT1.NE.1) GOTO 50
C----------WAS THERE A MEASUREMENT AT THIS BOUNDARY THIS TIME STEP?
        DO 40 NT = NT1, NT2
          IF (IQOB(NT).EQ.ITS .OR.
     &        (IQOB(NT).EQ.ITS-1.AND.TOFF(NHT+NT).GT.ZERO)) THEN
C----------ASSIGN VARIABLES ACCORDING TO BOUNDARY TYPE
            IRBOT = 0
            KRBOT = 0
            NBN = NRIVER
C----------LOOP THROUGH CELLS.
            NC1 = NC + 1
            NC2 = NC + NQCL(IQ)
            NB = 0
            DO 30 N = NC1, NC2
              K = QCLS(1,N)
              I = QCLS(2,N)
              J = QCLS(3,N)
C----------LOOP THROUGH DATA FILE TO FIND A MATCH.
              IFLAG = 0
              DO 10 MNB = 1, NBN
                NB = NB + 1
                IF (NB.GT.NBN) NB = 1
                KK = RIVR(1,NB)
                II = RIVR(2,NB)
                JJ = RIVR(3,NB)
C----------DO CALCULATIONS IF THIS IS A MATCH
                IF (I.EQ.II.AND.J.EQ.JJ.AND.K.EQ.KK) THEN
                  IFLAG = 1
                  IF (IBOUND(J,I,K).EQ.0) THEN
                    KRBOT = KRBOT + 1
                    GOTO 30
                  ENDIF
C-------------ASSIGN VARIABLE VALUES
                  HHNEW = HNEW(J,I,K)
                  HB = RIVR(4,NB)
                  C = RIVR(5,NB)
                  RBOT = RIVR(6,NB)
C-------------CALCULATE FLOWS
                  HH = C*(HB-HHNEW)
                  IF (HHNEW.LE.RBOT) THEN
                    HH = C*(HB-RBOT)
                    IF (JRBOT.EQ.0) WRITE (IOUT,500)
                    JRBOT = 1
                    IF (IRBOT.EQ.0) THEN
                      WRITE (IOUT,505) NHT + NT, OBSNAM(NHT+NT), ITS
                      WRITE (IOUT,510)
                    ENDIF
                    IRBOT = IRBOT + 1
                    WRITE (IOUT,520) K, I, J
                  ENDIF
                  GOTO 20
                ENDIF
   10         CONTINUE
              IF (IFLAG.EQ.0) THEN
                WRITE (IOUT,540) N, NHT + NT, OBSNAM(NHT+NT)
                CALL USTOP(' ')
              ENDIF
C-------------SUM VALUES FROM INDIVIDUAL CELLS.
C----------------CALCULATE FACTOR FOR TEMPORAL INTERPOLATION
   20         FACT = 1.0
              IF (TOFF(NHT+NT).GT.ZERO) THEN
                IF (IQOB(NT).EQ.ITS) FACT = 1. - TOFF(NHT+NT)
                IF (IQOB(NT).EQ.ITS-1) FACT = TOFF(NHT+NT)
              ENDIF
C---------------FLOWS
              H(NHT+NT) = H(NHT+NT) + HH*FACT*QCLS(4,N)
   30       CONTINUE
C-------PRINT NUMBER OF CELLS AT WHICH HEAD IS BELOW THE BOTTOM OF THE
C-------RIVER BED; CHECK FOR DISCONNECTED OBSERVATIONS.
            IF (IRBOT.GT.0) WRITE (IOUT,530) IRBOT, NQCL(IQ)
            IF (KRBOT.EQ.NQCL(IQ)) THEN
              WTQ(NT,NT) = -WTQ(NT,NT)
              WRITE (IOUT,535)
            ENDIF
          ENDIF
   40   CONTINUE
C-------UPDATE COUNTERS
   50   NC = NC + NQCL(IQ)
        NT1 = NT2 + 1
C
   60 CONTINUE
      RETURN
      END
C=======================================================================
      SUBROUTINE OBS1RIV6DR(NQ,NQOB,NQCL,IQOB,QCLS,IBT,MXRIVR,NRIVER,
     &                      RIVR,HNEW,IP,SNEW,NCOL,NROW,NLAY,IOUT,
     &                      IBOUND,NHT,X,OBSNAM,NPE,LN,TOFF,NPLIST,
     &                      ITS,NQAR,NQCAR,NQTAR,NRIVVL,IERR,IERRU,
     &                      ND)
C     VERSION 20000517 ERB
C     ******************************************************************
C     CALCULATE SENSITIVITIES FOR FLOW OBSERVATIONS FOR THE RIVER
C     PACKAGE
C     ******************************************************************
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      REAL C, FACT, HB, HHNEW, QCLS, RBOT, RIVR,
     &     TOFF, X, XX, ZERO
      INTEGER I, IBOUND, IBT, IBT1, IFLAG, II, IOUT, IP, IQ,
     &        IQOB, ITS, J, JJ, K, KK,
     &        LN, MXRIVR, N, NB, NBN, NC, NC1, NC2, NCOL, NHT,
     &        NLAY, NPE, NQ, NQCL, NQOB, NRIVER, NROW, NT,
     &        NT1, NT2
      CHARACTER*4 PIDTMP
      CHARACTER*12 OBSNAM(ND)
      DOUBLE PRECISION HNEW(NCOL,NROW,NLAY), SNEW(NCOL,NROW,NLAY)
      DIMENSION RIVR(NRIVVL,MXRIVR), IBOUND(NCOL,NROW,NLAY),
     &          X(NPE,ND), IBT(2,NQAR), LN(NPLIST), NQOB(NQAR),
     &          NQCL(NQAR), IQOB(NQTAR), QCLS(5,NQCAR), TOFF(ND)
      INCLUDE 'param.inc'
C     ------------------------------------------------------------------
  540 FORMAT (' CELL OR REACH # ',I6,
     &        ' OF HEAD-DEP. BOUNDARY GAIN OR LOSS OBS#',I5,' ID=',A,
     &        ' NOT FOUND IN CELLS LISTED FOR',/,' RIVER PACKAGE',
     &        ' -- STOP EXECUTION (OBS1RIV6DR)')
C
C-------INITIALIZE VARIABLES
      ZERO = 0.0
      NC = 0
      NT1 = 1
      IIPP = IPPTR(IP)
      PIDTMP = PARTYP(IIPP)
      IF (PIDTMP.EQ.'RIV ')
     &    CALL SOBS1RIV6QC(IBT,IIPP,NRIVVL,MXRIVR,NQ,NQCL,RIVR,QCLS,
     &                     NQAR,NQCAR)
C-------LOOP THROUGH BOUNDARY FLOWS
      DO 60 IQ = 1, NQ
        IBT1 = IBT(1,IQ)
        NT2 = NT1 + NQOB(IQ) - 1
        IF (IBT1.NE.1) GOTO 50
C----------WAS THERE A MEASUREMENT AT THIS BOUNDARY THIS TIME STEP?
        DO 40 NT = NT1, NT2
          IF (IQOB(NT).EQ.ITS .OR.
     &        (IQOB(NT).EQ.ITS-1.AND.TOFF(NHT+NT).GT.ZERO)) THEN
C----------ASSIGN VARIABLES ACCORDING TO BOUNDARY TYPE
            NBN = NRIVER
C----------LOOP THROUGH CELLS.
            NC1 = NC + 1
            NC2 = NC + NQCL(IQ)
            NB = 0
            DO 30 N = NC1, NC2
              K = QCLS(1,N)
              I = QCLS(2,N)
              J = QCLS(3,N)
C----------LOOP THROUGH DATA FILE TO FIND A MATCH.
              IFLAG = 0
              DO 10 MNB = 1, NBN
                NB = NB + 1
                IF (NB.GT.NBN) NB = 1
                KK = RIVR(1,NB)
                II = RIVR(2,NB)
                JJ = RIVR(3,NB)
C----------DO CALCULATIONS IF THIS IS A MATCH
                IF (I.EQ.II.AND.J.EQ.JJ.AND.K.EQ.KK) THEN
                  IF (IBOUND(J,I,K).LT.1) GOTO 30
                  IFLAG = 1
C-------------ASSIGN VARIABLE VALUES
                  HHNEW = HNEW(J,I,K)
                  HB = RIVR(4,NB)
                  C = RIVR(5,NB)
                  RBOT = RIVR(6,NB)
C-------------CALCULATE SENSITIVITIES
                  XX = ZERO
                  IF (HHNEW.GT.RBOT) XX = -C*SNEW(J,I,K)
                  IF (IIPP.EQ.IBT(2,IQ) .AND. (HHNEW.GT.RBOT))
     &                XX = XX + QCLS(5,N)*(HB-HHNEW)
                  IF (IIPP.EQ.IBT(2,IQ) .AND. HHNEW.LE.RBOT)
     &                XX = QCLS(5,N)*(HB-RBOT)
                  GOTO 20
                ENDIF
   10         CONTINUE
              IF (IFLAG.EQ.0) THEN
                WRITE (IOUT,540) N, NHT + NT, OBSNAM(NHT+NT)
                WRITE (IERRU,540) N, NHT + NT, OBSNAM(NHT+NT)
                IERR = 1
                RETURN
              ENDIF
C-------------SUM VALUES FROM INDIVIDUAL CELLS.
C----------------CALCULATE FACTOR FOR TEMPORAL INTERPOLATION
   20         FACT = 1.0
              IF (TOFF(NHT+NT).GT.ZERO) THEN
                IF (IQOB(NT).EQ.ITS) FACT = 1. - TOFF(NHT+NT)
                IF (IQOB(NT).EQ.ITS-1) FACT = TOFF(NHT+NT)
              ENDIF
C---------------SENSITIVITY-EQUATION SENSITIVITIES
              IF (LN(IIPP).GT.0) XX = XX*B(IIPP)
              X(IP,NHT+NT) = X(IP,NHT+NT) + XX*FACT*QCLS(4,N)
   30       CONTINUE
          ENDIF
   40   CONTINUE
C-------UPDATE COUNTERS
   50   NC = NC + NQCL(IQ)
        NT1 = NT2 + 1
C
   60 CONTINUE
      RETURN
      END
C=======================================================================
      SUBROUTINE OBS1RIV6PR(ITERSS,ITMXP,IUSS,SSRV)
C
C     VERSION 20010613 ERB
C     ******************************************************************
C     WRITE CONTRIBUTION TO SSWR OF FLOW OBSERVATIONS AT RIVER CELLS TO
C     _ss FILE
C     ******************************************************************
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      INTEGER ITMXP, IUSS
      LOGICAL LOP
      DIMENSION SSRV(ITMXP+1)
C     ------------------------------------------------------------------
  660 FORMAT(1X,'ITERATION',2X,A)
  670 FORMAT(1X,I5,6X,G14.7)
C
      INQUIRE(UNIT=IUSS,OPENED=LOP)
      IF (LOP) THEN
        WRITE (IUSS,660)'SSWR-(RIVER-PACKAGE FLOW OBSERVATIONS ONLY)'
C       WRITE CONTRIBUTION TO SSWR FOR EACH ITERATION
        DO 10 IT = 1, ITERSS
          WRITE(IUSS,670) IT,SSRV(IT)
   10   CONTINUE
      ENDIF
C
      RETURN
      END
C=======================================================================
      SUBROUTINE SOBS1RIV6OH(IO,IOWTQRV,IOUT,NHT,NQTRV,HOBS,H,WTQ,
     &                       OBSNAM,IDIS,WTQS,D,AVET,NPOST,NNEGT,NRUNS,
     &                       RSQ,ND,MPR,IPR,NDMH,WTRL,NRSO,IUGDO,OUTNAM,
     &                       IPLOT,IPLPTR,LCOBRIV,ISSWR,SSRV,ITMXP,
     &                       OTIME)
C     VERSION 19990423 ERB
C     ******************************************************************
C     CALCULATE AND PRINT WEIGHTED RESIDUALS FOR RIVER FLOW OBSERVATIONS
C     ******************************************************************
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      REAL AVE, AVET, D, H, HOBS, RES, RSQ, SWH, VMAX, VMIN, WT2, WTQ,
     &     WTQS, WTR, WTRL
      INTEGER IDIS, IO, IOUT, IOWTQRV, IPR,
     &        J, MPR, N, ND, NDMH, NHT, NMAX, NMIN,
     &        NNEG, NNEGT, NPOS, NPOST, NQ1, NQ2, NQTRV, NRSO, NRUNS
      INTEGER IUGDO(6), IPLOT(ND+IPR+MPR), IPLPTR(ND+IPR+MPR)
      CHARACTER*12 OBSNAM(ND)
      CHARACTER*200 OUTNAM
      DIMENSION H(ND), HOBS(ND), D(ND+MPR+IPR), WTQ(NDMH,NDMH),
     &          WTQS(NDMH,NDMH), SSRV(ITMXP+1)
      REAL OTIME(ND)
C     ------------------------------------------------------------------
C
  500 FORMAT (/,' DATA FOR FLOWS REPRESENTED USING THE RIVER PACKAGE',//
     &8X,'OBSERVATION      MEAS.      CALC.',14X,'WEIGHTED',/,
     &'   OBS#    NAME',10X,'FLOW',7X,'FLOW',4X,'RESIDUAL',3X,
     &        'RESIDUAL',/)
  505 FORMAT (/,' DATA FOR FLOWS REPRESENTED USING THE RIVER PACKAGE',//
     &8X,'OBSERVATION',6X,'MEAS.      CALC.',26X,'WEIGHTED',/,
     &3X,'OBS#    NAME',10X,'FLOW',7X,'FLOW',5X,'RESIDUAL',2X,
     &'WEIGHT**.5',2X,'RESIDUAL',/)
  510 FORMAT (1X,I6,1X,A12,1X,5(1X,G10.3))
  515 FORMAT (1X,I6,1X,A12,2X,G10.3,'   DISCONNECTED')
  520 FORMAT (/,' SUM OF SQUARED WEIGHTED RESIDUALS (RIVER FLOWS ONLY)',
     &        2X,G11.5)
  525 FORMAT (/,' STATISTICS FOR RIVER FLOW RESIDUALS :',/,
     &        ' MAXIMUM WEIGHTED RESIDUAL  :',G10.3,' OBS#',I7,/,
     &        ' MINIMUM WEIGHTED RESIDUAL  :',G10.3,' OBS#',I7,/,
     &        ' AVERAGE WEIGHTED RESIDUAL  :',G10.3,/,
     &        ' # RESIDUALS >= 0. :',I7,/,' # RESIDUALS < 0.  :',I7,/,
     &        ' NUMBER OF RUNS: ',I6,'  IN ',I6,' OBSERVATIONS')
  530 FORMAT (2G20.7)
  535 FORMAT (' ')
  540 FORMAT (2(G15.7,1X),I5,2X,A,2X,G15.7)
  550 FORMAT (G15.7,1X,I5,2X,A)
C
      IF (IO.EQ.1) THEN
        IF (IOWTQRV.GT.0) THEN
          WRITE (IOUT,500)
        ELSE
          WRITE (IOUT,505)
        ENDIF
      ENDIF
      IDISRV = 0
      NRESRV = 0
      NRUNSRV = 1
      RSQRV = 0.0
      NNEG = 0
      NPOS = 0
      VMAX = -1.E20
      VMIN = 1.E20
      AVE = 0.
      DO 20 N = LCOBRIV, LCOBRIV+NQTRV-1
        NQ1 = N - NHT
        IF (WTQ(NQ1,NQ1).LT.0.) THEN
          IF (IO.EQ.1) WRITE (IOUT,515) N, OBSNAM(N), HOBS(N)
          IDIS = IDIS + 1
          IDISRV = IDISRV + 1
          GOTO 20
        ENDIF
        NRSO = NRSO + 1
        NRESRV = NRESRV + 1
        IPLPTR(NRSO) = N
        RES = HOBS(N) - H(N)
        IF (IOWTQRV.GT.0) THEN
          WTR = 0.0
          SWH = 0.0
          OWH = 0.0
          DO 10 J = 1, NQTRV
            NQ2 = LCOBRIV + J - 1
            IF (WTQ(NQ2,NQ2).LT.0.0) GOTO 10
            WTR = WTR + WTQS(NQ1,J)*(HOBS(NQ2)-H(NQ2))
            SWH = SWH + WTQS(NQ1,J)*H(NQ2)
            OWH = OWH + WTQS(NQ1,J)*HOBS(NQ2)
   10     CONTINUE
          IF (IO.EQ.1) WRITE (IOUT,510) N, OBSNAM(N),
     &                                  HOBS(N), H(N), RES, WTR
        ELSE
          WT2 = WTQS(NQ1,NQ1)
          WTR = RES*WT2
          OWH = HOBS(N)*WT2
          SWH = H(N)*WT2
          IF (IO.EQ.1) WRITE (IOUT,510) N, OBSNAM(N),
     &                                  HOBS(N), H(N), RES, WT2, WTR
        ENDIF
        IF (IO.EQ.1) THEN
          IF (OUTNAM.NE.'NONE') THEN
            WRITE (IUGDO(1),540) H(N), HOBS(N), IPLOT(N), OBSNAM(N),
     &                           OTIME(N)
            WRITE (IUGDO(2),540) SWH, OWH, IPLOT(N), OBSNAM(N)
            WRITE (IUGDO(3),540) SWH, WTR, IPLOT(N), OBSNAM(N)
            WRITE (IUGDO(4),550) RES, IPLOT(N), OBSNAM(N)
            WRITE (IUGDO(5),550) WTR, IPLOT(N), OBSNAM(N)
            D(NRSO) = WTR
          ENDIF
        ENDIF
        RSQ = RSQ + (WTR**2)
        RSQRV = RSQRV + (WTR**2)
        IF (WTR.GT.VMAX) THEN
          VMAX = WTR
          NMAX = N
        ENDIF
        IF (WTR.LT.VMIN) THEN
          VMIN = WTR
          NMIN = N
        ENDIF
        IF (WTR.GE.0.) NPOS = NPOS + 1
        IF (WTR.LT.0.) NNEG = NNEG + 1
        IF (N.GT.1) THEN
          IF (WTRL*WTR.LT.0.) NRUNS = NRUNS + 1
        ENDIF
        IF (N.GT.LCOBRIV) THEN
          IF (WTRL*WTR.LT.0.) NRUNSRV = NRUNSRV + 1
        ENDIF
        WTRL = WTR
        AVE = AVE + WTR
   20 CONTINUE
      IF (ISSWR.GT.0) SSRV(ISSWR) = RSQRV
      IF (NQTRV.NE.IDISRV) THEN
        AVET = AVET + AVE
        NPOST = NPOST + NPOS
        NNEGT = NNEGT + NNEG
        AVE = AVE/REAL(NQTRV-IDISRV)
        IF (IO.EQ.1) THEN
          WRITE (IOUT,525) VMAX, NMAX, VMIN, NMIN, AVE, NPOS, NNEG,
     &                     NRUNSRV, NRESRV
          WRITE (IOUT,520) RSQRV
        ENDIF
      ENDIF
      RETURN
      END
C=======================================================================
      SUBROUTINE SOBS1RIV6QC(IBT,IIPP,NRIVVL,MXRIVR,NQ,NQCL,RIVR,QCLS,
     &                       NQAR,NQCAR)
C     VERSION 20011106 ERB
C     ******************************************************************
C     POPULATE QCLS ARRAY ELEMENT 5 FOR RIV FLOW-OBSERVATION CELLS ON
C     ONE RIV-PARAMETER-CONTROLLED HEAD-DEPENDENT BOUNDARY IF IT IS
C     ACTIVE THIS STRESS PERIOD
C     ******************************************************************
C     Modified 11/6/2001 to support parameter instances - ERB
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      DIMENSION IBT(2,NQAR), NQCL(NQAR), RIVR(NRIVVL,MXRIVR),
     &          QCLS(5,NQCAR)
      INCLUDE 'param.inc'
C     ------------------------------------------------------------------
      MCOLFAC = 5
      NC = 0
C     DEFINE NUMBER OF INSTANCES (NUMINST), NUMBER OF LIST ENTRIES PER
C     INSTANCE (NLST), AND, FOR THE ACTIVE INSTANCE, INSTANCE NUMBER
C     (NI) AND FIRST AND LAST POSITIONS IN PACKAGE ARRAY (IPL1 and IPL2)
      NI = IACTIVE(IIPP)
      NUMINST = IPLOC(3,IIPP)
      NLST = IPLOC(2,IIPP)-IPLOC(1,IIPP)+1
      IF (NUMINST.GT.1) NLST = NLST/NUMINST
      IPL1 = IPLOC(1,IIPP)+(NI-1)*NLST
      IPL2 = IPL1+NLST-1
C-------LOOP THROUGH ALL CELL GROUPS
      DO 120 IQ=1,NQ
C       IF THIS CELL GROUP INCLUDES A CELL ON BOUNDARY CONTROLLED BY
C       PARAMETER OF INTEREST, IFLAG WILL BE SET = 1
        IFLAG = 0
C---------DETERMINE STARTING AND ENDING CELLS IN THIS CELL GROUP
        NC1 = NC + 1
        NC2 = NC + NQCL(IQ)
        IF (NI.GT.0) THEN
C---------IF PID FOR THIS PARAMETER MATCHES BOUNDARY TYPE SPECIFIED
C         FOR THIS CELL GROUP
          IF (IBT(1,IQ).EQ.1) THEN
C-----------LOOP THROUGH CELLS IN THIS GROUP
            DO 100 JJ = NC1,NC2
              KO = QCLS(1,JJ)
              IO = QCLS(2,JJ)
              JO = QCLS(3,JJ)
C
C-------------LOOP THROUGH PARAMETER (INSTANCE) CELLS IN PACKAGE ARRAY
              DO 80 II = IPL1,IPL2
                K = RIVR(1,II)
                I = RIVR(2,II)
                J = RIVR(3,II)
C---------------IF OBSERVATION CELL IS ON THE PARAMETER-CONTROLLED
C               BOUNDARY, PUT PARAMETER FACTOR IN QCLS(5)
                IF (KO.EQ.K .AND. IO.EQ.I .AND. JO.EQ.J) THEN
                  IBT(2,IQ) = IIPP
                  IFLAG = 1
                  QCLS(5,JJ) = RIVR(MCOLFAC,II)
                  GOTO 100
                ENDIF
   80         CONTINUE
   90         CONTINUE
  100       CONTINUE
          ENDIF
        ENDIF
        IF (IFLAG.EQ.0 .AND. IBT(2,IQ).EQ.IIPP) IBT(2,IQ) = 0
        NC = NC2
  120 CONTINUE
C
      RETURN
      END

max_stars_repo_name: quanpands/wflow

max_issues_repo_name: quanpands/wflow

max_forks_repo_name: quanpands/wflow


text:
      SUBROUTINE DERQP(Q,CMX,ENERGY,P,CMV,CHTIME,DQ,DX,DE,DP,DV,DT)
*
*
*       Derivatives of chain variables.
*       -------------------------------
*
      INCLUDE 'commonc.h'
      LOGICAL  KSLOW,KCOLL
      REAL*8  KSCH
      COMMON/SLOW1/   TK2(0:NMX),EJUMP,KSCH(NMX),KSLOW,KCOLL
      COMMON/CPERT/  RGRAV,GPERT,IPERT,NPERT
      COMMON/CALLS/  TPR,TKPR,STEP,IDER,ICALL,NFN,NREG,ITER,IMCIRC
      COMMON/CCOLL2/  QK(NMX4),PK(NMX4),RIJ(NMX,NMX),SIZE(NMX),VSTAR1,
     &                ECOLL1,RCOLL,QPERI,ISTAR(NMX),ICOLL,ISYNC,NDISS1
      COMMON/BSSAVE/  EP(4),DSC,FACM,TFAC,ITFAC,JC
      COMMON/EBSAVE/  EBS
      COMMON/KSAVE/  K10,K20
      REAL*8  Q(NMX4),CMX(3),P(NMX4),CMV(3)
      REAL*8  DQ(NMX4),DX(3),DP(NMX4),DV(3),FCM(3),UPR(NMX4)
      REAL*8  W(NMX4),AK(NMX4),DK(NMX),FNC(NMX3),FXTNL(NMX3)
      REAL*8  TP(NMX4),TQ(NMX4),UQ(NMX4),FAUX(4),AUX(-2:2),XAUX(3)
      INCLUDE "mpif.h"
      INTEGER group,rank,ierr,isize,status(MPI_STATUS_SIZE)
      COMMON/MPIDAT/group,rank,ierr,isize,status
*
*
*       Mikkola & Aarseth 1990 eqs. (77) -> (80).
      UC=0.0
      RSUM=0.0
*       Eq. (77).
      DO I=1,N-1
      KS = 4*(I-1)
      K1 = KS + 1
      K2 = KS + 2
      K3 = KS + 3
      K4 = KS + 4
*       Obtain W = L^T(Q)P ; L^T = transpose of L times P.
      W(K1)=(Q(K1)*P(K1)-Q(K2)*P(K2)-Q(K3)*P(K3)+Q(K4)*P(K4))
      W(K2)=(Q(K2)*P(K1)+Q(K1)*P(K2)-Q(K4)*P(K3)-Q(K3)*P(K4))
      W(K3)=(Q(K3)*P(K1)+Q(K4)*P(K2)+Q(K1)*P(K3)+Q(K2)*P(K4))
      W(K4)=(Q(K4)*P(K1)-Q(K3)*P(K2)+Q(K2)*P(K3)-Q(K1)*P(K4))
      RIJL=Q(K1)**2+Q(K2)**2+Q(K3)**2+Q(K4)**2
*       Form RSUM for decision-making.
      RSUM=RSUM+RIJL
      RINV(I)=1./RIJL
      A=.5D0*RINV(I)
      UC=UC+MKK(I)*RINV(I)
      DO K=1,4
      W(KS+K)=A*W(KS+K)
      END DO
      END DO
      LRI=N-1
*       Eq. (78).
      TKIN=0.0
      DO I=1,N-1
      AUX(-2)=0.5D0*TK2(I-1)
      AUX(-1)=0.5D0*TK1(I)
      AUX( 0)=TKK(I)
      AUX(+1)=0.5D0*TK1(I+1)
      AUX(+2)=0.5D0*TK2(I+1)
      L=4*(I-1)
      DK(I)=0.0
      DO K=1,4
      AA=0.0
      DO J=-2,2
      LJ=L+4*J
      if (lj.ge.0.and.lj.le.4*n-8) then
      AA=AA+AUX(J)*W(LJ+K)
      end if
      END DO
      AK(L+K)=AA
*       Eq. (79).
      DK(I)=DK(I)+AA*W(L+K)
      END DO
*       Eq. (80).
      TKIN=TKIN+DK(I)
      END DO
*
*       Obtain physical coordinates.
      DO K=1,3
      XI(K)=0.0
      END DO
      DO I=1,N-1
      L=3*(I-1)
      KS=4*(I-1)
      XC(L+1)=Q(KS+1)**2-Q(KS+2)**2-Q(KS+3)**2+Q(KS+4)**2
      XC(L+2)=2.D0*(Q(KS+1)*Q(KS+2)-Q(KS+3)*Q(KS+4))
      XC(L+3)=2.D0*(Q(KS+1)*Q(KS+3)+Q(KS+2)*Q(KS+4))
      DO K=1,3
      XI(L+3+K)=XI(L+K)+XC(L+K)
      END DO
      END DO
*
*       Evaluate external force (FXTNL = W', not the 'usual' force!).
      IF (IPERT.GT.0) THEN
      ISKIP=0
      CALL XTF(FXTNL,FCM,CMX,CHTIME)
*       Note that perturbation limit must be consistent with chlist.f.
      IF (GPERT.LT.1.0D-06) THEN
      IPERT=0
      ELSE
      IPERT=1
      END IF
      ELSE
      ISKIP=1
      DO I=1,3*(N-1)
      FXTNL(I)=0.0D0
      END DO
      DO I=1,3
      FCM(I)=0.0D0
      END DO
      END IF
*
*	Form non-chained contributions.
      UNC=0.0
      DO I=1,3*(N-1)
      FNC(I)=FXTNL(I)
      END DO
*
      DO I=1,N-2
      LI=3*(I-1)
      DO J=I+2,N
      LJ=3*(J-1)
      RIJ2=0.0
      IF(J.GT.I+2)THEN
      DO K=1,3
      XAUX(K)=XI(LJ+K)-XI(LI+K)
      RIJ2=RIJ2+XAUX(K)**2
      END DO
      ELSE
      DO K=1,3
      XAUX(K)=XC(LI+K)+XC(LI+K+3)
      RIJ2=RIJ2+XAUX(K)**2
      END DO
      END IF
      RIJ2IN=1./RIJ2
*       Introduce the inverse distances.
      LRI=LRI+1
      RINV(LRI)=SQRT(RIJ2IN)
*
      FM=MIJ(I,J)*RINV(LRI)
      UNC=UNC+FM
      FM=FM*RIJ2IN
*       Fij attraction.
      DO K=1,3
      FAUX(K)=-FM*XAUX(K)
      END DO
*       Add the contribution to interactions depending on Rij.
      DO IK=I,J-1
      L=3*(IK-1)
      DO K=1,3
      FNC(L+K)=FNC(L+K)+FAUX(K)
      END DO
      END DO
      END DO
      END DO
*
*	Evaluate UQ & TP.
      DO I=1,N-1
      L1=3*(I-1)+1
      KS=4*(I-1)
      KS1=KS+1
      CALL QFORCE(Q(KS1),FNC(L1),UQ(KS1))
      CALL VECTOR(Q(KS1),AK(KS1),TP(KS1))
*       The * operation of eq. (84).
      AK(KS+4)=-AK(KS+4)
      CALL VECTOR(P(KS1),AK(KS1),TQ(KS1))
*
      DO K=1,4
      UQ(KS+K)=UQ(KS+K)-2.0D0*MKK(I)*Q(KS+K)*RINV(I)**2
      TQ(KS+K)=TQ(KS+K)-4.D0*DK(I)*Q(KS+K)
      END DO
      END DO
*	NOTE: The division by R above (in TP & TQ) is delayed.
*
*	Proceed to final evaluation of derivatives (90)->(94).
      UPOT=UC+UNC
      G=1./(TKIN+UPOT)
      H=TKIN-UPOT
*       Reset EJUMP after significant slow-down change (Seppo's procedure).
*     IF (KJUMP) THEN
*         EJUMP = H - ENERGY
*         KJUMP = .false.
*     END IF
      GAMMA=(H-(ENERGY+EJUMP))*G
*
      GT= (1.-GAMMA)*G
      GU=-(1.+GAMMA)*G
*
      DO I=1,N-1
      KS=4*(I-1)
*       Apply the division by R here (to TP & TQ).
      GToverR=GT*RINV(I)
*       NOTE: TP & TQ never get 'correct' values (thus TP = R*TPtrue).
      DO K=1,4
      DQ(KS+K)=GToverR*TP(KS+K)
      DP(KS+K)=-GToverR*TQ(KS+K)-GU*UQ(KS+K)
      END DO
      END DO
      DT=G
      TPR=G
      DO K=1,3
      DX(K)=CMV(K)*G
      DV(K)=FCM(K)*G
      END DO
*	Evaluate E'.
      DE=0.0
      IF (ISKIP.EQ.0) THEN
      DO I=1,N-1
      L1=3*(I-1)+1
      KS=4*(I-1)
      KS1=KS+1
      CALL QFORCE(Q(KS1),FXTNL(L1),FAUX)
      DO K=1,4
      DE=DE+DQ(KS+K)*FAUX(K)
      END DO
      END DO
      END IF
*
*       Copy the time derivative for step control and reset indicator.
*     IF (IDER.GT.0) THEN
*         TPR = G
*         IDER = 0
*     END IF
*
*       Check osculating pericentre of closest pair (first call only).
      IF (ICALL.EQ.0) GO TO 50
      IF (JC.GT.0) GO TO 10
*
*       Perform a fast pericentre search (saves unnecessary complications).
      RM = 0.0
      DO 5 I = 1,N-1
          IF (RINV(I).GT.RM) THEN
              RM = RINV(I)
              IM = I
          END IF
    5 CONTINUE
*
      K1 = INAME(IM)
      K2 = INAME(IM+1)
      K = 4*(IM - 1) + 1
      CALL PERI(Q(K),DQ(K),TPR/KSCH(IM),M(K1),M(K2),QPERI)
*
*       Switch off indicator and exit for large pericentre.
      IF (QPERI.GT.4.0*MAX(SIZE(K1),SIZE(K2))) THEN
          ICALL = 0
          GO TO 50
      END IF
*
*       Examine all close pair indices K1 & K2 in chain vector INAME.
   10 QPMIN = 100.0
      RPMIN = 100.0
      IM0 = 1
      DO 20 IM = 1,N-1
          K1 = INAME(IM)
          K2 = INAME(IM+1)
*
*       Set distance of nearest perturber (note extra test for N > 4).
          IF (IM.EQ.1) THEN
              RP = 1.0/RINV(2)
          ELSE IF (IM.EQ.2.AND.N.GT.3) THEN
              RP = MIN(1.0/RINV(1),1.0/RINV(3))
          ELSE IF (IM.EQ.3.AND.N.GT.4) THEN
              RP = MIN(1.0/RINV(2),1.0/RINV(4))
          ELSE
              RP = 1.0/RINV(IM-1)
          END IF
          RPMIN = MIN(RP,RPMIN)
*
*       Determine pericentre for small perturbations by Mikkola's algorithm.
          GI = (1.0/(RINV(IM)*RP))**3
          IF (GI.LT.0.005) THEN
              K = 4*(IM - 1) + 1
              CALL PERI(Q(K),DQ(K),TPR/KSCH(IM),M(K1),M(K2),QPERI)
          ELSE
              QPERI = 1.0/RINV(IM)
          END IF
*
*       Compare pericentre with previous mutual distances (note symmetry).
          RIJ(K1,K2) = MIN(RIJ(K1,K2),QPERI)
          RIJ(K2,K1) = MIN(RIJ(K2,K1),QPERI)
*
*       Save indices for smallest pericentre and switch off indicator.
          IF (QPERI.LT.QPMIN.AND.IM.NE.IMCIRC) THEN
              G0 = GI
              QPMIN = QPERI
              RP0 = RP
              K10 = K1
              K20 = K2
              IM0 = IM
              ICALL = 0
          END IF
   20 CONTINUE
*
*       Check smallest pericentre of current chain and copy latest value.
      RCOLL = MIN(RCOLL,QPMIN)
      QPERI = QPMIN
*
*       Specify KS index and closest pair indices.
      IM = IM0
      KS1 = 4*(IM - 1) + 1
      K1 = K10
      K2 = K20
*
*       Save TPR and current configuration (for EREL, TRANSK & CHAIN).
      TKPR = TPR
      DO 25 I = 1,N-1
          KS = 4*(I - 1)
          DO 24 J = 1,4
              QK(KS+J) = Q(KS+J)
              PK(KS+J) = P(KS+J)
   24     CONTINUE
   25 CONTINUE
*
*       Check for tidal two-body interaction or stellar collision.
      IF (QPMIN.LT.4.0*MAX(SIZE(K1),SIZE(K2)).OR.ITER.GT.0) THEN
*
*       Terminate after 25 iterations (convergence problem).
          ITER = ITER + 1
          IF (ITER.GE.25) THEN
              if(rank.eq.0)WRITE (6,26)  ITER, IM, G0, QPMIN, ECC
   26         FORMAT (' WARNING!    NO CONVERGENCE    # IM GI QP ECC ',
     &                                                I5,I4,1P,3E10.2)
              JC = 0
              KCOLL = .false.
              ICOLL = -1
              ITER = 0
              IMCIRC = 0
              GO TO 50
          END IF
*
*       Exit if collision candidate distance is not the smallest.
          RB = 1.0/RINV(IM)
*         IF (RB.GT.1.001*RPMIN) THEN
*             JC = 0
*             DSC = 1.0
*             GO TO 50
*         END IF
*
*       Convert Q' to standard KS with T' = R and form radial velocity R'.
          RPR = 0.0D0
          DO 30 J = KS1,KS1+3
              UPR(J) = DQ(J)*RB*KSCH(IM)/TPR
              RPR = RPR + 2.0D0*Q(J)*UPR(J)
   30     CONTINUE
*
*       Determine small semi-major axis from non-singular expressions.
          CALL EREL(IM,EB,SEMI)
*
*       Exclude circularized binaries and consider second smallest peri.
          ECC = 1.0 - QPMIN/SEMI
*       Skip exclusion in GR case (denoted by large VSTAR1).
          IF (ECC.LT.0.002.AND.IMCIRC.EQ.0.AND.VSTAR1.LT.100.0) THEN
              ITER = 0
              IMCIRC = IM
              GO TO 10
          END IF
*
*       Temporary exit because SEMI < 0 does not converge (tested OK 3/99).
*         IF (SEMI.LT.0.0) THEN
*             JC = 0
*             GO TO 50
*         END IF
*
*       Obtain pericentre time interval from elements & Kepler's equation.
          MB = M(K1) + M(K2)
          CALL TPERI(SEMI,Q(KS1),UPR(KS1),MB,TPER)
*
*       Activate collision indicator & B-S step selector (first time).
          IF (ICOLL.EQ.0) THEN
              ICOLL = -1
              JC = 1
              KCOLL = .true.
          ELSE
*
*       Check convergence: radial velocity < 1.0E-09 parabolic velocity.
              IF (ABS(RPR).LT.1.0E-09*SQRT(2.0D0*MB*RB)) THEN
*       Reset B-S step selector and copy chain index to collision indicator.
                  JC = 0
                  ICOLL = IM
                  EBS = EB
                  ITER = 0
                  IMCIRC = 0
              END IF
          END IF
*
*       Set regularized step for DIFSY1 using accelerated iteration rate.
*         IF (RB.GT.2.0*QPMIN) THEN
*             DSC = 2.0*ABS(TPER)/TPR
*         ELSE
*             DSC = ABS(TPER)/TPR
*         END IF
*
*       Evaluate regularized pericentre time (Stiefel & Scheifele, p. 85).
          HI = -0.5*MB/SEMI
*       Note use of Seppo's sign convention (TPER > 0 after peri).
          DSC = 2.0D0*(HI*TPER - 0.5D0*RPR)/MB
*       Scale from KS to chain time derivative.
          DSC = DSC*RB/TPR
*
*       Ensure negative step beyond pericentre (case RPR > 0 works OK).
          IF (JC.GT.0.AND.RPR.LT.0.0D0) THEN
              DSC = ABS(DSC)
          END IF
*       Restore step to dummy value at the end (not used).
          IF (JC.EQ.0) DSC = 1.0
*
*       Switch off iteration on large perturbation after five tries.
          IF (G0.GT.0.005.AND.ITER.GT.5) THEN
              JC = 0
              KCOLL = .false.
              ICOLL = -1
              ITER = 0
              IMCIRC = 0
*       Avoid apocentre region of secondary binary (algorithmic confusion).
          ELSE IF (RB.GT.SEMI.AND.IMCIRC.GT.0) THEN
              JC = 0
              ITER = 0
              IMCIRC = 0
              ICOLL = 0
              KCOLL = .false.
          END IF
      ELSE
          JC = 0
*       Enforce restart from saved variables on failed iteration.
          KCOLL = .false.
          ICOLL = -1
          ICALL = 0
          ITER = 0
      END IF
*
*       Increase function call counter.
   50 NFN = NFN + 1
*
      RETURN
*
      END

ext: f

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
!-----------------------------------------------------------------------
! Skeleton 2-1/2D Electromagnetic Vector PIC code
! written by Viktor K. Decyk, UCLA and Ricardo Fonseca, ISCTE
      program vbpic2
      use iso_c_binding
      use sselib2_c
      use ssebpush2_c
      use vbpush2_h
      implicit none
! indx/indy = exponent which determines grid points in x/y direction:
! nx = 2**indx, ny = 2**indy.
      integer, parameter :: indx =   9, indy =   9
! npx/npy = number of electrons distributed in x/y direction.
      integer, parameter :: npx =  3072, npy =   3072
! ndim = number of velocity coordinates = 3
      integer, parameter :: ndim = 4
! tend = time at end of simulation, in units of plasma frequency.
! dt = time interval between successive calculations.
! qme = charge on electron, in units of e.
      real, parameter :: tend = 10.0, dt = 0.04, qme = -1.0
! vtx/vty = thermal velocity of electrons in x/y direction
! vx0/vy0 = drift velocity of electrons in x/y direction.
      real, parameter :: vtx = 1.0, vty = 1.0, vx0 = 0.0, vy0 = 0.0
! vtx/vz0 = thermal/drift velocity of electrons in z direction
      real, parameter :: vtz = 1.0, vz0 = 0.0
! ax/ay = smoothed particle size in x/y direction
! ci = reciprocal of velocity of light.
      real :: ax = .912871, ay = .912871, ci = 0.1
! idimp = number of particle coordinates = 5
! ipbc = particle boundary condition: 1 = periodic
! sortime = number of time steps between standard electron sorting
! relativity = (no,yes) = (0,1) = relativity is used
      integer :: idimp = 5, ipbc = 1, sortime = 50, relativity = 1
! wke/we = particle kinetic/electrostatic field energy
! wf/wm/wt = magnetic field/transverse electric field/total energy
      real, target :: wke = 0.0, we = 0.0, wf = 0.0, wm = 0.0
      real :: wt = 0.0
! kvec = (1,2) = run (autovector,SSE2) version
      integer :: kvec = 1
!
! declare scalars for standard code
      integer :: np, nx, ny, nxh, nyh, nxe, nye, nxeh, nxyh, nxhy
      integer :: npe, ny1, ntime, nloop, isign
      integer :: irc = 0
      real :: qbme, affp, dth
!
! declare arrays for standard code:
! partt, partt2 = transposed particle arrays
      real, dimension(:,:), pointer :: partt => null()
      real, dimension(:,:), pointer :: partt2 => null()
      real, dimension(:,:), pointer :: tpartt => null()
! qe = electron charge density with guard cells
      real, dimension(:,:), pointer :: qe => null()
! cue = electron current density with guard cells
      real, dimension(:,:,:), pointer :: cue => null()
! fxyze/bxyze = smoothed electric/magnetic field with guard cells
      real, dimension(:,:,:), pointer :: fxyze => null()
      real, dimension(:,:,:), pointer :: bxyze => null()
! exyz/bxyz = transverse electric/magnetic field in fourier space
      complex, dimension(:,:,:), pointer :: exyz => null()
      complex, dimension(:,:,:), pointer :: bxyz => null()
! ffc = form factor array for poisson solver
      complex, dimension(:,:), pointer :: ffc => null()
! mixup = bit reverse table for FFT
      integer, dimension(:), pointer :: mixup => null()
! sct = sine/cosine table for FFT
      complex, dimension(:), pointer :: sct => null()
! npicy = scratch array for reordering particles
      integer, dimension(:), pointer :: npicy => null()
!
! declare and initialize timing data
      real :: time
      type, bind(C) :: timeval
         integer(c_long) :: tv_sec
         integer(c_long) :: tv_usec
      end type
      type (timeval) :: itime
      double precision :: dtime
      real :: tdpost = 0.0, tguard = 0.0, tfft = 0.0, tfield = 0.0
      real :: tdjpost = 0.0, tpush = 0.0, tsort = 0.0
!
      interface
         subroutine dtimer(time,itime,icntrl) bind(C,name='dtimer')
         use iso_c_binding
         import :: timeval
         implicit none
         real(c_double) :: time
         type (timeval) :: itime
         integer(c_int), value :: icntrl
         end subroutine
      end interface
!
! initialize scalars for standard code
! np = total number of particles in simulation
! nx/ny = number of grid points in x/y direction
      np = npx*npy; nx = 2**indx; ny = 2**indy; nxh = nx/2; nyh = ny/2
      nxe = nx + 2; nye = ny + 1; nxeh = nxe/2
      nxyh = max(nx,ny)/2; nxhy = max(nxh,ny); ny1 = ny + 1
! nloop = number of time steps in simulation
! ntime = current time step
      nloop = tend/dt + .0001; ntime = 0
      qbme = qme
      affp = real(nx*ny)/real(np)
      dth = 0.0
!
! allocate data for standard code
      allocate(mixup(nxhy),sct(nxyh))
!
! align memory for SSE
      npe = 4*((np - 1)/4 + 1)
      nxe = 4*((nxe - 1)/4 + 1)
      nxeh = nxe/2
      call sse_f2allocate(partt,npe,idimp,irc)
      if (sortime > 0) then
         call sse_f2allocate(partt2,npe,idimp,irc)
      endif
      call sse_f2allocate(qe,nxe,nye,irc)
      call sse_f3allocate(cue,ndim,nxe,nye,irc)
      call sse_f3allocate(fxyze,ndim,nxe,nye,irc)
      call sse_f3allocate(bxyze,ndim,nxe,nye,irc)
      call sse_c3allocate(exyz,ndim,nxeh,nye,irc)
      call sse_c3allocate(bxyz,ndim,nxeh,nye,irc)
      call sse_c2allocate(ffc,nxh,nyh,irc)
      call sse_i1allocate(npicy,ny1,irc)
      if (irc /= 0) then
         write (*,*) 'aligned allocation error: irc = ', irc
      endif
!
! prepare fft tables
      call WFFT2RINIT(mixup,sct,indx,indy,nxhy,nxyh)
! calculate form factors
      isign = 0
      call VPOIS23(qe,fxyze,isign,ffc,ax,ay,affp,we,nx,ny,nxeh,nye,nxh, &
     &nyh)
! initialize electrons
      call DISTR2HT(partt,vtx,vty,vtz,vx0,vy0,vz0,npx,npy,idimp,npe,nx, &
     &ny,ipbc)
!
! initialize transverse electromagnetic fields
      exyz = cmplx(0.0,0.0)
      bxyz = cmplx(0.0,0.0)
!
      if (dt > 0.45*ci) then
         write (*,*) 'Warning: Courant condition may be exceeded!'
      endif
!
! * * * start main iteration loop * * *
!
  500 if (nloop <= ntime) go to 2000
!     write (*,*) 'ntime = ', ntime
!
! deposit current with standard procedure: updates part, cue
      call dtimer(dtime,itime,-1)
      cue = 0.0
      if (relativity==1) then
         if (kvec==1) then
            call VGRJPOST2LT(partt,cue,qme,dth,ci,np,npe,idimp,nx,ny,nxe&
     &,nye,ipbc)
! SSE2 function
         else if (kvec==2) then
            call csse2grjpost2lt(c_loc(partt(1,1)),c_loc(cue(1,1,1)),qme&
     &,dth,ci,np,npe,idimp,nx,ny,nxe,nye,ipbc)
         endif
      else
         if (kvec==1) then
            call VGJPOST2LT(partt,cue,qme,dth,np,npe,idimp,nx,ny,nxe,nye&
     &,ipbc)
! SSE2 function
         else if (kvec==2) then
            call csse2gjpost2lt(c_loc(partt(1,1)),c_loc(cue(1,1,1)),qme,&
     &dth,np,npe,idimp,nx,ny,nxe,nye,ipbc)
         endif
      endif
      call dtimer(dtime,itime,1)
      time = real(dtime)
      tdjpost = tdjpost + time
!
! deposit charge with standard procedure: updates qe
      call dtimer(dtime,itime,-1)
      qe = 0.0
      if (kvec==1) then
         call VGPOST2LT(partt,qe,qme,np,npe,idimp,nxe,nye)
! SSE2 function
      else if (kvec==2) then
         call csse2gpost2lt(c_loc(partt(1,1)),c_loc(qe(1,1)),qme,np,npe,&
     &idimp,nxe,nye)
      endif
      call dtimer(dtime,itime,1)
      time = real(dtime)
      tdpost = tdpost + time
!
! add guard cells with standard procedure: updates cue, qe
      call dtimer(dtime,itime,-1)
      if (kvec==1) then
         call ACGUARD2L(cue,nx,ny,nxe,nye)
         call AGUARD2L(qe,nx,ny,nxe,nye)
! SSE2 function
      else if (kvec==2) then
         call csse2acguard2l(c_loc(cue(1,1,1)),nx,ny,nxe,nye)
         call csse2aguard2l(c_loc(qe(1,1)),nx,ny,nxe,nye)
      endif
      call dtimer(dtime,itime,1)
      time = real(dtime)
      tguard = tguard + time
!
! transform charge to fourier space with standard procedure: updates qe
      call dtimer(dtime,itime,-1)
      isign = -1
      if (kvec==1) then
         call WFFT2RVX(qe,isign,mixup,sct,indx,indy,nxeh,nye,nxhy,nxyh)
! SSE2 function
      else if (kvec==2) then
         call csse2wfft2rx(c_loc(qe(1,1)),isign,c_loc(mixup(1)),        &
     &c_loc(sct(1)),indx,indy,nxeh,nye,nxhy,nxyh)
      endif
      call dtimer(dtime,itime,1)
      time = real(dtime)
      tfft = tfft + time
!
! transform current to fourier space with standard procedure: update cue
      call dtimer(dtime,itime,-1)
      isign = -1
      if (kvec==1) then
         call WFFT2RV3(cue,isign,mixup,sct,indx,indy,nxeh,nye,nxhy,nxyh)
! SSE2 function
      else if (kvec==2) then
         call csse2wfft2r3(c_loc(cue(1,1,1)),isign,c_loc(mixup(1)),     &
     &c_loc(sct(1)),indx,indy,nxeh,nye,nxhy,nxyh)
      endif
      call dtimer(dtime,itime,1)
      time = real(dtime)
      tfft = tfft + time
!
! take transverse part of current with standard procedure: updates cue
      call dtimer(dtime,itime,-1)
      if (kvec==1) then
         call CUPERP2(cue,nx,ny,nxeh,nye)
! SSE2 function
      else if (kvec==2) then
         call csse2cuperp2(c_loc(cue(1,1,1)),nx,ny,nxeh,nye)
      endif
      call dtimer(dtime,itime,1)
      time = real(dtime)
      tfield = tfield + time
!
! calculate electromagnetic fields in fourier space with standard
! procedure: updates exyz, bxyz
      call dtimer(dtime,itime,-1)
      if (ntime==0) then
         if (kvec==1) then
            call VIBPOIS23(cue,bxyz,ffc,ci,wm,nx,ny,nxeh,nye,nxh,nyh)
! SSE2 function
         else if (kvec==2) then
            call csse2ibpois23(c_loc(cue(1,1,1)),c_loc(bxyz(1,1,1)),    &
     &c_loc(ffc(1,1)),ci,c_loc(wm),nx,ny,nxeh,nye,nxh,nyh)
         endif
         wf = 0.0
         dth = 0.5*dt
      else
         if (kvec==1) then
            call VMAXWEL2(exyz,bxyz,cue,ffc,ci,dt,wf,wm,nx,ny,nxeh,nye, &
     &nxh,nyh)
! SSE2 function
         else if (kvec==2) then
            call csse2maxwel2(c_loc(exyz(1,1,1)),c_loc(bxyz(1,1,1)),    &
     &c_loc(cue(1,1,1)),c_loc(ffc(1,1)),ci,dt,c_loc(wf),c_loc(wm),nx,ny,&
     &nxeh,nye,nxh,nyh)
         endif
      endif
      call dtimer(dtime,itime,1)
      time = real(dtime)
      tfield = tfield + time
!
! calculate force/charge in fourier space with standard procedure:
! updates fxyze
      call dtimer(dtime,itime,-1)
      isign = -1
      if (kvec==1) then
         call VPOIS23(qe,fxyze,isign,ffc,ax,ay,affp,we,nx,ny,nxeh,nye,  &
     &nxh,nyh)
! SSE2 function
      else if (kvec==2) then
         call csse2pois23(c_loc(qe(1,1)),c_loc(fxyze(1,1,1)),isign,     &
     &c_loc(ffc(1,1)),ax,ay,affp,c_loc(we),nx,ny,nxeh,nye,nxh,nyh)
      endif
      call dtimer(dtime,itime,1)
      time = real(dtime)
      tfield = tfield + time
!
! add longitudinal and transverse electric fields with standard
! procedure: updates fxyze
      call dtimer(dtime,itime,-1)
      isign = 1
      if (kvec==1) then
         call VEMFIELD2(fxyze,exyz,ffc,isign,nx,ny,nxeh,nye,nxh,nyh)
! SSE2 function
      else if (kvec==2) then
         call csse2emfield2(c_loc(fxyze(1,1,1)),c_loc(exyz(1,1,1)),     &
     &c_loc(ffc(1,1)),isign,nx,ny,nxeh,nye,nxh,nyh)
      endif
! copy magnetic field with standard procedure: updates bxyze
      isign = -1
      if (kvec==1) then
         call VEMFIELD2(bxyze,bxyz,ffc,isign,nx,ny,nxeh,nye,nxh,nyh)
! SSE2 function
      else if (kvec==2) then
         call csse2emfield2(c_loc(bxyze(1,1,1)),c_loc(bxyz(1,1,1)),     &
     &c_loc(ffc(1,1)),isign,nx,ny,nxeh,nye,nxh,nyh)
      endif
      call dtimer(dtime,itime,1)
      time = real(dtime)
      tfield = tfield + time
!
! transform electric force to real space with standard procedure:
! updates fxyze
      call dtimer(dtime,itime,-1)
      isign = 1
      if (kvec==1) then
         call WFFT2RV3(fxyze,isign,mixup,sct,indx,indy,nxeh,nye,nxhy,   &
     &nxyh)
! SSE2 function
      else if (kvec==2) then
         call csse2wfft2r3(c_loc(fxyze(1,1,1)),isign,c_loc(mixup(1)),   &
     &c_loc(sct(1)),indx,indy,nxeh,nye,nxhy,nxyh)
      endif
      call dtimer(dtime,itime,1)
      time = real(dtime)
      tfft = tfft + time
!
! transform magnetic force to real space with standard procedure:
! updates bxyze
      call dtimer(dtime,itime,-1)
      isign = 1
      if (kvec==1) then
         call WFFT2RV3(bxyze,isign,mixup,sct,indx,indy,nxeh,nye,nxhy,   &
     &nxyh)
! SSE2 function
      else if (kvec==2) then
         call csse2wfft2r3(c_loc(bxyze(1,1,1)),isign,c_loc(mixup(1)),   &
     &c_loc(sct(1)),indx,indy,nxeh,nye,nxhy,nxyh)
      endif
      call dtimer(dtime,itime,1)
      time = real(dtime)
      tfft = tfft + time
!
! copy guard cells with standard procedure: updates fxyze, bxyze
      call dtimer(dtime,itime,-1)
      if (kvec==1) then
         call BGUARD2L(fxyze,nx,ny,nxe,nye)
         call BGUARD2L(bxyze,nx,ny,nxe,nye)
! SSE2 function
      else if (kvec==2) then
         call csse2bguard2l(c_loc(fxyze(1,1,1)),nx,ny,nxe,nye)
         call csse2bguard2l(c_loc(bxyze(1,1,1)),nx,ny,nxe,nye)
      endif
      call dtimer(dtime,itime,1)
      time = real(dtime)
      tguard = tguard + time
!
! push particles with standard procedure: updates part, wke
      wke = 0.0
      call dtimer(dtime,itime,-1)
! GBPUSH23LT and GRBPUSH23LT give different results with ifort -O3
      if (relativity==1) then
         if (kvec==1) then
            call VGRBPUSH23LT(partt,fxyze,bxyze,qbme,dt,dth,ci,wke,idimp&
     &,np,npe,nx,ny,nxe,nye,ipbc)
! SSE2 function
         else if (kvec==2) then
            call csse2grbpush23lt(c_loc(partt(1,1)),c_loc(fxyze(1,1,1)),&
     &c_loc(bxyze(1,1,1)),qbme,dt,dth,ci,c_loc(wke),idimp,np,npe,nx,ny, &
     &nxe,nye,ipbc)
         endif
      else
         if (kvec==1) then
            call VGBPUSH23LT(partt,fxyze,bxyze,qbme,dt,dth,wke,idimp,np,&
     &npe,nx,ny,nxe,nye,ipbc)
! SSE2 function
         else if (kvec==2) then
            call csse2gbpush23lt(c_loc(partt(1,1)),c_loc(fxyze(1,1,1)), &
     &c_loc(bxyze(1,1,1)),qbme,dt,dth,c_loc(wke),idimp,np,npe,nx,ny,nxe,&
     &nye,ipbc)
         endif
      endif
      call dtimer(dtime,itime,1)
      time = real(dtime)
      tpush = tpush + time
!
! sort particles by cell for standard procedure
      if (sortime > 0) then
         if (mod(ntime,sortime)==0) then
            call dtimer(dtime,itime,-1)
            if (kvec==1) then
               call DSORTP2YLT(partt,partt2,npicy,idimp,np,npe,ny1)
! SSE2 function
            else if (kvec==2) then
               call csse2dsortp2ylt(c_loc(partt(1,1)),c_loc(partt2(1,1))&
     &,c_loc(npicy(1)),idimp,np,npe,ny1)
            endif
! exchange pointers
            tpartt => partt
            partt => partt2
            partt2 => tpartt
            call dtimer(dtime,itime,1)
            time = real(dtime)
            tsort = tsort + time
         endif
      endif
!
      if (ntime==0) then
         wt = we + wf + wm
         write (*,*) 'Initial Total Field, Kinetic and Total Energies:'
         write (*,'(3e14.7)') wt, wke, wke + wt
         write (*,*) 'Initial Electrostatic, Transverse Electric and Mag&
     &netic Field Energies:'
         write (*,'(3e14.7)') we, wf, wm
      endif
      ntime = ntime + 1
      go to 500
 2000 continue
!
! * * * end main iteration loop * * *
!
      write (*,*) 'ntime, relativity = ', ntime, relativity
      write (*,*) 'kvec = ', kvec
      wt = we + wf + wm
      write (*,*) 'Final Total Field, Kinetic and Total Energies:'
      write (*,'(3e14.7)') wt, wke, wke + wt
      write (*,*) 'Final Electrostatic, Transverse Electric and Magnetic&
     & Field Energies:'
      write (*,'(3e14.7)') we, wf, wm
!
      write (*,*)
      write (*,*) 'deposit time = ', tdpost
      write (*,*) 'current deposit time = ', tdjpost
      tdpost = tdpost + tdjpost
      write (*,*) 'total deposit time = ', tdpost
      write (*,*) 'guard time = ', tguard
      write (*,*) 'solver time = ', tfield
      write (*,*) 'fft time = ', tfft
      write (*,*) 'push time = ', tpush
      write (*,*) 'sort time = ', tsort
      tfield = tfield + tguard + tfft
      write (*,*) 'total solver time = ', tfield
      time = tdpost + tpush + tsort
      write (*,*) 'total particle time = ', time
      wt = time + tfield
      write (*,*) 'total time = ', wt
      write (*,*)
!
      wt = 1.0e+09/(real(nloop)*real(np))
      write (*,*) 'Push Time (nsec) = ', tpush*wt
      write (*,*) 'Deposit Time (nsec) = ', tdpost*wt
      write (*,*) 'Sort Time (nsec) = ', tsort*wt
      write (*,*) 'Total Particle Time (nsec) = ', time*wt
!
      call sse_deallocate(c_loc(npicy(1)),irc); nullify(npicy)
      call sse_deallocate(c_loc(ffc(1,1)),irc); nullify(ffc)
      call sse_deallocate(c_loc(bxyz(1,1,1)),irc); nullify(bxyz)
      call sse_deallocate(c_loc(exyz(1,1,1)),irc); nullify(exyz)
      call sse_deallocate(c_loc(bxyze(1,1,1)),irc); nullify(bxyze)
      call sse_deallocate(c_loc(fxyze(1,1,1)),irc); nullify(fxyze)
      call sse_deallocate(c_loc(cue(1,1,1)),irc); nullify(cue)
      call sse_deallocate(c_loc(qe(1,1)),irc); nullify(qe)
      if (sortime > 0) then
         call sse_deallocate(c_loc(partt2(1,1)),irc); nullify(partt2)
      endif
      call sse_deallocate(c_loc(partt(1,1)),irc); nullify(partt)
!
      stop
      end program

max_stars_repo_name: gcasabona/cuda

max_issues_repo_name: gcasabona/cuda

max_forks_repo_name: gcasabona/cuda


program test6
#if (_DP==0)
 use iso_fortran_env,only:int32,int64,wp=>real32
#else
 use iso_fortran_env,only:int32,int64,wp=>real64
#endif
 use modsparse
 implicit none
 integer(kind=int32)::nrow
 integer(kind=int32)::row
 integer(kind=int32)::col
 integer(kind=int32)::iunit, istat
 integer(kind=int32)::i,j
 integer(kind=int32),allocatable::iarray(:,:)
 integer(kind=int32),allocatable::perm(:)
 real(kind=wp)::val
 real(kind=wp),allocatable::x(:),y(:)
 logical::lup=.false.
 type(coosparse)::coo
 type(crssparse)::crs

 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !COO UPPER
 open(newunit=iunit,file='matrixija.ascii',status='old',action='read')
 !open(newunit=iunit,file='matkm.dat',status='old',action='read')
 read(iunit,*) nrow

 coo=coosparse(nrow,lupper=.true.)

 do
  read(iunit,*,iostat=istat) row,col,val
  if(istat.ne.0)exit
  call coo%add(row,col,val)
!  if(row.ne.col)call coo%add(col,row,val)
 end do
 close(iunit)

 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !CSR UPPER
 crs=coo

 call crs%printstats()

 call crs%printsquare()

#if (_METIS==1)
 call crs%setpermutation(crs%getordering(bglvl=0))
#else
 call crs%setpermutation((/(i,i=1,crs%getdim(1))/))
#endif
 call crs%printstats()

!!!!!!!!!!!!!!!!
 write(*,*)'Matrix'
 allocate(iarray(crs%getdim(1),crs%getdim(2)))
 iarray=0
 do i=1,crs%getdim(1)
  do j=1,crs%getdim(2)
   val=crs%get(i,j)
   if(val.ne.0.)iarray(i,j)=1
  enddo
 enddo

 do i=1,crs%getdim(1)
  write(*,'(10000(i2))')iarray(i,:)
 enddo

 write(*,*)'Permuted matrix'

#if (_METIS==1)
 perm=crs%getordering(bglvl=0)
#else
 perm=(/(i,i=1,crs%getdim(1))/)
#endif
 iarray=0
 do i=1,crs%getdim(1)
  do j=1,crs%getdim(2)
   val=crs%get(perm(i),perm(j))
   if(val.ne.0.)iarray(i,j)=1
  enddo
 enddo

 do i=1,crs%getdim(1)
  write(*,'(10000(i2))')iarray(i,:)
 enddo

!!!!!!

 call crs%spainv()

 call crs%printsquare()

!!!!!!
 crs=coo

 call crs%sort()

#if (_METIS==1)
 perm=crs%getordering(bglvl=0)
#else
 perm=(/(i,i=1,crs%getdim(1))/)
#endif
 call crs%setpermutation(perm)
 call crs%printstats()

 print*,'ord ',perm

 allocate(x(crs%getdim(1)),y(crs%getdim(1)))
 x=1.;y=1.

 call crs%solve(x,y)

print*,'x ',x
print*,'y ',y

end program


text:
      SUBROUTINE FLYBY(I,ITERM)
*
*
*      Flyby check of KS orbit.
*       -----------------------
*
      INCLUDE 'common6.h'
*
*
*       Set index of KS pair & first component of c.m. body #I.
      IPAIR = I - N
      I1 = 2*IPAIR - 1
      RJMIN2 = 1000.0
      ITERM = 0
      JCL = LIST(2,I1)
*
*       Find the closest body (JCL) and copy perturbers to JLIST.
      NNB = LIST(1,I1)
      DO 10 L = 2,NNB+1
          J = LIST(L,I1)
          JLIST(L-1) = J
          RIJ2 = (X(1,I) - X(1,J))**2 + (X(2,I) - X(2,J))**2 +
     &                                  (X(3,I) - X(3,J))**2
          IF (RIJ2.LT.RJMIN2) THEN
              RJMIN2 = RIJ2
              JCL = J
          END IF
   10 CONTINUE
*
*       Skip rare case of merged binary or chain c.m. (denoted by NAME <= 0).
      IF (NAME(I).LE.0.OR.NAME(JCL).LE.0) GO TO 20
*
      RDOT = (X(1,I) - X(1,JCL))*(XDOT(1,I) - XDOT(1,JCL)) +
     &       (X(2,I) - X(2,JCL))*(XDOT(2,I) - XDOT(2,JCL)) +
     &       (X(3,I) - X(3,JCL))*(XDOT(3,I) - XDOT(3,JCL))
*
      VREL2 = (XDOT(1,I) - XDOT(1,JCL))**2 +
     &        (XDOT(2,I) - XDOT(2,JCL))**2 +
     &        (XDOT(3,I) - XDOT(3,JCL))**2

*       Evaluate outer semi-major axis, eccentricity & impact parameter.
      RIJ = SQRT(RJMIN2)
      SEMI1 = 2.0/RIJ - VREL2/(BODY(I) + BODY(JCL))
      SEMI1 = 1.0/SEMI1
      ECC1 = SQRT((1.0D0 - RIJ/SEMI1)**2 +
     &                     RDOT**2/(SEMI1*(BODY(I) + BODY(JCL))))
      PMIN = SEMI1*(1.0D0 - ECC1)
*
*       Set semi-major axis & eccentricity of inner binary.
      SEMI = -0.5D0*BODY(I)/H(IPAIR)
      ECC2 = (1.0D0 - R(IPAIR)/SEMI)**2 + TDOT2(IPAIR)**2/(BODY(I)*SEMI)
      ECC = SQRT(ECC2)
*
      I2 = I1 + 1
      RJI1 = (X(1,JCL) - X(1,I1))**2 + (X(2,JCL) - X(2,I1))**2 +
     &                                 (X(3,JCL) - X(3,I1))**2
      RJI2 = (X(1,JCL) - X(1,I2))**2 + (X(2,JCL) - X(2,I2))**2 +
     &                                 (X(3,JCL) - X(3,I2))**2
*
*       Identify the dominant interaction (#I1 or #I2).
      FJI1 = (BODY(I1) + BODY(JCL))/RJI1
      FJI2 = (BODY(I2) + BODY(JCL))/RJI2
      NNB = NNB + 1
      IF (FJI1.GT.FJI2) THEN
          J = I1
          RJ2 = RJI1
*       Add the other component to the perturber list.
          JLIST(NNB) = I2
      ELSE
          J = I2
          RJ2 = RJI2
          JLIST(NNB) = I1
      END IF
*
*       Evaluate radial velocity (routine SEARCH requires RD < 0).
      RD = 0.0
      RJJ = 0.0
      VJJ = 0.0
      DO 15 K = 1,3
          RD = RD + (X(K,JCL) - X(K,J))*(XDOT(K,JCL) - XDOT(K,J))
          RJJ = RJJ + (X(K,JCL) - X(K,J))**2
          VJJ = VJJ + (XDOT(K,JCL) - XDOT(K,J))**2
   15 CONTINUE
*
*       Determine vectorial perturbation on intruder and closest component.
      CALL FPERT(JCL,J,NNB,PERT)
      GJ = PERT*RJ2/(BODY(JCL) + BODY(J))
      RJ = SQRT(RJJ)
*     WRITE (6,18) NAME(J),NAME(JCL),GAMMA(IPAIR),GJ,R(IPAIR),RJ,RD,
*    &             RDOT
*  18 FORMAT (' NMJ NMJC GI GJ R RJ RD RDI ',2I5,2F6.2,1P,4E9.1)
*
*       Terminate if perturber & I1/I2 provides dominant force.
      APO = ABS(SEMI)*(1.0 + ECC)
      RSUM = R(IPAIR) + SQRT(RJ2)
      XF = 2.0
*       Compare current radial velocity with #JCL & J.
      IF (TDOT2(IPAIR)/R(IPAIR).LT.RD/RJ) XF = 1.0
*
*       Increase the cross section for quadruples.
      IF (JCL.GT.N) THEN
          SEMI2 = -0.5*BODY(JCL)/H(JCL-N)
          APO = APO + ABS(SEMI2)
      END IF
*
      SEMIJ = 2.0/RJ - VJJ/(BODY(JCL) + BODY(J))
      SEMIJ = 1.0/SEMIJ
*       Check for chain regularization test or standard termination.
      IF (PMIN.LT.APO.AND.RSUM.LT.XF*RMIN.AND.RDOT.LT.0.0) THEN
          IF (SEMI.LT.0.0.AND.GAMMA(IPAIR).LT.0.5) THEN
              ITERM = 0
              GO TO 20
          END IF
          ITERM = 1
      ELSE IF (GJ.LT.0.7*GAMMA(IPAIR).AND.RD.LE.0.0) THEN
          ITERM = 2
      END IF
*
   20 RETURN
*
      END


text:
foreach([],_): 1.

foreach( head <~ [tail], F):
    F(head)
    foreach(tail, F).

max_stars_repo_name: diegovalverde/funk

max_issues_repo_name: diegovalverde/funk

max_forks_repo_name: diegovalverde/funk


text:
module cordriv
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!       NASA/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-----------------------------------------------------------------------
!
! !ROUTINE: cordriv - a F90 module of PSAS' correlation function drivers
!
! !INTERFACE:
!	use cordriv
!
! !DESCRIPTION:
!	cordriv() is a F90 module of PSAS' correlation function driver
!	subroutines.  It defines the interfaces of the following
!	correlation function driver subroutines,
!
!	  diagcorF: forming a diagonal block of the forecast error
!		    correlation matrix;
!	  diagcorO: forming a diagonal block of the observation error
!		    correlation matrix with correlated horizontal
!		    correlation functions;
!	  diagcorU: forming a diagonal block of the observation error
!		    correlation matrix with uncorrelated horizontal
!		    correlation functions;
!	  offdcorF: forming a off-diagonal block of the forecast error
!		    correlation matrix;
!	  offdcorO: forming a off-diagonal block of the observation
!		    error correlation matrix with correlation horizontal
!		    correlation functions.
!
! !BUGS:
!	Several changes are already planned for the interface of the
!	four existent subroutines, as while as a group of new entries
!	to be included in to this group.
!
!	Although it is convenient to use the module to insure the
!	correctness of the calling interface, it does not help much
!	if the actual subroutines/functions are coded different from
!	this interface.  So it is important to make sure the actual
!	interfaces are accurately recorded in the this module.  An
!	automatic process is desireble.
!
! !SEE ALSO:
!	diagcorF(), diagcorO(), diagcorU(), offdcorF(), and offdcorO().
!
! !REVISION HISTORY:
!	25Jan95 - J. Guo	- added diagcorU() interface.  Some
!		corrections also made according to the changes in other
!		include subroutines.
!
! 	15Dec95 - J. Guo	- programed and added the prolog
!_______________________________________________________________________
interface
!=======================================================================
!!$  subroutine diagcorF(	kind_cov,				&
!!$			kt,len,	qr_x,qr_y,qr_z,			&
!!$				qm_x,qm_y,qm_z,ql_x,ql_y,	&
!!$				ktab,				&
!!$			Mtyp,corF,istat				)
  subroutine diagcorF(	kind_cov,				&
			kt, iattr, rattr, &
			Mtyp,corF,istat				)
!-----------------------------------------------------------------------
    integer, intent(in)	:: kind_cov
    integer, intent(in)	:: kt
    real,    intent(in) :: rattr(:,:)
    integer, intent(in) :: iattr(:,:)
    character*1, intent(out)	:: Mtyp
!!$    real,    intent(out)	:: corF(len*(len+1)/2)
    real,    intent(out)	:: corF(size(iattr,2)*(size(iattr,2)+1)/2)
    integer, intent(out)	:: istat
  end subroutine diagcorF
!=======================================================================
  subroutine diagcorO(	kt, iattr, rattr,  &
			Mtyp,corO, istat	)
!-----------------------------------------------------------------------
    integer, intent(in)	:: kt
    real, intent(in)    :: rattr(:,:)
    integer, intent(in) :: iattr(:,:)
    character*1, intent(out)	:: Mtyp
!!$    real,    intent(out)	:: corO(len*(len+1)/2)
    real, intent(out) :: corO(size(iattr,2)*(size(iattr,2)+1)/2)
    integer, intent(out)	:: istat
  end subroutine diagcorO
!=======================================================================
!!$  subroutine diagcorU(	kt,len,	kx,ks,ktab,Mtyp,corU,istat)
  subroutine diagcorU(	kt, iattr, Mtyp, corU, istat)
!-----------------------------------------------------------------------
    integer, intent(in)	:: kt		! data type
    integer, intent(in) :: iattr(:,:)
    character*1, intent(out)	:: Mtyp
!!$    real,    intent(out)	:: corU(len*(len+1)/2)
    real, intent(out) :: corU(size(iattr,2)*(size(iattr,2)+1)/2)
    integer, intent(out)	:: istat
  end subroutine diagcorU
!=======================================================================
!!$  subroutine offdcorF(kind_cov,					&
!!$		      kti,leni,	ktabi, stride_i, coords_i,                &
!!$		      ktj,lenj,	ktabj, stride_j, coords_j,                &
!!$		      Mtyp,corF,istat				)
  subroutine offdcorF(kind_cov,					&
		      kti,rattr_i, iattr_i,                &
		      ktj,rattr_j, iattr_j,                &
		      Mtyp,corF,istat				)
!-----------------------------------------------------------------------
    integer, intent(in)	:: kind_cov

    integer, intent(in)	:: kti
    integer, intent(in) :: iattr_i(:,:)
    real,    intent(in) :: rattr_i(:,:)

    integer, intent(in)	:: ktj
    integer, intent(in) :: iattr_j(:,:)
    real,    intent(in) :: rattr_j(:,:)

    character*1, intent(out)	:: Mtyp
!!$    real,    intent(out)	:: corF(leni*lenj)
    real, intent(out) :: corf(:)
    integer, intent(out)	:: istat
  end subroutine offdcorF
!=======================================================================
!!$  subroutine offdcorO(kti,leni,kxi,qri_x,qri_y,qri_z,	&
!!$				ktabi,				&
!!$		      ktj,lenj,kxj,qrj_x,qrj_y,qrj_z,	&
!!$				ktabj,				&
!!$		      Mtyp,corO,istat				)
  subroutine offdcorO(kti,rattr_i, iattr_i, &
                      ktj,rattr_j, iattr_j, &
		      Mtyp,corO,istat				)
!-----------------------------------------------------------------------
    integer, intent(in)	:: kti
!!$    integer, intent(in)	:: leni
!!$    integer, intent(in)	::  kxi(leni)
!!$    real,    intent(in)	:: qri_x(leni)
!!$    real,    intent(in)	:: qri_y(leni)
!!$    real,    intent(in)	:: qri_z(leni)
!!$    integer, intent(in)	:: ktabi(leni)
    integer, intent(in) :: iattr_i(:,:)
    real, intent(in) :: rattr_i(:,:)

    integer, intent(in)	:: ktj
!!$    integer, intent(in)	:: lenj
!!$    integer, intent(in)	::  kxj(lenj)
!!$    real,    intent(in)	:: qrj_x(lenj)
!!$    real,    intent(in)	:: qrj_y(lenj)
!!$    real,    intent(in)	:: qrj_z(lenj)
!!$    integer, intent(in)	:: ktabj(lenj)
    integer, intent(in) :: iattr_j(:,:)
    real, intent(in) :: rattr_j(:,:)

    character*1, intent(out)	:: Mtyp
!!$    real,    intent(out)	:: corO(leni*lenj)
    real, intent(out) :: corO(size(iattr_i,2)*size(iattr_j,2))
    integer, intent(out)	:: istat
  end subroutine offdcorO
!=======================================================================
end interface
end module cordriv

max_stars_repo_name: GEOS-ESM/AeroApps

max_forks_repo_name: GEOS-ESM/AeroApps


text:
This has multiple meanings.  You are probably looking for information on one of the following:

  UC Davis Italian Program  Language program offered by UC Davis
  Restaurants Italian Restaurants  Local Italian restaurants

This is a Disambiguation disambiguation page  a navigational aid which lists other pages that might otherwise share the same title. If an article link referred you here, you might want to go back and fix it to point directly to the intended page.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
! :: ----------------------------------------------------------
! :: Volume-weight average the fine grid data onto the coarse
! :: grid.  Overlap is given in coarse grid coordinates.
! ::
! :: INPUTS / OUTPUTS:
! ::  crse      <=  coarse grid data
! ::  clo,chi    => index limits of crse array interior
! ::  ngc        => number of ghost cells in coarse array
! ::  nvar	 => number of components in arrays
! ::  fine       => fine grid data
! ::  flo,fhi    => index limits of fine array interior
! ::  ngf        => number of ghost cells in fine array
! ::  rfine      => (ignore) used in 2-D RZ calc
! ::  lo,hi      => index limits of overlap (crse grid)
! ::  lrat       => refinement ratio
! ::
! :: NOTE:
! ::  Assumes all data cell centered
! :: ----------------------------------------------------------
! ::
      subroutine avgdown(crse,c_l1,c_l2,c_h1,c_h2,nvar, &
                         cv,cv_l1,cv_l2,cv_h1,cv_h2, &
                         fine,f_l1,f_l2,f_h1,f_h2, &
                         fv,fv_l1,fv_l2,fv_h1,fv_h2,lo,hi,lrat)
      implicit none
      integer c_l1,c_l2,c_h1,c_h2
      integer cv_l1,cv_l2,cv_h1,cv_h2
      integer f_l1,f_l2,f_h1,f_h2
      integer fv_l1,fv_l2,fv_h1,fv_h2
      integer lo(2), hi(2)
      integer nvar, lrat(2)
      double precision crse(c_l1:c_h1,c_l2:c_h2,nvar)
      double precision cv(cv_l1:cv_h1,cv_l2:cv_h2)
      double precision fine(f_l1:f_h1,f_l2:f_h2,nvar)
      double precision fv(fv_l1:fv_h1,fv_l2:fv_h2)

      integer i, j, n, ic, jc, ioff, joff
      integer lenx, leny, mxlen
      integer lratx, lraty

      lratx = lrat(1)
      lraty = lrat(2)
      lenx = hi(1)-lo(1)+1
      leny = hi(2)-lo(2)+1
      mxlen = max(lenx,leny)

      if (lenx .eq. mxlen) then
         do n = 1, nvar

!           Set coarse grid to zero on overlap
            do jc = lo(2), hi(2)
               do ic = lo(1), hi(1)
                  crse(ic,jc,n) = 0.d0
               enddo
            enddo

!           Sum fine data
            do joff = 0, lraty-1
               do jc = lo(2), hi(2)
                  j = jc*lraty + joff
                  do ioff = 0, lratx-1
                     do ic = lo(1), hi(1)
                        i = ic*lratx + ioff
                        crse(ic,jc,n) = crse(ic,jc,n) + fv(i,j) * fine(i,j,n)
                     enddo
                  enddo
               enddo
            enddo

!           Divide out by volume weight
            do jc = lo(2), hi(2)
               do ic = lo(1), hi(1)
                  crse(ic,jc,n) = crse(ic,jc,n) / cv(ic,jc)
               enddo
            enddo

         enddo

      else

         do n = 1, nvar

!           Set coarse grid to zero on overlap
            do ic = lo(1), hi(1)
               do jc = lo(2), hi(2)
                  crse(ic,jc,n) = 0.d0
               enddo
            enddo

!           Sum fine data
            do ioff = 0, lratx-1
               do ic = lo(1), hi(1)
                  i = ic*lratx + ioff
                  do joff = 0, lraty-1
                     do jc = lo(2), hi(2)
                        j = jc*lraty + joff
                        crse(ic,jc,n) = crse(ic,jc,n) + fv(i,j) * fine(i,j,n)
                     enddo
                  enddo
               enddo
            enddo

!           Divide out by volume weight
            do ic = lo(1), hi(1)
               do jc = lo(2), hi(2)
                  crse(ic,jc,n) = crse(ic,jc,n) / cv(ic,jc)
               enddo
            enddo

         enddo

      end if

      end subroutine avgdown

! ::
! :: ----------------------------------------------------------
! ::

      subroutine normalize_species_fluxes(  &
                        flux1,flux1_l1,flux1_l2,flux1_h1,flux1_h2, &
                        flux2,flux2_l1,flux2_l2,flux2_h1,flux2_h2, &
                        lo,hi)

      use network, only : nspec
      use meth_params_module, only : NVAR, URHO, UFS

      implicit none

      integer          :: lo(2),hi(2)
      integer          :: flux1_l1,flux1_l2,flux1_h1,flux1_h2
      integer          :: flux2_l1,flux2_l2,flux2_h1,flux2_h2
      double precision :: flux1(flux1_l1:flux1_h1,flux1_l2:flux1_h2,NVAR)
      double precision :: flux2(flux2_l1:flux2_h1,flux2_l2:flux2_h2,NVAR)

      ! Local variables
      integer          :: i,j,n
      double precision :: sum,fac

      do j = lo(2),hi(2)
            do i = lo(1),hi(1)+1
               sum = 0.d0
               do n = UFS, UFS+nspec-1
                  sum = sum + flux1(i,j,n)
      	       end do
               if (sum .ne. 0.d0) then
                  fac = flux1(i,j,URHO) / sum
               else
                  fac = 1.d0
               end if
               do n = UFS, UFS+nspec-1
                  flux1(i,j,n) = flux1(i,j,n) * fac
      	       end do
            end do
      end do
      do j = lo(2),hi(2)+1
            do i = lo(1),hi(1)
               sum = 0.d0
               do n = UFS, UFS+nspec-1
                  sum = sum + flux2(i,j,n)
      	       end do
               if (sum .ne. 0.d0) then
                  fac = flux2(i,j,URHO) / sum
               else
                  fac = 1.d0
               end if
               do n = UFS, UFS+nspec-1
                  flux2(i,j,n) = flux2(i,j,n) * fac
      	       end do
            end do
      end do

      end subroutine normalize_species_fluxes

! :::
! ::: ------------------------------------------------------------------
! :::

      subroutine enforce_nonnegative_species(uout,uout_l1,uout_l2,uout_h1,uout_h2,lo,hi)

      use network, only : nspec
      use meth_params_module, only : NVAR, URHO, UFS

      implicit none

      integer          :: lo(2), hi(2)
      integer          :: uout_l1,uout_l2,uout_h1,uout_h2
      double precision :: uout(uout_l1:uout_h1,uout_l2:uout_h2,NVAR)

      ! Local variables
      integer          :: i,j,n,nn
      integer          :: int_dom_spec
      logical          :: any_negative
      double precision :: dom_spec,x,eps

      eps = -1.0d-16

      do j = lo(2),hi(2)
      do i = lo(1),hi(1)

         any_negative = .false.

         ! First deal with tiny undershoots by just setting them to zero
         do n = UFS, UFS+nspec-1
           if (uout(i,j,n) .lt. 0.d0) then
              x = uout(i,j,n)/uout(i,j,URHO)
              if (x .gt. eps) then
                 uout(i,j,n) = 0.d0
              else
                 any_negative = .true.
              end if
           end if
         end do

         ! We know there are one or more undershoots needing correction
         if (any_negative) then

            ! Find the dominant species
            int_dom_spec  = UFS
            dom_spec      = uout(i,j,int_dom_spec)

            do n = UFS,UFS+nspec-1
              if (uout(i,j,n) .gt. dom_spec) then
                dom_spec = uout(i,j,n)
                int_dom_spec = n
              end if
            end do

           ! Now take care of undershoots greater in magnitude than 1e-16.
           do n = UFS, UFS+nspec-1

              if (uout(i,j,n) .lt. 0.d0) then

                 x = uout(i,j,n)/uout(i,j,URHO)

                 ! Here we only print the bigger negative values
                 if (x .lt. -1.d-2) then
                    print *,'At cell (i,j) = ',i,j
                    print *,'... Fixing negative species ',n           ,' with X = ',x
                    print *,'...   from dominant species ',int_dom_spec,' with X = ',&
                             uout(i,j,int_dom_spec) / uout(i,j,URHO)
                 end if

                 ! Take enough from the dominant species to fill the negative one.
                 uout(i,j,int_dom_spec) = uout(i,j,int_dom_spec) + uout(i,j,n)

                 ! Test that we didn't make the dominant species negative
                 if (uout(i,j,int_dom_spec) .lt. 0.d0) then
                    print *,'Just made dominant species negative ',int_dom_spec,' at ',i,j
                    print *,'... We were fixing species ',n,' which had value ',x
                    print *,'... Dominant species became ',uout(i,j,int_dom_spec) / uout(i,j,URHO)
                    call bl_error("Error:: ACTReact_2d.f90 :: enforce_nonnegative_species")
                 end if

                 ! Now the negative species to zero
                 uout(i,j,n) = 0.d0

              end if

           enddo
         end if

      enddo
      enddo

      ! Test again to make sure all species are between 0 and 1.
      do j = lo(2),hi(2)
      do i = lo(1),hi(1)

         any_negative = .false.

         do n = UFS, UFS+nspec-1
           if (uout(i,j,n) .lt. eps) then
              print *,'Species ',n,' still negative at (i,j) = ',i,j,uout(i,j,n)
              any_negative = .true.
           end if
           if ( uout(i,j,n) .gt. 1.d0 .and. abs(uout(i,j,n)-1.d0) .gt. 1.e-15) then
              print *,'Species ',n,' still overshoots at (i,j) = ',i,j,(uout(i,j,n)-1.d0)
              any_negative = .true.
           end if
         end do
      enddo
      enddo

      end subroutine enforce_nonnegative_species

! :::
! ::: ------------------------------------------------------------------
! :::

      subroutine normalize_new_species(u,u_l1,u_l2,u_h1,u_h2,lo,hi)

      use network, only : nspec
      use meth_params_module, only : NVAR, URHO, UFS

      implicit none

      integer          :: lo(2), hi(2)
      integer          :: u_l1,u_l2,u_h1,u_h2
      double precision :: u(u_l1:u_h1,u_l2:u_h2,NVAR)

      ! Local variables
      integer          :: i,j,n
      double precision :: fac,sum

      do j = lo(2),hi(2)
         do i = lo(1),hi(1)
            sum = 0.d0
            do n = UFS, UFS+nspec-1
               sum = sum + u(i,j,n)
            end do
            if (sum .ne. 0.d0) then
               fac = u(i,j,URHO) / sum
            else
               fac = 1.d0
            end if
            do n = UFS, UFS+nspec-1
               u(i,j,n) = u(i,j,n) * fac
            end do
         end do
      end do

      end subroutine normalize_new_species

max_stars_repo_name: memmett/BoxLib

max_issues_repo_name: memmett/BoxLib

max_forks_repo_name: memmett/BoxLib


text:
This editor can edit this entry and tell us a bit about themselves by clicking the Edit icon.

Hey, I saw the trouble you were having, so I fixed the photo of Monkey.  Best of luck! Users/JabberWokky Evan JabberWokky Edwards