text:
C     ------------------------------------------------------------------
      SUBROUTINE EXANGLES (EANGLES, WHICH)
C     ------------------------------------------------------------------
C     Calls EXANGLES2, which does the actual energy calculation
C
C     By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003)
C     ------------------------------------------------------------------
      IMPLICIT NONE

      INCLUDE  "cns.inc"
      INCLUDE  "heap.inc"
      INCLUDE  "xangle.inc"

      DOUBLE PRECISION EANGLES
      CHARACTER*7 WHICH

      CALL EXANGLES2(EANGLES,HEAP(XANGLESJPTR),HEAP(XANGLESKPTR),
     &                       HEAP(XANGLESJLST),HEAP(XANGLESKLST),
     &                       HEAP(XANGLESOBS1PTR),HEAP(XANGLESOBS2PTR),
     &                       HEAP(XANGLESERRPTR),
     &                       HEAP(CALCXANGLESPTR),WHICH)

      RETURN
      END
C     ------------------------------------------------------------------
      SUBROUTINE EXANGLES2 (EANGLES,ATOMJPTR,ATOMKPTR,ATOMJLST,ATOMKLST,
     &                      XANGLESOBS1,XANGLESOBS2,XANGLESERR,
     &                      XANGLESCALC,WHICH)
C     ------------------------------------------------------------------
C     Calculates RDC-ANGLES energies
C
C     Energies are of the form:
C        E = K1*(SIN(ANG)**2)
C     where:
C        K1 = Force constant,
C        ANG = The angle between the target vector and the actual vector
C
C     Note that:
C        Atom J = Detected nucleus (e.g. H)
C        Atom K = Coupled nucleus (e.g. N)
C
C     By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003)
C     ------------------------------------------------------------------
      IMPLICIT NONE

      INCLUDE  "cns.inc"
      INCLUDE  "mtf.inc"
      INCLUDE  "coord.inc"
      INCLUDE  "numbers.inc"
      INCLUDE  "deriv.inc"
      INCLUDE  "xangle.inc"
      INCLUDE  "consta.inc"

      INTEGER ATOMJPTR(*),ATOMKPTR(*)
      INTEGER ATOMJLST(*),ATOMKLST(*)
      DOUBLE PRECISION XANGLESOBS1(*),XANGLESOBS2(*),XANGLESERR(*)
      DOUBLE PRECISION XANGLESCALC(*),EANGLES
      CHARACTER*7 WHICH
C     ------------------------------------------------------------------
C     Local variables
C     ------------------------------------------------------------------
      INTEGER COUNT,CLASS,MM,NN,II,I,J,K,L,NA,M,N
      PARAMETER (NA=20)
      INTEGER DMM,DNN,OUTFLAG(NA),MCOUNT,NCOUNT,PCOUNT
      DOUBLE PRECISION XJ,XK,YJ,YK,ZJ,ZK,CALCXANGLES,E,
     &                 DFXJ(NA),DFYJ(NA),DFZJ(NA),
     &                 DFXK(NA),DFYK(NA),DFZK(NA),DF1,
     &                 O1,O2,O3,O4,O5,O6,O7,O8,O9,O10,O11,O12,
     &                 PS,DJK,DJK2,DJK3
      DOUBLE PRECISION OBSXANGLES1,OBSXANGLES2,ERRXANGLES,K1,
     &                 A,B,
     &                 DT,DP,DELTA,DELTAXANGLES,
     &                 DD
C     ------------------------------------------------------------------
C     Zero out partial energy
C     ------------------------------------------------------------------
      EANGLES = ZERO

      CLASS = 1
      K1=XANGLESFORCES(1)

      COUNT=0

      DO WHILE(COUNT.LT.ANGLESNUM)
         COUNT = COUNT+1
C     ------------------------------------------------------------------
C     Reset individual E to zero
C     ------------------------------------------------------------------
         E=0

         DO WHILE ((XANGLESASSNDX(CLASS).LT.COUNT).OR.
     &             (XANGLESASSNDX(CLASS).EQ.0))
                    CLASS=CLASS+1
         END DO

         IF (XANGLESASSNDX(CLASS).EQ.XANGLESASSNDX(CLASS-1)) THEN
             COUNT=COUNT-1
         ENDIF

         K1=XANGLESFORCES(CLASS)
C     ------------------------------------------------------------------
C     Note there should only be one atom for J and K
C     ------------------------------------------------------------------
         J=ATOMJPTR(COUNT)+1
         K=ATOMKPTR(COUNT)+1

         XJ=X(ATOMJLST(J))
         XK=X(ATOMKLST(K))
         YJ=Y(ATOMJLST(J))
         YK=Y(ATOMKLST(K))
         ZJ=Z(ATOMJLST(J))
         ZK=Z(ATOMKLST(K))

         OBSXANGLES1=XANGLESOBS1(COUNT)
         OBSXANGLES2=XANGLESOBS2(COUNT)
         ERRXANGLES=XANGLESERR(COUNT)
C     ------------------------------------------------------------------
C     Initialize calculated RDC-angles and counter
C     ------------------------------------------------------------------
         CALCXANGLES=0
         II=0
         MCOUNT=0
         NCOUNT=0
C     ------------------------------------------------------------------
C     Check for correct permutations of paired atoms to get the nucleus-
C     nucleus vector. This depends solely on the order of the assigned
C     atoms. OUTFLAG=1 indicates that the permutation is not allowed.
C     ------------------------------------------------------------------
         DMM=ATOMJPTR(COUNT+1)-ATOMJPTR(COUNT)
         DNN=ATOMKPTR(COUNT+1)-ATOMKPTR(COUNT)
         DO MM=ATOMJPTR(COUNT)+1,ATOMJPTR(COUNT+1)
            MCOUNT=MCOUNT+1
            NCOUNT=0
            DO NN=ATOMKPTR(COUNT)+1,ATOMKPTR(COUNT+1)
               NCOUNT=NCOUNT+1
               II=II+1
               IF ((DMM.GT.1).AND.(DNN.GT.1)) THEN
                    IF (MCOUNT.EQ.NCOUNT) THEN
                        OUTFLAG(II)=0
                    ELSE
                        OUTFLAG(II)=1
                    ENDIF
               ELSE
                    OUTFLAG(II)=0
               END IF
            END DO
         END DO

         II=0
         PCOUNT=0
         DO MM=ATOMJPTR(COUNT)+1,ATOMJPTR(COUNT+1)
            DO NN=ATOMKPTR(COUNT)+1,ATOMKPTR(COUNT+1)
               II=II+1
               IF (OUTFLAG(II).NE.1) THEN
                   PCOUNT=PCOUNT+1
C     ------------------------------------------------------------------
C     Calculate the function and its derivatives
C     ------------------------------------------------------------------
                   OBSXANGLES1=OBSXANGLES1*PI/180.
                   IF (OBSXANGLES1.LT.0.) OBSXANGLES1=OBSXANGLES1+2.*PI
                   OBSXANGLES2=OBSXANGLES2*PI/180.
                   IF (OBSXANGLES2.LT.0.) OBSXANGLES2=OBSXANGLES2+2.*PI

                   O1=SIN(OBSXANGLES1)*COS(OBSXANGLES2)
                   O2=SIN(OBSXANGLES1)*SIN(OBSXANGLES2)
                   O3=COS(OBSXANGLES1)
                   O4=XJ-XK
                   O5=YJ-YK
                   O6=ZJ-ZK
                   O7=O4**2
                   O8=O5**2
                   O9=O6**2
                   DJK2=O7+O8+O9
                   DJK=SQRT(DJK2)
                   DJK3=DJK**3
                   PS=O1*O4+O2*O5+O3*O6
                   O10=O1*DJK-PS*O4/DJK
                   O11=O2*DJK-PS*O5/DJK
                   O12=O3*DJK-PS*O6/DJK

                   DFXJ(II)=-2.*PS*O10/DJK3
                   DFYJ(II)=-2.*PS*O11/DJK3
                   DFZJ(II)=-2.*PS*O12/DJK3
                   DFXK(II)=-DFXJ(II)
                   DFYK(II)=-DFYJ(II)
                   DFZK(II)=-DFZJ(II)
C     ------------------------------------------------------------------
C     Calculate the cosine square of the angle
C     ------------------------------------------------------------------
                   CALCXANGLES=(PS/DJK)**2
               END IF
            END DO
         END DO

         IF (WHICH.EQ.'ANALYZE') THEN
             XANGLESCALC(COUNT)=CALCXANGLES
         END IF
C     ------------------------------------------------------------------
C     The deviation is the sine square of the angle
C     ------------------------------------------------------------------
         DELTA=(1.-CALCXANGLES)
C     ------------------------------------------------------------------
C     Adjust the deviation based on the error range (expressed as the
C     accepted deviation of the cosine square from 1)
C     ------------------------------------------------------------------
         IF ((DELTA.GT.0.000).AND.(DELTA.GT.ERRXANGLES)) THEN
              DELTAXANGLES=DELTA-ERRXANGLES
         ELSE
              DELTAXANGLES=0.0
         END IF

         DD=DELTAXANGLES

         II=0
         PCOUNT=0
         DO MM=ATOMJPTR(COUNT)+1,ATOMJPTR(COUNT+1)
            DO NN=ATOMKPTR(COUNT)+1,ATOMKPTR(COUNT+1)
               II=II+1
               IF (OUTFLAG(II).NE.1) THEN
                   PCOUNT=PCOUNT+1
                   IF (DD.EQ.0.0) THEN
                       E=E+0.0
                       DF1=0.0
                   ELSE
                       E=E+K1*DELTAXANGLES
                       DF1=K1
                   END IF
               END IF
            END DO
         END DO
C     ------------------------------------------------------------------
C     Accumulate energy
C     ------------------------------------------------------------------
         EANGLES=EANGLES+E
C     ------------------------------------------------------------------
C     Now update forces if in energy/force mode
C     ------------------------------------------------------------------
         IF (WHICH.NE.'ANALYZE') THEN
             II=0
             DO MM=ATOMJPTR(COUNT)+1,ATOMJPTR(COUNT+1)
                DO NN=ATOMKPTR(COUNT)+1,ATOMKPTR(COUNT+1)
                   II=II+1
                   IF (OUTFLAG(II).NE.1) THEN
                       DX(ATOMJLST(J))=DX(ATOMJLST(J))+DF1*DFXJ(II)
                       DY(ATOMJLST(J))=DY(ATOMJLST(J))+DF1*DFYJ(II)
                       DZ(ATOMJLST(J))=DZ(ATOMJLST(J))+DF1*DFZJ(II)
                       DX(ATOMKLST(K))=DX(ATOMKLST(K))+DF1*DFXK(II)
                       DY(ATOMKLST(K))=DY(ATOMKLST(K))+DF1*DFYK(II)
                       DZ(ATOMKLST(K))=DZ(ATOMKLST(K))+DF1*DFZK(II)
                   END IF
                END DO
             END DO
         END IF
      END DO

      RETURN
      END
C     ------------------------------------------------------------------
      SUBROUTINE READXANGLES
C     ------------------------------------------------------------------
C     Reads in RDC-ANGLES information
C
C     By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003)
C     ------------------------------------------------------------------
      IMPLICIT NONE

      INCLUDE  "cns.inc"
      INCLUDE  "comand.inc"
      INCLUDE  "xangle.inc"
      INCLUDE  "funct.inc"
      INCLUDE  "mtf.inc"
      INCLUDE  "heap.inc"
      INCLUDE  "numbers.inc"

      INTEGER COUNT,SPTR,OLDCLASS,OLDMAXXANGLES
      DOUBLE PRECISION K1,CUTOFF
      CHARACTER*4 THENAME

      SPTR=ALLHP(INTEG4(NATOM))
      CALL PUSEND('XANGLE>')
862   CONTINUE
      CALL NEXTWD('XANGLE>')
      CALL MISCOM('XANGLE>',USED)
      IF (.NOT.USED) THEN
C     ------------------------------------------------------------------
C     Documentation
C     ------------------------------------------------------------------
          IF (WD(1:4).EQ.'HELP') THEN
              WRITE(DUNIT,'(10X,A)')
     &              ' XANGLE {<ANGLE-STATEMENT>} END ',
     &              ' <ANGLE-STATEMENT>:== ',
     &              ' ASSIGN <SEL> <SEL> <REAL> <REAL> <REAL>',
     &              ' * Restraint: Detected_nucleus Coupled_nucleus',
     &              ' Theta Phi Err *',
     &              ' CLASsification <NAME>',
     &              ' * Starts a new class *',
     &              ' FORCeconstant <REAL>',
     &              ' * Force constant for the current class *',
     &              ' NREStraints <INTEGER>',
     &              ' * Number of slots for restraints to allocate *',
     &              ' PRINt THREshold <REAL>',
     &              ' * Prints violations larger than the THREshold *',
     &              ' RESEt',
     &              ' * Erases the restraint table, but keeps NRES *'
C     ------------------------------------------------------------------
C     Get class name. Determine if it's an already-defined class.
C     Insert a new class if it's not.
C     ------------------------------------------------------------------
          ELSE IF (WD(1:4).EQ.'CLAS') THEN
              OLDCLASS=CURCLASS
              CALL NEXTA4('Class name =',THENAME)
              MODE=NEW
              DO COUNT=1,NCLASSES
                 IF (XANGLESCLASSNAMES(COUNT).EQ.THENAME) THEN
                     MODE=UPDATE
                     CURCLASS=COUNT
                  END IF
              END DO
              IF (MODE.EQ.NEW) THEN
C     ------------------------------------------------------------------
C     Make sure you can't add more than the maximum number of classes
C     ------------------------------------------------------------------
                  IF (OLDCLASS.EQ.MAXXANGLESCLASSES) THEN
                      CALL DSPERR('XANGLE','Too many classes.')
                      CALL DSPERR('XANGLE',
     &                     'Increase MAXXANGLESCLASSES and recompile.')
                      CALL WRNDIE(-5, 'READXANGLES',
     &                     'Too many RDC-angles classes.')
                   END IF
                   NCLASSES=NCLASSES+1
                   CURCLASS=NCLASSES
                   XANGLESCLASSNAMES(CURCLASS)=THENAME
C     ------------------------------------------------------------------
C     If this isn't the first class, close off the old class
C     ------------------------------------------------------------------
                   IF (NCLASSES.GT.1) THEN
                       XANGLESASSNDX(OLDCLASS)=ANGLESNUM
                   END IF
              END IF
C     ------------------------------------------------------------------
C     Set force constant for current class
C     ------------------------------------------------------------------
          ELSE IF (WD(1:4).EQ.'FORC') THEN
C     ------------------------------------------------------------------
C     Start a default class if there isn't one defined
C     ------------------------------------------------------------------
              IF (CURCLASS.EQ.0) THEN
                  NCLASSES = 1
                  CURCLASS = 1
              END IF
              CALL NEXTF('Force constant =',K1)
              WRITE(DUNIT,'(A,A,A,F8.3)')
     &              'Setting force constant for class ',
     &              XANGLESCLASSNAMES(CURCLASS),' to ',K1
                    XANGLESFORCES(CURCLASS)=K1
C     ------------------------------------------------------------------
C     Reset RDC-angles database
C     ------------------------------------------------------------------
          ELSE IF (WD(1:4).EQ.'RESE') THEN
               CALL XANGLESDEFAULTS
               CALL ALLOCXANGLES(0,MAXXANGLES)
C     ------------------------------------------------------------------
C     Change number of assignment slots
C     ------------------------------------------------------------------
          ELSE IF (WD(1:4).EQ.'NRES') THEN
               OLDMAXXANGLES=MAXXANGLES
               CALL NEXTI('Number of slots =',MAXXANGLES)
               CALL ALLOCXANGLES(OLDMAXXANGLES,MAXXANGLES)
               WRITE(DUNIT,'(A,I8,A)')
     &            'XANGLE: Allocating space for',MAXXANGLES,
     &            'number of restraints.'
C     ------------------------------------------------------------------
C     Read in an assignment
C     ------------------------------------------------------------------
          ELSE IF (WD(1:4).EQ.'ASSI') THEN
C     ------------------------------------------------------------------
C     Make sure you can't add more assignments than you have slots for
C     ------------------------------------------------------------------
               IF (XANGLESMX.EQ.MAXXANGLES) THEN
                   CALL DSPERR('XANGLE','Too many assignments.')
                   CALL DSPERR('XANGLE',
     &                         'Increasing NREStraints by 100.')
                   OLDMAXXANGLES=MAXXANGLES
                   MAXXANGLES=MAXXANGLES+100
                   CALL ALLOCXANGLES(OLDMAXXANGLES,MAXXANGLES)
               END IF
C     ------------------------------------------------------------------
C     If there isn't a class specified, start a default class
C     ------------------------------------------------------------------
               IF (CURCLASS.EQ.0) THEN
                   NCLASSES=1
                   CURCLASS=1
               END IF
               CALL READXANGLES2(HEAP(XANGLESJPTR),HEAP(XANGLESKPTR),
     &                           HEAP(XANGLESJLST),HEAP(XANGLESKLST),
     &                           HEAP(XANGLESOBS1PTR),
     &                           HEAP(XANGLESOBS2PTR),
     &                           HEAP(XANGLESERRPTR),HEAP(SPTR))
C     ------------------------------------------------------------------
C     Print violations
C     ------------------------------------------------------------------
          ELSE IF (WD(1:4).EQ.'PRIN') THEN
               CALL NEXTWD('PRINT>')
               IF (WD(1:4).NE.'THRE') THEN
                   CALL DSPERR('XANGLE',
     &                  'PRINt expects THREshold parameter.')
               ELSE
                   CALL NEXTF('THREshold =',CUTOFF)
                   IF (CUTOFF.LT.ZERO) THEN
                       CALL DSPERR('XANGLE',
     &                             'Cutoff must be positive.')
                       CUTOFF = ABS(CUTOFF)
                   END IF
                   CALL NEXTA4('ALL OR CLASs>',THENAME)
                   IF (THENAME(1:3).EQ.'ALL') THEN
                       DO COUNT=1,NCLASSES
                          PRINTCLASS(COUNT)=.TRUE.
                       END DO
                   ELSE IF (THENAME(1:4).EQ.'CLAS') THEN
                       CALL NEXTA4('CLASS NAME =',THENAME)
                       DO COUNT=1,NCLASSES
                          IF (XANGLESCLASSNAMES(COUNT).EQ.THENAME) THEN
                              PRINTCLASS(COUNT)=.TRUE.
                          ELSE
                              PRINTCLASS(COUNT)=.FALSE.
                          END IF
                       END DO
                    ELSE
                       DO COUNT=1,NCLASSES
                          PRINTCLASS(COUNT)=.TRUE.
                       END DO
                    END IF
                    CALL PRINTXANGLES(CUTOFF,HEAP(CALCXANGLESPTR),
     &                                HEAP(XANGLESOBS1PTR),
     &                                HEAP(XANGLESOBS2PTR),
     &                                HEAP(XANGLESERRPTR),
     &                                HEAP(XANGLESJPTR),
     &                                HEAP(XANGLESKPTR),
     &                                HEAP(XANGLESJLST),
     &                                HEAP(XANGLESKLST))
               END IF
C     ------------------------------------------------------------------
C     Check for END statement
C     ------------------------------------------------------------------
          ELSE
              CALL CHKEND('XANGLE>', DONE)
          END IF
      END IF
      IF (.NOT.(DONE)) GOTO 862
      DONE=.FALSE.
      CALL FREHP(SPTR,INTEG4(NATOM))

      RETURN
      END
C     ------------------------------------------------------------------
      SUBROUTINE ALLOCXANGLES (OLDSIZE, NEWSIZE)
C     ------------------------------------------------------------------
C     Resets RDC-angles arrays to hold size entries
C
C     By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003)
C     ------------------------------------------------------------------
      IMPLICIT NONE

      INCLUDE  "funct.inc"
      INCLUDE  "xangle.inc"

      INTEGER OLDSIZE,NEWSIZE

      IF (OLDSIZE.NE.0) THEN
           CALL FREHP(XANGLESJPTR,INTEG4(OLDSIZE))
           CALL FREHP(XANGLESKPTR,INTEG4(OLDSIZE))

           CALL FREHP(XANGLESJLST,INTEG4(OLDSIZE))
           CALL FREHP(XANGLESKLST,INTEG4(OLDSIZE))

           CALL FREHP(XANGLESOBS1PTR,IREAL8(OLDSIZE))
           CALL FREHP(XANGLESOBS2PTR,IREAL8(OLDSIZE))
           CALL FREHP(XANGLESERRPTR,IREAL8(OLDSIZE))
           CALL FREHP(CALCXANGLESPTR,IREAL8(OLDSIZE))
      END IF

      XANGLESJPTR=ALLHP(INTEG4(NEWSIZE))
      XANGLESKPTR=ALLHP(INTEG4(NEWSIZE))

      XANGLESJLST=ALLHP(INTEG4(NEWSIZE))
      XANGLESKLST=ALLHP(INTEG4(NEWSIZE))

      XANGLESOBS1PTR=ALLHP(IREAL8(NEWSIZE))
      XANGLESOBS2PTR=ALLHP(IREAL8(NEWSIZE))
      XANGLESERRPTR=ALLHP(IREAL8(NEWSIZE))
      CALCXANGLESPTR=ALLHP(IREAL8(NEWSIZE))

      RETURN
      END
C     ------------------------------------------------------------------
      SUBROUTINE XANGLESDEFAULTS
C     ------------------------------------------------------------------
C     Sets up defaults
C
C     By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003)
C     ------------------------------------------------------------------
      IMPLICIT NONE

      INCLUDE  "xangle.inc"

      INTEGER COUNT

      MODE=NEW
      MAXXANGLES=200
      XANGLESMX=200
      ANGLESNUM=0
      NCLASSES=0
      CURCLASS=0
      DO COUNT=1,MAXXANGLESCLASSES
           XANGLESCLASSNAMES(COUNT)='DEFAULT'
           XANGLESASSNDX(COUNT)=0
           XANGLESFORCES(COUNT)=5.0
      END DO

      RETURN
      END
C     ------------------------------------------------------------------
      SUBROUTINE READXANGLES2(ATOMJPTR,ATOMKPTR,ATOMJLST,ATOMKLST,
     &                        XANGLESOBS1,XANGLESOBS2,XANGLESERR,SEL)
C     ------------------------------------------------------------------
C     Reads actual RDC-angles assignments into arrays
C
C     By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003)
C     ------------------------------------------------------------------
      IMPLICIT NONE

      INCLUDE  "cns.inc"
      INCLUDE  "mtf.inc"
      INCLUDE  "coord.inc"
      INCLUDE  "xangle.inc"

      INTEGER ATOMJPTR(*),ATOMKPTR(*),ATOMJLST(*),ATOMKLST(*),SEL(*)
      INTEGER JTMP,KTMP,LTMP,II
      DOUBLE PRECISION XANGLESOBS1(*),XANGLESOBS2(*),XANGLESERR(*)
C     ------------------------------------------------------------------
C     Local variables
C     ------------------------------------------------------------------
      INTEGER NSEL,INSERTPOS,COUNT,CURSTOP,OTHERSTOP,NFLAG
      INTEGER I
      DOUBLE PRECISION XANGLESO1,XANGLESO2,XANGLESE
C     ------------------------------------------------------------------
C     If we're in UPDATE mode, make a space for the new line
C     ------------------------------------------------------------------
      NFLAG = 0
      IF (MODE.EQ.UPDATE) THEN
          DO COUNT=ANGLESNUM+1,XANGLESASSNDX(CURCLASS)+1,-1
             ATOMJPTR(COUNT)=ATOMJPTR(COUNT-1)
             ATOMKPTR(COUNT)=ATOMKPTR(COUNT-1)
             XANGLESOBS1(COUNT)=XANGLESOBS1(COUNT-1)
             XANGLESOBS2(COUNT)=XANGLESOBS2(COUNT-1)
             XANGLESERR(COUNT)=XANGLESERR(COUNT-1)
           END DO
           CURSTOP=XANGLESASSNDX(CURCLASS)
           DO COUNT=1,NCLASSES
                OTHERSTOP=XANGLESASSNDX(COUNT)
                IF (OTHERSTOP.GT.CURSTOP) THEN
                     XANGLESASSNDX(COUNT)=OTHERSTOP+1
                END IF
           END DO
           XANGLESASSNDX(CURCLASS)=CURSTOP+1
           INSERTPOS=CURSTOP
           ANGLESNUM=ANGLESNUM + 1
      ELSE
           ANGLESNUM=ANGLESNUM + 1
           INSERTPOS=ANGLESNUM
           XANGLESASSNDX(CURCLASS)=INSERTPOS
      END IF
C     ------------------------------------------------------------------
C     Reading in the atom selection in the constraint table
C     ------------------------------------------------------------------
      CALL SELCTA(SEL,NSEL,X,Y,Z,.TRUE.)
      IF (NSEL.GT.1) THEN
          CALL DSPERR('XANGLE',
     &                'More than 1 atom in SEL for atom J. Using first')
      END IF
      IF (NSEL.EQ.0) THEN
          CALL DSPERR('XANGLE','Error with atom selection')
          NFLAG=1
      END IF
      CALL MAKIND(SEL, NATOM, NSEL)
      IF (INSERTPOS.EQ.1) ATOMJPTR(INSERTPOS)=0
      II=ATOMJPTR(INSERTPOS)
      II=II+1
      IF (II.GT.XANGLESMX) XANGLESMX=II
      ATOMJLST(II)=SEL(1)
      ATOMJPTR(INSERTPOS+1)=II

      CALL SELCTA(SEL,NSEL,X,Y,Z,.TRUE.)
      IF (NSEL.GT.1) THEN
          CALL DSPERR('XANGLE',
     &                'More than 1 atom in SEL for atom K. Using first')
      END IF
      IF (NSEL.EQ.0) THEN
          CALL DSPERR('XANGLE','Error with atom selection')
          NFLAG=1
      END IF
      CALL MAKIND(SEL,NATOM,NSEL)
      IF (INSERTPOS.EQ.1) ATOMKPTR(INSERTPOS)=0
      II=ATOMKPTR(INSERTPOS)
      II=II+1
      IF (II.GT.XANGLESMX) XANGLESMX=II
      ATOMKLST(II)=SEL(1)
      ATOMKPTR(INSERTPOS+1)=II
C     ------------------------------------------------------------------
C     Reading in the observed Theta and Phi
C     ------------------------------------------------------------------
      CALL NEXTF('Observed THETA =',XANGLESO1)
      CALL NEXTF('Observed PHI =',XANGLESO2)
      CALL NEXTF('Error on sine square =',XANGLESE)

      XANGLESOBS1(INSERTPOS)=XANGLESO1
      XANGLESOBS2(INSERTPOS)=XANGLESO2
      XANGLESERR(INSERTPOS)=XANGLESE
C     ------------------------------------------------------------------
C     Check for error in atom selection. If there is one then reset the
C     counter for restraint
C     ------------------------------------------------------------------
      IF (NFLAG.EQ.1) THEN
          ANGLESNUM=ANGLESNUM-1
      END IF

      RETURN
      END
C     ------------------------------------------------------------------
      SUBROUTINE XANGLESINIT
C     ------------------------------------------------------------------
C     Initializes RDC-angles stuff
C
C     By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003)
C     ------------------------------------------------------------------
      IMPLICIT NONE

      INCLUDE  "xangle.inc"

      CALL XANGLESDEFAULTS
      CALL ALLOCXANGLES(0,MAXXANGLES)

      RETURN
      END
C     ------------------------------------------------------------------
      SUBROUTINE PRINTXANGLES (CUTOFF,XANGLESCALC,XANGLESOBS1,
     &                         XANGLESOBS2,XANGLESERR,
     &                         ATOMJPTR,ATOMKPTR,ATOMJLST,ATOMKLST)
C     ------------------------------------------------------------------
C     Prints RDC-angles with deviation of cosine square from 1 larger
C     than cutoff, calculates RMSD and puts it into $RESULT
C
C     By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003)
C     ------------------------------------------------------------------
      IMPLICIT NONE

      INCLUDE  "cns.inc"
      INCLUDE  "xangle.inc"
      INCLUDE  "comand.inc"
      INCLUDE  "mtf.inc"
      INCLUDE  "numbers.inc"

      DOUBLE PRECISION CUTOFF,XANGLESCALC(*),XANGLESOBS1(*),
     &                        XANGLESOBS2(*),XANGLESERR(*)
      INTEGER ATOMJLST(*),ATOMKLST(*)
      INTEGER ATOMJPTR(*),ATOMKPTR(*)
C     ------------------------------------------------------------------
C     Local variables
C     ------------------------------------------------------------------
      DOUBLE PRECISION CALCXANGLES,OBSXANGLES1,OBSXANGLES2,
     &                 DELTAXANGLES,DELTA,DP
      INTEGER COUNT,CLASS,J,K,NUM,II
      DOUBLE PRECISION RMS,VIOLS,ERRXANGLES,XANGLESENERGY
      DOUBLE COMPLEX DUMMY2
      LOGICAL PRINTTHISCLASS

      RMS=ZERO
      VIOLS=ZERO
      NUM=0
C     ------------------------------------------------------------------
C     Make sure that the array of calculated RDC-angles is up to date
C     ------------------------------------------------------------------
      CALL EXANGLES(XANGLESENERGY,'ANALYZE')
      WRITE (PUNIT,'(A)') 'The following RDC-angles have'
      WRITE (PUNIT,'(A)') 'deviation of cosine square from 1'
      WRITE (PUNIT,'(A)') 'larger than or equal to the cutoff:'
C     ------------------------------------------------------------------
C     Write out first class heading
C     ------------------------------------------------------------------
      CLASS=1
      PRINTTHISCLASS=PRINTCLASS(CLASS)
      IF (PRINTTHISCLASS) THEN
           WRITE (PUNIT,'(A,A)') 'Class ',XANGLESCLASSNAMES(1)
      END IF
C     ------------------------------------------------------------------
C     For every RDC-angles entry...
C     ------------------------------------------------------------------
      COUNT=0
      DO WHILE(COUNT.LT.ANGLESNUM)
          COUNT=COUNT+1
C     ------------------------------------------------------------------
C     Is this the start of a new class?
C     ------------------------------------------------------------------
          IF (XANGLESASSNDX(CLASS).LT.COUNT) THEN
              CLASS=CLASS+1
              PRINTTHISCLASS=PRINTCLASS(CLASS)
              IF (XANGLESASSNDX(CLASS).EQ.XANGLESASSNDX(CLASS-1)) THEN
                  COUNT=COUNT-1
              END IF
              IF (PRINTTHISCLASS) THEN
                  WRITE(PUNIT,'(A,A)') 'Class ',XANGLESCLASSNAMES(CLASS)
              END IF
          END IF
C     ------------------------------------------------------------------
C     If this assignment is in a class to be printed
C     and make sure there is an entry for that class
C     ------------------------------------------------------------------
          IF((PRINTTHISCLASS).AND.
     &       (XANGLESASSNDX(CLASS).NE.XANGLESASSNDX(CLASS-1))) THEN
C     ------------------------------------------------------------------
C     Always update RMSD
C     ------------------------------------------------------------------
              CALCXANGLES=XANGLESCALC(COUNT)
              OBSXANGLES1=XANGLESOBS1(COUNT)
              OBSXANGLES2=XANGLESOBS2(COUNT)
              ERRXANGLES=XANGLESERR(COUNT)
              DP=(1.-CALCXANGLES)

              IF ((DP.GT.0.000).AND.(DP.GT.ERRXANGLES)) THEN
                   DELTAXANGLES=DP-ERRXANGLES
              ELSE
                   DELTAXANGLES=0.0
              END IF

              RMS=RMS+DELTAXANGLES**2
              NUM=NUM+1
C     ------------------------------------------------------------------
C     Print out deviations larger than cutoff
C     and update number of violations
C     ------------------------------------------------------------------
              IF (ABS(DELTAXANGLES).GE.CUTOFF) THEN
                  J=ATOMJLST(ATOMJPTR(COUNT)+1)
                  K=ATOMKLST(ATOMKPTR(COUNT)+1)
                  WRITE(PUNIT,'(A,A)') '==============================',
     &                                 '=============================='
                  WRITE(PUNIT,'(A)') ' Set-J-atoms'
                  DO II=ATOMJPTR(COUNT)+1,ATOMJPTR(COUNT+1)
                     J=ATOMJLST(II)
                     WRITE(PUNIT,'(9X,4(1X,A))') SEGID(J),RESID(J),
     &                                           RES(J),TYPE(J)
                  END DO
                  WRITE(PUNIT,'(A)') ' Set-K-atoms'
                  DO II=ATOMKPTR(COUNT)+1,ATOMKPTR(COUNT+1)
                     K=ATOMKLST(II)
                     WRITE(PUNIT,'(9X,4(1X,A))') SEGID(K),RESID(K),
     &                                           RES(K),TYPE(K)
                  END DO
                  WRITE(PUNIT,'(2(2X,A,1X,F8.3))')
     &                  'Calc: ',CALCXANGLES,'Obs: ',1.
                  WRITE(PUNIT,'(2X,A,1X,F8.3,2X,A,1X,F8.3)')
     &                  'Error: ',ERRXANGLES,'Delta: ',DELTAXANGLES
                  VIOLS=VIOLS+ONE
              END IF
          END IF
      END DO

      IF (NUM.GT.0) THEN
          RMS=SQRT(RMS/NUM)
      ELSE
          RMS=0.0
      ENDIF

      CALL DECLAR('RESULT','DP',' ',DUMMY2,RMS)
      CALL DECLAR('VIOLATIONS','DP',' ',DUMMY2,VIOLS)

      RETURN
      END


text:
program test_products
  use mod_orrb
  implicit none

  character(len=*), parameter :: fmt="(A40, 'Time: ',ES8.2,', error: ',ES8.2, ', ', A10)"
  character(len=40) :: test_name

  real(kind=dp) :: t0, t1
  integer(kind=int32) :: ma, na, lbwb, ubwb
  type(error_info) :: error
  real(kind=dp), parameter :: tol=1e-15, c=7.0
  !
  real(kind=dp), dimension(:,:), allocatable :: b_d, a_d, c0_d, c1_d
  complex(kind=dp), dimension(:,:), allocatable :: b_z, a_z, c0_z, c1_z

  type(d_ub), allocatable :: ub_d
  type(d_bv), allocatable :: bv_d
  type(d_wb), allocatable :: wb_d
  type(d_bt), allocatable :: bt_d
  type(d_ubt), allocatable :: ubt_d
  type(d_wbv), allocatable :: wbv_d

  type(z_ub), allocatable :: ub_z
  type(z_bv), allocatable :: bv_z
  type(z_wb), allocatable :: wb_z
  type(z_bt), allocatable :: bt_z
  type(z_ubt), allocatable :: ubt_z
  type(z_wbv), allocatable :: wbv_z

  call initialize_errors
  print *
  print *, "--------------------------------"
  print *
  print *, "Tests of Real Products"
  print *

  ma=40; na=10; lbwb=5; ubwb=7
  a_d=d_random_matrix(ma,na)
  ub_d=d_random_ub(ma,lbwb,ubwb,error=error)
  b_d = general(ub_d,error)
  call cpu_time(t0)
  c0_d=d_product_of_ub_and_general(ub_d,a_d,error)
  call cpu_time(t1)
  c1_d=matmul(b_d,a_d)
  test_name = "Real UB Times Gen. (40x10)"
  call d_output_result(test_name,c0_d,c1_d,t0,t1,c*tol,error)

  ma=40; na=10; lbwb=5; ubwb=7
  a_d=d_random_matrix(ma,na)
  bv_d=d_random_bv(ma,lbwb,ubwb,error=error)
  b_d = general(bv_d,error)
  call cpu_time(t0)
  c0_d=d_product_of_bv_and_general(bv_d,a_d,error)
  call cpu_time(t1)
  c1_d=matmul(b_d,a_d)
  test_name = "Real BV Times Gen. (40x10)"
  call d_output_result(test_name,c0_d,c1_d,t0,t1,c*tol,error)

  ma=40; na=10; lbwb=5; ubwb=7
  a_d=d_random_matrix(ma,na)
  wb_d=d_random_wb(ma,lbwb,ubwb,error=error)
  b_d = general(wb_d,error)
  call cpu_time(t0)
  c0_d=d_product_of_wb_and_general(wb_d,a_d,error)
  call cpu_time(t1)
  c1_d=matmul(b_d,a_d)
  test_name = "Real WB Times Gen. (40x10)"
  call d_output_result(test_name,c0_d,c1_d,t0,t1,c*tol,error)

  ma=40; na=10; lbwb=5; ubwb=7
  a_d=d_random_matrix(ma,na)
  bt_d=d_random_bt(ma,lbwb,ubwb,error=error)
  b_d = general(bt_d,error)
  call cpu_time(t0)
  c0_d=d_product_of_bt_and_general(bt_d,a_d,error)
  call cpu_time(t1)
  c1_d=matmul(b_d,a_d)
  test_name = "Real BT Times Gen. (40x10)"
  call d_output_result(test_name,c0_d,c1_d,t0,t1,c*tol,error)

  ma=40; na=10; lbwb=5; ubwb=7
  a_d=d_random_matrix(ma,na)
  ubt_d=d_random_ubt(ma,lbwb,ubwb,error=error)
  b_d = general(ubt_d,error)
  call cpu_time(t0)
  c0_d=d_product_of_ubt_and_general(ubt_d,a_d,error)
  call cpu_time(t1)
  c1_d=matmul(b_d,a_d)
  test_name = "Real UBT Times Gen. (40x10)"
  call d_output_result(test_name,c0_d,c1_d,t0,t1,c*tol,error)

  ma=40; na=10; lbwb=5; ubwb=7
  a_d=d_random_matrix(ma,na)
  wbv_d=d_random_wbv(ma,lbwb,ubwb,error=error)
  b_d = general(wbv_d,error)
  call cpu_time(t0)
  c0_d=d_product_of_wbv_and_general(wbv_d,a_d,error)
  call cpu_time(t1)
  c1_d=matmul(b_d,a_d)
  test_name = "Real WBV Times Gen. (40x10)"
  call d_output_result(test_name,c0_d,c1_d,t0,t1,c*tol,error)

  ma=10; na=40; lbwb=5; ubwb=7
  a_d=d_random_matrix(ma,na)
  ub_d=d_random_ub(na,lbwb,ubwb,error=error)
  b_d = general(ub_d,error)
  call cpu_time(t0)
  c0_d=d_product_of_general_and_ub(a_d, ub_d,error)
  call cpu_time(t1)
  c1_d=matmul(a_d,b_d)
  test_name = "Real Gen. Times UB. (10x40)"
  call d_output_result(test_name,c0_d,c1_d,t0,t1,c*tol,error)

  print *
  print *, "--------------------------------"
  print *
  print *, "Tests of Complex Products"
  print *

  ma=40; na=10; lbwb=5; ubwb=7
  a_z=z_random_matrix(ma,na)
  ub_z=z_random_ub(ma,lbwb,ubwb,error=error)
  b_z = general(ub_z,error)
  call cpu_time(t0)
  c0_z=z_product_of_ub_and_general(ub_z,a_z,error)
  call cpu_time(t1)
  c1_z=matmul(b_z,a_z)
  test_name = "Complex UB Times Gen. (40x10)"
  call z_output_result(test_name,c0_z,c1_z,t0,t1,c*tol,error)

  ma=40; na=10; lbwb=5; ubwb=7
  a_z=z_random_matrix(ma,na)
  bv_z=z_random_bv(ma,lbwb,ubwb,error=error)
  b_z = general(bv_z,error)
  call cpu_time(t0)
  c0_z=z_product_of_bv_and_general(bv_z,a_z,error)
  call cpu_time(t1)
  c1_z=matmul(b_z,a_z)
  test_name = "Complex BV Times Gen. (40x10)"
  call z_output_result(test_name,c0_z,c1_z,t0,t1,c*tol,error)

  ma=40; na=10; lbwb=5; ubwb=7
  a_z=z_random_matrix(ma,na)
  wb_z=z_random_wb(ma,lbwb,ubwb,error=error)
  b_z = general(wb_z,error)
  call cpu_time(t0)
  c0_z=z_product_of_wb_and_general(wb_z,a_z,error)
  call cpu_time(t1)
  c1_z=matmul(b_z,a_z)
  test_name = "Complex WB Times Gen. (40x10)"
  call z_output_result(test_name,c0_z,c1_z,t0,t1,c*tol,error)

  ma=40; na=10; lbwb=5; ubwb=7
  a_z=z_random_matrix(ma,na)
  bt_z=z_random_bt(ma,lbwb,ubwb,error=error)
  b_z = general(bt_z,error)
  call cpu_time(t0)
  c0_z=z_product_of_bt_and_general(bt_z,a_z,error)
  call cpu_time(t1)
  c1_z=matmul(b_z,a_z)
  test_name = "Complex BT Times Gen. (40x10)"
  call z_output_result(test_name,c0_z,c1_z,t0,t1,c*tol,error)

  ma=40; na=10; lbwb=5; ubwb=7
  a_z=z_random_matrix(ma,na)
  ubt_z=z_random_ubt(ma,lbwb,ubwb,error=error)
  b_z = general(ubt_z,error)
  call cpu_time(t0)
  c0_z=z_product_of_ubt_and_general(ubt_z,a_z,error)
  call cpu_time(t1)
  c1_z=matmul(b_z,a_z)
  test_name = "Complex UBT Times Gen. (40x10)"
  call z_output_result(test_name,c0_z,c1_z,t0,t1,c*tol,error)

  ma=40; na=10; lbwb=5; ubwb=7
  a_z=z_random_matrix(ma,na)
  wbv_z=z_random_wbv(ma,lbwb,ubwb,error=error)
  b_z = general(wbv_z,error)
  call cpu_time(t0)
  c0_z=z_product_of_wbv_and_general(wbv_z,a_z,error)
  call cpu_time(t1)
  c1_z=matmul(b_z,a_z)
  test_name = "Complex WBV Times Gen. (40x10)"
  call z_output_result(test_name,c0_z,c1_z,t0,t1,c*tol,error)

  ma=10; na=40; lbwb=5; ubwb=7
  a_z=z_random_matrix(ma,na)
  ub_z=z_random_ub(na,lbwb,ubwb,error=error)
  b_z = general(ub_z,error)
  call cpu_time(t0)
  c0_z=z_product_of_general_and_ub(a_z, ub_z,error)
  call cpu_time(t1)
  c1_z=matmul(a_z,b_z)
  test_name = "Complex Gen. Times UB. (10x40)"
  call z_output_result(test_name,c0_z,c1_z,t0,t1,c*tol,error)

contains

  subroutine d_output_result(name,a0,a1,t0,t1,bnd,error)
    character(len=*) :: name
    real(kind=dp), dimension(:,:) :: a0, a1
    real(kind=dp) :: bnd, t0, t1
    type(error_info) :: error

    real(kind=dp) :: err
    character(len=10) :: test_result

    if (error%code > 0) then
       print *, "Calling error in test: ", name
    else
       err = maxabs(a1-a0)
       if (err < bnd) then
          test_result="PASSED"
       else
          test_result="    FAILED"
       end if
       write (*,fmt) name, t1-t0, err, test_result
    end if
  end subroutine d_output_result

  subroutine z_output_result(name,a0,a1,t0,t1,bnd,error)
    character(len=*) :: name
    complex(kind=dp), dimension(:,:) :: a0, a1
    real(kind=dp) :: bnd, t0, t1
    type(error_info) :: error

    real(kind=dp) :: err
    character(len=10) :: test_result

    if (error%code > 0) then
       print *, "Calling error in test: ", name
    else
       err = maxabs(a1-a0)
       if (err < bnd) then
          test_result="PASSED"
       else
          test_result="    FAILED"
       end if
       write (*,fmt) name, t1-t0, err, test_result
    end if
  end subroutine z_output_result

end program test_products

max_stars_repo_name: m-a-stewart/BandGivensWeight

max_issues_repo_name: m-a-stewart/BandGivensWeight

max_forks_repo_name: m-a-stewart/BandGivensWeight


text:
      double precision function dnrm2(n,x,incx)
      integer n, incx
      double precision x(n)
c     **********
c
c     Function dnrm2
c
c     Given a vector x of length n, this function calculates the
c     Euclidean norm of x with stride incx.
c
c     The function statement is
c
c       double precision function dnrm2(n,x,incx)
c
c     where
c
c       n is an integer variable.
c         On entry n is the length of x.
c         On exit n is unchanged.
c
c       x is a double precision array of dimension n.
c         On entry x specifies the vector x.
c         On exit x is unchanged.
c
c       incx is an integer variable.
c         On entry incx specifies the stride.
c         On exit incx is unchanged.
c
c     MINPACK-2 Project. October 1993.
c     Argonne National Laboratory and University of Minnesota.
c     Brett M. Averick and Jorge J. More'.
c
c     **********
      double precision zero
      parameter (zero=0.0d0)

      integer i
      double precision scale

      scale = zero
      do 10 i = 1, n, incx
         scale = max(scale,abs(x(i)))
   10 continue

      dnrm2 = zero
      if (scale .eq. zero) return

      do 20 i = 1, n, incx
         dnrm2 = dnrm2 + (x(i)/scale)**2
   20 continue
      dnrm2 = scale*sqrt(dnrm2)

      end

ext: f

max_forks_repo_name: dynamics-of-stellar-systems/dynamite


text:
************************************************************************
*
*     Subroutine INIT                          Called by: Main program
*
*     initialize, open files, determine the length of the calibration
*     and estimation files, and read the first part of the header file
*
************************************************************************
      SUBROUTINE INIT(LDOPT,NOBSCI,NOBSE,NSEAS,PRTOPT,SEOPT,TITLE)
*
*     subroutine args
*
      INTEGER*4 LDOPT,NOBSCI,NOBSE,NSEAS,PRTOPT,SEOPT
      CHARACTER*80 TITLE
*
*     initialize Logical Device Assignments (LDAs)
*
      CALL LDAINIT
*
*     open echo.out and all the input files; determine the length of the
*     calibration and estimation files.
*
      CALL OPENFILES(NOBSCI,NOBSE)
*
*     read the title, print option, standard error option, etc. from the
*     header file.
*
      CALL INPUTH(TITLE,LDOPT,PRTOPT,SEOPT,NSEAS)

      RETURN
      END

max_stars_repo_name: lthiamodelers/baseflow-coefficients

max_issues_repo_name: lthiamodelers/baseflow-coefficients

max_forks_repo_name: lthiamodelers/baseflow-coefficients


     subroutine header_wetland

     use basin_module
     use reservoir_module
     use hydrograph_module

     implicit none

    !! RESERVOIR/WETLAND - DAILY
      if (pco%res%d == "y") then
        open (2548,file="wetland_day.txt",recl=1500)
        write (2548,*) bsn%name, prog
        write (9000,*) "RES_WET                   wetland_day.txt"
        write (2548,*) ch_wbod_hdr, hyd_stor_hdr, hyd_in_hdr, hyd_out_hdr
        write (2548,*) ch_wbod_hdr_units, hyd_hdr_units1, hyd_hdr_units1, hyd_hdr_units1
          if (pco%csvout == "y") then
            open (2552,file="wetland_day.csv",recl=1500)
            write (2552,*) bsn%name, prog
            write (2552,'(*(G0.3,:","))') ch_wbod_hdr, hyd_stor_hdr, hyd_in_hdr, hyd_out_hdr
            write (2552,'(*(G0.3,:","))') ch_wbod_hdr_units, hyd_hdr_units1, hyd_hdr_units1, hyd_hdr_units1
            write (9000,*) "RES_WET                   wetland_day.csv"
          end if
      end if

    !! RESERVOIR/WETLAND - MONTHLY
      if (pco%res%m == "y") then
        open (2549,file="wetland_mon.txt",recl=1500)
        write (2549,*) bsn%name, prog
        write (9000,*) "RES_WET                   wetland_mon.txt"
        write (2549,*) ch_wbod_hdr, hyd_stor_hdr, hyd_in_hdr, hyd_out_hdr
        write (2549,*) ch_wbod_hdr_units, hyd_hdr_units1, hyd_hdr_units1, hyd_hdr_units1
          if (pco%csvout == "y") then
            open (2553,file="wetland_mon.csv",recl=1500)
            write (2553,*) bsn%name, prog
            write (2553,'(*(G0.3,:","))') ch_wbod_hdr, hyd_stor_hdr, hyd_in_hdr, hyd_out_hdr
            write (2553,'(*(G0.3,:","))') ch_wbod_hdr_units, hyd_hdr_units1, hyd_hdr_units1, hyd_hdr_units1
            write (9000,*) "RES_WET                   wetland_mon.csv"
          end if
      end if

   !! RESERVOIR/WETLAND YEARLY
     if (pco%res%y == "y") then
        open (2550,file="wetland_yr.txt",recl=1500)
        write (2550,*) bsn%name, prog
        write (9000,*) "RES_WET                   wetland_yr.txt"
        write (2550,*) ch_wbod_hdr, hyd_stor_hdr, hyd_in_hdr, hyd_out_hdr
        write (2550,*) ch_wbod_hdr_units, hyd_hdr_units1, hyd_hdr_units1, hyd_hdr_units1
          if (pco%csvout == "y") then
            open (2554,file="wetland_yr.csv",recl=1500)
            write (2554,*) bsn%name, prog
            write (2554,'(*(G0.3,:","))') ch_wbod_hdr, hyd_stor_hdr, hyd_in_hdr, hyd_out_hdr
            write (2554,'(*(G0.3,:","))') ch_wbod_hdr_units, hyd_hdr_units1, hyd_hdr_units1, hyd_hdr_units1
            write (9000,*) "RES_WET                   wetland_yr.csv"
          end if
     end if

    !! RESERVOIR/WETLAND - AVERAGE ANNUAL
      if (pco%res%a == "y") then
        open (2551,file="wetland_aa.txt",recl=1500)
        write (2551,*) bsn%name, prog
        write (2551,*) ch_wbod_hdr, hyd_stor_hdr, hyd_in_hdr, hyd_out_hdr
        write (2551,*) ch_wbod_hdr_units, hyd_hdr_units1, hyd_hdr_units1, hyd_hdr_units1
        write (9000,*) "RES_WET                   wetland_aa.txt"
          if (pco%csvout == "y") then
            open (2555,file="wetland_aa.csv",recl=1500)
            write (2555,*) bsn%name, prog
            write (2555,'(*(G0.3,:","))') ch_wbod_hdr, hyd_stor_hdr, hyd_in_hdr, hyd_out_hdr
            write (2555,'(*(G0.3,:","))') ch_wbod_hdr_units, hyd_hdr_units1, hyd_hdr_units1, hyd_hdr_units1
            !write (2555,'(*(G0.3,:","))') res_hdr_add
            write (9000,*) "RES_WET                   wetland_aa.csv"
          end if
      end if

      return
     end subroutine header_wetland


text:
MODULE m_umtx
   !*********************************************************************
   !* The calculation of the "U"-contribution to Hartree-Fock matrix.   *
   !*-------------------------------------------------------------------*
   !* Extension to multiple U per atom type by G.M. 2017                *
   !*********************************************************************
   USE m_juDFT
   USE m_constants
   USE m_sgaunt
   USE m_types

   IMPLICIT NONE

   CONTAINS

   SUBROUTINE umtx(u_in,n_u,f0,f2,f4,f6,u)

      INTEGER,       INTENT(IN)  :: n_u
      TYPE(t_utype), INTENT(IN)  :: u_in(:)
      REAL,          INTENT(IN)  :: f0(:),f2(:),f4(:),f6(:)
      REAL,          INTENT(OUT) :: u(-lmaxU_const:,-lmaxU_const:,-lmaxU_const:,-lmaxU_const:,:)

      INTEGER, PARAMETER :: lmmaxw1=(2*lmaxU_const+2)**2
      REAL,    PARAMETER :: tol=1e-14

      INTEGER :: i,j,k,l,m,mk,nfk,itype,i_u
      INTEGER :: m1,m2,m3,m4,lm1,lm2,lm3,lm4,kf
      REAL    :: uk,uq,avu,avj,cgk1,cgk2
      REAL    :: fk(lmaxU_const+1,n_u)
      REAL,   ALLOCATABLE :: c(:,:,:)

      !
      ! transformation to Hr-units:
      !
      DO i_u = 1, n_u
         itype = u_in(i_u)%atomType
         l = u_in(i_u)%l
         fk(1,i_u) = f0(i_u) / hartree_to_ev_const
         fk(2,i_u) = f2(i_u) / hartree_to_ev_const
         fk(3,i_u) = f4(i_u) / hartree_to_ev_const
         IF (l.EQ.3) THEN
            fk(4,i_u) = f6(i_u) / hartree_to_ev_const
         ELSE IF (l.GT.3) THEN
            CALL juDFT_error("LDA+U for p, d or f-states!", calledby="umtx")
         END IF
      END DO

      !
      ! evaluate Gaunt parameter
      !
      ALLOCATE(c(0:2*lmaxU_const+1,lmmaxw1,lmmaxw1),source=0.0)
      CALL sgaunt(lmaxU_const,c)

      DO i_u = 1, n_u
         l = u_in(i_u)%l
         kf = 2*l
         DO m1 = -l,l
            lm1 = l*(l+1)+m1+1
            DO m2 = -l,l
               lm2 = l*(l+1)+m2+1
               DO m3 = -l,l
                  lm3 = l*(l+1)+m3+1
                  DO m4 = -l,l
                     lm4 = l*(l+1)+m4+1
                     uk = 0.0e0
                     DO k=0,kf,2
                        uq = 0.e0
                        DO mk=-k,k
                           IF (mk.NE.m1-m3)  CYCLE
                           cgk1 = c(k/2,lm1,lm3)
                           IF (ABS(cgk1).LT.tol) CYCLE
                           IF (mk.NE.m4-m2)  CYCLE
                           cgk2 = c(k/2,lm4,lm2)
                           IF (ABS(cgk2).LT.tol) CYCLE
                           uq = uq+cgk1*cgk2
                        END DO                   ! mk
                        IF (ABS(uq).LT.tol) CYCLE
                        nfk=k/2+1
                        uk=uk+uq*fk(nfk,i_u)*4*pi_const/(2*k+1)
                     END DO                     ! k
                     u(m1,m2,m3,m4,i_u)=uk
                  END DO                       ! m4 etc.
               END DO
            END DO
         END DO
         avu=0.e0
         avj=0.e0

         DO i = -l,l
            DO j = -l,l
               avu = avu+u(i,j,i,j,i_u)
               avj = avj+(u(i,j,i,j,i_u)-u(i,j,j,i,i_u))
            END DO
         END DO
         avu = avu/(2*l+1)/(2*l+1)
         avj = avj/(2*l+1)/(2*l)
         avj = avu-avJ
         !        WRITE (oUnit,*) 'U-matr:'
         !        IF (l.eq.2) WRITE (oUnit,111) ((u(i,j,i,j,i_u),i=-l,l),j=-l,l)
         !        IF (l.eq.3) WRITE (oUnit,211) ((u(i,j,i,j,i_u),i=-l,l),j=-l,l)
         !        WRITE (oUnit,*) 'J-matr:'
         !        IF (l.eq.2) WRITE (oUnit,111) ((u(i,j,j,i,i_u),i=-l,l),j=-l,l)
         !        IF (l.eq.3) WRITE (oUnit,211) ((u(i,j,j,i,i_u),i=-l,l),j=-l,l)
         !         PRINT*,'U-av:',avu*hartree_to_ev_const
         !         PRINT*,'J-av:',avj*hartree_to_ev_const
111      FORMAT (5f8.4)
211      FORMAT (7f8.4)
112      FORMAT (10e20.10)
         !c         WRITE (9,112) ((((u(i,j,k,m,i_u),i=-l,l),j=-l,l),k=-l,l),m=-l,l)
      END DO
      DEALLOCATE (c)

   END SUBROUTINE umtx
END MODULE m_umtx

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR


text:
module unum_size_limit_mod
  implicit none

  integer, parameter :: MIN_ESIZESIZE = 0
  integer, parameter :: MAX_ESIZESIZE = 4

  integer, parameter :: MIN_FSIZESIZE = 0
  integer, parameter :: MAX_FSIZESIZE = 11

end module unum_size_limit_mod

max_stars_repo_name: CrashBurnRepeat/fUnums

max_issues_repo_name: CrashBurnRepeat/fUnums

max_forks_repo_name: CrashBurnRepeat/fUnums


text:
      subroutine sub(num1, *, num2, *)
      if (num1+num2 .eq. 2) return 1
      if (num1 .eq. 3) return 2
      if (num1 .eq. 4) return 0
      if (num1 .eq. 5) return 3
      return
      end

      program prog
      do i=0,5
        call sub(i, *10, 1, *20)
        write(6, '(a)') 'regular'
        goto 30
   10   write(6, '(a)') 'goto 10'
        goto 30
   20   write(6, '(a)') 'goto 20'
   30   continue
      enddo
      end


text:
      subroutine wetlan

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine simulates wetlands

!!    ~ ~ ~ INCOMING VARIABLES ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    bw1(:)      |none          |1st shape parameter for wetland surface area
!!                               |equation
!!    bw2(:)      |none          |2nd shape parameter for the wetland surface
!!                               |area equation
!!    chlaw(:)    |none          |chlorophyll-a production coefficient for
!!                               |wetland
!!    hru_dafr(:) |none          |fraction of watershed area in HRU
!!    hru_ha(:)   |ha            |area of HRU in hectares
!!    ihru        |none          |HRU number
!!    ipnd1(:)    |none          |beginning month of 2nd "season" of nutrient
!!                               |settling
!!    ipnd2(:)    |none          |ending month of 2nd "season" of nutrient
!!                               |settling
!!    minpgw(:)   |kg P/ha       |soluble P loading to reach in groundwater
!!    no3gw(:)    |kg N/ha       |nitrate loading to reach in groundwater
!!    nsetlw(1,:) |m/day         |nitrogen settling rate for 1st season
!!    nsetlw(2,:) |m/day         |nitrogen settling rate for 2nd season
!!    pet_day     |mm H2O        |potential evapotranspiration for day in HRU
!!    psetlw(1,:) |m/day         |phosphorus settling rate for 1st season
!!    psetlw(2,:) |m/day         |phosphorus settling rate for 2nd season
!!    qdr(:)      |mm H2O        |net water loading from HRU to main channel
!!    secciw(:)   |none          |water clarity coefficient for wetland
!!    sed_stl(:)  |kg/kg         |fraction of sediment remaining suspended in
!!                               |impoundment after settling for one day
!!    sedminpa(:) |kg P/ha       |amount of active mineral phosphorus sorbed to
!!                               |sediment in surface runoff in HRU for day
!!    sedminps(:) |kg P/ha       |amount of stable mineral phosphorus sorbed to
!!                               |sediment in surface runoff in HRU for day
!!    sedorgn(:)  |kg N/ha       |amount of organic nitrogen in surface runoff
!!                               |in HRU for the day
!!    sedorgp(:)  |kg P/ha       |amount of organic phosphorus in surface runoff
!!                               |in HRU for the day
!!    sedyld(:)   |metric tons   |daily soil loss caused by water erosion in HRU
!!    shallst(:)  |mm H2O        |depth of water in shallow aquifer
!!    subp(:)     |mm H2O        |precipitation for the day in HRU
!!    surqno3(:)  |kg N/ha       |amount of NO3-N in surface runoff in HRU for
!!                               |the day
!!    surqsolp(:) |kg P/ha       |amount of soluble phosphorus in surface runoff
!!                               |in HRU for the day
!!    wet_fr(:)   |none          |fraction of HRU/subbasin area that drains
!!                               |into wetlands
!!    wet_k(:)    |mm/hr         |hydraulic conductivity of bottom of wetlands
!!    wet_mxvol(:)|m^3 H2O       |volume of water required to fill wetlands to
!!                               |maximum water level
!!    wet_no3(:)  |kg N          |amount of nitrate originating from surface
!!                               |runoff in wetland at beginning of day
!!    wet_no3g(:) |kg N          |amount of nitrate originating from
!!                               |groundwater in wetland at beginning of day
!!    wet_no3s(:) |kg N          |amount of nitrate originating from lateral
!!                               |flow in wetland at beginning of day
!!    wet_nsed(:) |kg/L          |normal sediment concentration in wetland water
!!    wet_nvol(:) |m^3 H2O       |volume of water needed to
!!                               |fill wetlands to normal water level
!!    wet_orgn(:) |kg N          |amount of organic N originating from
!!                               |surface runoff in wetland at beginning of day
!!    wet_orgp(:) |kg P          |amount of organic P originating from
!!                               |surface runoff in wetland at beginning of day
!!    wet_psed(:) |kg P          |amount of mineral P attached to sediment
!!                               |originating from surface runoff in wetland at
!!                               |beginning of day
!!    wet_sed(:)  |kg/L          |sediment concentration in wetland water
!!    wet_solp(:) |kg P          |amount of soluble P originating from surface
!!                               |runoff in wetland at beginning of day
!!    wet_solpg(:)|kg P          |amount of soluble P originating from
!!                               |groundwater in wetland at beginning of day
!!    wet_vol(:)  |m^3 H2O       |volume of water in wetlands
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ OUTGOING VARIABLES ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    qdr(:)      |mm H2O        |net water loading from HRU to main channel
!!    sedyld(:)   |metric tons   |daily soil loss caused by water erosion in HRU
!!    shallst(:)  |mm H2O        |depth of water in shallow aquifer
!!    twlwet      |mm H2O        |water lost through seepage from wetlands on
!!                               |day in HRU
!!    wet_chla(:) |kg chla       |amount of chlorophyll-a in wetland at end
!!                               |of day
!!    wet_no3(:)  |kg N          |amount of nitrate originating from surface
!!                               |runoff in wetland at end of day
!!    wet_no3g(:) |kg N          |amount of nitrate originating from
!!                               |groundwater in wetland at end of day
!!    wet_no3s(:) |kg N          |amount of nitrate originating from lateral
!!                               |flow in wetland at end of day
!!    wet_orgn(:) |kg N          |amount of organic N originating from
!!                               |surface runoff in wetland at end of day
!!    wet_orgp(:) |kg P          |amount of organic P originating from
!!                               |surface runoff in wetland at end of day
!!    wet_psed(:) |kg P          |amount of mineral P attached to sediment
!!                               |originating from surface runoff in wetland at
!!                               |end of day
!!    wet_seci(:) |m             |secchi-disk depth in wetland at end of day
!!    wet_sed(:)  |kg/L          |sediment concentration in wetland water
!!    wet_solp(:) |kg P          |amount of soluble P originating from surface
!!                               |runoff in wetland at end of day
!!    wet_solpg(:)|kg P          |amount of soluble P originating from
!!                               |groundwater in wetland at end of day
!!    wet_vol(:)  |m^3 H2O       |volume of water in wetlands
!!    wetev       |m^3 H2O       |evaporation from wetland for day
!!    wetflwi     |m^3 H2O       |volume of water flowing in wetland on day
!!    wetflwo     |m^3 H2O       |volume of water flowing out wetland on day
!!    wetpcp      |m^3 H2O       |precipitation on wetland for day
!!    wetsedc     |metric tons   |net change in sediment in wetland on day
!!    wetsedi     |metric tons   |sediment loading to wetland for day
!!    wetsedo     |metric tons   |sediment loading from wetland for day
!!    wetsep      |m^3 H2O       |seepage from wetland bottom for day
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ LOCAL DEFINITIONS ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    cnv         |none          |conversion factor (mm/ha => m^3)
!!    iseas       |none          |nutrient settling rate season
!!    j           |none          |HRU number
!!    nitrok      |none          |fraction of nitrogen in wetland removed by
!!                               |settling
!!    phosk       |none          |fraction of phosphorus in wetland removed by
!!                               |settling
!!    sed         |kg/kg         |sediment concentration in wetland at beginning
!!                               |of day
!!    tpco        |ppb (ug/L)    |concentration of phosphorus in pond water
!!                               |on day
!!    vol         |m^3 H2O       |volume of wetland at beginning of day
!!    wetsa       |ha            |surface area of wetland on current day
!!    xx          |none          |variable to hold intermediate calculation
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

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

      use parm

      integer :: j, iseas
      real*8 :: vol, cnv, sed, wetsa, xx, phosk, nitrok, tpco
      real*8 :: wetsani, wetsili, wetclai, wetsagi, wetlagi
	real*8 :: san, sil, cla, sag, lag, inised, finsed,setsed,remsetsed
      real*8 :: wetsano, wetsilo, wetclao, wetsago, wetlago
      real*8 :: qdayi, latqi

      j = 0
      j = ihru

      if (wet_fr(j) > 0.) then
        cnv = 0.
        cnv = hru_ha(j) * 10.               !conversion factor

        !! store initial values
        vol = 0.
        sed = 0.
        san = 0.
        sil = 0.
        cla = 0.
        sag = 0.
        lag = 0.
        inised = 0.
        finsed = 0.
        setsed = 0.
        remsetsed = 0.

        vol = wet_vol(j)
        sed = wet_sed(j)
        san = wet_san(j)
        sil = wet_sil(j)
        cla = wet_cla(j)
        sag = wet_sag(j)
        lag = wet_lag(j)

        !! calculate water balance for day
        wetsa = 0.
        wetsa = bw1(j) * wet_vol(j) ** bw2(j)

        wetev = 10. * evwet(j) * pet_day * wetsa
        wetsep = wet_k(j) * wetsa * 240.
        wetpcp = subp(j) * wetsa * 10.

        !! calculate water flowing into wetland from HRU
        fr_cur = wet_fr(j) * ((hru_ha(j) - wetsa) / hru_ha(j))
        fr_cur = Max(0.,fr_cur)
        wetflwi = qday + latq(j)

        if (iwetile(j) == 1) then
          wetflwi = wetflwi + qtile
          qtile = qtile * (1. - fr_cur)
        end if
        if (iwetgw(j) == 1) then
          wetflwi = wetflwi + gw_q(j)
          gw_q(j) = gw_q(j) * (1. - fr_cur)
        end if

        wetflwi = wetflwi * 10. * (hru_ha(j) * wet_fr(j) - wetsa)
        qdayi = qday
        latqi = latq(j)
        qday = qday * (1. - wet_fr(j))
        latq(j) = latq(j) * (1. - wet_fr(j))
        wetloss = qdayi - qday
        lwetloss = latqi - latq(j)

        qdr(j) = qdr(j) - wetloss - lwetloss
!       qdr(j) = qdr(j) - qdr(j) * wet_fr(j)

        !! sediment loading to wetland from HRU
        wetsedi = sedyld(j) * (wet_fr(j) - (wetsa / hru_ha(j)))

        wetsani = sanyld(j) * (wet_fr(j) - (wetsa / hru_ha(j)))
        wetsili = silyld(j) * (wet_fr(j) - (wetsa / hru_ha(j)))
        wetclai = clayld(j) * (wet_fr(j) - (wetsa / hru_ha(j)))
        wetsagi = sagyld(j) * (wet_fr(j) - (wetsa / hru_ha(j)))
        wetlagi = lagyld(j) * (wet_fr(j) - (wetsa / hru_ha(j)))

        sedyld(j) = sedyld(j) - sedyld(j) * wet_fr(j)
        sanyld(j) = sanyld(j) - sanyld(j) * wet_fr(j)
        silyld(j) = silyld(j) - silyld(j) * wet_fr(j)
        clayld(j) = clayld(j) - clayld(j) * wet_fr(j)
        sagyld(j) = sagyld(j) - sagyld(j) * wet_fr(j)
        lagyld(j) = lagyld(j) - lagyld(j) * wet_fr(j)

        !! compute nitrogen and phosphorus levels in wetland at beginning
        !! of day
        !! equation 29.1.1 in SWAT manual
        xx = 0.
        xx = wet_fr(j) * hru_ha(j)
        wet_solp(j) = wet_solp(j) + (surqsolp(j) + sedminpa(j)) * xx
        wet_psed(j) = wet_psed(j) + sedminps(j) * xx
        wet_orgp(j) = wet_orgp(j) + sedorgp(j) * xx
        wet_solpg(j) = wet_solpg(j) + minpgw(j) * xx
        wet_orgn(j) = wet_orgn(j) + sedorgn(j) * xx
        wet_no3(j) = wet_no3(j) + surqno3(j) * xx
        wet_no3s(j) = wet_no3s(j) + latno3(j) * xx
        wet_no3g(j) = wet_no3g(j) + no3gw(j) * xx

        !! remove nutrients entering wetlands from HRU loadings
        xx = 0.
        xx = 1. - wet_fr(j)
        sedorgn(j) = sedorgn(j) * xx
        surqno3(j) = surqno3(j) * xx
        latno3(j) = latno3(j) * xx
        no3gw(j) = no3gw(j) * xx
        sedorgp(j) = sedorgp(j) * xx
        sedminpa(j) = sedminpa(j) * xx
        sedminps(j) = sedminps(j) * xx
        surqsolp(j) = surqsolp(j) * xx
        minpgw(j) = minpgw(j) * xx

        !! new water volume for day
        wet_vol(j) = wet_vol(j) - wetsep - wetev + wetpcp + wetflwi

        if (wet_vol(j) < 0.001) then

          !! check for volume deficit in wetland
          !! reduce seepage so that the wetland volume is zero
          wetsep = wetsep + wet_vol(j)
          wet_vol(j) = 0.
          !! if seepage is less than the volume deficit, take the
          !! remainder from evaporation
          if (wetsep < 0.) then
            wetev = wetev + wetsep
            wetsep = 0.
          end if
          wet_sed(j) = 0.

          wet_san(j) = 0.
          wet_sil(j) = 0.
          wet_cla(j) = 0.
          wet_sag(j) = 0.
          wet_lag(j) = 0.

          wet_solp(j) = 0.
          wet_psed(j) = 0.
          wet_orgp(j) = 0.
          wet_solpg(j) = 0.
          wet_orgn(j) = 0.
          wet_no3(j) = 0.
          wet_no3s(j) = 0.
          wet_no3g(j) = 0.
          wet_chla(j) = 0.
          wet_seci(j) = 0.

        else

          !! compute new sediment concentration
          wet_sed(j) = (sed * vol + wetsedi) / wet_vol(j)

          wet_san(j) = (san * vol + wetsani) / wet_vol(j)
          wet_sil(j) = (sil * vol + wetsili) / wet_vol(j)
          wet_cla(j) = (cla * vol + wetclai) / wet_vol(j)
          wet_sag(j) = (sag * vol + wetsagi) / wet_vol(j)
          wet_lag(j) = (lag * vol + wetlagi) / wet_vol(j)

          !! compute outflow if wetland water volume > 0
          if (wet_vol(j) <= wet_nvol(j)) then
            wetflwo = 0.
          else
            if (wet_vol(j) <= wet_mxvol(j)) then
              wetflwo = (wet_vol(j) - wet_nvol(j)) / 10.
              wet_vol(j) = wet_vol(j) - wetflwo
            else
              wetflwo = wet_vol(j) - wet_mxvol(j)
              wet_vol(j) = wet_mxvol(j)
            end if
          end if
          qday= qday + wetflwo / cnv
          qdr(j) = qdr(j) + wetflwo / cnv

          !! compute sediment settling
	    if (sed_stl(j) < 1.e-6) sed_stl(j) = 0.0
          inised = wet_sed(j)
          if (wet_sed(j) > wet_nsed(j)) then
            wet_sed(j) = (wet_sed(j) - wet_nsed(j)) * sed_stl(j) +
     &                                                       wet_nsed(j)
          end if
          finsed = wet_sed(j)
	    setsed = inised - finsed
          setsed = Max(0.,setsed)

	    if (wet_lag(j) >= setsed) then
	      wet_lag(j) = wet_lag(j) - setsed
	      remsetsed = 0.
	    else
	      remsetsed = setsed - wet_lag(j)
	      wet_lag(j) = 0.
	      if (wet_san(j) >= remsetsed) then
	        wet_san(j) = wet_san(j) - remsetsed
	        remsetsed = 0.
	      else
	        remsetsed = remsetsed - wet_san(j)
	        wet_san(j) = 0.
              if (wet_sag(j) >= remsetsed) then
	          wet_sag(j) = wet_sag(j) - remsetsed
	          remsetsed = 0.
	        else
	          remsetsed = remsetsed - wet_sag(j)
	          wet_sag(j) = 0.
                if (wet_sil(j) >= remsetsed) then
  	            wet_sil(j) = wet_sil(j) - remsetsed
	            remsetsed = 0.
	          else
	            remsetsed = remsetsed - wet_sil(j)
	            wet_sil(j) = 0.
                  if (wet_cla(j) >= remsetsed) then
	              wet_cla(j) = wet_cla(j) - remsetsed
	              remsetsed = 0.
	            else
	              remsetsed = remsetsed - wet_cla(j)
	              wet_cla(j) = 0.
	            end if
                end if
	        end if
	      end if
	    end if

          !! compute sediment leaving wetland
          wetsedo = wet_sed(j) * wetflwo

          wetsano = wet_san(j) * wetflwo
          wetsilo = wet_sil(j) * wetflwo
          wetclao = wet_cla(j) * wetflwo
          wetsago = wet_sag(j) * wetflwo
          wetlago = wet_lag(j) * wetflwo

          sedyld(j) = sedyld(j) + wetsedo

          sanyld(j) = sanyld(j) + wetsano
          silyld(j) = silyld(j) + wetsilo
          clayld(j) = clayld(j) + wetclao
          sagyld(j) = sagyld(j) + wetsago
          lagyld(j) = lagyld(j) + wetlago

          !! net change in amount of sediment in wetland for day
          wetsedc = vol * sed + wetsedi - wetsedo - wet_sed(j) *
     &                                                        wet_vol(j)
          !! determine settling rate for nutrients
          !! part of equation 29.1.3 in SWAT manual
          if (i_mo >= ipnd1(j) .and. i_mo <= ipnd2(j)) then
            iseas = 1
          else
            iseas = 2
          endif
          phosk = 0.
          nitrok = 0.
          phosk = psetlw(iseas,j) * wetsa * 10000. / wet_vol(j)
          phosk = Min(phosk, 1.)
          nitrok = nsetlw(iseas,j) * wetsa * 10000. / wet_vol(j)
          nitrok = Min(nitrok, 1.)

          !! remove nutrients by settling
          !! other part of equation 29.1.3 in SWAT manual
          wet_solp(j) = wet_solp(j) * (1. - phosk)
          wet_psed(j) = wet_psed(j) * (1. - phosk)
          wet_orgp(j) = wet_orgp(j) * (1. - phosk)
          wet_solpg(j) = wet_solpg(j) * (1. - phosk)
          wet_orgn(j) = wet_orgn(j) * (1. - nitrok)
          wet_no3(j) = wet_no3(j) * (1. - nitrok)
          wet_no3s(j) = wet_no3s(j) * (1. - nitrok)
          wet_no3g(j) = wet_no3g(j) * (1. - nitrok)

          if (wet_vol(j) < 1.e-6) wet_vol(j) = 0.0
          if (wet_orgn(j) < 1.e-6) wet_orgn(j) = 0.0
          if (wet_no3(j) < 1.e-6) wet_no3(j) = 0.0
          if (wet_no3s(j) < 1.e-6) wet_no3s(j) = 0.0
          if (wet_no3g(j) < 1.e-6) wet_no3g(j) = 0.0
          if (wet_orgp(j) < 1.e-6) wet_orgp(j) = 0.0
          if (wet_psed(j) < 1.e-6) wet_psed(j) = 0.0
          if (wet_solp(j) < 1.e-6) wet_solp(j) = 0.0
          if (wet_solpg(j) < 1.e-6) wet_solpg(j) = 0.0

          tpco = 0.
          tpco = 1.e+6 * (wet_solp(j) + wet_orgp(j) + wet_psed(j) +
     &                            wet_solpg(j)) / (wet_vol(j) + wetflwo)
          chlaco = 0.
          wet_chla(j) = 0.
          wet_seci(j) = 0.
          if (tpco > 1.e-4) then
            !! equation 29.1.6 in SWAT manual
            chlaco = chlaw(j) * 0.551 * (tpco**0.76)
            wet_chla(j) = chlaco * (wet_vol(j) + wetflwo) * 1.e-6
          endif
          if (chlaco > 1.e-4) then
            !! equation 29.1.8 in SWAT manual
            wet_seci(j) = secciw(j) * 6.35 * (chlaco**(-0.473))
          endif

          !! compute nutrients leaving wetland
          yy = 0.
          yy = wetflwo / (wet_vol(j) + wetflwo)
          sedorgn(j) = sedorgn(j) + wet_orgn(j) * yy / hru_ha(j)
          surqno3(j) = surqno3(j) + wet_no3(j) * yy / hru_ha(j)
          latno3(j) = latno3(j) + wet_no3s(j) * yy / hru_ha(j)
          no3gw(j) = no3gw(j) + wet_no3g(j) * yy / hru_ha(j)
          sedorgp(j) = sedorgp(j) + wet_orgp(j) * yy / hru_ha(j)
          sedminps(j) = sedminps(j) + wet_psed(j) * yy / hru_ha(j)
          surqsolp(j) = surqsolp(j) + wet_solp(j) * yy / hru_ha(j)
          minpgw(j) = minpgw(j) + wet_solpg(j) * yy / hru_ha(j)

          !!update nutrient pools in wetlands
          wet_orgn(j) = wet_orgn(j) * (1. - yy)
          wet_no3(j) = wet_no3(j) * (1. - yy)
          wet_no3s(j) = wet_no3s(j) * (1. - yy)
          wet_no3g(j) = wet_no3g(j) * (1. - yy)
          wet_orgp(j) = wet_orgp(j) * (1. - yy)
          wet_psed(j) = wet_psed(j) * (1. - yy)
          wet_solp(j) = wet_solp(j) * (1. - yy)
          wet_solpg(j) = wet_solpg(j) * (1. - yy)
          wet_chla(j) = wet_chla(j) * (1. - yy)

        end if

        !! add impoundment seepage to shallow aquifer
        !! This will be done in gwmod
!        shallst(j) = shallst(j) + wetsep / cnv

        !! compute seepage depth for HRU water balance
        twlwet(j) = wetsep / cnv

      end if

      if (qdr(j) < 0.) qdr(j) = 0.
      if (sedyld(j) < 0.) then
	    sedyld(j) = 0.0
        sanyld(j) = 0.0
        silyld(j) = 0.0
        clayld(j) = 0.0
        sagyld(j) = 0.0
        lagyld(j) = 0.0
	end if

      return
      end

max_stars_repo_name: GISWAT/erosion-sediment

max_issues_repo_name: GISWAT/erosion-sediment

max_forks_repo_name: GISWAT/erosion-sediment


text:
! { dg-do run }
      subroutine test_lower
      implicit none
      character(3), dimension(3) :: zsymel,zsymelr
      common /xx/ zsymel, zsymelr
      integer :: znsymelr
      zsymel = (/ 'X', 'Y', ' ' /)
      zsymelr= (/ 'X', 'Y', ' ' /)
      znsymelr=2
      call check_zsymel(zsymel,zsymelr,znsymelr)

      contains

      subroutine check_zsymel(zsymel,zsymelr,znsymelr)
        implicit none
        integer znsymelr, isym
        character(*) zsymel(*),zsymelr(*)
        character(len=80) buf
        zsymel(3)(lenstr(zsymel(3))+1:)='X'
        write (buf,10) (trim(zsymelr(isym)),isym=1,znsymelr)
10      format(3(a,:,','))
        if (trim(buf) /= 'X,Y') STOP 1
      end subroutine check_zsymel

      function lenstr(s)
        character(len=*),intent(in) :: s
        integer :: lenstr
        if (len_trim(s) /= 0) STOP 2
        lenstr = len_trim(s)
      end function lenstr

      end subroutine test_lower

      subroutine test_upper
      implicit none
      character(3), dimension(3) :: zsymel,zsymelr
      common /xx/ zsymel, zsymelr
      integer :: znsymelr
      zsymel = (/ 'X', 'Y', ' ' /)
      zsymelr= (/ 'X', 'Y', ' ' /)
      znsymelr=2
      call check_zsymel(zsymel,zsymelr,znsymelr)

      contains

      subroutine check_zsymel(zsymel,zsymelr,znsymelr)
        implicit none
        integer znsymelr, isym
        character(*) zsymel(*),zsymelr(*)
        character(len=80) buf
        zsymel(3)(:lenstr(zsymel(3))+1)='X'
        write (buf,20) (trim(zsymelr(isym)),isym=1,znsymelr)
20      format(3(a,:,','))
        if (trim(buf) /= 'X,Y') STOP 3
      end subroutine check_zsymel

      function lenstr(s)
        character(len=*),intent(in) :: s
        integer :: lenstr
        if (len_trim(s) /= 0) STOP 4
        lenstr = len_trim(s)
      end function lenstr

      end subroutine test_upper

      program test
        call test_lower
        call test_upper
      end program test


text:
The Davis Enforcer is an Campus oncampus student publication for the Christian Democratic Party, and is now part of AS Papers.  Two issues have been released thus far, and both were printed at the Davis Enterprise. It is a four page newspaper meant to detail mainly ASUCD news and information. Articles are written by staff and various other members of the campus community, some anonymously. It is printed in black and white and currently free of advertisements.  An attempt was being made to hire people to find advertisers.  It has been widely ridiculed as ASUCD election propaganda, and issues typically contain many spelling errors and grammatical mistakes.

The publication was founded by a $750 grant to Steve Ostrowski by http://www.leadershipinstitute.org/ The Leadership Institute: Conservative Training for Tomorrows Leaders, a national organization whose mission is to identify, recruit, train, and place conservatives in politics, government, and media. The Institute sent a packet in the mail to the staff on how to make noise and how to take over ASUCD. Most ASUCD insiders and observers believe that this noisemaking is the entire mission of the paper.
   LI is back lol  GregWebb

 Distribution

The first issue was distributed on on November 15th 2006, the same day as the Fall ASUCD election. The first edition was intentionally placed in AS Papers and The California Aggie Aggie newsstands, violating various regulations as the publication was not affiliated with AS Papers at the time.

On Tuesday February 13, 2007, three months after the initial publication, the Enforcer published their second issue. This time around, the publication affiliated with AS PAPERS and utilized their distribution infrastructure legitimately.

The Davis Enforcer also used Internet distribution through Facebook ads (only around ASUCD elections, of course).

 ASUCD Court Case
The paper was the defendant in the ASUCD Court Case #41, Laabs v. AS PAPERs (Davis Enforcer). The second issue contained glowing articles about certain ASUCD candidates and not others. Brent Laabs alleged, inter alia, that these stories were tantamount to endorsements. While AS PAPERS publications are free to make political endorsements, they are required to invite each candidate running to an interview, which they did not do. The Enforcer, of course, argued that the stories did not constitute endorsements. The ASUCD Court found that while the stories were unquestionably endorsements in the usual sense and that the DE intended to make endorsements, they were not endorsements in the context of ASUCD. The ASUCD bylaws require that certain forms be filled out in order for the endorsement to be official. No such forms were ever filed by the Enforcer and they won the case on what most considered to be a technicality. After the case, a committee was formed to close the loophole that gave the Enforcer their victory. No proposed changes were made due to free speech concerns.

Staff

  Ronald Payne  Chief Editor, Liaison to the Leadership Institute.
  Steve Ostrowski  Chief Editor, Liaison to  AS Papers.

20061115 22:08:24 nbsp All the Davis Enforcer enforced was my disdain for individuals who cant write.  If youre going to publish propaganda, try writing it about something people are actually concerned with.  You didnt even spell the names of Student Focus candidates correctly.  What a shining endorsement. Users/JillWeinstein

20061116 10:35:05 nbsp The only effect the Davis Enforcer will have will be adding more weight to the recycling bins.  Users/JamesSchwab

20061116 13:23:05 nbsp I picked one up yesterday (out of the ASPAPERs rack at the MU) and let me just say that this is nothing more than a soft money donation to the campaigns of Student Focus and Steve Ostrowski ordained independent candidates. Steve, how come you didnt go to the the Green party looking for writers. Oh wait, because this is propaganda... the Green party tries to stick to the inconvenient truths of the world and politics, which is probably why we lose so often. Anyway the Enforcer is nothing but fluff. Heck, I think Ill pick up a few extra copies and stuff a pillow with them  that is if I wasnt worried that any of the ridiculous ideas would seep into my head through osmosis. Steve, I advise you to stop trying to help with people campaigns and instead simply pray. Let Jesus determine the winner. Jesus certainly does not believe in the separation of church and state. Oh wait, as Christ said in Matthew 15:21 Give to Caesar what is Caesars, and to God what is Gods. Therefore allowing for a separation between governments and the divine. So that means leave God out of the pledge of allegiance. And on that matter, the pledge of allegiance violates the second commandment anyways so I guess as Evangelicals, youre damned if you do and damned if you dont. Ill stop rambling. Users/RobRoy

20061117 03:20:26 nbsp I love propaganda. Writing it is fun, and reading it can be hillarious. Why? Because it doesnt stand up to the truth. And, Steve, sorry bout that last night, and if ya want to, you can remove this comment. Users/MaxMikalonis

20061118 00:21:01 nbsp The only problem was that this was written so poorly.  I think that all laughter was largely a product of terrible composition...and certainly not reflective of quality propaganda.  Also, I think about 20 people on campus actually read it. Users/JillWeinstein

20061118 22:27:46 nbsp The humor in the enforcer is rather subtle, but quite ROTFL inducing. It appears to be geared more towards the ASUCD savy consumers, however. I suggest that the Supreme Editors tone down the amount of insider language so that non student government people know what its talking about, or revamp it into a proLEADing Focus publication. Users/JosephBleckman

20061215 11:55:36 nbsp In the PDF version available on this page Steve has illegal reprinted photos of Rob Roy and Dan Xie.  Users/PaulHarms
   True, interestingly the quoted material from the Wiki is properly attributed.  All that has to be done is put a Photo from DavisWiki as a photocredit caption and hed be set (at least the Rob Roy one, not sure about the other).  Users/JabberWokky jw

20070125 22:19:38 nbsp Is this a joke? I havent seen the publication on campus. Users/TusharRawat

20070125 23:15:40 nbsp We are involved in AS Papers new advertisement deal.   Users/SteveOstrowski

20070210 20:26:40 nbsp Its real in that it had one briefly released volume. I had to get one through Joseph though. Users/CarlosOverstreet

20070212 23:30:30 nbsp I liked the Suyat and Peake quotations, I pray that they didnt really say that stuff. Users/GregWebb

20070212 23:34:23 nbsp How come the Aggie Spoof is not called the Davis Enforcer? Users/JamesSchwab

20070213 07:20:17 nbsp i was thinking the same thing....they should really think about calling it that this year Users/JillWeinstein

20070220 16:42:21 nbsp I picked up issue 2 without knowing anything about the Enforcer, and found that Ostrowski is both the subject of the frontpage article and the chief editor. Theres not much I can add to that. Oh, wait. Make that TWO front page articles that mention Ostrowski.Users/GeoffJohnson
   One has to make news to report on it.  Users/SteveOstrowski
   Actually, I have to correct myself here. Its THREE articles. Steve, dont you see writing about yourself as a conflict of interest? Users/GeoffJohnson
     Nah, seeing as how I am not paid.  Users/SteveOstrowski

20070220 16:55:32 nbsp News is something relevant and has an affect on ones life. Ostrowski is not relevant nor does he have an affect on anyones life. Users/JamesSchwab

20070220 16:58:38 nbsp But theres no indication, save for the bottom of the corner of the back page, that you wrote this. No author is attributed. If people cant trust your journalism, how can they trust your politics?   Let me rephrase that. What youve written is clearly meant to appear to the untrained eye as a newspaper, but by definition, it cant be objective if the author is the subject. Thats subversive and unfair to the people whose vote you are asking for. Users/GeoffJohnson
    I am not asking for anyone to vote for anyone or for anything.  The authors are left off because I have no intention of getting others sued in ASUCD Court.  Users/SteveOstrowski

20070220 17:09:43 nbsp I dont know what you guys are talking about, the Davis Enforcer is as fine as a newspaper publication as the http://www.weeklyworldnews.com/ Weekly World News or other such publications.  People walk by, pick it up, flip through a pages, read the news, have a giggle, and walk away.  Long live both!  (And besides, the Comic Press News always runs out quickly, and I think this paper helps.)  Users/EdwinS ES

20070220 17:12:44 nbsp Its good you dont look into the other AS PAPERS publications.   Users/SteveOstrowski

20070220 19:54:28 nbsp none of the current Davis Enforcer members have any experience in journalism or Englishsic which may hinder the affectivenesssic of the first issue.  So, I guess journalism is yet another field that evangelicals consider expertise to be a hinderance. Users/JesseSingh

20070220 19:58:58 nbsp First of all none of us are evangelicals, in fact most of us are Catholic. And most AS PAPERS publications are not journalism students.  Users/SteveOstrowski

20070220 20:59:27 nbsp The lawsuit only had a minor effect on AS PAPERs decision not to distribute the Davis Enforcer.  The main reason is that the Enforcer was already on probation for failing to volunteer enough at the office. Users/BrentLaabs

20070220 21:10:51 nbsp There is a conflict of interests though.  The Davis Enforcer is under the editorial control of the Christian Slate, and of Steve Ostrowski.  The Christian Slate is under the control of Steve Ostrowski.  The Christian Slate has repeatedly said (on fliers and facebook) that their webpage is this Davis Wiki one.   Interestingly enough, the Christian Democratic Party page does list endorsements on their so seeming official site.  And the link between the two is obviously very clear and stated on multiple pages, much less its all in Steve Os hands.  And the two that the Christian Democratic Party endorsed coincidentally had a segment about them in the issue (the only prez/VP combo with a segment!).  Just reading the rest of the issue makes it obvious to everyone its not exactly objective reporting of news.   Since you didnt use the word endorse, then obviously its not one. >.>   I guess this court case is a huge thing in letter of the law vs the spirit of it.    Users/EdwinS ES
  Actually, the letter of the law supports this one too.  Since the endorsement clause is the only thing that allows AS PAPERs to bypass the requirement that an ASUCD unit may not make a public statement on elections, he still has to follow that clause whether he makes endorsements or not.   Users/BrentLaabs
  Is this really necessary, you are only slowing me down a little and giving me a lot of information to send to other places.  Isnt this more helpful to me than to you?  Users/SteveOstrowski
   You keep on hurting someone I care about.  There is nothing I wont do to stop you from doing that, and I will work against you in any endeavor no matter how much it hurts me.    If you keep destroying those around you, then I will destroy you.   Users/BrentLaabs
  Well said, Brent! Users/KevinPowers
   Never assume that the information for the DE is only comming from me.  Look around in ASUCD and there will always be someone who is assisting me with information. As for Pham it was accurate without a doubt.  Users/SteveOstrowski

20070220 22:32:20 nbsp If the paper were another vessel of LEAD (even a puppet conservative paper), I doubt half of you would have anything to say about it... Users/CarlosOverstreet
   Youre right.  Because no one from LEAD would have violated the bylaws (which were written before LEAD existed, by the way).  They would have done the interviews, because they are professional.  And they wouldnt use a paper with ambiguous authorship and minimal production values to spread lies about their opponents.  Hell, Student Focus wouldnt do that either.  Only the Christian Slate has so little honor.   Users/BrentLaabs

20070220 22:44:17 nbsp Brent... are you stealing DE papers? Users/CarlosOverstreet
  No, because SJA has specifically told me not to take multiple papers out of the rack and I fear the wrath of Don Dudley.  See California Aggie Printing for more.  Users/BrentLaabs
  A small percentage of Enforcers are going through AS PAPERS on a daily basis to catch thieves, dont try it.  Users/SteveOstrowski
   Well, thanks for the admission of guilt!   Users/BrentLaabs
   It was always anticipated that LEAD would try to steal some of them.   Users/SteveOstrowski
    So wait, as small percentage were placed in stands, and they somehow disappeared?  zOMG, someone might be taking your paper and reading it!!!!!11!!   Users/BrentLaabs
    Not when all of the racks are empty in multiple places.  I will find out soon enough, but this has happened all over the country and thus I am in suspicious mode.  Users/SteveOstrowski

20070220 23:32:41 nbsp I think everyone here, especially Steve, needs to stop being such little kids fighting over nothing. I seriously think we should just delete this page since it is such a source of argument on the wiki and detracts from other pages getting viewed or editted as all the editors seem to be focused on this insignificant battle here. I also think that if the creator of this publication wanted it to be taken seriously he would not have given it the name enforcer (whos gonna read something like that?), he wouldnt have made its symbol the Imperial emblem from Star Wars (possible copyright infringement?), and he wouldnt have put a crybaby that looks for conspiracy under every bush in charge as the Editor in Chief. Users/JamesHaile
  Uhmm no because deletion of a page that has info on it is bad.  Its a real thing at Davis and the reason there are so many edits is because there seems to be some resistance from some ASUCDers.  Users/SteveOstrowski

20070220 23:37:23 nbsp Also, where can I get a copy of the Enforcer? Users/JamesHaile

     It should be available for download at the PDF link above. Users/CarlosOverstreet

20070221 17:40:46 nbsp I love how the enforcer logo is a star wars empire logo spinoff, it fits nicely with chief editors aim of empirestv... Users/StevenDaubert
   Spinoff? Im pretty sure its just the wiki:WikiPedia:Galactic_Empire_(Star_Wars) Galactic Empires logo with the letters D E MSPainted into it. Users/TravisGrathwell
     as corny as MSPainting DE is, its still a change, and not the origional logo... Hence why I used spinoff. Users/ThomasLloyd has offered his un lawer semi educated oppinion below and I will take his statement at face value.

20070223 10:26:35 nbsp Steve!  You should have told me you were running an article on the ASUCD Supreme Court!  I would have loved to give an interview. Users/TimCoady

20070223 11:24:12 nbsp I dont pretend to be a real lawyer, but Ive had enough Intro to IP to say that the use of the above logo is probably not a sufficient transformative use to avoid copyright enforcement. Users/ThomasLloyd
   It is probably not legal for DEs use, but since the Wiki is about DE, not actually DE itself, it is an illustration of a logo rather than a logo usage.  I.e., we can have a screen shot of UCDavisLife.com with the Aggie logo as an example of what they were before they became UDavisLife.com, even though it has been asserted UCDavisLife and Aggie Logo was an illegal use of the UC Davis trademark (as opposed to copyright, which that entry erroneously asserts). Users/JabberWokky

20070223 18:36:34 nbsp I saw the enforcer in the silo today... It had a full page Christian Slate endorsment! Turns out there man running is also an editor! I speed read some of the article and then threw it in the trash. In hindsight I should have used the recyle next to the trash. Users/StevenDaubert

20070502 19:03:46 nbsp You people are assaulting the Enforcer for spelling issues and lack of journalism experience of the editor(s). Most of the attacks on the Davis Enforcer I’ve read on this page are simpleminded and dismissable, despite how intricate these people are trying to sound with their arguments. Your attacks are pathetic only because they address that which is well outside of the crucial point that the Enforcer is trying to convey. Your attention to irrelevant details illuminates the way you draw attention away from what is actually being said in the Enforcer. It seems like people care a lot about who YOU are, Steve, and how they can label YOU as unimportant, so that less people read the Enforcer. That is a classic defense, using the high level rhetoric that you can find by scrolling up a bit, instead of putting your point simply, the way a real, goodnatured and educated human being would do.

Bashing the Enforcer the way you are makes you a disgrace not only to our campus, but to supporters of free speech everywhere. It is as much propaganda as the issue of the California Aggie that i read today, May 2, 2007, which quoted officials criticizing DCR members for their capture the flag game on the QUAD, labeling it as racist and disrespectful while including NO derogatory quotes regarding the protestors down the block that made countless students late for class and hundreds ordinary citizens frustrated. Lets also disregard that the same person wrote both articles. Disrespect? Your attacks on the Enforcer show your hate and fear for people just like Steve that you so deeply conceal. And for good reason, because what he writes reeks of truth and reality, and it makes your blood boil how much of a wakeup call it is, unlike the sugar coated crap you are fed for breakfast every morning. If you read thus far, thank you for your time. I can only hope that there are more of you out there that see this the way I do.

And another thing....that logo is probably reminding you of the dark side for a reason. Steve, touchè. Users/boristheblade

20070509 18:59:32 nbsp I find this publication extra difficult to read knowing the Chief Editor and some of his thoughts. In regards to an article printed in the California Aggie about an anonymous, racist midterm turned into an AfricanAmerican studies course, the Chief Editor wrote If you want to destroy an idea you need to silence them in the press, find those who did it, and disappear them.  (Free Speech! Right Boris?)  Seeing that ideology expressed by a supposed newspaper, the Chief Editor and founder no less, makes it difficult for many, myself included, to attempt to take this publication as anything other than a skewed ASUCD political propaganda.  Users/EdWins
   My blog entries are of course separate from the Davis Enforcer and I am correct about that midterm.  Because it is out in the open it gave that racist idiot a forum and it gave him publicity he didnt deserve.  If you want to play the media game right you silence extremists.   Users/SteveOstrowski
     And who gets to decide whos extreme enough to silence?  Theres a difference between silencing extremists (whatever that means), and as you suggested, a newspaper deliberately not reporting on a story because it makes UC Davis look bad! This makes UC Davis and the city of Davis look bad. Because of one person our reputation is going to be ruined because this had to go to the press. Why should the entire community pay for the idiocy of one person.  Yeah, newspapers shouldnt report any news or issues if it only involves one person. Lets keep all that stuff quiet, shred it.  Then again, maybe one day newspapers will http://www.spiegel.de/international/0,1518,479023,00.html have to figure out the balance of good and bad.  Users/EdwinS ES

20070510 02:34:51 nbsp I demand a copy of the first Davis Enforcer.   Users/GregWebb

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

!
!  *********************************************************************
!  *   SUBROUTINE INPUT1 TO READ TITLE AND SOME LAKE-SPECIFIC INPUTS   *
!  *********************************************************************
!
!  0895 MvdV
!  read time series biomass for several grazer types
!
      SUBROUTINE INPUT1 (NDEC,INPU)
      IMPLICIT REAL*8 (A-H,O-Z)
      INCLUDE 'blmdim.inc'
      INCLUDE 'putin1.inc'
      INCLUDE 'putin2.inc'
!
!  Note: this subroutine is NOT used by the coupled BLOOM II/CHARON
!  program.
!
!  Read year number and comments for particular run.
!
      READ (INPU,99999) IYEAR,(CASE(I),I=1,9)
      READ (INPU,99998) COM
!
!  Read number of weeks in the input data set
!
      READ (INPU,99997) NDEC
!
!  Read data for all weeks of the year:
!  Read: DAYLGT--day length; DEATH--death rate for weeks of year
!  Read: zooplankton dryweigth for weeks; depth for weeks.
!
!  0895 MvdV
!  read time series biomass for several grazer types
!
      DO 10 I=1,NDEC
      READ (INPU,99996) DAYLGT(I),DEATH(I),ZOOD(0,I),DEPTH(I),
     1     (ZOOD(J,I),J=1,MG)
   10 CONTINUE

99999 FORMAT (I4,1X,9A8)
99998 FORMAT (9A8,8X)
99997 FORMAT (I5,10X,F5.2)
99996 FORMAT (40F10.2)
      RETURN
      END

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


 program ifs

      implicit none

!     Program to demonstate if statements.

!     Declarations
      integer*4 i, j
      integer*4 ierr    ! Error when reading from screen
      logical result
      character*80 line

!**   Use some simple input to get value
      line = 'Program to demonsrate if statements'
      write(*,'(a)') line    ! Demonstrates another way to output

!***  Ask the user for a value
      write(*,'(a)') 'Enter two values between 1 and 10'
      read(*,*) i,j

!***  Output values from user
      write(*,100) i,j
 100  format('User input values ',I3,' and ',i3)

!     See if error generated on read (! used for comment, so /- for NE)
      if( ierr /= 0 ) then   ! If error is not zero
         write(*,200) ierr
 200     format('IOSTAT Error ',i5,' occurred reading i and j')
         print *,'Values ',i,j
!        If this is servere we could uncomment line below
!        stop 'Error reading i and j'
      endif

!***  Now test the values.  Notice the indenting (not needed but
!     good to allow easy reading of code).
      if( i.gt.10 .and. j.gt.10 ) then
          print *,'User gave values too large'
      elseif( i.lt.1 .and. j.le.0 ) then
!         Notice in line above that test is the same for i and j
!         (because they are integer)
          print *,'User gave values too small'
      elseif( i.lt.1 .or. i.gt. 10 ) then
          print *,'User gave invalid value of I'
      elseif( j.le.0 .or. j.ge. 11 ) then
          print *,'User gave invalid value of J'
      else
          write(*,'(a)') 'Users input was good'
          if( i.gt.j ) then
             write(*,310) i,i-j
 310         format('I=',i3,' is greater than J. Difference is ',i3)
          else
             write(*,320) J,i-j
 320         format('J=',i3,' is greater than or equal I. Difference is ',i3)
          endif
      endif

 end


text:
	SUBROUTINE TG_DUAL  ( time, clntm, iret )
C************************************************************************
C* TG_DUAL								*
C*									*
C* This subroutine converts a time(2) array into a colon separated dual	*
C* grid time stamp.							*
C*									*
C* TG_DUAL  ( TIME, CLNTM, IRET )					*
C*									*
C* Input parameters:							*
C*	TIME (2)	CHAR*		Grid time array			*
C*									*
C* Output parameters:							*
C*	CLNTM		CHAR*36		Colon separated dual grid time	*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C**									*
C* Log:									*
C* K. Brill/HPC		 2/04						*
C************************************************************************
	CHARACTER*(*)	time (*), clntm
C*
C------------------------------------------------------------------------
	iret  = 0
C*
	CALL ST_LSTR ( time (1), lt1, ier )
	IF ( time (2) .ne. ' ' ) THEN
	    clntm = time (1)(1:lt1) // ':' // time (2)
	ELSE
	    clntm = time (1)(1:lt1)
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
!Autor: Claudio Iván Esparza Castañeda
!Título: Pi
!Descripción: Calcula Pi por medio de la serie de Leibinz
!Fecha: 15/03/2020

Program CPi !Inicio del programa
  implicit none !Quitar variables implícitas
  REAL*16::pi, s !Declaración de variables reales de 16 bits
  INTEGER::i !Variable entera iterativa
  s=0.0 !inicializar memoria auxiliar
  do i=0, 1000000000 !Ciclo que contiene a la serie
     s=(-1.0)**i/(2.0*i+1.0)+s !Forma matemática de la serie
  end do
  pi=4.0*s !Valor de pi igual a cuatro veces el valor de la serie
  write(*, *) pi !Imprimir pi
end Program CPi !Finalizar programa

max_stars_repo_name: clivan/Fortran

max_issues_repo_name: clivan/Fortran

max_forks_repo_name: clivan/Fortran