text:
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                             c
c                  Calculating the Peak Likelihood of the ZDCF                c
c                                                                             c
c        A special version for processing the ZDCF output file *.dcf          c
c                                                                             c
c                                Version 3.0                                  c
c                                                                             c
c    (Bayesian prior assumption of uniform distribution in z-space)           c
c                                                                             c
c                                                                             c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                             c
c SAMPLE RUN (user input marked by "<-USER". See explanations below)          c
c ----------                                                                  c
c                                                                             c
c PLIKE V3.0 begins.                                                          c
c Enter dcf file name:                                                        c
c test.dcf       <-USER                                                       c
c Enter lower bound on peak location:                                         c
c -200           <-USER                                                       c
c Enter upper bound on peak location:                                         c
c 200            <-USER                                                       c
c Calculating ML in range  -200.  to  200.                                    c
c                                                                             c
c ZDCF peak at  +0.133     r_max =  +0.881     , ZDCF C.O.M at   -4.31        c
c                                                                             c
c  #      lag         r         -dr        +dr    likelihood                  c
c --- ---------- ---------- ---------- ---------- ----------                  c
c  1  -76.8      0.355      0.314      0.275      2.701E-03                   c
c  2  -19.8      0.780      0.121      9.680E-02  0.111                       c
c  3  0.133      0.881      8.032E-02  6.053E-02  0.458                       c
c  4   5.62      0.821      0.124      9.322E-02  0.205                       c
c  5   24.1      0.392      0.306      0.265      3.669E-03                   c
c  6   87.5     -0.162      0.319      0.339      2.365E-05                   c
c  7   200.     -0.435      0.251      0.295      7.652E-07                   c
c                                                                             c
c ML Peak at  +0.133     , L at peak =  +0.458                                c
c 1 sigma ML interval =  +0.133       +10.0       -23.7                       c
c                     = (  -23.6     ,   +10.2    )                           c
c                                                                             c
c Program ended.                                                              c
c                                                                             c
c EXPLANATION                                                                 c
c -----------                                                                 c
c                                                                             c
c o The program uses the *.dcf output file of the ZDCF.                       c
c                                                                             c
c o The user has to give lower and upper bounds on the peak location          c
c                                                                             c
c o The output consists of the likelihood function and the fiducial interval, c
c   given both as t_peak + dt_up - dt_low, and as the interval (t_min, t_max) c
c                                                                             c
c FORTRAN ISSUES                                                              c
c --------------                                                              c
c                                                                             c
c o Compile with static variables flag:                                       c
c      Intel Fortran: ifort -save -o plike plike_v3.0.f                       c
c      g77 Fortran: g77 -fno-automatic  -o plike plike_v3.0.f                 c
c      g95 Fortran: g95 -fstatic  -o plike plike_v3.0.f                       c
c      Other compilers: -static (?)                                           c
c                                                                             c
c o The program can be VERY slow - don't panic!                               c
c                                                                             c
c o This program incorporates the functions QROMO, POLINT and MIDPNT          c
c   from Numerical Recipes / Press, Flannery, Teukolsky & Vetterling.         c
c                                                                             c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer i,np,maxp
      parameter (maxp = 200)
      integer n(maxp)
      real t(maxp),dtm(maxp),dtp(maxp),r(maxp),drm(maxp),drp(maxp)
      real liklhd(maxp),mxlik,tmxlik
      real tlow,thigh,dtlikl,dtliku,tcom
      character*72 dcf
      integer oerr,cerr
      real com
c
      print *,'PLIKE V3.0 begins.'
c
      print *,'Enter dcf file name:'
      read '(a72)',dcf
      open (unit=11,file=dcf,status='OLD',iostat=oerr,err=901)
c
      print *,'Enter lower bound on peak location:'
      read *,tlow
      print *,'Enter upper bound on peak location:'
      read *,thigh
      print *,'Calculating ML in range ',tlow,' to ', thigh
      np = 0
1     continue
         if (np .ge. maxp) then
            print *,'Only 1st ',maxp,' points were read!'
            goto 2
         endif
         read (11,*,end=2)
     >      t(np+1),dtm(np+1),dtp(np+1),
     >      r(np+1),drm(np+1),drp(np+1),n(np+1)
         if (t(np+1) .lt. tlow .or. t(np+1) .gt. thigh) goto 1
         np = np+1
         goto 1
2     continue
      close (unit=11,status='KEEP',iostat=cerr,err=902)
c
      print *
      tcom = com(t,r,np,i)
      print '(''ZDCF peak at '',sp,1p,g11.3,'' r_max = '',g11.3,
     >        '' , ZDCF C.O.M at '',g11.3)',t(i),r(i),tcom
      call plike3(t,dtm,dtp,r,drm,drp,n,np,tlow,thigh,
     >            liklhd,tmxlik,mxlik,dtlikl,dtliku,2)
c
      print *
      print *,
     >' #      lag         r         -dr        +dr    likelihood '
      print *,
     >'--- ---------- ---------- ---------- ---------- ---------- '
      print '(i3,'' '',1p,g10.3,1x,g10.3,1x,g10.3,1x,g10.3,1x,g10.3)',
     >   (i,t(i),r(i),drm(i),drp(i),liklhd(i),i=1,np)
c
      print *
      print '(''ML Peak at '',sp,1p,g11.3,'' , L at peak = '',g11.3)',
     >   tmxlik,mxlik
      print '(''1 sigma ML interval = '',sp,1p,g11.3,1x,g11.3,1x,
     >        g11.3,/,19x,'' = ('',g11.3,'' , '',g11.3,'')'')',
     >   tmxlik,dtliku,-dtlikl,tmxlik-dtlikl,tmxlik+dtliku
      print *
      print *,'Program ended.'
      stop
 901  print *,'ERROR ON OPENING FILE - ERRCOD=',oerr
      stop
 902  print *,'ERROR ON CLOSING FILE - ERRCOD=',cerr
      stop
      end
c//////////////////////////////////////////////////////////////////////////////
c
c com calculates the center-of-mass (centroid) around the ZDCF peak.
c The sum is taken over all the points lying between the two nearest
c the maximum with r < r_max/2.
c
      real function com(t,r,n,imax)
      implicit none
      integer n,imax,nerr
      real t(n),r(n)
      real sumr,sumtr,rhalf
      integer i,i1,i2
c Locating the maximum
      imax = 1
      do i = 1,n
         if (r(imax) .lt. r(i)) imax = i
      enddo
c
      rhalf = r(imax)/2.
c Locating the peak
      i1 = imax-1
      do while (i1.ge.1 .and. r(i1).ge.rhalf)
         i1 = i1-1
      enddo
      i1 = i1+1
c Warning
      if (i1.eq.1 .and. r(1).ge.rhalf) then
         print *,
     >      'COM: Warning - Peak extends beyond range to the left...'
      endif
c
      i2 = imax+1
      do while (i2.le.n .and. r(i2).ge.rhalf)
         i2 = i2+1
      enddo
      i2 = i2-1
c Warning
      if (i2.eq.n .and. r(n).ge.rhalf) then
         print *,
     >      'COM: Warning - Peak extends beyond range to the right...'
         nerr = nerr+1
      endif
c Calculating the center of mass
      sumtr = 0.0
      sumr = 0.0
      do i = i1,i2
         sumtr = sumtr+t(i)*r(i)
         sumr = sumr+r(i)
      enddo
      if (sumr .le. 0.0) then
         print *,'COM: Problems. n = ',n,' imax = ',imax,
     >           ' rmax = ',r(imax),
     >           ' i1 = ',i1,' i2 = ',i2,' sumr = ',sumr
         print '(1x,f10.3,1x,f10.3)',(t(i),r(i),i=1,n)
         stop
      else
         com = sumtr/sumr
      endif
      return
      end
c//////////////////////////////////////////////////////////////////////////////
cold      include "/home/tal/ZDCF/src/plike_func3.f"
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                             c
c                  Calculating the Peak Likelihood of the ZDCF                c
c                                                                             c
c        A special version for processing the ZDCF output file *.dcf          c
c                                                                             c
c                                Version 3.0                                  c
c                                                                             c
c    (Bayesian prior assumption of uniform distribution in z-space)           c
c                                                                             c
c                                                                             c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                             c
c SAMPLE RUN (user input marked by "<-USER". See explanations below)          c
c ----------                                                                  c
c                                                                             c
c PLIKE V3.0 begins.                                                          c
c Enter dcf file name:                                                        c
c test.dcf      <-USER                                                        c
c Enter lower bound on peak location:                                         c
c -200           <-USER                                                       c
c Enter upper bound on peak location:                                         c
c 200            <-USER                                                       c
c Working on   7 points                                                       c
c                                                                             c
c   i       tau        z          s_z       likelihood                        c
c -----------------------------------------------------                       c
c    1  -76.77000    0.39168    0.35056     3.2769E-03                        c
c    2  -19.79000    1.07518    0.28485     0.1270                            c
c    3    0.13330    1.42651    0.32573     0.5923                            c
c    4    5.61800    1.20773    0.34551     0.2729                            c
c    5   24.10000    0.43595    0.35034     4.4965E-03                        c
c    6   87.48000   -0.17234    0.35143     2.2088E-05                        c
c    7  199.89999   -0.49122    0.35009     4.8453E-07                        c
c                                                                             c
c Peak at 0.133     , ML(peak) = 0.592                                        c
c 1 sigma ML interval =  0.133    + 13.4    -  23.9                           c
c                     = ( -23.8     ,   13.5    )                             c
c                                                                             c
c Program ended.                                                              c
c                                                                             c
c EXPLANATION                                                                 c
c -----------                                                                 c
c                                                                             c
c o The program uses the *.dcf output file of the ZDCF.                       c
c                                                                             c
c o The user has to give lower and upper bounds on the peak location          c
c                                                                             c
c o The output consists of the likelihood function and the fiducial interval, c
c   given both as t_peak + dt_up - dt_low, and as the interval (t_min, t_max) c
c                                                                             c
c FORTRAN ISSUES                                                              c
c --------------                                                              c
c                                                                             c
c o The program is VERY slow - don't panic!                                   c
c                                                                             c
c o This program incorporates the functions QROMO, POLINT and MIDPNT          c
c   from Numerical Recipes / Press, Flannery, Teukolsky & Vetterling.         c
c                                                                             c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine plike3(t,dtm,dtp,r,drm,drp,n,np0,tlow,thigh,
c                      i i   i   i i  i  i i i   i    i
     >                 liklhd,tmxlik,mxlik,dtlikl,dtliku,zrflag)
c                      o      o      o     o      o      i
      implicit none
      integer np0,zrflag
      real tlow,thigh,tmxlik,mxlik,dtlikl,dtliku
      real t(np0),dtm(np0),dtp(np0),r(np0),drm(np0),drp(np0),liklhd(np0)
      integer n(np0)
c
      integer i,maxp,np,ilow
      parameter (maxp = 200)
      integer pnt(maxp)
      real z(maxp),pcum(maxp)
      real p,rpeak
      integer imxlik,ipeak
      real eps,epslik,tiny,maxsig
      parameter (eps = 1e-6, epslik = 1e-10, tiny = 1e-30, maxsig = 5)
      real likint,likin0,likin1
      external likint,midpt1
c
      np = 0
      ilow = 0
      rpeak = -1.0
c     Input given in z-space (zrflag = 1): vector r holds z
      if (zrflag .eq. 1) then
         do i = 1,np0
            if (t(i) .ge. tlow .and. t(i) .le. thigh) then
               np = np+1
               if (np .gt. maxp) then
                  print *,'PLIKE: Work area too small ! ',
     >                 'Increase maxp to >= ',maxp
                  stop
               endif
               z(np) = r(i)
               liklhd(np) = 0.0
               if (r(i) .gt. rpeak) then
                  ipeak = np
                  rpeak = r(i)
               endif
            else if (t(i) .lt. tlow) then
               ilow = i
            endif
         enddo
c     Input given in r-space (rzflag <> 1)
      else
         do i = 1,np0
            if (t(i) .ge. tlow .and. t(i) .le. thigh) then
               np = np+1
               if (np .gt. maxp) then
                  print *,'PLIKE: Work area too small ! ',
     >                 'Increase maxp to >= ',maxp
                  stop
               endif
               z(np)  = log((1.+r(np))/(1.-r(np)))/2
               liklhd(np) = 0.0
               if (r(i) .gt. rpeak) then
                  ipeak = np
                  rpeak = r(i)
               endif
            else if (t(i) .lt. tlow) then
               ilow = i
            endif
         enddo
      endif
c Creating a vector of pointers to cells in t, dtm, dtp which relate to z
      do i = 1,np
         pnt(i) = ilow+i
      enddo
      ilow = ilow+1
c
c Calculating the likelihood function
c
      p = likin0(z,n,np)
      do i = ipeak-1,1,-1
         p = likin1(i)
         print '(''.'',$)'
c Note small modifications of original qromo
         call qromo1(likint,-1d0,1d0,liklhd(i),midpt1)
cdbg(
c         print '(1x,i3,1x,0p,f10.5,1x,f10.5,5x,i3,7x,1p,g10.4)',
c     >         i,t(i),z(i),n(i),liklhd(i)
cdbg)
         if (liklhd(i) .lt. epslik) goto 301
      enddo
 301  continue
      do i = ipeak,np
         p = likin1(i)
         print '(''.'',$)'
c Note small modifications of original qromo
         call qromo1(likint,-1d0,1d0,liklhd(i),midpt1)
cdbg(
c         print '(1x,i3,1x,0p,f10.5,1x,f10.5,5x,i3,7x,1p,g10.4)',
c     >         i,t(i),z(i),n(i),liklhd(i)
cdbg)
         if (liklhd(i) .lt. epslik) goto 302
      enddo
 302  continue
      print '(/)'
c
c Finding the maximal likelihood value and error interval
c
      call conf1(liklhd,t(ilow),dtm(ilow),dtp(ilow),np,
     >           pcum,imxlik,tmxlik,mxlik,dtlikl,dtliku)
c
cdbg      print *
cdbg      print '(''Peak at '',1p,g9.3,'' , ML(peak) = '',g9.3)',
cdbg     >   tmxlik,mxlik
cdbg      print '(''1 sigma ML interval = '',1p,g10.3,''+'',g9.3,''- '',g9.3,
cdbg     >        /,19x,'' = ('',g10.3,'' , '',g10.3,'')'')',
cdbg     >   tmxlik,t1sigu-tmxlik,tmxlik-t1sigl,t1sigl,t1sigu
c
      return
      end
c//////////////////////////////////////////////////////////////////////////////
c
c The integrand of the likelyhood function
c
      real function likint(rho)
      implicit none
      integer maxp
      parameter (maxp = 200)
      real*8 rho
      real z0(*),z(maxp),cn,zbar,sz,dzdr
      integer n0(*),n(maxp)
      integer i0,i,np0,np
c cpi = 1/sqrt(2*Pi)
      integer j
      real cpi,tiny
      parameter (cpi = 0.398942280402, tiny = 1e-10)
      real clike,exp,log
      real likin0,likin1
c
cdbg(
      if (abs(rho) .ge. 1.0) then
	print *,'LIKINT: |rho| >= 1'
	stop
      endif
cdbg)
      zbar = log((1.+rho)/(1.-rho))/2.+rho/2./(n(i)-1)*
     >         (1.+(5.+rho**2)/4./(n(i)-1.)+
     >             (11.+2.*rho**2+3.*rho**4)/8./(n(i)-1)**2)
      sz = sqrt(1./(n(i)-1)*(1.+(4-rho**2)/2./(n(i)-1)+
     >                       (22.-6*rho**2-3*rho**4)/6./(n(i)-1)**2))
      dzdr = 1./(1.+rho)/(1.-rho)
      likint = log(dzdr)+log(cpi)-log(sz)-((z(i)-zbar)/sz)**2/2
      do j = 1,i-1
         cn = clike(rho,z(j),n(j))
         if (cn .le. tiny) then
            likint = 0.0
            return
         endif
         likint = likint+log(cn)
      enddo
      do j = i+1,np
         cn = clike(rho,z(j),n(j))
         if (cn .le. tiny) then
            likint = 0.0
            return
         endif
         likint = likint+log(cn)
      enddo
      likint = exp(likint)
      return
c
      entry likin0(z0,n0,np0)
      if (np .gt. maxp) then
         print *,'LIKINT: maxp too small. np = ',np
         stop
      endif
      np = np0
      do j = 1,np
         z(j) = z0(j)
         n(j) = n0(j)
      enddo
      return
c
      entry likin1(i0)
      i = i0
      return
      end
c///////////////////////////////////////////////////////////////////////
c Integrating over the likelihood function from rho=-1 to rho=rho0 for a
c given value of the estimator r evaluated on n pairs.
      real function clike(rho0,z,n)
      implicit none
      real*8 rho0
      real z
      integer n
c
      real phiz0,phiz1
      external midpt2,phiz1
c
      clike = phiz0(n,z)
cdbg(
      if (abs(rho0) .gt. 1.0) then
	print *,'CLIKE: |rho| > 1'
	stop
      endif
cdbg)

c Note small modifications of original qromo
      call qromo2(phiz1,-1d0,rho0,clike,midpt2)
      return
      end
c///////////////////////////////////////////////////////////////////////
c The integrand of the cumulative likelihood function
c exp(-(z(r)-zbar(rho))/sz(rho))**2/2)/sqrt(2Pi)/sz(rho)
      real function phiz1(rho)
      implicit none
      integer n0,n
      real z0,z
      real*8 rho
      real dzdr,zbar,sz
c sqrpi2 = 1/sqrt(2*Pi)
      real sqrpi2
      parameter (sqrpi2 = .39894228040143267794)
      real sqrt,exp,phiz0
c
cdbg(
      if (abs(rho) .ge. 1.0) then
	print *,'PHIZ1: |rho| >= 1'
	stop
      endif
cdbg)
      zbar = log((1.+rho)/(1.-rho))/2.+rho/2./(n-1)*
     >         (1.+(5.+rho**2)/4./(n-1.)+
     >             (11.+2.*rho**2+3.*rho**4)/8./(n-1)**2)
      sz = sqrt(1./(n-1)*(1.+(4-rho**2)/2./(n-1)+
     >                       (22.-6*rho**2-3*rho**4)/6./(n-1)**2))
      dzdr = 1./(1.+rho)/(1.-rho)
      phiz1 = dzdr*sqrpi2/sz*exp(-((z-zbar)/sz)**2/2)
cdbg      phiz1 = log(dzdr*sqrpi2)-log(sz)-((z-zbar)/sz)**2/2
cdbg      phiz1 = exp(phiz1)
      return
c
      entry phiz0(n0,z0)
      n = n0
      z = z0
      phiz0 = 0.0
      return
      end
c//////////////////////////////////////////////////////////////////////////////
c Finding the +/- 1 sigma intervals on the interpolated likelihood function
c which (after normalization) is the "Fiducial Probability Function" (Kendall
c & Stuart, Vol 2, 3rd ed. p. 141)
      subroutine conf1(liklhd,t,dtm,dtp,np,pcum,imxlik,tmxlik,mxlik,
c                      i      i i   i   i  o    o      o      o
     >                 dtlikl,dtliku)
c                      o      o
      implicit none
      integer np,imxlik
      real liklhd(np),t(np),dtm(np),dtp(np),pcum(np)
      real tmxlik,mxlik,dtlikl,dtliku
      integer i
      real cn,s,p,t1sigl,t1sigu
c
      mxlik = -1.0
      do i = 1,np
         if (mxlik .lt. liklhd(i)) then
            tmxlik = t(i)
            imxlik = i
            mxlik = liklhd(i)
         endif
      enddo
      cn = 0.0
      pcum(1) = 0.0
      do i = 2,np
         cn = cn+(liklhd(i)+liklhd(i-1))/2.0*
     >           (t(i)-t(i-1))
         pcum(i) = cn
      enddo
      do i = 1,np
         pcum(i) = pcum(i)/cn
      enddo
c
      p = pcum(imxlik)+0.3413*(pcum(np)-pcum(imxlik))*2
      do i = imxlik,np
         if (pcum(i) .ge. p) then
            s = (liklhd(i)-liklhd(i-1))/(t(i)-t(i-1))
            t1sigu = t(i-1)+
     >             (-liklhd(i-1)+
     >              sqrt(liklhd(i-1)**2+2*(p-pcum(i-1))*s*cn))/s
            goto 100
         endif
      enddo
cdbg
      print *,'t1sigu = end point...'
      t1sigu = t(np)
 100  continue
c
      p = pcum(imxlik)-0.3413*(pcum(imxlik)-pcum(1))*2
      do i = imxlik,1,-1
         if (pcum(i) .le. p) then
            s = (liklhd(i+1)-liklhd(i))/(t(i+1)-t(i))
            t1sigl = t(i)+
     >             (-liklhd(i)+
     >              sqrt(liklhd(i)**2+2*(p-pcum(i))*s*cn))/s
            goto 200
         endif
      enddo
cdbg
      print *,'t1sigl = end point...'
      t1sigl = t(1)
 200  continue
c
c???      dtliku = amax0(t1sigu-tmxlik,dtp(imxlik))
c???      dtlikl = amax0(tmxlik-t1sigl,dtm(imxlik))
      dtliku = amax1(t1sigu-tmxlik,dtp(imxlik))
      dtlikl = amax1(tmxlik-t1sigl,dtm(imxlik))
cdbg
cdbg      print *,'CONF1: t1sigl = ',t1sigl,' t1sigu = ',t1sigu
cdbg
      return
      end
c///////////////////////////////////////////////////////////////////////
      SUBROUTINE QROMO1(FUNC,A,B,SS,CHOOSE)
      save all
      PARAMETER (EPS=1.E-6,JMAX=14,JMAXP=JMAX+1,KM=4,K=KM+1)
      DIMENSION S(JMAXP),H(JMAXP)
cTal(
      real*8 a,b
cTal)
      H(1)=1.
      DO 11 J=1,JMAX
        CALL CHOOSE(FUNC,A,B,S(J),J)
        IF (J.GE.K) THEN
          CALL POLIN1(H(J-KM),S(J-KM),K,0.0,SS,DSS)
cTal          IF (ABS(DSS).LT.EPS*ABS(SS)) RETURN
          IF (ABS(DSS).LE.EPS*ABS(SS)) RETURN
        ENDIF
        S(J+1)=S(J)
        H(J+1)=H(J)/9.
11    CONTINUE
      print *,'QROMO1: Too many steps...'
      END
c///////////////////////////////////////////////////////////////////////
      SUBROUTINE MIDPT1(FUNC,A,B,S,N)
cTal(
      save all
      real*8 a,b,tnm,del,ddel,x,sum
cTal)
      IF (N.EQ.1) THEN
        S=(B-A)*FUNC(0.5*(A+B))
        IT=1
      ELSE
        TNM=IT
        IF(TNM.EQ.0.)PAUSE 'MIDPT1: attempting to divide by den = 0'
        DEL=(B-A)/(3.*TNM)
        DDEL=DEL+DEL
        X=A+0.5*DEL
        SUM=0.
        DO 11 J=1,IT
          SUM=SUM+FUNC(X)
          X=X+DDEL
          SUM=SUM+FUNC(X)
          X=X+DEL
11      CONTINUE
        S=(S+(B-A)*SUM/TNM)/3.
        IT=3*IT
      ENDIF
      RETURN
      END
c//////////////////////////////////////////////////////////////////////////////
      SUBROUTINE POLIN1(XA,YA,N,X,Y,DY)
      save all
      PARAMETER (NMAX=10)
      DIMENSION XA(N),YA(N),C(NMAX),D(NMAX)
      NS=1
      DIF=ABS(X-XA(1))
      DO 11 I=1,N
        DIFT=ABS(X-XA(I))
        IF (DIFT.LT.DIF) THEN
          NS=I
          DIF=DIFT
        ENDIF
        C(I)=YA(I)
        D(I)=YA(I)
11    CONTINUE
      Y=YA(NS)
      NS=NS-1
      DO 13 M=1,N-1
        DO 12 I=1,N-M
          HO=XA(I)-X
          HP=XA(I+M)-X
          W=C(I+1)-D(I)
          DEN=HO-HP
          IF(DEN.EQ.0.)PAUSE 'POLIN1: attempting to divide by den = 0'
          DEN=W/DEN
          D(I)=HP*DEN
          C(I)=HO*DEN
12      CONTINUE
        IF (2*NS.LT.N-M)THEN
          DY=C(NS+1)
        ELSE
          DY=D(NS)
          NS=NS-1
        ENDIF
        Y=Y+DY
13    CONTINUE
      RETURN
      END
c///////////////////////////////////////////////////////////////////////
      SUBROUTINE QROMO2(FUNC,A,B,SS,CHOOSE)
      save all
cTal      PARAMETER (EPS=1.E-6,JMAX=14,JMAXP=JMAX+1,KM=4,K=KM+1)
      PARAMETER (EPS=1.E-6,JMAX=15,JMAXP=JMAX+1,KM=4,K=KM+1)
      DIMENSION S(JMAXP),H(JMAXP)
cTal(
      real*8 a,b
cTal)
      H(1)=1.
      DO 11 J=1,JMAX
        CALL CHOOSE(FUNC,A,B,S(J),J)
        IF (J.GE.K) THEN
          CALL POLIN2(H(J-KM),S(J-KM),K,0.0,SS,DSS)
cTal          IF (ABS(DSS).LT.EPS*ABS(SS)) RETURN
          IF (ABS(DSS).LE.EPS*ABS(SS)) RETURN
        ENDIF
        S(J+1)=S(J)
        H(J+1)=H(J)/9.
11    CONTINUE
      print *,'QROMO2: Too many steps...'
      END
c///////////////////////////////////////////////////////////////////////
      SUBROUTINE MIDPT2(FUNC,A,B,S,N)
cTal(
      save all
      real*8 a,b,tnm,del,ddel,x,sum
cTal)
      IF (N.EQ.1) THEN
        S=(B-A)*FUNC(0.5*(A+B))
        IT=1
      ELSE
        TNM=IT
        IF(TNM.EQ.0.)PAUSE 'MIDPT2: attempting to divide by den = 0'
        DEL=(B-A)/(3.*TNM)
        DDEL=DEL+DEL
        X=A+0.5*DEL
        SUM=0.
        DO 11 J=1,IT
          SUM=SUM+FUNC(X)
          X=X+DDEL
          SUM=SUM+FUNC(X)
          X=X+DEL
11      CONTINUE
        S=(S+(B-A)*SUM/TNM)/3.
        IT=3*IT
      ENDIF
      RETURN
      END
c//////////////////////////////////////////////////////////////////////////////
      SUBROUTINE POLIN2(XA,YA,N,X,Y,DY)
      save all
      PARAMETER (NMAX=10)
      DIMENSION XA(N),YA(N),C(NMAX),D(NMAX)
      NS=1
      DIF=ABS(X-XA(1))
      DO 11 I=1,N
        DIFT=ABS(X-XA(I))
        IF (DIFT.LT.DIF) THEN
          NS=I
          DIF=DIFT
        ENDIF
        C(I)=YA(I)
        D(I)=YA(I)
11    CONTINUE
      Y=YA(NS)
      NS=NS-1
      DO 13 M=1,N-1
        DO 12 I=1,N-M
          HO=XA(I)-X
          HP=XA(I+M)-X
          W=C(I+1)-D(I)
          DEN=HO-HP
          IF(DEN.EQ.0.)PAUSE 'POLIN2: attempting to divide by den = 0'
          DEN=W/DEN
          D(I)=HP*DEN
          C(I)=HO*DEN
12      CONTINUE
        IF (2*NS.LT.N-M)THEN
          DY=C(NS+1)
        ELSE
          DY=D(NS)
          NS=NS-1
        ENDIF
        Y=Y+DY
13    CONTINUE
      RETURN
      END
c//////////////////////////////////////////////////////////////////////////////

max_stars_repo_name: orelgueta/blazar-variability-study

max_issues_repo_name: orelgueta/blazar-variability-study

max_forks_repo_name: orelgueta/blazar-variability-study


MODULE Soilevaporationmod
    IMPLICIT NONE
CONTAINS

    SUBROUTINE model_soilevaporation(diffusionLimitedEvaporation, &
        energyLimitedEvaporation, &
        soilEvaporation)
        IMPLICIT NONE
        REAL, INTENT(IN) :: diffusionLimitedEvaporation
        REAL, INTENT(IN) :: energyLimitedEvaporation
        REAL, INTENT(OUT) :: soilEvaporation
        !- Description:
    !            * Title: SoilEvaporation Model
    !            * Author: Pierre Martre
    !            * Reference: Modelling energy balance in the wheat crop model SiriusQuality2:
    !            Evapotranspiration and canopy and soil temperature calculations
    !            * Institution: INRA Montpellier
    !            * Abstract: Starting from a soil at field capacity, soil evaporation  is assumed to
    !                be energy limited during the first phase of evaporation and diffusion limited thereafter.
    !                Hence, the soil evaporation model considers these two processes taking the minimum between
    !                the energy limited evaporation (PtSoil) and the diffused limited
    !                evaporation
        !- inputs:
    !            * name: diffusionLimitedEvaporation
    !                          ** description : diffusion Limited Evaporation
    !                          ** variablecategory : state
    !                          ** datatype : DOUBLE
    !                          ** default : 6605.505
    !                          ** min : 0
    !                          ** max : 10000
    !                          ** unit : g m-2 d-1
    !                          ** uri : http://www1.clermont.inra.fr/siriusquality/?page_id=547
    !                          ** inputtype : variable
    !            * name: energyLimitedEvaporation
    !                          ** description : energy Limited Evaporation
    !                          ** variablecategory : state
    !                          ** datatype : DOUBLE
    !                          ** default : 448.240
    !                          ** min : 0
    !                          ** max : 1000
    !                          ** unit : g m-2 d-1
    !                          ** uri : http://www1.clermont.inra.fr/siriusquality/?page_id=547
    !                          ** inputtype : variable
        !- outputs:
    !            * name: soilEvaporation
    !                          ** description : soil Evaporation
    !                          ** variablecategory : auxiliary
    !                          ** datatype : DOUBLE
    !                          ** min : 0
    !                          ** max : 5000
    !                          ** unit : g m-2 d-1
    !                          ** uri : http://www1.clermont.inra.fr/siriusquality/?page_id=547
        soilEvaporation = min(diffusionLimitedEvaporation, energyLimitedEvaporation)
    END SUBROUTINE model_soilevaporation

END MODULE


text:
bil = 4*Pgas*cgs_mhydr*cgs_stef*miu*(T - Trad)*(4*Trad**4*cgs_boltz*F(1, &
      2)/(cgs_c**2*cgs_mel) + (T + Trad)*(T**2 + Trad**2)*F(1, 1))/(T* &
      cgs_boltz*heat) - 1
bil_T = 4*Pgas*cgs_mhydr*cgs_stef*miu*(T**5*cgs_c**2*cgs_mel*F(2, 1) + 3 &
      *T**4*cgs_c**2*cgs_mel*F(1, 1) + 4*T**2*Trad**4*cgs_boltz*F(2, 2 &
      ) - 4*T*Trad**5*cgs_boltz*F(2, 2) - T*Trad**4*cgs_c**2*cgs_mel*F( &
      2, 1) + 4*Trad**5*cgs_boltz*F(1, 2) + Trad**4*cgs_c**2*cgs_mel*F( &
      1, 1))/(T**2*cgs_boltz*cgs_c**2*cgs_mel*heat)

max_stars_repo_name: gronki/pydiskvert

max_issues_repo_name: gronki/pydiskvert

max_forks_repo_name: gronki/pydiskvert


text:
c     DRSLESUM
c>> 1995-05-28 DRSLESUM Krogh  Changes to use M77CON
c>> 1994-08-09 DRSLESUM WVS  Set up for CHGTYP
c>> 1994-07-14 DRSLESUM CLL
c>> 1992-04-29 DRSLESUM CAO Replaced '1' in format.
c>> 1991-11-19 DRSLESUM CLL
c>> 1987-12-09 DRSLESUM Lawson  Initial Code.
c--S replaces "?": ?LESUM, DR?LESUM
c
c     Demonstration driver for evaluation of a Legendre series.
c     ------------------------------------------------------------------
      integer j
      real             x,a(0:5),y,w,z
      data a/0.07e0, 0.27e0, 0.20e0, 0.28e0, 0.08e0, 0.08e0/
c     ------------------------------------------------------------------
      print '(1x,3x,a1,14x,a1,17x,a1/)','x','y','z'
      do 20 j = -10,10,2
        x = real(j) /10.e0
        call slesum (x,5,a,y)
        w = 0.35e0 * (x**4) + 0.63e0 * (x**5)
        z = y - w
        print '(1x,f5.2,5x,g15.7,g15.2)',x,y,z
   20 continue
      end

ext: f


text:
write(*,*)"a"
end

max_stars_repo_name: xevolver/xevxml

max_issues_repo_name: xevolver/xevxml

max_forks_repo_name: xevolver/xevxml


text:
c
C === dummy routines with ADIFOR with SparsLinC option
C     used when you cannot access the SparsLinC library from ADIFOR
C
      subroutine xspini
      return
      end
      subroutine dspxsq(indvec,ajac,num,iwk,joutlen,ires)
      integer indvec, num, iwk, joutlen, ires
      real*8  ajac
      return
      end
      subroutine  dspsd(indvec, i, a, n)
      integer indvec, i, n
      real*8 a
      return
      end
      subroutine ehrpt
      return
      end

max_stars_repo_name: davidbrochart/daetk

max_issues_repo_name: davidbrochart/daetk

max_forks_repo_name: davidbrochart/daetk


text:
*DECK PJAC
      SUBROUTINE PJAC (NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM, F,
     +   JAC, RPAR, IPAR)
C***BEGIN PROLOGUE  PJAC
C***SUBSIDIARY
C***PURPOSE  Subsidiary to DEBDF
C***LIBRARY   SLATEC
C***TYPE      SINGLE PRECISION (PJAC-S, DPJAC-D)
C***AUTHOR  Watts, H. A., (SNLA)
C***DESCRIPTION
C
C   PJAC sets up the iteration matrix (involving the Jacobian) for the
C   integration package DEBDF.
C
C***SEE ALSO  DEBDF
C***ROUTINES CALLED  SGBFA, SGEFA, VNWRMS
C***COMMON BLOCKS    DEBDF1
C***REVISION HISTORY  (YYMMDD)
C   800901  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900328  Added TYPE section.  (WRB)
C   910722  Updated AUTHOR section.  (ALS)
C   920422  Changed DIMENSION statement.  (WRB)
C***END PROLOGUE  PJAC
C
CLLL. OPTIMIZE
      INTEGER NEQ, NYH, IWM, I, I1, I2, IER, II, IOWND, IOWNS, J, J1,
     1   JJ, JSTART, KFLAG, L, LENP, MAXORD, MBA, MBAND, MEB1, MEBAND,
     2   METH, MITER, ML, ML3, MU, N, NFE, NJE, NQ, NQU, NST
      EXTERNAL F, JAC
      REAL Y, YH, EWT, FTEM, SAVF, WM,
     1   ROWND, ROWNS, EL0, H, HMIN, HMXI, HU, TN, UROUND,
     2   CON, DI, FAC, HL0, R, R0, SRUR, YI, YJ, YJJ, VNWRMS
      DIMENSION         Y(*), YH(NYH,*), EWT(*), FTEM(*), SAVF(*),
     1   WM(*), IWM(*), RPAR(*), IPAR(*)
      COMMON /DEBDF1/ ROWND, ROWNS(210),
     1   EL0, H, HMIN, HMXI, HU, TN, UROUND, IOWND(14), IOWNS(6),
     2   IER, JSTART, KFLAG, L, METH, MITER, MAXORD, N, NQ, NST, NFE,
     3   NJE, NQU
C-----------------------------------------------------------------------
C PJAC IS CALLED BY STOD  TO COMPUTE AND PROCESS THE MATRIX
C P = I - H*EL(1)*J , WHERE J IS AN APPROXIMATION TO THE JACOBIAN.
C HERE J IS COMPUTED BY THE USER-SUPPLIED ROUTINE JAC IF
C MITER = 1 OR 4, OR BY FINITE DIFFERENCING IF MITER = 2, 3, OR 5.
C IF MITER = 3, A DIAGONAL APPROXIMATION TO J IS USED.
C J IS STORED IN WM AND REPLACED BY P.  IF MITER .NE. 3, P IS THEN
C SUBJECTED TO LU DECOMPOSITION IN PREPARATION FOR LATER SOLUTION
C OF LINEAR SYSTEMS WITH P AS COEFFICIENT MATRIX. THIS IS DONE
C BY SGEFA IF MITER = 1 OR 2, AND BY SGBFA IF MITER = 4 OR 5.
C
C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION
C WITH PJAC USES THE FOLLOWING..
C Y    = ARRAY CONTAINING PREDICTED VALUES ON ENTRY.
C FTEM = WORK ARRAY OF LENGTH N (ACOR IN STOD ).
C SAVF = ARRAY CONTAINING F EVALUATED AT PREDICTED Y.
C WM   = REAL WORK SPACE FOR MATRICES.  ON OUTPUT IT CONTAINS THE
C        INVERSE DIAGONAL MATRIX IF MITER = 3 AND THE LU DECOMPOSITION
C        OF P IF MITER IS 1, 2 , 4, OR 5.
C        STORAGE OF MATRIX ELEMENTS STARTS AT WM(3).
C        WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA..
C        WM(1) = SQRT(UROUND), USED IN NUMERICAL JACOBIAN INCREMENTS.
C        WM(2) = H*EL0, SAVED FOR LATER USE IF MITER = 3.
C IWM  = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING AT
C        IWM(21), IF MITER IS 1, 2, 4, OR 5.  IWM ALSO CONTAINS THE
C        BAND PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER IS 4 OR 5.
C EL0  = EL(1) (INPUT).
C IER  = OUTPUT ERROR FLAG,  = 0 IF NO TROUBLE, .NE. 0 IF
C        P MATRIX FOUND TO BE SINGULAR.
C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, TN, UROUND,
C MITER, N, NFE, AND NJE.
C-----------------------------------------------------------------------
C***FIRST EXECUTABLE STATEMENT  PJAC
      NJE = NJE + 1
      HL0 = H*EL0
      GO TO (100, 200, 300, 400, 500), MITER
C IF MITER = 1, CALL JAC AND MULTIPLY BY SCALAR. -----------------------
 100  LENP = N*N
      DO 110 I = 1,LENP
 110    WM(I+2) = 0.0E0
      CALL JAC (TN, Y, WM(3), N, RPAR, IPAR)
      CON = -HL0
      DO 120 I = 1,LENP
 120    WM(I+2) = WM(I+2)*CON
      GO TO 240
C IF MITER = 2, MAKE N CALLS TO F TO APPROXIMATE J. --------------------
 200  FAC = VNWRMS (N, SAVF, EWT)
      R0 = 1000.0E0*ABS(H)*UROUND*N*FAC
      IF (R0 .EQ. 0.0E0) R0 = 1.0E0
      SRUR = WM(1)
      J1 = 2
      DO 230 J = 1,N
        YJ = Y(J)
        R = MAX(SRUR*ABS(YJ),R0*EWT(J))
        Y(J) = Y(J) + R
        FAC = -HL0/R
        CALL F (TN, Y, FTEM, RPAR, IPAR)
        DO 220 I = 1,N
 220      WM(I+J1) = (FTEM(I) - SAVF(I))*FAC
        Y(J) = YJ
        J1 = J1 + N
 230    CONTINUE
      NFE = NFE + N
C ADD IDENTITY MATRIX. -------------------------------------------------
 240  J = 3
      DO 250 I = 1,N
        WM(J) = WM(J) + 1.0E0
 250    J = J + (N + 1)
C DO LU DECOMPOSITION ON P. --------------------------------------------
      CALL SGEFA (WM(3), N, N, IWM(21), IER)
      RETURN
C IF MITER = 3, CONSTRUCT A DIAGONAL APPROXIMATION TO J AND P. ---------
 300  WM(2) = HL0
      IER = 0
      R = EL0*0.1E0
      DO 310 I = 1,N
 310    Y(I) = Y(I) + R*(H*SAVF(I) - YH(I,2))
      CALL F (TN, Y, WM(3), RPAR, IPAR)
      NFE = NFE + 1
      DO 320 I = 1,N
        R0 = H*SAVF(I) - YH(I,2)
        DI = 0.1E0*R0 - H*(WM(I+2) - SAVF(I))
        WM(I+2) = 1.0E0
        IF (ABS(R0) .LT. UROUND*EWT(I)) GO TO 320
        IF (ABS(DI) .EQ. 0.0E0) GO TO 330
        WM(I+2) = 0.1E0*R0/DI
 320    CONTINUE
      RETURN
 330  IER = -1
      RETURN
C IF MITER = 4, CALL JAC AND MULTIPLY BY SCALAR. -----------------------
 400  ML = IWM(1)
      MU = IWM(2)
      ML3 =  3
      MBAND = ML + MU + 1
      MEBAND = MBAND + ML
      LENP = MEBAND*N
      DO 410 I = 1,LENP
 410    WM(I+2) = 0.0E0
      CALL JAC (TN, Y, WM(ML3), MEBAND, RPAR, IPAR)
      CON = -HL0
      DO 420 I = 1,LENP
 420    WM(I+2) = WM(I+2)*CON
      GO TO 570
C IF MITER = 5, MAKE MBAND CALLS TO F TO APPROXIMATE J. ----------------
 500  ML = IWM(1)
      MU = IWM(2)
      MBAND = ML + MU + 1
      MBA = MIN(MBAND,N)
      MEBAND = MBAND + ML
      MEB1 = MEBAND - 1
      SRUR = WM(1)
      FAC = VNWRMS (N, SAVF, EWT)
      R0 = 1000.0E0*ABS(H)*UROUND*N*FAC
      IF (R0 .EQ. 0.0E0) R0 = 1.0E0
      DO 560 J = 1,MBA
        DO 530 I = J,N,MBAND
          YI = Y(I)
          R = MAX(SRUR*ABS(YI),R0*EWT(I))
 530      Y(I) = Y(I) + R
        CALL F (TN, Y, FTEM, RPAR, IPAR)
        DO 550 JJ = J,N,MBAND
          Y(JJ) = YH(JJ,1)
          YJJ = Y(JJ)
          R = MAX(SRUR*ABS(YJJ),R0*EWT(JJ))
          FAC = -HL0/R
          I1 = MAX(JJ-MU,1)
          I2 = MIN(JJ+ML,N)
          II = JJ*MEB1 - ML + 2
          DO 540 I = I1,I2
 540        WM(II+I) = (FTEM(I) - SAVF(I))*FAC
 550      CONTINUE
 560    CONTINUE
      NFE = NFE + MBA
C ADD IDENTITY MATRIX. -------------------------------------------------
 570  II = MBAND + 2
      DO 580 I = 1,N
        WM(II) = WM(II) + 1.0E0
 580    II = II + MEBAND
C DO LU DECOMPOSITION OF P. --------------------------------------------
      CALL SGBFA (WM(3), MEBAND, N, ML, MU, IWM(21), IER)
      RETURN
C----------------------- END OF SUBROUTINE PJAC -----------------------
      END


text:
subroutine rdwaqpar(lundia, error, kmax, lsed, dt, itcomf, itcomi, itcoml, gdp)
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  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 version 3.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: rdwaqpar.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/flow2d3d/packages/io/src/input/rdwaqpar.f90 $
!!--description-----------------------------------------------------------------
!
! Read waq parameters from mdf file
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    use precision
    use properties
    use globaldata
    use system_utils, only: exifil
    !
    implicit none
    !
    type(globdat),target :: gdp
    !
    ! The following list of pointer parameters is used to point inside the gdp structure
    !
    integer                       , pointer :: aggre
    integer                       , pointer :: itwqff
    integer                       , pointer :: itwqfi
    integer                       , pointer :: itwqfl
    integer                       , pointer :: itstrt
    integer                       , pointer :: itfinish
    integer       , dimension(:)  , pointer :: ilaggr    ! layer aggregation pointer
    real(fp)      , dimension(:,:), pointer :: quwaq     ! Cumulative qxk
    real(fp)      , dimension(:,:), pointer :: qvwaq     ! Cumulative qyk
    real(fp)      , dimension(:,:), pointer :: qwwaq     ! Cumulative qzk
    real(fp)      , dimension(:,:), pointer :: cumsedflx ! Cumulative sedimentation flux
    real(fp)      , dimension(:,:), pointer :: cumresflx ! Cumulative resuspension flux
    real(fp)      , dimension(:)  , pointer :: discumwaq ! Cumulated sources m3/s*nstep
    logical                       , pointer :: waqfil
    logical                       , pointer :: waqol
    character(256)                , pointer :: flaggr
!
! Global variables
!
    integer   , intent(in)  :: kmax   !  Description and declaration in esm_alloc_int.f90
    integer   , intent(in)  :: lsed   !!  Number of sediment fractions
    integer                 :: lundia !  Description and declaration in inout.igs
    integer                 :: itcomf !  Description and declaration in inttim.igs
    integer                 :: itcomi !  Description and declaration in inttim.igs
    integer                 :: itcoml !  Description and declaration in inttim.igs
    logical   , intent(out) :: error  !!  Flag=TRUE if an error is encountered
    real(fp)                :: dt     !  Description and declaration in esm_alloc_real.f90
!
! Local variables
!
    integer                                             :: i
    integer                                             :: k
    integer                                             :: ilaggrIndex
    integer        ,dimension(:), allocatable           :: ilaggrInput
    integer                                             :: istat
    real(fp)                                            :: ftwqff
    real(fp)                                            :: ftwqfi
    real(fp)                                            :: ftwqfl
    real(fp)       ,dimension(3)                        :: times
    logical                                  , external :: dtn
    character(3000)                                     :: txtput2      ! usage depends on kmax
    character(30)                                       :: txtput1
    character(300)                                      :: message
!
!! executable statements -------------------------------------------------------
!
    aggre      => gdp%gdwaqpar%aggre
    itwqff     => gdp%gdwaqpar%itwqff
    itwqfi     => gdp%gdwaqpar%itwqfi
    itwqfl     => gdp%gdwaqpar%itwqfl
    itstrt     => gdp%gdinttim%itstrt
    itfinish   => gdp%gdinttim%itfinish
    ilaggr     => gdp%gdwaqpar%ilaggr
    quwaq      => gdp%gdwaqpar%quwaq
    qvwaq      => gdp%gdwaqpar%qvwaq
    qwwaq      => gdp%gdwaqpar%qwwaq
    discumwaq  => gdp%gdwaqpar%discumwaq
    waqfil     => gdp%gdwaqpar%waqfil
    waqol      => gdp%gdwaqpar%waqol
    flaggr     => gdp%gdwaqpar%flaggr
    !
    error = .false.
    !
    ! If keyword Flwq is in the mdf-file, waqfil = true
    !
    txtput1 = ' '
    call prop_get_string(gdp%mdfile_ptr, '*', 'Flwq', txtput1)
    if (txtput1 == ' ') then
       waqfil = .false.
    else
       waqfil = .true.
    endif
    !
    if (.not.waqfil) return
    !
    ! read output times
    ! Flwq  = 1.4400000e+003  6.0000000e+001  2.8800000e+003
    !
    times = 0.0_fp
    call prop_get(gdp%mdfile_ptr, '*', 'Flwq', times, 3)
    !
    ftwqff = times(1)
    ftwqfi = times(2)
    ftwqfl = times(3)
    !
    ! read vertical aggregation array ilaggr
    !
                  allocate(gdp%gdwaqpar%ilaggr     (kmax), stat=istat)
    if (istat==0) allocate(             ilaggrInput(kmax), stat=istat)
    if (istat/=0) then
       call prterr(lundia, 'U021', 'rdwaqpar: memory alloc error')
       call d3stop(1, gdp)
    endif
    !
    ! Update local pointer to point to newly allocated memory
    !
    ilaggr     => gdp%gdwaqpar%ilaggr
    !
    ! Read value of keyword WaqOL in the mdf-file
    ! When WaqOL = true, FLOW and WAQ communicate online via WAQ input files
    ! Default WaqOL = false
    !
    call prop_get_logical(gdp%mdfile_ptr, '*', 'WaqOL', waqol)
    !
    ilaggrInput = 0
    call prop_get(gdp%mdfile_ptr, '*', 'ilaggr', ilaggrInput, kmax)
    if (ilaggrInput(1) == 0) then
       ilaggrInput = 1 ! set ilaggrInput to "all ones" for echo to tri-diag file.
       do i = 1, kmax
          ilaggr(i) = i
       enddo
    else
       !
       ! Test whether the sum of ilaggrInput is equal to kmax
       !
       k = 0
       do i=1,kmax
          if (ilaggrInput(i) < 0) then
             write(message, '(a,i0,a)') 'ilaggr element ', i, &
                 & ' is negative.'
             call prterr(lundia, 'P004', trim(message))
             error = .true.
             call d3stop(1, gdp)
          endif
          k = k + ilaggrInput(i)
       enddo
       if (k /= kmax) then
          write(message, '(a,i0,a,i0,a)') 'The sum of all ilaggr elements (', k, &
              & ') is not equal to the total number of layers (',kmax, ').'
          call prterr(lundia, 'P004', trim(message))
          error = .true.
          call d3stop(1, gdp)
       endif
       !
       ! Transform ilaggrInput into ilaggr
       !
       k           = 0
       ilaggrIndex = 1
       do i=1,kmax
          k = k + 1
          do ilaggrIndex = ilaggrIndex, ilaggrIndex + ilaggrInput(i) - 1
             ilaggr(ilaggrIndex) = k
          enddo
       enddo
    endif
    !
    ! calculate integer multiples of DT and test calculated values
    !
    itwqff = nint(ftwqff/dt)
    itwqfi = nint(ftwqfi/dt)
    itwqfl = nint(ftwqfl/dt)
    !
    if (dtn(itwqff, ftwqff, dt)) then
       error = .true.
       call prterr(lundia    ,'U044'    ,'Waq store start time'          )
    endif
    if (dtn(itwqfi, ftwqfi, dt)) then
       error = .true.
       call prterr(lundia    ,'U044'    ,'Waq store time interval'       )
    endif
    if (itwqfi<=0) then
       error = .true.
       call prterr(lundia    ,'U021'    ,'Waq store time interval should be > 0' )
    endif
    if (dtn(itwqfl, ftwqfl, dt)) then
       error = .true.
       call prterr(lundia    ,'U044'    ,'Waq store stop time')
    endif
    if (mod(itwqfl-itwqff,max(itwqfi,1))/=0) then
       error = .true.
       call prterr(lundia    ,'U021'    ,'Waq store period should be multiple of interval' )
    endif
    if (itwqfl > itfinish) then
       itwqfl = itfinish
       call prterr(lundia, 'U190', 'Waq store end time is set to simulation end time')
    endif
    !
    ! read aggregation file
    ! WaqAgg= #example.dwq#
    !
    flaggr = ' '
    call prop_get_string(gdp%mdfile_ptr, '*', 'WaqAgg', flaggr)
    if (flaggr /= ' ' .and. flaggr /= 'active only') then
       if (.not. exifil(flaggr, lundia)) then
          error = .true.
          call d3stop(1, gdp)
       endif
    endif
    !
    ! echo input in diagnose-file
    !
    write (lundia, '(a)') '*** Start  of waq parameters input'
    !
    txtput1 = 'Write waq output'
    txtput2 = '       Yes'
    write (lundia, '(3a)') txtput1, ':  ', trim(txtput2)
    !
    txtput1 = 'Start at time step'
    write (lundia, '(2a,i12)') txtput1, ':', itwqff
    txtput1 = 'Output every N time steps'
    write (lundia, '(2a,i12)') txtput1, ':', itwqfi
    txtput1 = 'Stop at time step'
    write (lundia, '(2a,i12)') txtput1, ':', itwqfl
    txtput1 = 'Aggregation (horizontal)'
    txtput2 = ' '
    message = ' '
    if (flaggr == ' ') then
       write (message, '(a)') ':        none'
    else
       !
       ! flaggr == 'active only' or flaggr == [filename]
       !
       write (message, '(a)') trim(flaggr)
    endif
    txtput2 = trim(message)
    write (lundia, '(3a)') txtput1, ': ', trim(txtput2)
    txtput1 = 'Aggregation (vertical)'
    txtput2 = ' '
    do i=1,kmax
       if (ilaggr(i) /= 0) then
          message = ' '
          if (i == 1) then
             write(message,'(i0)') ilaggrInput(i)
          elseif (ilaggrInput(i)>0) then ! don't show zeroes at end
             write(message,'(a,i0)') ', ', ilaggrInput(i)
          endif
          txtput2 = trim(txtput2)//trim(message)
       endif
    enddo
    write (lundia, '(3a)') txtput1, ': ', trim(txtput2)
    write (lundia, '(a)') '*** End    of waq parameters input'
    !
    ! When no communication file is written:
    ! Force the communication file to be written for one time step
    ! This is needed by WAQ
    !
    if (itcomi == 0) then
       itcomf = itstrt
       itcomi = 1
       itcoml = itstrt + 1
       call prterr(lundia, 'G051', 'The communication file is written during the first time step only. Needed by WAQ.')
    endif
    write (lundia, *)
    !
    ! Writing binary waq file enabled: allocate the required memory
    !
    istat = 0
    if (.not. associated(gdp%gdwaqpar%quwaq)) then
       !
       ! allocate all arrays needed for writing binary waq files
       !
       if (istat==0) allocate (gdp%gdwaqpar%quwaq    (gdp%d%nmlb:gdp%d%nmub,  kmax)     , stat = istat)
       if (istat==0) allocate (gdp%gdwaqpar%qvwaq    (gdp%d%nmlb:gdp%d%nmub,  kmax)     , stat = istat)
       if (istat==0) allocate (gdp%gdwaqpar%qwwaq    (gdp%d%nmlb:gdp%d%nmub,0:kmax)     , stat = istat)
       if (istat==0) allocate (gdp%gdwaqpar%discumwaq(gdp%d%nsrc)                       , stat = istat)
       if (istat==0) allocate (gdp%gdwaqpar%cumsedflx(gdp%d%nmlb:gdp%d%nmub,  lsed)     , stat = istat)
       if (istat==0) allocate (gdp%gdwaqpar%cumresflx(gdp%d%nmlb:gdp%d%nmub,  lsed)     , stat = istat)
       !
       if (istat /= 0) then
          call prterr(lundia, 'U021', 'RDWAQPAR: memory alloc error')
          call d3stop(1, gdp)
       endif
       !
       ! Update local pointers to point to newly allocated memory
       !
       quwaq      => gdp%gdwaqpar%quwaq
       qvwaq      => gdp%gdwaqpar%qvwaq
       qwwaq      => gdp%gdwaqpar%qwwaq
       discumwaq  => gdp%gdwaqpar%discumwaq
       cumsedflx  => gdp%gdwaqpar%cumsedflx
       cumresflx  => gdp%gdwaqpar%cumresflx
    endif
    !
    ! Initialize arrays allocated
    !
    quwaq     = 0.0_fp
    qvwaq     = 0.0_fp
    qwwaq     = 0.0_fp
    discumwaq = 0.0_fp
    cumsedflx = 0.0_fp
    cumresflx = 0.0_fp
    deallocate(ilaggrInput)
end subroutine rdwaqpar

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
c
c     version 1  Oct 1995
c
c     ----------------------------------------------------------------
c     |                                                              |
c     |           *******  S M O O T H L S Q R  *******              |
c     |                                                              |
c     |                 Least-squares solution of                    |
c     |                traveltime inversion problem                  |
c     |               using smoothing regularisation                 |
c     |                                                              |
c     |                   Written by C. A. Zelt                      |
c     |                 Modified by N. J. Bruguier                   |
c     |                                                              |
c     |                   Bullard Laboratories                       |
c     |                  Cambridge, CB3 0EZ, UK                      |
c     |                                                              |
c     ----------------------------------------------------------------
c
c
c     I/O units:
c
c        10 -- input:  damping parameters
c
c        11 -- input:  matrix of partial derivatives and vector of
c                      traveltime residuals
c
c        12 -- output: model parameter adjustments and resolution and
c                      covariance estimates and updated model
c
c        13 -- input/output: initial and updated velocity model
c
c        30 -- input:  initial and updated floating reflectors
c
c
c     ----------------------------------------------------------------
c
c
      program main
c
      include 'rayinvr.par'
c
      real apart(prayi,pnvar),tres(prayi),ata(pnvar,pnvar),
     +     att(pnvar),atad(pnvar,pnvar),atadi(pnvar,pnvar),
     +     dx(pnvar),res(pnvar,pnvar),var(pnvar),parorg(pnvar),
     +     xm(pncntr,ppcntr),zm(pncntr,ppcntr),grad(player*ppvel),
     +     xvel(player,ppvel,2),vf(player,ppvel,2),
     +     tunc(prayi),thick(player*ppcntr),parunc(pnvar),
     +     xfrefl(pfrefl,ppfref),zfrefl(pfrefl,ppfref),
     +     d(pnvar,pnvar),dtd(pnvar,pnvar),dtdm(pnvar),
     +     parbsf(pnvar),velw(9),bndw(9)
      integer partyp(pnvar),nzed(pncntr),nvel(player,2),
     +        ivarz(player,ppcntr),ivarv(player,ppvel,2),
     +        mpinch(pncntr,ppcntr,2),igrad(player*ppvel),
     +        ivarf(pfrefl,ppfref),npfref(pfrefl),parw(pnvar),
     +        iparz(player,ppcntr),iparv(player,ppvel,2)
c
      namelist /dmppar/ iscrn,ilayer,dmpfct,velunc,bndunc,xmax,
     +                  velw,bndw
c
      open(unit=10, file='d.in', status='old')
      open(unit=11, file='i.out', status='old')
      open(unit=12, file='d.out')
      open(unit=13, file='v.in', status='old')
      open(unit=14, file='v.bak')
c
c default parameter values
c
      iscrn=1
      ilayer=0
      dmpfct=1.0
      velunc=-1.
      bndunc=-1.
      xmax=-99999.
      ifrbnd=0
c
      read(10,dmppar)
c
      if(xmax.lt.-99998) then
        write(6,105)
105     format(/'***  xmax not specified  ***'/)
        stop
      end if
c
c read in matrix of partial derivatives
c and vector of traveltime residuals
c
      read(11,1)
1     format(' ')
      read(11,115) narinv,nvar
115   format(2i5)
      read(11,1)
      do 6 i=1,nvar
         read(11,5) partyp(i),parorg(i),parw(i)
5        format(i5,f15.5,i10)
6     continue
      read(11,1)
      do 10 i=1,narinv
         read(11,15) (apart(i,j),j=1,nvar)
15       format(5e12.5)
10    continue
      read(11,1)
      read(11,15) (tres(i),i=1,narinv)
      read(11,1)
      read(11,15) (tunc(i),i=1,narinv)
c
c parameter uncertainties
c
      do 70 i=1,nvar
         if(partyp(i).eq.1) then
            if(bndunc.le.0.) then
               parunc(i)=bndw(parw(i))**2
            else
               parunc(i)=bndunc**2
            end if
         end if
         if(partyp(i).eq.2) then
            if(velunc.le.0.) then
               parunc(i)=velw(parw(i))**2
            else
               parunc(i)=velunc**2
            end if
         end if
         if(partyp(i).eq.3) then
            if(bndunc.le.0.) then
               parunc(i)=bndw(parw(i))**2
            else
               parunc(i)=bndunc**2
            end if
            ifrbnd=1
         end if
70    continue
c
c read in velocity model
c
      ncont=1
      nrzmax=ppcntr/10
      nrvmax=ppvel/10
      do 170 icont=1,player+1
         nrz=1
         j1=1
         j2=10
11       if(nrz.gt.nrzmax) go to 211
         read(13,55,end=999) ilyr,(xm(icont,j),j=j1,j2)
         read(13,55,end=999) icnt,(zm(icont,j),j=j1,j2)
         read(13,65,end=99) (ivarz(icont,j),j=j1,j2)
55       format(i2,1x,10f7.2)
65       format(3x,10i7)
         nrz=nrz+1
         if(icnt.ne.1) go to 211
         j1=j1+10
         j2=j2+10
         go to 11
211      nrv=1
         j1=1
         j2=10
21       if(nrv.gt.nrvmax) go to 311
         read(13,55,end=999) ilyr,(xvel(icont,j,1),j=j1,j2)
         read(13,55,end=999) icnt,(vf(icont,j,1),j=j1,j2)
         read(13,65,end=999) (ivarv(icont,j,1),j=j1,j2)
         nrv=nrv+1
         if(icnt.ne.1) go to 311
         j1=j1+10
         j2=j2+10
         go to 21
311      nrv=1
         j1=1
         j2=10
31       if(nrv.gt.nrvmax) go to 411
         read(13,55,end=999) ilyr,(xvel(icont,j,2),j=j1,j2)
         read(13,55,end=999) icnt,(vf(icont,j,2),j=j1,j2)
         read(13,65,end=999) (ivarv(icont,j,2),j=j1,j2)
         nrv=nrv+1
         if(icnt.ne.1) go to 411
         j1=j1+10
         j2=j2+10
         go to 31
411      ncont=ncont+1
170   continue
c
99    nlayer=ncont-1
c
c calculate number of parameters per layer
c
      do 171 i=1,ncont
         nzed(i)=1
171   continue
      do 172 i=1,nlayer
         nvel(i,1)=1
         nvel(i,2)=1
172   continue
c
      do 180 i=1,ncont
         do 190 j=1,ppcntr
            if(abs(xm(i,j)-xmax).lt..0001) go to 180
            nzed(i)=nzed(i)+1
190      continue
180   continue
c
      do 210 i=1,nlayer
         do 220 j=1,ppvel
            if(abs(xvel(i,j,1)-xmax).lt..0001) go to 212
            nvel(i,1)=nvel(i,1)+1
220      continue
212      if(nvel(i,1).eq.1.and.vf(i,1,1).eq.0.) nvel(i,1)=0
210   continue
c
      do 240 i=1,nlayer
         do 250 j=1,ppvel
            if(abs(xvel(i,j,2)-xmax).lt..0001) go to 260
            nvel(i,2)=nvel(i,2)+1
250      continue
260      if(nvel(i,2).eq.1.and.vf(i,1,2).eq.0.) nvel(i,2)=0
240   continue
c
c read in floating reflectors
c
      if(ifrbnd.eq.1) then
        open(unit=30, file='f.in', status='old')
        open(unit=31, file='f.bak')
        nfrefl=0
690     read(30,545,end=595) nfrefr
        if(nfrefr.lt.2.or.nfrefr.gt.ppfref) then
          write(6,585)
585       format(/'***  error in f.in file  ***'/)
          stop
        end if
        nfrefl=nfrefl+1
        npfref(nfrefl)=nfrefr
545     format(i2)
        read(30,555) (xfrefl(nfrefl,i),i=1,npfref(nfrefl))
        read(30,555) (zfrefl(nfrefl,i),i=1,npfref(nfrefl))
        read(30,575) (ivarf(nfrefl,i),i=1,npfref(nfrefl))
555     format(3x,<ppfref>f7.2)
575     format(3x,<ppfref>i7)
        go to 690
595     continue
      end if
c
c change depths to layer thicknesses if required
c
      np=0
      if (ilayer.eq.1) then
        do 1281 i=2,nlayer
          do 1282 j=1,nzed(i)
            if(ivarz(i,j).gt.0) then
              np=np+1
              do 1284 k=1,nzed(i-1)
                if(abs(xm(i,j)-xm(i-1,k)).lt..005) then
                  parbsf(np)=zm(i,j)-zm(i-1,k)
                  goto 1282
                end if
1284          continue
            endif
1282      continue
c
          if (nvel(i,1).gt.0) then
            do 1285 j=1,nvel(i,1)
              if(ivarv(i,j,1).gt.0) then
                np=np+1
                parbsf(np)=parorg(np)
              end if
1285        continue
          end if
c
          if (nvel(i,2).gt.0) then
            do 1286 j=1,nvel(i,2)
              if(ivarv(i,j,2).gt.0) then
                np=np+1
                parbsf(np)=parorg(np)
              end if
1286        continue
          end if
1281   continue
c
       if(ifrbnd.eq.1) then
         do 1287 i=1,nfrefl
           do 1288 j=1,npfref(j)
             if(ivarf(i,j).eq.1) then
               np=np+1
               parbsf(np)=parorg(np)
             end if
1288       continue
1287     continue
      end if
c
      if(np.ne.nvar) goto 999
      endif
c
c calculate difference matrix D
c
      do 1250 i=1,nvar
         do 1260 j=1,nvar
            d(i,j)=0.
1260     continue
1250  continue
c
      np=0
c
      do 1270 i=1,nlayer
         do 1280 j=1,nzed(i)
            if(ivarz(i,j).gt.0) then
              np=np+1
              if(j.gt.1) then
                dist=xm(i,j)-xm(i,j-1)
                d(np,np)=1./dist
                d(np,np-1)=-1./dist
              endif
            end if
1280     continue
c
         if (nvel(i,1).gt.0) then
           do 1310 j=1,nvel(i,1)
             if(ivarv(i,j,1).gt.0) then
               np=np+1
               if(j.gt.1) then
                 dist=xvel(i,j,1)-xvel(i,j-1,1)
                 d(np,np)=1./dist
                 d(np,np-1)=-1./dist
               end if
             end if
1310       continue
         end if
c
         if (nvel(i,2).gt.0) then
           do 1430 j=1,nvel(i,2)
             if(ivarv(i,j,2).gt.0) then
               np=np+1
               if(j.gt.1) then
                 dist=xvel(i,j,2)-xvel(i,j-1,2)
                 d(np,np)=1./dist
                 d(np,np-1)=-1./dist
               end if
             end if
1430       continue
         end if
1270  continue
c
      if(ifrbnd.eq.1) then
        do 1880 i=1,nfrefl
          do 1890 j=1,npfref(j)
            if(ivarf(i,j).eq.1) then
              np=np+1
              if(j.gt.1) then
                dist=xfrefl(i,j)-xfrefl(i,j-1)
                d(np,np)=1./dist
                d(np,np-1)=-1./dist
              end if
            end if
1890      continue
1880    continue
      end if
c
      if(np.ne.nvar) goto 999
c
c calculate ATA
c
      do 20 i=1,nvar
         do 30 j=1,i
            ata(i,j)=0.
            do 40 k=1,narinv
               ata(i,j)=ata(i,j)+apart(k,i)*apart(k,j)/tunc(k)**2
40          continue
            if(i.ne.j) ata(j,i)=ata(i,j)
30       continue
20    continue
c
c calculate ATT
c
      do 50 i=1,nvar
         att(i)=0.
         do 60 j=1,narinv
            att(i)=att(i)+apart(j,i)*tres(j)/tunc(j)**2
60       continue
50    continue
c
c calculate DTD and (ATA+DTD)
c
      do 1520 i=1,nvar
         do 1530 j=1,nvar
           dtd(i,j)=0.
           do 1540 k=1,nvar
             dtd(i,j)=dtd(i,j)+d(k,i)*d(k,j)*dmpfct/parunc(k)
1540       continue
	   atad(i,j)=ata(i,j)+dtd(i,j)
1530     continue
1520  continue
c
c calculate DTDm
c
      if (ilayer.eq.1) then
        do 1550 i=1,nvar
          dtdm(i)=0.
          do 1560 j=1,nvar
             dtdm(i)=dtdm(i)+dtd(i,j)*parbsf(j)
1560      continue
1550    continue
      else
        do 1551 i=1,nvar
          dtdm(i)=0.
          do 1561 j=1,nvar
             dtdm(i)=dtdm(i)+dtd(i,j)*parorg(j)
1561      continue
1551    continue
      endif
c
c invert matrix (ATA+DTD)
c
      call matinv(atad,atadi,nvar)
c
c calculate parameter adjustments
c
      do 90 i=1,nvar
         dx(i)=0.
         do 100 j=1,nvar
            dx(i)=dx(i)+atadi(i,j)*(att(j)-dtdm(j))
100      continue
90    continue
c
c calculate resolution matrix
c
      do 110 i=1,nvar
         do 120 j=1,nvar
            res(i,j)=0.
            do 130 k=1,nvar
               res(i,j)=res(i,j)+atadi(i,k)*ata(k,j)
130         continue
120      continue
110   continue
c
c calculate variance
c
      do 140 i=1,nvar
         var(i)=var(i)+(1.0-res(i,i))*parunc(i)
140   continue
c
c output results
c
      write(12,25) dmpfct
25    format(/'overall damping factor: ',f10.2)
      write(12,35)
35    format(/'type  orig. val.  uncert.    adjust.   new val.',
     +        ' resolution std. error')
      write(12,45) (partyp(i),parorg(i),parunc(i)**.5,dx(i),
     +              parorg(i)+dx(i),res(i,i),sqrt(var(i)),i=1,nvar)
45    format(i3,6f11.4)
c
      if(iscrn.eq.1) then
        write(6,25) dmpfct
        write(6,35)
        write(6,45) (partyp(i),parorg(i),parunc(i)**.5,dx(i),parorg(i)+
     +               dx(i),res(i,i),sqrt(var(i)),i=1,nvar)
      end if
c
c write out original velocity model
c
      do 1570 i=1,nlayer
         nstart=1
1590     j1=nstart
         j2=j1+9
         if(j2.gt.nzed(i)) j2=nzed(i)
         if(j2.lt.nzed(i)) then
           icnt=1
         else
           icnt=0
         end if
         write(14,55) i,(xm(i,j),j=j1,j2)
         write(14,55) icnt,(zm(i,j),j=j1,j2)
         write(14,65) (ivarz(i,j),j=j1,j2)
         if(j2.eq.nzed(i)) go to 1600
         nstart=j2+1
         go to 1590
1600     if(nvel(i,1).le.0) then
           write(14,55) i,xmax
           write(14,55) 0,0.
           write(14,65) 0
           go to 1630
         end if
         nstart=1
1620     j1=nstart
         j2=j1+9
         if(j2.gt.nvel(i,1)) j2=nvel(i,1)
         if(j2.lt.nvel(i,1)) then
           icnt=1
         else
           icnt=0
         end if
         write(14,55) i,(xvel(i,j,1),j=j1,j2)
         write(14,55) icnt,(vf(i,j,1),j=j1,j2)
         write(14,65) (ivarv(i,j,1),j=j1,j2)
         if(j2.eq.nvel(i,1)) go to 1630
         nstart=j2+1
         go to 1620
1630     if(nvel(i,2).le.0) then
           write(14,55) i,xmax
           write(14,55) 0,0.
           write(14,65) 0
           go to 1570
         end if
         nstart=1
1650     j1=nstart
         j2=j1+9
         if(j2.gt.nvel(i,2)) j2=nvel(i,2)
         if(j2.lt.nvel(i,2)) then
           icnt=1
         else
           icnt=0
         end if
         write(14,55) i,(xvel(i,j,2),j=j1,j2)
         write(14,55) icnt,(vf(i,j,2),j=j1,j2)
         write(14,65) (ivarv(i,j,2),j=j1,j2)
         if(j2.eq.nvel(i,2)) go to 1570
         nstart=j2+1
         go to 1650
1570  continue
c
      write(14,55) ncont,(xm(ncont,j),j=1,nzed(ncont))
      write(14,55) 0,(zm(ncont,j),j=1,nzed(ncont))
c
c write out original floating reflectors
c
      if(ifrbnd.eq.1) then
        do 1790 j=1,nfrefl
           write(31,545) npfref(j)
           write(31,565) j,(xfrefl(j,i),i=1,npfref(j))
           write(31,555) (zfrefl(j,i),i=1,npfref(j))
           write(31,575) (ivarf(j,i),i=1,npfref(j))
1790    continue
      end if
c
c calculate layer thickness if held constant
c
      nlyr=0
      do 470 i=2,ncont
         do 480 j=1,nzed(i)
            if(ivarz(i,j).eq.-1) then
              nlyr=nlyr+1
              thick(nlyr)=zm(i,j)-zm(i-1,j)
            end if
480      continue
470   continue
c
c check for fixed velocity gradients
c
      ngrad=0
      do 370 i=1,nlayer
         iflagg=0
         if(nvel(i,2).gt.0) then
           do 380 j=1,nvel(i,2)
              if(ivarv(i,j,2).eq.-1) then
                ngrad=ngrad+1
                xbndc=xvel(i,j,2)
                if(nzed(i).gt.1) then
                  do 390 k=1,nzed(i)-1
                     if(xbndc.ge.xm(i,k).and.xbndc.le.xm(i,k+1))
     +               then
                       zu=(zm(i,k+1)-zm(i,k))*(xbndc-xm(i,k))/
     +                    (xm(i,k+1)-xm(i,k))+zm(i,k)
                       go to 400
                     end if
390               continue
                else
                  zu=zm(i,1)
                end if
400             if(nvel(i,1).gt.0) then
                  vu=vf(i,j,1)
                else
                  do 401 k=i-1,1,-1
                     if(nvel(k,2).gt.0) then
                       vu=vf(k,j,2)
                       go to 402
                     end if
                     if(nvel(k,1).gt.0) then
                       vu=vf(k,j,1)
                       go to 402
                     end if
401               continue
                end if
402             if(nzed(i+1).gt.1) then
                  do 410 k=1,nzed(i+1)-1
                     if(xbndc.ge.xm(i+1,k).and.xbndc.le.xm(i+1,k+1))
     +               then
                       zl=(zm(i+1,k+1)-zm(i+1,k))*(xbndc-xm(i+1,k))/
     +                    (xm(i+1,k+1)-xm(i+1,k))+zm(i+1,k)
                       go to 420
                     end if
410               continue
                else
                  zl=zm(i+1,1)
                end if
420             vl=vf(i,j,2)
                if(abs(vu-vl).gt..001) then
                  if(abs(zl-zu).gt..001) then
                    grad(ngrad)=(vl-vu)/(zl-zu)
                    igrad(ngrad)=1
                  else
                    grad(ngrad)=vl-vu
                    igrad(ngrad)=0
                    iflagg=1
                  end if
                else
                  grad(ngrad)=0.
                  igrad(ngrad)=1
                end if
              end if
380        continue
         end if
c
         if(iflagg.eq.1) then
           write(6,85) i
           write(12,85) i
85         format(/'***  check gradient in layer ',i2,'  ***'/)
         end if
c
370   continue
c
      rewind(13)
c
      write(12,75)
      if(iscrn.eq.1) write(6,75)
75    format(/'velocity model:'/)
c
c check for pinchouts
c
      if(nlayer.gt.1) then
        ipinch=0
        do 281 i=2,nlayer
           do 282 j=1,nzed(i)
              mpinch(i,j,1)=0
              do 283 k=1,i-1
                 do 284 l=1,nzed(k)
                    if(abs(xm(i,j)-xm(k,l)).lt..005.and.
     +                 abs(zm(i,j)-zm(k,l)).lt..005) then
                       ipinch=+1
                       mpinch(i,j,1)=k
                       mpinch(i,j,2)=l
                       go to 282
                    end if
284              continue
283           continue
282        continue
281     continue
      end if
c
c add the parameter adjustments
c
      nvarw=0
      nlyr=0
      do 270 i=1,nlayer
         do 280 j=1,nzed(i)
            if(ipinch.eq.1.and.i.gt.1) then
              if(mpinch(i,j,1).gt.0) then
                zm(i,j)=zm(mpinch(i,j,1),mpinch(i,j,2))
                go to 280
              end if
            end if
            if(ivarz(i,j).gt.0) then
              nvarw=nvarw+1
              zm(i,j)=zm(i,j)+dx(nvarw)
            end if
            if(ivarz(i,j).eq.-1) then
              nlyr=nlyr+1
              zm(i,j)=zm(i-1,j)+thick(nlyr)
            end if
280      continue
c
         if(nvel(i,1).gt.0) then
           do 310 j=1,nvel(i,1)
              if(ivarv(i,j,1).gt.0) then
                nvarw=nvarw+1
                vf(i,j,1)=vf(i,j,1)+dx(nvarw)
              end if
310        continue
         else
           nvel(i,1)=1
           xvel(i,1,1)=xmax
           vf(i,1,1)=0.
         end if
c
         if(nvel(i,2).gt.0) then
           do 430 j=1,nvel(i,2)
              if(ivarv(i,j,2).gt.0) then
                nvarw=nvarw+1
                vf(i,j,2)=vf(i,j,2)+dx(nvarw)
              end if
430        continue
         else
           nvel(i,2)=1
           xvel(i,1,2)=xmax
           vf(i,1,2)=0.
         end if
270   continue
c
      if(ifrbnd.eq.1) then
        do 880 j=1,nfrefl
           do 890 i=1,npfref(j)
              if(ivarf(j,i).eq.1) then
                nvarw=nvarw+1
                zfrefl(j,i)=zfrefl(j,i)+dx(nvarw)
              end if
890        continue
880     continue
      end if
c
c maintain fixed velocity gradients
c
      ngrad=0
      do 271 i=1,nlayer
         do 272 j=1,nvel(i,2)
            if(ivarv(i,j,2).eq.-1) then
              ngrad=ngrad+1
              xbndc=xvel(i,j,2)
              if(nzed(i).gt.1) then
                do 440 k=1,nzed(i)-1
                   if(xbndc.ge.xm(i,k).and.xbndc.le.xm(i,k+1))
     +             then
                     zu=(zm(i,k+1)-zm(i,k))*(xbndc-xm(i,k))/
     +                  (xm(i,k+1)-xm(i,k))+zm(i,k)
                     go to 450
                   end if
440             continue
              else
                zu=zm(i,1)
              end if
450           if(nvel(i,1).gt.0.and.vf(i,1,1).gt.0.) then
                vu=vf(i,j,1)
              else
                do 451 k=i-1,1,-1
                   if(nvel(k,2).gt.0.and.vf(k,1,2).gt.0.) then
                     vu=vf(k,j,2)
                     go to 452
                   end if
                   if(nvel(k,1).gt.0.and.vf(k,1,1).gt.0.) then
                     vu=vf(k,j,1)
                     go to 452
                   end if
451             continue
              end if
              if(igrad(ngrad).eq.1) then
452             if(nzed(i+1).gt.1) then
                  do 460 k=1,nzed(i+1)-1
                     if(xbndc.ge.xm(i+1,k).and.xbndc.le.xm(i+1,k+1))
     +               then
                       zl=(zm(i+1,k+1)-zm(i+1,k))*(xbndc-xm(i+1,k))/
     +                    (xm(i+1,k+1)-xm(i+1,k))+zm(i+1,k)
                       vf(i,j,2)=vu+grad(ngrad)*(zl-zu)
                       go to 272
                     end if
460               continue
                else
                  zl=zm(i+1,1)
                  vf(i,j,2)=vu+grad(ngrad)*(zl-zu)
                end if
              else
                vf(i,j,2)=vu+grad(ngrad)
              end if
            end if
272      continue
271   continue
c
c write out velocity model
c
      do 570 i=1,nlayer
         nstart=1
590      j1=nstart
         j2=j1+9
         if(j2.gt.nzed(i)) j2=nzed(i)
         if(j2.lt.nzed(i)) then
           icnt=1
         else
           icnt=0
         end if
         write(12,55) i,(xm(i,j),j=j1,j2)
         write(12,55) icnt,(zm(i,j),j=j1,j2)
         write(12,65) (ivarz(i,j),j=j1,j2)
         write(13,55) i,(xm(i,j),j=j1,j2)
         write(13,55) icnt,(zm(i,j),j=j1,j2)
         write(13,65) (ivarz(i,j),j=j1,j2)
         if(iscrn.eq.1) then
           write(6,55) i,(xm(i,j),j=j1,j2)
           write(6,55) icnt,(zm(i,j),j=j1,j2)
           write(6,65) (ivarz(i,j),j=j1,j2)
         end if
         if(j2.eq.nzed(i)) go to 600
         nstart=j2+1
         go to 590
600      nstart=1
620      j1=nstart
         j2=j1+9
         if(j2.gt.nvel(i,1)) j2=nvel(i,1)
         if(j2.lt.nvel(i,1)) then
           icnt=1
         else
           icnt=0
         end if
         write(12,55) i,(xvel(i,j,1),j=j1,j2)
         write(12,55) icnt,(vf(i,j,1),j=j1,j2)
         write(12,65) (ivarv(i,j,1),j=j1,j2)
         write(13,55) i,(xvel(i,j,1),j=j1,j2)
         write(13,55) icnt,(vf(i,j,1),j=j1,j2)
         write(13,65) (ivarv(i,j,1),j=j1,j2)
         if(iscrn.eq.1) then
           write(6,55) i,(xvel(i,j,1),j=j1,j2)
           write(6,55) icnt,(vf(i,j,1),j=j1,j2)
           write(6,65) (ivarv(i,j,1),j=j1,j2)
         end if
         if(j2.eq.nvel(i,1)) go to 630
         nstart=j2+1
         go to 620
630      nstart=1
650      j1=nstart
         j2=j1+9
         if(j2.gt.nvel(i,2)) j2=nvel(i,2)
         if(j2.lt.nvel(i,2)) then
           icnt=1
         else
           icnt=0
         end if
         write(12,55) i,(xvel(i,j,2),j=j1,j2)
         write(12,55) icnt,(vf(i,j,2),j=j1,j2)
         write(12,65) (ivarv(i,j,2),j=j1,j2)
         write(13,55) i,(xvel(i,j,2),j=j1,j2)
         write(13,55) icnt,(vf(i,j,2),j=j1,j2)
         write(13,65) (ivarv(i,j,2),j=j1,j2)
         if(iscrn.eq.1) then
           write(6,55) i,(xvel(i,j,2),j=j1,j2)
           write(6,55) icnt,(vf(i,j,2),j=j1,j2)
           write(6,65) (ivarv(i,j,2),j=j1,j2)
         end if
         if(j2.eq.nvel(i,2)) go to 570
         nstart=j2+1
         go to 650
570   continue
c
      write(12,55) ncont,(xm(ncont,j),j=1,nzed(ncont))
      write(12,55) 0,(zm(ncont,j),j=1,nzed(ncont))
      write(13,55) ncont,(xm(ncont,j),j=1,nzed(ncont))
      write(13,55) 0,(zm(ncont,j),j=1,nzed(ncont))
      if(iscrn.eq.1) then
        write(6,55) ncont,(xm(ncont,j),j=1,nzed(ncont))
        write(6,55) 0,(zm(ncont,j),j=1,nzed(ncont))
      end if
c
c write out floating reflectors
c
      if(ifrbnd.eq.1) then
        write(12,175)
        if(iscrn.eq.1) write(6,175)
175     format(/'floating reflectors:'/)
        rewind(30)
        do 790 j=1,nfrefl
           write(30,545) npfref(j)
           write(12,545) npfref(j)
           write(30,565) j,(xfrefl(j,i),i=1,npfref(j))
           write(12,565) j,(xfrefl(j,i),i=1,npfref(j))
565        format(i2,1x,<ppfref>f7.2)
           write(30,555) (zfrefl(j,i),i=1,npfref(j))
           write(12,555) (zfrefl(j,i),i=1,npfref(j))
           write(30,575) (ivarf(j,i),i=1,npfref(j))
           write(12,575) (ivarf(j,i),i=1,npfref(j))
           if(iscrn.eq.1) then
             write(6,545) npfref(j)
             write(6,565) j,(xfrefl(j,i),i=1,npfref(j))
             write(6,555) (zfrefl(j,i),i=1,npfref(j))
             write(6,575) (ivarf(j,i),i=1,npfref(j))
           end if
790     continue
      end if
c
      stop
c
999   write(6,95)
95    format(/'***  error in velocity model  ***'/)
      stop
c
      end
c
c     ----------------------------------------------------------------
c
      subroutine matinv(a,y,n)
c
c     invert the nxn matrix a
c
      include 'rayinvr.par'
c
      real a(pnvar,pnvar),y(pnvar,pnvar)
      integer indx(pnvar)
c
      do 10 i=1,n
         do 20 j=1,n
            y(i,j)=0.
20       continue
         y(i,i)=1.
10    continue
c
      call ludcmp(a,n,indx,d)
c
      do 30 j=1,n
         call lubksb(a,n,indx,y(1,j))
30    continue
c
      return
      end
c
c     ----------------------------------------------------------------
c
      subroutine ludcmp(a,n,indx,d)
c
c     replace a by its LU decomposition
c
      include 'rayinvr.par'
c
      real a(pnvar,pnvar),vv(pnvar)
      integer indx(n)
c
      tiny=1.0e-20
c
      d=1.
      do 10 i=1,n
         aamax=0.
         do 20 j=1,n
            if(abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))
20       continue
         if(aamax.eq.0.) then
           write(6,5)
5          format(/'***  singular matrix  ***'/)
           stop
         end if
         vv(i)=1./aamax
10    continue
      do 30 j=1,n
         if(j.gt.1) then
           do 40 i=1,j-1
              sum=a(i,j)
              if(i.gt.1) then
                do 50 k=1,i-1
                   sum=sum-a(i,k)*a(k,j)
50              continue
                a(i,j)=sum
              end if
40         continue
         end if
         aamax=0.
         do 60 i=j,n
            sum=a(i,j)
            if(j.gt.1) then
              do 70 k=1,j-1
                 sum=sum-a(i,k)*a(k,j)
70            continue
              a(i,j)=sum
            end if
            dum=vv(i)*abs(sum)
            if(dum.ge.aamax) then
              imax=i
              aamax=dum
            end if
60       continue
         if(j.ne.imax) then
           do 80 k=1,n
              dum=a(imax,k)
              a(imax,k)=a(j,k)
              a(j,k)=dum
80         continue
           d=-d
           vv(imax)=vv(j)
         end if
         indx(j)=imax
         if(j.ne.n) then
           if(a(j,j).eq.0.) a(j,j)=tiny
           dum=1./a(j,j)
           do 90 i=j+1,n
              a(i,j)=a(i,j)*dum
90         continue
         end if
30    continue
      if(a(n,n).eq.0.) a(n,n)=tiny
      return
      end
c
c     ----------------------------------------------------------------
c
      subroutine lubksb(a,n,indx,b)
c
c     solve the system of n linear equations ax=b
c
      include 'rayinvr.par'
c
      real a(pnvar,pnvar),b(n)
      integer indx(n)
c
      ii=0
      do 10 i=1,n
         ll=indx(i)
         sum=b(ll)
         b(ll)=b(i)
         if(ii.ne.0) then
           do 20 j=ii,i-1
              sum=sum-a(i,j)*b(j)
20         continue
         else if(sum.ne.0.) then
           ii=i
         end if
         b(i)=sum
10    continue
      do 30 i=n,1,-1
         sum=b(i)
         if(i.lt.n) then
           do 40 j=i+1,n
              sum=sum-a(i,j)*b(j)
40         continue
         end if
         b(i)=sum/a(i,i)
30    continue
      return
      end


text:
subroutine tem_solite(a,itask)
   use typre
   use Mod_Temperature
   use def_parame
   implicit none
   class(TemperatureProblem) :: a
   integer(ip) :: itask

   !Nothing to do here

end subroutine

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


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

20090525 16:55:30 nbsp Welcome to the Wiki. Allow me to introduce you to the Bike Recovery System. Please move your postings about your missing bike there, and good luck getting it back. Users/JasonAller

20100801 00:10:23 nbsp Geez... you had another bike stolen?!  Ill keep an eye out. Users/Aaron.Curtin

20100801 08:25:10 nbsp Thats a shame about your bike, but you should post lost bikes only on the Bike Recovery System. Users/BrianPakpour

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
*DECK DPCHFE
      SUBROUTINE DPCHFE (N, X, F, D, INCFD, SKIP, NE, XE, FE, IERR)
C***BEGIN PROLOGUE  DPCHFE
C***PURPOSE  Evaluate a piecewise cubic Hermite function at an array of
C            points.  May be used by itself for Hermite interpolation,
C            or as an evaluator for DPCHIM or DPCHIC.
C***LIBRARY   SLATEC (PCHIP)
C***CATEGORY  E3
C***TYPE      DOUBLE PRECISION (PCHFE-S, DPCHFE-D)
C***KEYWORDS  CUBIC HERMITE EVALUATION, HERMITE INTERPOLATION, PCHIP,
C             PIECEWISE CUBIC EVALUATION
C***AUTHOR  Fritsch, F. N., (LLNL)
C             Lawrence Livermore National Laboratory
C             P.O. Box 808  (L-316)
C             Livermore, CA  94550
C             FTS 532-4275, (510) 422-4275
C***DESCRIPTION
C
C          DPCHFE:  Piecewise Cubic Hermite Function Evaluator
C
C     Evaluates the cubic Hermite function defined by  N, X, F, D  at
C     the points  XE(J), J=1(1)NE.
C
C     To provide compatibility with DPCHIM and DPCHIC, includes an
C     increment between successive values of the F- and D-arrays.
C
C ----------------------------------------------------------------------
C
C  Calling sequence:
C
C        PARAMETER  (INCFD = ...)
C        INTEGER  N, NE, IERR
C        DOUBLE PRECISION  X(N), F(INCFD,N), D(INCFD,N), XE(NE), FE(NE)
C        LOGICAL  SKIP
C
C        CALL  DPCHFE (N, X, F, D, INCFD, SKIP, NE, XE, FE, IERR)
C
C   Parameters:
C
C     N -- (input) number of data points.  (Error return if N.LT.2 .)
C
C     X -- (input) real*8 array of independent variable values.  The
C           elements of X must be strictly increasing:
C                X(I-1) .LT. X(I),  I = 2(1)N.
C           (Error return if not.)
C
C     F -- (input) real*8 array of function values.  F(1+(I-1)*INCFD) is
C           the value corresponding to X(I).
C
C     D -- (input) real*8 array of derivative values.  D(1+(I-1)*INCFD)
C           is the value corresponding to X(I).
C
C     INCFD -- (input) increment between successive values in F and D.
C           (Error return if  INCFD.LT.1 .)
C
C     SKIP -- (input/output) logical variable which should be set to
C           .TRUE. if the user wishes to skip checks for validity of
C           preceding parameters, or to .FALSE. otherwise.
C           This will save time in case these checks have already
C           been performed (say, in DPCHIM or DPCHIC).
C           SKIP will be set to .TRUE. on normal return.
C
C     NE -- (input) number of evaluation points.  (Error return if
C           NE.LT.1 .)
C
C     XE -- (input) real*8 array of points at which the function is to
C           be evaluated.
C
C          NOTES:
C           1. The evaluation will be most efficient if the elements
C              of XE are increasing relative to X;
C              that is,   XE(J) .GE. X(I)
C              implies    XE(K) .GE. X(I),  all K.GE.J .
C           2. If any of the XE are outside the interval [X(1),X(N)],
C              values are extrapolated from the nearest extreme cubic,
C              and a warning error is returned.
C
C     FE -- (output) real*8 array of values of the cubic Hermite
C           function defined by  N, X, F, D  at the points  XE.
C
C     IERR -- (output) error flag.
C           Normal return:
C              IERR = 0  (no errors).
C           Warning error:
C              IERR.GT.0  means that extrapolation was performed at
C                 IERR points.
C           "Recoverable" errors:
C              IERR = -1  if N.LT.2 .
C              IERR = -2  if INCFD.LT.1 .
C              IERR = -3  if the X-array is not strictly increasing.
C              IERR = -4  if NE.LT.1 .
C             (The FE-array has not been changed in any of these cases.)
C               NOTE:  The above errors are checked in the order listed,
C                   and following arguments have **NOT** been validated.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  DCHFEV, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   811020  DATE WRITTEN
C   820803  Minor cosmetic changes for release 1.
C   870707  Corrected XERROR calls for d.p. name(s).
C   890206  Corrected XERROR calls.
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890831  Modified array declarations.  (WRB)
C   891006  Cosmetic changes to prologue.  (WRB)
C   891006  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C***END PROLOGUE  DPCHFE
C  Programming notes:
C
C     1. To produce a single precision version, simply:
C        a. Change DPCHFE to PCHFE, and DCHFEV to CHFEV, wherever they
C           occur,
C        b. Change the double precision declaration to real,
C
C     2. Most of the coding between the call to DCHFEV and the end of
C        the IR-loop could be eliminated if it were permissible to
C        assume that XE is ordered relative to X.
C
C     3. DCHFEV does not assume that X1 is less than X2.  thus, it would
C        be possible to write a version of DPCHFE that assumes a
C        decreasing X-array by simply running the IR-loop backwards
C        (and reversing the order of appropriate tests).
C
C     4. The present code has a minor bug, which I have decided is not
C        worth the effort that would be required to fix it.
C        If XE contains points in [X(N-1),X(N)], followed by points .LT.
C        X(N-1), followed by points .GT.X(N), the extrapolation points
C        will be counted (at least) twice in the total returned in IERR.
C
C  DECLARE ARGUMENTS.
C
      INTEGER  N, INCFD, NE, IERR
      DOUBLE PRECISION  X(*), F(INCFD,*), D(INCFD,*), XE(*), FE(*)
      LOGICAL  SKIP
C
C  DECLARE LOCAL VARIABLES.
C
      INTEGER  I, IERC, IR, J, JFIRST, NEXT(2), NJ
C
C  VALIDITY-CHECK ARGUMENTS.
C
C***FIRST EXECUTABLE STATEMENT  DPCHFE
      IF (SKIP)  GO TO 5
C
      IF ( N.LT.2 )  GO TO 5001
      IF ( INCFD.LT.1 )  GO TO 5002
      DO 1  I = 2, N
         IF ( X(I).LE.X(I-1) )  GO TO 5003
    1 CONTINUE
C
C  FUNCTION DEFINITION IS OK, GO ON.
C
    5 CONTINUE
      IF ( NE.LT.1 )  GO TO 5004
      IERR = 0
      SKIP = .TRUE.
C
C  LOOP OVER INTERVALS.        (   INTERVAL INDEX IS  IL = IR-1  . )
C                              ( INTERVAL IS X(IL).LE.X.LT.X(IR) . )
      JFIRST = 1
      IR = 2
   10 CONTINUE
C
C     SKIP OUT OF LOOP IF HAVE PROCESSED ALL EVALUATION POINTS.
C
         IF (JFIRST .GT. NE)  GO TO 5000
C
C     LOCATE ALL POINTS IN INTERVAL.
C
         DO 20  J = JFIRST, NE
            IF (XE(J) .GE. X(IR))  GO TO 30
   20    CONTINUE
         J = NE + 1
         GO TO 40
C
C     HAVE LOCATED FIRST POINT BEYOND INTERVAL.
C
   30    CONTINUE
         IF (IR .EQ. N)  J = NE + 1
C
   40    CONTINUE
         NJ = J - JFIRST
C
C     SKIP EVALUATION IF NO POINTS IN INTERVAL.
C
         IF (NJ .EQ. 0)  GO TO 50
C
C     EVALUATE CUBIC AT XE(I),  I = JFIRST (1) J-1 .
C
C       ----------------------------------------------------------------
        CALL DCHFEV (X(IR-1),X(IR), F(1,IR-1),F(1,IR), D(1,IR-1),D(1,IR)
     *              ,NJ, XE(JFIRST), FE(JFIRST), NEXT, IERC)
C       ----------------------------------------------------------------
         IF (IERC .LT. 0)  GO TO 5005
C
         IF (NEXT(2) .EQ. 0)  GO TO 42
C        IF (NEXT(2) .GT. 0)  THEN
C           IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(2) TO THE
C           RIGHT OF X(IR).
C
            IF (IR .LT. N)  GO TO 41
C           IF (IR .EQ. N)  THEN
C              THESE ARE ACTUALLY EXTRAPOLATION POINTS.
               IERR = IERR + NEXT(2)
               GO TO 42
   41       CONTINUE
C           ELSE
C              WE SHOULD NEVER HAVE GOTTEN HERE.
               GO TO 5005
C           ENDIF
C        ENDIF
   42    CONTINUE
C
         IF (NEXT(1) .EQ. 0)  GO TO 49
C        IF (NEXT(1) .GT. 0)  THEN
C           IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(1) TO THE
C           LEFT OF X(IR-1).
C
            IF (IR .GT. 2)  GO TO 43
C           IF (IR .EQ. 2)  THEN
C              THESE ARE ACTUALLY EXTRAPOLATION POINTS.
               IERR = IERR + NEXT(1)
               GO TO 49
   43       CONTINUE
C           ELSE
C              XE IS NOT ORDERED RELATIVE TO X, SO MUST ADJUST
C              EVALUATION INTERVAL.
C
C              FIRST, LOCATE FIRST POINT TO LEFT OF X(IR-1).
               DO 44  I = JFIRST, J-1
                  IF (XE(I) .LT. X(IR-1))  GO TO 45
   44          CONTINUE
C              NOTE-- CANNOT DROP THROUGH HERE UNLESS THERE IS AN ERROR
C                     IN DCHFEV.
               GO TO 5005
C
   45          CONTINUE
C              RESET J.  (THIS WILL BE THE NEW JFIRST.)
               J = I
C
C              NOW FIND OUT HOW FAR TO BACK UP IN THE X-ARRAY.
               DO 46  I = 1, IR-1
                  IF (XE(J) .LT. X(I)) GO TO 47
   46          CONTINUE
C              NB-- CAN NEVER DROP THROUGH HERE, SINCE XE(J).LT.X(IR-1).
C
   47          CONTINUE
C              AT THIS POINT, EITHER  XE(J) .LT. X(1)
C                 OR      X(I-1) .LE. XE(J) .LT. X(I) .
C              RESET IR, RECOGNIZING THAT IT WILL BE INCREMENTED BEFORE
C              CYCLING.
               IR = MAX(1, I-1)
C           ENDIF
C        ENDIF
   49    CONTINUE
C
         JFIRST = J
C
C     END OF IR-LOOP.
C
   50 CONTINUE
      IR = IR + 1
      IF (IR .LE. N)  GO TO 10
C
C  NORMAL RETURN.
C
 5000 CONTINUE
      RETURN
C
C  ERROR RETURNS.
C
 5001 CONTINUE
C     N.LT.2 RETURN.
      IERR = -1
      CALL XERMSG ('SLATEC', 'DPCHFE',
     +   'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1)
      RETURN
C
 5002 CONTINUE
C     INCFD.LT.1 RETURN.
      IERR = -2
      CALL XERMSG ('SLATEC', 'DPCHFE', 'INCREMENT LESS THAN ONE', IERR,
     +   1)
      RETURN
C
 5003 CONTINUE
C     X-ARRAY NOT STRICTLY INCREASING.
      IERR = -3
      CALL XERMSG ('SLATEC', 'DPCHFE',
     +   'X-ARRAY NOT STRICTLY INCREASING', IERR, 1)
      RETURN
C
 5004 CONTINUE
C     NE.LT.1 RETURN.
      IERR = -4
      CALL XERMSG ('SLATEC', 'DPCHFE',
     +   'NUMBER OF EVALUATION POINTS LESS THAN ONE', IERR, 1)
      RETURN
C
 5005 CONTINUE
C     ERROR RETURN FROM DCHFEV.
C   *** THIS CASE SHOULD NEVER OCCUR ***
      IERR = -5
      CALL XERMSG ('SLATEC', 'DPCHFE',
     +   'ERROR RETURN FROM DCHFEV -- FATAL', IERR, 2)
      RETURN
C------------- LAST LINE OF DPCHFE FOLLOWS -----------------------------
      END