text:
C FUNCTION ISTRLN (STRING)  Returns index of last non-blank
C                           character.  Returns zero if string is
C                           null or all blank.

      FUNCTION ISTRLN (STRING)
      CHARACTER*(*)  STRING
      CHARACTER BLANK, TAB
      PARAMETER (BLANK = ' ', TAB = '	')
C     there is a tab character here  ^

C  -- If null string or blank string, return length zero.
      ISTRLN = 0
      IF (STRING (1:1) .EQ. CHAR(0))  RETURN
      IF (STRING .EQ. ' ')  RETURN

C  -- Find rightmost non-blank character.
      ILEN = LEN (STRING)
      DO 20  I = ILEN, 1, -1
         IF (STRING(I:I).NE.BLANK .AND. STRING(I:I).NE.TAB
     &        .AND. STRING(I:I).NE.CHAR(0))  GOTO 30
   20 CONTINUE
   30 ISTRLN = I

      RETURN
      END
C SUBROUTINE TRIML (STRING)  Removes leading blanks.

      SUBROUTINE TRIML (STRING)
      CHARACTER*(*)  STRING
      CHARACTER*200  TMP
      CHARACTER BLANK, TAB
      PARAMETER (BLANK = ' ', TAB = '	')
C     there is a tab character here  ^

      JLEN = ISTRLN (STRING)

C  -- All blank and null strings are special cases.
      IF (JLEN .EQ. 0)  RETURN

C  -- FInd first non-blank char
      DO 10  I = 1, JLEN
         IF (STRING(I:I).NE.BLANK .AND. STRING(I:I).NE.TAB)  GOTO 20
   10 CONTINUE
   20 CONTINUE

C  -- If I is greater than JLEN, no non-blanks were found.
      IF (I .GT. JLEN)  RETURN

C  -- Remove the leading blanks.
      TMP = STRING (I:)
      STRING = TMP
      RETURN
      END
C SUBROUTINE UPPER (STRING)  Changes a-z to upper case.

      SUBROUTINE UPPER (STRING)
      CHARACTER*(*)  STRING

      JLEN = ISTRLN (STRING)

      DO 10  I = 1, JLEN
         IC = ICHAR (STRING (I:I))
         IF ((IC .LT. 97)  .OR.  (IC .GT. 122))  GOTO 10
         STRING (I:I) = CHAR (IC - 32)
   10 CONTINUE

      RETURN
      END
C SUBROUTINE LOWER (STRING)  Changes A-Z to lower case.

      SUBROUTINE LOWER (STRING)
      CHARACTER*(*)  STRING

      JLEN = ISTRLN (STRING)

      DO 10  I = 1, JLEN
         IC = ICHAR (STRING (I:I))
         IF ((IC .LT. 65) .OR.  (IC .GT. 90))  GOTO 10
         STRING (I:I) = CHAR (IC + 32)
   10 CONTINUE

      RETURN
      END
C***********************************************************************
C
      SUBROUTINE BWORDS (S, NWORDS, WORDS)
C
C     Breaks string into words.  Words are seperated by one or more
C     blanks or tabs, or a comma and zero or more blanks.
C
C     ARGS        I/O      DESCRIPTION
C     ----        ---      -----------
C     S            I       CHAR*(*)  String to be broken up
C     NWORDS      I/O      Input:  Maximum number of words to get
C                          Output: Number of words found
C     WORDS(NWORDS) O      CHAR*(*) WORDS(NWORDS)
C                          Contains words found.  WORDS(J), where J is
C                          greater then NWORDS found, are undefined on
C                          output.
C
C      Written by:  Steven Zabinsky, September 1984
C      Tab char added July 1994.
C
C**************************  Deo Soli Gloria  **************************

C  -- No floating point numbers in this routine.
      IMPLICIT INTEGER (A-Z)

      CHARACTER*(*) S, WORDS(NWORDS)

      CHARACTER BLANK, COMMA, TAB
      PARAMETER (BLANK = ' ', COMMA = ',', TAB = '	')
C     there is a tab character here               ^.

C  -- BETW    .TRUE. if between words
C     COMFND  .TRUE. if between words and a comma has already been found
      LOGICAL BETW, COMFND

C  -- Maximum number of words allowed
      WORDSX = NWORDS

C  -- SLEN is last non-blank character in string
      SLEN = ISTRLN (S)

C  -- All blank string is special case
      IF (SLEN .EQ. 0)  THEN
         NWORDS = 0
         RETURN
      ENDIF

C  -- BEGC is beginning character of a word
      BEGC = 1
      NWORDS = 0

      BETW   = .TRUE.
      COMFND = .TRUE.

      DO 10  I = 1, SLEN
         IF (S(I:I) .EQ. BLANK .OR. S(I:I) .EQ. TAB)  THEN
            IF (.NOT. BETW)  THEN
               NWORDS = NWORDS + 1
               WORDS (NWORDS) = S (BEGC : I-1)
               BETW = .TRUE.
               COMFND = .FALSE.
            ENDIF
         ELSEIF (S(I:I) .EQ. COMMA)  THEN
            IF (.NOT. BETW)  THEN
               NWORDS = NWORDS + 1
               WORDS (NWORDS) = S(BEGC : I-1)
               BETW = .TRUE.
            ELSEIF (COMFND)  THEN
               NWORDS = NWORDS + 1
               WORDS (NWORDS) = BLANK
            ENDIF
            COMFND = .TRUE.
         ELSE
            IF (BETW)  THEN
               BETW = .FALSE.
               BEGC = I
            ENDIF
         ENDIF

         IF (NWORDS .GE. WORDSX)  RETURN

   10 CONTINUE

      IF (.NOT. BETW  .AND.  NWORDS .LT. WORDSX)  THEN
         NWORDS = NWORDS + 1
         WORDS (NWORDS) = S (BEGC :SLEN)
      ENDIF

      RETURN
      END
      SUBROUTINE UNTAB (STRING)
C REPLACE TABS WITH BLANKS :    TAB IS ASCII DEPENDENT
      INTEGER        ITAB , I, ILEN, ISTRLN
      PARAMETER      (ITAB = 9)
      CHARACTER*(*)  STRING, TAB*1
      EXTERNAL ISTRLN
      TAB  = CHAR(ITAB)
      ILEN = MAX(1, ISTRLN(STRING))
 10   CONTINUE
        I = INDEX(STRING(:ILEN), TAB )
        IF (I .NE. 0) THEN
            STRING(I:I) = ' '
            GO TO 10
        END IF
      RETURN
C END SUBROUTINE UNTAB
      END

      logical function iscomm (line)
c     returns true if line is a comment or blank line, false otherwise
c#mn{ rewritten to allow ";*%#" as comment characters
       character*(*) line
       iscomm = ((line.eq.' ').or.(index(';*%#',line(1:1)).ne.0))
c#mn}
      return
      end


text:
************************************************************************
*
*     Subroutine OUTEQN2                              Called by: CALIBR
*
*     print out concentration regression equation
*
************************************************************************
      SUBROUTINE OUTEQN2(BESTMOD,CCONC,CCONCAML,CUNITSTR,DVNAME,LDOUT,
     &                   NEXPL,NOBSC,NPAR,PARAMLC,PVALC,RSQC,SSECONC,
     &                   SSQCOBS,STDDEV,SUMCAML,SUMCOBS,TRANS)
*
*     dimensional parameters
*
      INCLUDE 'fmodules.inc'
*
*     subroutine arguments
*
      INTEGER*4 BESTMOD,LDOUT,NEXPL,NOBSC,NPAR
      DOUBLE PRECISION RSQC,SSECONC,SSQCOBS,SUMCAML,SUMCOBS
      DOUBLE PRECISION CCONC(*),CCONCAML(*),PARAMLC(*),PVALC(*)
      DOUBLE PRECISION STDDEV(MAXPARMS+1,*)
      CHARACTER*(*) CUNITSTR
      CHARACTER*5 DVNAME(*),TRANS(*)
*
*     local vars
*
      INTEGER*4 I
      CHARACTER*2 COEFSTR(MAXPARMS) /'a0','a1','a2','a3','a4','a5','a6'/
*
*     format statements
*
 1000 FORMAT(///,1X,70('-'),//,1X,'Constituent Output File Part Ib: ',
     &       'Calibration (Concentration Regression)',//,1X,70('-'),/)
 1200 FORMAT(//,' AMLE Regression Statistics',/,1X,26('-'))
 1400 FORMAT(/,' Model #',I2,' was selected for the load regression ',
     &       '(PART Ia) and is used here:',/)
 2001 FORMAT(' Ln(Conc) = a0 + a1 LnQ')
 2002 FORMAT(' Ln(Conc) = a0 + a1 LnQ + a2 LnQ^2')
 2003 FORMAT(' Ln(Conc) = a0 + a1 LnQ + a2 dtime')
 2004 FORMAT(' Ln(Conc) = a0 + a1 LnQ + a2 Sin(2 pi dtime) + a3 ',
     &       'Cos(2 pi dtime)')
 2005 FORMAT(' Ln(Conc) = a0 + a1 LnQ + a2 LnQ^2 + a3 ',
     &       'dtime')
 2006 FORMAT(' Ln(Conc) = a0 + a1 LnQ + a2 LnQ^2 + a3 ',
     &       'Sin(2 pi dtime) + a4 Cos(2 pi dtime)')
 2007 FORMAT(' Ln(Conc) = a0 + a1 LnQ + a2 Sin( pi dtime) + a3 ',
     &       'Cos(2 pi dtime) + a4 dtime')
 2008 FORMAT(' Ln(Conc) = a0 + a1 LnQ + a2 LnQ^2 + a3 ',
     &       'Sin(2 pi dtime) + a4 Cos(2 pi dtime)',//,12X,'+ a5 dtime')
 2009 FORMAT(' Ln(Conc) = a0 + a1 LnQ + a2 LnQ^2 + a3 Sin(2 pi dtime)',
     &       ' + a4 Cos(2 pi dtime)',//,12X,'+ a5 dtime ',
     &       '+ a6 dtime^2')
 2010 FORMAT(' Ln(Conc) = a0 + a1 per + a2 LnQ + a3 ',
     &       'LnQ*per')
 2011 FORMAT(' Ln(Conc) = a0 + a1 per + a2 LnQ + a3 ',
     &       'LnQ*per + a4 LnQ^2',//,12X,'+ a5 ',
     &       'LnQ^2*per')
 2012 FORMAT(' Ln(Conc) =   a0')
 2020 FORMAT(12X,'+ a',I1,1X,A5)
 2021 FORMAT(12X,'+ a',I1,1X,A5,'^2')
 2022 FORMAT(12X,'+ a',I1,1X,'sqrt(',A5,')')
 2023 FORMAT(12X,'+ a',I1,1X,'LnQ')
 2024 FORMAT(12X,'+ a',I1,1X,'ln(',A5,')')
 2025 FORMAT(12X,'+ a',I1,1X,'LnQ^2')
 2026 FORMAT(12X,'+ a',I1,1X,'ln(',A5,')^2')
 2027 FORMAT(12X,'+ a',I1,1X,'sin(2*Pi*',A5,')')
 2028 FORMAT(12X,'+ a',I1,1X,'sin(4*Pi*',A5,')')
 2029 FORMAT(12X,'+ a',I1,1X,'sin(6*Pi*',A5,')')
 2030 FORMAT(12X,'+ a',I1,1X,'cos(2*Pi*',A5,')')
 2031 FORMAT(12X,'+ a',I1,1X,'cos(4*Pi*',A5,')')
 2032 FORMAT(12X,'+ a',I1,1X,'cos(6*Pi*',A5,')')
 2100 FORMAT(/' where:',/,7X,'Conc  = constituent concentration',/,7X,
     &       'LnQ   = Ln(Q) - center of Ln(Q)')
 2110 FORMAT(7X,'dtime = decimal time - center of decimal time')
 2120 FORMAT(7X,'per   = period; 1 or 0 depending on defined period')
 3000 FORMAT(//,' Concentration Regression Results',/,1X,32('-'),/,
     &       ' R-Squared [%]                  : ',F5.2,/,
     &       ' Residual Variance              : ',F6.4,//,
     &       ' Coeff.',4X,'Value',9X,'Std.Dev.',5X,'t-ratio',5X,
     &       'P Value',/,1X,58('-'))
 3100 FORMAT(1X,A2,7X,0PF7.4,8X,0PF6.4,5X,0PF7.2,7X,1PE9.3)
 3200 FORMAT(1X,A2,7X,0PF7.4,8X,0PF6.4,5X,0PF7.2,7X,'See WARNING')
 4600 FORMAT(//,' Comparison of Observed and Estimated Concentrations',
     &       /,1X,51('-'),
     &       /,'   The summary statistics and bias diagnostics ',
     &         'presented below are',
     &       /,' based on a comparison of observed and estimated ',
     &         'concentrations for',
     &       /,' all dates/times within the calibration data set. ',
     &         'Although this'
     &       /,' comparison does not directly address errors in ',
     &         'concentration estimation',
     &       /,' for unsampled dates/times, large discrepancies',
     &         ' between observed and',
     &       /,' estimated concentrations are indicative of a poor ',
     &         'model fit. Additional',
     &       /,' details and warnings are ',
     &         'provided below.',
     &       //,' Note: The comparison that follows uses a ',
     &         'concentration equal to 1/2 the',
     &       /,' detection limit when an observation ',
     &         'is censored. The summary stats and',
     &       /,' bias diagnostics are therefore slightly inaccurate ',
     &         'for censored data sets.')
 4700 FORMAT(//,' Summary Stats: Est. and Obs. Concentrations in ',A9,/,
     &       1X,59('-'),/,17X,'25th',14X,'75th',5X,'90th',5X,'95th',5X,
     &       '99th',/,8X,'Min.',5X,'Pct',6X,'Med.',5X,'Pct',
     &       3(6X,'Pct'),6X,'Max.',/,7X,72('-'))
 4800 FORMAT(//,' Bias Diagnostics',/,1X,16('-'),/,' Bp [%] ',F9.3,
     &       /,' PCR',4X,F9.3,/,' E',6X,F9.3,//,' where:',//,4X,
     &       'Bp',T9,'Concentration Bias in Percent',/,T9,'Positive '
     &       '(negative) values indicate over (under) ',
     &       'estimation.',/,T9,'***The model should not be ',
     &       'used when the + or - bias exceeds 25%***',
     &       //,4X,'PCR',T9,'Partial Concentration Ratio',/,T9,'Sum of',
     &       ' est. concentrations divided by sum of obs. ',
     &       'concentrations.',/,T9,'Values > 1 indicate over',
     &       'estimation; values < 1 indicate underestimation.',/,T9,
     &       'PCR = (Bp + 100) / 100',//,4X,'E',T9,
     &       'Nash Sutcliffe Efficiency Index',/,T9,'E ranges',
     &       ' from -infinity to 1.0',/,T9,'E = 1; a perfect',
     &       ' fit to observed data.',/,T9,'E = 0; model ',
     &       'estimates are as accurate as the mean of observed ',
     &       'data.',/,T9,'E < 0; the observed mean is a',
     &       ' better estimate than the model estimates.',/)
 4900 FORMAT(' IMPORTANT WARNING:',//,' Concentration Bias (Bp) ',
     &       'Exceeds + or - 25%',/,' THE CALIBRATED MODEL SHOULD NOT ',
     &       'BE USED FOR CONCENTRATION ESTIMATION',//)
 5000 FORMAT(' IMPORTANT WARNING:',//,' Nash Sutcliffe Efficiency ',
     &       'Index (E) is less than zero.',/,' The observed mean is ',
     &       'a better estimate than the model estimates.',/,' THE ',
     &       'CALIBRATED MODEL SHOULD NOT BE USED FOR CONCENTRATION ',
     &       'ESTIMATION',//)
 5100 FORMAT(/,' NOTE: Additional information on model calibration is ',
     &       'included in the',/,7X,'residual output file. LOADEST ',
     &       'users should conduct a thorough',/,7X,'residuals ',
     &       'analysis using the data contained therein (checks for',/,
     &       7X,'heteroscedasticity and non-normality).  Example ',
     &       'residual plots',/,7X,'are shown in Figures 7, 8, 9, and',
     &       ' 17 of the LOADEST documentation',/,7X,
     &       '(Runkel et al., 2004).')
*
*     output selected model number; If the regression model is user
*     defined (MODNO=99), BESTMOD will equal NMODELS, as a user defined
*     model is the last model - output '99' rather than BESTMOD in this
*     case.
*
      WRITE(LDOUT,1000)
      WRITE(LDOUT,1200)
      IF (BESTMOD .EQ. NMODELS) THEN
         WRITE(LDOUT,1400) 99
      ELSE
         WRITE(LDOUT,1400) BESTMOD
      ENDIF
*
*     output model equation
*
      IF (BESTMOD.EQ.1) THEN
         WRITE(LDOUT,2001)
         WRITE(LDOUT,2100)
      ELSEIF (BESTMOD.EQ.2) THEN
         WRITE(LDOUT,2002)
         WRITE(LDOUT,2100)
      ELSEIF (BESTMOD.EQ.3) THEN
         WRITE(LDOUT,2003)
         WRITE(LDOUT,2100)
         WRITE(LDOUT,2110)
      ELSEIF (BESTMOD.EQ.4) THEN
         WRITE(LDOUT,2004)
         WRITE(LDOUT,2100)
         WRITE(LDOUT,2110)
      ELSEIF (BESTMOD.EQ.5) THEN
         WRITE(LDOUT,2005)
         WRITE(LDOUT,2100)
         WRITE(LDOUT,2110)
      ELSEIF (BESTMOD.EQ.6) THEN
         WRITE(LDOUT,2006)
         WRITE(LDOUT,2100)
         WRITE(LDOUT,2110)
      ELSEIF (BESTMOD.EQ.7) THEN
         WRITE(LDOUT,2007)
         WRITE(LDOUT,2100)
         WRITE(LDOUT,2110)
      ELSEIF (BESTMOD.EQ.8) THEN
         WRITE(LDOUT,2008)
         WRITE(LDOUT,2100)
         WRITE(LDOUT,2110)
      ELSEIF (BESTMOD.EQ.9) THEN
         WRITE(LDOUT,2009)
         WRITE(LDOUT,2100)
         WRITE(LDOUT,2110)
      ELSEIF (BESTMOD.EQ.10) THEN
         WRITE(LDOUT,2010)
         WRITE(LDOUT,2100)
         WRITE(LDOUT,2120)
      ELSEIF (BESTMOD.EQ.11) THEN
         WRITE(LDOUT,2011)
         WRITE(LDOUT,2100)
         WRITE(LDOUT,2120)
      ELSEIF (BESTMOD.EQ.NMODELS) THEN
         WRITE(LDOUT,2012)
         DO 20 I=1,NEXPL
            IF (TRANS(I) .EQ. 'NONE') THEN
               WRITE(LDOUT,2020) I,DVNAME(I)
            ELSEIF(TRANS(I) .EQ. 'SQ') THEN
               WRITE(LDOUT,2021) I,DVNAME(I)
            ELSEIF(TRANS(I) .EQ. 'SQRT') THEN
               WRITE(LDOUT,2022) I,DVNAME(I)
            ELSEIF(TRANS(I) .EQ. 'LN') THEN
               IF (DVNAME(I) .EQ. 'Q') THEN
                  WRITE(LDOUT,2023) I
               ELSE
                  WRITE(LDOUT,2024) I,DVNAME(I)
               ENDIF
            ELSEIF(TRANS(I) .EQ. 'LNSQ') THEN
               IF (DVNAME(I) .EQ. 'Q') THEN
                  WRITE(LDOUT,2025) I
               ELSE
                  WRITE(LDOUT,2026) I,DVNAME(I)
               ENDIF
            ELSEIF(TRANS(I) .EQ. 'SIN2P') THEN
               WRITE(LDOUT,2027) I,DVNAME(I)
            ELSEIF(TRANS(I) .EQ. 'SIN4P') THEN
               WRITE(LDOUT,2028) I,DVNAME(I)
            ELSEIF(TRANS(I) .EQ. 'SIN6P') THEN
               WRITE(LDOUT,2029) I,DVNAME(I)
            ELSEIF(TRANS(I) .EQ. 'COS2P') THEN
               WRITE(LDOUT,2030) I,DVNAME(I)
            ELSEIF(TRANS(I) .EQ. 'COS4P') THEN
               WRITE(LDOUT,2031) I,DVNAME(I)
            ELSEIF(TRANS(I) .EQ. 'COS6P') THEN
               WRITE(LDOUT,2032) I,DVNAME(I)
            ENDIF
 20      CONTINUE
         WRITE(LDOUT,2100)
         WRITE(LDOUT,2110)
      ENDIF
*
*     output AMLE regression statistics
*
      WRITE(LDOUT,3000) RSQC,PARAMLC(NPAR+1)

      DO 30 I=1,NPAR
         IF (PVALC(I) .LT. 9999999.99) THEN
            WRITE(LDOUT,3100) COEFSTR(I),PARAMLC(I),STDDEV(I,BESTMOD),
     &                        PARAMLC(I)/STDDEV(I,BESTMOD),PVALC(I)
         ELSE
            WRITE(LDOUT,3200) COEFSTR(I),PARAMLC(I),STDDEV(I,BESTMOD),
     &                           PARAMLC(I)/STDDEV(I,BESTMOD)
         ENDIF
 30   CONTINUE
*
*     start a new section w/i Constituent Output File Part Ib; this new
*     section compares observed and estimated concentrations on the
*     days/times for which observations are available (the calibration
*     data set).
*
      WRITE(LDOUT,4600)
*
*     print summary statistics for estimated concentrations
*
      WRITE(LDOUT,4700) CUNITSTR
      CALL OUTSUMMC(NOBSC,CCONCAML,CCONC,LDOUT)
*
*     print bias diagnostics, including:
*
*     Bp   Concentration Bias in % (email from B Hirsch, 16Jan2013)
*     PCR  Partial Concentration Ratio (Stenback et al., 2011)
*     E    Nash-Sutcliffe Efficiency Index (Nash & Sutcliffe, 1970)
*
*     Issue WARNINGS if |Bp| > 25% or  E < 0
*
      WRITE(LDOUT,4800) 100.D0*(SUMCAML/SUMCOBS - 1.D0),
     &                  SUMCAML/SUMCOBS,
     &                  1.D0-SSECONC/(SSQCOBS - SUMCOBS*SUMCOBS/NOBSC)

      IF (ABS(100.D0*(SUMCAML/SUMCOBS-1.D0)).GT.25.D0) THEN
         CALL JROGER(LDOUT)
         WRITE(LDOUT,4900)
      ENDIF

      IF (1.D0-SSECONC/(SSQCOBS - SUMCOBS*SUMCOBS/NOBSC).LT.0.D0) THEN
         CALL JROGER(LDOUT)
         WRITE(LDOUT,5000)
      ENDIF
*
*     write final note in regard to residual analyses
*
      WRITE(LDOUT,5100)

      RETURN
      END

ext: f

max_stars_repo_name: lthiamodelers/baseflow-coefficients

max_issues_repo_name: lthiamodelers/baseflow-coefficients

max_forks_repo_name: lthiamodelers/baseflow-coefficients


text:
      SUBROUTINE ZERO
*
*
*       Initialization of global scalars.
*       ---------------------------------
*
      INCLUDE 'common6.h'
*
*
*       Initialize parameters & counters and set useful constants.
      TIME = 0.0D0
      TOFF = 0.0D0
      TTOT = 0.0D0
      TADJ = 0.0D0
      TBLIST = 0.0D0
      TPLOT = 0.0D0
      DTB = 0.0D0
      TNEXT = 0.0D0
      TMDOT = 0.0D0
      ZMDOT = 0.0D0
      TLASTS = 0.0D0
      TLASTT = 0.0D0
      RMAX = 0.0
      EPOCH0 = 0.0
      DMIN1 = 100.0
      DMIN2 = 100.0
      DMIN3 = 100.0
      DMIN4 = 100.0
      DMINC = 100.0
      CPUTOT = 0.0
      WTOT = 0.0D0
      ERROR = 0.0D0
      ERRTOT = 0.0D0
      DETOT = 0.0D0
      ETCORR = 0.0
      EBIN0 = 0.0D0
      ESUB = 0.0D0
      EMERGE = 0.0D0
      EDISS = 0.0D0
      ECOLL = 0.0D0
      ECDOT = 0.0D0
      EMDOT = 0.0D0
      EKICK = 0.0
      ESYNC = 0.0D0
      ETIDE = 0.0D0
      ESESC = 0.0D0
      EBESC = 0.0D0
      EMESC = 0.0D0
      EGRAV = 0.0
      SBCOLL = 0.0D0
      BBCOLL = 0.0D0
      CHCOLL = 0.0D0
      NPAIRS = 0
      NBPREV = 0
      NNTB = 0
      NMERGE = 0
      NEWHI = 0
      NSUB = 0
      NCH = 0
      NBIN0 = 0
      NHI0 = 0
      NDISR = 0
      IFIRST = 1
      IPHASE = 0
      ICOMP = 0
      JCOMP = 0
      KSPAIR = 0
      NPRINT = 0
      MODEL = 0
      NDUMP = 0
      NSTEPS = 0
      NTIMER = 0
      LISTR(1) = 0
      LISTD(1) = 0
      LISTV(1) = 0
      NTAIL = 0
      NTTOT = 0
      NRSAVE = 0
      NNBOPT = 0  ! NBODY6++ compatibility.
*       Jpetts - initialise DYNFRI vars
      COULOG = 0.0D0
      BDENS = 0.0D0
      BVDISP = 0.0D0
      BRATIO = 0.0D0
      BDIST = 0.0D0
      EDYNFRI = 0.0D0
      DEDYNFRI = 0.0D0
      RGDENSMAG = 0.0D0
      VGCOREMAG = 0.0D0
      MASSCL = 1.0D0
      NBOUND = 0
      DYNFTPREV = 0.0D0
      DFOFF = 0
*
      DO 10 J = 1,KMAX
          KVEC(2*J-1) = J
          KVEC(2*J) = J
   10 CONTINUE
*
      DO 20 J = 1,12
          E(J) = 0.0D0
   20 CONTINUE
*
      DO 30 J = 1,9
          TLASTB(J) = 0.0D0
   30 CONTINUE
*
      DO 40 J = 1,NMAX
          RADIUS(J) = 0.0D0
          BODY0(J) = 0.0
          EPOCH(J) = 0.0
          TEV(J) = 0.0
          TEV0(J) = 0.0
          KSTAR(J) = 0
*       Jpetts - initialise DYNFRI vars
          DYNFPOT(J) = 0.0D0
          DYNFPOTBG(J) = 0.0D0
          BOUND(J) = 0.0D0
   40 CONTINUE
*
      DO 50 K = 1,3
          RDENS(K) = 0.0D0
          TIDAL(K) = 0.0
*       Jpetts - initialise DYNFRI vars
          RGDENS(K) = 0.0D0
          VCORE(K) = 0.0D0
          VGCORE(K) = 0.0D0
          VPREV(K) = 0.0D0
          VDOT(K) = 0.0D0
   50 CONTINUE
      TIDAL(4) = 0.0
*
*       Set fractional constants & two PI.
      ONE3 = 1.0D0/3.0D0
      ONE6 = 1.0D0/6.0D0
      ONE9 = 1.0D0/9.0D0
      ONE12 = 1.0D0/12.0D0
      TWOPI = 8.0D0*ATAN(1.0D0)
*       Jpetts - define 1xpi
      ONEPI = 0.5d0*TWOPI
      TINY = 1.0D-14
*
      RETURN
*
      END


text:
C The contents of this file are copyright protected:
C --------------------------------------------------
C NCAR Graphics - UNIX Version 3.1.2
C Copyright (C) 1987, 1988, 1989, 1991
C University Corporation for Atmospheric Research
C All Rights Reserved
C
C Modified for high resolution plots as part
C of the shallow water model STSWM by R. Jakob
C
C-----------------------------------------------------------------------
C
      SUBROUTINE CPCNRC (ZDAT,KZDT,MZDT,NZDT,FINC,NDSH)
C
      INTEGER KZDT,MZDT,NZDT,NDSH
      REAL ZDAT, FINC
      DIMENSION ZDAT(KZDT,*)
C
C This routine simulates the old routine CONREC.
C
C Define some needed dimensions.
C
      INTEGER LRWK, LIWK, LAWK
      PARAMETER (LRWK=4000,LIWK=40000,LAWK=300000)
      INTEGER LIRA, LCRA
      PARAMETER (LIRA=1000,LCRA=30000)
C
C     Flag to turn on labels
C
      LOGICAL PLABEL
      PARAMETER (PLABEL = .FALSE.)
C
C Declare arrays needed by ARSCAM
C
      INTEGER IARA,IGRA
      REAL XCRA,YCRA
      DIMENSION IARA(LIRA),IGRA(LIRA)
      DIMENSION XCRA(LCRA),YCRA(LCRA)
C
C Define required workspace arrays.
C
      REAL RWRK
      INTEGER IWRK,IAMA
      DIMENSION RWRK(LRWK),IWRK(LIWK),IAMA(LAWK)
C
C Contour level value
C
      REAL CLEV
C
C Temporaries
C
      INTEGER IDSH,ICFF,IWU,IRWU,NCLO,ICLO,ICLU
C
C     ROUTINE FOR DRAWING CONTOUR LINES
C
      EXTERNAL DRAWCL
C
C     ROUTINE FOR SHADING (NEGATIVE VALUES)
C
      EXTERNAL SHADER
C
C     BIT OPERATIONS
C
      INTEGER IAND, IOR, ISHIFT
      EXTERNAL IAND, IOR, ISHIFT
C
C GKS aspect source flags
C
      INTEGER IASF
      DIMENSION IASF(13)
      DATA IASF / 13*1 /
C
C-----------------------------------------------------------------------
C
C Turn off clipping, so that the informational label won't disappear.
C
      CALL GSCLIP (0)
C
C     Set GKS aspect source flags to 'individual'
C
      CALL GSASF (IASF)
C
C Arrange for the selection of contour levels as desired by the user.
C
      IF (FINC.LT.0.) THEN
        CALL CPSETI ('CLS - CONTOUR LEVEL SELECTOR',MAX(1,INT(-FINC)))
        CALL CPSETR ('CIS - CONTOUR INTERVAL SPECIFIER',0.)
      ELSE IF (FINC.EQ.0.) THEN
        CALL CPSETI ('CLS - CONTOUR LEVEL SELECTOR',16)
        CALL CPSETR ('CIS - CONTOUR INTERVAL SPECIFIER',0.)
      ELSE
        CALL CPSETI ('CLS - CONTOUR LEVEL SELECTOR',1)
        CALL CPSETR ('CIS - CONTOUR INTERVAL SPECIFIER',FINC)
      END IF
C
C     Turn on positioning of labels by penalty scheme
C
      CALL CPSETI ('LLP - LINE LABEL POSITIONING',2)
C
C Decide what dash pattern to use.
C
      IDSH=ABS(NDSH)
      IF ((IDSH .EQ. 0) .OR. (IDSH .EQ. 1) .OR.
     &    (IDSH .EQ. 1023)) THEN
        IDSH=IOR(ISHIFT(32767,1),1)
      ELSE
        IDSH=IOR(ISHIFT(IDSH,6),IAND(ISHIFT(IDSH,-4),63))
      END IF
C
      IF (PLABEL) THEN
C
C Label highs and lows.
C
         CALL CPSETC ('HLT - HIGH/LOW LABEL TEXT',
     &               'H:B:$ZDV$:E:''L:B:$ZDV$:E:')
         CALL CPSETC ('HIT - HIGH/LOW LABEL TEXT',
     &               'H:B:$ZDV$:E:')
         CALL CPSETC ('LOT - HIGH/LOW LABEL TEXT',
     &               'L:B:$ZDV$:E:')
C
C And put box around it
C
         CALL CPSETI('HLB - HIGH/LOW LABEL BOX FLAG',1)
         CALL CPSETI('HLO - HIGH/LOW LABEL OVERLAP FLAG',11)
      ELSE
C
C No high/low labels
C
         CALL CPSETC ('HLT - HIGH/LOW LABEL TEXT','''')
      ENDIF
C
C Initialize CONPACK and give it all array dimensions.
C
      CALL CPRECT (ZDAT,KZDT,MZDT,NZDT,RWRK,LRWK,IWRK,LIWK)
C
C If the field was constant, skip some of the following code.
C
      CALL CPGETI ('CFF - CONSTANT FIELD FLAG',ICFF)
C
      IF (ICFF.EQ.0) THEN
C
C Pick contour levels.
C
        CALL CPPKCL (ZDAT,RWRK,IWRK)
C
C     Check work space used
C
      CALL CPGETI ('IWU',IWU)
      CALL CPGETI ('RWU',IRWU)
      IF ((IWU .GT. LIWK) .OR. (IRWU .GT. LRWK)) THEN
        WRITE(0,*)'CPPKCL WORKSPACE USED/AVAILABLE'
        WRITE(0,*)' INTEGER =',IWU,LIWK
        WRITE(0,*)' REAL =', IRWU,LRWK
      ENDIF
C
C Use fat lines each fifth contour interval, mark zero contour
C for later shading, and set dash pattern for negative contours.
C
        CALL CPGETI ('NCL - NUMBER OF CONTOUR LEVELS',NCLO)
C
        DO 10001 ICLO=1,NCLO
          CALL CPSETI ('PAI - PARAMETER ARRAY INDEX',ICLO)
          CALL CPGETI ('CLU - CONTOUR LEVEL USE FLAG',ICLU)
          IF (ICLU .EQ. 3) THEN
             CALL CPSETI ('CLL - CONTOUR LINE WIDTH',2)
          ENDIF
          CALL CPGETR ('CLV - CONTOUR LEVEL VALUE',CLEV)
          IF (CLEV .NE. 0.0) THEN
             CALL CPSETI ('AIA - AREA IDENTIFIER ABOVE LINE',0)
             CALL CPSETI ('AIB - AREA IDENTIFIER BELOW LINE',0)
          ELSE
             CALL CPSETI ('AIA - AREA IDENTIFIER ABOVE LINE',2)
             CALL CPSETI ('AIB - AREA IDENTIFIER BELOW LINE',1)
          ENDIF
          IF (NDSH .GT. 0 .OR.
     &        ((NDSH .LT. 0.0) .AND. (CLEV .LT. 0.0))) THEN
             CALL CPSETI ('CLD - CONTOUR LINE DASH PATTERN',IDSH)
          ENDIF
10001   CONTINUE
C
      END IF
C
C Draw the contour lines, masking them if necessary.
C
C     CALL CPCLDR (ZDAT,RWRK,IWRK)
C
C
C     INITIALIZE AREA MAP
C
      CALL ARINAM(IAMA,LAWK)
C
C     PUT LABEL BOXES INTO AREA MAP
C
      CALL CPLBAM(ZDAT,RWRK,IWRK,IAMA)
C
C     Check work space used
C
      CALL CPGETI ('IWU',IWU)
      CALL CPGETI ('RWU',IRWU)
      IF ((IWU .GT. LIWK) .OR. (IRWU .GT. LRWK)) THEN
        WRITE(0,*)'CPLBAM WORKSPACE USED/AVAILABLE'
        WRITE(0,*)' INTEGER =',IWU,LIWK
        WRITE(0,*)' REAL =', IRWU,LRWK
      ENDIF
C
C     DRAW CONTOUR LINES
C
      CALL CPCLDM(ZDAT,RWRK,IWRK,IAMA,DRAWCL)
C
C     Check work space used
C
      CALL CPGETI ('IWU',IWU)
      CALL CPGETI ('RWU',IRWU)
      IF ((IWU .GT. LIWK) .OR. (IRWU .GT. LRWK)) THEN
        WRITE(0,*)'CPCLDM WORKSPACE USED/AVAILABLE'
        WRITE(0,*)' INTEGER =',IWU,LIWK
        WRITE(0,*)' REAL =', IRWU,LRWK
      ENDIF
C
C Plot labels.
C
      CALL CPLBDR (ZDAT,RWRK,IWRK)
C
C     Check work space used
C
      CALL CPGETI ('IWU',IWU)
      CALL CPGETI ('RWU',IRWU)
      IF ((IWU .GT. LIWK) .OR. (IRWU .GT. LRWK)) THEN
        WRITE(0,*)'CPLBDR WORKSPACE USED/AVAILABLE'
        WRITE(0,*)' INTEGER =',IWU,LIWK
        WRITE(0,*)' REAL =', IRWU,LRWK
      ENDIF
C
C     Add contour lines to area map
C
      CALL CPCLAM(ZDAT,RWRK,IWRK,IAMA)
C
C     Check work space used
C
      CALL CPGETI ('IWU',IWU)
      CALL CPGETI ('RWU',IRWU)
      IF ((IWU .GT. LIWK) .OR. (IRWU .GT. LRWK)) THEN
        WRITE(0,*)'CPLBDR WORKSPACE USED/AVAILABLE'
        WRITE(0,*)' INTEGER =',IWU,LIWK
        WRITE(0,*)' REAL =', IRWU,LRWK
      ENDIF
C
C     SCAN AREA MAP FOR SHADING OF NEGATIVE VALUES
C
      CALL ARSCAM(IAMA,XCRA,YCRA,LCRA,IARA,IGRA,LIRA,SHADER)
C
      RETURN
C
      END


text: A department in the College of Letters and Science. Offers programs in UC Davis Chinese Program Chinese and UC Davis Japanese Program Japanese.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      SUBROUTINE PLTMEQ                                                 ABSH0045
C***********************************************************************ABSH0046
C*****PLTMEQ GENERATES A DATA FILE FOR THE EQUILIBRIUM QUANTITIES.     *ABSH0047
C***********************************************************************ABSH0048
      PARAMETER MNMX=61                                                 ABSH0049
      PARAMETER MTHETA=11,MM1=10,MNMOMS=3                               ABSH0050
      PARAMETER MNEQ=6                                                  ABSH0051
      PARAMETER MM1T=20                                                 ABSH0052
      COMMON /CD02/ XP(MM1), C(MM1,MNEQ), PARAM(MNEQ)                   ABSH0053
      COMMON /LOCGPA/ BETAT0, F0, FWALL, PSIC0, Q0, ARGEOM,             ABSH0054
     *     BETAPB, BETATB                                               ABSH0055
      COMMON /LOCMEQ/ XL0, DX, R0, A0, BT0, P0, AI0, E0, E1,            ABSH0056
     *     D0, ARC, SLPE, SLPR0, BETAP0, KMHD, IUNDER                   ABSH0057
      COMMON /NIO/ NIN, NOUT, NPLOT, NDEBUG, NTTY                       ABSH0058
      COMMON /SWORK1/ PRES(MNMX), CURT(MNMX), RHO(MNMX),                ABSH0059
     *     BPRES(MNMX), CPRES(MNMX), DPRES(MNMX), BCURT(MNMX),          ABSH0060
     *     CCURT(MNMX), DCURT(MNMX), NRHO                               ABSH0061
      COMMON /SWORK3/ PH(MM1), PPH(MM1), AIH(MM1), AIPH(MM1),           ABSH0062
     *     PSI(MM1), FPOL(MM1), SMQ(MM1), AGSQRT(MM1),                  ABSH0063
     *     AGTTGS(MM1), ATAUOR(MM1), AGRHO(MM1), ARIF(MM1),             ABSH0064
     *     AR2IF(MM1), FFP(MM1), AJORF(MM1), AJTM(MM1T),                ABSH0065
     *     RMID(MM1T), GTTGSQ(MTHETA), TAUOR(MTHETA),                   ABSH0066
     *     GRXX(MTHETA), R2I(MTHETA)                                    ABSH0067
      COMMON /VOUPUT/ SHIFT(MNMX), ELONG(MNMX), TRIANG(MNMX)            ABSH0068
      M1 = MM1                                                          ABSH0069
      WRITE (NPLOT,99999) M1, NRHO                                      ABSH0070
      WRITE (NPLOT,99998) ARGEOM, A0, D0, PSIC0, PARAM                  ABSH0071
      DO 10 I=1,MM1                                                     ABSH0072
           WRITE (NPLOT,99998) XP(I), C(I,1), C(I,3), C(I,5)            ABSH0073
           WRITE (NPLOT,99998) PSI(I), AGSQRT(I), SMQ(I),               ABSH0074
     *          AJORF(I), FFP(I), AJTM(I), AJTM(MM1+I),                 ABSH0075
     *          FPOL(I), RMID(I), RMID(MM1+I), AIH(I)                   ABSH0076
           WRITE (NPLOT,99998) PH(I), PPH(I)                            ABSH0077
   10 CONTINUE                                                          ABSH0078
      DO 20 I=1,NRHO                                                    ABSH0079
           WRITE (NPLOT,99998) RHO(I), SHIFT(I), ELONG(I),              ABSH0080
     *          TRIANG(I)                                               ABSH0081
   20 CONTINUE                                                          ABSH0082
      RETURN                                                            ABSH0083
99999 FORMAT (1X, 2I5)                                                  ABSH0084
99998 FORMAT (1X, 10E12.4)                                              ABSH0085
      END                                                               ABSH0086

max_stars_repo_name: jonathanschilling/VMOMS

max_issues_repo_name: jonathanschilling/VMOMS

max_forks_repo_name: jonathanschilling/VMOMS


text:
      subroutine runCY_00lllli(k,l,i1,Xtwiddle,Gtwiddle,Shat7,N0)
      implicit none
C---  Expression for extension of Eq. 5.60b
C---  Calculates C00llli, requires C00lllll
C---  Small terms of order Xtwiddle(0,k)*Ciiiiii,Xtwiddle(0,0)*Ciiiiiii
C---  Denominator Gtwiddle(k,l)
      include 'pvCnames.f'
      include 'pvCv.f'
      include 'Carraydef.f'
      include 'Carrays.f'
      integer ep,N0,k,l,i1,np
      parameter(np=2)
      double precision Xtwiddle(0:np,0:np),Gtwiddle(np,np)
      double complex Shat7(np,z6max,-2:0)

      if ((i1 .eq. l) .or. (i1 .eq. 0)) then
      return
      endif

      do ep=-2,0

      Cv(czziiiii(z5(l,l,l,l,i1))+N0,ep)=
     . (-2d0*Gtwiddle(k,i1)*Cv(czziiiii(z5(l,l,l,l,l))+N0,ep)
     . +Gtwiddle(k,1)*Shat7(1,z6(l,l,l,l,l,i1),ep)
     . +Gtwiddle(k,2)*Shat7(2,z6(l,l,l,l,l,i1),ep)
     . +Xtwiddle(k,0)*Cv(ciiiiii(z6(l,l,l,l,l,i1))+N0,ep)
     . -Xtwiddle(0,0)*Cv(ciiiiiii(z7(k,l,l,l,l,l,i1))+N0,ep))
     . /(10d0*Gtwiddle(k,l))

      enddo

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
The Deadly Davis Dash is a series of  Bicycling Advocacy bicycle races that has included intoxication and may include jungle gyms, singing, music, costumes ... oh yeah ... and Bicycling bicycling. Riders will not be subject to blood tests following the race (you should enter, wiki:WikiPedia:Floyd_Landis Floyd). Past sponsors have included: Bike Church, Crossbill Records, Delta of Venus, and Primary Concepts.

DDD5

3d5, the Day of the Deadly Derby, is planned for October 31st, 4pm start, 3:30 register at NW corner of 4th and L streets.

Moondash
A stripped down fast version of the Dash, this is all about speed and gumption. Inspired by the events of 3D4, this is a criteriumesque race by moonlight.

The Ministry presents...
5/8/09 11:30pm
Northwest corner of 4th and L
Bring helmet, lights, and cell phone.

DDD4
DAVIS DEADLY DERBY (3D4) happened Friday, April 10th, 2009. 5:30pm12am @ the Delta of Venus DOV

The Somali Coast Guard and Bike Church DBC presents:
The Deadly Davis Derby.  A team bicycle scavenge race. 36 person teams. Costumes highly encouraged. Team Registration at 5:30pm @ Delta of Venus DOV, race 68pm $5 All skill levels welcome. Fun>Speed.

and Bike Church Davis Bicycle Collective Benefit Dance Party with Dj’s TAO, Nix, Dog Tonez, and Jordie La Forge and a silent auction or raffle for the Bike Church DBC. @ Delta of Venus DOV (8pm12pm)
$35 sliding scale. Free admission for Derby racers. All Ages.

and the Deadliest Davis Dash, 3D4 “The Alley Trak” 12 AM MIDNIGHT Delta of Venus DOV
(Not for the amateur, squeamish, or easily injured) FREE

DDD3
The third Deadly Davis Dash was a roaming music festival with hella participants.

DDD2
The second Deadly Davis Dash, subtitled Race and Pillage, was held on September 16th, 2006.  More than 40 racers made their way through six checkpoints to two frosty kegs at the Delta of Venus.  The route was given in pieces to the riders as they arrived at the various checkpoints, which included the tops of two campus parking garages and some really outoftheway neighborhood parks.  Those who raced in costume had one less checkpoint, and were eligible for the Best Costume award.\

DDD
The first Deadly Davis Dash was held on July 8th, 2006.  Approximately 45 racers left Central Park at about 7 pm to seek out the five checkpoints around town.  At each checkpoint the riders had to complete various challenges such as climbing through a jungle gym, drinking a shot of whiskey, or singing a rendition of Im a little tea pot.  The race ended at Delta of Venus with a keg of Pabst waiting for the tired and thirsty riders.

See Bicycling Activities and Clubs for other bikerelated events.

20060717 11:40:29 nbsp who won? Users/BenTaylor

20061211 11:28:20 nbsp Is anyone else ready for #3? Users/RobiPochapin

20070316 19:59:52 nbsp Ill plan the Deadliesd Davis Dash when I get back to davis in May. wanna help? Users/PxlAted

20070316 20:59:44 nbsp How does this not get people BUIs? Is the amount of booze consumed in the course insignificant? Users/WilliamLewis
     Drinking is afterward.

20070617 04:29:20 nbsp im sorry for getting sick :( Users/KatieDelwiche

20070617 11:51:42 nbsp I went to DDD3 yesterday  it was a lot of fun. Good exercise too. Users/JoseBleckman

20070621 14:57:30 nbsp the not so deadly dash was fucking epic. Users/PxlAted

20080206 21:38:09 nbsp DDD4. March. Yes? Users/Darach

20090129 12:08:45 nbsp Yes, March 7th. Users/Darach

Make that April 10th. Darach

20090324 21:20:12 nbsp What is the Davis Bike Collective? What is the Somali Coast Guard? Are these real organizations? Is the Davis Bike Club somehow involved? This is confusing. Users/NickSchmalenberger

20090324 21:43:35 nbsp The Davis Bike Church is changing its name, still DBC, just now that stands for Davis Bike Collective. The Somali coastguard are as real as the Somali Pirates and more real than the Somali fisheries which are gone for years now or the Somali government! Another kiss from the global North to the global South. The Bike race event is real, even if one of the sponsors is never going to hear of this event. Users/PxlAted

20090328 14:33:06 nbsp The Davis Bicycle Club is running the bike jousting, dude. I hear theyve got some Cervelo tall bikes. Users/Darach

20090401 20:37:19 nbsp The third Deadly Davis Dash was a roaming music festival with hella participants.

This made me giggle. Users/BrookeB

20090402 13:38:27 nbsp THE PLAN flyer  3rd image is DDD3, where is the Plan for DDD4? The times on the flyer are going to confuse everyone who doesnt realize thats for #3 Users/GarrettGallegos


text:
      SUBROUTINE UPSURF (COORD)
C This routine is called alternatively to CONSTS during geometry
C optimization. It updates the segments on the SCAS , i.e. their
C areas and centers, according to a new geometry, but it doesnot create
C a new segmentation. The reason for using UPSURF is to smooth the total
C energy with respect to fluctations arisig from a fluctuating segmentation.
C Be aware that the use of UPSURF brings some 'hysteresis' into the
C Hamiltonian.
C UPSURF fills the A-matrix and calls the inversion routine.
C written by A. Klamt, Burscheider Str. 524, 5090 Leverkusen 3, Germany
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
      DIMENSION COORD(3,*)
      DIMENSION XX(3),XA(3),COBAS(3,NPPA,NUMATM)
      DIMENSION XSP(4,MAXNSS),NSET(MAXNSS,MAXNSS)
C     Added to satisfy the LAPACK call to DGETRI. SJC 1/10/93
      DIMENSION IPIV(MAXNSS)
      LOGICAL DIN(NPPA)
      COMMON / SOLVI / NSPA,NSS,IATSP(MAXNSS),NAR(MAXNSS),NN(2,NUMATM)
      COMMON / SOLVR / DIPL,FEPSI,RDS,DISEX2,TM(3,3,NUMATM),
     &                AAMAT(MAXNSS,MAXNSS),ADMAT(MAXNSS,MAXNSS),
     &                CCMAT(MAXORB,MAXORB), BBMAT(MAXORB,MAXNSS),
     &                COSURF(3,MAXNSS), SRAD(NUMATM)
      COMMON /DIRVEC/ DIRVEC(3,NPPA), DIRTM(3,NPPA)
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      COMMON /CHANEL/ IFILES(30)
      EQUIVALENCE (IW,IFILES(6))
      EQUIVALENCE (ADMAT,COBAS),(BBMAT,XSP)
      FDIAG=1.05D0*SQRT(NPPA+0.)
      AREA=0.0D0
      DS=SQRT(4./NSPA)
      C2DS=COS(2.0D0*DS)
      NSS1=0
      DO 1000 I=1,NUMAT
        DO 4 IX=1,3
4       XA(IX)=COORD(IX,I)
        R=SRAD(I)
        RM=R-RDS
        NSS0=NSS1+1
        IF (IATSP(NSS0) .NE. I) GO TO 1000
        IPS=NSS0
304     IF (IPS .GT. NSS .OR. IATSP(IPS) .NE. I) GO TO 306
        NAR(IPS)=0
        XSP(1,IPS)=0.D0
        XSP(2,IPS)=0.D0
        XSP(3,IPS)=0.D0
        XSP(4,IPS)=ADMAT(IPS,IPS)
        IPS=IPS+1
        GO TO 304
306     NSS1=IPS-1
C TRANSFORM COSURF ACCORDING TO TM(INV)
        DO 5 J=NSS0,NSS1
          XX(1)=COSURF(1,J)
          XX(2)=COSURF(2,J)
          XX(3)=COSURF(3,J)
          COSURF(1,J)=XX(1)*TM(1,1,I)+XX(2)*TM(1,2,I)+XX(3)*TM(1,3,I)
          COSURF(2,J)=XX(1)*TM(2,1,I)+XX(2)*TM(2,2,I)+XX(3)*TM(2,3,I)
          COSURF(3,J)=XX(1)*TM(3,1,I)+XX(2)*TM(3,2,I)+XX(3)*TM(3,3,I)
5       CONTINUE
C UPDATE TM
        NN1=NN(1,I)
        NN2=NN(2,I)
        IF (NN1 .EQ. 0) THEN
          TM(1,1,I)=1.D0
          TM(1,2,I)=0.D0
          TM(1,3,I)=0.D0
        ELSE
          DIST1=0.D0
          DO 6 IX=1,3
6         DIST1=DIST1+(XA(IX)-COORD(IX,NN1))**2
          DIST=1.D0/DSQRT(DIST1)
          TM(1,1,I)=(COORD(1,NN1)-XA(1))*DIST
          TM(1,2,I)=(COORD(2,NN1)-XA(2))*DIST
          TM(1,3,I)=(COORD(3,NN1)-XA(3))*DIST
        END IF
        IF (NN2 .EQ. 0) THEN
          TM(2,1,I)=-TM(1,2,I)
          TM(2,2,I)=TM(1,1,I)
          TM(2,3,I)=0.D0
        ELSE
          DIST2=0.D0
          DO 7 IX=1,3
7         DIST2=DIST2+(XA(IX)-COORD(IX,NN2))**2
          DIST=1.D0/DSQRT(DIST2)
          XX(1)=(COORD(1,NN2)-XA(1))*DIST
          XX(2)=(COORD(2,NN2)-XA(2))*DIST
          XX(3)=(COORD(3,NN2)-XA(3))*DIST
          SP=XX(1)*TM(1,1,I)+XX(2)*TM(1,2,I)+XX(3)*TM(1,3,I)
          SININV=1.D0/DSQRT(1.D0-SP*SP)
          TM(2,1,I)=(XX(1)-SP*TM(1,1,I))*SININV
          TM(2,2,I)=(XX(2)-SP*TM(1,2,I))*SININV
          TM(2,3,I)=(XX(3)-SP*TM(1,3,I))*SININV
        END IF
        TM(3,1,I)=TM(1,2,I)*TM(2,3,I)-TM(2,2,I)*TM(1,3,I)
        TM(3,2,I)=TM(1,3,I)*TM(2,1,I)-TM(2,3,I)*TM(1,1,I)
        TM(3,3,I)=TM(1,1,I)*TM(2,2,I)-TM(2,1,I)*TM(1,2,I)
C TRANSFORM COSURF ACCORDING TO TM
        DO 35 J=NSS0,NSS1
          XX(1)=COSURF(1,J)
          XX(2)=COSURF(2,J)
          XX(3)=COSURF(3,J)
          COSURF(1,J)=XX(1)*TM(1,1,I)+XX(2)*TM(2,1,I)+XX(3)*TM(3,1,I)
          COSURF(2,J)=XX(1)*TM(1,2,I)+XX(2)*TM(2,2,I)+XX(3)*TM(3,2,I)
          COSURF(3,J)=XX(1)*TM(1,3,I)+XX(2)*TM(2,3,I)+XX(3)*TM(3,3,I)
35      CONTINUE
        DO 45 J=1,NPPA
          XX(1)=DIRVEC(1,J)
          XX(2)=DIRVEC(2,J)
          XX(3)=DIRVEC(3,J)
          DO 45 IX=1,3
            X=XX(1)*TM(1,IX,I)+XX(2)*TM(2,IX,I)+XX(3)*TM(3,IX,I)
            DIRTM(IX,J)=X
            COBAS(IX,J,I)=RM*X+XA(IX)
45      CONTINUE
        DO 2080 J = 1,NPPA
          DIN(J)=.FALSE.
          DO 2020 IX=1,3
            XX(IX) = XA(IX) + DIRTM(IX,J)* R
2020      CONTINUE
          DO 2040 K = 1, NUMAT
            IF (K . EQ. I) GO TO 2040
            DIST=0.
            DO 2030 IX=1,3
              DIST = DIST + (XX(IX) - COORD(IX,K))**2
2030        CONTINUE
            DIST=SQRT(DIST)-SRAD(K)
            IF (DIST .LT. 0.D0) GO TO 2080
2040      CONTINUE
          DIN(J)=.TRUE.
2080    CONTINUE
        DO 350 J=1,NPPA
          IF (.NOT. DIN(J)) GO TO 350
          SPM=-1.D0
          X1=DIRTM(1,J)
          X2=DIRTM(2,J)
          X3=DIRTM(3,J)
          DO 340 IPS=NSS0,NSS1
            SP=X1*COSURF(1,IPS)+X2*COSURF(2,IPS)+X3*COSURF(3,IPS)
            IF (SP .LT. SPM) GO TO 340
            SPM=SP
            IPM=IPS
340       CONTINUE
          IF (SPM .LT. C2DS) GO TO 350
          NARA=NAR(IPM)+1
          NSET(IPM,NARA)=J
          NAR(IPM)=NARA
          DO 345 IX=1,3
345       XSP(IX,IPM)=XSP(IX,IPM)+DIRTM(IX,J)
350     CONTINUE
        DO 400 IPS=NSS0,NSS1
          DIST=0.
          IF (NAR(IPS) .EQ. 0) THEN
            GO TO 400
          END IF
          DO 390 IX=1,3
            X=XSP(IX,IPS)
            DIST=DIST+X*X
390       CONTINUE
          DIST=1.D0/SQRT(DIST)
          DO 391 IX=1,3
391       COSURF(IX,IPS)=XSP(IX,IPS)*DIST
400     CONTINUE
        DO 490 IPS=NSS0,NSS1
          DO 405 IX=1,3
405       XSP(IX,IPS)=XA(IX)+COSURF(IX,IPS)*RM
490     CONTINUE
1000  CONTINUE
C
C filling AAMAT
      DO 170 IPS=1,NSS
        I=IATSP(IPS)
        RI=SRAD(I)-RDS
        NARI=NAR(IPS)
        AA=0.D0
        DO 140 K=1,NARI
          J1=NSET(IPS,K)
          AA=AA+.5*FDIAG
          X1=DIRVEC(1,J1)
          X2=DIRVEC(2,J1)
          X3=DIRVEC(3,J1)
          DO 140 L=1,K-1
            J2=NSET(IPS,L)
            AA=AA+1./SQRT((X1-DIRVEC(1,J2))**2+
     &             (X2-DIRVEC(2,J2))**2+(X3-DIRVEC(3,J2))**2)
140     CONTINUE
        AA=2*AA/RI
        IF (NARI .EQ. 0) THEN
          AAMAT(IPS,IPS)=XSP(4,IPS)
        ELSE
          AAMAT(IPS,IPS)=AA
        ENDIF
        DO 141 IX=1,3
141     XA(IX)=XSP(IX,IPS)
        DO 169 JPS=1,IPS-1
          NARJ=NAR(JPS)
          DIST=0.
          DO 143 IX=1,3
143       DIST=DIST+(XSP(IX,JPS)-XA(IX))**2
          IF (DIST .LT. DISEX2) THEN
            J=IATSP(JPS)
            AIJ=0.D0
            DO 149 K=1,NARI
              J1=NSET(IPS,K)
              X1=COBAS(1,J1,I)
              X2=COBAS(2,J1,I)
              X3=COBAS(3,J1,I)
              DO 149 L=1,NARJ
                J2=NSET(JPS,L)
                AIJ=AIJ+1.D0/SQRT((COBAS(1,J2,J)-X1)**2+
     &          (COBAS(2,J2,J)-X2)**2+(COBAS(3,J2,J)-X3)**2)
149         CONTINUE
          ELSE
            AIJ=NARI*NARJ/SQRT(DIST)
          END IF
          AAMAT(IPS,JPS)=AIJ
          AAMAT(JPS,IPS)=AIJ
169     CONTINUE
170   CONTINUE
      DO 175 I=1,NSS
        DO 175 J=1,NSS
175   ADMAT(I,J)=AAMAT(I,J)
C invert AAMAX
C     CALL MA22BD(AAMAT,MAXNSS,NSS,BBMAT,ERROR)
C     CALL INVAA(AAMAT,MAXNSS,NSS)
      CALL DGETRF(NSS,NSS,AAMAT,MAXNSS,IPIV,INFO)
      IF (INFO .NE. 0) THEN
         WRITE (IW,*) 'UPSURF: Factorization failed. INFO=',INFO
         STOP
      ENDIF
      CALL DGETRI(NSS,AAMAT,MAXNSS,IPIV,BBMAT,MAXORB*MAXNSS,INFO)
      IF (INFO .NE. 0) THEN
         WRITE (IW,*) 'UPSURF: Inversion failed. INFO=',INFO
         STOP
      ENDIF
      RETURN
      END

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive


text:
!-------------------------------------------------------------------------------------------------!
!> Processig command line option
!!
!! @copyright
!! Copyright (c) 2019 Takuto Maeda. All rights reserved.
!!
!! @license
!! This software is released under the MIT license. See LICENSE for details.
!<
!--
module m_getopt

  use iso_fortran_env
  use m_system
  implicit none
  private

  public :: getopt

  !-----------------------------------------------------------------------------------------------!
  !> Obtain command-line option
  !!
  !! @par Usage
  !! getopt( OPTNAME, exist, var, default_bar )
  !!   - var, default_var are optional.
  !!   - types of var can be double, single integer or character
  !--
  interface getopt

   module procedure getopt_d, getopt_s, getopt_i, getopt_c

  end interface getopt
  !-----------------------------------------------------------------------------------------------!

contains

  !-----------------------------------------------------------------------------------------------!
  !> Obtain command line as characters
  !!
  !! @par Example
  !!   For the command
  !!   hoge -A -T 01
  !!     call getopt( 'A', isA, Aval ) -> isA=.true.,  Aval=''
  !!     call getopt( 'T', isT, Tval ) -> isT=.true.,  Tval='01'
  !!     call getopt( 'C', isC, Cval ) -> isC=.false., Cval=''
  !!     call getopt( 'A', isA )       -> isA=.true.
  !--
  subroutine getopt_c(opt, isExist, val, default)
    character(*),           intent(in)  :: opt
    logical,                intent(out) :: isExist
    character(*), optional, intent(out) :: val
    character(*), optional, intent(in)  :: default
    !--
    integer :: narg
    character(256), allocatable :: argv(:)
    integer :: i
    character(256) :: optkey
    !----

    narg = system__iargc()
    allocate( argv(1:narg) )

    do i=1, narg
       call system__getarg( i, argv(i)(:) )
    end do

    if( present( val ) ) then
       val = ''
    end if

    isExist = .false.
    optkey = '-'//trim(adjustl( opt ) )

    do i=1, narg

       if( trim(optkey) == trim(argv(i)) ) then

          if( isExist ) then
             write(error_unit,*) 'getopt: ', &
                  'option '//trim(optkey)//' is multiplly defined. '
          end if

          isExist = .true.

          if( present( val ) ) then

             if( i==narg ) then
                val = ''
             else
                val = argv(i+1)
             end if

          end if
       end if
    end do

    deallocate( argv )

    if( .not. isExist .and. present( default ) ) then
       val = default
    end if

  end subroutine getopt_c
  !-----------------------------------------------------------------------------------------------!

  !-----------------------------------------------------------------------------------------------!
  !> Obtain command line option for double precision number
  !--
  subroutine getopt_d( opt, isExist, val, default )

    character(*),       intent(in)  :: opt
    logical,            intent(out) :: isExist
    real(real64),           intent(out) :: val
    real(real64), optional, intent(in)  :: default
    !--
    character(1024) :: aval
    !----

    call getopt_c( opt, isExist, aval )

    if( .not. isExist ) then
       if( present( default ) ) then
          val = default
       else
          val = -99999.9_real64
       end if
       return
    end if

    read( aval, * ) val

  end subroutine getopt_d
  !-----------------------------------------------------------------------------------------------!

  !-----------------------------------------------------------------------------------------------!
  !> Obtain command line option for single precision number
  !--
  subroutine getopt_s( opt, isExist, val, default )

    character(*),   intent(in)  :: opt
    logical,        intent(out) :: isExist
    real,           intent(out) :: val
    real, optional, intent(in)  :: default
    !--
    character(1024) :: aval
    !----

    call getopt_c( opt, isExist, aval )

    if( .not. isExist ) then
       if( present( default ) ) then
          val = default
       else
          val = -99999.9
       end if
       return
    end if

    read(aval, *) val

  end subroutine getopt_s
  !-----------------------------------------------------------------------------------------------!

  !-----------------------------------------------------------------------------------------------!
  !> Obtain command line option for integer number
  !--
  subroutine getopt_i( opt, isExist, val, default )

    character(*),      intent(in)  :: opt
    logical,           intent(out) :: isExist
    integer,           intent(out) :: val
    integer, optional, intent(in)  :: default
    !--
    character(1024) :: aval
    !----

    call getopt_c( opt, isExist, aval )

    if( .not. isExist ) then
       if( present( default ) ) then
          val = default
       else
          val = -99999
       end if
       return
    end if

    read(aval, *) val

  end subroutine getopt_i
  !-----------------------------------------------------------------------------------------------!

end module m_getopt
!-------------------------------------------------------------------------------------------------!

max_stars_repo_name: tktmyd/fwin

max_issues_repo_name: tktmyd/fwin

max_forks_repo_name: tktmyd/fwin


text:
      subroutine bacteria

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine calculates bacteria growth, transport with runoff and
!!    loss due to percolation into soil

!!    ~ ~ ~ INCOMING VARIABLES ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    bactkdq     |none          |Bacteria partition coefficient.
!!                               |Partition coefficient for bacteria between
!!                               |soluble and sorbed phase in surface runoff.
!!    bactlp_plt(:)|# cfu/m^2    |less persistent bacteria on foliage
!!    bactlpq(:)  |# cfu/m^2     |less persistent bacteria in soil solution
!!    bactlps(:)  |# cfu/m^2     |less persistent bacteria attached to soil
!!                               |particles
!!    bactminlp   |# cfu/m^2     |Threshold detection level for less persistent
!!                               |bacteria
!!                               |when bacteria levels drop to this amount the
!!                               |model considers bacteria in the soil to be
!!                               |insignificant and sets the levesl to zero
!!    bactminp    |# cfu/m^2     |Threshold detection level for persistent
!!                               |bacteria
!!                               |when bacteria levels drop to this amount the
!!                               |model considers bacterial in the soil to be
!!                               |insignificant and sets the levels to zero
!!    bactmx      |none          |bacteria percolation coefficient
!!                               |Ratio of solution bacteria in surface layer
!!                               |to solution bacteria in percolate
!!    bactp_plt(:)|# cfu/m^2     |persistent bacteria on foliage
!!    bactpq(:)   |# cfu/m^2     |persistent bacteria in soil solution
!!    bactps(:)   |# cfu/m^2     |persistent bacteria attached to soil particles
!!    curyr       |none          |current year of simulation
!!    enratio     |none          |enrichment ratio calculated for current day
!!                               |in HRU
!!    filterw(:)  |m             |filter strip width for bacteria transport
!!    hru_dafr(:) |none          |fraction of watershed area in HRU
!!    ihru        |none          |HRU number
!!    nyskip      |none          |number of years to skip output summarization
!!                               |and printing
!!    precipday   |mm H2O        |precipitation for the day in HRU
!!    sbactlchlp  |# cfu/m^2     |average annual number of less persistent
!!                               |bacteria lost from soil surface layer by
!!                               |percolation
!!    sbactlchp   |# cfu/m^2     |average annual number of persistent bacteria
!!                               |lost from soil surface layer by percolation
!!    sdiegrolpq  |# cfu/m^2     |average annual change in the number of
!!                               |less persistent bacteria colonies in soil
!!                               |solution in watershed
!!    sdiegrolps  |# cfu/m^2     |average annual change in the number of
!!                               |less persistent bacteria colonies on soil
!!                               |particles in watershed
!!    sdiegropq   |# cfu/m^2     |average annual change in the number of
!!                               |persistent bacteria colonies in soil solution
!!                               |in watershed
!!    sdiegrops   |# cfu/m^2     |average annual change in the number of
!!                               |persistent bacteria colonies on soil particles
!!                               |in watershed
!!    sedyld(:)   |metric tons   |daily soil loss caused by water erosion
!!    sol_bd(:,:) |Mg/m**3       |bulk density of the soil
!!    sol_z(:,:)  |mm            |depth to bottom of soil layer
!!    surfq(:)    |mm H2O        |surface runoff generated on day in HRU
!!    thbact      |none          |temperature adjustment factor for bacteria
!!                               |die-off/growth
!!    tmpav(:)    |deg C         |average air temperature on current day in HRU
!!    wlpq20      |1/day         |Overall rate change for less persistent
!!                               |bacteria in soil solution.
!!    wlps20      |1/day         |Overall rate change for less persistent
!!                               |bacteria adsorbed to soil particles.
!!    wof_lp      |none          |fraction of less persistent bacteria on foliage
!!                               |that is washed off by a rainfall event
!!    wof_p       |none          |fraction of persistent bacteria on foliage that
!!                               |is washed off by a rainfall event
!!    wp20lp_plt  |1/day         |Overall rate change for less persistent bacteria
!!                               |on foliage
!!    wp20p_plt   |1/day         |Overall rate change for persistent bacteria on
!!                               |foliage
!!    wpq20       |1/day         |Overall rate change for persistent bacteria in
!!                               |soil solution.
!!    wps20       |1/day         |Overall rate change for persistent bacteria
!!                               |adsorbed to soil particles.
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ OUTGOING VARIABLES ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    bactlchlp   |# cfu/m^2     |less persistent bacteria removed from soil
!!                               |surface layer by percolation
!!    bactlchp    |# cfu/m^2     |persistent bacteria removed from soil surface
!!                               |layer by percolation
!!    bactlp_plt(:)|# cfu/m^2    |less persistent bacteria on foliage
!!    bactlpq(:)  |# cfu/m^2     |less persistent bacteria in soil solution
!!    bactlps(:)  |# cfu/m^2     |less persistent bacteria attached to soil
!!                               |particles
!!    bactp_plt(:)|# cfu/m^2     |persistent bacteria on foliage
!!    bactpq(:)   |# cfu/m^2     |persistent bacteria in soil solution
!!    bactps(:)   |# cfu/m^2     |persistent bacteria attached to soil particles
!!    bactrolp    |# cfu/m^2     |less persistent bacteria transported to main
!!                               |channel with surface runoff
!!    bactrop     |# cfu/m^2     |persistent bacteria transported to main
!!                               |channel with surface runoff
!!    bactsedlp   |# cfu/m^2     |less persistent bacteria transported with
!!                               |sediment in surface runoff
!!    bactsedp    |# cfu/m^2     |persistent bacteria transported with
!!                               |sediment in surface runoff
!!    sbactlchlp  |# cfu/m^2     |average annual number of less persistent
!!                               |bacteria lost from soil surface layer by
!!                               |percolation
!!    sbactlchp   |# cfu/m^2     |average annual number of persistent bacteria
!!                               |lost from soil surface layer by percolation
!!    sdiegrolpq  |# cfu/m^2     |average annual change in the number of
!!                               |less persistent bacteria colonies in soil
!!                               |solution in watershed
!!    sdiegrolps  |# cfu/m^2     |average annual change in the number of
!!                               |less persistent bacteria colonies on soil
!!                               |particles in watershed
!!    sdiegropq   |# cfu/m^2     |average annual change in the number of
!!                               |persistent bacteria colonies in soil solution
!!                               |in watershed
!!    sdiegrops   |# cfu/m^2     |average annual change in the number of
!!                               |persistent bacteria colonies on soil particles
!!                               |in watershed
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ LOCAL DEFINITIONS ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    bactmn      |
!!    bactmx      |
!!    blpq        |# cfu/m^2     |less persistent bacteria in soil solution at
!!                               |beginning of day
!!    blps        |# cfu/m^2     |less persistent bacteria attached to soil
!!                               |particles at beginning of day
!!    bpq         |# cfu/m^2     |persistent bacteria in soil solution at
!!                               |beginning of day
!!    bps         |# cfu/m^2     |persistent bacteria attached to soil particles
!!                               |at beginning of day
!!    cbact       |
!!    j           |none          |HRU number
!!    wt1         |none          |conversion factor to convert kg/ha to g/t(ppm)
!!    xx          |
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ SUBROUTINES/FUNCTIONS CALLED ~ ~ ~
!!    Intrinsic: Exp, Min, Max
!!    SWAT: Theta

!!    ~ ~ ~ ~ ~ ~ END SPECIFICATIONS ~ ~ ~ ~ ~ ~

      use parm

      integer :: j
      real :: bpq, blpq, bps, blps, wt1, cbact, xx

      j = 0
      j = ihru

      if (bactlps(j) < 1.e-6) bactlps(j) = 0.0
      if (bactlpq(j) < 1.e-6) bactlpq(j) = 0.0
      if (bactpq(j) < 1.e-6) bactpq(j) = 0.0
      if (bactps(j) < 1.e-6) bactps(j) = 0.0
      if (bactp_plt(j) < 1.e-6) bactp_plt(j) = 0.0

!! compute bacteria wash off
      if (precipday >= 2.54) then
        xx = 0.
        xx = wof_p * bactp_plt(j)
        if (xx > bactp_plt(j)) xx = bactp_plt(j)
        bactpq(j) = bactpq(j) + xx
        bactp_plt(j) = bactp_plt(j) - xx
        xx = 0.
        xx = wof_lp * bactlp_plt(j)
        if (xx > bactlp_plt(j)) xx = bactlp_plt(j)
        bactlpq(j) = bactlpq(j) + xx
        bactlp_plt(j) = bactlp_plt(j) - xx
      end if

!! compute bacteria die-off and re-growth on foilage
      if (tmpav(j) > 1.e-6) then
        bactp_plt(j) = bactp_plt(j) * Exp(-Theta(wp20p_plt,thbact,
     &    tmpav(j))) - bactminp
        bactp_plt(j) = Max(0., bactp_plt(j))
        if (bactp_plt(j) < bactminp) bactp_plt(j) = 0.
        bactlp_plt(j) = bactlp_plt(j) * Exp(-Theta(wp20lp_plt,thbact,
     &    tmpav(j))) - bactminlp
        bactlp_plt(j) = Max(0., bactlp_plt(j))
        if (bactlp_plt(j) < bactminlp) bactlp_plt(j) = 0.
      endif

!! compute bacteria die-off and re-growth in surface soil layer
      bpq = 0.
      blpq = 0.
      bps = 0.
      blps = 0.
      bpq = bactpq(j)
      bactpq(j) = bactpq(j) * Exp(-Theta(wpq20,thbact,tmpav(j))) -
     &    bactminp
      bactpq(j) = Max(0., bactpq(j))
      if (bactpq(j) < bactminp) bactpq(j) = 0.
      blpq = bactlpq(j)
      bactlpq(j) = bactlpq(j) * Exp(-Theta(wlpq20,thbact,tmpav(j))) -
     &    bactminlp
      bactlpq(j) = Max(0., bactlpq(j))
      if (bactlpq(j) < bactminlp) bactlpq(j) = 0.
      bps = bactps(j)
      bactps(j) = bactps(j) * Exp(-Theta(wps20,thbact,tmpav(j))) -
     &    bactminp
      bactps(j) = Max(0., bactps(j))
      if (bactps(j) < bactminp) bactps(j) = 0.
      blps = bactlps(j)
      bactlps(j) = bactlps(j) * Exp(-Theta(wlps20,thbact,tmpav(j))) -
     &    bactminlp
      bactlps(j) = Max(0., bactlps(j))
      if (bactlps(j) < bactminlp) bactlps(j) = 0.

!! compute bacteria in the runoff
      bactrop = bactpq(j) * surfq(j) /
     &                        (sol_bd(1,j) * sol_z(1,j) * bactkdq)
      bactrop = Min(bactrop, bactpq(j))
      bactrop = Max(bactrop, 0.)
      bactpq(j) = bactpq(j) - bactrop

      bactrolp = bactlpq(j) * surfq(j) /
     &                          (sol_bd(1,j) * sol_z(1,j) * bactkdq)
      bactrolp = Min(bactrolp, bactlpq(j))
      bactrolp = Max(bactrolp, 0.)
      bactlpq(j) = bactlpq(j) - bactrolp

!! compute bacteria transported with sediment
      if (enratio > 0.) then
        wt1 = 0.
        wt1 = sol_bd(1,j) * sol_z(1,j) / 1000.

        cbact = 0.
        cbact = bactps(j) * enratio / wt1
        bactsedp = .0001 * cbact * sedyld(j) / hru_ha(j)
        bactsedp = Min(bactsedp, bactps(j))
        bactps(j) = bactps(j) - bactsedp

        cbact = 0.
        cbact = bactlps(j) * enratio / wt1
        bactsedlp = .0001 * cbact * sedyld(j) / hru_ha(j)
        bactsedlp = Min(bactsedlp, bactlps(j))
        bactlps(j) = bactlps(j) - bactsedlp
      end if

!! compute bacteria incorporated into the soil
      bactlchp = bactpq(j) * sol_prk(1,j) / ((conv_wt(1,j) / 1000.)
     &                                                       * bactmx)
      bactlchp = Min(bactlchp, bactpq(j))
      bactlchp = Max(bactlchp, 0.)
      bactpq(j) = bactpq(j) - bactlchp

      bactlchlp = bactlpq(j) * sol_prk(1,j) / ((conv_wt(1,j) / 1000.)
     &                                                       * bactmx)
      bactlchlp = Min(bactlchlp, bactlpq(j))
      bactlchlp = Max(bactlchlp, 0.)
      bactlpq(j) = bactlpq(j) - bactlchlp

!! summary calculations
      if (curyr > nyskip) then
        sdiegropq = sdiegropq + (bpq - bactpq(j)) * hru_dafr(j)
        sdiegrolpq = sdiegrolpq + (blpq - bactlpq(j)) * hru_dafr(j)
        sdiegrops = sdiegrops + (bps - bactps(j)) * hru_dafr(j)
        sdiegrolps = sdiegrolps + (blps - bactlps(j)) * hru_dafr(j)
!! added 4 here
        sbactrop = sbactrop + bactrop * hru_dafr(j)
        sbactrolp = sbactrolp + bactrolp * hru_dafr(j)
        sbactsedp = sbactsedp + bactsedp * hru_dafr(j)
        sbactsedlp = sbactsedlp +  bactsedlp * hru_dafr(j)
!! added 4 here
        sbactlchp = sbactlchp + bactlchp * hru_dafr(j)
        sbactlchlp = sbactlchlp + bactlchlp * hru_dafr(j)
      end if

! 1 is HRU number!
!      xx = bactpq(1) + bactps(1) + bactp_plt(1)
!      yy = bactlpq(1) + bactlps(1) + bactlp_plt(1)
!      write (17,100) iida, xx, yy,
!     &     bactpq(1), bactps(1), bactlpq(1), bactlps(1),
!     &     bactrop, bactrolp, bactsedp, bactsedlp, bactlchp, bactlchlp,
!     &     bactp_plt(1), bactlp_plt(1)
! 100  format (i4,14f10.7)

      return
      end


! chk_dir_separator.f90 --
!     Check if a program can use the UNIX directory separator
!
!     Background:
!     If you work on Windows and UNIX/Linux-like operating systems, it
!     is annoying to have to distinguish between forward and backward slashes
!     in file names. If the forward slash is accepted, then you only need
!     a single version.
!
program chk_dir_separator
    implicit none

    integer :: lun, ierr

    lun = 10

    open( lun, file = 'chk_dir_separator.inp' )
    write( lun, '(a)' ) 'Text'
    close( lun )

    open( lun, file = './chk_dir_separator.inp', status = 'old', iostat = ierr )
    if ( ierr /= 0 ) then
        write( *, '(a)' ) 'Open statement on chk_dir_separator.inp failed - ', &
                          '    apparently the forward slash is not seen as a directory separator'
    else
        write( *, '(a)' ) 'Open statement on chk_dir_separator.inp succeeded -'
        write( *, '(a)' ) '    forward slashes acceptable as directory separator'
    endif
end program chk_dir_separator


text:
 Venues for conferences, bands, meetings, gatherings
 This is meant to be a comprehensive list of spaces for you to hold your event, be it a weddings wedding, conference, meeting, trade show, banquet, musical performance, film screening, or other gathering.
Also see the Yolo County Visitors Bureau.

 City of Davis venues are administered by the http://cityofdavis.org/cs/facility City of Davis Community Services

|| The Attendance Office

 UC Davis venues are administered by UC Davis Conference and Event Services

Capacity, unless otherwise stated, is for seated (i.e. for lectures). Increase for standing (reception), decrease for dining (seated at tables).
Definitions used here:
  Auditorium: Large room with a downward slope and fixed seating
  Multipurpose room: Large empty room to which chairs, tables may be brought. Good for bands, dances, movies, presentations, banquets.

 Rural Venues near Davis
  Plainfield Station
  Yolo Land & Cattle Local cattle ranch rents places for weddings, meetings, banquets, etc.
  http://www.taberranch.com/ Taber Ranch Vineyards and Event Center

 Other facilities to consider
  Athletic fields
  Churches (may have restrictions)
  Hotel conference rooms (may be costly)
  Parks (possible rain, mosquitos, harsh sun, no kitchen, usually no electricity)
  Private homes (free, few restrictions, comfortable atmosphere, small, but risk from strangers)
  Restaurants (loud, limited time, pressure to purchase food)

 Other things to consider, depending on the type of event
Calendar conflicts, food (catering? potluck? brown bag? none?), transportation, lodging, deposits & insurance, labor (volunteer? paid?), permits, noise permits, A/V equipment/stage/lights, electric power charges, security, cleanup, waste & garbage.


text:
*----------------------------------------------------------------------*
      subroutine mem_iget(ffbuf,ibuf,idxst,idxnd)
*----------------------------------------------------------------------*
*     get integer_array(idxst:idxnd) on ibuf(1:)
*----------------------------------------------------------------------*
      use memman
      implicit none
c      include 'def_filinf.h'

      type(filinf) ::
     &     ffbuf
      integer, intent(in) ::
     &     idxst, idxnd
      integer, intent(out) ::
     &     ibuf(*)

      call memman_vbuffer_get_int(ffbuf%buf_id,idxst,idxnd,ibuf)

      return
      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


! { dg-do run }
! { dg-options "-frepack-arrays" }
!
! Check that arrays marked with TARGET attribute are not repacked.
!
program test2
   use iso_c_binding
   implicit none
   real, target :: x(7)
   type(c_ptr) cp1, cp2

   x = 42
   if (.not. c_associated(c_loc(x(3)),point(x(::2)))) STOP 1
contains
  function point(x)
    use iso_c_binding
    real, intent(in), target :: x(:)
    type(c_ptr) point
    real, pointer :: p

    p => x(2)
    point = c_loc(p)
  end function point
end program test2


text:
c-------------------------------------------------------------------------
      subroutine H_ADJUST(ZUs,ZTs,ZQs,U_hs,T_hs,Q_hs,IHUMID)
c
c This subroutine adjusts the U,T,Q variables to the specified
c standard height (ZUs,ZTs,ZQs) using the loglayer profiles.
c The DELTA correction (adjustment) is relative to the surface
c measurement.             Cronin 4/13/94
c
      real*8 ZUs,ZTs,ZQs,U_hs,T_hs,Q_hs,ZUsL,ZTsL,ZQsL
      real*8 U,T,Q,TS,QS,rns,rnl,ZU,ZT,ZQ,zi,P,PUZs,PTZs,PQZs
      real*8 U_wg_hs,Rho_hs,Rho_avg,QA,Rho,P_hs,ee
      real*8 USR,TSR,QSR,ZO,zot,zoq,ZL,RR,RT,RQ,RI,dter,dqer,tkt
      real*8 al,beta,cpa,cpw,grav,xlv,rhoa,rhow,rgas,toK,visa
      real*8 visw,von,fdg,DU_Wg,Wg,S,D
      integer ID,ihumid
c
      COMMON/PIN/U,T,Q,TS,QS,rns,rnl,ZU,ZT,ZQ,zi,P,ID
      COMMON/POUT/USR,TSR,QSR,ZO,zot,zoq,ZL,RR,RT,RQ,RI,dter,dqer,tkt
      COMMON/const/al,beta,cpa,cpw,grav,xlv,rhoa,rhow,rgas,toK,
     &             visa,visw,von,fdg
      COMMON/wgust/DU_Wg,Wg

      call ZETA(T,Q,USR,TSR,QSR,ZUs,ZUsL)
      call ZETA(T,Q,USR,TSR,QSR,ZTs,ZTsL)
      call ZETA(T,Q,USR,TSR,QSR,ZQs,ZQsL)
      PUZs= PSI(1,ZUsL)
      PTZs= PSI(2,ZTsL)
      PQZs= PSI(2,ZQsL)

      S = (dlog(ZTs*USR/VISA/RT)-PTZs)/(von*fdg)
      D = (dlog(ZQs*USR/VISA/RQ)-PQZs)/(von*fdg)
      T_hs =TSR*S +TS - dter -.0098*ZTs
      Q_hs =(QSR*D + QS - dqer)*1000
      U_wg_hs = USR*(dlog(ZUs/ZO) - PUZs)/0.4
      if(U_wg_hs.ge.Wg) then
         U_hs = SQRT(U_wg_hs**2 - Wg**2)
      else
         U_hs = U_wg_hs
      endif
c
      if(IHUMID.eq.1) then    ! then need to convert sp hum into rh
         Q_hs = Q_hs/1000     ! sh kg/kg
         RHO=1./(287.*(T+273.16)*(1.+.61*Q))*P*100.
         P_hs = P - (RHO*grav*(ZTs - ZT))/100 !Approx hydrost.Pressure mb
         RHO_hs=1./(287.*(T_hs+273.16)*(1.+.61*Q_hs))*P_hs*100
         RHO_avg = (RHO + RHO_hs)/2
         P_hs = P -(RHO_avg*grav*(ZTs - ZT))/100 !hydrostatic Pressure
         call humidity(T_hs,P_hs,QA)         !Teten's formula for Pvap,sat
         ee=Q_hs*P_hs/(.62197 + .378*Q_hs)   !to get vapor pressure
         Q_hs = ee/QA                        !to get relative humidity
      endif
      return
      end


text:
!------------------------------------------------------------------------------!
!  The Community Multiscale Air Quality (CMAQ) system software is in           !
!  continuous development by various groups and is based on information        !
!  from these groups: Federal Government employees, contractors working        !
!  within a United States Government contract, and non-Federal sources         !
!  including research institutions.  These groups give the Government          !
!  permission to use, prepare derivative works of, and distribute copies       !
!  of their work in the CMAQ system to the public and to permit others         !
!  to do so.  The United States Environmental Protection Agency                !
!  therefore grants similar permission to use the CMAQ system software,        !
!  but users are requested to provide copies of derivative works or            !
!  products designed to operate in the CMAQ system to the United States        !
!  Government without restrictions as to use by others.  Software              !
!  that is used with the CMAQ system but distributed under the GNU             !
!  General Public License or the GNU Lesser General Public License is          !
!  subject to their copyright restrictions.                                    !
!------------------------------------------------------------------------------!

FUNCTION ptemp (temp, press)

!-------------------------------------------------------------------------------
! Name:     Potential Temperature
! Purpose:  Returns potential temperature [K] as a function of temperature [K]
!           and pressure [Pa] using Poisson's equation.
! Revised:  ?? ??? ????  Original version as a statement function in MCIP
!                        routine getpblht.f90.
!           23 Feb 2011  Converted to independent routine.  (T. Otte)
!           07 Sep 2011  Updated disclaimer.  (T. Otte)
!-------------------------------------------------------------------------------

  IMPLICIT NONE

  REAL,          INTENT(IN)    :: press      ! pressure [Pa]
  REAL                         :: ptemp      ! virtual temperature [K]
  REAL,          PARAMETER     :: rdovcp     = 2.0 / 7.0  ! Rd / cP
  REAL,          INTENT(IN)    :: temp       ! temperature [K]

  ptemp = temp * ( 100000.0 / press )**rdovcp

END FUNCTION ptemp

max_stars_repo_name: Simeng-unique/CMAQ-changed

max_issues_repo_name: Simeng-unique/CMAQ-changed

max_forks_repo_name: Simeng-unique/CMAQ-changed


module dlg

contains
    function dlg_pDeriv ( x, y, array, dir ) result ( dArray )
        implicit none

        real, dimension (:,:), allocatable :: dArray
        integer :: nX, nY, i, j
        integer, intent(IN) :: dir
        real, intent(IN) :: array (:,:), x(:), y(:)
        real :: dS

        nX  = size ( array, 1 )
        nY  = size ( array, 2 )

        allocate ( dArray ( nX, nY ) )
        dArray  = 0.0

        if ( dir == 2 ) then

            do i=1,nX
                do j=1,nY

                    if ( j > 1 .AND. j < nY ) then

                        dS = ((y(j+1)-y(j)) + (y(j)-y(j-1))) / 2.0

                        dArray(i,j) = 1.0 / ( 2.0 * dS ) * &
                            ( array(i,j+1) - array(i,j-1) )

                    elseif ( j == 1 ) then

                        dS = ((y(j+1)-y(j)) + (y(j+2)-y(j+1))) / 2.0

                        dArray(i,j) = 1.0 / ( 2.0 * dS ) * &
                            ( &
                            -3.0 * array(i,j) + 4.0 * array(i,j+1) &
                            - array(i,j+2) )

                    elseif ( j == nY ) then

                        dS = ((y(j-1)-y(j-2)) + (y(j)-y(j-1))) / 2.0

                        dArray(i,j) = 1.0 / ( 2.0 * dS ) * &
                            ( &
                            3.0 * array(i,j) - 4.0 * array(i,j-1) &
                            + array(i,j-2) )

                    endif

                enddo
            enddo

        elseif ( dir == 1 ) then

            do i=1,nX
                do j=1,nY

                    if ( i > 1 .AND. i < nX ) then

                        dS = ((x(i+1)-x(i)) + (x(i)-x(i-1))) / 2.0

                        dArray(i,j) = 1.0 / ( 2.0 * dS ) * &
                            ( array(i+1,j) - array(i-1,j) )

                    elseif ( i == 1 ) then

                        dS = ((x(i+1)-x(i)) + (x(i+2)-x(i+1))) / 2.0

                        dArray(i,j) = 1.0 / ( 2.0 * dS ) * &
                            ( &
                            -3.0 * array(i,j) + 4.0 * array(i+1,j) &
                            - array(i+2,j) )

                    elseif ( i == nX ) then

                        dS = ((x(i-1)-x(i-2)) + (x(i)-x(i-1))) / 2.0

                        dArray(i,j) = 1.0 / ( 2.0 * dS ) * &
                            ( &
                            3.0 * array(i,j) - 4.0 * array(i-1,j) &
                            + array(i-1,j) )

                    endif

                enddo
            enddo

        end if

    end function dlg_pDeriv

end module dlg


text:
	SUBROUTINE GG_TCDA ( rlat, rlon, f34kt, ocean, radii, elat,
     +			     elon, iret )
C************************************************************************
C* GG_TCDA								*
C*									*
C* This subroutine plots the danger area for a single tropical system.  *
C* The danger area is determined by adding 50-100-150-200-300 nm (for   *
C* the TPC) or 50-80-125-180-250 nm (for the CPHC) to the 12-24-36-48-72*
C* hour forecast position and maximum 34kt wind radii respectively, and *
C* by the 34kt wind radii at the current position.                      *
C*									*
C* GG_TCDA ( RLAT, RLON, F34KT, OCEAN, RADII, ELAT, ELON, IRET )        *
C*									*
C* Input parameters:							*
C*	RLAT (6)    	REAL	    Track point latitudes               *
C*	RLON (6)   	REAL	    Track point longitudes              *
C*	F34KT (6)	CHAR*	    Current and fcst 34kt wind radii    *
C*	OCEAN		CHAR*	    Basin indicator - AL, EP or CP      *
C*	RADII (6)	REAL	    Values to be added to 34kt winds    *
C*									*
C* Output parameters:							*
C*	ELAT (4)	REAL	    Latitudes of danger area extrema    *
C*	ELON (4)	REAL	    Longitudes of danger area extrema   *
C*	IRET		INTEGER	    Return code			        *
C*									*
C**									*
C* Log:									*
C* D. Kidwell/NCEP	 7/01	                                        *
C* D. Kidwell/NCEP	 2/02	Rewrote to improve interpolation &      *
C*				smoothing; added CPHC and extrema       *
C* D. Kidwell/NCEP	 5/03	Added argument iptsm to GG_TCSM call    *
C* D. Kidwell/NCEP	 5/04	Initialized elat, elon                  *
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	CHARACTER*(*)	f34kt (*), ocean
	REAL		rlat (*), rlon (*), elat (*), elon (*),
     +			radii (*)
C*
C*	MXPTS should have the same value as MAXPTS in gemprm.h
C*
	PARAMETER	( ITRACK = 6, MXINTP = 1000, MXWORK = 10000,
     +			  MXPTS = 500, IANG1 = 0, IANG2 = 90 ,
     +			  IPOLY = ( (IANG2-IANG1+1) * 4 ),
     +			  RADCUR = 40. )
C*
	CHARACTER	ddev*2, sys*2
	INTEGER		intarr (5), ino12 (IPOLY), ino21 (IPOLY),
     +			ipt1 (2,100), ipt2 (2,100), iptsv (ITRACK)
	LOGICAL		good (ITRACK), out
	REAL		tlat (ITRACK), tlon (ITRACK), quad (4,ITRACK),
     +			angle (4), dist (5),
     +			xlat (IPOLY,ITRACK), xlon (IPOLY,ITRACK),
     +			dlat (MXINTP), dlon (MXINTP),
     +			dse (MXINTP), dne (MXINTP), dnw (MXINTP),
     +			dsw (MXINTP), wlat (MXWORK), wlon (MXWORK),
     +			olat (MXWORK), olon (MXWORK),
     +			blat (IPOLY), blon (IPOLY)
C*
	INCLUDE		'ERMISS.FNC'
C*
	DATA		angle / 45., 135., 225., 315. /
C-----------------------------------------------------------------------
	iret = 0
	DO ii = 1, 4
	    elat ( ii ) = RMISSD
	    elon ( ii ) = RMISSD
	END DO
C
	CALL GQLINE ( jltyp, jlthw, jwidth, jlwhw, ier )
C
	DO ii = 1, ITRACK
	    tlat ( ii ) = rlat ( ii )
	    tlon ( ii ) = rlon ( ii )
	END DO
C
C*	Check for systems above 60N (Atlantic) or 40N (Pacific).
C
	IF ( ocean .eq. 'AL' ) THEN
	    xmxlat = 60.
	  ELSE
	    xmxlat = 40.
	END IF
	IF ( tlat ( ITRACK ) .gt. xmxlat ) tlat ( ITRACK ) = RMISSD
	IF ( ERMISS ( tlat ( ITRACK ) ) ) THEN
	    out = .true.
	    ii  = ITRACK - 1
	    DO WHILE ( out )
		IF ( tlat ( ii ) .gt. xmxlat ) THEN
		    tlat ( ii ) = RMISSD
		  ELSE
		    IF ( .not. ERMISS ( tlat ( ii ) ) ) out = .false.
		END IF
		ii = ii - 1
		IF ( ii .lt. 1 ) out = .false.
	    END DO
	END IF
C
C*	Check for systems west of the dateline and do not plot the
C*	danger area.
C
	IF ( ocean .ne. 'AL' ) THEN
	    ii = 1
	    out = .true.
	    DO WHILE ( out )
		IF ( tlon ( ii ) .gt. 0. ) THEN
		    DO jj = ii, ITRACK
			tlon ( jj ) = RMISSD
		    END DO
		    out = .false.
		  ELSE
		    ii = ii + 1
		    IF ( ii .gt. ITRACK ) out = .false.
		END IF
	    END DO
	END IF
C
C*	Get the quadrant radii for the danger area.
C
	out = .false.
	DO ii = 1, ITRACK
	    CALL ST_C2I ( f34kt ( ii ), 5, intarr, num, ier )
	    IF ( ( ier .lt. 0 ) .or. ( num .lt. 5 ) .or.
     +		 ERMISS ( tlat (ii) ) .or. ERMISS ( tlon (ii) ) ) THEN
		good ( ii ) = .false.
	      ELSE
		good ( ii ) = .true.
		DO jj = 1, 4
		    IF ( intarr ( jj + 1 ) .ge. 0 ) THEN
		        quad ( jj, ii ) = FLOAT ( intarr ( jj+1 ) ) +
     +					  radii ( ii )
		      ELSE
			quad ( jj, ii ) = radii ( ii )
			IF ( ii .eq. 2 ) out = .true.
		    END IF
		    IF ( ii .eq. 1 ) THEN
			IF ( quad ( jj, ii ) .lt. RADCUR )
     +			     quad ( jj, ii ) = RADCUR
		    END IF
		END DO
		IF ( ( ii .gt. 2 ) .and. out ) good ( 2 ) = .false.
	    END IF
	END DO
C
C*	Loop to get the bounds for each track point.
C*      This section of code was adapted from routine GH_SWLN.
C
       	DO ii = 1, ITRACK
            IF ( good ( ii ) ) THEN
		dist ( 1 ) = quad ( 2, ii )
		dist ( 2 ) = quad ( 1, ii )
		dist ( 3 ) = quad ( 4, ii )
		dist ( 4 ) = quad ( 3, ii )
		dist ( 5 ) = dist ( 1 )
C
C*	 	Loop through each quadrant of a track point.
C
                nn = 0
                DO jj = 1, 4
                    DO iangle = IANG1, IANG2
                        nn = nn + 1
			theta = FLOAT ( iangle )
			rtheta = DTR * theta
                        rad = COS (rtheta) * COS (rtheta) * dist (jj)
     +                      + SIN (rtheta) * SIN (rtheta) * dist (jj+1)
C
			rtheta = DTR * ( theta + 270.0 + angle ( jj ) )
                        y1 = rad * SIN ( rtheta ) / 60.0
                        x1 = (-rad) * COS ( rtheta ) / 60.0
     +                       / COS ( DTR * tlat ( ii ) )
C
                        xlat ( nn, ii ) = tlat ( ii ) + y1
                        xlon ( nn, ii ) = tlon ( ii ) - x1
                    END DO
                END DO
	    END IF
	END DO
C
C*      Determine if the bounds for one track point lie completely
C*	inside the bounds for another.  If so, the inner bounds can be
C*	eliminated.
C
        sys = 'M'
        CALL ST_NULL ( sys, sys, lens, ier )
C
	DO ii = 1, ITRACK - 1
	    IF ( good ( ii ) ) THEN
		DO jj = ii + 1, ITRACK
		    IF ( good ( jj ) ) THEN
                        CALL CGR_INPOLY ( sys, IPOLY, xlat (1,ii),
     +				   xlon (1,ii), sys, IPOLY, xlat (1,jj),
     +				   xlon (1,jj), ino12, ier )
        		CALL GG_TCEP ( xlat (1,ii), xlon (1,ii), IPOLY,
     +				       ino12, ipt1, nseg1, ier )
        		IF ( nseg1 .eq. 0 ) THEN
            		    good ( ii ) = .false.
          		  ELSE
            		    CALL CGR_INPOLY ( sys, IPOLY, xlat (1,jj),
     +			           xlon (1,jj), sys, IPOLY, xlat (1,ii),
     +                             xlon (1,ii), ino21, ier )
            		    CALL GG_TCEP ( xlat (1,jj), xlon (1,jj),
     +			                IPOLY, ino21, ipt2, nseg2, ier )
            		    IF ( nseg2 .eq. 0 ) THEN
                		good ( jj ) = .false.
            		    END IF
        		END IF
		    END IF
		END DO
	    END IF
	END DO
C
	CALL GQDEV ( ddev, junit, jatyp, ier )
	IF ( ddev .ne. 'PS' ) THEN
	    ilwid = 4
	  ELSE
	    ilwid = 8
	END IF
 	CALL GSLINE ( 0, 0, ilwid, 0, ier )
C
C*	Get all the "good" track points consecutively for ease of
C*	interpolation and plotting.
C
	igood = 0
	DO ii = 1, ITRACK
	    IF ( good ( ii ) ) THEN
		igood = igood + 1
		tlat ( igood ) = tlat ( ii )
		tlon ( igood ) = tlon ( ii )
		DO jj = 1, 4
		    quad ( jj, igood ) = quad ( jj, ii )
		END DO
		DO jj = 1, IPOLY
		    xlat ( jj, igood ) = xlat ( jj, ii )
		    xlon ( jj, igood ) = xlon ( jj, ii )
		END DO
	    END IF
	END DO
C
	IF ( igood .eq. 0 ) THEN
	    RETURN
	END IF
C
C*	Interpolate the radii between the polygons for a smooth
C*	boundary.
C
	ipt = 1
	DO ip = 1, igood
	    dlat ( ipt ) = tlat ( ip )
	    dlon ( ipt ) = tlon ( ip )
	    dse ( ipt )  = quad ( 2, ip )
	    dne ( ipt )  = quad ( 1, ip )
	    dnw ( ipt )  = quad ( 4, ip )
	    dsw ( ipt )  = quad ( 3, ip )
	    iptsv ( ip ) = ipt
C
	    IF ( ip .lt. igood ) THEN
		CALL CLO_DIST ( tlat (ip), tlon (ip), 1, tlat (ip+1),
     +				tlon (ip+1), dd, ier )
		IF ( dd .gt. 5000. ) THEN
		    radmin = AMIN1 ( quad (1,ip), quad (2,ip),
     +				     quad (3,ip), quad (4,ip),
     +				     quad (1,ip+1), quad (2,ip+1),
     +				     quad (3,ip+1), quad (4,ip+1) )
		    radavg = NINT ( ( quad (1,ip) + quad (2,ip) +
     +				      quad (3,ip) + quad (4,ip) +
     +				      quad (1,ip+1) + quad (2,ip+1) +
     + 				      quad (3,ip+1) + quad (4,ip+1) )
     +				      / 8. )
		    space = PR_HGNM ( radavg )
                    intrp = NINT ( dd / space + .6 )
                    incr  = NINT ( radavg - radmin ) / 9
                    intrp = intrp + incr
		    IF ( ( ipt + intrp ) .gt. MXINTP )
     +				 intrp = MXINTP - ipt
                    intvl = intrp + 1
                    delta = dd / FLOAT ( intvl )
                    CALL CLO_DIRECT ( tlat (ip+1), tlon (ip+1),
     +                             tlat ( ip ), tlon ( ip ), ddir, ier )
                    ddist = delta
                    DO ii = 1, intrp
                        ipt = ipt + 1
                        CALL CLO_DLTLN ( tlat (ip), tlon (ip), ddist,
     +                             ddir, dlat (ipt), dlon (ipt), ier )
                        wt2 = FLOAT ( ii ) / FLOAT ( intvl )
                        wt1 = 1. - wt2
                        dse ( ipt ) = quad ( 2,ip ) * wt1 +
     +                                quad ( 2, ip+1 ) * wt2
                        dne ( ipt ) = quad ( 1,ip ) * wt1 +
     +                                quad ( 1, ip+1 ) * wt2
                        dnw ( ipt ) = quad ( 4,ip ) * wt1 +
     +                                quad ( 4, ip+1 ) * wt2
                        dsw ( ipt ) = quad ( 3,ip ) * wt1 +
     +                                quad ( 3, ip+1 ) * wt2
                        ddist = ddist + delta
                    END DO
                END IF
                ipt = ipt + 1
		IF ( ipt .gt. MXINTP ) THEN
		    ipt = MXINTP
		END IF
            END IF
        END DO
C
C*      Calculate polygon bounds along the interpolated track.
C
	jj = 1
        DO ip = 1, ipt
C
	    IF ( ip .eq. iptsv ( jj ) ) THEN
C
C*		Get the bounds for the track point.
C
		DO ii = 1, IPOLY
		    blat ( ii ) = xlat ( ii, jj )
		    blon ( ii ) = xlon ( ii, jj )
		END DO
		jj = jj + 1
	      ELSE
C
C*		Calculate the bounds for an interpolated track point.
C
                dist ( 1 ) = dse ( ip )
                dist ( 2 ) = dne ( ip )
                dist ( 3 ) = dnw ( ip )
                dist ( 4 ) = dsw ( ip )
                dist ( 5 ) = dist (1)
C
C*	        Loop through each quadrant of the polygon.  The radii
C*              are currently calculated to generate a smooth contour
C*              when sweeping from one quadrant to the next in 1 degree
C*	        azimuthal increments.
C
                nn = 0
                DO ii = 1, 4
                    DO iangle  = IANG1, IANG2
                        nn     = nn + 1
		        theta  = FLOAT ( iangle )
		        rtheta = DTR * theta
                        rad    = COS (rtheta) * COS (rtheta) * dist (ii)
     +                     + SIN (rtheta) * SIN (rtheta) * dist (ii + 1)
C
		        rtheta = DTR * ( theta + 270.0 + angle(ii) )
                        y1     = rad * SIN ( rtheta ) / 60.0
                        x1     = (-rad) * COS ( rtheta ) / 60.0
     +                           / COS ( DTR * dlat ( ip ) )
C
                        blat ( nn ) = dlat ( ip ) + y1
                        blon ( nn ) = dlon ( ip ) - x1
                    END DO
                END DO
	    END IF
C
	    IF ( ip .eq. 1 ) THEN
C
C*		Save the first polygon.
C
		DO ii = 1, IPOLY
		    wlat ( ii ) = blat ( ii )
		    wlon ( ii ) = blon ( ii )
		END DO
		npt = IPOLY
	      ELSE
C
C*	        Combine the new polygon, defined by (blat, blon, IPOLY),
C*	        with the preceding one (wlat, wlon, npt).  The points
C*	        for the combined polygons are defined by (olat, olon,
C*		npts).
C
	        CALL GG_TCSH ( wlat, wlon, npt, blat, blon, IPOLY,
     +			       olat, olon, npts, ier )
		IF ( ier .eq. 0 ) THEN
C
C*		    The polygons intersect.  Move the combined
C*		    polygons to the first polygon array.
C
		    iskip = 1
		    IF ( npts .ge. ( MXWORK - IPOLY ) ) iskip = 2
		    nn = 0
		    DO ii = 1, npts, iskip
			nn = nn + 1
			wlat ( nn ) = olat ( ii )
			wlon ( nn ) = olon ( ii )
		    END DO
		    npt = nn
	    	    wlat ( npt ) = wlat ( 1 )
	    	    wlon ( npt ) = wlon ( 1 )
		  ELSE
C
C*		    The polygons do not intersect.
C
		    CALL ER_WMSG ( 'GG', 4, ' ', ier )
		END IF
	    END IF
	END DO
C
	IF ( npt .gt. MXPTS ) THEN
C
C*	    Thin the points for vgf.
C
	    iskip = ( npt - 1 ) / MXPTS + 1
	    nn    = 0
	    DO ii = 1, npt, iskip
	        nn = nn + 1
		wlat ( nn ) = wlat ( ii )
		wlon ( nn ) = wlon ( ii )
	    END DO
	    IF ( nn .gt. MXPTS ) nn = MXPTS
	    npt = nn
	    wlat ( npt ) = wlat ( 1 )
	    wlon ( npt ) = wlon ( 1 )
	END IF
C
	nsmth = 50
	iptsm = 5
	CALL GG_TCSM ( nsmth, npt, iptsm, wlat, wlon, ier )
	CALL GFILL ( 'M', npt, wlat, wlon, ier )
	IF ( ddev .ne. 'VG' ) THEN
C
C*	    The above GFILL call generates bounding lines for the
C*	    polygon for vgf, but not for other devices.
C
   	    CALL GLINE ( 'M', npt, wlat, wlon, ier )
	END IF
	CALL GSLINE ( jltyp, jlthw, jwidth, jlwhw, ier )
C
C*	Get the bottommost, rightmost, topmost and leftmost points.
C
	elat ( 1 ) = wlat ( 1 )
	elat ( 3 ) = wlat ( 1 )
	elon ( 2 ) = wlon ( 1 )
	IF ( elon ( 2 ) .gt. 0. ) elon ( 2 ) = -180.
	elon ( 4 ) = wlon ( 1 )
	i1 = 1
	i2 = 1
	i3 = 1
	i4 = 1
	DO ii = 2, npt
	    IF ( wlat ( ii ) .lt. elat ( 1 ) ) THEN
		elat ( 1 ) = wlat ( ii )
		i1 = ii
	    END IF
	    IF ( wlat ( ii ) .gt. elat ( 3 ) ) THEN
		elat ( 3 ) = wlat ( ii )
		i3 = ii
	    END IF
	    IF ( ( wlon ( ii ) .gt. elon ( 2 ) ) .and.
     +		 ( wlon ( ii ) .le. 0. ) ) THEN
		elon ( 2 ) = wlon ( ii )
		i2 = ii
	    END IF
	    IF ( wlon ( ii ) .lt. elon ( 4 ) ) THEN
		elon ( 4 ) = wlon ( ii )
		i4 = ii
	    END IF
	END DO
	elat ( 2 ) = wlat ( i2 )
	elat ( 4 ) = wlat ( i4 )
	elon ( 1 ) = wlon ( i1 )
	elon ( 3 ) = wlon ( i3 )
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
C     -*- fortran -*-
C     This file is autogenerated with f2py (version:1.21.1)
C     It contains Fortran 77 wrappers to fortran functions.

      subroutine f2pywraperf_c (erf_cf2pywrap, xx)
      external erf_c
      real xx
      real erf_cf2pywrap, erf_c
      erf_cf2pywrap = erf_c(xx)
      end

      subroutine f2pywrapfermid (fermidf2pywrap, xx)
      external fermid
      real xx
      real fermidf2pywrap, fermid
      fermidf2pywrap = fermid(xx)
      end

      subroutine f2pywrapdelthm (delthmf2pywrap, xx)
      external delthm
      real xx
      real delthmf2pywrap, delthm
      delthmf2pywrap = delthm(xx)
      end

      subroutine f2pywrapspline (splinef2pywrap, x)
      external spline
      real x
      real splinef2pywrap, spline
      splinef2pywrap = spline(x)
      end

      subroutine f2pywrapposhm (poshmf2pywrap, x)
      external poshm
      real x
      real poshmf2pywrap, poshm
      poshmf2pywrap = poshm(x)
      end

      subroutine f2pywrapposhm2 (poshm2f2pywrap, x)
      external poshm2
      real x
      real poshm2f2pywrap, poshm2
      poshm2f2pywrap = poshm2(x)
      end

      subroutine f2pyinitwarn(setupfunc)
      external setupfunc
      integer iwarn
      common /warn/ iwarn
      call setupfunc(iwarn)
      end

max_stars_repo_name: aiidateam/aiida-sssp-workflow

max_issues_repo_name: aiidateam/aiida-sssp-workflow

max_forks_repo_name: aiidateam/aiida-sssp-workflow


text:
        SUBROUTINE GDGHIR ( inttyp, glat, glon, np, maxx, maxy, grid,
     +			    grdo, iret )
C************************************************************************
C* GDGHIR								*
C*									*
C* This subroutine interplates the data in the input grid to the	*
C* output grid whose positions relative to the input grid are given	*
C* in GLAT and GLON as I and J positions.				*
C*									*
C* GRID and GRDO may be the same array.					*
C*									*
C* GDGHIR ( INTTYP, GLAT, GLON, NP, MAXX, MAXY, GRID, GRDO, IRET )	*
C*									*
C* Input parameters:							*
C*	INTTYP		INTEGER		Interpolation type		*
C*					  0 = linear interpolation	*
C*					  1 = nearest point assignment	*
C*	GLAT (NP)	REAL		GRID relative I positions	*
C*	GLON (NP)	REAL		GRID relative J positions	*
C*	NP		INTEGER		Number of output points		*
C*	MAXX		INTEGER		Maximum I dimension of GRID	*
C*	MAXY		INTEGER		Maximum J dimension of GRID	*
C*	GRID (MAXX,MAXY)REAL		Grid of data			*
C*									*
C* Output parameters:							*
C*	GRDO (*)	REAL		Output grid of interplated data *
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					-18 = invalied interp type	*
C**									*
C* Log:									*
C* K. Brill/HPC		 2/00						*
C* T. Lee/SAIC		 7/03		Changed LLMXGD -> LLMXTG	*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	REAL 		grid (*), grdo (*), glat (*), glon (*)
C*
	REAL		gout (LLMXTG)
C*
	INCLUDE		'ERMISS.FNC'
	INDXQ ( lx, ly, kx ) = ( ly - 1 ) * kx + lx
C------------------------------------------------------------------------
	iret = 0
C
C*	Interpolate to output grid at valid points.
C
	IF ( inttyp .eq. 0 ) THEN
C
C*	    Bilinear interpolation.  Note:  GR_INTP is not used
C*	    because it does no check for missing coordinates.
C
	    DO ii = 1, np
		IF ( .not.
     +		  ( ERMISS ( glat (ii) ) .or. ERMISS ( glon (ii) )
     +		 .or. glat (ii) .lt. 1 .or. glon (ii) .lt. 1
     +		 .or. glat (ii) .gt. maxx .or. glon (ii) .gt. maxy ) )
     +		THEN
		    icol = glat (ii)
		    irow = glon (ii)
		    IF ( icol .eq. maxx ) icol = icol - 1
		    IF ( irow .eq. maxy ) irow = irow - 1
		    icolp1 = icol + 1
		    irowp1 = irow + 1
		    c = glat (ii) - FLOAT (icol)
		    r = glon (ii) - FLOAT (irow)
		    omc = 1. - c
		    omr = 1. - r
		    IF ( ERMISS (grid(INDXQ(icol,irow,maxx)))
     +			 .or.
     +			 ERMISS (grid(INDXQ(icolp1,irow,maxx)))
     +			 .or.
     +			 ERMISS (grid(INDXQ(icol,irowp1,maxx)))
     +			 .or.
     +			 ERMISS (grid(INDXQ(icolp1,irowp1,maxx))) )
     +		    THEN
			gout (ii) = RMISSD
		    ELSE
			gout (ii) =
     +			( grid (INDXQ(icol,irow,maxx)) * omc +
     +			  grid (INDXQ(icolp1,irow,maxx)) * c ) * omr +
     +		        ( grid (INDXQ(icol,irowp1,maxx)) * omc +
     +		          grid (INDXQ(icolp1,irowp1,maxx) ) * c ) * r
		    END IF
		ELSE
		    gout (ii) = RMISSD
		END IF
	    END DO
	ELSE IF ( inttyp .eq. 1 ) THEN
C
C*	    Nearest value interpolation (NVI).
C
	    DO ii = 1, np
		IF ( .not.
     +		  ( ERMISS ( glat (ii) ) .or. ERMISS ( glon (ii) )
     +		 .or. glat (ii) .lt. 1 .or. glon (ii) .lt. 1
     +		 .or. glat (ii) .gt. maxx .or. glon (ii) .gt. maxy ) )
     +		THEN
		    nrx = NINT ( glat (ii) )
		    nry = NINT ( glon (ii) )
		    gout (ii) = grid (INDXQ(nrx,nry,maxx))
		ELSE
		    gout (ii) = RMISSD
		END IF
	    END DO
	ELSE
	    iret = -18
	    RETURN
	END IF
	DO ii = 1, np
	    grdo (ii) = gout (ii)
	END DO
C*
	RETURN
	END