text:
C$Procedure  ZZEKTRES ( Private: EK, resolve times in encoded query )

      SUBROUTINE ZZEKTRES ( QUERY,  EQRYI,  EQRYC,   EQRYD,
     .                              ERROR,  ERRMSG,  ERRPTR  )

C$ Abstract
C
C     SPICE Private routine intended solely for the support of SPICE
C     routines.  Users should not call this routine directly due
C     to the volatile nature of this routine.
C
C     Resolve time values in an encoded EK query.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     EK
C
C$ Keywords
C
C     EK
C     PRIVATE
C
C$ Declarations

      INCLUDE 'ekattdsc.inc'
      INCLUDE 'ekbool.inc'
      INCLUDE 'ekcnamsz.inc'
      INCLUDE 'ekopcd.inc'
      INCLUDE 'ekqlimit.inc'
      INCLUDE 'ekquery.inc'
      INCLUDE 'ektnamsz.inc'
      INCLUDE 'ektype.inc'

      INTEGER               LBCELL
      PARAMETER           ( LBCELL = -5 )

      CHARACTER*(*)         QUERY
      INTEGER               EQRYI  ( LBCELL : * )
      CHARACTER*(*)         EQRYC
      DOUBLE PRECISION      EQRYD  ( * )
      LOGICAL               ERROR
      CHARACTER*(*)         ERRMSG
      INTEGER               ERRPTR

C$ Brief_I/O
C
C     Variable  I/O  Description
C     --------  ---  --------------------------------------------------
C     EQRYI     I-O  Integer component of query.
C     EQRYC     I-O  Character component of query.
C     EQRYD     I-O  Numeric component of query.
C     ERROR      O   Error flag.
C     ERRMSG     O   Error message.
C     ERRPTR     O   Position in query where error was detected.
C
C$ Detailed_Input
C
C     QUERY          is the original query from which EQRYI and EQRYC
C                    were obtained.  QUERY is used only for
C                    construction of error messages.
C
C     EQRYI          is the integer portion of an encoded EK query.
C                    The query must have been parsed.
C
C     EQRYC          is the character portion of an encoded EK query.
C
C     EQRYD          is the numeric portion of an encoded EK query.
C
C$ Detailed_Output
C
C     EQRYI          is the integer portion of an encoded EK query.
C                    On output, all valid time values will have been
C                    converted from strings to equivalent numeric
C                    values which represent times as ephemeris
C                    seconds past J2000 (TDB).
C
C     EQRYC          is the character portion of an encoded EK query.
C
C     ERROR          is a logical flag indicating whether an error was
C                    detected.  Note that a time string might be
C                    syntactically valid, but incapable of being
C                    converted to ET if the appropriate time kernels
C                    (Leapseconds or SCLK) are not loaded.
C
C     ERRMSG         is an error message describing an error in the
C                    input query, if one was detected.  If ERROR is
C                    returned .FALSE., then ERRPTR is undefined.
C
C     ERRPTR         is the character position in the original query
C                    at which an error was detected, if an error was
C                    found.  This index refers to the offending lexeme's
C                    position in the original query represented by the
C                    input encoded query.  If ERROR is returned .FALSE.,
C                    ERRPTR is undefined.
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     1)  If the input query is not initialized, the error will be
C         diagnosed by routines called by this routine.  The outputs
C         will not be modified.
C
C     2)  If names have not been resolved in the input query, the error
C         SPICE(NAMESNOTRESOLVED) will be signalled.  The outputs
C         will not be modified.
C
C     3)  If any sort of time conversion error occurs, the output flag
C         ERROR is set, and an error message is returned.
C
C$ Files
C
C     None.
C
C$ Particulars
C
C     Strings representing time values are interpreted as follows:
C
C        1)  The string is first examined to see whether it's an
C            SCLK string for a recognized clock; if it is, the
C            string is converted to the equivalent ET.
C
C        2)  If the string is not a SCLK string, it is expected
C            to be some sort of UTC representation.  The string is
C            checked to see whether it's an ISO format UTC time that
C            ISO2UTC can handle.
C
C        3)  If the string does not conform to an ISO format, the
C            last chance is to try to get the string through
C            TPARSE.  If TPARSE can't deal with it, it's considered
C            to be invalid.
C
C
C     This routine assumes that encoded EK query architecture version
C     1 is to be used with the query to be initialized; this routine
C     will not work with any other architecture version.
C
C$ Examples
C
C     See EKFIND.
C
C$ Restrictions
C
C     1) A leapseconds kernel must be loaded at the time this routine
C        is called.
C
C     2) In order to convert SCLK strings, an appropriate SCLK kernel
C        must be loaded at the time this routine is called.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     N.J. Bachman       (JPL)
C
C$ Version
C
C-    Beta Version 1.0.0, 11-OCT-1995 (NJB)
C
C-&

C
C     SPICELIB functions
C
      LOGICAL               FAILED

C
C     Local parameters
C
      INTEGER               SHORT
      PARAMETER           ( SHORT  = 32 )

C
C     Local variables
C
      CHARACTER*(TNAMSZ)    ALIAS  ( MAXTAB )
      CHARACTER*(CNAMSZ)    COLNAM
      CHARACTER*(TNAMSZ)    TABLE  ( MAXTAB )
      CHARACTER*(SHORT)     TIMSTR
      CHARACTER*(1)         TOUCHC

      DOUBLE PRECISION      ET

      INTEGER               ATTDSC ( ADSCSZ )
      INTEGER               BASE
      INTEGER               CNSTYP
      INTEGER               COLIDX
      INTEGER               DESCR  ( EQVDSZ )
      INTEGER               DTYPE
      INTEGER               I
      INTEGER               IRSOLV
      INTEGER               LXB
      INTEGER               LXE
      INTEGER               NCNS
      INTEGER               NTAB
      INTEGER               OPCODE
      INTEGER               SB
      INTEGER               SE
      INTEGER               TABIDX

C
C     No error to start with.
C
      ERROR  =  .FALSE.
      ERRMSG =  ' '
      ERRPTR =  0

      TOUCHC =  QUERY(1:1)

C
C     The query must have had names resolved at this point, or it's no
C     go.
C
      CALL ZZEKREQI ( EQRYI, 'NAMES_RESOLVED', IRSOLV )

      IF ( FAILED() ) THEN
         RETURN
      END IF

      IF ( IRSOLV .EQ. IFALSE ) THEN

         CALL CHKIN  ( 'ZZEKTRES'                                  )
         CALL SETMSG ( 'Names are not resolved in encoded query.'  )
         CALL SIGERR ( 'SPICE(NAMESNOTRESOLVED)'                   )
         CALL CHKOUT ( 'ZZEKTRES'                                  )
         RETURN

      END IF

C
C     Time strings occur only on the right sides of constraints.
C     Examine each constraint that compares a column and a value.
C
      CALL ZZEKREQI ( EQRYI, 'NUM_TABLES',      NTAB )
      CALL ZZEKREQI ( EQRYI, 'NUM_CONSTRAINTS', NCNS )

      DO I = 1, NCNS
C
C        Calculate the base address of the constraint.
C
         BASE    =  EQVBAS + (NTAB * 2 * EQVDSZ) + (I-1)*EQCDSZ

C
C        Obtain the constraint type.  If the RHS is not a value or if
C        the RHS is null (as indicated by the opcode), we can skip it.
C
         CNSTYP  =  EQRYI ( BASE + EQCTYP )
         OPCODE  =  EQRYI ( BASE + EQOPCD )

         IF (        ( CNSTYP .EQ. EQVAL  )
     .        .AND.  ( OPCODE .NE. ISNULL )
     .        .AND.  ( OPCODE .NE. NOTNUL )  ) THEN
C
C           Get the index of the table containing the LHS column, and
C           get the index of this column within that table.  Get the
C           table name, then get the column's attributes.
C
            TABIDX  =  EQRYI ( BASE + EQLTAB - 1 + EQTORD )
            COLIDX  =  EQRYI ( BASE + EQLCOL - 1 + EQCIDX )

            CALL ZZEKQTAB ( EQRYI, EQRYC, TABIDX, TABLE, ALIAS )

            CALL EKCII ( TABLE, COLIDX, COLNAM, ATTDSC )

            DTYPE  =  ATTDSC(2)

            IF ( DTYPE .EQ. TIME ) THEN
C
C              The RHS points to a string representing a time
C              value.
C
               LXB    =  EQRYI ( BASE + EQBVAL - 1 + EQBLEX )
               LXE    =  EQRYI ( BASE + EQBVAL - 1 + EQBLEX )
               SB     =  EQRYI ( BASE + EQBVAL - 1 + EQBSTR )
               SE     =  EQRYI ( BASE + EQBVAL - 1 + EQESTR )
               TIMSTR =  EQRYC ( SB:SE )

C
C              Convert the time to ET, if possible.
C
               CALL ZZEKTCNV ( TIMSTR, ET, ERROR, ERRMSG )

               IF ( ERROR ) THEN

                  ERRPTR = SB
                  RETURN

               END IF
C
C              Insert the ET value into the query, and replace the
C              value descriptor for the time string.
C
               CALL ZZEKINQN ( ET, TIME, LXB, LXE, EQRYI, EQRYD, DESCR )

               CALL MOVEI ( DESCR, EQVDSZ, EQRYI(BASE+EQBVAL) )

            END IF
C
C           We've parsed a time string, if the current column's type
C           was TIME.
C

         END IF
C
C        We've examined the current constraint, if it compares a
C        column with a value.
C

      END DO

C
C     Indicate completion of time resolution.
C
      CALL ZZEKWEQI (  'TIMES_RESOLVED',  ITRUE,  EQRYI  )

      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
C
C     ZOEK IS NIET AFHANKELIJK VAN UPPERCASE/LOWERCASE
C
C     De routine maakt gebruik van ICHAR()
C     a t/m z hebben codes 97 t/m 122
C     A t/m Z hebben codes 65 t/m 90
C
      SUBROUTINE ZOEK (NAAM,NOTOT,SYNAME,NZOEK,INDEX)
      INTEGER NOTOT, NZOEK, INDEX
      CHARACTER*(*) NAAM
      CHARACTER*(*) SYNAME(NOTOT)
      INTEGER I,K,I1,I2
      INDEX = -1
c     WRITE (*,'(A)') NAAM(1:NZOEK)
      DO 100 I = 1,NOTOT
c         WRITE (*,'(I5,A)') I,SYNAME(I)(1:NZOEK)
          DO 50 K = 1,NZOEK
              I1 = ICHAR(NAAM     (K:K))
              I2 = ICHAR(SYNAME(I)(K:K))
              IF (I1.GE. 97.AND.I1.LE.122) THEN
                  IF (I1.NE.I2.AND.I1.NE.(I2+32)) GOTO 100
              ELSEIF (I1.GE. 65.AND.I1.LE. 90) THEN
                  IF (I1.NE.I2.AND.I1.NE.(I2-32)) GOTO 100
              ELSE
                  IF (I1.NE.I2                  ) GOTO 100
              ENDIF
   50     CONTINUE
          INDEX = I
          GOTO 200
C         IF (NAAM(1:NZOEK).EQ.SYNAME(I)(1:NZOEK)) GOTO 200
  100 CONTINUE
c     WRITE (*,*) ' ZOEK:',INDEX
      RETURN
  200 CONTINUE
c     WRITE (*,*) ' ZOEK:',INDEX
      RETURN
      END

ext: f

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
subroutine nsi_reaous(a,itask)
   use typre
   use Mod_iofile
   use Mod_NavierStokes
   implicit none

   class(NavierStokesProblem) :: a
   integer(ip) :: itask

   !For Listener
   real(rp), pointer     :: param(:) => NULL()
   character(5), pointer :: words(:) => NULL()
   integer(ip), pointer  :: nnpar => NULL(),nnwor => NULL()
   integer(ip) :: ndime,ipara

   !Initializations

   call a%Listener%getarrs(words,param,nnpar,nnwor)

   if (itask == 0) then         !Initializations.
      a%kfl_dispa = 0
      a%kfl_postBtract = 0
      a%kfl_outfm = 0

   elseif (itask == 1) then

      ! Step or time for post-process
      if(words(1)=='POSTP') then
         if(words(2)=='VELOC') then
            if(words(3)=='STEPS') then
               a%npp_stepi(1) = a%Listener%getint('STEPS',1,'#Postprocess step interval for velocities')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,1)=param(4:13)
            end if
         else if(words(2)=='PRESS') then
            if(words(3)=='STEPS') then
               a%npp_stepi(2) = a%Listener%getint('STEPS',1,'#Postprocess step interval for pressure')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,2)=param(4:13)
            end if
         else if(words(2)=='STREA') then
            if(words(3)=='STEPS') then
               a%npp_stepi(3) = a%Listener%getint('STEPS',1,'#Postprocess step interval for streamlines')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,3)=param(4:13)
            end if
         else if(words(2)=='DENSI') then
            if(words(3)=='STEPS') then
               a%npp_stepi(4) = a%Listener%getint('STEPS',1,'#Postprocess step interval for Density')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,4)=param(4:13)
            end if
         else if(words(2)=='VISCO') then
            if(words(3)=='STEPS') then
               a%npp_stepi(5) = a%Listener%getint('STEPS',1,'#Postprocess step interval for Viscosity')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,5)=param(4:13)
            end if
         else if(words(2)=='STRES') then
            if(words(3)=='STEPS') then
               a%npp_stepi(6) = a%Listener%getint('STEPS',1,'#Postprocess step interval for Stress')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,6)=param(4:13)
            end if
         else if(words(2)=='TURBU') then
            if(words(3)=='STEPS') then
               a%npp_stepi(7) = a%Listener%getint('STEPS',1,'#Postprocess step interval for Turbulent viscosity')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,7)=param(4:13)
            end if
         else if(words(2)=='PERME') then
            if(words(3)=='STEPS') then
               a%npp_stepi(8) = a%Listener%getint('STEPS',1,'#Postprocess step interval for Permeability')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,8)=param(4:13)
            end if
         else if(words(2)=='TAUSM') then
            if(words(3)=='STEPS') then
               a%npp_stepi(9) = a%Listener%getint('STEPS',1,'#Postprocess step interval for Tausmoothing')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,9)=param(4:13)
            end if
         else if(words(2)=='VORTI') then
            if(words(3)=='STEPS') then
               a%npp_stepi(10) = a%Listener%getint('STEPS',1,'#Postprocess step interval for Vorticity')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,10)=param(4:13)
            end if
         else if(words(2)=='DIVER') then
            if(words(3)=='STEPS') then
               a%npp_stepi(11) = a%Listener%getint('STEPS',1,'#Postprocess step interval for Velocity divergence')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,11)=param(4:13)
            end if
         else if(words(2)=='KEFEM') then
            if(words(3)=='STEPS') then
               a%npp_stepi(12) = a%Listener%getint('STEPS',1,'#Postprocess step interval for fem Kinetic energy')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,12)=param(4:13)
            end if
         else if(words(2)=='KESGS') then
            if(words(3)=='STEPS') then
               a%npp_stepi(13) = a%Listener%getint('STEPS',1,'#Postprocess step interval for sgs Kinetic energy')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,13)=param(4:13)
            end if
         else if(words(2)=='VESGS') then
            if(words(3)=='STEPS') then
               a%npp_stepi(14) = a%Listener%getint('STEPS',1,'#Postprocess step interval for vesgs ')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,14)=param(4:13)
            end if
         else if(words(2)=='PRSGS') then
            if(words(3)=='STEPS') then
               a%npp_stepi(26) = a%Listener%getint('STEPS',1,'#Postprocess step interval for prsgs ')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,26)=param(4:13)
            end if
         else if(words(2)=='DISSI') then
            if(words(3)=='STEPS') then
               a%npp_stepi(15) = a%Listener%getint('STEPS',1,'#Postprocess step interval for dissipation')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,15)=param(4:13)
               if (a%npp_stepi(15) /= 0) a%kfl_dispa=1
            end if
         elseif(words(2)=='AVGVE') then
            if (words(3) == 'STEPS') then
               a%npp_stepi(16) = a%Listener%getint('STEPS',1,'#Postprocess step interval for avg1D_velocity')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,16)=param(4:13)
            endif

         elseif(words(2)=='AVGDI') then
            if(words(3)=='STEPS') then
               a%npp_stepi(17) = a%Listener%getint('STEPS',1,'#Postprocess step interval for avg1D_Dissipation')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,17)=param(4:13)
            endif
         !Residual
         elseif(words(2)=='RESID') then
            if(words(3)=='STEPS') then
               a%npp_stepi(18) = a%Listener%getint('STEPS',1,'#Postprocess step interval for avg1D_Dissipation')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,18)=param(4:13)
            endif

         elseif(words(2)=='REPRO') then
            if(words(3)=='STEPS') then
               a%npp_stepi(19) = a%Listener%getint('STEPS',1,'#Postprocess step interval for avg1D_Dissipation')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,19)=param(4:13)
            endif

         elseif(words(2)=='GPDIS') then
            if(words(3)=='STEPS') then
               a%npp_stepi(20) = a%Listener%getint('STEPS',1,'#Postprocess step interval for avg1D_Dissipation')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,20)=param(4:13)
            endif

         else if(words(2)=='QFACT') then
            if(words(3)=='STEPS') then
               a%npp_stepi(21) = a%Listener%getint('STEPS',1,'#Postprocess step interval for Q-factor')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,21)=param(4:13)
            end if

         elseif(words(2)=='BOUND') then
            if(words(3)=='STEPS') then
               a%npp_stepi(22) = a%Listener%getint('STEPS',1,'#Postprocess step interval for boundary traction')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,22)=param(4:13)
            end if
            if(a%npp_stepi(22).gt.0 ) then
                a%kfl_postBtract=1
            end if

          elseif(words(2)=='RELAX') then
            if(words(3)=='STEPS') then
               a%npp_stepi(23) = a%Listener%getint('STEPS',1,'#Postprocess step interval for boundary traction')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,23)=param(4:13)
            end if

          else if(words(2)=='STATI') then
            if(words(3)=='STEPS') then
               a%npp_stepi(24) = a%Listener%getint('STEPS',1,'#Postprocess step interval for Statistics')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,24)=param(4:13)
            end if
            if(words(4)=='START') then
               a%StatisticsStartTime = a%Listener%getint('START',0,'#Postprocess start time for Statistics')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,24)=param(4:13)
            end if

          !External Forces
          else if(words(2)=='EXTER') then
            if(words(3)=='STEPS') then
               a%npp_stepi(25) = a%Listener%getint('STEPS',1,'#Postprocess step interval for Statistics')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,25)=param(4:13)
            end if
            if(words(4)=='START') then
               a%StatisticsStartTime = a%Listener%getint('START',0,'#Postprocess start time for Statistics')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,25)=param(4:13)
            end if

          !Exact solution
          elseif(words(2) == 'EXACT') then
            if(words(3) == 'STEPS') then
               a%npp_stepi(27) = a%Listener%getint('STEPS',1,'#Postprocess step interval for Statistics')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,27)=param(4:13)
            endif
            if(words(4)=='START') then
               a%StatisticsStartTime = a%Listener%getint('START',0,'#Postprocess start time for Statistics')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,27)=param(4:13)
            endif

          elseif(words(2)=='TAUST') then
            if(words(3)=='make STEPS') then
               a%npp_stepi(28) = a%Listener%getint('STEPS',1,'#Postprocess step interval for stabilization taus')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,27)=param(4:13)
            end if

          elseif(words(2)=='PSIVA') then
            if(words(3)=='STEPS') then
               a%npp_stepi(29) = a%Listener%getint('STEPS',1,'#Postprocess step interval for internal variable psi')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,28)=param(4:13)
            end if

            elseif(words(2)=='PSIRE') then
            if(words(3)=='STEPS') then
               a%npp_stepi(30) = a%Listener%getint('STEPS',1,'#Postprocess step interval for real variable psi')
               if(a%Listener%exists('ATTIM')) a%pos_times(1:10,29)=param(4:13)
            end if

         end if

      elseif(words(1)=='AVG1D') then
         if (a%Listener%exists('X    ')) then
            a%Avg1DIdime = 1
         elseif (a%Listener%exists('Y    ')) then
            a%Avg1DIdime = 2
         elseif (a%Listener%exists('Z    ')) then
            a%Avg1DIdime = 3
         endif

      ! Step or time for post-process
      !Force dissipation computation
      elseif(words(1)=='DISCO') then
         if (a%Listener%exists('ON   ')) a%kfl_dispa = 2

      ! Forces and Moments.
      elseif(words(1)=='FORCE') then
         a%kfl_outfm = 1
         a%adimf(1) = a%Listener%getrea('FX   ',0.0_rp,'#x-component of CF')
         a%adimf(2) = a%Listener%getrea('FY   ',0.0_rp,'#y-component of CF')
         a%adimf(3) = a%Listener%getrea('FZ   ',0.0_rp,'#z-component of CF')
         a%adimm(1) = a%Listener%getrea('MX   ',0.0_rp,'#x-component of CM')
         a%adimm(2) = a%Listener%getrea('MY   ',0.0_rp,'#y-component of CM')
         a%adimm(3) = a%Listener%getrea('MZ   ',0.0_rp,'#z-component of CM')
         a%origm(1) = a%Listener%getrea('OX   ',0.0_rp,'#x-component of ori')
         a%origm(2) = a%Listener%getrea('OY   ',0.0_rp,'#y-component of ori')
         a%origm(3) = a%Listener%getrea('OZ   ',0.0_rp,'#z-component of ori')
      end if

   elseif (itask == 100) then
      if (a%nptra/=0) then
         if (a%MPIrank == a%MPIroot) then
            a%fil_trap = adjustl(trim(a%namda))//'.nsi.tp'
            call iofile(zero,a%lun_trap,a%fil_trap,'NSI TRACKING OF POINTS')
         endif
      endif
   endif

   1 format(a)
   5 format(a,1x,i2)

   101 format(5x,'WARNING: ',a)

end subroutine nsi_reaous

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


   subroutine  getcount_bufr(inpfile,nmsg,nsub)

!$$$  subprogram documentation block
!  ............................................................
!   subprogram:  getcount
!   prgmmr:     Woollen, Jack, Su, Xiujuan
!  abstract:   this subroutine is to read the bufr file to get information on the
!              counts of message and subset.

!  program history log
!  2015-03-27  Su    Modify original code from Jack Woollen

!  input argument list
!     inpfile       - input bufr file

!  output argument list:
!     nmsg         - messge count from input bufr file
!     nsub         - subset count from input bufr file
!$$$

use kinds, only: i_kind
use file_utility, only : get_lun

   implicit none
!  Declare passed variables
   character(len=*)                      ,intent(in ) :: inpfile
   integer(i_kind)                       ,intent(out) :: nmsg,nsub

!  Declare local parameters

   character(len=8)  :: subset
   integer(i_kind)   lunit,ireadmg,nmsub,idate

   lunit=get_lun()
   nsub=0;nmsg=0
   call closbf(lunit)
   open(lunit,file=trim(inpfile),form='unformatted')
   call openbf(lunit,'IN',lunit)
   do while(ireadmg(lunit,subset,idate) >=0)
      nmsg = nmsg+1; nsub = nsub+nmsub(lunit)
   enddo
   call closbf(lunit)

   return
   end


text:
subroutine	convol(x,nx,lsfarr,y)

!
! 'Convolving' a flux array x with a 1D or 2D (wavelength-dependent) lsf
!

use share, only: dp,lsf,mlsf,nlsf,nsynth,hs

implicit none

!input/output
real(dp),intent(in) 	:: x(nx)	  			!input data array
integer, intent(in)     :: nx 					!size of input data array
real(dp), intent(in)    :: lsfarr(mlsf,nlsf)    ! lsfarray
real(dp),intent(out) 	:: y(nx)	  			!output data array

!locals
integer			:: i,j,k,ii,n,range
real(dp) 		:: tot
real(dp)     	:: y2(nx)	  					!tmp holder for output data

!write(*,*)'lsfarr in convol=', lsfarr(1,:)
!write(*,*)'x in convol=',x

range=(nlsf-1)/2
!write(*,*)'nlsf,range=',nlsf,range

!initialize y2
y2(1:nx)=0.0_dp

do j=1,nsynth
	if (abs(hs(j)%res-0._dp) > 1e-6_dp) then	!skip photometry
		n=hs(j)%pixend-hs(j)%pixbegin+1
		do i=hs(j)%pixbegin,hs(j)%pixend

			ii=i-hs(j)%pixbegin+1

			!write(*,*)hs(j)%pixbegin,hs(j)%pixend

			k=1 							!1D
			if (lsf == 2 .or. lsf == 4 .or. lsf == 12 .or. lsf == 14) k=i !2D

			!write(*,*)'i,ii=',i,ii
			if (ii <= range) then
				tot=sum(lsfarr(k,2+range-ii:2*range+1))
				y2(i)=sum( lsfarr(k,2+range-ii:2*range+1)/tot * x(hs(j)%pixbegin:i+range) )
				!write(*,*)'ii <= range, lsf:',2+range-ii,2*range+1
				!write(*,*)'             x  :', hs(j)%pixbegin, i+range
			elseif (ii >= n-range) then
				tot=sum(lsfarr(k,1:range+n-ii+1 ))
				y2(i) = sum( lsfarr(k,1:range+n-ii+1 )/tot * x(i-range:hs(j)%pixend) )
				!write(*,*)'ii <= range, lsf:',1,range+n-ii+1
				!write(*,*)'             x  :', i-range,hs(j)%pixend
			else
				!write(*,*)'i,ii,range,i-range,i+range=',i,ii,range,i-range,i+range
				!write(*,*)'x(i-range:i+range)=',x(i-range:i+range)
				y2(i)=sum( lsfarr(k,:) * x(i-range:i+range) )
				!write(*,*)'lsfarr=', lsfarr(k,:)
				!write(*,*)'y2(i)=',y2(i)
			endif
		enddo
	endif
enddo

!now safely copy to output array
y=y2

end subroutine convol

max_stars_repo_name: dnidever/apogee

max_issues_repo_name: dnidever/apogee

max_forks_repo_name: dnidever/apogee


text: 20070420 23:19:04 nbsp hi, and welcome to the wiki! Users/JessicaRockwell

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
\ I2C temperature measuring temperature with LM75CIM
\ AUTHOR      : Willem Ouwerkerk, November 25, 1999
\ LAST CHANGE : Willem Ouwerkerk, Februari 1, 2022
\ Added compensation value for the LM75 named #COMP
\ Used special character B0 that displays the degrees sign

decimal
 -50 constant #COMP              \ Correction value in half degrees please adjust!
: {REGISTER     ( reg +n -- )   {i2c-write  bus! ;   \ Select any LM75 register
: TEMP-REGISTER ( -- )          0 1 {register i2c} ; \ Select temp. register

hex
: @REG8         ( reg -- b )    \ Read 8-bit register
    1 {register i2c}  1 {i2c-read  bus@ i2c} ;

: @REG16        ( reg -- x )    \ Read 16-bit register
    1 {register i2c}  2 {i2c-read  bus@ bus@ i2c} swap b+b ;

\ Temperature is given in half degrees for each bit, where zero
\ is 0 degrees celcius, the range goes from -55 tot 125 degrees.
: TEMPERATURE)   ( -- tr )          \ Read raw temperature
    2 {i2c-read  bus@ bus@ i2c} swap b+b ;

: TEMPERATURE   ( -- n )            \ Read corrected temperature
    temperature)  dup 0F rshift     \ Shift 15 bits to right
    if  -1  FFFF xor  or  then      \ Convert sign to systems word width
    2/ 2/ 2/ 2/ 2/ 2/ 2/  #comp + ; \ Divide by 128, keep sign & compensation value

: CONFIGURATION ( b  -- )           \ Set LM75 configuration
    1 1 {register  i2c}             \ Select config. reg.
    temp-register ;                 \ Select temp. reg. again

\ The temperatute T needs to be given in whole degrees here
: >TEMPERATURE  ( t -- tl th )      \ Convert temp. to two bytes
    0  swap FF and ;                \ for LM75 thermostatic function

: LOW-LIMIT     ( t -- )
    >temperature  2 3 {register     \ Select THYST reg.
    bus! bus! i2c}  temp-register ; \ Set thermostat lower boundary

: HIGH-LIMIT    ( t -- )
    >temperature  3 3 {register     \ Select TOS reg. (over temperature)
    bus! bus! bus! i2c}  temp-register ; \ Set thermostat higher boundary

i2c-on  4C device!  \ Set LM75 device address

\ Example programs

\ Show temperature in whole degrees celcius
: .CELCIUS1         ( n -- )        2/ 3 .r  BA emit ." C " ;

: TEMPERATURE1      ( -- )
    i2c-on  4C device!
    base @ >r decimal
    begin
        temperature .celcius1  100 ms
    key? until  r> base ! ;

\ Show temperature in half degrees celsius
\ temperature dm 10 dm 2 */ dm .
: .CELCIUS2         ( n -- )
    0A * 2 / s>d <# # ch . hold #s #> type  BA emit ." C " ;

: TEMPERATURE2      ( -- )
    i2c-on  4C device!
    base @ >r  decimal
    begin
        temperature .celcius2  100 ms
    key? until  r> base ! ;

\ End ;;;

max_stars_repo_name: embeddingforth/embeddinForth

max_issues_repo_name: embeddingforth/embeddinForth

max_forks_repo_name: embeddingforth/embeddinForth


text:
C     HELLO WORLD IN FORTRAN IV
      PROGRAM HELLO
      WRITE (*,100)
      STOP
  100 FORMAT (' HELLO WORLD! ' /)
      END

max_stars_repo_name: barrettotte/punch

max_issues_repo_name: barrettotte/punch

max_forks_repo_name: barrettotte/punch


text:
!------------------------------------------------------------------------------
!    This code is part of the MondoSCF suite of programs for linear scaling
!    electronic structure theory and ab initio molecular dynamics.
!
!    Copyright (2004). The Regents of the University of California. This
!    material was produced under U.S. Government contract W-7405-ENG-36
!    for Los Alamos National Laboratory, which is operated by the University
!    of California for the U.S. Department of Energy. The U.S. Government has
!    rights to use, reproduce, and distribute this software.  NEITHER THE
!    GOVERNMENT NOR THE UNIVERSITY MAKES ANY WARRANTY, EXPRESS OR IMPLIED,
!    OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.
!
!    This program is free software; you can redistribute it and/or modify
!    it under the terms of the GNU General Public License as published by the
!    Free Software Foundation; either version 2 of the License, or (at your
!    option) any later version. Accordingly, 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 at www.gnu.org for details.
!
!    While you may do as you like with this software, the GNU license requires
!    that you clearly mark derivative software.  In addition, you are encouraged
!    to return derivative works to the MondoSCF group for review, and possible
!    disemination in future releases.
!------------------------------------------------------------------------------
!    PACKAGE FOR THE SPARSE APPROXIMATE MATRIX MULTIPLY (SPAMMPACK)
!    Matt Challacombe and Nick Bock
!------------------------------------------------------------------------------
MODULE SpAMM_TYPES
  IMPLICIT NONE
  !
  INTEGER, PARAMETER :: INT1=SELECTED_INT_KIND(2)  !--Integer*1
  INTEGER, PARAMETER :: INT2=SELECTED_INT_KIND(4)  !--Integer*2
  INTEGER, PARAMETER :: INT4=SELECTED_INT_KIND(9)  !--Integer*4
  INTEGER, PARAMETER :: INT8=SELECTED_INT_KIND(18) !--Integer*8
  !
  INTEGER, PARAMETER :: SINGLE=KIND(0.0)           !--Real*4
  INTEGER, PARAMETER :: DOUBLE=KIND(0.0D0)         !--Real*8
  REAL(DOUBLE),PARAMETER :: Zero=0D0,Half=5D-1,One=1D0,Two=2D0,Eight=8D0
  !
  TYPE QuTree
     INTEGER                      :: Siz
     INTEGER                      :: Lev     ! Level of this box
     INTEGER                      :: Num     ! Box number
     REAL(DOUBLE)                 :: Norm
     INTEGER,      DIMENSION(2,2) :: Box
     TYPE(QuTree), POINTER        :: Quad00,Quad01,Quad10,Quad11
     REAL(DOUBLE), DIMENSION(:,:), ALLOCATABLE :: Blok
  END TYPE QuTree
  !
  TYPE SpAMM_cubes
     TYPE(SpAMM_cubes), POINTER  :: Next
     INTEGER                     :: Lev
     INTEGER                     :: Hash
     REAL(DOUBLE)                :: Bound
     INTEGER,     DIMENSION(3,2) :: Box
  END TYPE SpAMM_cubes
  !
  !----------------------------------------------------------------
  REAL(DOUBLE)       :: SpAMM_tolerance, SpAMM_multiplies, SpAMM_average_multiplies, &
                        SpAMM_dimension
  INTEGER            :: SpAMM_tiles, SpAMM_levels, SpAMM_quadcount
  TYPE(SpAMM_cubes),POINTER :: SpAMM_stream
!  INTEGER, PARAMETER :: SpAMM_BLOCK_SIZE=4
  INTEGER :: SpAMM_BLOCK_SIZE

  REAL(DOUBLE) :: GlobalDropped,GlobalError

  !----------------------------------------------------------------
  !===============================================================================
  !  Interface blocks for generic linear algebra routines
  !===============================================================================
  INTERFACE Trace
     MODULE PROCEDURE Trace_qutree, Trace_dns
  END INTERFACE

  INTERFACE Copy
     MODULE PROCEDURE Copy_qutree_eq_qutree
  END INTERFACE

  INTERFACE Add
     MODULE PROCEDURE Add_qutree_pls_qutree
  END INTERFACE

  INTERFACE Multiply
     MODULE PROCEDURE Multiply_qutree_tms_qutree,Multiply_qutree_tms_scalar
  END INTERFACE
  !
CONTAINS
  !
  FUNCTION Dense2Quad(A) RESULT(qA)
    REAL(DOUBLE),DIMENSION(:,:) :: A
    TYPE(QuTree),POINTER        :: qA
    !
    SpAMM_quadcount=0
    !
    CALL NewQuNode(qA,init=.TRUE.)
    CALL SpAMM_Mat2Quad(A,qA)
    qA%Norm=SQRT(SpAMM_NormReduce(qA))
    CALL Filter_qutree(qA,1D-16)
  END FUNCTION Dense2Quad

  RECURSIVE SUBROUTINE SpAMM_Mat2Quad(A,qA)
    REAL(DOUBLE),DIMENSION(:,:) :: A
    TYPE(QuTree),POINTER        :: qA
    INTEGER                     :: I,J,Delta,Tier
    REAL(DOUBLE)                :: Norm
    !
    I=SIZE(A,1)
    J=SIZE(A,2)
    IF(I<=SpAMM_BLOCK_SIZE.AND.J<=SpAMM_BLOCK_SIZE)THEN
       IF(I<SpAMM_BLOCK_SIZE)THEN
          WRITE(*,*)' I = ',I,' J = ',J,' Blok Siz = ',SpAMM_BLOCK_SIZE
          STOP ' LOGIC ERROR IN SpAMM: padding error '
       ELSE
          qA%Siz=SpAMM_BLOCK_SIZE
          ALLOCATE(qA%Blok(SpAMM_BLOCK_SIZE,SpAMM_BLOCK_SIZE))
          qA%Blok(1:I,1:J)=A(1:I,1:J)
!          qA%Norm=SUM(A(1:I,1:J)**2)
          NULLIFY(qA%Quad00)
          NULLIFY(qA%Quad01)
          NULLIFY(qA%Quad10)
          NULLIFY(qA%Quad11)
       ENDIF
       RETURN
    ELSE
       !
       ALLOCATE(qA%Quad00)
       ALLOCATE(qA%Quad01)
       ALLOCATE(qA%Quad10)
       ALLOCATE(qA%Quad11)

       Tier=qA%Lev+1
       qA%Quad00%Lev=Tier
       qA%Quad01%Lev=Tier
       qA%Quad10%Lev=Tier
       qA%Quad11%Lev=Tier
       !
       qA%Quad00%Num=SpAMM_quadcount+1
       qA%Quad01%Num=SpAMM_quadcount+2
       qA%Quad10%Num=SpAMM_quadcount+3
       qA%Quad11%Num=SpAMM_quadcount+4
       SpAMM_quadcount=SpAMM_quadcount+4
       !
       qA%Siz=(qA%Box(1,2)-qA%Box(1,1))+1
       Delta=(qA%Siz-1)/2
       !
       qA%Quad00%Box(:,1)=qA%Box(:,1)
       qA%Quad00%Box(:,2)=qA%Box(:,1)+(/Delta,Delta/)
       !
       qA%Quad10%Box(:,1)=(/qA%Box(1,1)+Delta+1,qA%Box(2,1)/)
       qA%Quad10%Box(:,2)=(/qA%Box(1,2),qA%Box(2,1)+Delta/)
       !
       qA%Quad01%Box(:,1)=(/qA%Box(1,1),qA%Box(2,1)+Delta+1/)
       qA%Quad01%Box(:,2)=(/qA%Box(1,1)+Delta,qA%Box(2,2)/)
       !
       qA%Quad11%Box(:,1)=(/qA%Box(1,1)+Delta+1,qA%Box(2,1)+Delta+1/)
       qA%Quad11%Box(:,2)=qA%Box(:,2)
       !
       CALL SpAMM_Mat2Quad(A(1:I/2  ,1:J/2  )  , qA%Quad00 )
       CALL SpAMM_Mat2Quad(A(1:I/2  ,J/2+1:J)  , qA%Quad01 )
       CALL SpAMM_Mat2Quad(A(I/2+1:I,1:J/2  )  , qA%Quad10 )
       CALL SpAMM_Mat2Quad(A(I/2+1:I,J/2+1:J)  , qA%Quad11 )
!       !
!       qA%Norm=qA%Quad00%Norm+qA%Quad01%Norm+qA%Quad10%Norm+qA%Quad11%Norm
!       !
!       qA%Quad00%Norm=SQRT(qA%Quad00%Norm)
!       qA%Quad01%Norm=SQRT(qA%Quad01%Norm)
!       qA%Quad10%Norm=SQRT(qA%Quad10%Norm)
!       qA%Quad11%Norm=SQRT(qA%Quad11%Norm)
    ENDIF
    !
  END SUBROUTINE SpAMM_Mat2Quad

  FUNCTION Quad2Dense(qA) RESULT(A)
    REAL(DOUBLE),DIMENSION(SpAMM_tiles*SpAMM_BLOCK_SIZE,SpAMM_tiles*SpAMM_BLOCK_SIZE) :: A
    TYPE(QuTree),POINTER :: qA
    A=Zero
    CALL SpAMM_Quad2Mat(A,qA)
  END FUNCTION Quad2Dense

  RECURSIVE SUBROUTINE SpAMM_Quad2Mat(A,qA)
    INTEGER :: I1,I2,J1,J2
    REAL(DOUBLE),DIMENSION(:,:) :: A
    TYPE(QuTree),POINTER :: qA

    IF(.NOT.ASSOCIATED(qA))RETURN
    IF(qA%Norm==Zero)RETURN
    !
    IF(ALLOCATED(qA%Blok))THEN
       I1=qA%Box(1,1)
       I2=qA%Box(1,2)
       J1=qA%Box(2,1)
       J2=qA%Box(2,2)
       A(I1:I2,J1:J2)=qA%Blok
       RETURN
    ELSE
       CALL SpAMM_Quad2Mat(A,qA%Quad00)
       CALL SpAMM_Quad2Mat(A,qA%Quad01)
       CALL SpAMM_Quad2Mat(A,qA%Quad10)
       CALL SpAMM_Quad2Mat(A,qA%Quad11)
    ENDIF
  END SUBROUTINE SpAMM_Quad2Mat

  SUBROUTINE DeleteQuad(qA)
    TYPE(QuTree),POINTER :: qA
    CALL SpAMM_DeleteQuad(qA)
    DEALLOCATE(qA)
    NULLIFY(qA)
  END SUBROUTINE DeleteQuad

  RECURSIVE SUBROUTINE SpAMM_DeleteQuad(qA)
    TYPE(QuTree),POINTER  :: qA
    INTEGER :: Status
    IF(ALLOCATED(qA%Blok))THEN
       DEALLOCATE(qA%Blok,STAT=Status)
       NULLIFY(qA%Blok)
       SpAMM_quadcount=SpAMM_quadcount-1
    ENDIF
    IF(ASSOCIATED(qA%Quad00))THEN
       CALL SpAMM_DeleteQuad(qA%Quad00)
       DEALLOCATE(qA%Quad00)
       NULLIFY(qA%Quad00)
       SpAMM_quadcount=SpAMM_quadcount-1
    ENDIF
    IF(ASSOCIATED(qA%Quad01))THEN
       CALL SpAMM_DeleteQuad(qA%Quad01)
       DEALLOCATE(qA%Quad01)
       NULLIFY(qA%Quad01)
       SpAMM_quadcount=SpAMM_quadcount-1
    ENDIF
    IF(ASSOCIATED(qA%Quad10))THEN
       CALL SpAMM_DeleteQuad(qA%Quad10)
       DEALLOCATE(qA%Quad10)
       NULLIFY(qA%Quad10)
       SpAMM_quadcount=SpAMM_quadcount-1
    ENDIF
    IF(ASSOCIATED(qA%Quad11))THEN
       CALL SpAMM_DeleteQuad(qA%Quad11)
       DEALLOCATE(qA%Quad11)
       NULLIFY(qA%Quad11)
       SpAMM_quadcount=SpAMM_quadcount-1
    ENDIF
  END SUBROUTINE SpAMM_DeleteQuad

  SUBROUTINE NewQuNode(qA,count,init)
    LOGICAL,OPTIONAL :: init
    INTEGER,OPTIONAL :: count
    TYPE(QuTree), POINTER :: qA
    IF(PRESENT(init))THEN
       ALLOCATE(qA)
       qA%Num=1
       qA%Lev=1
       qA%Siz=SpAMM_tiles*SpAMM_BLOCK_SIZE
       qA%Box(:,1)=(/1,1/)
       qA%Box(:,2)=(/SpAMM_tiles,SpAMM_tiles/)*SpAMM_BLOCK_SIZE
    ELSE
       IF(.NOT.ASSOCIATED(qA))THEN
          ALLOCATE(qA)
       ENDIF
       IF(PRESENT(count))THEN
          count=count+1
          qA%Num=count
       ELSE
          SpAMM_quadcount=SpAMM_quadcount+1
          qA%Num=SpAMM_quadcount
       ENDIF
    ENDIF
    qA%Norm=Zero
    NULLIFY(qA%Quad00)
    NULLIFY(qA%Quad01)
    NULLIFY(qA%Quad10)
    NULLIFY(qA%Quad11)
  END SUBROUTINE NewQuNode

  RECURSIVE SUBROUTINE Print_Quad(qA)
    TYPE(QuTree),POINTER :: qA

    IF(ASSOCIATED(qA))THEN
       WRITE(*,111)qA%Box(:,1),qA%Box(:,2),qA%Num,qA%Siz,qA%Norm
    ENDIF

    IF(ASSOCIATED(qA%Quad00))&
         WRITE(*,111)qA%Quad00%Box(:,1),qA%Quad00%Box(:,2),qA%Quad00%Num,qA%Quad00%Siz,qA%Quad00%Norm
    IF(ASSOCIATED(qA%Quad01))&
         WRITE(*,111)qA%Quad01%Box(:,1),qA%Quad01%Box(:,2),qA%Quad01%Num,qA%Quad01%Siz,qA%Quad01%Norm
    IF(ASSOCIATED(qA%Quad10))&
         WRITE(*,111)qA%Quad10%Box(:,1),qA%Quad10%Box(:,2),qA%Quad10%Num,qA%Quad10%Siz,qA%Quad10%Norm
    IF(ASSOCIATED(qA%Quad11))&
         WRITE(*,111)qA%Quad11%Box(:,1),qA%Quad11%Box(:,2),qA%Quad11%Num,qA%Quad11%Siz,qA%Quad11%Norm
    WRITE(*,*)' ========================================'
    !
    IF(ASSOCIATED(qA%Quad00)) &
         CALL Print_Quad(qA%Quad00)
    IF(ASSOCIATED(qA%Quad01)) &
         CALL Print_Quad(qA%Quad01)
    IF(ASSOCIATED(qA%Quad10)) &
         CALL Print_Quad(qA%Quad10)
    IF(ASSOCIATED(qA%Quad11)) &
         CALL Print_Quad(qA%Quad11)

111 FORMAT("Cuboid[{",I8,",",I8,":",I8,",",I8,"}, (* ",I6,", ",I6,", ",F12.6,"*)")

  END SUBROUTINE Print_Quad
  !=================================================================
  ! QuTree trace routines: a=trace[A]
  !=================================================================
  RECURSIVE FUNCTION Trace_qutree(qA) RESULT(Trc)
    IMPLICIT NONE
    TYPE(QuTree), POINTER  :: qA
    REAL(DOUBLE) :: trc
    INTEGER :: I

    IF(qA%Siz==SpAMM_BLOCK_SIZE)THEN
       trc=Zero
       IF(.NOT.ALLOCATED(qA%Blok))RETURN
       DO I=1,SpAMM_BLOCK_SIZE
          trc=trc+qA%Blok(I,I)
       ENDDO
    ELSEIF(.NOT.ASSOCIATED(qA%Quad00).AND..NOT.ASSOCIATED(qA%Quad11))THEN
       Trc=Zero
   ELSEIF(.NOT.ASSOCIATED(qA%Quad11))THEN
       trc=Trace_qutree(qA%Quad00)
    ELSEIF(.NOT.ASSOCIATED(qA%Quad00))THEN
       trc=Trace_qutree(qA%Quad11)
    ELSE
       trc=Trace_qutree(qA%Quad00)+Trace_qutree(qA%Quad11)
    ENDIF
  END FUNCTION Trace_qutree
  !=================================================================
  ! QuTree filter routines: B=filter[A,tau]
  !=================================================================
  RECURSIVE SUBROUTINE Filter_qutree(qA,tau)
    IMPLICIT NONE
    TYPE(QuTree), POINTER  :: qA
    REAL(DOUBLE)           :: tau
    IF(.NOT.ASSOCIATED(qA))RETURN
    IF(qA%Norm<tau)THEN
       CALL DeleteQuad(qA)
    ELSE
       CALL Filter_qutree(qA%Quad00,tau)
       CALL Filter_qutree(qA%Quad01,tau)
       CALL Filter_qutree(qA%Quad10,tau)
       CALL Filter_qutree(qA%Quad11,tau)
    ENDIF
  END SUBROUTINE Filter_qutree
  !=================================================================
  ! QuTree add routines: C=A+B
  !=================================================================
  SUBROUTINE Add_qutree_pls_qutree(qA,qB,qC)
    IMPLICIT NONE
    TYPE(QuTree), POINTER       :: qA,qB,qC
    REAL(DOUBLE) :: norm
    INTEGER :: add_quadcount
    IF(.NOT.ASSOCIATED(qC))THEN
       CALL NewQuNode(qC,init=.TRUE.)
    ELSE
       CALL DeleteQuad(qC)
       CALL NewQuNode(qC,init=.TRUE.)
    ENDIF
    add_quadcount=1
    CALL SpAMM_Add(qA,qB,qC,count=add_quadcount)
    qC%Norm=SQRT(SpAMM_NormReduce(qC))
    CALL Filter_qutree(qC,1D-16)
  END SUBROUTINE Add_qutree_pls_qutree

  RECURSIVE SUBROUTINE SpAMM_Add(qA,qB,qC,count)
    IMPLICIT NONE
    TYPE(QuTree), POINTER    :: qA,qB,qC
    INTEGER,OPTIONAL :: count
    !
    IF(ASSOCIATED(qA).AND.ASSOCIATED(qB))THEN
       IF(.NOT.ASSOCIATED(qC))THEN
          CALL NewQuNode(qC,count=count)
       ELSE
          count=count+1
          qC%Num=count
       ENDIF
       qC%Lev=qA%Lev
       qC%Siz=qA%Siz
       qC%Box=qA%Box
       !
       IF(qA%Siz==SpAMM_BLOCK_SIZE)THEN
          IF(.NOT.ALLOCATED(qC%Blok))THEN
             ALLOCATE(qC%Blok(1:SpAMM_BLOCK_SIZE,1:SpAMM_BLOCK_SIZE))
             qC%Blok=0D0
          ENDIF
          ! Add
          IF(ALLOCATED(qA%Blok).AND.ALLOCATED(qB%Blok))THEN
             qC%Blok=qA%Blok+qB%Blok
             qC%Lev=qA%Lev
          ELSEIF(ALLOCATED(qA%Blok))THEN
             qC%Blok=qA%Blok
             qC%Lev=qA%Lev
          ELSEIF(ALLOCATED(qB%Blok))THEN
             qC%Blok=qB%Blok
             qC%Lev=qB%Lev
          ELSE
             CALL DeleteQuad(qC)
          ENDIF
       ELSE
          CALL SpAMM_Add(qA%Quad00,qB%Quad00,qC%Quad00,count)
          CALL SpAMM_Add(qA%Quad01,qB%Quad01,qC%Quad01,count)
          CALL SpAMM_Add(qA%Quad10,qB%Quad10,qC%Quad10,count)
          CALL SpAMM_Add(qA%Quad11,qB%Quad11,qC%Quad11,count)
       ENDIF
    ELSEIF(.NOT.ASSOCIATED(qA))THEN
       CALL SpAMM_Copy(qB,qC,count=count)
    ELSEIF(.NOT.ASSOCIATED(qB))THEN
       CALL SpAMM_Copy(qA,qC,count=count)
    ELSE
       CALL DeleteQuad(qC)
    ENDIF
  END SUBROUTINE SpAMM_Add
  !=================================================================
  ! QuTree copy routines: C=A
  !=================================================================
  FUNCTION Copy_qutree_eq_qutree(qA) RESULT(qC)
    IMPLICIT NONE
    TYPE(QuTree), POINTER       :: qA,qC
    INTEGER :: node_count
    CALL NewQuNode(qC,init=.TRUE.)
    qC%Num=1
    qC%Lev=1
    node_count=1
    qC%Box(:,1)=(/1,1/)
    qC%Box(:,2)=(/SpAMM_tiles,SpAMM_tiles/)*SpAMM_BLOCK_SIZE
    CALL SpAMM_Copy(qA,qC,count=node_count)
  END FUNCTION Copy_qutree_eq_qutree

  RECURSIVE SUBROUTINE SpAMM_Copy(qA,qC,count)
    IMPLICIT NONE
    TYPE(QuTree), POINTER    :: qA,qC
    INTEGER, OPTIONAL :: count
    !
    IF(.NOT.ASSOCIATED(qA))RETURN
    !
    IF(.NOT.ASSOCIATED(qC))THEN
       CALL NewQuNode(qC,count=count)
    ELSEIF(PRESENT(count))THEN
       count=count+1
       qC%Num=count
    ELSE
       SpAMM_quadcount=SpAMM_quadcount+1
       qC%Num=SpAMM_quadcount
    ENDIF
    qC%Siz=qA%Siz
    qC%Lev=qA%Lev
    qC%Box=qA%Box
    qC%Norm=qA%Norm
    IF(qA%Siz==SpAMM_BLOCK_SIZE.AND.ALLOCATED(qA%Blok))THEN
       IF(.NOT.ALLOCATED(qC%Blok)) &
          ALLOCATE(qC%Blok(1:SpAMM_BLOCK_SIZE,1:SpAMM_BLOCK_SIZE))
       qC%Blok=qA%Blok
    ELSE
       IF(ASSOCIATED(qA%Quad00))THEN
          CALL SpAMM_Copy(qA%Quad00,qC%Quad00,count=count)
       ELSEIF(ASSOCIATED(qC%Quad00))THEN
          CALL DeleteQuad(qC%Quad00)
       ENDIF
       IF(ASSOCIATED(qA%Quad01))THEN
          CALL SpAMM_Copy(qA%Quad01,qC%Quad01,count=count)
       ELSEIF(ASSOCIATED(qC%Quad01))THEN
          CALL DeleteQuad(qC%Quad01)
       ENDIF
       IF(ASSOCIATED(qA%Quad10))THEN
          CALL SpAMM_Copy(qA%Quad10,qC%Quad10,count=count)
       ELSEIF(ASSOCIATED(qC%Quad10))THEN
          CALL DeleteQuad(qC%Quad10)
       ENDIF
       IF(ASSOCIATED(qA%Quad11))THEN
          CALL SpAMM_Copy(qA%Quad11,qC%Quad11,count=count)
       ELSEIF(ASSOCIATED(qC%Quad11))THEN
          CALL DeleteQuad(qC%Quad11)
       ENDIF
    ENDIF
  END SUBROUTINE SpAMM_Copy
  !=================================================================
  ! QuTree multiply routines
  !=================================================================
  ! A=A*B
  !-----------------------------------------------------------------
  RECURSIVE SUBROUTINE Multiply_qutree_tms_scalar(qA,B)
    IMPLICIT NONE
    TYPE(QuTree), POINTER    :: qA
    REAL(DOUBLE)             :: B
    IF(.NOT.ASSOCIATED(qA))RETURN
    IF(qA%Siz==SpAMM_BLOCK_SIZE.AND.ALLOCATED(qA%Blok))THEN
       qA%Norm=qA%Norm*B
       qA%Blok=qA%Blok*B
    ELSE
       CALL Multiply_qutree_tms_scalar(qA%Quad00,B)
       CALL Multiply_qutree_tms_scalar(qA%Quad01,B)
       CALL Multiply_qutree_tms_scalar(qA%Quad10,B)
       CALL Multiply_qutree_tms_scalar(qA%Quad11,B)
       qA%Norm=qA%Norm*B
    ENDIF
  END SUBROUTINE Multiply_qutree_tms_scalar
  !-----------------------------------------------------------------
  ! C=A*B
  !-----------------------------------------------------------------
  SUBROUTINE Multiply_qutree_tms_qutree(qA,qB,qC,tolerance)
    IMPLICIT NONE
    TYPE(QuTree), POINTER       :: qA,qB,qC
    REAL(DOUBLE),OPTIONAL       :: tolerance
    REAL(DOUBLE)                :: Saved_tolerance,norm,AbsError,RelError,TargetError,ABNorm
    INTEGER                     :: SpAMM_count
    REAL(DOUBLE),ALLOCATABLE,DIMENSION(:,:) :: dA,dB,dC

    IF(PRESENT(tolerance))THEN
       Saved_tolerance=SpAMM_tolerance
       SpAMM_tolerance=tolerance
    ENDIF
    CALL DeleteQuad(qC)
    CALL NewQuNode(qC,init=.TRUE.)
    SpAMM_count=1
    SpAMM_multiplies=0

    GlobalDropped=Zero
    GlobalError=Zero

    CALL SpAMM_Multiply2(qC,qA,qB,count=SpAMM_count)
    qC%Norm=SQRT(SpAMM_NormReduce(qC))
    CALL Filter_qutree(qc,1D-16)

!!$    ABNorm=qA%Norm*qB%Norm
!!$    TargetError=SpAMM_tolerance/SQRT(Eight)
!!$    CALL SpAMM_Multiply(qC,qA,qB,TargetError,count=SpAMM_count)
!!$    qC%Norm=SQRT(SpAMM_NormReduce(qC))
!!$    CALL Filter_qutree(qc,1D-16)
!!$    ALLOCATE(dA(1:SpAMM_tiles*SpAMM_BLOCK_SIZE,1:SpAMM_tiles*SpAMM_BLOCK_SIZE))
!!$    ALLOCATE(dB(1:SpAMM_tiles*SpAMM_BLOCK_SIZE,1:SpAMM_tiles*SpAMM_BLOCK_SIZE))
!!$    ALLOCATE(dC(1:SpAMM_tiles*SpAMM_BLOCK_SIZE,1:SpAMM_tiles*SpAMM_BLOCK_SIZE))
!!$
!!$    dA=Quad2Dense(qA)
!!$    dB=Quad2Dense(qB)
!!$    dC=MATMUL(dA,dB)
!!$    AbsError=NormReduce_dns(dC-Quad2Dense(qC))
!!$    RelError=AbsError/ABNorm !NormReduce_dns(dC)
!!$    DEALLOCATE(dA)
!!$    DEALLOCATE(dB)
!!$    DEALLOCATE(dC)
!!$!    WRITE(*,*)' True  Error          = ',TrueError
!!$    WRITE(*,33)SpAMM_tolerance,GlobalDropped,SpAMM_multiplies,AbsError,RelError
!!$33  FORMAT('Tol = ',D12.6,', Drops = ',D12.6,', Mults = ',D12.6,', Abs Error = ',D12.6,', Rel Error = ',D12.6)

    !
    IF(PRESENT(tolerance))THEN
       SpAMM_tolerance=Saved_tolerance
    ENDIF
  END SUBROUTINE Multiply_qutree_tms_qutree
  !
  RECURSIVE FUNCTION SpAMM_NormReduce(qC) RESULT(norm)
    IMPLICIT NONE
    TYPE(QuTree), POINTER :: qC
    REAL(DOUBLE) :: norm
    ! Associated
    IF(.NOT.ASSOCIATED(qC))THEN
       norm=Zero
       RETURN
    ELSEIF(qC%Siz==SpAMM_BLOCK_SIZE.AND.ALLOCATED(qC%Blok))THEN
       norm=SUM(qC%Blok(1:SpAMM_BLOCK_SIZE,1:SpAMM_BLOCK_SIZE)**2)
       qC%Norm=SQRT(norm)
    ELSE
       norm=SpAMM_NormReduce(qC%Quad00) + &
            SpAMM_NormReduce(qC%Quad01) + &
            SpAMM_NormReduce(qC%Quad10) + &
            SpAMM_NormReduce(qC%Quad11)
       qC%Norm=SQRT(norm)
    ENDIF
  END FUNCTION SpAMM_NormReduce

  FUNCTION NormReduce_dns(A) RESULT(norm)
    REAL(DOUBLE), DIMENSION(:,:) :: A
    REAL(DOUBLE) :: norm
    INTEGER :: I,J
    norm=zero
    DO I=1,SIZE(A,1)
       DO J=1,SIZE(A,2)
          norm=norm+A(I,J)**2
       ENDDO
    ENDDO
    norm=SQRT(norm)
  END FUNCTION NormReduce_dns

  RECURSIVE SUBROUTINE SpAMM_Multiply(qC,qA,qB,TotalError,count)
    IMPLICIT NONE
    TYPE(QuTree), POINTER :: qC,qA,qB
    INTEGER     :: count,LocalCount,LocalDropped
    REAL(DOUBLE) :: TotalError,LocalError,ErrorEstimate
    REAL(DOUBLE):: norm00x00,norm00x01,norm01x11, &
                   norm10x00,norm10x01,norm11x11,norm01x10,norm11x10
    REAL(DOUBLE):: frac00x00,frac01x10,frac00x01,frac01x11, &
                   frac10x00,frac11x10,frac10x01,frac11x11
    LOGICAL     :: do00x00,do01x10,do10x00,do11x10, &
                   do00x01,do01x11,do10x01,do11x11
    ! Bounds
    qC%Lev=qA%Lev
    qC%Siz=qA%Siz
    qC%Box(1,:)=qA%Box(1,:)
    qC%Box(2,:)=qB%Box(2,:)
    ! Blocks
    IF(qA%Siz==SpAMM_BLOCK_SIZE)THEN
       ! Allocate
       IF(.NOT.ALLOCATED(qC%Blok))THEN
          ALLOCATE(qC%Blok(1:SpAMM_BLOCK_SIZE,1:SpAMM_BLOCK_SIZE))
          qC%Blok=Zero
       END IF
       ! Count
       SpAMM_multiplies=SpAMM_multiplies+1
       ! Accumulate
       qC%Blok=qC%Blok+MATMUL(qA%Blok,qB%Blok)
    ELSE
       ! Nodes exist?
       do00x00=ASSOCIATED(qA%Quad00).AND.ASSOCIATED(qB%Quad00)
       do01x10=ASSOCIATED(qA%Quad01).AND.ASSOCIATED(qB%Quad10)
       do00x01=ASSOCIATED(qA%Quad00).AND.ASSOCIATED(qB%Quad01)
       do01x11=ASSOCIATED(qA%Quad01).AND.ASSOCIATED(qB%Quad11)
       do10x00=ASSOCIATED(qA%Quad10).AND.ASSOCIATED(qB%Quad00)
       do11x10=ASSOCIATED(qA%Quad11).AND.ASSOCIATED(qB%Quad10)
       do10x01=ASSOCIATED(qA%Quad10).AND.ASSOCIATED(qB%Quad01)
       do11x11=ASSOCIATED(qA%Quad11).AND.ASSOCIATED(qB%Quad11)
       !
       norm00x00=Zero
       norm01x10=Zero
       norm00x01=Zero
       norm01x11=Zero
       norm10x00=Zero
       norm11x10=Zero
       norm10x01=Zero
       norm11x11=Zero
       !
       IF(do00x00)norm00x00=qA%Quad00%Norm*qB%Quad00%Norm
       IF(do01x10)norm01x10=qA%Quad01%Norm*qB%Quad10%Norm
       IF(do00x01)norm00x01=qA%Quad00%Norm*qB%Quad01%Norm
       IF(do01x11)norm01x11=qA%Quad01%Norm*qB%Quad11%Norm
       IF(do10x00)norm10x00=qA%Quad10%Norm*qB%Quad00%Norm
       IF(do11x10)norm11x10=qA%Quad11%Norm*qB%Quad10%Norm
       IF(do10x01)norm10x01=qA%Quad10%Norm*qB%Quad01%Norm
       IF(do11x11)norm11x11=qA%Quad11%Norm*qB%Quad11%Norm
       ! SpAMM criteria
       IF(do00x00)do00x00=do00x00.AND.norm00x00>TotalError
       IF(do01x10)do01x10=do01x10.AND.norm01x10>TotalError
       IF(do00x01)do00x01=do00x01.AND.norm00x01>TotalError
       IF(do01x11)do01x11=do01x11.AND.norm01x11>TotalError
       IF(do10x00)do10x00=do10x00.AND.norm10x00>TotalError
       IF(do11x10)do11x10=do11x10.AND.norm11x10>TotalError
       IF(do10x01)do10x01=do10x01.AND.norm10x01>TotalError
       IF(do11x11)do11x11=do11x11.AND.norm11x11>TotalError
       ! Counting the next generation
       LocalCount=0
       IF(do00x00)LocalCount=LocalCount+1
       IF(do01x10)LocalCount=LocalCount+1
       IF(do00x01)LocalCount=LocalCount+1
       IF(do01x11)LocalCount=LocalCount+1
       IF(do10x00)LocalCount=LocalCount+1
       IF(do11x10)LocalCount=LocalCount+1
       IF(do10x01)LocalCount=LocalCount+1
       IF(do11x11)LocalCount=LocalCount+1
       ! Check for done by division
       IF(LocalCount==0)RETURN
       ! global debug/development counters
       GlobalDropped=GlobalDropped+(8-LocalCount)
       GlobalError=GlobalError+LocalError
       ! Statistical estimate of error partitioning
       TotalError=TotalError/SQRT(DBLE(LocalCount))
       ! 00=00*00+01*10
       IF(do00x00.OR.do01x10)THEN
          IF(.NOT.ASSOCIATED(qC%Quad00))CALL NewQuNode(qC%Quad00,count=count)
          IF(do00x00)THEN
             frac00x00=TotalError
             CALL SpAMM_Multiply(qC%Quad00,qA%Quad00,qB%Quad00,frac00x00,count)
          ENDIF
          IF(do01x10)THEN
             frac01x10=TotalError
             CALL SpAMM_Multiply(qC%Quad00,qA%Quad01,qB%Quad10,frac01x10,count)
          ENDIF
       ENDIF
       ! 01=00*01+01*11
       IF(do00x01.OR.do01x11)THEN
          IF(.NOT.ASSOCIATED(qC%Quad01))CALL NewQuNode(qC%Quad01,count=count)
          IF(do00x01)THEN
             frac00x01=TotalError
             CALL SpAMM_Multiply(qC%Quad01,qA%Quad00,qB%Quad01,frac00x01,count)
          ENDIF
          IF(do01x11)THEN
             frac01x11=TotalError
             CALL SpAMM_Multiply(qC%Quad01,qA%Quad01,qB%Quad11,frac01x11,count)
          ENDIF
       ENDIF
       ! 10=10*00+11*10
       IF(do10x00.OR.do11x10)THEN
          IF(.NOT.ASSOCIATED(qC%Quad10))CALL NewQuNode(qC%Quad10,count=count)
          IF(do10x00)THEN
             frac10x00=TotalError
             CALL SpAMM_Multiply(qC%Quad10,qA%Quad10,qB%Quad00,frac10x00,count)
          ENDIF
          IF(do11x10)THEN
             frac11x10=TotalError
             CALL SpAMM_Multiply(qC%Quad10,qA%Quad11,qB%Quad10,frac11x10,count)
          ENDIF
       ENDIF
       ! 11=10*01+11*11
       IF(do10x01.OR.do11x11)THEN
          IF(.NOT.ASSOCIATED(qC%Quad11))CALL NewQuNode(qC%Quad11,count=count)
          IF(do10x01)THEN
             frac10x01=TotalError
             CALL SpAMM_Multiply(qC%Quad11,qA%Quad10,qB%Quad01,frac10x01,count)
          ENDIF
          IF(do11x11)THEN
             frac11x11=TotalError
             CALL SpAMM_Multiply(qC%Quad11,qA%Quad11,qB%Quad11,frac11x11,count)
          ENDIF
       ENDIF
    ENDIF
  END SUBROUTINE SpAMM_Multiply

  RECURSIVE SUBROUTINE SpAMM_Multiply2(qC,qA,qB,count)
    IMPLICIT NONE
    TYPE(QuTree), POINTER :: qC,qA,qB
    INTEGER :: count
    REAL(DOUBLE) :: norm
    REAL(DOUBLE), DIMENSION(1:SpAMM_BLOCK_SIZE,1:SpAMM_BLOCK_SIZE) :: CTmp
    ! Associated
    IF(ASSOCIATED(qA).AND.ASSOCIATED(qB))THEN
    ! Bounds
    IF(qA%Norm*qB%Norm<SpAMM_tolerance)RETURN
    ! Allocate
    IF(.NOT.ASSOCIATED(qC))CALL NewQuNode(qC,count=count)
    ! Boxing
     qC%Siz=qA%Siz
     qC%Lev=qA%Lev
     qC%Box(1,:)=qA%Box(1,:)
     qC%Box(2,:)=qB%Box(2,:)
       ! Blocks
       IF(qA%Siz==SpAMM_BLOCK_SIZE)THEN
          ! Allocate
          IF(.NOT.ALLOCATED(qC%Blok))THEN
             ALLOCATE(qC%Blok(1:SpAMM_BLOCK_SIZE,1:SpAMM_BLOCK_SIZE))
             qC%Blok=Zero
          END IF
          ! Count
          SpAMM_multiplies=SpAMM_multiplies+1
          ! Accumulate
          qC%Blok=qC%Blok+MATMUL(qA%Blok,qB%Blok)
       ELSE
          ! 00=00*00+01*10
          CALL SpAMM_Multiply2(qC%Quad00,qA%Quad00,qB%Quad00,count=count)
          CALL SpAMM_Multiply2(qC%Quad00,qA%Quad01,qB%Quad10,count=count)
          ! 01=00*01+01*11
          CALL SpAMM_Multiply2(qC%Quad01,qA%Quad00,qB%Quad01,count=count)
          CALL SpAMM_Multiply2(qC%Quad01,qA%Quad01,qB%Quad11,count=count)
          ! 10=10*00+11*10
          CALL SpAMM_Multiply2(qC%Quad10,qA%Quad10,qB%Quad00,count=count)
          CALL SpAMM_Multiply2(qC%Quad10,qA%Quad11,qB%Quad10,count=count)
          ! 11=10*01+11*11
          CALL SpAMM_Multiply2(qC%Quad11,qA%Quad10,qB%Quad01,count=count)
          CALL SpAMM_Multiply2(qC%Quad11,qA%Quad11,qB%Quad11,count=count)
       ENDIF
    ENDIF
  END SUBROUTINE SpAMM_Multiply2

  FUNCTION Trace_dns(A) RESULT(trc)
    REAL(DOUBLE) :: trc
    REAL(DOUBLE),DIMENSION(:,:) :: A
    INTEGER I
    trc=Zero
    DO I=1,SIZE(A,1)
       trc=trc+A(I,I)
    ENDDO
  END FUNCTION Trace_dns

  SUBROUTINE TC2_test(P,P2,Tmp1,Norm,TrP,I)
    TYPE(QuTree),POINTER :: P,P2,Tmp1 , Tmp2
    REAL(DOUBLE) :: Norm, CR1, CR2, TrP, TrP2,    TrdP,TrdP2, sum
    INTEGER      :: I,  II,J,KK

    INTEGER, PARAMETER:: NN=1024
    REAL(DOUBLE),DIMENSION(NN,NN) :: dP,dP2,dTmp1, dTmp2,dP2b
    !-------------------------------------------------------------------------------

    dP=Quad2Dense(P)

    IF (I.EQ.1) THEN
       TrP=Trace(P)
       TrdP=Trace(dP)
       ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       !       TrP=TrdP
       IF(ABS(TrP-TrdP)>SpAMM_tolerance*1D1)THEN
          WRITE(*,*)' 1TrP = ',TrP,TrdP
          STOP
       ENDIF
       ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    ENDIF
    !
    CALL Multiply(P,P,P2)             ! The only multiplication is a square
    dP2=MATMUL(dP,dP)

!    P2=>Dense2Quad(Quad2Dense(P2))
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    dP2b=Quad2Dense(P2)
    DO II=1,NN
       DO J=1,NN
          IF(ABS(dP2(II,J)-dP2b(II,J))>SpAMM_tolerance*1D1)THEN
             WRITE(*,*)'++++++++++++++++++++++++++++++++++++++++++++'
             WRITE(*,*)' sparse P*P'

             sum=Zero
             DO KK=1,NN
                IF(ABS(dP(II,KK)*dP(KK,J))>1D-1*SpAMM_tolerance)THEN
                   WRITE(*,*)dP(II,KK)*dP(KK,J),sum
                   sum=sum+dP(II,KK)*dP(KK,J)
                ENDIF
             ENDDO
             WRITE(*,*)' SpAMM_tolerance = ',SpAMM_tolerance
            WRITE(*,*)'ARAW',II,J,dP2(II,J),dP2b(II,J),dP2(II,J)-dP2b(II,J)
            WRITE(*,*)'ASUM',II,J,sum,dP2b(II,J),dP2(II,J)-sum

!             CALL Print_quad(P2)
 !            P=>Dense2Quad(dP2)
             WRITE(*,*)'++++++++++++++++++++++++++++++++++++++++++++'
             WRITE(*,*)' dense P*P'

 !            CALL Print_quad(P)
             STOP
          ENDIF
       ENDDO
    ENDDO
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    TrP2=Trace(P2)
    TrdP2=Trace(dP2)
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !    TrP2=TrdP2
    IF(ABS(TrP2-TrdP2)>SpAMM_tolerance)THEN
       WRITE(*,*)'TrP = ',TrP2,TrdP2
       STOP "P2 STOP "
    ENDIF
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    CR1 = ABS(TrP2-Norm)              ! CR1 = Occupation error criteria
    CR2 = ABS(2.D0*TrP - TrP2 - Norm) ! CR2 = Occupation error criteria

    WRITE(*,33)I,NORM,TrP,SpAMM_multiplies/SpAMM_dimension**3
33  FORMAT(I4,", N =",F12.5,", Tr(P)=",F12.5,", O(N)/N^3 = ",F8.5)

    IF (CR1 < CR2) THEN               ! Too many states
      P=>Copy(P2)
      dP=dP2
      ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      dP2b=Quad2Dense(P)
       DO II=1,NN
          DO J=1,NN
             IF(ABS(dP(II,J)-dP2b(II,J))>SpAMM_tolerance*1D1)THEN
                WRITE(*,*)' B ',II,J,dP(II,J)-dP2b(II,J)
                STOP
             ENDIF
          ENDDO
       ENDDO
      ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    ELSE
       CALL Multiply(P,Two)
       CALL Multiply(P2,-One)
       CALL Add(P,P2,Tmp1)             ! P = 2P-P^2
       P=>Copy(Tmp1)
       dP=Two*dP-dP2
       ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       dP2b=Quad2Dense(Tmp1)
       DO II=1,NN
          DO J=1,NN
             IF(ABS(dP(II,J)-dP2b(II,J))>SpAMM_tolerance*1D1)THEN
                WRITE(*,*)II,J,dP(II,J)-dP2b(II,J)
                STOP
             ENDIF
          ENDDO
       ENDDO
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    ENDIF
    TrP=Trace(P)
    TrdP=Trace(dP)

!    TrP=TrdP

    IF(ABS(TrP-TrdP)>SpAMM_tolerance*1D1)THEN
       WRITE(*,*)'TrP = ',TrP,TrdP
       STOP
    ENDIF

  END SUBROUTINE TC2_TEST

  SUBROUTINE TC2(P,P2,Tmp1,Norm,TrP,I)
    TYPE(QuTree),POINTER :: P,P2,Tmp1 , Tmp2
    REAL(DOUBLE) :: Norm, CR1, CR2, TrP, TrP2,    TrdP,TrdP2
    INTEGER      :: I,  II,J
    !-------------------------------------------------------------------------------
    IF (I.EQ.1)TrP=Trace(P)
    CALL Multiply(P,P,P2)             ! The only multiplication is a square
    TrP2=Trace(P2)
    CR1 = ABS(TrP2-Norm)              ! CR1 = Occupation error criteria
    CR2 = ABS(2.D0*TrP - TrP2 - Norm) ! CR2 = Occupation error criteria
!    WRITE(*,33)I,NORM,TrP,SpAMM_multiplies/SpAMM_tiles**3
33  FORMAT(I4,", N =",F18.8,", Tr(P)=",F18.8,", O(N)/N^3 = ",F8.5)
    IF (CR1 < CR2) THEN               ! Too many states
      P=>Copy(P2)
    ELSE
       CALL Multiply(P,Two)
       CALL Multiply(P2,-One)
       CALL Add(P,P2,Tmp1)             ! P = 2P-P^2
       P=>Copy(Tmp1)
    ENDIF
    TrP=Trace(P)
!    CALL Filter_qutree(P,SpAMM_tolerance*1D-2)
    !
  END SUBROUTINE TC2

  SUBROUTINE TC2_filter(P,P2,Tmp1,Norm,TrP,I)
    TYPE(QuTree),POINTER :: P,P2,Tmp1 , Tmp2
    REAL(DOUBLE) :: Norm, CR1, CR2, TrP, TrP2,    TrdP,TrdP2
    INTEGER      :: I,  II,J
    !-------------------------------------------------------------------------------
    IF(I.EQ.1)THEN
       CALL Filter_qutree(P,SpAMM_tolerance)
       TrP=Trace(P)
    ENDIF
    CALL Multiply_qutree_tms_qutree(P,P,P2,tolerance=1D-100)
    TrP2=Trace(P2)
    CR1 = ABS(TrP2-Norm)              ! CR1 = Occupation error criteria
    CR2 = ABS(2.D0*TrP - TrP2 - Norm) ! CR2 = Occupation error criteria
!    WRITE(*,33)I,NORM,TrP,SpAMM_multiplies/SpAMM_tiles**3
33  FORMAT(I4,", N =",F12.8,", Tr(P)=",F12.8,", O(N)/N^3 = ",F8.5)
    IF (CR1 < CR2) THEN               ! Too many states
       P=>Copy(P2)
    ELSE
       CALL Multiply(P,Two)
       CALL Multiply(P2,-One)
       CALL Add(P,P2,Tmp1)             ! P = 2P-P^2
       P=>Copy(Tmp1)
    ENDIF
    CALL Filter_qutree(P,SpAMM_tolerance)
    TrP=Trace(P)
  END SUBROUTINE TC2_FILTER

END MODULE SpAMM_TYPES

PROGRAM SpAMM_TEST
  USE SpAMM_TYPES
  IMPLICIT NONE
  INTEGER           ::  L,M,N_OLD,N,Nel,I,J,K,II,JJ,TC2_cycles
  REAL(DOUBLE),ALLOCATABLE, DIMENSION(:,:) ::  A, B, C,C2,CDiff,A_NOPADDING
  TYPE(QuTree),POINTER  :: qC,qP,qF,qTmp1,qTmp2
  TYPE(SpAMM_cubes),POINTER  :: Stream
  REAL(DOUBLE)          :: CNorm,CDiffNorm,CErrBound,CElementErr,Opacity,Occ0,Occ1,Occ2,Occ3,TrE, &
                           RelativeErrorE,RelativeErrorN,TargetTrE
  CHARACTER *100        :: Name,Buffer
  LOGICAL :: DoFilter

  CALL GETARG(1,NAME)
  NAME=TRIM(ADJUSTL(NAME))
  CALL GETARG(2,BUFFER)
  READ(BUFFER,*)N
  CALL GETARG(3,BUFFER)
  READ(BUFFER,*)Nel
  CALL GETARG(4,BUFFER)
  READ(BUFFER,*)TC2_cycles
  CALL GETARG(5,BUFFER)
  READ(BUFFER,*)SpAMM_tolerance
  CALL GETARG(6,BUFFER)
  READ(BUFFER,*)DoFilter
  CALL GETARG(7,BUFFER)
  READ(BUFFER,*)TargetTrE
  !--------------------------------------------------
  K=CEILING(LOG10(DBLE(N))/LOG10(2D0))
  N_OLD=N
  N=2**K
  SpAMM_dimension=N

  SpAMM_BLOCK_SIZE=4
!!!  SpAMM_BLOCK_SIZE=N

  SpAMM_tiles=CEILING(DBLE(N)/SpAMM_BLOCK_SIZE)
  SpAMM_levels=CEILING(LOG(DBLE(SpAMM_tiles))/LOG(2D0))+1
  ! ~ ~
  ALLOCATE(A_NOPADDING(1:N_OLD,1:N_OLD))
  !--------------------------------------------------
  OPEN(UNIT=66,FILE=TRIM(NAME)//".f_ortho")
  DO I=1,N_OLD
     DO J=1,N_OLD
        READ(66,*)L,M,A_NOPADDING(L,M)
     ENDDO
  ENDDO
  CLOSE(66)
  ALLOCATE(A(1:N,1:N))
  A=Zero
  A(1:N_OLD,1:N_OLD)=A_NOPADDING(1:N_OLD,1:N_OLD)
  qF=>Dense2Quad(A)
  !--------------------------------------------------
  OPEN(UNIT=66,FILE=TRIM(NAME)//".p_guess_ortho")
  DO I=1,N_OLD
     DO J=1,N_OLD
        READ(66,*)L,M,A_NOPADDING(L,M)
     ENDDO
  ENDDO
  CLOSE(66)
  A=Zero
  A(1:N_OLD,1:N_OLD)=A_NOPADDING(1:N_OLD,1:N_OLD)
  qP=>Dense2Quad(A)
  !  CALL Print_quad(qP)
  !--------------------------------------------------
  DEALLOCATE(A)
  DEALLOCATE(A_NOPADDING)
  !
  NULLIFY(qTmp1)
  NULLIFY(qTmp2)
  CALL NewQuNode(qTmp1,init=.TRUE.)
  CALL NewQuNode(qTmp2,init=.TRUE.)
  !--------------------------------------------------
  SpAMM_average_multiplies=Zero
  IF(.NOT.DoFilter)THEN
!     SpAMM_tolerance=SpAMM_tolerance*DBLE(Nel)/Two
  ENDIF
  !--------------------------------------------------
  Occ0 = 0.D0;Occ1 = 0.D0;Occ2 = 0.D0;Occ3 = 0.D0
  DO I=1,TC2_cycles
     IF(DoFilter)THEN
        CALL TC2_filter(qP,qTmp1,qTmp2,Half*FLOAT(NEl),Occ0,I)
     ELSE
        CALL TC2(qP,qTmp1,qTmp2,Half*FLOAT(NEl),Occ0,I)
     ENDIF
     Occ3 = Occ2;Occ2 = Occ1;Occ1 = Occ0
     SpAMM_average_multiplies=SpAMM_average_multiplies+SpAMM_multiplies
  ENDDO

!  SpAMM_average_multiplies=SpAMM_average_multiplies/DBLE(TC2_cycles)
!  WRITE(*,*)' AV MULTIPLIES = ',SpAMM_average_multiplies

  !

  !
  IF(DoFilter)THEN
     CALL Filter_qutree(qF,SpAMM_tolerance)
     CALL Multiply_qutree_tms_qutree(qP,qF,qTmp1,tolerance=1D-100)
  ELSE
!     CALL Multiply_qutree_tms_qutree(qP,qF,qTmp1,tolerance=1D-100)
     CALL Multiply(qP,qF,qTmp1)
  ENDIF
  TrE=Trace(qTmp1)
  RelativeErrorE=ABS(TrE-TargetTrE)/ABS(TargetTrE)
  RelativeErrorN=ABS(Occ0-Half*DBLE(Nel))/(Half*DBLE(Nel))

  SpAMM_average_multiplies=SpAMM_average_multiplies+SpAMM_multiplies
  SpAMM_average_multiplies=SpAMM_average_multiplies/DBLE(TC2_cycles+1)

!!!  WRITE(*,*)"TrE = ",TrE

  !
  WRITE(*,66)N_old,N,SpAMM_BLOCK_SIZE,SpAMM_tolerance,RelativeErrorE,RelativeErrorN, &
             SpAMM_average_multiplies
66 FORMAT(I6,", ",I6,", ",I3,", ",E8.2,7(", ",E18.8))
  !
END PROGRAM SpAMM_TEST

max_stars_repo_name: FreeON/spammpack

max_issues_repo_name: FreeON/spammpack

max_forks_repo_name: FreeON/spammpack


text:
Address(La Paz Drive) is a residential street in South Davis.

 Intersecting Streets

  El Macero Drive
  Tallow Place
  Hackberry Place


text:
      program prog
      character buf*6
      num = -2
      buf = 'AbCdEf'
      write(buf(2:4), '(i3)') num
      write(6, '(''num = ('',a,'')'')') buf
      end


text:
        SUBROUTINE STEPQF(N,NFE,IFLAG,START,CRASH,HOLD,H,
     $            WK,RELERR,ABSERR,S,Y,YP,YOLD,YPOLD,A,QT,R,
     $            F0,F1,Z0,DZ,W,T,SSPAR,PAR,IPAR)
C
C SUBROUTINE  STEPQF  TAKES ONE STEP ALONG THE ZERO CURVE OF THE
C HOMOTOPY MAP  RHO(LAMBDA,X)  USING A PREDICTOR-CORRECTOR ALGORITHM.
C THE PREDICTOR USES A HERMITE CUBIC INTERPOLANT, AND THE CORRECTOR
C RETURNS TO THE ZERO CURVE USING A QUASI-NEWTON ALGORITHM, REMAINING
C IN A HYPERPLANE PERPENDICULAR TO THE MOST RECENT TANGENT VECTOR.
C  STEPQF  ALSO ESTIMATES A STEP SIZE  H  FOR THE NEXT STEP ALONG THE
C ZERO CURVE.
C
C
C ON INPUT:
C
C N = DIMENSION OF  X.
C
C NFE = NUMBER OF JACOBIAN MATRIX EVALUATIONS.
C
C IFLAG = -2, -1, OR 0, INDICATING THE PROBLEM TYPE.
C
C START = .TRUE. ON FIRST CALL TO  STEPQF, .FALSE. OTHERWISE.
C         SHOULD NOT BE MODIFIED BY THE USER AFTER THE FIRST CALL.
C
C HOLD = ||Y - YOLD|| ; SHOULD NOT BE MODIFIED BY THE USER.
C
C H = UPPER LIMIT ON LENGTH OF STEP THAT WILL BE ATTEMPTED.  H  MUST
C    BE SET TO A POSITIVE NUMBER ON THE FIRST CALL TO  STEPQF.
C    THEREAFTER,  STEPQF  CALCULATES AN OPTIMAL VALUE FOR  H, AND  H
C    SHOULD NOT BE MODIFIED BY THE USER.
C
C WK = APPROXIMATE CURVATURE FOR THE LAST STEP (COMPUTED BY PREVIOUS
C    CALL TO  STEPQF).  UNDEFINED ON FIRST CALL.  SHOULD NOT BE
C    MODIFIED BY THE USER.
C
C RELERR, ABSERR = RELATIVE AND ABSOLUTE ERROR VALUES.  THE ITERATION
C    IS CONSIDERED TO HAVE CONVERGED WHEN A POINT  Z=(LAMBDA,X)  IS
C    FOUND SUCH THAT
C       ||DZ|| .LE. RELERR*||Z|| + ABSERR,
C    WHERE  DZ  IS THE LAST QUASI-NEWTON STEP.
C
C S  = (APPROXIMATE) ARC LENGTH ALONG THE HOMOTOPY ZERO CURVE UP TO
C    Y(S) = (LAMBDA(S), X(S)).
C
C Y(1:N+1) = PREVIOUS POINT (LAMBDA(S),X(S)) FOUND ON THE ZERO CURVE
C    OF THE HOMOTOPY MAP.
C
C YP(1:N+1) = UNIT TANGENT VECTOR TO THE ZERO CURVE OF THE HOMOTOPY
C    MAP AT  Y.  INPUT IN THIS VECTOR IS NOT USED ON THE FIRST CALL
C    TO  STEPQF.
C
C YOLD(1:N+1) = A POINT BEFORE  Y  ON THE ZERO CURVE OF THE HOMOTOPY
C    MAP.  INPUT IN THIS VECTOR IS NOT USED ON THE FIRST CALL TO
C    STEPQF.
C
C YPOLD(1:N+1) = UNIT TANGENT VECTOR TO THE ZERO CURVE OF THE
C    HOMOTOPY MAP AT  YOLD.
C
C A(1:N) = PARAMETER VECTOR IN THE HOMOTOPY MAP.
C
C QT(1:N+1,1:N+1) = HOLDS  Q TRANSPOSE  OF THE QR FACTORIZATION OF
C    THE AUGMENTED JACOBIAN MATRIX AT  Y.
C
C R((N+1)*(N+2)/2) = HOLDS THE UPPER TRIANGLE OF  R  OF THE QR
C    FACTORIZATION, STORED BY ROWS.
C
C F0(1:N+1), F1(1:N+1), Z0(1:N+1), DZ(1:N+1), W(1:N+1), T(1:N+1) ARE
C    WORK ARRAYS.
C
C SSPAR(1:4) = PARAMETERS USED FOR COMPUTATION OF THE OPTIMAL STEP SIZE.
C    SSPAR(1) = HMIN, SSPAR(2) = HMAX, SSPAR(3) = BMIN, SSPAR(4) = BMAX.
C    THE OPTIMAL STEP  H  IS RESTRICTED SUCH THAT
C       HMIN .LE. H .LE. HMAX, AND  BMIN*HOLD .LE. H .LE. BMAX*HOLD.
C
C PAR(1:*)  AND  IPAR(1:*)  ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS,
C    WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES
C    RHO, RHOJAC.
C
C
C ON OUTPUT:
C
C NFE HAS BEEN UPDATED.
C
C IFLAG
C
C    = -2, -1, OR 0 (UNCHANGED) ON A NORMAL RETURN.
C
C    = 4 IF A JACOBIAN MATRIX WITH RANK <  N  HAS OCCURRED.  THE
C        ITERATION WAS NOT COMPLETED.
C
C    = 6 IF THE ITERATION FAILED TO CONVERGE.
C
C START = .FALSE. ON A NORMAL RETURN.
C
C CRASH
C
C    = .FALSE. ON A NORMAL RETURN.
C
C    = .TRUE. IF THE STEP SIZE  H  WAS TOO SMALL.  H  HAS BEEN
C      INCREASED TO AN ACCEPTABLE VALUE, WITH WHICH  STEPQF  MAY BE
C      CALLED AGAIN.
C
C    = .TRUE. IF  RELERR  AND/OR  ABSERR  WERE TOO SMALL.  THEY HAVE
C      BEEN INCREASED TO ACCEPTABLE VALUES, WITH WHICH  STEPQF  MAY
C      BE CALLED AGAIN.
C
C HOLD = ||Y-YOLD||.
C
C H = OPTIMAL VALUE FOR NEXT STEP TO BE ATTEMPTED.  NORMALLY  H  SHOULD
C     NOT BE MODIFIED BY THE USER.
C
C WK = APPROXIMATE CURVATURE FOR THE STEP TAKEN BY  STEPQF.
C
C S = (APPROXIMATE) ARC LENGTH ALONG THE ZERO CURVE OF THE HOMOTOPY
C     MAP UP TO THE LATEST POINT FOUND, WHICH IS RETURNED IN  Y.
C
C RELERR, ABSERR  ARE UNCHANGED ON A NORMAL RETURN.  THEY ARE POSSIBLY
C     CHANGED IF  CRASH  = .TRUE. (SEE DESCRIPTION OF  CRASH  ABOVE).
C
C Y, YP, YOLD, YPOLD  CONTAIN THE TWO MOST RECENT POINTS AND TANGENT
C     VECTORS FOUND ON THE ZERO CURVE OF THE HOMOTOPY MAP.
C
C QT, R  STORE THE QR FACTORIZATION OF THE AUGMENTED JACOBIAN MATRIX
C     EVALUATED AT  Y.
C
C
C CALLS  D1MACH, DAXPY, DCOPY, DDOT, DNRM2, DSCAL, F (OR RHO), FJAC
C     (OR RHOJAC), QRFAQF, QRSLQF, TANGQF, UPQRQF.
C
C ***** DECLARATIONS *****
C
C     FUNCTION DECLARATIONS
C
        DOUBLE PRECISION D1MACH, DDOT, DNRM2, QOFS
C
C     LOCAL VARIABLES
C
        DOUBLE PRECISION ALPHA, DD001, DD0011, DD01, DD011, DELS, ETA,
     $    FOURU, GAMMA, HFAIL, HTEMP, IDLERR, ONE, P0, P1, PP0, PP1,
     $    TEMP, TWOU, WKOLD
        INTEGER I, ITCNT, LITFH, J, JP1, NP1
        LOGICAL FAILED
C
C     SCALAR ARGUMENTS
C
        INTEGER N, NFE, IFLAG
        LOGICAL START, CRASH
        DOUBLE PRECISION HOLD, H, WK, RELERR, ABSERR, S
C
C     ARRAY DECLARATIONS
C
        DOUBLE PRECISION Y(N+1), YP(N+1), YOLD(N+1), YPOLD(N+1),
     $    A(N), QT(N+1,N+1), R((N+1)*(N+2)/2), F0(N+1), F1(N+1),
     $    Z0(N+1), DZ(N+1), W(N+1), T(N+1), SSPAR(4), PAR(1)
        INTEGER IPAR(1)
C
        SAVE
C
C ***** END OF DECLARATIONS *****
C
C DEFINITION OF HERMITE CUBIC INTERPOLANT VIA DIVIDED DIFFERENCES.
C
        DD01(P0,P1,DELS) = (P1-P0)/DELS
        DD001(P0,PP0,P1,DELS) = (DD01(P0,P1,DELS)-PP0)/DELS
        DD011(P0,P1,PP1,DELS) = (PP1-DD01(P0,P1,DELS))/DELS
        DD0011(P0,PP0,P1,PP1,DELS) = (DD011(P0,P1,PP1,DELS) -
     $                                DD001(P0,PP0,P1,DELS))/DELS
        QOFS(P0,PP0,P1,PP1,DELS,S) = ((DD0011(P0,PP0,P1,PP1,DELS)*
     $    (S-DELS) + DD001(P0,PP0,P1,DELS))*S + PP0)*S + P0
C
C ***** FIRST EXECUTABLE STATEMENT *****
C
C
C ***** INITIALIZATION *****
C
C ETA = PARAMETER FOR BROYDEN'S UPDATE.
C LITFH = MAXIMUM NUMBER OF QUASI-NEWTON ITERATIONS ALLOWED.
C
        ONE = 1.0
        TWOU = 2.0*D1MACH(4)
        FOURU = TWOU + TWOU
        NP1 = N+1
        FAILED = .FALSE.
        CRASH = .TRUE.
        ETA = 50.0*TWOU
        LITFH = 2*(INT(-LOG10(ABSERR+RELERR*DNRM2(NP1,Y,1)))+1)
C
C CHECK THAT ALL INPUT PARAMETERS ARE CORRECT.
C
C     THE ARCLENGTH  S MUST BE NONNEGATIVE.
C
        IF (S .LT. 0.0) RETURN
C
C     IF STEP SIZE IS TOO SMALL, DETERMINE AN ACCEPTABLE ONE.
C
        IF (H .LT. FOURU*(1.0+S)) THEN
          H=FOURU*(1.0 + S)
          RETURN
        END IF
C
C     IF ERROR TOLERANCES ARE TOO SMALL, INCREASE THEM TO ACCEPTABLE
C     VALUES.
C
        TEMP=DNRM2(NP1,Y,1) + 1.0
        IF (.5*(RELERR*TEMP+ABSERR) .LT. TWOU*TEMP) THEN
          IF (RELERR .NE. 0.0) THEN
            RELERR = FOURU*(1.0+FOURU)
            TEMP = 0.0
            ABSERR = MAX(ABSERR,TEMP)
          ELSE
            ABSERR=FOURU*TEMP
          END IF
          RETURN
        END IF
C
C     INPUT PARAMETERS WERE ALL ACCEPTABLE.
C
        CRASH = .FALSE.
C
C COMPUTE  YP  ON FIRST CALL.
C NOTE:  DZ  IS USED SIMPLY AS A WORK ARRAY HERE.
C
        IF (START) THEN
          CALL TANGQF(Y,YP,YPOLD,A,QT,R,W,DZ,T,N,IFLAG,NFE,PAR,IPAR)
          IF (IFLAG .GT. 0) RETURN
        END IF
C
C F0 = (RHO(Y), YP*Y) TRANSPOSE (DIFFERENT FOR EACH PROBLEM TYPE).
C
         IF (IFLAG .EQ. -2) THEN
C
C          CURVE TRACKING PROBLEM.
C
           CALL RHO(A,Y(1),Y(2),F0,PAR,IPAR)
         ELSE IF (IFLAG .EQ. -1) THEN
C
C          ZERO FINDING PROBLEM.
C
           CALL F(Y(2),F0)
           DO 5 I=1,N
             F0(I) = Y(1)*F0(I) + (1.0-Y(1))*(Y(I+1)-A(I))
  5        CONTINUE
         ELSE
C
C          FIXED POINT PROBLEM.
C
           CALL F(Y(2),F0)
           DO 10 I=1,N
             F0(I) = Y(1)*(A(I)-F0(I))+Y(I+1)-A(I)
  10       CONTINUE
         END IF
C
C        DEFINE LAST ROW OF F0  =  YP*Y.
C
           F0(NP1) = DDOT(NP1,YP,1,Y,1)
C
C ***** END OF INITIALIZATION *****
C
C ***** COMPUTE PREDICTOR POINT Z0 *****
C
  20    IF (START) THEN
C
C         COMPUTE Z0 WITH LINEAR PREDICTOR USING Y, YP --
C         Z0 = Y+H*YP.
C
          CALL DCOPY(NP1,Y,1,Z0,1)
          CALL DAXPY(NP1,H,YP,1,Z0,1)
C
        ELSE
C
C         COMPUTE Z0 WITH CUBIC PREDICTOR.
C
          DO 30 I=1,NP1
            Z0(I) = QOFS(YOLD(I),YPOLD(I),Y(I),YP(I),HOLD,HOLD+H)
  30      CONTINUE
C
        END IF
C
C F1 = (RHO(Z0), YP*Z0) TRANSPOSE.
C
        IF (IFLAG .EQ. -2) THEN
          CALL RHO(A,Z0(1),Z0(2),F1,PAR,IPAR)
        ELSE IF (IFLAG .EQ. -1) THEN
          CALL F(Z0(2),F1)
          DO 40 I=1,N
            F1(I) = Z0(1)*F1(I) + (1.0-Z0(1))*(Z0(I+1)-A(I))
  40      CONTINUE
        ELSE
          CALL F(Z0(2),F1)
          DO 50 I=1,N
            F1(I) = Z0(1)*(A(I)-F1(I))+Z0(I+1)-A(I)
  50      CONTINUE
        END IF
        F1(NP1) = DDOT(NP1,YP,1,Z0,1)
C
C ***** END OF PREDICTOR SECTION *****
C
C ***** SET-UP FOR QUASI-NEWTON ITERATION *****
C
        IF (FAILED) THEN
C
C GENERATE QT = AUGMENTED JACOBIAN MATRIX FOR POINT Z0=(LAMBDA,X).
C
          IF (IFLAG .EQ. -2) THEN
C
C           CURVE TRACKING PROBLEM:
C           D(RHO) = (D RHO(A,LAMBDA,X)/D LAMBDA, D RHO(A,LAMBDA,X)/DX).
C
            DO 60 J = 1,NP1
              CALL RHOJAC(A,Z0(1),Z0(2),QT(1,J),J,PAR,IPAR)
  60        CONTINUE
          ELSE IF (IFLAG .EQ. -1) THEN
C
C           ZERO FINDING PROBLEM:
C           D(RHO) = (F(X) - X + A, LAMBDA*DF(X) + (1-LAMBDA)*I).
C
            CALL F(Z0(2),QT(1,1))
            DO 70 I=1,N
              QT(I,1) = A(I) - Z0(I+1) + QT(I,1)
  70        CONTINUE
            DO 80 J= 1,N
              JP1 = J+1
              CALL FJAC(Z0(2),QT(1,JP1),J)
              CALL DSCAL(N, Z0(1), QT(1,JP1), 1)
              QT(J,JP1) = 1.0 - Z0(1) + QT(J,JP1)
  80        CONTINUE
          ELSE
C
C           FIXED POINT PROBLEM:
C           D(RHO) = (A - F(X), I - LAMBDA*DF(X)).
C
            CALL F(Z0(2),QT(1,1))
            CALL DSCAL(N,-ONE,QT(1,1),1)
            CALL DAXPY(N,ONE,A,1,QT(1,1),1)
            DO 90 J=1,N
              JP1 = J+1
              CALL FJAC(Z0(2),QT(1,JP1),J)
              CALL DSCAL(N, -Z0(1), QT(1,JP1), 1)
              QT(J,JP1) = 1.0 + QT(J,JP1)
  90        CONTINUE
          END IF
C
C       DEFINE LAST ROW OF QT = YP.
C
          CALL DCOPY(NP1, YP, 1, QT(NP1,1), NP1)
C
C       COUNT JACOBIAN EVALUATION.
C
          NFE = NFE+1
C
C DO FIRST QUASI NEWTON STEP.
C
C       FACTOR AUG.
C
          CALL QRFAQF(QT,R,NP1,IFLAG)
          IF (IFLAG .GT. 0) RETURN
C
C       COMPUTE NEWTON STEP.
C
          CALL DCOPY(N,F1,1,DZ,1)
          CALL DSCAL(N,-ONE,DZ,1)
          DZ(NP1) = 0.0
          CALL QRSLQF(QT,R,DZ,W,NP1)
C
C       TAKE STEP AND SET F0 = F1.
C
          CALL DAXPY(NP1, ONE, DZ, 1, Z0, 1)
          CALL DCOPY(NP1, F1, 1, F0, 1)
C
C       F1 = (RHO(Z0), YP*Z0) TRANSPOSE.
C
          IF (IFLAG .EQ. -2) THEN
            CALL RHO(A,Z0(1),Z0(2),F1,PAR,IPAR)
          ELSE IF (IFLAG .EQ. -1) THEN
            CALL F(Z0(2),F1)
            DO 100 I=1,N
              F1(I) = Z0(1)*F1(I) + (1.0-Z0(1))*(Z0(I+1)-A(I))
  100       CONTINUE
          ELSE
            CALL F(Z0(2),F1)
            DO 110 I=1,N
              F1(I) = Z0(1)*(A(I)-F1(I))+Z0(I+1)-A(I)
  110       CONTINUE
          END IF
          F1(NP1) = DDOT(NP1,YP,1,Z0,1)
C
        ELSE
C
C IF NOT FAILED THEN DEFINE  DZ=Z0-Y  PRIOR TO MAIN LOOP.
C
          CALL DCOPY(NP1,Z0,1,DZ,1)
          CALL DAXPY(NP1,-ONE,Y,1,DZ,1)
        END IF
C
C ***** END OF PREPARATION FOR QUASI-NEWTON ITERATION *****
C
C ***** QUASI-NEWTON ITERATION *****
C
        DO 140 ITCNT = 1,LITFH
C
C PERFORM UPDATE FOR NEWTON STEP JUST TAKEN.
C
          CALL UPQRQF(NP1,ETA,DZ,F0,F1,QT,R,W,T)
C
C COMPUTE NEXT NEWTON STEP.
C
          CALL DCOPY(N,F1,1,DZ,1)
          CALL DSCAL(N,-ONE,DZ,1)
          DZ(NP1) = 0.0
          CALL QRSLQF(QT,R,DZ,W,NP1)
C
C TAKE STEP.
C
          CALL DAXPY(NP1, ONE, DZ, 1, Z0, 1)
C
C CHECK FOR CONVERGENCE.
C
          IF (DNRM2(NP1,DZ,1) .LE. RELERR*DNRM2(NP1,Z0,1)+ABSERR) THEN
             GO TO 160
          END IF
C
C IF NOT CONVERGED, PREPARE FOR NEXT ITERATION.
C
C       F0 = F1.
C
          CALL DCOPY(NP1, F1, 1, F0, 1)
C
C       F1 = (RHO(Z0), YP*Z0) TRANSPOSE.
C
          IF (IFLAG .EQ. -2) THEN
            CALL RHO(A,Z0(1),Z0(2),F1,PAR,IPAR)
          ELSE IF (IFLAG .EQ. -1) THEN
            CALL F(Z0(2),F1)
            DO 120 I=1,N
              F1(I) = Z0(1)*F1(I) + (1.0-Z0(1))*(Z0(I+1)-A(I))
  120       CONTINUE
          ELSE
            CALL F(Z0(2),F1)
            DO 130 I=1,N
              F1(I) = Z0(1)*(A(I)-F1(I))+Z0(I+1)-A(I)
  130       CONTINUE
          END IF
          F1(NP1) = DDOT(NP1,YP,1,Z0,1)
C
  140   CONTINUE
C
C ***** END OF QUASI-NEWTON LOOP *****
C
C ***** DIDN'T CONVERGE OR TANGENT AT NEW POINT DID NOT MAKE
C       AN ACUTE ANGLE WITH YPOLD -- TRY AGAIN WITH A SMALLER H *****
C
  150   FAILED = .TRUE.
        HFAIL = H
        IF (H .LE. FOURU*(1.0 + S)) THEN
          IFLAG = 6
          RETURN
        ELSE
          H = .5 * H
        END IF
        GO TO 20
C
C ***** END OF CONVERGENCE FAILURE SECTION *****
C
C ***** CONVERGED -- MOP UP AND RETURN *****
C
C COMPUTE TANGENT & AUGMENTED JACOBIAN AT  Z0.
C NOTE:  DZ  AND  F1  ARE USED SIMPLY AS WORK ARRAYS HERE.
C
  160   CALL TANGQF(Z0,T,YP,A,QT,R,W,DZ,F1,N,IFLAG,NFE,PAR,IPAR)
        IF (IFLAG .GT. 0) RETURN
C
C CHECK THAT COMPUTED TANGENT  T  MAKES AN ANGLE NO LARGER THAN
C 60 DEGREES WITH CURRENT TANGENT  YP.  (I.E. COS OF ANGLE < .5)
C IF NOT, STEP SIZE WAS TOO LARGE, SO THROW AWAY Z0, AND TRY
C AGAIN WITH A SMALLER STEP.
C
        ALPHA = DDOT(NP1,T,1,YP,1)
        IF (ALPHA .LT. 0.5) GOTO 150
        ALPHA = ACOS(ALPHA)
C
C SET UP VARIABLES FOR NEXT CALL.
C
        CALL DCOPY(NP1,Y,1,YOLD,1)
        CALL DCOPY(NP1,Z0,1,Y,1)
        CALL DCOPY(NP1,YP,1,YPOLD,1)
        CALL DCOPY(NP1,T,1,YP,1)
C
C UPDATE ARCLENGTH   S = S + ||Y-YOLD||.
C
        HTEMP = HOLD
        CALL DAXPY(NP1,-ONE,YOLD,1,Z0,1)
        HOLD = DNRM2(NP1,Z0,1)
        S = S+HOLD
C
C COMPUTE OPTIMAL STEP SIZE.
C   IDLERR = DESIRED ERROR FOR NEXT PREDICTOR STEP.
C   WK = APPROXIMATE CURVATURE = 2*SIN(ALPHA/2)/HOLD  WHERE
C        ALPHA = ARCCOS(YP*YPOLD).
C   GAMMA = EXPECTED CURVATURE FOR NEXT STEP, COMPUTED BY
C        EXTRAPOLATING FROM CURRENT CURVATURE  WK, AND LAST
C        CURVATURE  WKOLD.  GAMMA IS FURTHER REQUIRED TO BE
C        POSITIVE.
C
        WKOLD = WK
        IDLERR = SQRT(SQRT(ABSERR + RELERR*DNRM2(NP1,Y,1)))
C
C     IDLERR SHOULD BE NO BIGGER THAN 1/2 PREVIOUS STEP.
C
        IDLERR = MIN(.5*HOLD,IDLERR)
        WK = 2.0*ABS(SIN(.5*ALPHA))/HOLD
        IF (START) THEN
           GAMMA = WK
        ELSE
           GAMMA = WK + HOLD/(HOLD+HTEMP)*(WK-WKOLD)
        END IF
        GAMMA = MAX(GAMMA, 0.01*ONE)
        H = SQRT(2.0*IDLERR/GAMMA)
C
C     ENFORCE RESTRICTIONS ON STEP SIZE SO AS TO ENSURE STABILITY.
C        HMIN <= H <= HMAX, BMIN*HOLD <= H <= BMAX*HOLD.
C
        H = MIN(MAX(SSPAR(1),SSPAR(3)*HOLD,H),SSPAR(4)*HOLD,SSPAR(2))
        IF (FAILED) H = MIN(HFAIL,H)
        START = .FALSE.
C
C ***** END OF MOP UP SECTION *****
C
        RETURN
C
C ***** END OF SUBROUTINE STEPQF *****
        END

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg


text:
Four Seasons Gourmet, formerly known as Gold Town Chinese Restaurant, is a Chinese food Chinese restaurant. They opened sometime in September 2013, replacing Gold Town.

As of September 2013 they are owned and operated by a new set of folks who dont appear to have any relationship with any of the previous owners, and no longer have Taiwanese specific menu offerings.

20130915 16:04:48 nbsp Im just excited they have sour spicy cobbage. It sounds even more delicious than sour spicy cabbage! Users/JaneBF

20131002 10:24:17 nbsp The mongolian chicken over rice is excellent. They also have dim sum offerings, but they take awhile to cook, so its best to call ahead.  Users/ChrisLambertus

20131124 14:00:07 nbsp I hit it up today. The place is kinda endearing with the typos. One on the sign that notify the problem to the chevron people if the bathroom needs work, and the other was MondayFried. Teehee.

That being said the food was pretty damn good. I perceived a massive portion of the general chicken. It was also all nice with its sticky candied goodness. BBQ pork fried rice was decent. I will report back when I try the walnut shrimp Users/StevenDaubert

reporting back on the walnut shrimp. The shrimp themselves were good, they split nicely and the sauce was slathered, could have been slightly more copious but it was good. The whole meal would have been great if it wasnt for the dealbreaker: 5 pieces of glazed walnuts total. 3 of those 5 pieces werent even full halves. How am I supposed to get some walnut on all 13 of the shrimp with only 5 partial pieces... What was even more irksome was the fact that the candied walnuts were fantastic and I really really desired more. ~SD
Post Script: There is a CC minimum of $10. This is slight miffing but I was intending on tipping a dollar anyways so it evened out on the $9 meal


text:
module test_filesystem
    use testsuite, only : new_unittest, unittest_t, error_t, test_failed
    use fpm_filesystem, only: canon_path, is_dir, mkdir, os_delete_dir, &
                              join_path
    use fpm_environment, only: OS_WINDOWS, get_os_type, os_is_unix
    implicit none
    private

    public :: collect_filesystem

contains

    !> Collect all exported unit tests
    subroutine collect_filesystem(testsuite)

        !> Collection of tests
        type(unittest_t), allocatable, intent(out) :: testsuite(:)

        testsuite = [ &
            & new_unittest("canon-path", test_canon_path), &
            & new_unittest("create-delete-directory", test_mkdir_rmdir) &
            ]

    end subroutine collect_filesystem

    subroutine test_canon_path(error)

        !> Error handling
        type(error_t), allocatable, intent(out) :: error

        call check_string(error, &
            & canon_path("git/project/src/origin"), "git/project/src/origin")
        if (allocated(error)) return

        call check_string(error,  &
            & canon_path("./project/src/origin"), "project/src/origin")
        if (allocated(error)) return

        call check_string(error, &
            & canon_path("./project/src///origin/"), "project/src/origin")
        if (allocated(error)) return

        call check_string(error, &
            & canon_path("../project/./src/origin/"), "../project/src/origin")
        if (allocated(error)) return

        call check_string(error, &
            & canon_path("/project//src/origin/"), "/project/src/origin")
        if (allocated(error)) return

        call check_string(error, &
            & canon_path("/project/src/../origin/"), "/project/origin")
        if (allocated(error)) return

        call check_string(error, &
            & canon_path("/project/src/../origin/.."), "/project")
        if (allocated(error)) return

        call check_string(error, &
            & canon_path("/project/src//../origin/."), "/project/origin")
        if (allocated(error)) return

        call check_string(error, &
            & canon_path("../project/src/./../origin/."), "../project/origin")
        if (allocated(error)) return

        call check_string(error, &
            & canon_path("../project/src/../../../origin/."), "../../origin")
        if (allocated(error)) return

        call check_string(error, &
            & canon_path("/../.."), "/")
        if (allocated(error)) return

        call check_string(error, &
            & canon_path("././././././/////a/b/.///././////.///c/../../../"), ".")
        if (allocated(error)) return

        call check_string(error, &
            & canon_path("/./././././/////a/b/.///././////.///c/../../../"), "/")
        if (allocated(error)) return

    end subroutine test_canon_path

    !> Check a character variable against a reference value
    subroutine check_string(error, actual, expected)

        !> Error handling
        type(error_t), allocatable, intent(out) :: error

        !> Actual string value
        character(len=*), intent(in) :: actual

        !> Expected string value
        character(len=*), intent(in) :: expected

        if (actual /= expected) then
            call test_failed(error, &
                "Character value mismatch "//&
                "expected '"//expected//"' but got '"//actual//"'")
        end if

    end subroutine check_string

    subroutine test_mkdir_rmdir(error)

        !> Error handling
        type(error_t), allocatable, intent(out) :: error

        call check_mkdir(error, join_path("tmpdir","subdir"))
        if (allocated(error)) return

        call check_rmdir(error, "tmpdir")
        if (allocated(error)) return

    end subroutine test_mkdir_rmdir

    !> Create a directory and verify its existence
    subroutine check_mkdir(error, path)

        !> Error handling
        type(error_t), allocatable, intent(out) :: error

        !> Directory path
        character(len=*), intent(in) :: path

        ! Directory shouldn't exist before it's created
        if (is_dir(path)) then
            call test_failed(error, &
                "Directory path "//path//" already exists before its creation")
            return
        end if

        ! Create directory
        call mkdir(path)

        ! Check that directory is indeed created
        if (.not.is_dir(path)) then
            call test_failed(error, &
                "Directory path "//path//" cannot be created")
        end if

      end subroutine check_mkdir

    !> Create a directory and verify its existence
    subroutine check_rmdir(error, path)

        !> Error handling
        type(error_t), allocatable, intent(out) :: error

        !> Directory path
        character(len=*), intent(in) :: path

        ! Directory should exist before it's deleted
        if (.not. is_dir(path)) then
            call test_failed(error, &
                "Directory path "//path//" doesn't exist before its deletion")
            return
        end if

        ! Delete directory
        call os_delete_dir(os_is_unix(),path)

        ! Check that directory is indeed deleted
        if (is_dir(path)) then
            call test_failed(error, &
                "Directory path "//path//" cannot be deleted")
        end if

      end subroutine check_rmdir

end module test_filesystem

max_stars_repo_name: freevryheid/fpm

max_issues_repo_name: freevryheid/fpm

max_forks_repo_name: freevryheid/fpm


! { dg-do compile }

SUBROUTINE check_dnucint_ana (dcore)
  IMPLICIT NONE
  INTEGER, PARAMETER :: dp=8
  REAL(dp), DIMENSION(10, 2), INTENT(IN),&
       OPTIONAL                            :: dcore
  INTEGER                                  :: i, j
  REAL(dp)                                 :: delta, nssss, od, rn, ssssm, &
       ssssp
  REAL(dp), DIMENSION(10, 2)               :: corem, corep, ncore
  LOGICAL                                  :: check_value

  delta = 1.0E-8_dp
  od = 0.5_dp/delta
  ncore = od * (corep - corem)
  nssss = od * (ssssp - ssssm)
  IF (PRESENT(dcore)) THEN
     DO i = 1, 2
        DO j = 1, 10
           IF (.NOT.check_value(ncore(j,i), dcore(j,i), delta, 0.1_dp)) THEN
           END IF
        END DO
     END DO
  END IF
END SUBROUTINE check_dnucint_ana


text: The California Animal Health and Food Safety Lab is a cooperative effort between UC Davis School of Veterinary Medicine, The California Department of Food and Agriculture, and The California Animal Agriculture Industries. Both Maddy Lab and Thurman Laboratory Thurman Lab are part of the Lab.


text:
subroutine Hilm(cilm, ba, gridglq, lmax, nmax, mass, r0, rho, w, plx, zero, filter_type, filter_deg)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!	This routine will compute the next estimate of Moho coefficients
!	given an initial estimate of the moho in a gridded data file. This
!	is simply Equation 18 in Wieczorek and Phillips (1998). Note that the
!	degree-0 topography must be included in the gridded relief. Note that
!	the array plx is optional, and should not be precomputed when memory
!	is an issue (i.e., lmax>360).
!
!	Calling Parameters:
!		IN
!			ba:		Bouguer Anomaly spherical harmonic coefficients.
!			gridglq:	Initial estimate of Moho refief, gridded according
!					to a call to MakeGridGLQ.
!			lmax:		Maxmimum spherical harmonic degree to compute.
!			nmax:		Order of potential coefficient expansion.
!			mass:		Mass of planet.
!			R0:		The radius that the coefficients are referenced
!					to.
!			rho:		Density contrast between the mantle and crust.
!			w:		Gauss-Legendre points used in integrations
!					(determined from a call to PreCompute).
!		OUT
!			cilm:		Estimate of Moho relief spherical harmonic coefficients,
!					with dimensions (2, lmax+1, lmax+1).
!
!		OPTIONAL
!			filter_deg:	If specified, each interation will be filtered according to
!					equation 18 in Wieczorek and Phillips (1998), where the value
!					of filter_deg corresponds to the spherical harmonic degree where
!					the filter is 1/2.
!			filter_type:	If filter_deg is specified, this must be as well.
!					A value of (1) corresponds to the minimum amplitude filter in
!					Wieczorek and Phillips (1998), whereas as value of (2) corresponds
!					to a minimum curvature filter.
!			plx:		Input array of Associated Legendre Polnomials computed
!					at the Gauss points (determined from a call to
!					PreCompute). If this is not included, then the optional
!					array zero MUST be inlcuded.
!			zero		Array of dimension lmax+1 that contains the latitudinal
!					gridpoints used in the Gauss-Legendre quadrature integration
!					scheme. Only needed if plx is not included.
!
!	All units assumed to be SI.
!
!	Dependencies:		NGLQSH, SHExpandGLQ, wl, wl_curv
!
!	Written by Mark Wieczorek 2003
!	September 3, 2005. Modifed so that the array plx is now optional.
!
!	Copyright (c) 2005, Mark A. Wieczorek
!	All rights reserved.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	use SHTOOLS, only: NGLQSH, SHExpandGLQ, Wl, WlCurv

	implicit none

	real*8, intent(out) :: 			cilm(:,:,:)
	real*8, intent(in) ::			ba(:,:,:), gridglq(:,:), mass, r0, rho, w(:)
	real*8, intent(in), optional ::		plx(:,:), zero(:)
	integer, intent(in) ::			lmax, nmax
	integer, intent(in), optional :: 	filter_type, filter_deg
	real*8 ::				prod, pi, d, depth, filter(lmax+1)
	real*8, allocatable ::			cilmn(:, :, :), grid2(:,:)
	integer  ::				j, l, n, nlong, nlat, astat(2)

	if (size(cilm(:,1,1)) < 2 .or. size(cilm(1,:,1)) < lmax+1 .or. size(cilm(1,1,:)) < lmax+1) then
		print*, "Error --- HilmGLQ"
		print*, "CILM must be dimensioned as (2, LMAX+1, LMAX+1) where LMAX is ", lmax
		print*, "Input dimension is ", size(cilm(:,1,1)), size(cilm(1,:,1)), size(cilm(1,1,:))
		stop
	elseif (size(ba(:,1,1)) < 2 .or. size(ba(1,:,1)) < lmax+1 .or. size(ba(1,1,:)) < lmax+1) then
		print*, "Error --- HilmGLQ"
		print*, "BA must be dimensioned as (2, LMAX+1, LMAX+1) where LMAX is ", lmax
		print*, "Input dimension is ", size(ba(:,1,1)), size(ba(1,:,1)), size(ba(1,1,:))
		stop
	elseif(size(gridglq(1,:)) < 2*lmax+1 .or. size(gridglq(:,1)) <  lmax +1) then
		print*, "Error --- HilmGLQ"
		print*, "GRIDGLQ must be dimensioned as (LMAX+1, 2*LMAX+1) where LMAX is ", lmax
		print*, "Input dimension is ", size(gridglq(1,:)), size(gridglq(:,1))
		stop
	elseif(size(w) < lmax+1) then
		print*, "Error --- HilmGLQ"
		print*, "W must be dimensioned as (LMAX+1) where LMAX is ", lmax
		print*, "Input dimension is ", size(w)
		stop
	endif

	if (present(zero)) then
		if (size(zero) < lmax + 1) then
			print*, "Error --- HilmGLQ"
			print*, "ZERO must be dimensioned as (LMAX+1) where LMAX is ", lmax
			print*, "Input dimension is ", size(zero)
			stop
		endif
	endif

	if (present(plx)) then
		if (size(plx(:,1)) < lmax+1 .or. size(plx(1,:)) < (lmax+1)*(lmax+2)/2) then
			print*, "Error --- HilmGLQ"
			print*, "PLX must be dimensioned as (LMAX+1, (LMAX+1)*(LMAX+2)/2) where LMAX is ", lmax
			print*, "Input dimension is ", size(plx(:,1)), size(plx(1,:))
			stop
		endif
	endif

	if(present(filter_type)) then
		if (filter_type /=0 .and. filter_type /=1 .and. filter_type /=2) then
			print*, "Error --- HilmGLQ"
			print*, "FILTER_TYPE must be either 0 (none), 1 (minimum amplitude), or 2 (minimum curvature)"
			print*, "Input value is ", filter_type
			stop
		endif
	endif

	if (present(plx)) then
		continue
	else
		if (present(zero)) then
			continue
		else
			print*, "Error --- Hilm"
			print*, "If the optional array PLX is not included, it is necessary to include the optional array ZERO."
			stop
		endif
	endif

	allocate(cilmn(2, lmax+1, lmax+1), stat = astat(1))
	allocate(grid2(lmax+1,2*lmax+1), stat = astat(2))
	if (astat(1) /= 0 .or. astat(2) /= 0) then
		print*, "Error --- HilmGLQ"
		print*, "Problem allocating arrays CILMN and GRID2", astat(1), astat(2)
		stop
	endif

	cilm = 0.0d0
	cilmn = 0.0d0
	pi = acos(-1.0d0)
	nlat = NGLQSH(lmax)
	nlong = 2*lmax+1
	grid2(1:nlat,1:nlong) = gridglq(1:nlat,1:nlong)

	if (present(plx)) then
		call SHExpandGLQ(cilmn, lmax, grid2(1:nlat,1:nlong), w(1:lmax+1), plx = plx(1:lmax+1, 1:(lmax+1)*(lmax+2)/2), &
			norm = 1, csphase = 1)
	else
		call SHExpandGLQ(cilmn, lmax, grid2(1:nlat,1:nlong), w(1:lmax+1), zero = zero(1:lmax+1), norm = 1, csphase = 1)
	endif

	d = cilmn(1,1,1)
	depth = r0-d
	print*, "Average depth of Moho (km) = ", depth/1.d3

	cilm(1,1,1) = d
	grid2(1:nlat,1:nlong) = grid2(1:nlat,1:nlong) - d

	filter = 1.0d0
	if (present(filter_type) .and. present(filter_deg)) then
		do l=1, lmax
			if (filter_type ==1) then
				filter(l+1) = Wl(l, filter_deg, r0, d)
			elseif(filter_type==2) then
				filter(l+1) = WlCurv(l, filter_deg, r0, d)
			endif
		enddo
	endif

	do l=1, lmax
		cilm(:,l+1,1:l+1) = filter(l+1)*ba(:,l+1,1:l+1) * mass * dble(2*l+1) * ( (r0/d)**l) / &
			(4.0d0*pi*rho*d**2)
	enddo

	do n=2, nmax

		if (present(plx)) then
			call SHExpandGLQ(cilmn, lmax, (grid2(1:nlat,1:nlong)/d)**n, w(1:lmax+1), &
				plx = plx(1:lmax+1, 1:(lmax+1)*(lmax+2)/2), norm = 1, csphase = 1)
		else
			call SHExpandGLQ(cilmn, lmax, (grid2(1:nlat,1:nlong)/d)**n, w(1:lmax+1), zero = zero(1:lmax+1), &
				norm = 1, csphase = 1)
		endif

		do l = 1, lmax
			prod = 1.0d0
			do j=1, n
				prod = prod * dble(l+4-j)
			enddo
			prod = d * prod/( dble(l+3) * dble(fact(n)) )

			cilm(:,l+1,1:l+1) = cilm(:,l+1,1:l+1) - filter(l+1)*cilmn(:,l+1,1:l+1)*prod
		enddo
	enddo

	deallocate(cilmn)
	deallocate(grid2)

	contains

		function fact(i)
		!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
		!
		!	This function computes the factorial of an integer
		!
		!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
			implicit none
			integer ::	i, j, fact

			if (i == 0) then
				fact = 1
			elseif (i .lt. 0) then
				print*, "Argument to FACT must be positive"
				stop
			else
				fact = 1
				do j = 1, i
					fact = fact * j
				enddo
			endif

		end function fact

end subroutine Hilm

max_stars_repo_name: pletnes/SHTOOLS

max_issues_repo_name: pletnes/SHTOOLS

max_forks_repo_name: pletnes/SHTOOLS