text:
c
c This a driver to test fastlap V2.0 as a callable procedure from fortran
c for computing a desingularized solution of a single time step of
c Steve Scorpio's free surface problem.
c Input is from Steve's dumped data.
c
c NB you can precondition these problems with OL which has been extended to
c handle non-square local problems, but it is not very effective.  Use the SP
c preconditioner which was designed for the desingularized case.
c
c Written in Modern Fortran by Korsmeyer, winter 1993.  (c) MIT 1993.
c
      implicit none
      integer bsize, fssize, size
      parameter (bsize=510, fssize=2486)
      parameter (size=bsize+fssize)
      double precision a(size,size), r, rx, ry, rz
      integer ipvt(size), info, idum
      double precision x(3,4,size), xnrm(3,size),
     $     xcoll(3,size), lhsvect(size), rhsvect(size)
      integer type(size), shape(size), dtype(size), job, fljob
      integer rhstype(size), lhstype(size), rhsindex(4,size),
     $     lhsindex(4,size)
      integer NULL, POINT_SOURCE, CONSTANT_SOURCE, CONSTANT_DIPOLE,
     $     LINEAR_SOURCE, LINEAR_DIPOLE
      data NULL, POINT_SOURCE, CONSTANT_SOURCE, CONSTANT_DIPOLE,
     $     LINEAR_SOURCE, LINEAR_DIPOLE /0, 1, 11, 12, 21, 22/
      double precision scorp_sol(size)
      integer nummom, numlev, i, j, k, nlhs, nrhs, nsing
      double precision error, max_a_error, ave_a_error,
     $     max_r_error, ave_r_error
      character*1 dumchar
      integer fastlap
      external fastlap
      double precision tol
      integer numitr, maxitr
      data maxitr /1024/
c
c Query for run parameters and data file name.
c
c
      write(*,'(//''We solve a single time step of the Scorpio''/
     $     ''free-sruface problem  '')')
      write(*,'(//''You should use the SP preconditioner for''/
     $     ''this problem.  Set PRECOND = SP in mulGlobal.h. '')')
      write(*,'(''Select expansion order and tree depth:  ''$)')
      read(*,*) nummom, numlev
      write(*,'(''Select the iteration tolerance:  ''$)')
      read(*,*) tol
c
c Read desingularized source, collocation points, etc.
c (There are 510 body points and 2486 free-surface points.)
c
      open(unit=20,file='bdbvp-1.dat',status='old')
      open(unit=21,file='bdbvp-2.dat',status='old')
      open(unit=22,file='fsbvp.dat',status='old')
c
c Get past the headers.
c
      do i = 1,3
         read(20,'(a)') dumchar
         read(21,'(a)') dumchar
         read(22,'(a)') dumchar
      end do
      do i = 1,4
         read(20,'(a)') dumchar
         read(22,'(a)') dumchar
      end do
      print *, ' Reading',bsize+fssize,' source and field points.'
      nsing = size
      nlhs = nsing
      nrhs = nsing
      do i = 1, bsize
         read(20,*) (xcoll(k,i),k=1,3),(x(k,1,i),k=1,3),
     $        rhsvect(i),scorp_sol(i)
         rhsindex(1,i) = i
         lhsindex(1,i) = i
         shape(i) = 1
         dtype(i) = 1
         rhstype(i) = POINT_SOURCE
         lhstype(i) = POINT_SOURCE
         read(21,*) (xnrm(k,i),k=1,3)
      end do
      do i = bsize+1, size
         read(22,*) (xcoll(k,i),k=1,3),(x(k,1,i),k=1,3),
     $        rhsvect(i),scorp_sol(i)
         rhsindex(1,i) = i
         lhsindex(1,i) = i
         shape(i) = 1
         dtype(i) = 0
         rhstype(i) = POINT_SOURCE
         lhstype(i) = POINT_SOURCE
      end do
c
c Call fastlap to compute a solution or field at the vector of field points.
c
      fljob = 2
      print *, "Calling fastlap."
      numitr = fastlap(nlhs,nrhs,nsing,x,shape,dtype,lhstype,rhstype,
     $     lhsindex,rhsindex,lhsvect,rhsvect,xcoll,xnrm,numlev,nummom,
     $     maxitr,tol,fljob)
      write(*,'(//i3,
     $     '' iterations knocked down residual to:'',e14.8)')
     $     numitr, tol
c
c Compute the solution by LU, and call it exact.  Or possibly
c go around this block and read a previous solution.
c
      go to 100
      do j = 1,size
         do i = 1,bsize
            r = (xcoll(1,i) - x(1,1,j))**2
     $           + (xcoll(2,i) - x(2,1,j))**2
     $           + (xcoll(3,i) - x(3,1,j))**2
            rx = -(xcoll(1,i) - x(1,1,j)) / r**1.5
            ry = -(xcoll(2,i) - x(2,1,j)) / r**1.5
            rz = -(xcoll(3,i) - x(3,1,j)) / r**1.5
            a(i,j) = xnrm(1,i) * rx + xnrm(2,i) * ry + xnrm(3,i) * rz
         end do
         do i = bsize+1,size
            r = (xcoll(1,i) - x(1,1,j))**2
     $           + (xcoll(2,i) - x(2,1,j))**2
     $           + (xcoll(3,i) - x(3,1,j))**2
            a(i,j) = 1./sqrt(r)
         end do
      end do
      call DGETRF(size,size,a,size,ipvt,INFO)
      call DGETRS('N',size,1,a,size,ipvt,
     $     rhsvect,size,INFO)
 100  continue
      open(unit=24,file='oldLU.dat',status='old',form='formatted')
      do i = 1,size
         read(24,*) idum, rhsvect(i)
      end do
      close(24)
c
c Compute the error in the multi and the scorp solution.
c

      ave_a_error = 0.
      ave_r_error = 0.
      max_a_error = 0.
      max_r_error = 0.
      open(unit=25,file='solutions.dat',
     $     status='unknown',form='formatted')
      write(25,'(''Index   LU solution   Dom-Decomp   fastlap'')')
      do i = 1,nlhs
         write(25,999) i,rhsvect(i),scorp_sol(i),lhsvect(i)
         error = abs(rhsvect(i) - lhsvect(i))
         write(98,999) i,error,error/rhsvect(i)
         ave_a_error = ave_a_error + error
         ave_r_error = ave_r_error + error/rhsvect(i)
         max_a_error = max(max_a_error,error)
         max_r_error = max(max_r_error,error/rhsvect(i))
         if (error/rhsvect(i) .gt. 100) write(88,*) 'multi:  ',
     $        i, rhsvect(i), lhsvect(i)
      end do
 999  format(i6,3e20.10)
      ave_a_error = ave_a_error / float(nlhs)
      ave_r_error = ave_r_error / float(nlhs)
      write(*,'(//''Average absolute multi error ='',
     $     f13.8)') ave_a_error
      write(*,'(''Maximum absolute multi error ='',
     $     f13.8)') max_a_error
      write(*,'(//''Average relative multi error ='',
     $     f13.8)') ave_r_error
      write(*,'(''Maximum relative multi error ='',
     $     f13.8)') max_r_error
      ave_a_error = 0.
      ave_r_error = 0.
      max_a_error = 0.
      max_r_error = 0.
      do i = 1,nlhs
         error = abs(rhsvect(i) - scorp_sol(i))
         write(97,999) i,error,error/rhsvect(i)
         ave_a_error = ave_a_error + error
         ave_r_error = ave_r_error + error/rhsvect(i)
         max_a_error = max(max_a_error,error)
         max_r_error = max(max_r_error,error/rhsvect(i))
         if (error/rhsvect(i) .gt. 100) write(88,*) 'scorp:  ',
     $        i, rhsvect(i), scorp_sol(i)
      end do
      ave_a_error = ave_a_error / float(nlhs)
      ave_r_error = ave_r_error / float(nlhs)
      write(*,'(//''Average absolute scorpio error ='',
     $     f13.8)') ave_a_error
      write(*,'(''Maximum absolute scorpio error ='',
     $     f13.8)') max_a_error
      write(*,'(//''Average relative scorpio error ='',
     $     f13.8)') ave_r_error
      write(*,'(''Maximum relative scorpio error ='',
     $     f13.8)') max_r_error
      stop
      end

max_stars_repo_name: OregonIons/bem

max_issues_repo_name: OregonIons/bem

max_forks_repo_name: OregonIons/bem


text:
     SUBROUTINE SYMTRY
     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
     INCLUDE 'SIZES/NOLIST'
     COMMON /GMETRY/ GEO(3,NUMATM)
     COMMON /GEOSYM/ NDEP, LOCPAR(200), IDEPFN(200), LOCDEP(200)
**********************************************************************

  SYMTRY COMPUTES THE BOND LENGTHS AND ANGLES THAT ARE FUNCTIONS OF
         OTHER BOND LENGTHS AND ANGLES.

 ON INPUT GEO     = KNOWN INTERNAL COORDINATES
          NDEP    = NUMBER OF DEPENDENCY FUNCTIONS.
          IDEPFN  = ARRAY OF DEPENDENCY FUNCTIONS.
          LOCDEP  = ARRAY OF LABELS OF DEPENDENT ATOMS.
          LOCPAR  = ARRAY OF LABELS OF REFERENCE ATOMS.

  ON OUTPUT THE ARRAY "GEO" IS FILLED
***********************************************************************

     NOW COMPUTE THE DEPENDENT PARAMETERS.

     DO 30 I=1,NDEP
        CALL HADDON (VALUE,LOCN,IDEPFN(I),LOCPAR(I),GEO)
        J=LOCDEP(I)
  30 GEO(LOCN,J)=VALUE
     RETURN
     END

ext: for

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive


text:
-makelib ies/dist_mem_gen_v8_0_9 \
  "../../../ipstatic/dist_mem_gen_v8_0_9/simulation/dist_mem_gen_v8_0.vhd" \
-endlib
-makelib ies/xil_defaultlib \
  "../../../../Pipelining.srcs/sources_1/ip/InstructionROM_1/sim/InstructionROM.vhd" \
-endlib

ext: f

max_stars_repo_name: CrossingX/LPCO

max_issues_repo_name: CrossingX/LPCO

max_forks_repo_name: CrossingX/LPCO


text:
      SUBROUTINE INTEGRATE( TIN, TOUT )

      IMPLICIT KPP_REAL (A-H,O-Z)
      INCLUDE 'KPP_ROOT_params.h'
      INCLUDE 'KPP_ROOT_global.h'

C TIN - Start Time
      KPP_REAL TIN
C TOUT - End Time
      KPP_REAL TOUT
      INTEGER i

      PARAMETER (LWORK=2*NVAR*NVAR+14*NVAR+20,LIWORK=NVAR+20)
      KPP_REAL WORK(LWORK)
      INTEGER IWORK(LIWORK)
      EXTERNAL FUNC_CHEM,JAC_CHEM

      ITOL=1     ! --- VECTOR TOLERANCES
      IJAC=1     ! --- COMPUTE THE JACOBIAN ANALYTICALLY
      MLJAC=NVAR ! --- JACOBIAN IS A FULL MATRIX
      MUJAC=NVAR ! --- JACOBIAN IS A FULL MATRIX
      IMAS=0     ! --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
      IOUT=0     ! --- OUTPUT ROUTINE IS NOT USED DURING INTEGRATION
      IDFX=0     ! --- INTERNAL TIME DERIVATIVE

      DO i=1,20
        IWORK(i) = 0
        WORK(i) = 0.D0
      ENDDO

      IWORK(3) = 1

      CALL ATMRODAS(NVAR,FUNC_CHEM,Autonomous,TIN,VAR,TOUT,
     &                  STEPMIN,RTOL,ATOL,ITOL,
     &                  JAC_CHEM,IJAC,MLJAC,MUJAC,FUNC_CHEM,IDFX,
     &                  FUNC_CHEM,IMAS,
     &                  WORK,LWORK,IWORK,LIWORK,IDID)

      IF (IDID.LT.0) THEN
        print *,'ATMRODAS: Unsucessfull exit at T=',
     &          TIN,' (IDID=',IDID,')'
      ENDIF

      RETURN
      END

      SUBROUTINE ATMRODAS(N,FCN,IFCN,X,Y,XEND,H,
     &                  RelTol,AbsTol,ITOL,
     &                  JAC ,IJAC,MLJAC,MUJAC,DFX,IDFX,
     &                  MAS ,IMAS,
     &                  WORK,LWORK,IWORK,LIWORK,IDID)
C ----------------------------------------------------------
C     NUMERICAL SOLUTION OF A STIFF (OR DIFFERENTIAL ALGEBRAIC)
C     SYSTEM OF FIRST 0RDER ORDINARY DIFFERENTIAL EQUATIONS  MY'=F(X,Y).
C     THIS IS AN EMBEDDED ROSENBROCK METHOD OF ORDER (3)4
C     (WITH STEP SIZE CONTROL).
C     C.F. SECTIONS IV.7  AND VI.3
C
C     AUTHORS: E. HAIRER AND G. WANNER
C              UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES
C              CH-1211 GENEVE 24, SWITZERLAND
C              E-MAIL:  HAIRER@DIVSUN.UNIGE.CH,  WANNER@DIVSUN.UNIGE.CH
C ---------------------------------------------------------
C *** *** *** *** *** *** *** *** *** *** *** *** ***
C          DECLARATIONS
C *** *** *** *** *** *** *** *** *** *** *** *** ***
      IMPLICIT KPP_REAL (A-H,O-Z)
      DIMENSION Y(N),AbsTol(*),RelTol(*),WORK(LWORK),IWORK(LIWORK)
      LOGICAL AUTNMS,IMPLCT,JBAND,ARRET,PRED
      EXTERNAL FCN,JAC,DFX,MAS
      COMMON /STATISTICS/ NFCN,NACCPT,NREJCT,NSTEP,NJAC,NDEC,NSOL
C *** *** *** *** *** *** ***
C        SETTING THE PARAMETERS
C *** *** *** *** *** *** ***
      ARRET=.FALSE.
      METH = 1
      NMAX=100000
C -------- PRED   STEP SIZE CONTROL
      IF(IWORK(3).LE.1)THEN
         PRED=.TRUE.
      ELSE
         PRED=.FALSE.
      END IF
      UROUND=1.D-16
      NM1 = N
      M1  = N
      M2  = N
C -------- MAXIMAL STEP SIZE
      IF(WORK(2).EQ.0.D0)THEN
         HMAX=XEND-X
      ELSE
         HMAX=WORK(2)
      END IF
C -------  FAC1,FAC2     PARAMETERS FOR STEP SIZE SELECTION
      IF(WORK(3).EQ.0.D0)THEN
         FAC1=5.D0
      ELSE
         FAC1=1.D0/WORK(3)
      END IF
      IF(WORK(4).EQ.0.D0)THEN
         FAC2=1.D0/6.0D0
      ELSE
         FAC2=1.D0/WORK(4)
      END IF
      IF (FAC1.LT.1.0D0.OR.FAC2.GT.1.0D0) THEN
            WRITE(6,*)' CURIOUS INPUT WORK(3,4)=',WORK(3),WORK(4)
            ARRET=.TRUE.
         END IF
C --------- SAFE     SAFETY FACTOR IN STEP SIZE PREDICTION
      SAFE=0.9D0
C --------- CHECK IF TOLERANCES ARE O.K.
      IF (ITOL.EQ.0) THEN
          IF (AbsTol(1).LE.0.D0.OR.RelTol(1).LE.10.D0*UROUND) THEN
              WRITE (6,*) ' TOLERANCES ARE TOO SMALL'
              ARRET=.TRUE.
          END IF
      ELSE
          DO I=1,N
             IF (AbsTol(I).LE.0.D0.OR.RelTol(I).LE.10.D0*UROUND) THEN
                WRITE (6,*) ' TOLERANCES(',I,') ARE TOO SMALL'
                ARRET=.TRUE.
             END IF
          END DO
      END IF

      IF (ARRET) STOP
      NM1 = N
C *** *** *** *** *** *** *** *** *** *** *** *** ***
C         COMPUTATION OF ARRAY ENTRIES
C *** *** *** *** *** *** *** *** *** *** *** *** ***
C ---- AUTONOMOUS, IMPLICIT, BANDED OR NOT ?
      AUTNMS=IFCN.EQ.0
      IMPLCT=IMAS.NE.0
      JBAND=MLJAC.LT.NM1
C -------- COMPUTATION OF THE ROW-DIMENSIONS OF THE 2-ARRAYS ---
C -- JACOBIAN AND MATRIX E
         MLJAC=NM1
         MUJAC=NM1
         LDJAC=NM1
         LDE=NM1
C -- MASS MATRIX
      IF (IMPLCT) THEN
         print *, 'Implicit 1'
      ELSE
         LDMAS=0
         IJOB=1
      END IF
      LDMAS2=MAX(1,LDMAS)
C ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK -----
      IEYNEW=21
      IEDY1=IEYNEW+N
      IEDY=IEDY1+N
      IEAK1=IEDY+N
      IEAK2=IEAK1+N
      IEAK3=IEAK2+N
      IEAK4=IEAK3+N
      IEAK5=IEAK4+N
      IEAK6=IEAK5+N
      IEFX =IEAK6+N
      IECON=IEFX+N
      IEJAC=IECON+4*N
      IEMAS=IEJAC+N*LDJAC
      IEE  =IEMAS+NM1*LDMAS
C ------ TOTAL STORAGE REQUIREMENT -----------
      ISTORE=IEE+NM1*LDE-1
      IF(ISTORE.GT.LWORK)THEN
         WRITE(6,*)' INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=',ISTORE
         ARRET=.TRUE.
      END IF
C ------- ENTRY POINTS FOR INTEGER WORKSPACE -----
      IEIP=21
      ISTORE=IEIP+NM1-1
      IF(ISTORE.GT.LIWORK)THEN
         WRITE(6,*)' INSUFF. STORAGE FOR IWORK, MIN. LIWORK=',ISTORE
         ARRET=.TRUE.
      END IF
C ------ WHEN A FAIL HAS OCCURED, WE RETURN WITH IDID=-1
      IF (ARRET) THEN
         IDID=-1
         RETURN
      END IF
C -------- CALL TO CORE INTEGRATOR ------------
      CALL ROSCOR(N,FCN,X,Y,XEND,HMAX,H,RelTol,AbsTol,ITOL,JAC,IJAC,
     &   MLJAC,MUJAC,DFX,IDFX,MAS,MLMAS,MUMAS,IOUT,IDID,NMAX,
     &   UROUND,METH,IJOB,FAC1,FAC2,SAFE,AUTNMS,IMPLCT,JBAND,PRED,LDJAC,
     &   LDE,LDMAS2,WORK(IEYNEW),WORK(IEDY1),WORK(IEDY),WORK(IEAK1),
     &   WORK(IEAK2),WORK(IEAK3),WORK(IEAK4),WORK(IEAK5),WORK(IEAK6),
     &   WORK(IEFX),WORK(IEJAC),WORK(IEE),WORK(IEMAS),IWORK(IEIP),
     &   WORK(IECON),
     &   M1,M2,NM1,NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL)
      IWORK(14)=NFCN
      IWORK(15)=NJAC
      IWORK(16)=NSTEP
      IWORK(17)=NACCPT
      IWORK(18)=NREJCT
      IWORK(19)=NDEC
      IWORK(20)=NSOL
C ----------- RETURN -----------
      RETURN
      END
C
C
C
C  ----- ... AND HERE IS THE CORE INTEGRATOR  ----------
C
      SUBROUTINE ROSCOR(N,FCN,X,Y,XEND,HMAX,H,RelTol,AbsTol,
     &  ITOL,JAC,IJAC,
     &  MLJAC,MUJAC,DFX,IDFX,MAS,MLMAS,MUMAS,IOUT,IDID,NMAX,
     &  UROUND,METH,IJOB,FAC1,FAC2,SAFE,AUTNMS,IMPLCT,BANDED,
     &  PRED,LDJAC,
     &  LDE,LDMAS,YNEW,DY1,DY,AK1,AK2,AK3,AK4,AK5,AK6,
     &  FX,FJAC,E,FMAS,IP,CONT,
     &  M1,M2,NM1,NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL)
C ----------------------------------------------------------
C     CORE INTEGRATOR FOR RODAS
C     PARAMETERS SAME AS IN RODAS WITH WORKSPACE ADDED
C ----------------------------------------------------------
C         DECLARATIONS
C ----------------------------------------------------------
      IMPLICIT KPP_REAL (A-H,O-Z)
      INCLUDE 'KPP_ROOT_params.h'
      INCLUDE 'KPP_ROOT_sparse.h'
      DIMENSION Y(N),YNEW(N),DY1(N),DY(N),AK1(N),
     *  AK2(N),AK3(N),AK4(N),AK5(N),AK6(N),FX(N),
     *  FJAC(LU_NONZERO),E(LDE,NM1),FMAS(LDMAS,NM1),
     *  AbsTol(*),RelTol(*)
      DIMENSION CONT(4*N)
      INTEGER IP(NM1)
      LOGICAL REJECT,AUTNMS,IMPLCT,BANDED
      LOGICAL ONE,LAST,PRED,SINGULAR
      EXTERNAL FCN, MAS, JAC, DFX
      COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
      COMMON /CONROS/XOLD,HOUT,NN
C *** *** *** *** *** *** ***
C  INITIALISATIONS
C *** *** *** *** *** *** ***
      NN=N
      NN2=2*N
      NN3=3*N
      LRC=4*N
C ------- COMPUTE MASS MATRIX FOR IMPLICIT CASE ----------
      IF (IMPLCT) CALL MAS (NM1,FMAS,LDMAS)
C ------ SET THE PARAMETERS OF THE METHOD -----
      CALL ROCOE(METH,A21,A31,A32,A41,A42,A43,A51,A52,A53,A54,
     &    C21,C31,C32,C41,C42,C43,C51,C52,C53,C54,C61,
     &    C62,C63,C64,C65,GAMMA,C2,C3,C4,D1,D2,D3,D4,
     &    D21,D22,D23,D24,D25,D31,D32,D33,D34,D35)
C --- INITIAL PREPARATIONS
      IF (M1.GT.0) IJOB=IJOB+10
      POSNEG=SIGN(1.D0,XEND-X)
      HMAXN=DMIN1(DABS(HMAX),DABS(XEND-X))
      IF (DABS(H).LE.10.D0*UROUND) H=1.0D-6
      H=DMIN1(DABS(H),HMAXN)
      H=SIGN(H,POSNEG)
      HACC = H
      ERRACC = 1.0d0
      REJECT=.FALSE.
      LAST=.FALSE.
      NSING=0
      IRTRN=1
      IF (AUTNMS) THEN
         HD1=0.0D0
         HD2=0.0D0
         HD3=0.0D0
         HD4=0.0D0
      END IF
C -------- PREPARE BAND-WIDTHS --------
      MBDIAG=MUMAS+1

C --- BASIC INTEGRATION STEP
      LAST = .FALSE.
      DO WHILE (.NOT.LAST)
      IF (.NOT. REJECT) THEN
      IF (NSTEP.GT.NMAX) CALL FAIL_EXIT(3,X,IDID,H,NMAX)
      IF ( 0.1D0*DABS(H) .LE. DABS(X)*UROUND )
     *       CALL FAIL_EXIT(2,X,IDID,H,NMAX)
      HOPT=H
      IF ((X+H*1.0001D0-XEND)*POSNEG.GE.0.D0) THEN
         H=XEND-X
         LAST=.TRUE.
      END IF
C *** *** *** *** *** *** ***
C  COMPUTATION OF THE JACOBIAN
C *** *** *** *** *** *** ***
      CALL FCN(N,X,Y,DY1)
      CALL JAC(N,X,Y,FJAC,LDJAC)
      NFCN=NFCN+1
      NJAC=NJAC+1

      IF (.NOT.AUTNMS) THEN
C --- COMPUTE NUMERICALLY THE DERIVATIVE WITH RESPECT TO X
            DELT=DSQRT(UROUND*DMAX1(1.D-5,DABS(X)))
            XDELT=X+DELT
            CALL FCN(N,XDELT,Y,AK1)
            DO J=1,N
               FX(J)=(AK1(J)-DY1(J))/DELT
            END DO
      END IF
      END IF

C *** *** *** *** *** *** ***
C  COMPUTE THE STAGES
C *** *** *** *** *** *** ***
      SINGULAR = .TRUE.
      DO WHILE (SINGULAR)
        FAC=1.D0/(H*GAMMA)
        CALL DECOMR(N,FJAC,LDJAC,FMAS,LDMAS,MLMAS,MUMAS,
     &            M1,M2,NM1,FAC,E,LDE,IP,IER,IJOB,IMPLCT,IP)
        SINGULAR = IER.NE.0
        IF (SINGULAR) THEN
          NSING=NSING+1
          IF (NSING.GE.5) CALL FAIL_EXIT(1,X,IDID,H,NMAX)
          H=H*0.5D0
          REJECT=.TRUE.
          LAST=.FALSE.
          ONE = .FALSE.
        END IF
      END DO

      NDEC=NDEC+1
C --- PREPARE FOR THE COMPUTATION OF THE 6 STAGES
      HC21=C21/H
      HC31=C31/H
      HC32=C32/H
      HC41=C41/H
      HC42=C42/H
      HC43=C43/H
      HC51=C51/H
      HC52=C52/H
      HC53=C53/H
      HC54=C54/H
      HC61=C61/H
      HC62=C62/H
      HC63=C63/H
      HC64=C64/H
      HC65=C65/H
      IF (.NOT.AUTNMS) THEN
         HD1=H*D1
         HD2=H*D2
         HD3=H*D3
         HD4=H*D4
      END IF
C --- THE STAGES
      CALL SLVROD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
     &    M1,M2,NM1,FAC,E,LDE,IP,DY1,AK1,FX,YNEW,HD1,IJOB,.FALSE.)
      DO  I=1,N
         YNEW(I)=Y(I)+A21*AK1(I)
      END DO
      CALL FCN(N,X+C2*H,YNEW,DY)
      DO I=1,N
         YNEW(I)=HC21*AK1(I)
      END DO
      CALL SLVROD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
     &    M1,M2,NM1,FAC,E,LDE,IP,DY,AK2,FX,YNEW,HD2,IJOB,.TRUE.)
      DO I=1,N
         YNEW(I)=Y(I)+A31*AK1(I)+A32*AK2(I)
      END DO
      CALL FCN(N,X+C3*H,YNEW,DY)
      DO I=1,N
         YNEW(I)=HC31*AK1(I)+HC32*AK2(I)
      END DO
      CALL SLVROD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
     &    M1,M2,NM1,FAC,E,LDE,IP,DY,AK3,FX,YNEW,HD3,IJOB,.TRUE.)
      DO I=1,N
         YNEW(I)=Y(I)+A41*AK1(I)+A42*AK2(I)+A43*AK3(I)
      END DO
      CALL FCN(N,X+C4*H,YNEW,DY)
      DO I=1,N
         YNEW(I)=HC41*AK1(I)+HC42*AK2(I)+HC43*AK3(I)
      END DO
      CALL SLVROD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
     &    M1,M2,NM1,FAC,E,LDE,IP,DY,AK4,FX,YNEW,HD4,IJOB,.TRUE.)
      DO I=1,N
         YNEW(I)=Y(I)+A51*AK1(I)+A52*AK2(I)+A53*AK3(I)+A54*AK4(I)
      END DO
      CALL FCN(N,X+H,YNEW,DY)
      DO I=1,N
         AK6(I)=HC52*AK2(I)+HC54*AK4(I)+HC51*AK1(I)+HC53*AK3(I)
      END DO
      CALL SLVROD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
     &    M1,M2,NM1,FAC,E,LDE,IP,DY,AK5,FX,AK6,0.D0,IJOB,.TRUE.)
C ------------ EMBEDDED SOLUTION ---------------
      DO I=1,N
         YNEW(I)=YNEW(I)+AK5(I)
      END DO
      CALL FCN(N,X+H,YNEW,DY)
      DO I=1,N
         AK5(I)=HC61*AK1(I)+HC62*AK2(I)+HC65*AK5(I)
     &              +HC64*AK4(I)+HC63*AK3(I)
      END DO
      CALL SLVROD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
     &    M1,M2,NM1,FAC,E,LDE,IP,DY,AK6,FX,AK5,0.D0,IJOB,.TRUE.)
C ------------ NEW SOLUTION ---------------
      DO  I=1,N
         YNEW(I)=YNEW(I)+AK6(I)
      END DO
      NSOL=NSOL+6
      NFCN=NFCN+5

C *** *** *** *** *** *** ***
C  ERROR ESTIMATION
C *** *** *** *** *** *** ***
      NSTEP=NSTEP+1
C ------------ COMPUTE ERROR ESTIMATION ----------------
      ERR=0.0D0
      DO I=1,N
         IF (ITOL.EQ.0) THEN
            SK=AbsTol(1)+RelTol(1)*DMAX1(DABS(Y(I)),DABS(YNEW(I)))
         ELSE
            SK=AbsTol(I)+RelTol(I)*DMAX1(DABS(Y(I)),DABS(YNEW(I)))
         END IF
         ERR=ERR+(AK6(I)/SK)**2
c2           ERR = DMAX1(ERR, AK6(I)/SK)
      END DO
      ERR=DSQRT(ERR/N)

C --- COMPUTATION OF HNEW
C --- WE REQUIRE .2<=HNEW/H<=6.
      FAC=DMAX1(FAC2,DMIN1(FAC1,(ERR)**0.25D0/SAFE))
      HNEW=DMAX1(H/FAC, STEPMIN)

C *** *** *** *** *** *** ***
C  IS THE ERROR SMALL ENOUGH ?
C *** *** *** *** *** *** ***

      IF ( (ERR.LE.1.D0).or.(H.LE.STEPMIN) ) THEN
C --- STEP IS ACCEPTED
         NACCPT=NACCPT+1
         IF (PRED) THEN
C       --- PREDICTIVE CONTROLLER OF GUSTAFSSON
            IF (NACCPT.GT.1) THEN
               FACGUS=(HACC/H)*(ERR**2/ERRACC)**0.25D0/SAFE
               FACGUS=DMAX1(FAC2,DMIN1(FAC1,FACGUS))
               FAC=DMAX1(FAC,FACGUS)
               HNEW=DMAX1(H/FAC, STEPMIN)
            END IF
            HACC=H
            ERRACC=DMAX1(1.0D-2,ERR)
         END IF
         DO I=1,N
            Y(I)=YNEW(I)
         END DO
         XOLD=X
         X=X+H
         IF (DABS(HNEW).GT.HMAXN) HNEW=POSNEG*HMAXN
         IF (REJECT) HNEW=POSNEG*DMIN1(DABS(HNEW),DABS(H))
         REJECT=.FALSE.
         H=HNEW
      ELSE
C --- STEP IS REJECTED
         REJECT=.TRUE.
         LAST=.FALSE.
         H=HNEW
         IF (NACCPT.GE.1) NREJCT=NREJCT+1
      END IF
      END DO
      RETURN
      END
C
      SUBROUTINE FAIL_EXIT(NERR,X,IDID,H,NMAX)
      INTEGER NERR, NMAX
      KPP_REAL X, H
      GO TO (1,2,3,4) NERR
 1    CONTINUE
      WRITE(6,979)X
      WRITE(6,*) ' MATRIX IS REPEATEDLY SINGULAR, IER=',IER
      IDID=-4
      STOP
 2    CONTINUE
      WRITE(6,979)X
      WRITE(6,*) ' STEP SIZE TOO SMALL, H=',H
      IDID=-3
      STOP
 3    CONTINUE
      WRITE(6,979)X
      WRITE(6,*) ' MORE THAN NMAX =',NMAX,'STEPS ARE NEEDED'
      IDID=-2
      STOP
C --- EXIT CAUSED BY solout
 4    CONTINUE
      WRITE(6,979)X
 979  FORMAT(' EXIT OF RODAS AT X=',E18.4)
      IDID=2
      RETURN
      END

      SUBROUTINE ROCOE(METH,A21,A31,A32,A41,A42,A43,A51,A52,A53,A54,
     &  C21,C31,C32,C41,C42,C43,C51,C52,C53,C54,C61,
     &  C62,C63,C64,C65,GAMMA,C2,C3,C4,D1,D2,D3,D4,
     &  D21,D22,D23,D24,D25,D31,D32,D33,D34,D35)
      IMPLICIT KPP_REAL (A-H,O-Z)

      if (METH.ne.1) print *, 'WRONG CHOICE OF METHOD'
        C2=0.386D0
        C3=0.21D0
        C4=0.63D0
        BET2P=0.0317D0
        BET3P=0.0635D0
        BET4P=0.3438D0
       D1= 0.2500000000000000D+00
       D2=-0.1043000000000000D+00
       D3= 0.1035000000000000D+00
       D4=-0.3620000000000023D-01
       A21= 0.1544000000000000D+01
       A31= 0.9466785280815826D+00
       A32= 0.2557011698983284D+00
       A41= 0.3314825187068521D+01
       A42= 0.2896124015972201D+01
       A43= 0.9986419139977817D+00
       A51= 0.1221224509226641D+01
       A52= 0.6019134481288629D+01
       A53= 0.1253708332932087D+02
       A54=-0.6878860361058950D+00
       C21=-0.5668800000000000D+01
       C31=-0.2430093356833875D+01
       C32=-0.2063599157091915D+00
       C41=-0.1073529058151375D+00
       C42=-0.9594562251023355D+01
       C43=-0.2047028614809616D+02
       C51= 0.7496443313967647D+01
       C52=-0.1024680431464352D+02
       C53=-0.3399990352819905D+02
       C54= 0.1170890893206160D+02
       C61= 0.8083246795921522D+01
       C62=-0.7981132988064893D+01
       C63=-0.3152159432874371D+02
       C64= 0.1631930543123136D+02
       C65=-0.6058818238834054D+01
       GAMMA= 0.2500000000000000D+00

       D21= 0.1012623508344586D+02
       D22=-0.7487995877610167D+01
       D23=-0.3480091861555747D+02
       D24=-0.7992771707568823D+01
       D25= 0.1025137723295662D+01
       D31=-0.6762803392801253D+00
       D32= 0.6087714651680015D+01
       D33= 0.1643084320892478D+02
       D34= 0.2476722511418386D+02
       D35=-0.6594389125716872D+01
      RETURN
      END
C

C ******************************************
C     VERSION OF SEPTEMBER 18, 1995
C ******************************************
C
      SUBROUTINE DECOMR(N,FJAC,LDJAC,FMAS,LDMAS,MLMAS,MUMAS,
     &            M1,M2,NM1,FAC1,E1,LDE1,IP1,IER,IJOB,CALHES,IPHES)
      IMPLICIT KPP_REAL (A-H,O-Z)
      INCLUDE 'KPP_ROOT_params.h'
      INCLUDE 'KPP_ROOT_sparse.h'
      DIMENSION FJAC(LU_NONZERO),FMAS(LDMAS,NM1),E1(LU_NONZERO),
     &          IP1(NM1),IPHES(N)
      LOGICAL CALHES
      COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
C

C ---  B=IDENTITY, JACOBIAN A FULL MATRIX
      DO J=1,LU_NONZERO
         E1(J) = -FJAC(J)
      END DO
      DO J=1,N
         E1(LU_DIAG(J)) = E1(LU_DIAG(J)) + FAC1
      END DO
      CALL KppDecomp (E1,IER)
      RETURN
      END
C
C     END OF SUBROUTINE DECOMR
C
C ***********************************************************
C
C
C
C
      SUBROUTINE SLVROD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
     &          M1,M2,NM1,FAC1,E,LDE,IP,DY,AK,FX,YNEW,HD,IJOB,STAGE1)
      IMPLICIT KPP_REAL (A-H,O-Z)
      INCLUDE 'KPP_ROOT_params.h'
      INCLUDE 'KPP_ROOT_sparse.h'
      DIMENSION FJAC(LU_NONZERO),FMAS(LDMAS,NM1),E(LU_NONZERO),
     &          IP(NM1),DY(N),AK(N),FX(N),YNEW(N)
      LOGICAL STAGE1
      COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
C
      IF (HD.EQ.0.D0) THEN
         DO  I=1,N
           AK(I)=DY(I)
         END DO
      ELSE
         DO I=1,N
            AK(I)=DY(I)+HD*FX(I)
         END DO
      END IF

C ---  B=IDENTITY, JACOBIAN A FULL MATRIX
      IF (STAGE1) THEN
         DO I=1,N
            AK(I)=AK(I)+YNEW(I)
         END DO
      END IF
      CALL KppSolve (E,AK)
      RETURN
      END
C
C     END OF SUBROUTINE SLVROD
C
C
C ***********************************************************

      SUBROUTINE FUNC_CHEM(N, T, Y, P)
      INCLUDE 'KPP_ROOT_params.h'
      INCLUDE 'KPP_ROOT_global.h'
      INTEGER N
      KPP_REAL   T, Told
      KPP_REAL   Y(NVAR), P(NVAR)
      Told = TIME
      TIME = T
      CALL Update_SUN()
      CALL Update_RCONST()
      CALL Fun( Y,  FIX, RCONST, P )
      TIME = Told
      RETURN
      END

      SUBROUTINE JAC_CHEM(N, T, Y, J)
      INCLUDE 'KPP_ROOT_params.h'
      INCLUDE 'KPP_ROOT_global.h'
      INTEGER N
      KPP_REAL   Told, T
      KPP_REAL   Y(NVAR), J(LU_NONZERO)
      Told = TIME
      TIME = T
      CALL Update_SUN()
      CALL Update_RCONST()
      CALL Jac_SP( Y,  FIX, RCONST, J )
      TIME = Told
      RETURN
      END

max_stars_repo_name: matzegoebel/WRF-fluxavg

max_issues_repo_name: teb-model/wrf-teb

max_forks_repo_name: teb-model/wrf-teb


module fapi_read_tests

    use odc
    use odc_config
    use, intrinsic :: iso_fortran_env
    implicit none

contains

    subroutine check_call(err, msg, success)
        integer, intent(in) :: err
        character(*), intent(in) :: msg
        logical, intent(inout) :: success

        if (err /= ODC_SUCCESS) then
            write(error_unit, *) 'Failed API call: ', msg
            write(error_unit, *) 'Error: ', odc_error_string(err)
            success = .false.
        end if
    end subroutine

    function test_count_lines() result(success)

        ! Test that we obtain the expected version number

        type(odc_reader) :: reader
        type(odc_frame) :: frame
        integer(8) :: frame_count, row_count, tmp_8
        integer :: err, tmp_4
        logical :: success

        success = .true.
        call check_call(reader%open_path("../2000010106.odb"), "open ODB", success)
        call check_call(frame%initialise(reader), "initialise frame", success)

        frame_count = 0
        row_count = 0

        err = frame%next()
        do while (err == ODC_SUCCESS)

            call check_call(frame%row_count(tmp_8), "row count", success)

            frame_count = frame_count + 1
            row_count = row_count + tmp_8

            call check_call(frame%column_count(tmp_4), "column count", success)
            if (tmp_4 /= 51) then
                write(error_unit, *) 'Unexpected column count: ', tmp_4, ' /= 51'
                success = .false.
            endif

            err = frame%next()
        end do

        if (err /= ODC_ITERATION_COMPLETE) call check_call(err, "next frame", success)

        if (frame_count /= 333) then
            write(error_unit, *) 'Unexpected frame count: ', frame_count, ' /= 333'
            success = .false.
        endif

        if (row_count /= 3321753) then
            write(error_unit, *) 'Unexpected row count: ', row_count, ' /= 3321753'
            success = .false.
        endif

        call check_call(reader%close(), "close reader", success)

    end function

    function test_column_details() result(success)

        type(odc_reader) :: reader
        type(odc_frame) :: frame
        character(:), allocatable :: column_name, field_name
        integer :: ncols, col, column_type, field, field_size, expected_offset, field_offset
        integer :: element_size, element_size_doubles, bitfield_count
        logical :: success

        character(23), parameter :: example_column_names(*) = [ character(23) :: &
            "expver@desc", "andate@desc", "antime@desc", "seqno@hdr", "obstype@hdr", &
            "obschar@hdr", "subtype@hdr", "date@hdr", "time@hdr", "rdbflag@hdr", &
            "status@hdr", "event1@hdr", "blacklist@hdr", "sortbox@hdr", "sitedep@hdr", &
            "statid@hdr", "ident@hdr", "lat@hdr", "lon@hdr", "stalt@hdr", &
            "modoro@hdr", "trlat@hdr", "trlon@hdr", "instspec@hdr", "event2@hdr", &
            "anemoht@hdr", "baroht@hdr", "sensor@hdr", "numlev@hdr", "varno_presence@hdr", &
            "varno@body", "vertco_type@body", "rdbflag@body", "anflag@body", "status@body", &
            "event1@body", "blacklist@body", "entryno@body", "press@body", "press_rl@body", &
            "obsvalue@body", "aux1@body", "event2@body", "ppcode@body", "level@body", &
            "biascorr@body", "final_obs_error@errstat", "obs_error@errstat", "repres_error@errstat", &
            "pers_error@errstat", "fg_error@errstat"]

        integer, parameter :: example_column_types(*) = [ &
            ODC_STRING, ODC_INTEGER, ODC_INTEGER, ODC_INTEGER, ODC_INTEGER, ODC_BITFIELD , &
            ODC_INTEGER, ODC_INTEGER, ODC_INTEGER, ODC_BITFIELD , ODC_BITFIELD , ODC_BITFIELD , &
            ODC_BITFIELD , ODC_INTEGER, ODC_INTEGER, ODC_STRING, ODC_INTEGER, ODC_REAL, &
            ODC_REAL, ODC_REAL, ODC_REAL, ODC_REAL, ODC_REAL, ODC_INTEGER, &
            ODC_INTEGER, ODC_REAL, ODC_REAL, ODC_INTEGER, ODC_INTEGER, ODC_BITFIELD , &
            ODC_INTEGER, ODC_INTEGER, ODC_BITFIELD , ODC_BITFIELD , ODC_BITFIELD , ODC_BITFIELD , &
            ODC_BITFIELD , ODC_INTEGER, ODC_REAL, ODC_REAL, ODC_REAL, ODC_REAL, &
            ODC_INTEGER, ODC_INTEGER, ODC_BITFIELD , ODC_REAL, ODC_REAL, ODC_REAL, &
            ODC_REAL, ODC_REAL, ODC_REAL]

        character(14), parameter :: column_10_bitfield_names(*) = [ character(14) :: &
            "lat_humon", "lat_qcsub", "lat_override", "lat_flag", "lat_hqc_flag", "lon_humon", "lon_qcsub", &
            "lon_override", "lon_flag", "lon_hqc_flag", "date_humon", "date_qcsub", "date_override", &
            "date_flag", "date_hqc_flag", "time_humon", "time_qcsub", "time_override", "time_flag", &
            "time_hqc_flag", "stalt_humon", "stalt_qcsub", "stalt_override", "stalt_flag", "stalt_hqc_flag" &
        ]

        integer, parameter :: column_10_bitfield_sizes(*) = [ &
            1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1]

        success = .true.
        call check_call(reader%open_path("../2000010106.odb"), "open reader", success)
        call check_call(frame%initialise(reader), "initialise frame", success)

        call check_call(frame%next(), "get the first frame", success)

        call check_call(frame%column_count(ncols), "column count", success)
        if (ncols /= 51) then
            write(error_unit, *) 'Expected 51 columns'
            success = .false.
        endif

        ! n.b. -- 1-based indexing!
        do col = 1, ncols

            call check_call(frame%column_attributes(col, &
                                                    name=column_name, &
                                                    type=column_type, &
                                                    element_size=element_size, &
                                                    element_size_doubles=element_size_doubles, &
                                                    bitfield_count=bitfield_count), "column attrs", success)

            if (column_name /= trim(example_column_names(col))) then
                write(error_unit,'(3a,i2,3a)') 'Unexpected column name ', column_name, &
                           ' for column ', col, ' (expected ', trim(example_column_names(col)), ')'
                success = .false.
            end if

            if (column_type /= example_column_types(col)) then
                write(error_unit, '(a,i1,a,i2,a,i1,a)') 'Unexpected column type ', column_type, &
                        ' for column ', col, ' (expected ', example_column_types(col), ')'
                success = .false.
            end if

            if (element_size /= 8) then
                write(error_unit, '(a,i1,a,i2,a)') 'Unexpected column data size ', element_size, &
                        ' for column ', col, ' (expected 8)'
                success = .false.
            end if

            if (element_size_doubles /= 1) then
                write(error_unit, '(a,i1,a,i2,a)') 'Unexpected column doubles data size ', &
                        element_size_doubles, ' for column ', col, ' (expected 1)'
                success = .false.
            end if

            if (column_type == ODC_BITFIELD) then
                if (bitfield_count <= 0) then
                    write(error_unit, *) "Bitfields expected for bitfield column"
                    success = .false.
                end if
            else
                if (bitfield_count /= 0) then
                    write(error_unit, *) "Unexpected bitfields for non-bitfield column"
                    success = .false.
                end if
            end if

        end do

        ! Test bitfields for column 10

        call check_call(frame%column_attributes(10, bitfield_count=bitfield_count), "bitfield count", success)
        if (bitfield_count /= 25) then
            write(error_unit, *) "Expected 25 bitfield fields for column 10. Got ", bitfield_count
            success = .false.
        end if

        expected_offset = 0
        do field = 1, 25

            ! Look at column 10
            call check_call(frame%bitfield_attributes(10, field, &
                                                      name=field_name, &
                                                      offset=field_offset, &
                                                      size=field_size), "bitfield attrs", success)

            if (field_name /= trim(column_10_bitfield_names(field))) then
                write(error_unit, '(3a,i2,3a)') 'Unexpected field name ', field_name, ' for field ', &
                        field, ' (expected ', trim(column_10_bitfield_names(field)), ')'
                success = .false.
            end if

            if (field_size /= column_10_bitfield_sizes(field)) then
                write(error_unit, '(a,i2,a,i2,a,i2,a)') 'Unexpected field size ', field_size, &
                        ' for field ', field, ' (expected ', column_10_bitfield_sizes(field), ')'
                success = .false.
            end if

            if (field_offset /= expected_offset) then
                write(error_unit, '(a,i2,a,i2,a,i2,a)') 'Unexpected field offset ', field_offset, &
                        ' for field ', field, ' (expected ', expected_offset, ')'
                success = .false.
            end if

            expected_offset = expected_offset + field_size
        end do

        call check_call(frame%free(), "free frame", success)
        call check_call(reader%close(), "close reader", success)
    end function

    function check_frame_2_values(array_data) result(success)

        real(8) :: array_data(:,:)
        logical :: success

        integer :: row, i
        integer, parameter :: expected_seqno(*) = [6106691, 6002945, 6003233, 6105819]
        integer, parameter :: expected_obschar(*) = [537918674, 135265490, 135265490, 537918674]
        integer(8) :: missing_integer
        real(8) :: missing_double

        success = .true.

        call check_call(odc_missing_integer(missing_integer), "missing integer", success)
        call check_call(odc_missing_double(missing_double), "missing double", success)

        do i = 1, 4
            row = 1 + ((i-1) * 765)

            ! Expver
            if (trim(transfer(array_data(row, 1), "        ")) /= "0018") then
                write(error_unit, *) 'unexpected expver in row ', row, ' (expected 0018, got ', &
                                     transfer(array_data(row, 1), "        ") ,')'
                success = .false.
            end if

            ! Test seqno (INTEGER)
            if (int(array_data(row, 4)) /= expected_seqno(i)) then
                write(error_unit, *) 'Unexpected seqno value. row=', row, ", expected=", &
                                     expected_seqno(i), ", got=", int(array_data(row, 4))
                success = .false.
            end if

            ! obschar (BITFIELD)
            if (int(array_data(row, 6)) /= expected_obschar(i)) then
                write(error_unit, *) 'Unexpected obschar value. row=', row, ", expected=", &
                                     expected_obschar(i), ", got=", int(array_data(row, 6))
                success = .false.
            end if

            ! Sortbox (INTEGER, missing)
            if (int(array_data(row, 14)) /= missing_integer) then
                write(error_unit, *) 'Expected value with set missing value. Got ', int(array_data(row, 14)), ', &
                                     &expected ', missing_integer
                success = .false.
            end if

            ! repres_error (REAL, missing)
            if (array_data(row, 49) /= missing_double) then
                write(error_unit, *) 'Expected value with set missing value. Got ', array_data(row, 49), ', &
                                     &expected ', missing_double
                success = .false.
            end if
        end do

    end function

    function test_decode_columns_allocate() result(success)

        type(odc_reader) :: reader
        type(odc_frame) :: frame
        type(odc_decoder) :: decoder
        integer(8) :: nrows, nrows2
        integer :: ncols
        logical :: success, column_major
        real(8), pointer :: array_data(:,:)

        success =.true.
        call check_call(reader%open_path("../2000010106.odb"), "open reader", success)
        call check_call(frame%initialise(reader), "initialise frame", success)

        ! Read the second frame, because why not.
        call check_call(frame%next(), "get first frame", success)
        call check_call(frame%next(), "get second frame", success)

        call check_call(decoder%initialise(), "initialise decoder", success)
        call check_call(decoder%defaults_from_frame(frame), "decoder from frame", success)
        call check_call(decoder%decode(frame, nrows), "do decode", success)

        if (nrows /= 10000) then
            write(error_unit, *) 'Unexpected number of rows decoded'
            success = .false.
        end if

        call check_call(decoder%row_count(nrows2), "decoder row count", success)
        if (nrows2 /= 10000) then
            write(error_unit, *) 'Got row count ', nrows, ' not 10000'
            success = .false.
        end if

        call check_call(decoder%column_count(ncols), "decoder column count", success)
        if (ncols /= 51) then
            write(error_unit, *) 'Got column count ', ncols, ' not 51'
            success = .false.
        end if

        call check_call(decoder%data(array_data, column_major), "get decoded data", success)

        if (any(shape(array_data) /= [10000, 51])) then
            write(error_unit, *) 'Unexpected data dimensions'
            success = .false.
        end if

        if (.not. column_major) then
            write(error_unit, *) 'Expected column major by default'
            success = .false.
        end if

        success = success .and. check_frame_2_values(array_data)

        call check_call(decoder%free(), "free decoder", success)
        call check_call(reader%close(), "free reader", success)

    end function

    function test_decode_array_reuse() result(success)

        use, intrinsic :: iso_c_binding

        type(odc_reader) :: reader
        type(odc_frame) :: frame
        type(odc_decoder) :: decoder
        integer(8) :: rows_decoded, nrows
        integer :: ncols
        logical :: success
        real(8), target :: array_data(11000, 51)

        success = .true.
        call check_call(reader%open_path("../2000010106.odb"), "open reader", success)
        call check_call(frame%initialise(reader), "initialise frame", success)

        call check_call(frame%next(), "get first frame", success)

        call check_call(decoder%initialise(), "initialise decoder", success)
        call check_call(decoder%defaults_from_frame(frame), "decoder frame defaults", success)
        call check_call(decoder%set_data(array_data), "set array data", success)
        call check_call(decoder%decode(frame, rows_decoded), "decode first frame", success)

        if (rows_decoded /= 10000) then
            write(error_unit, *) 'Unexpected number of rows decoded'
            success = .false.
        end if

        call check_call(frame%next(), "get second frame", success)

        call check_call(decoder%decode(frame, rows_decoded), "decode second frame", success)

        if (rows_decoded /= 10000) then
            write(error_unit, *) 'Unexpected number of rows decoded'
            success = .false.
        end if

        call check_call(decoder%row_count(nrows), "decoder row count", success)
        if (nrows /= 11000) then
            write(error_unit, *) 'Got row count ', nrows, ' not 11000'
            write(error_unit, *) 'Row count should be related to the size of the decode target, not the decode data'
            success = .false.
        end if

        call check_call(decoder%column_count(ncols), "decoder column count", success)
        if (ncols /= 51) then
            write(error_unit, *) 'Got column count ', ncols, ' not 51'
            success = .false.
        end if

        success = success .and. check_frame_2_values(array_data)

        call check_call(decoder%free(), "free decoder", success)
        call check_call(reader%close(), "close reader", success)

    end function

    function test_decode_aggregate() result(success)

        type(odc_reader) :: reader
        type(odc_frame) :: frame
        type(odc_decoder) :: decoder
        integer(8) :: rows_decoded, nrows
        integer :: ncols
        logical :: success
        real(8), pointer :: array_data(:,:)

        success = .true.
        call check_call(reader%open_path("../2000010106.odb"), "open reader", success)
        call check_call(frame%initialise(reader), "initialise frame", success)

        call check_call(frame%next(maximum_rows=99999_8), "get first (aggregate) frame", success)

        call check_call(decoder%initialise(), "initialise decoder", success)
        call check_call(decoder%defaults_from_frame(frame), "decoder frame defaults", success)
        call check_call(decoder%decode(frame, rows_decoded, nthreads=4), "decode threaded", success)

        if (rows_decoded /= 90000) then
            write(error_unit, *) 'Unexpected number of rows decoded'
            success = .false.
        end if

        call check_call(decoder%row_count(nrows), "decoder row count", success)
        if (nrows /= 90000) then
            write(error_unit, *) 'Got row count ', nrows, ' not 90000'
            success = .false.
        end if

        call check_call(decoder%column_count(ncols), "decoder column count", success)
        if (ncols /= 51) then
            write(error_unit, *) 'Got column count ', ncols, ' not 51'
            success = .false.
        end if

        call check_call(decoder%data(array_data), "get array data", success)

        if (any(shape(array_data) /= [90000, 51])) then
            write(error_unit, *) 'Unexpected data dimensions'
            success = .false.
        end if

        call check_call(decoder%free(), "free decoder", success)
        call check_call(reader%close(), "close reader", success)

    end function

end module

program fapi_general

    use fapi_read_tests
    implicit none

    logical :: success

    success = .true.
    call check_call(odc_initialise_api(), "initialise api", success)

    success = test_count_lines() .and. success
    success = test_column_details() .and. success
    success = test_decode_columns_allocate() .and. success
    success = test_decode_array_reuse() .and. success
    success = test_decode_aggregate() .and. success

    if (.not. success) stop -1

end program


text:
!
! Copyright (c) 2009           HLRS. All rights reserved.
! $COPYRIGHT$
!
! Additional copyrights may follow
!
! $HEADER$
!

module auxdata3df
   use adcl
contains

!***************************************************************************
subroutine set_data_3D_cont ( data, rank, dims, hwidth, cart_comm)
! numbers the data array continously over all processes

    implicit none
    include 'ADCL.inc'

    integer, intent(in) :: rank, hwidth, cart_comm
    integer, dimension(3), intent(in) :: dims
    double precision, intent(inout) :: data(dims(1),dims(2),dims(3))

    integer :: i, j, k, ierr
    integer, dimension(3) :: coords, cart_dims, dims_wo_halos
    logical, dimension(3) :: period

    call MPI_Cart_get(cart_comm, 3, cart_dims, period, coords, ierr)
    ! cube: cart_dims with index coords
    ! inside cube: dims with index (i,j,k)

    do i=1, dims(1)
        do j=1, dims(2)
            do k=1, hwidth
                data(i,j,k) = -1
            end do
            do k=dims(3)-hwidth+1, dims(3)
                data(i,j,k) = -1
           end do
       end do
   end do

    do i=1, dims(1)
        do k=1, dims(3)
            do j=1, hwidth
                data(i,j,k) = -1
           end do
            do j=dims(2)-hwidth+1, dims(2)
                data(i,j,k) = -1
           end do
       end do
   end do

    do j=1, dims(2)
        do k=1, dims(3)
            do i=1, hwidth
                data(i,j,k) = -1
           end do
            do i=dims(1)-hwidth+1, dims(1)
                data(i,j,k) = -1
           end do
       end do
   end do

    !be aware of the change in i and j in the equation compare to the c code. (i-1) and (j-1)
    ! without hwidth:
    !    in x direction offset dims(1)*coords(1) + i
    !    in y direction offset dims(1)*cart_dims(1) * ( dims(2)*coords(2) + j-1 )
    !    in z direction offset dims(1)*cart_dims(1) * dims(2)*cart_dims(2) * ( dims(3)*coords(3) + z-1 )
    dims_wo_halos = dims - 2*hwidth

    do  i=hwidth+1, dims(1)-hwidth
        do j=hwidth+1, dims(2)-hwidth
            do k=hwidth+1, dims(3)-hwidth
                !print *, "rank =", rank, ", coords = ( ", coords, " )"
                data(i,j,k) = i-hwidth + dims_wo_halos(1)*coords(1) +                                                     &
                              dims_wo_halos(1)*cart_dims(1) * ( ( dims_wo_halos(2)*coords(2) + j-hwidth-1 ) +             &
                                             dims_wo_halos(2)*cart_dims(2) * ( dims_wo_halos(3)*coords(3) + k-hwidth-1 ) )
           end do
       end do
    end do

    return
end subroutine set_data_3D_cont

!****************************************************************************
function check_data_3D_cont ( data, rank, dims, hwidth, neighbors, cart_comm ) result(isok)

   implicit none
   include 'ADCL.inc'

   integer, intent(in) :: rank, hwidth, cart_comm
   integer, dimension(3), intent(in) :: dims
   integer, dimension(6), intent(in) :: neighbors
   double precision, intent(in) :: data(dims(1),dims(2), dims(3))
   logical :: isok

   integer :: i, j, k, ierr, lres=1, gres, prod, ierror
   integer, dimension(3) :: coords, n_coords, c_coords, cart_size
   double precision should_be

   integer :: x_direction, y_direction, z_direction

   ! check for each of the 27 possible locations
   do z_direction = -1, 1
        do y_direction = -1, 1
            do x_direction = -1, 1
                prod = x_direction * y_direction * z_direction
#ifdef INCCORNER
                if ( prod .ne. 0 ) then
                    ! corner
                    lres = check_region_3D (x_direction, y_direction, z_direction, data, rank, cart_comm, &
                           dims, hwidth, neighbors)
                endif
#endif
                if ( prod == 0) then
                    ! edge, face or inside
                    lres = check_region_3D (x_direction, y_direction, z_direction, data, rank, cart_comm, &
                           dims, hwidth, neighbors )
                end if
                call MPI_Allreduce ( lres, gres, 1, MPI_INTEGER, MPI_MIN, cart_comm, ierror )
                if ( gres .ne. 1 ) then
                   isok = .false.
                   return
                endif
            end do
        end do
    end do

    !if ( gres .eq. 1 ) then
    !    if ( rank == 0 )  then
    !        write(*,'(1x,a,i0,a)') "3-D C testsuite: hwidth = ", hwidth, ", nc = 0 passed"
    !    end if
    !else
    !    if ( rank == 0 ) then
    !        write(*,'(1x,a,i0,a)') "3-D C testsuite: hwidth = ", hwidth, ", nc = 0 failed"
    !    end if
    !    call dump_vector_3D_mpi_dp ( data, dims, cart_comm )
    !endif

   isok = .true.
   return

end function check_data_3D_cont

!****************************************************************************

function check_region_3D (x_direction, y_direction, z_direction, data, rank, cart_comm, dims, hwidth, neighbors) result(lres)
! checks all data in one direction i.e. on one face, edge (and corner if defined)
    integer, intent(in) :: x_direction, y_direction, z_direction, rank, cart_comm
    integer, intent(in) :: dims(3)                          ! size of one cube
    integer, intent(in) :: hwidth, neighbors(6)
    double precision, intent(in) :: data (dims(1), dims(2), dims(3) )
    integer :: lres

    integer :: i, j, k, ierr
    double precision :: should_be
    integer :: cart_dims(3), dims_wo_halos(3)
    integer :: coords(3), n_coords(3), c_coords(3) ! coords, coords of neighbor and corrected coords of MPI process
    logical, dimension(3) :: period ! not really used
    integer :: compensate(3)  ! what do I have to add / substract to my coordinate on the
                              ! neighboring process to compare the values in the halo cell and
                              ! in the domain of the neighboring process
                              ! be aware of shift (hwidth+1,hwidth+1,hwidth+1) is no. 1 !
    integer :: loopstart(3), loopend(3)  ! defines part of data to check
    logical :: neighbor_cond(3)          ! is there a neighbor in x,y,z-direction?

    lres = 1

    call MPI_Cart_get (cart_comm, 3, cart_dims, period, coords, ierr)
    neighbor_cond = .false.

    select case (x_direction)
       case (0)
            loopstart(1)     = hwidth+1
            loopend(1)       = dims(1)-hwidth
            compensate(1)    = -hwidth
            c_coords(1)      = coords(1)
            neighbor_cond(1) = .true.
       case (-1)
            loopstart(1) = 1
            loopend(1) = hwidth
            compensate(1) = dims(1) - 3*hwidth ! 2*hwidth for positioning and hwidth for shift of numbering
            if (neighbors(1) .ne. MPI_PROC_NULL ) then
                call MPI_Cart_coords (cart_comm, neighbors(1), 3, n_coords, ierror )
                c_coords(1) = n_coords(1)
                neighbor_cond(1) = .true.
            endif
       case (1)
            loopstart(1)     = dims(1) - hwidth + 1
            loopend(1)       = dims(1)
            compensate(1)    = - dims(1) + hwidth
            if (neighbors(2) .ne. MPI_PROC_NULL ) then
                call MPI_Cart_coords (cart_comm, neighbors(2), 3, n_coords, ierror )
                c_coords(1) = n_coords(1)
                neighbor_cond(1) = .true.
            endif
    end select

    select case (y_direction)
        case (0)
            loopstart(2) = hwidth+1
            loopend(2) = dims(2) - hwidth
            compensate(2) = -hwidth
            c_coords(2) = coords(2)
            neighbor_cond(2) = .true.
        case (-1)
            loopstart(2) = 1
            loopend(2)   = hwidth
            compensate(2) = dims(2) - 3*hwidth
            if (neighbors(3) .ne. MPI_PROC_NULL ) then
                call MPI_Cart_coords (cart_comm, neighbors(3), 3, n_coords, ierror )
                c_coords(2) = n_coords(2)
                neighbor_cond(2) = .true.
            endif
        case (1)
            loopstart(2)  = dims(2) - hwidth + 1
            loopend(2)    = dims(2)
            compensate(2) = - dims(2) + hwidth
            if (neighbors(4) .ne. MPI_PROC_NULL ) then
                call MPI_Cart_coords (cart_comm, neighbors(4), 3, n_coords, ierror )
                c_coords(2) = n_coords(2)
                neighbor_cond(2) = .true.
            endif
    end select

    select case (z_direction)
        case (0)
            loopstart(3)  = hwidth + 1
            loopend(3)    = dims(3) - hwidth
            compensate(3) = -hwidth
            c_coords(3)   = coords(3)
            neighbor_cond(3) = .true.
        case (-1)
            loopstart(3)  = 1
            loopend(3)    = hwidth
            compensate(3) = dims(3) - 3*hwidth
            if (neighbors(5) .ne. MPI_PROC_NULL ) then
                call MPI_Cart_coords (cart_comm, neighbors(5), 3, n_coords, ierror )
                c_coords(3) = n_coords(3)
                neighbor_cond(3) = .true.
            endif
        case (1)
            loopstart(3)  = dims(3) - hwidth + 1
            loopend(3)    = dims(3)
            compensate(3) = - dims(3) + hwidth
            if (neighbors(6) .ne. MPI_PROC_NULL ) then
                call MPI_Cart_coords (cart_comm, neighbors(6), 3, n_coords, ierror )
                c_coords(3) = n_coords(3)
                neighbor_cond(3) = .true.
            endif
    end select

   dims_wo_halos = dims - 2*hwidth

    do k = loopstart(3), loopend(3)
       do  j= loopstart(2), loopend(2)
           do  i = loopstart(1), loopend(1)
               if (neighbor_cond(1) .and. neighbor_cond(2) .and. neighbor_cond(3)) then
                  should_be = i + compensate(1) + dims_wo_halos(1)*c_coords(1) +               &
                       dims_wo_halos(1)*cart_dims(1) *  &
                                ( ( dims_wo_halos(2)*c_coords(2) + j-1 + compensate(2) ) +                      &
                                    dims_wo_halos(2)*cart_dims(2) * ( dims_wo_halos(3)*c_coords(3) + k-1 + compensate(3) ) )
               else
                   should_be = -1
               end if
               if ( data(i,j,k) .ne. should_be ) then
                   lres = 0
                   write(*,'(i4,a,3I4,a,f12.4,a,f12.4,a,3i3)') rank, ": data(",i,j,k,") = ", data(i,j,k), &
                      ", should_be, ", should_be, ", direction =", x_direction, y_direction, z_direction
               end if
            end do
        end do
    end do
    return

end function check_region_3D

!****************************************************************************
      subroutine DUMP_VECTOR ( arr, rank, dims )

        implicit none
        include 'ADCL.inc'
        character(100) ::  fmt

        integer rank, dims(3)
        DATATYPE arr(dims(1), dims(2), dims(3))
        integer i, j, k

        write (fmt,'(a,i0,a)') "(i3,a1,2i5,a1,", dims(1), "f12.5)"
        do k = 1, dims(3)
           do j = 1, dims(2)
              write (*,fmt) rank, ":", k, j, ":", (arr(i,j,k), i=1,dims(1))
           end do
        end do

        return
      end subroutine DUMP_VECTOR

!***************************************************************************
      subroutine DUMP_VECTOR_MPI ( arr, dims, comm )

        implicit none
        include 'ADCL.inc'

        integer dims(3), comm
        DATATYPE arr(dims(1), dims(2), dims(3))
        integer i, j, k, iproc, rank, size, ierror
        character(100) ::  fmt

        call MPI_Comm_rank ( comm, rank, ierror )
        call MPI_Comm_size ( comm, size, ierror )

        write (fmt,'(a,i0,a)') "(i3,a1,2i5,a1,", dims(1), "f12.5)"

        do iproc = 0, size-1
           if ( iproc .eq. rank) then
              do k = 1, dims(3)
                 do j = 1, dims(2)
                    write (6,fmt) rank, ":", k, j, ":", (arr(i,j,k), i=1,dims(1))
                 end do
              end do
           end if
           call MPI_Barrier ( comm, ierror )
        end do
        return
      end subroutine DUMP_VECTOR_MPI

end module auxdata3df

max_stars_repo_name: edgargabriel/ADCL

max_issues_repo_name: PSTL-UH/ADCL

max_forks_repo_name: PSTL-UH/ADCL