text:
      integer :: m = 12, clock,inirand,i,j,k,l,p,dimn,i1,ntm,mtc,q,r
      integer counter,i2,nrow,tend,count,tspc,flag,indi,indj,sampsz
      integer ii,jj
      integer, parameter ::Totrun=10
      integer, parameter ::burntim=0
      double precision, parameter ::stepszu=1.0d0
      integer N,NRHS,LDA,INFO,LDB,tag,marker,Montrun
      integer, allocatable,dimension(:) :: iseed,IPIV
      double precision x,chi1,chi2,time,sum1,sigma,ratio,chiint,prob,y
      double precision pert,sumM,stepsz,sigmal,alpha
      double precision, allocatable, dimension(:) ::Delta,thetasave
      double precision, allocatable, dimension(:) ::mu2diff,Indic
      double precision, allocatable, dimension(:) :: mutreal,mu0real
      double precision, allocatable, dimension(:) :: mu1t,mu1diff,mu2t
      double precision, allocatable, dimension(:,:) :: Jtreal,J0real,M1
      double precision, allocatable, dimension(:,:) :: M2,J1diff,J2diff
      integer, allocatable, dimension(:,:) :: M1ind,M2ind
      integer, allocatable, dimension(:) :: Accpt,List
      double precision, allocatable, dimension(:,:) :: J1,J2,A1,A2
      double precision, allocatable, dimension(:,:) :: AVG
      double precision, allocatable, dimension(:) :: chisave,EqCORR
      double precision, allocatable, dimension(:,:) ::Areal,M1ini
      double precision, dimension(10,Totrun) :: chiequli,stdchiequli
      integer num_times
cc      integer :: pos1 = 1
cc      character (len=100) :: t_string
      character(len=5),allocatable,dimension(:) :: t_i
      integer :: t_i_int

cc       call system("wc -m <time.txt > t_len")
cc       call system("wc -l <time.txt > num_times")
cc       open(unit=11,file='t_len',status='old')
cc       read(11,*)t_len
cc       close(11)
cc       open(unit=11,file='num_times',status='old')
cc       read(11,*)num_times
cc       close(11)
       num_times = 2

       ALLOCATE(t_i(num_times))
cc       ALLOCATE(t_string(t_len))
       open(unit=11,file='time.txt',status='old',action='read')
cc       read(unit=11,fmt='(A)')t_string
       do j=1,num_times
          read(11,*)t_i(j)
       end do
       close(11)

cc This segment of code below was taken from rosettacode.org.  The section is "Tokenize a string", subsection "Fortran"
cc       DO z = 1, num_times
cc          pos2 = INDEX(t_string(pos1:), " ")
cc          IF (pos2 == 0) THEN
cc             t_i(z) = t_string(pos1:)
cc             EXIT
cc          END IF
cc          t_i(z) = t_string(pos1:pos1+pos2-2)
cc          pos1 = pos2+pos1
cc       END DO

cc BEGIN SAYAK's CODE cc

      call system("wc -l <avg_new.txt>
     & inpl")

      call system("awk 'NR==2{print NF}'
     & avg_new.txt> inpc")

      open(unit=11,file='inpc',status='old')
       read(11,*)dimn
      close(11)

      open(unit=12,file='inpl',status='old')
       read(12,*)nrow
      close(12)
      write(*,*)"row =",nrow

      allocate(EqCORR((nrow)*(dimn)**2))
      allocate(AVG(nrow,dimn))

      open(unit=10,file='equaltime.txt',status='old')
         do i=1,((nrow)*(dimn)**2)
           read(10,*)EqCORR(i)
c            write(*,*)'i =',i
c           write(*,*)'EqCorr = ',EqCorr(i)
c            call sleep(3)
         enddo
      close(10)

      open(unit=11,file='avg.txt',
     &          status='old')
         do i=1,(nrow)
           read(11,*)(AVG(i,j),j=1,dimn)
         enddo
      close(11)

cc      open(unit=12,file='tend.txt',status='old')
cc           read(12,*)tend,tspc
cc           write(*,*)"tend=",tend,"tspc=",tspc
cc      close(12)
       tend=1
       tspc=1

c This is the part Sayak pointed out needs changing

       read(t_i(tend),*) t_i_int1
       read(t_i(tend+1),*) t_i_int2
       time=t_i_int2-t_i_int1

cc       if(tend.eq.1) then
cc       time=24.0d0
cc       else if(tend.eq.2) then
cc       time = 32.0d0
cc       else if(tend.eq.3) then
cc       time = 64.0d0
cc       else if(tend.eq.4) then
cc       time = 128.0d0
cc       endif
       write(*,*)"time =",time

c *******    Initialize the random number generator ********

      allocate(iseed(m))
      call random_seed(size = m)

      call system_clock(COUNT=clock)

      iseed = clock + 37 * [(i1, i1 = 0,m-1)]

      call random_seed(PUT = iseed)

c **************************************************************

c **************    Allocate arrays  **************************

      allocate(mutreal(dimn),mu0real(dimn),mu1t(dimn),mu1diff(dimn))
      allocate(Jtreal(dimn,dimn),J0real(dimn,dimn),M1(dimn,dimn))
      allocate(M2(dimn,dimn),A1(dimn,dimn),A2(dimn,dimn),J1(dimn,dimn))
      allocate(J2(dimn,dimn),M1ind(dimn,dimn),M2ind(dimn,dimn))
      allocate(J1diff(dimn,dimn),mu2diff(dimn),J2diff(dimn,dimn))
      allocate(chisave(Totrun-burntim),mu2t(dimn))
      allocate(Areal(dimn,dimn),Delta(dimn),IPIV(dimn))
      allocate(M1ini(dimn,dimn),Accpt(Totrun))
      allocate(Indic(Totrun))

c ******************************************************************

c *******  Initialize the arrays given from data *******************

      mutreal=0.0d0
      mu0real=0.0d0
      Jtreal=0.0d0
      J0real=0.0d0
      IPIV=0

c ******************************************************************

c **************** Read initial data *******************************

      p=(tend-1)*(dimn)**2+1

      do i=1,dimn
         mu0real(i)=AVG(tend,i)
         mutreal(i)=AVG(tend+tspc,i)
         do j=1,dimn
          J0real(i,j)=EqCORR(p)
c          write(*,*)'p =',p
c          write(*,*)'J0real =',J0real(i,j)
          Jtreal(i,j)=EqCORR(tspc*(dimn)**2+p)
c          write(*,*)'dimn =',dimn
c          write(*,*)'Jtreal =',Jtreal(i,j)
c          call sleep(5)
          p=p+1
         enddo
      enddo

c      write(*,*)"Mmat of J0="
c      do i=1,dimn
c          write(*,'(1X,512F10.4)')(J0real(i,j),j=1,dimn)
c      enddo
c
c      write(*,*)"Mmat of Jt="
c      do i=1,dimn
c          write(*,'(1X,512F10.4)')(Jtreal(i,j),j=1,dimn)
c      enddo

c       write(*,*)" C0="
c       write(*,'(1X,512F10.5)')(mu0real(j),j=1,dimn)

c       write(*,*)"C1="
c       write(*,'(1X,512F10.5)')(mutreal(j),j=1,dimn)

      M1=0.0d0

      write(*,*)"dimn=",dimn

       open(unit=21,file='Minitial.txt',status='old')
         do i=1,dimn
           read(21,*)(M1(i,j),j=1,dimn)
         enddo
      close(21)

cc      open(unit=15,file='prob.txt',status='old')
cc           read(15,*)prob
cc      close(15)
       prob = 0.0d0

c      open(unit=15,file='step.txt',status='old')
c           read(15,*)stepsz
c      close(15)
       stepsz = 0.01d0
       write(*,*)"stepsz= ",stepsz
      do i=1,dimn
       do j=1,dimn
          if(i.ne.j) then
             if(M1(i,j).lt.prob) then
                 M1(i,j)=0.0d0*x
                endif
             endif
        enddo
       enddo

        do j=1,dimn
             sum1=0.0d0
             M1(j,j)=0.0d0
             do i=1,dimn
                 sum1=sum1+M1(i,j)
             enddo
             M1(j,j)=-1.0*sum1
         enddo

      write(*,*)"Mmat of choice="
      do i=1,dimn
          write(*,'(1X,512F10.5)')(M1(i,j),j=1,dimn)
      enddo

      Allocate(List(dimn*dimn))
      tag=0
      List=0
      marker=0
      do i=1,dimn
         do j=1,dimn
            marker=marker+1
c            if(i.ne.j) then
                if(M1(i,j).ne.0.0) then
                  tag=tag+1
                  List(tag)=marker
                endif
c             endif
          enddo
       enddo

c       Montrun=10000*tag         !assigning Monte Carlo Steps
        Montrun=1000000
c        allocate(thetasave(Montrun))
c  *******************************************************************

c **************  Initial distance calculation  **********************************

      chi1=0.0d0
      A1=0.0d0
      A1=expm(time,M1,dimn)
      mu1t=0.0d0
      mu1diff=0.0d0
      J1=0.0d0
      J1diff=0.0d0

      do i=1,dimn
         do j=1,dimn
c          write(*,*)'A1 =',A1(i,j)
          mu1t(i)=mu1t(i)+A1(i,j)*mu0real(j)
c          write(*,*)'mu1t =',mu1t(i)
c          call sleep(5)
         enddo
         mu1diff(i)=1.0d0-mu1t(i)/mutreal(i)
c          write(*,*)'mu1diff =',mu1diff(i)
c          call sleep(5)
       enddo

      do i=1,dimn
         do j=i,dimn
            count=count+1
               do l=1,dimn
                    do p=1,dimn
                      J1(i,j)=J1(i,j)+A1(i,l)*J0real(l,p)*A1(j,p)
c                      write(*,*)'J1'
c                      write(*,*)J1(i,j)
                    enddo
                enddo
c            write(*,*)'Jtreal =',Jtreal(i,j)
c            write(*,*)Jtreal(i,j)
c            write(*,*)'J1diff'
             J1diff(i,j)=1.0d0-J1(i,j)/Jtreal(i,j)
c            write(*,*)'J1diff =',J1diff(i,j)
c            call sleep(5)
         enddo
      enddo

      do i=1,dimn
         do j=i,dimn
           chi1=chi1+J1diff(i,j)**2
c            write(*,*)'part 1'
c            write(*,*)chi1
         enddo

           chi1=chi1+mu1diff(i)**2
c           write(*,*)'part2'
c           write(*,*)chi1
      enddo

      sigma=10*chi1
      write(*,*)"sigma=",sigma

      write(*,*)"chi1 =",chi1

      alpha=dexp((dlog(0.000001/sigma))/(1.0d0*Totrun))
      write(*,*)"alpha= ",alpha

c      call sleep(300)

c      write(*,*)"tag=",tag
c      open(unit=20,file='prob_new.txt',status='old')
c      read(20,*)prob
c      close(20)

c      open(unit=21,file='pert.txt',status='old')
c      read(21,*)pert
c      close(21)

c ***********************************************************************************

c *********************  Temperature loop *********************************************

c HERE'S SOME THINGS I ALLOCATED IN A DEBUGGING ATTEMPT
c      allocate(Indic(1))
      allocate(thetasave(1))
c      allocate(Accpt(1))
c      allocate(chisave(1))
c      allocate(

      chisave=0.0d0
      Accpt=0
      chiequli=0.0d0
      stdchiequli=0.0d0
      Indic=0.0d0

c     Beginning of temperture loop
      do ntm=1,Totrun
c	   prog=ntm/Totrun
c	   write(prog_str,’(F3.1)’)prog
        write(*,'(a)',advance='no')"Run "
        write(*,'(i4)',advance='no')ntm
        write(*,'(a)',advance='no')" out of "
        write(*,'(i4)',advance='no')Totrun
          thetasave=0.0d0
          Indic(ntm)=sigma
         if(MOD(ntm,1000)==0) then
            write(*,*) ntm
         endif
         flag=0
         sampsz=0
c  Beginning of Monte Carlo loop
         do mtc=1,Montrun
               M2=0.0d0
               call random_number(x)

               if(MOD(List(int(1.0+1.0d0*(tag)*x)),dimn).eq.0)  then
                    indi=int((1.0d0*List(int(1.0+1.0d0*(tag)*x)))/
     &                                                 (dimn*1.0d0))
                    indj=dimn
               else
                    indi= int((1.0d0*List(int(1.0+1.0d0*(tag)*x)))/
     &                                               (dimn*1.0d0))+1
                    indj= int(MOD(List(int(1.0+1.0d0*(tag)*x)),dimn))
               endif
c               write(*,*)"indi=",indi,"indj=",indj

               M2=M1
               call random_number(x)
               if(indi.eq.indj) then
               M2(indi,indj)=M1(indi,indj)+stepsz*(2.0d0*x-1.0d0)
               else
c                write(*,*)M1(indi,indj)+ stepsz*(2.0d0*x-1.0d0)
                  if(((M1(indi,indj)+ stepsz*(2.0d0*x-1.0d0)).lt.0.0)
     &          .OR.((M1(indi,indj)+ stepsz*(2.0d0*x-1.0d0)).gt.1.0))
     &                 then
                      M2(indi,indj)=M1(indi,indj)
                   else
                      M2(indi,indj)=M1(indi,indj)+stepsz*(2.0d0*x-1.0d0))
                  endif

               endif

c                write(*,*)"Mmat of M2="
c                do i=1,dimn
c                   write(*,'(1X,512F10.5)')(M2(i,j),j=1,dimn)
c                 enddo

                chi2=0.0d0
                A2=0.0d0
                A2=expm(time,M2,dimn)
                mu2t=0.0d0
                mu2diff=0.0d0
                J2=0.0d0
                J2diff=0.0d0

                do i=1,dimn
                  do j=1,dimn
                     mu2t(i)=mu2t(i)+A2(i,j)*mu0real(j)
                  enddo
                 mu2diff(i)=1.0d0-mu2t(i)/mutreal(i)
                enddo

                do i=1,dimn
                     do j=i,dimn
                       count=count+1
                          do l=1,dimn
                              do p=1,dimn
                                J2(i,j)=J2(i,j)+A2(i,l)*J0real(l,p)
     &                                                     *A2(j,p)
                              enddo
                          enddo
                       J2diff(i,j)=1.0d0-J2(i,j)/Jtreal(i,j)
                     enddo
                enddo

                do i=1,dimn
cccTHIS IS THE SPOT YOU CHANGED TO REMOVE THE COVARIANCE FROM THE COST FUNCTION

                   do j=i,dimn
                        chi2=chi2+J2diff(i,j)**2
                   enddo

                      chi2=chi2+mu2diff(i)**2
                enddo

                do j=1,dimn
                  sumM=0.0d0
                  do i=1,dimn
                    sumM=sumM+M2(i,j)
                  enddo
                  chi2=chi2+(sumM**2)
                enddo
c                write(*,*)"chi2= ",chi2
                ratio=dexp((-chi2+chi1)/(2.0d0*sigma))

                call random_number(x)
c                write(*,*)"ratio=",ratio,"x=",x,"chi2=",chi2
                if(x.lt.ratio) then
c                   write(*,*) "yeah"
                   flag=flag+1
                   M1=M2
                   chi1=chi2
                 endif

                 if(mtc.gt.Montrun-100000) then
                    chiequli(int(sampsz/10000)+1,ntm)=
     &                  chiequli(int(sampsz/10000)+1,ntm)+chi1
                    stdchiequli(int(sampsz/10000)+1,ntm)=
     &                  stdchiequli(int(sampsz/10000)+1,ntm)+chi1**2
                    sampsz=sampsz+1

                 endif

c             thetasave(mtc)=chi1
         enddo   !end of Monte Carlo run
         sigma=alpha*sigma
         Accpt(ntm)=flag
             if(ntm.gt.burntim) then
c                   counter=0
c                    do i=1,dimn
c                         do j=1,dimn
c                            if(i.ne.j) then
c                            counter=counter+1
c                             thetasave(ntm-burntim,counter)=M1(i,j)
c                            endif
c                         enddo
c                    enddo
                chisave(ntm-burntim)=chi1
             endif

      write(*,'(a)',advance='no')char(13)
      enddo    ! end of temperature loop
      write(*,*)"sigma=",sigma

      open(unit=41,status='unknown')
        do i=1,dimn
           write(41,'(1X,512F20.10)')(M1(i,j),j=1,dimn)
        enddo
      close(41)

      i2=Totrun-burntim

      open(unit=42,status='unknown')
        do i=1,i2
c           write(*,*)chisave(i)
           write(42,'(1X,512F30.10)')chisave(i)
        enddo
      close(42)

      open(unit=40,status='unknown')
        do i=1,Totrun
           write(40,'(1X,512I7)')Accpt(i)
         enddo
       close(40)

      open(unit=49,status='unknown')
        do i=1,10
           write(49,'(1X,100000F30.10)')(chiequli(i,j)/10000.0d0,
     &                                        j=1,Totrun)
         enddo
       close(49)

      open(unit=59,status='unknown')
        do i=1,10
          write(59,'(1X,100000F30.10)')(sqrt(stdchiequli(i,j)/10000.0d0
     &                      -(chiequli(i,j)/10000.0d0)**2),j=1,Totrun)
         enddo
       close(59)

c      open(unit=50,status='unknown')
c        do i=1,Montrun
c           write(50,'(1X,512F20.10)')thetasave(i)
c        enddo
c       close(50)

      open(unit=51,status='unknown')
        do i=1,Totrun
           write(51,'(1X,512F30.15)')Indic(i)
        enddo
       close(51)

      contains
         function expm(t,H,N1) result(expH)
         double precision, intent(in):: t
         integer, intent(in)::N1
         double precision, dimension(N1,N1), intent(in)::H
         double precision, dimension(size(H,1),size(H,2)) :: expH
         external :: DGPADM
          integer, parameter :: ideg = 6
          double precision, dimension(4*size(H,1)*size(H,2) + ideg + 1)
     &       :: wsp
          integer, dimension(size(H,1))  :: iwsp
          integer :: iexp, ns, iflag, n
          if (size(H,1) /= size(H,2)) then
          stop 'expm: matrix must be square'
          end if
          n = size(H,1)
          call DGPADM(ideg, n, t, H, n, wsp, size(wsp,1), iwsp, iexp,
     &       ns,iflag)
          expH = reshape(wsp(iexp:iexp+n*n-1), shape(expH))
         end function expm

 10    end

c----------------------------------------------------------------------|
      subroutine DGPADM( ideg,m,t,H,ldh,wsp,lwsp,ipiv,iexph,ns,iflag )

      implicit none
      integer ideg, m, ldh, lwsp, iexph, ns, iflag, ipiv(m)
      double precision t, H(ldh,m), wsp(lwsp)

c-----Purpose----------------------------------------------------------|
c
c     Computes exp(t*H), the matrix exponential of a general matrix in
c     full, using the irreducible rational Pade approximation to the
c     exponential function exp(x) = r(x) = (+/-)( I + 2*(q(x)/p(x)) ),
c     combined with scaling-and-squaring.
c
c-----Arguments--------------------------------------------------------|
c
c     ideg      : (input) the degre of the diagonal Pade to be used.
c                 a value of 6 is generally satisfactory.
c
c     m         : (input) order of H.
c
c     H(ldh,m)  : (input) argument matrix.
c
c     t         : (input) time-scale (can be < 0).
c
c     wsp(lwsp) : (workspace/output) lwsp .ge. 4*m*m+ideg+1.
c
c     ipiv(m)   : (workspace)
c
c>>>> iexph     : (output) number such that wsp(iexph) points to exp(tH)
c                 i.e., exp(tH) is located at wsp(iexph ... iexph+m*m-1)
c                       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
c                 NOTE: if the routine was called with wsp(iptr),
c                       then exp(tH) will start at wsp(iptr+iexph-1).
c
c     ns        : (output) number of scaling-squaring used.
c
c     iflag     : (output) exit flag.
c                      0 - no problem
c                     <0 - problem
c
c----------------------------------------------------------------------|
c     Roger B. Sidje (rbs@maths.uq.edu.au)
c     EXPOKIT: Software Package for Computing Matrix Exponentials.
c     ACM - Transactions On Mathematical Software, 24(1):130-156, 1998
c----------------------------------------------------------------------|
c
      integer mm,i,j,k,ih2,ip,iq,iused,ifree,iodd,icoef,iput,iget
      double precision hnorm,scale,scale2,cp,cq

      intrinsic INT,ABS,DBLE,LOG,MAX

c---  check restrictions on input parameters ...
      mm = m*m
      iflag = 0
      if ( ldh.lt.m ) iflag = -1
      if ( lwsp.lt.4*mm+ideg+1 ) iflag = -2
      if ( iflag.ne.0 ) stop 'bad sizes (in input of DGPADM)'
c
c---  initialise pointers ...
c
      icoef = 1
      ih2 = icoef + (ideg+1)
      ip  = ih2 + mm
      iq  = ip + mm
      ifree = iq + mm
c
c---  scaling: seek ns such that ||t*H/2^ns|| < 1/2;
c     and set scale = t/2^ns ...
c
      do i = 1,m
         wsp(i) = 0.0d0
      enddo
      do j = 1,m
         do i = 1,m
            wsp(i) = wsp(i) + ABS( H(i,j) )
         enddo
      enddo
      hnorm = 0.0d0
      do i = 1,m
         hnorm = MAX( hnorm,wsp(i) )
      enddo
      hnorm = ABS( t*hnorm )
      if ( hnorm.eq.0.0d0 ) stop 'Error - null H in input of DGPADM.'
      ns = MAX( 0,INT(LOG(hnorm)/LOG(2.0d0))+2 )
      scale = t / DBLE(2**ns)
      scale2 = scale*scale
c
c---  compute Pade coefficients ...
c
      i = ideg+1
      j = 2*ideg+1
      wsp(icoef) = 1.0d0
      do k = 1,ideg
         wsp(icoef+k) = (wsp(icoef+k-1)*DBLE( i-k ))/DBLE( k*(j-k) )
      enddo
c
c---  H2 = scale2*H*H ...
c
      call DGEMM( 'n','n',m,m,m,scale2,H,ldh,H,ldh,0.0d0,wsp(ih2),m )
c
c---  initialize p (numerator) and q (denominator) ...
c
      cp = wsp(icoef+ideg-1)
      cq = wsp(icoef+ideg)
      do j = 1,m
         do i = 1,m
            wsp(ip + (j-1)*m + i-1) = 0.0d0
            wsp(iq + (j-1)*m + i-1) = 0.0d0
         enddo
         wsp(ip + (j-1)*(m+1)) = cp
         wsp(iq + (j-1)*(m+1)) = cq
      enddo
c
c---  Apply Horner rule ...
c
      iodd = 1
      k = ideg - 1
 100  continue
      iused = iodd*iq + (1-iodd)*ip
      call DGEMM( 'n','n',m,m,m, 1.0d0,wsp(iused),m,
     .             wsp(ih2),m, 0.0d0,wsp(ifree),m )
      do j = 1,m
         wsp(ifree+(j-1)*(m+1)) = wsp(ifree+(j-1)*(m+1))+wsp(icoef+k-1)
      enddo
      ip = (1-iodd)*ifree + iodd*ip
      iq = iodd*ifree + (1-iodd)*iq
      ifree = iused
      iodd = 1-iodd
      k = k-1
      if ( k.gt.0 )  goto 100
c
c---  Obtain (+/-)(I + 2*(p\q)) ...
c
      if ( iodd .eq. 1 ) then
         call DGEMM( 'n','n',m,m,m, scale,wsp(iq),m,
     .                H,ldh, 0.0d0,wsp(ifree),m )
         iq = ifree
      else
         call DGEMM( 'n','n',m,m,m, scale,wsp(ip),m,
     .                H,ldh, 0.0d0,wsp(ifree),m )
         ip = ifree
      endif
      call DAXPY( mm, -1.0d0,wsp(ip),1, wsp(iq),1 )
      call DGESV( m,m, wsp(iq),m, ipiv, wsp(ip),m, iflag )
      if ( iflag.ne.0 ) stop 'Problem in DGESV (within DGPADM)'
      call DSCAL( mm, 2.0d0, wsp(ip), 1 )
      do j = 1,m
         wsp(ip+(j-1)*(m+1)) = wsp(ip+(j-1)*(m+1)) + 1.0d0
      enddo
      iput = ip
      if ( ns.eq.0 .and. iodd.eq.1 ) then
         call DSCAL( mm, -1.0d0, wsp(ip), 1 )
         goto 200
      endif
c
c--   squaring : exp(t*H) = (exp(t*H))^(2^ns) ...
c
      iodd = 1
      do k = 1,ns
         iget = iodd*ip + (1-iodd)*iq
         iput = (1-iodd)*ip + iodd*iq
         call DGEMM( 'n','n',m,m,m, 1.0d0,wsp(iget),m, wsp(iget),m,
     .                0.0d0,wsp(iput),m )
         iodd = 1-iodd
      enddo
 200  continue
      iexph = iput
      END
c----------------------------------------------------------------------|


! { dg-do compile }
!
! PR fortran/51605
!

contains
  subroutine foo
    BLOCK_NAME: block
    end block BLOCK_NAME
  end subroutine foo

  subroutine BLOCK_NAME()
  end subroutine BLOCK_NAME

  subroutine bar()
  end subroutine bar
end

subroutine test()
contains
  subroutine BLOCK_NAME()
  end subroutine BLOCK_NAME

  subroutine foobar()
  end subroutine foobar

  subroutine foo
    BLOCK_NAME: block
    end block BLOCK_NAME
  end subroutine foo

  subroutine bar()
  end subroutine bar
end


text:
  !! Local Data
  LOGICAL :: bubble
  LOGICAL :: swap_occured
  INTEGER, DIMENSION(:), ALLOCATABLE :: values_per_row
  INTEGER, DIMENSION(:), ALLOCATABLE :: offset_array
  INTEGER, DIMENSION(:), ALLOCATABLE :: inserted_per_row
  !! Counters and temporary variables
  INTEGER :: counter
  INTEGER :: temp_index
  INTEGER :: alloc_stat
  INTEGER :: list_length

  IF (PRESENT(bubble_in)) THEN
     bubble = bubble_in
  ELSE
     bubble = .TRUE.
  END IF

  list_length = input_list%CurrentSize

  IF (bubble .AND. list_length .GT. matrix_rows*matrix_columns*0.1) THEN
     CALL SortDenseTripletList(input_list, matrix_columns, matrix_rows, &
          & sorted_list)
  ELSE
     !! Data Allocation
     CALL ConstructTripletList(sorted_list, list_length)
     ALLOCATE(values_per_row(matrix_columns), stat=alloc_stat)
     ALLOCATE(offset_array(matrix_columns), stat=alloc_stat)
     ALLOCATE(inserted_per_row(matrix_columns), stat=alloc_stat)

     !! Initial one dimensional sort
     values_per_row = 0
     inserted_per_row = 0

     !! Do a first pass bucket sort
     DO counter = 1, input_list%CurrentSize
        values_per_row(input_list%data(counter)%index_column) = &
             & values_per_row(input_list%data(counter)%index_column) + 1
     END DO
     offset_array(1) = 1
     DO counter = 2, UBOUND(offset_array,dim=1)
        offset_array(counter) = offset_array(counter-1) + &
             & values_per_row(counter-1)
     END DO
     DO counter = 1, input_list%CurrentSize
        temp_index = input_list%data(counter)%index_column
        sorted_list%data(offset_array(temp_index)+inserted_per_row(temp_index))=&
             & input_list%data(counter)
        inserted_per_row(temp_index) = inserted_per_row(temp_index) + 1
     END DO

     !! Finish with bubble sort
     !! Not necessary for transposing or unpacking.
     swap_occured = .TRUE.
     IF (bubble) THEN
        DO WHILE (swap_occured .EQV. .TRUE.)
           swap_occured = .FALSE.
           DO counter = 2, sorted_list%CurrentSize
              IF (CompareTriplets(sorted_list%data(counter-1), &
                   & sorted_list%data(counter))) THEN
                 temporary = sorted_list%data(counter)
                 sorted_list%data(counter) = sorted_list%data(counter-1)
                 sorted_list%data(counter-1) = temporary
                 swap_occured = .TRUE.
              END IF
           END DO
        END DO
     END IF

     !! Cleanup
     DEALLOCATE(values_per_row)
     DEALLOCATE(offset_array)
     DEALLOCATE(inserted_per_row)
  END IF

max_stars_repo_name: shoshijak/NTPoly

max_issues_repo_name: omni-compiler/NTPoly

max_forks_repo_name: omni-compiler/NTPoly


text:
!Crown Copyright 2012 AWE.
!
! This file is part of CloverLeaf.
!
! CloverLeaf is free software: you can redistribute it and/or modify it under
! the terms of the GNU General Public License as published by the
! Free Software Foundation, either version 3 of the License, or (at your option)
! any later version.
!
! CloverLeaf 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
! CloverLeaf. If not, see http://www.gnu.org/licenses/.

!>  @brief String manipulation utilities
!>  @author Wayne Gaudin
!>  @details Provides utilities to manipulate and parse Fortran strings.

MODULE clover_case_change

USE data_module

CONTAINS

FUNCTION tolower(string) RESULT (tolower_result)

  IMPLICIT NONE
  CHARACTER (LEN=*), INTENt(IN) :: string
  CHARACTER (LEN=LEN(string)) :: tolower_result
  INTEGER         :: i,ii

  DO i=1,len(string)
    ii=IACHAR(string(i:i))
    SELECT CASE (ii)
      CASE (65:90)
        tolower_result(i:i)=ACHAr(ii+32)
      CASE DEFAULT
        tolower_result(i:i)=string(i:i)
    END SELECT
  ENDDO

END FUNCTION TOLOWER

FUNCTION toupper(string) RESULT (toupper_result)

  IMPLICIT NONE
  CHARACTER (LEN=*), INTENT(IN) :: string
  CHARACTER (LEN=LEN(string)) :: toupper_result
  INTEGER         :: i,ii

  DO i=1,LEN(string)
    ii=IACHAR(string(i:i))
    SELECT CASE (ii)
      CASE (97:122)
        toupper_result(i:i)=ACHAR(ii-32)
      CASE DEFAULT
        toupper_result(i:i)=string(i:i)
    END SELECT
  ENDDO

END FUNCTION toupper

END MODULE clover_case_change

MODULE clover_isitanint_mod

CONTAINS

FUNCTION isitanint(instring) RESULT(isitanint_result)

  IMPLICIT NONE

  CHARACTER(LEN=*), INTENT(IN) :: instring
  LOGICAL :: isitanint_result

  INTEGER         :: i,ii

  isitanint_result=.TRUE.

  DO i=1,LEN(instring)
    ii=IACHAR(instring(i:i))
    SELECT CASE(ii)
      CASE (43,45,48:57)
        IF(i.NE.1) THEN
          IF(ii.EQ.43 .OR.ii.EQ.45) isitanint_result=.FALSE.
        ENDIF
      CASE DEFAULT
        isitanint_result=.FALSE.
    END SELECT
  ENDDO

END FUNCTION isitanint

END MODULE clover_isitanint_mod

MODULE parse_module

  USE data_module
  USE report_module

  IMPLICIT NONE

  INTEGER, PARAMETER :: len_max=g_len_max &
                       ,dummy=0

  INTEGER         :: iu ! Current unit number.

  CHARACTER(LEN=len_max) :: mask   &
                           ,line   &
                           ,here   &
                           ,sel    &
                           ,rest
CONTAINS

FUNCTION parse_init(iunit,cmask)

  ! Initialise for new set of reads.

  IMPLICIT NONE

  INTEGER :: parse_init

  INTEGER         :: iunit
  CHARACTER(LEN=*)  :: cmask

  INTEGER :: ios

  iu=iunit
  mask=cmask ! Set mask for which part of the file we are interested in.
  line=''
  here=''
  rest=''

  REWIND(UNIT=iunit,IOSTAT=ios)

  parse_init=ios

END FUNCTION parse_init

FUNCTION parse_getline(dummy)

  IMPLICIT NONE

  INTEGER :: parse_getline
  INTEGER        ,INTENT(IN) :: dummy

  INTEGER                :: s,ios,i,parse_out
  CHARACTER(LEN=len_max) :: l,nugget,string_temp1,string_temp2

  DO
    READ(UNIT=iu,IOSTAT=ios,FMT='(a100)') l ! Read in next line.

    parse_getline=ios

    IF(parse_getline.NE.0) RETURN

      DO i=1,len_trim(l)
        if (IACHAR(l(i:)).LT.32.OR.IACHAR(l(i:i)).GT.128) l(i:i)=' '
      ENDDO

      l=TRIM(ADJUSTL(l))

      s=SCAN(l,'!')
      IF(s.GT.0) l=TRIM(l(1:s-1))
      s=scan(l,';')
      IF(s.GT.0) l=TRIM(l(1:s-1))

      IF( ABS(dummy).NE.1) THEN
         DO i=1,LEN(l)
            IF(IACHAR(l(i:i)).GT.64.AND.IACHAR(l(i:i)).LT.91) &
                 l(i:i)=ACHAR(IACHAR(l(i:i))+32)
         ENDDO
      ENDIF

      IF(dummy.LT.0)THEN
        line=l
      ELSE

        DO i=1,LEN(l)
          IF(l(i:i).EQ.'='.OR.l(i:i).EQ.',') l(i:i)=' '
       ENDDO

       IF(l(1:8).NE.'*select:')THEN
         IF(l(1:1).EQ.'*')THEN
           s=SCAN(l,' ')
           IF(s.EQ.0) s=LEN_TRIM(l)+1
             nugget=TRIM(l(1:s-1))
             IF(nugget(2:4).EQ.'end')THEN
               nugget='*'//TRIM(nugget(5:))

               s=SCAN(here,'*',.TRUE.)

               parse_out=parse_scan(sel,nugget)
               IF(sel.NE.''.AND.parse_out.GT.0)THEN
                 nugget=''
               ENDIF

               IF(nugget.GT.'')THEN
                 string_temp1=TRIM(here(s:))
                 string_temp2=TRIM(nugget(1:))
                 IF(string_temp1.NE.string_temp2)THEN
                   WRITE(*,*) 'l:      ',trim(adjustl(l))
                   WRITE(*,*) 'nugget: ',trim(adjustl(nugget))
                   WRITE(*,*) 'here:   ',trim(adjustl(here))
                   WRITE(*,*) 'sel:    ',trim(adjustl(sel))
                   CALL report_error('parse_getline','Unmatched */*end pair.')
                 ELSE
                   IF(mask.EQ.here) THEN
                     here=here(1:s-1)
                     parse_getline=-1
                     rest='returned_before_eof'
                     RETURN
                   ENDIF
                     here=here(1:s-1)
                   ENDIF
                 ENDIF
               ELSE
                 parse_out=parse_scan(sel,nugget)
                 IF(sel.NE.''.AND.parse_out.GT.0)THEN
                   nugget=''
                 ENDIF

                 here=TRIM(TRIM(here)//nugget)
                 IF(rest.EQ.'returned_before_eof') rest=''
                 rest=TRIM(TRIM(rest)//l(s:))
              ENDIF
                line=''
            ELSE
              IF(here.EQ.mask.OR.mask.EQ.'')THEN
                line=TRIM(l)
              ELSE
                line=''
              ENDIF
            ENDIF
          ELSE
            sel=TRIM(l(9:))
            line=''
          ENDIF
        ENDIF

        IF(line.GT.'') EXIT
     ENDDO

  END FUNCTION parse_getline

  RECURSIVE FUNCTION parse_getword(wrap) RESULT(getword)

    INTEGER                :: stat
    LOGICAL                :: wrap
    CHARACTER(LEN=len_max) :: getword,temp

    INTEGER :: s

    DO WHILE(line(1:1).EQ.' '.AND.LEN_TRIM(line).GT.0)
      line=TRIM(line(2:))
    ENDDO

    s=SCAN(line,' ')
    IF(s.EQ.0) s=LEN_TRIM(line)

    temp=TRIM(line(1:s))

    IF(temp.EQ.ACHAR(92).OR.(temp.EQ.''.AND.wrap))THEN
       temp=''
       stat=parse_getline(dummy)
       IF(stat.EQ.0) temp=parse_getword(wrap)
       getword=temp
    ELSE
       getword=temp

       line=TRIM(line(s+1:))
    ENDIF

  END FUNCTION parse_getword

  FUNCTION parse_getival(word)

    USE clover_module

    CHARACTER(LEN=*)  :: word
    INTEGER         :: temp,parse_getival

    INTEGER :: ios

    READ(UNIT=word,FMT="(I7)",IOSTAT=ios) temp

    IF(ios.NE.0)THEN
       CALL report_error('parse_getival','Error attempting to convert to integer:'//word)
       CALL clover_abort
    ENDIF

    parse_getival=temp

  END FUNCTION parse_getival

  FUNCTION parse_getlval(word)

    USE clover_module

    CHARACTER(LEN=*)  :: word
    LOGICAL :: temp,parse_getlval

    INTEGER :: ios

    ios=0

    SELECT CASE(word)
       CASE('on')
          temp=.TRUE.
       CASE('true')
          temp=.TRUE.
       CASE('off')
          temp=.FALSE.
       CASE('false')
          temp=.FALSE.
       CASE DEFAULT
          ios=99999
    END SELECT

    IF(ios.NE.0)THEN
       CALL report_error('parse_getlval','Error attempting to convert to logical:'//word)
       CALL clover_abort
    ENDIF

    parse_getlval=temp

  END FUNCTION parse_getlval

  FUNCTION parse_getrval(word)

    USE clover_module

    CHARACTER(LEN=*) :: word
    REAL(KIND=8)   :: temp,parse_getrval

    INTEGER :: ios

    ! Make an integer into a float if necessary.

    IF(SCAN(word,'.').EQ.0) word=TRIM(TRIM(word)//'.0')

    READ(UNIT=word,FMT="(E27.20)",IOSTAT=ios) temp

    IF(ios.NE.0)THEN
       CALL report_error('parse_getrval','Error attempting to convert to real:'//word)
       CALL clover_abort
    ENDIF

    parse_getrval=temp

  END FUNCTION parse_getrval

  FUNCTION parse_scan(string,set)

    ! Improved version of F90 SCAN.

    INTEGER                     :: parse_scan
    CHARACTER(LEN=*),INTENT(IN) :: string,set
    CHARACTER(LEN=LEN_MAX)      :: set_temp

    INTEGER :: i,l,temp

    l=LEN_TRIM(set)-1

    temp=0
    set_temp=TRIM(set)
    DO i=1,LEN_TRIM(string)-l
       IF(string(i:i+l).EQ.set_temp)THEN
          temp=i
          EXIT
       ENDIF
    ENDDO

    parse_scan=temp

  END FUNCTION parse_scan

END MODULE parse_module

max_stars_repo_name: uwplse/stng

max_issues_repo_name: uwplse/stng

max_forks_repo_name: uwplse/stng


text:
subroutine hus_vertical_interpolation(hus_on_model_level,pres, &
	pressure_full_level,m,n,o,p,req,hus_on_press_level)
implicit none
integer :: m,n,o,p,req
integer :: x,y,s,t,plev
double precision :: pres(req),hus_on_model_level(m,n,o,p)
double precision :: pressure_full_level(m,n,o,p)
double precision :: delta_hus,delta_p,grad_hus_p,diff_p
double precision, intent(out) :: hus_on_press_level(m,n,req,p)
real :: arg = -1.0,NaN
NaN= sqrt(arg)
do plev=1,req
 do t=1,p
  do x=1,m
   do y=1,n
    do s=1,o
!above uppest level
    if(pres(plev) .LT. pressure_full_level(x,y,1,t)) then
     hus_on_press_level(x,y,plev,t) = NaN
    end if
! in between levels
    if(pres(plev) .GE. pressure_full_level(x,y,s,t) .AND. pres(plev) .LE. &
     pressure_full_level(x,y,s+1,t) ) then
     delta_hus = hus_on_model_level(x,y,s,t) - hus_on_model_level(x,y,s+1,t)
     delta_p = log(pressure_full_level(x,y,s,t))&
     - log(pressure_full_level(x,y,s+1,t))
     grad_hus_p = delta_hus /delta_p
     diff_p = log(pres(plev)) - log(pressure_full_level(x,y,s,t))
     hus_on_press_level(x,y,plev,t) = hus_on_model_level(x,y,s,t)&
     + grad_hus_p * diff_p
    end if
! extrapolation below the ground
    if(pres(plev) .GT. pressure_full_level(x,y,o,t)) then
     hus_on_press_level(x,y,plev,t) = hus_on_model_level(x,y,o,t)
    end if
    end do
   end do
  end do
 end do
end do
end subroutine hus_vertical_interpolation

ext: f


text:
      SUBROUTINE HAFTON (Z,L,M,N,FLO,HI,NLEV,NOPT,NPRM,ISPV,SPVAL)
C
C
C +-----------------------------------------------------------------+
C |                                                                 |
C |                Copyright (C) 1986 by UCAR                       |
C |        University Corporation for Atmospheric Research          |
C |                    All Rights Reserved                          |
C |                                                                 |
C |                 NCARGRAPHICS  Version 1.00                      |
C |                                                                 |
C +-----------------------------------------------------------------+
C
C
C     SUBROUTINE HAFTON (Z,L,M,N,FLO,HI,NLEV,NOPT,NPRM,ISPV,SPVAL)
C
C
C     DIMENSION OF           Z(L,M)
C     ARGUMENTS
C
C     LATEST REVISION        JULY,1984
C
C     PURPOSE                HAFTON DRAWS A HALF-TONE PICTURE FROM DATA
C                            STORED IN A RECTANGULAR ARRAY WITH THE
C                            INTENSITY IN THE PICTURE PROPORTIONAL TO
C                            THE DATA VALUE.
C
C     USAGE                  IF THE FOLLOWING ASSUMPTIONS ARE MET, USE
C
C                                 CALL EZHFTN (Z,M,N)
C
C                              ASSUMPTIONS:
C                                  .ALL OF THE ARRAY IS TO BE DRAWN.
C                                  .LOWEST VALUE IN Z WILL BE AT LOWEST
C                                   INTENSITY ON READER/PRINTER OUTPUT.
C                                  .HIGHEST VALUE IN Z WILL BE AT
C                                   HIGHEST INTENSITY.
C                                  .VALUES IN BETWEEN WILL APPEAR
C                                   LINEARLY SPACED.
C                                  .MAXIMUM POSSIBLE NUMBER OF
C                                   INTENSITIES ARE USED.
C                                  .THE PICTURE WILL HAVE A PERIMETER
C                                   DRAWN.
C                                  .FRAME WILL BE CALLED AFTER THE
C                                   PICTURE IS DRAWN.
C                                  .Z IS FILLED WITH NUMBERS THAT SHOULD
C                                   BE USED (NO MISSING VALUES).
C
C                            IF THESE ASSUMPTIONS ARE NOT MET, USE
C
C                                 CALL HAFTON (Z,L,M,N,FLO,HI,NLEV,
C                                              NOPT,NPRM,ISPV,SPVAL)
C
C     ARGUMENTS
C
C     ON INPUT               Z
C     FOR EZHFTN               M BY N ARRAY TO BE USED TO GENERATE A
C                              HALF-TONE PLOT.
C
C                            M
C                              FIRST DIMENSION OF Z.
C
C                            N
C                              SECOND DIMENSION OF Z.
C
C     ON OUTPUT              ALL ARGUMENTS ARE UNCHANGED.
C     FOR EZHFTN
C
C     ON INPUT               Z
C     FOR HAFTON               THE ORIGIN OF THE ARRAY TO BE PLOTTED.
C
C                            L
C                              THE FIRST DIMENSION OF Z IN THE CALLING
C                              PROGRAM.
C
C                            M
C                              THE NUMBER OF DATA VALUES TO BE PLOTTED
C                              IN THE X-DIRECTION (THE FIRST SUBSCRIPT
C                              DIRECTION).  WHEN PLOTTING ALL OF AN
C                              ARRAY, L = M.
C
C                            N
C                              THE NUMBER OF DATA VALUES TO BE PLOTTED
C                              IN THE Y-DIRECTION (THE SECOND SUBSCRIPT
C                              DIRECTION).
C
C                            FLO
C                              THE VALUE OF Z THAT CORRESPONDS TO THE
C                              LOWEST INTENSITY.  (WHEN NOPT.LT.0, FLO
C                              CORRESPONDS TO THE HIGHEST INTENSITY.)
C                              IF FLO=HI=0.0, MIN(Z) WILL BE USED FOR FLO.
C
C                            HI
C                              THE VALUE OF Z THAT CORRESPONDS TO THE
C                              HIGHEST INTENSITY.  (WHEN NOPT.LT.0, HI
C                              CORRESPONDS TO THE LOWEST INTENSITY.)  IF
C                              HI=FLO=0.0, MAX(Z) WILL BE USED FOR HI.
C
C                            NLEV
C                              THE NUMBER OF INTENSITY LEVELS DESIRED.
C                              16 MAXIMUM.  IF NLEV = 0 OR 1, 16 LEVELS
C                              ARE USED.
C
C                            NOPT
C                              FLAG TO CONTROL THE MAPPING OF Z ONTO THE
C                              INTENSITIES.  THE SIGN OF NOPT CONTROLS
C                              THE DIRECTNESS OR INVERSENESS OF THE
C                              MAPPING.
C
C                              . NOPT POSITIVE YIELDS DIRECT MAPPING.
C                                THE LARGEST VALUE OF Z PRODUCES THE
C                                MOST DENSE DOTS.  ON MECHANICAL PLOTTERS,
C                                LARGE VALUES OF Z WILL PRODUCE A DARK
C                                AREA ON THE PAPER.  WITH THE FILM
C                                DEVELOPMENT METHODS USED AT NCAR,
C                                LARGE VALUES OF Z WILL PRODUCE MANY
C                                (WHITE) DOTS ON THE FILM, ALSO
C                                RESULTING IN A DARK AREA ON
C                                READER-PRINTER PAPER.
C                              . NOPT NEGATIVE YIELDS INVERSE MAPPING.
C                                THE SMALLEST VALUES OF Z PRODUCE THE
C                                MOST DENSE DOTS RESULTING IN DARK
C                                AREAS ON THE PAPER.
C
C                              THE ABSOLUTE VALUE OF NOPT DETERMINES THE
C                              MAPPING OF Z ONTO THE INTENSITIES.  FOR
C                              IABS(NOPT)
C                              = 0  THE MAPPING IS LINEAR.  FOR
C                                   EACH INTENSITY THERE IS AN EQUAL
C                                   RANGE IN Z VALUE.
C                              = 1  THE MAPPING IS LINEAR.  FOR
C                                   EACH INTENSITY THERE IS AN EQUAL
C                                   RANGE IN Z VALUE.
C                              = 2  THE MAPPING IS EXPONENTIAL.  FOR
C                                   LARGER VALUES OF Z, THERE IS A
C                                   LARGER DIFFERENCE IN INTENSITY FOR
C                                   RELATIVELY CLOSE VALUES OF Z.  DETAILS
C                                   IN THE LARGER VALUES OF Z ARE DISPLAYED
C                                   AT THE EXPENSE OF THE SMALLER VALUES
C                                   OF Z.
C                              = 3  THE MAPPING IS LOGRITHMIC, SO
C                                   DETAILS OF SMALLER VALUES OF Z ARE SHOWN
C                                   AT THE EXPENSE OF LARGER VALUES OF Z.
C                              = 4  SINUSOIDAL MAPPING, SO MID-RANGE VALUES
C                                   OF Z SHOW DETAILS AT THE EXPENSE OF
C                                   EXTREME VALUES OF Z.
C                              = 5  ARCSINE MAPPING, SO EXTREME VALUES OF
C                                   Z ARE SHOWN AT THE EXPENSE OF MID-RANGE
C                                   VALUES OF Z.
C
C                            NPRM
C                              FLAG TO CONTROL THE DRAWING OF A
C                              PERIMETER AROUND THE HALF-TONE PICTURE.
C
C                              . NPRM=0:  THE PERIMETER IS DRAWN WITH
C                                TICKS POINTING AT DATA LOCATIONS.
C                                (SIDE LENGTHS ARE PROPORTIONAL TO NUMBER
C                                OF DATA VALUES.)
C                              . NPRM POSITIVE:  NO PERIMETER IS DRAWN.  THE
C                                PICTURE FILLS THE FRAME.
C                              . NPRM NEGATIVE:  THE PICTURE IS WITHIN THE
C                                CONFINES OF THE USER'S CURRENT VIEWPORT
C                                SETTING.
C
C                            ISPV
C                              FLAG TO TELL IF THE SPECIAL VALUE FEATURE
C                              IS BEING USED.  THE SPECIAL VALUE FEATURE
C                              IS USED TO MARK AREAS WHERE THE DATA IS
C                              NOT KNOWN OR HOLES ARE WANTED IN THE
C                              PICTURE.
C
C                              . ISPV = 0:  SPECIAL VALUE FEATURE NOT IN
C                                USE.  SPVAL IS IGNORED.
C                              . ISPV NON-ZERO:  SPECIAL VALUE FEATURE
C                                IN USE.  SPVAL DEFINES THE SPECIAL
C                                VALUE.  WHERE Z CONTAINS THE SPECIAL
C                                VALUE, NO HALF-TONE IS DRAWN.  IF ISPV
C                                = 0  SPECIAL VALUE FEATURE NOT IN USE.
C                                     SPVAL IS IGNORED.
C                                = 1  NOTHING IS DRAWN IN SPECIAL VALUE
C                                     AREA.
C                                = 2  CONTIGUOUS SPECIAL VALUE AREAS ARE
C                                     SURROUNDED BY A POLYGONAL LINE.
C                                = 3  SPECIAL VALUE AREAS ARE FILLED
C                                     WITH X(S).
C                                = 4  SPECIAL VALUE AREAS ARE FILLED IN
C                                     WITH THE HIGHEST INTENSITY.
C
C                            SPVAL
C                              THE VALUE USED IN Z TO DENOTE MISSING
C                              VALUES.  THIS ARGUMENT IS IGNORED IF
C                              ISPV = 0.
C
C     ON OUTPUT              ALL ARGUMENTS ARE UNCHANGED.
C     FOR HAFTON
C
C     NOTE                   THIS ROUTINE PRODUCES A HUGE NUMBER OF
C                            PLOTTER INSTRUCTIONS PER PICTURE, AVERAGING
C                            OVER 100,000 LINE-DRAWS PER FRAME WHEN M = N.
C
C
C     ENTRY POINTS           EZHFTN, HAFTON, ZLSET, GRAY, BOUND, HFINIT
C
C     COMMON BLOCKS          HAFT01, HAFT02, HAFT03, HAFT04
C
C     REQUIRED LIBRARY       GRIDAL, THE ERPRT77 PACKAGE AND THE SPPS.
C     ROUTINES
C
C     I/O                    PLOTS HALF-TONE PICTURE.
C
C     PRECISION              SINGLE
C
C     LANGUAGE               FORTRAN
C
C     HISTORY                REWRITE OF PHOMAP ORIGINALLY WRITTEN BY
C                            M. PERRY OF HIGH ALTITUDE OBSERVATORY,
C                            NCAR.
C
C     ALGORITHM              BI-LINEAR INTERPOLATION ON PLOTTER
C                            (RESOLUTION-LIMITED) GRID OF NORMALIZED
C                            REPRESENTATION OF DATA.
C
C     PORTABILITY            ANSI FORTRAN 77.
C
C
C
C INTERNAL PARAMTERSS
C VALUES SET IN BLOCK DATA
C         NAME  DEFAULT                   FUNCTION
C         ----  -------                   ________
C
C         XLT     0.1    LEFT-HAND EDGE OF THE PLOT WHEN NSET=0.  (0.0=
C                        LEFT EDGE OF FRAME, 1.0=RIGHT EDGE OF FRAME.)
C         YBT     0.1    BOTTOM EDGE OF THE PLOT WHEN NSET=0.  (0.0=
C                        BOTTOM OF FRAME, 1.0=TOP OF FRAME.)
C         SIDE    0.8    LENGTH OF LONGER EDGE OF PLOT (SEE ALSO EXT).
C         EXT     .25    LENGTHS OF THE SIDES OF THE PLOT ARE PROPOR-
C                        TIONAL TO M AND N (WHEN NSET=0) EXCEPT IN
C                        EXTREME CASES, NAMELY, WHEN MIN(M,N)/MAX(M,N)
C                        IS LESS THAN EXT.  THEN A SQUARE PLOT IS PRO-
C                        DUCED.  WHEN A RECTANGULAR PLOT IS PRODUCED,
C                        THE PLOT IS CENTERED ON THE FRAME (AS LONG AS
C                        SIDE+2*XLT = SIDE+2*YBT=1., AS WITH THE
C                        DEFAULTS.)
C         ALPHA   1.6    A PARAMETER TO CONTROL THE EXTREMENESS OF THE
C                        MAPPING FUNCTION SPECIFIED BY NOPT.  (FOR
C                        IABS(NOPT)=0 OR 1, THE MAPPING FUNCTION IS
C                        LINEAR AND INDEPENDENT OF ALPHA.)  FOR THE NON-
C                        LINEAR MAPPING FUNCTIONS, WHEN ALPHA IS CHANGED
C                        TO A NUMBER CLOSER TO 1., THE MAPPING FUNCTION
C                        BECOMES MORE LINEAR; WHEN ALPHA IS CHANGED TO
C                        A LARGER NUMBER, THE MAPPING FUNCTION BECOMES
C                        MORE EXTREME.
C         MXLEV   16     MAXIMUM NUMBER OF LEVELS.  LIMITED BY PLOTTER.
C         NCRTG   8      NUMBER OF CRT UNITS PER GRAY-SCALE CELL.
C                        LIMITED BY PLOTTER.
C         NCRTF   1024   NUMBER OF PLOTTER ADDRESS UNITS PER FRAME.
C         IL    (BELOW)  AN ARRAY DEFINING WHICH OF THE AVAILABLE IN-
C                        TENSITIES ARE USED WHEN LESS THAN THE MAXIMUM
C                        NUMBER OF INTENSITIES ARE REQUESTED.
C
C
C   NLEV         INTENSITIES USED
C   ____         ________________
C    2           5,11,
C    3           4, 8,12,
C    4           3, 6,10,13,
C    5           2, 5, 8,11,14,
C    6           1, 4, 7, 9,12,15,
C    7           1, 4, 6, 8,10,12,15,
C    8           1, 3, 5, 7, 9,11,13,15,
C    9           1, 3, 4, 6, 8,10,12,13,15
C   10           1, 3, 4, 6, 7, 9,10,12,13,15,
C   11           1, 2, 3, 5, 6, 8,10,11,13,14,15,
C   12           1, 2, 3, 5, 6, 7, 9,10,11,13,14,15,
C   13           1, 2, 3, 4, 6, 7, 8, 9,10,12,13,14,15
C   14           1, 2, 3, 4, 5, 6, 7, 9,10,11,12,13,14,15,
C   15           1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13,14,15,
C   16           0, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13,14,15
C
C

      SAVE
      DIMENSION       Z(L,N)      ,PX(2)       ,PY(2)
      DIMENSION       ZLEV(16)    ,VWPRT(4)    ,WNDW(4)
      DIMENSION                    VWPR2(4)    ,WND2(4)
      CHARACTER*11    IDUMMY
C
C
      COMMON /HAFTO1/ I          ,J          ,INTEN
      COMMON /HAFTO2/ GLO        ,HA         ,NOPTN      ,ALPHA      ,
     1                NSPV       ,SP         ,ICNST
      COMMON /HAFTO3/ XLT        ,YBT        ,SIDE       ,EXT        ,
     1                IOFFM      ,ALPH       ,MXLEV      ,NCRTG      ,
     2                NCRTF      ,IL(135)
      COMMON /HAFTO4/ NPTMAX     ,NPOINT     ,XPNT(50)  ,YPNT(50)
C +NOAO - Blockdata rewritten as run time initialization subroutine
C
C     EXTERNAL HFINIT
      call hfinit
C -NOAO
C
C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR
C
      CALL Q8QST4 ('GRAPHX','HAFTON','HAFTON','VERSION  1')
C
      NPOINT = 0
      ALPHA = ALPH
      GLO = FLO
      HA = HI
      NLEVL = MIN0(IABS(NLEV),MXLEV)
      IF (NLEVL .LE. 1) NLEVL = MXLEV
      NOPTN = NOPT
      IF (NOPTN .EQ. 0) NOPTN = 1
      NPRIM = NPRM
      NSPV = MAX0(MIN0(ISPV,4),0)
      IF (NSPV .NE. 0) SP = SPVAL
      MX = L
      NX = M
      NY = N
      CRTF = NCRTF
      MSPV = 0
C
C SET INTENSITY BOUNDARY LEVELS
C
      CALL ZLSET (Z,MX,NX,NY,ZLEV,NLEVL)
C
C SET UP PERIMETER
C
      X3 = NX
      Y3 = NY
      CALL GQCNTN (IERR,NTORIG)
      CALL GETUSV('LS',IOLLS)
      IF (NPRIM.LT.0) THEN
         CALL GQNT (NTORIG,IERR,WNDW,VWPRT)
         X1 = VWPRT(1)
         X2 = VWPRT(2)
         Y1 = VWPRT(3)
         Y2 = VWPRT(4)
      ELSE IF (NPRIM.EQ.0) THEN
         X1 = XLT
         X2 = XLT+SIDE
         Y1 = YBT
         Y2 = YBT+SIDE
         IF (AMIN1(X3,Y3)/AMAX1(X3,Y3) .GE. EXT) THEN
            IF (NX-NY.LT.0) THEN
               X2 =SIDE*X3/Y3+XLT
               X2 = (AINT(X2*CRTF/FLOAT(NCRTG))*FLOAT(NCRTG))/CRTF
            ELSE IF (NX-NY.GT.0) THEN
               Y2 = SIDE*Y3/X3+YBT
               Y2 = (AINT(Y2*CRTF/FLOAT(NCRTG))*FLOAT(NCRTG))/CRTF
            END IF
         END IF
      ELSE IF (NPRIM.GT.0) THEN
         X1 = 0.0
         X2 = 1.0
         Y1 = 0.0
         Y2 = 1.0
      END IF
      MX1 = X1*CRTF
      MX2 = X2*CRTF
      MY1 = Y1*CRTF
      MY2 = Y2*CRTF
      IF (NPRIM.GT.0) THEN
         MX1 = 1
         MY1 = 1
         MX2 = NCRTF
         MY2 = NCRTF
      END IF
C
C SAVE NORMALIZATION TRANS 1
C
      CALL GQNT (1,IERR,WNDW,VWPRT)
C
C DEFINE NORMALIZATION TRANS 1 AND LOG SCALING FOR USE WITH PERIM
C DRAW PERIMETER IF NPRIM EQUALS 0
C
      CALL SET(X1,X2,Y1,Y2,1.0,X3,1.0,Y3,1)
      IF (NPRIM .EQ. 0) CALL PERIM (NX-1,1,NY-1,1)
      IF (ICNST .NE. 0) THEN
        CALL GSELNT (0)
        CALL WTSTR(XLT*1.1,0.5,'CONSTANT FIELD',2,0,0)
        GO TO 132
      END IF
C
C FIND OFFSET FOR REFERENCE TO IL, WHICH IS TRIANGULAR
C
      IOFFST = NLEVL*((NLEVL-1)/2)+MOD(NLEVL-1,2)*(NLEVL/2)-1
C
C OUTPUT INTENSITY SCALE
C
      IF (NPRIM .GT. 0) GO TO 112
      LEV = 0
      KX = (1.1*XLT+SIDE)*CRTF
      KY = YBT*CRTF
      NNX = KX/NCRTG
  109 LEV = LEV+1
C +NOAO
C The following statement moved from after statement label 111 (CONTINUE) to
C here.  Otherwise an extra (unlabelled) grayscale box was being drawn.
C This was (eventually) causing a [floating operand error] on a Sun-3.
      IF (LEV .GT. NLEVL) GO TO 112
C -NOAO
      ISUB = IOFFST+LEV
      INTEN = IL(ISUB)
      IF (NOPTN .LT. 0) INTEN = MXLEV-INTEN
      NNY = KY/NCRTG
      DO 111 JJ=1,3
         DO 110 II=1,10
            I = NNX+II
            J = NNY+JJ
            CALL GRAY
  110    CONTINUE
  111 CONTINUE
C +NOAO - FTN internal write rewritten as call to encode.
C       WRITE(IDUMMY,'(G11.4)') ZLEV(LEV)
	call encode (11, '(g11.4)', idummy, zlev(lev))
C -NOAO
      TKX = KX
      TKY = KY+38
      CALL GQNT(1,IERR,WND2,VWPR2)
      CALL SET(0.,1.,0.,1.,0.,1023.,0.,1023.,1)
      CALL WTSTR (TKX,TKY,IDUMMY,0,0,-1)
      CALL SET(VWPR2(1),VWPR2(2),VWPR2(3),VWPR2(4),
     -         WND2(1),WND2(2),WND2(3),WND2(4),1)
C
C ADJUST 38 TO PLOTTER.
C
      KY = KY+52
C
C ADJUST 52 TO PLOTTER.
C
      GO TO 109
C
C STEP THROUGH PLOTTER GRID OF INTENSITY CELLS.
C
  112 IMIN = (MX1-1)/NCRTG+1
      IMAX = (MX2-1)/NCRTG
      JMIN = (MY1-1)/NCRTG+1
      JMAX = (MY2-1)/NCRTG
      XL = IMAX-IMIN+1
      YL = JMAX-JMIN+1
      XN = NX
      YN = NY
      LSRT = NLEVL/2
      DO 130 J=JMIN,JMAX
C
C FIND Y FOR THIS J AND Z FOR THIS Y.
C
       YJ = (FLOAT(J-JMIN)+.5)/YL*(YN-1.)+1.
         LOWY = YJ
         YPART = YJ-FLOAT(LOWY)
         IF (LOWY .NE. NY) GO TO 113
         LOWY = LOWY-1
         YPART = 1.
  113    IPEN = 0
         ZLFT = Z(1,LOWY)+YPART*(Z(1,LOWY+1)-Z(1,LOWY))
         ZRHT = Z(2,LOWY)+YPART*(Z(2,LOWY+1)-Z(2,LOWY))
         IF (NSPV .EQ. 0) GO TO 114
         IF (Z(1,LOWY).EQ.SP .OR. Z(2,LOWY).EQ.SP .OR.
     1       Z(1,LOWY+1).EQ.SP .OR. Z(2,LOWY+1).EQ.SP) IPEN = 1
  114    IF (IPEN .EQ. 1) GO TO 117
C
C FIND INT FOR THIS Z.
C
         IF (ZLFT .GT. ZLEV(LSRT+1)) GO TO 116
  115    IF (ZLFT .GE. ZLEV(LSRT)) GO TO 117
C
C LOOK LOWER
C
         IF (LSRT .LE. 1) GO TO 117
         LSRT = LSRT-1
         GO TO 115
C
C LOOK HIGHER
C
  116    IF (LSRT .GE. NLEVL) GO TO 117
         LSRT = LSRT+1
         IF (ZLFT .GT. ZLEV(LSRT+1)) GO TO 116
C
C OK
C
  117    IRHT = 2
         LAST = LSRT
         DO 129 I=IMIN,IMAX
C
C FIND X FOR THIS I AND Z FOR THIS X AND Y.
C
            IADD = 1
            XI = (FLOAT(I-IMIN)+.5)/XL*(XN-1.)+1.
            LOWX = XI
            XPART = XI-FLOAT(LOWX)
            IF (LOWX .NE. NX) GO TO 118
            LOWX = LOWX-1
            XPART = 1.
C
C TEST FOR INTERPOLATION POSITIONING
C
  118       IF (LOWX .LT. IRHT) GO TO 119
C
C MOVE INTERPOLATION ONE CELL TO THE RIGHT
C
            ZLFT = ZRHT
            IRHT = IRHT+1
            ZRHT = Z(IRHT,LOWY)+YPART*(Z(IRHT,LOWY+1)-Z(IRHT,LOWY))
            IF (NSPV .EQ. 0) GO TO 118
            IPEN = 0
            IF (Z(IRHT-1,LOWY).EQ.SP .OR. Z(IRHT,LOWY).EQ.SP .OR.
     1          Z(IRHT-1,LOWY+1).EQ.SP .OR. Z(IRHT,LOWY+1).EQ.SP)
     2          IPEN = 1
            GO TO 118
  119       IF (IPEN .NE. 1) GO TO 123
C
C SPECIAL VALUE AREA
C
            GO TO (129,120,121,122),NSPV
  120       MSPV = 1
            GO TO 129
  121       PX(1) = I*NCRTG
            PY(1) = J*NCRTG
            PX(2) = PX(1)+NCRTG-1
            PY(2) = PY(1)+NCRTG-1
            CALL GPL (2,PX,PY)
            PYTMP = PY(1)
            PY(1) = PY(2)
            PY(2) = PYTMP
            CALL GPL (2,PX,PY)
C
            GO TO 129
  122       INTEN = MXLEV
            GO TO 128
  123       ZZ = ZLFT+XPART*(ZRHT-ZLFT)
C
C TEST FOR SAME INT AS LAST TIME.
C
            IF (ZZ .GT. ZLEV(LAST+1)) GO TO 126
  124       IF (ZZ .GE. ZLEV(LAST)) GO TO 127
C
C LOOK LOWER
C
            IF (LAST .LE. 1) GO TO 125
            LAST = LAST-1
            GO TO 124
  125       IF (ZZ .LT. ZLEV(LAST)) IADD = 0
            GO TO 127
C
C LOOK HIGHER
C
  126       IF (LAST .GE. NLEVL) GO TO 127
            LAST = LAST+1
            IF (ZZ .GE. ZLEV(LAST+1)) GO TO 126
C
C OK
C
  127       ISUB = LAST+IOFFST+IADD
            INTEN = IL(ISUB)
            IF (NOPTN .LT. 0) INTEN = MXLEV-INTEN
  128       CALL GRAY
  129    CONTINUE
  130 CONTINUE
C
C PUT OUT ANY REMAINING BUFFERED POINTS.
C
      IF (NPOINT.GT.0) THEN
         CALL GQNT(1,IERR,WND2,VWPR2)
         CALL SET(0.,1.,0.,1.,0.,1023.,0.,1023.,1)
         CALL POINTS(XPNT,YPNT,NPOINT,0,0)
         CALL SET(VWPR2(1),VWPR2(2),VWPR2(3),VWPR2(4),
     -            WND2(1),WND2(2),WND2(3),WND2(4),1)
      ENDIF
C
C CALL BOUND IF ISPV=2 AND SPECIAL VALUES WERE FOUND.
C
      IF (MSPV .EQ. 1) THEN
        CALL SET(X1,X2,Y1,Y2,1.0,X3,1.0,Y3,1)
        CALL BOUND (Z,MX,NX,NY,SP)
      END IF
  132 CONTINUE
C
C RESTORE NORMALIZATION TRANS 1 AND ORIGINAL NORMALIZATION NUMBER
C
      CALL SET(VWPRT(1),VWPRT(2),VWPRT(3),VWPRT(4),
     -         WNDW(1),WNDW(2),WNDW(3),WNDW(4),IOLLS)
      CALL SETUSV('LS',IOLLS)
      CALL GSELNT (NTORIG)
      RETURN
C
      END
      SUBROUTINE ZLSET (Z,MX,NX,NY,ZL,NLEVL)
      SAVE
C
      DIMENSION       Z(MX,NY)   ,ZL(NLEVL)
C
      COMMON /HAFTO2/ GLO        ,HA         ,NOPTN      ,ALPHA      ,
     1                NSPV       ,SP         ,ICNST
C
      BIG = R1MACH(2)
C
C ZLSET PUTS THE INTENSITY LEVEL BREAK POINTS IN ZL.
C ALL ARGUMENTS ARE AS IN HAFTON.
C
      LX = NX
      LY = NY
      NLEV = NLEVL
      NOPT = IABS(NOPTN)
      RALPH = 1./ALPHA
      ICNST = 0
      IF (GLO.NE.0. .OR. HA.NE.0.) GO TO 106
C
C FIND RANGE IF NOT KNOWN.
C
      GLO = BIG
      HA = -GLO
      IF (NSPV .NE. 0) GO TO 103
      DO 102 J=1,LY
         DO 101 I=1,LX
            ZZ = Z(I,J)
            GLO = AMIN1(ZZ,GLO)
            HA = AMAX1(ZZ,HA)
  101    CONTINUE
  102 CONTINUE
      GO TO 106
  103 DO 105 J=1,LY
         DO 104 I=1,LX
            ZZ = Z(I,J)
            IF (ZZ .EQ. SP) GO TO 104
            GLO = AMIN1(ZZ,GLO)
            HA = AMAX1(ZZ,HA)
  104    CONTINUE
  105 CONTINUE
C
C FILL ZL
C
  106 DELZ = HA-GLO
      IF (DELZ .EQ. 0.) GO TO 115
      DZ = DELZ/FLOAT(NLEV)
      NLEVM1 = NLEV-1
      DO 114 K=1,NLEVM1
         ZNORM = FLOAT(K)/FLOAT(NLEV)
         GO TO (107,108,109,110,111),NOPT
C
C NOPT=1
C
  107    ZL(K) = GLO+FLOAT(K)*DZ
         GO TO 114
C
C NOPT=2
C
  108    ONORM = (1.-(1.-ZNORM)**ALPHA)**RALPH
         GO TO 113
C
C NOPT=3
C
  109    ONORM = 1.-(1.-ZNORM**ALPHA)**RALPH
         GO TO 113
C
C NOPT=4
C
  110    ONORM = .5*(1.-(ABS(ZNORM+ZNORM-1.))**ALPHA)**RALPH
         GO TO 112
C
C NOPT=5
C
  111    ZNORM2 = ZNORM+ZNORM
         IF (ZNORM .GT. .5) ZNORM2 = 2.-ZNORM2
         ONORM = .5*(1.-(1.-ABS(ZNORM2)**ALPHA)**RALPH)
  112    IF (ZNORM .GT. .5) ONORM = 1.-ONORM
  113    ZL(K) = GLO+DELZ*ONORM
  114 CONTINUE
      ZL(NLEV) = BIG
      RETURN
  115 ICNST = 1
      RETURN
      END
      SUBROUTINE GRAY
C
C SUBROUTINE GRAY COLORS HALF-TONE CELL (I,J) WITH INTENSITY INTEN.
C THE ROUTINE ASSUMES 8X8 CELL SIZE ON A VIRTUAL SCREEN 1024X1024.
C
      DIMENSION       IFOT(16)   ,JFOT(16)
      DIMENSION       WNDW(4)    ,VWPRT(4)
CCC      DIMENSION       MX(16)  ,MY(16)
      COMMON /HAFTO1/ I          ,J          ,INTEN
      COMMON /HAFTO4/ NPTMAX     ,NPOINT     ,XPNT(50)  ,YPNT(50)
      SAVE
C
      DATA
     1 IFOT(1),IFOT(2),IFOT(3),IFOT(4),IFOT(5),IFOT(6),IFOT(7),IFOT(8)/
     2   1,      5,      1,      5,      3,      7,      3,      7   /
      DATA
     1 IFOT(9),IFOT(10),IFOT(11),IFOT(12),IFOT(13),IFOT(14),IFOT(15)/
     2   3,      7,       3,       7,       1,       5,       1/,
     3 IFOT(16)/
     4   5    /
C
      DATA
     1 JFOT(1),JFOT(2),JFOT(3),JFOT(4),JFOT(5),JFOT(6),JFOT(7),JFOT(8)/
     2   1,      5,      5,      1,      3,     7,      7,      3   /
      DATA
     1 JFOT(9),JFOT(10),JFOT(11),JFOT(12),JFOT(13),JFOT(14),JFOT(15)/
     2   1,      5,       5,       1,       3,       7,       7/,
     3 JFOT(16)/
     4   3    /
C
      IF (INTEN) 103,103,101
  101 I1 = I*8
      J1 = J*8
      IF ((NPOINT+INTEN) .LE.NPTMAX)  GO TO 1015
      CALL GQNT(1,IERR,WNDW,VWPRT)
      CALL SET(0.,1.,0.,1.,0.,1023.,0.,1023.,1)
      CALL POINTS(XPNT,YPNT,NPOINT,0,0)
      CALL SET(VWPRT(1),VWPRT(2),VWPRT(3),VWPRT(4),
     -         WNDW(1),WNDW(2),WNDW(3),WNDW(4),1)
      NPOINT = 0
 1015 DO 102 I2=1,INTEN
      NPOINT = NPOINT + 1
      XPNT(NPOINT) = I1+IFOT(I2)
      YPNT(NPOINT) = J1+JFOT(I2)
  102 CONTINUE
  103 RETURN
      END
      SUBROUTINE BOUND (Z,MX,NNX,NNY,SSP)
      DIMENSION       Z(MX,NNY)   ,PX(2)   ,PY(2)
C
C BOUND DRAWS A POLYGONAL BOUNDRY AROUND ANY SPECIAL-VALUE AREAS IN Z.
C
      SAVE
      NX = NNX
      NY = NNY
C
C VERTICAL LINES
C
      SP = SSP
      DO 103 IP1=3,NX
         I = IP1-1
         PX(1) = I
         PX(2) = I
         IM1 = I-1
         DO 102 JP1=2,NY
            PY(2) = JP1
            J = JP1-1
            PY(1) = J
            KLEFT = 0
            IF (Z(IM1,J).EQ.SP .OR. Z(IM1,JP1).EQ.SP) KLEFT = 1
            KCENT = 0
            IF (Z(I,J).EQ.SP .OR. Z(I,JP1).EQ.SP) KCENT = 1
            KRIGT = 0
            IF (Z(IP1,J).EQ.SP .OR. Z(IP1,JP1).EQ.SP) KRIGT = 1
            JUMP = KLEFT*4+KCENT*2+KRIGT+1
            GO TO (102,101,102,102,101,102,102,102,102),JUMP
  101       CALL GPL (2,PX,PY)
  102    CONTINUE
  103 CONTINUE
C
C HORIZONTAL
C
      DO 106 JP1=3,NY
         J = JP1-1
         PY(1) = J
         PY(2) = J
         JM1 = J-1
         DO 105 IP1=2,NX
            PX(2) = IP1
            I = IP1-1
            PX(1) = I
            KLOWR = 0
            IF (Z(I,JM1).EQ.SP .OR. Z(IP1,JM1).EQ.SP) KLOWR = 1
            KCENT = 0
            IF (Z(I,J).EQ.SP .OR. Z(IP1,J).EQ.SP) KCENT = 1
            KUPER = 0
            IF (Z(I,JP1).EQ.SP .OR. Z(IP1,JP1).EQ.SP) KUPER = 1
            JUMP = KLOWR*4+KCENT*2+KUPER+1
            GO TO (105,104,105,105,104,105,105,105,105),JUMP
  104       CALL GPL (2,PX,PY)
  105    CONTINUE
  106 CONTINUE
      RETURN
      END
      SUBROUTINE EZHFTN (Z,M,N)
C
      DIMENSION       Z(M,N)
      SAVE
C
C HALF-TONE PICTURE VIA SHORTEST ARGUMENT LIST.
C ASSUMPTIONS--
C     ALL OF THE ARRAY IS TO BE DRAWN,
C     LOWEST VALUE IN Z WILL BE AT LOWEST INTENSITY ON READER/PRINTER
C     OUTPUT, HIGHEST VALUE IN Z WILL BE AT HIGHEST INTENSITY, VALUES IN
C     BETWEEN WILL APPEAR LINEARLY SPACED, MAXIMUM POSSIBLE NUMBER OF
C     INTENSITIES ARE USED, THE PICTURE WILL HAVE A PERIMETER DRAWN,
C     FRAME WILL BE CALLED AFTER THE PICTURE IS DRAWN, Z IS FILLED WITH
C     NUMBERS THAT SHOULD BE USED (NO UNKNOWN VALUES).
C IF THESE CONDITIONS ARE NOT MET, USE HAFTON.
C EZHFTN ARGUMENTS--
C     Z   2 DIMENSIONAL ARRAY TO BE USED TO GENERATE A HALF-TONE PLOT.
C     M   FIRST DIMENSION OF Z.
C     N   SECOND DIMENSION OF Z.
C
      DATA FLO,HI,NLEV,NOPT,NPRM,ISPV,SPV/0.0,0.0,0,0,0,0,0.0/
C
C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR
C
      CALL Q8QST4 ('GRAPHX','HAFTON','EZHFTN','VERSION  1')
C
      CALL HAFTON (Z,M,M,N,FLO,HI,NLEV,NOPT,NPRM,ISPV,SPV)
C
C +NOAO - EZHFTN no longer calls frame.
C     CALL FRAME
C -NOAO
      RETURN
      END
C
C-----------------------------------------------------------------------
C
C REVISION HISTORY---
C
C JULY     1984    CONVERTED TO FORTAN 77 AND GKS
C
C MARCH    1983    INSTITUTED BUFFERING OF POINTS WITHIN ROUTINE GRAY,
C                  WHICH DRAMATICALLY REDUCES SIZE OF OUTPUT PLOT CODE,
C                  METACODE.  THIS IN TURN GENERALLY IMPROVES THROUGHPUT
C                  OF METACODE INTERPRETERS.
C
C FEBRUARY 1979    MODIFIED CODE TO CONFORM TO FORTRAN 66 STANDARD
C
C JANUARY  1978    DELETED REFERENCES TO THE  *COSY  CARDS AND
C                  ADDED REVISION HISTORY
C
C-----------------------------------------------------------------------
C

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


!
!MIT License
!
!Copyright (c) 2020 flavianowilliams
!
!Permission is hereby granted, free of charge, to any person obtaining a copy
!of this software and associated documentation files (the "Software"), to deal
!in the Software without restriction, including without limitation the rights
!to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
!copies of the Software, and to permit persons to whom the Software is
!furnished to do so, subject to the following conditions:
!
!THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
!IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
!FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
!AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
!LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
!OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
!SOFTWARE.
!
module interaction_module
  !*******************************************************************************
  !*******************************************************************************

  use neighbourlist_module
  use bonds_module
  use angles_module
  use coulomb_module
  use dihedrals_module
  use vanderwaals_module

  implicit none

  private
  public :: interaction

  type, extends(neighbourlist) :: interaction
     type(bonds)                   :: bnd
     type(angles)                  :: tht
     type(coulomb)                 :: coul
     type(dihedrals)               :: dih
     type(vanderwaals)             :: vdw
     real(8), private              :: eintra
     real(8), private              :: einter
     real(8), private              :: enpot
     real(8), private              :: virtot
     real(8), private              :: encorr
     real(8), private              :: vircorr
     real(8), private              :: fbi(3)
     real(8), private              :: fbj(3)
     real(8), private              :: fbk(3)
     real(8), private              :: fbl(3)
   contains
     procedure :: interaction_init
     procedure :: set_forcefield
     procedure :: set_intraff
     procedure :: set_interff
     procedure :: set_force2
     procedure :: set_force3
     procedure :: set_force4
     generic   :: set_force => set_force2, set_force3, set_force4
     procedure :: set_enpot
     procedure :: get_enpot
     procedure :: set_eintra
     procedure :: get_eintra
     procedure :: set_einter
     procedure :: get_einter
     procedure :: set_virtot
     procedure :: get_virtot
     procedure :: set_vdwcorr
     procedure :: get_encorr
     procedure :: get_vircorr
  end type interaction

contains

  subroutine interaction_init(this)
    implicit none
    class(interaction), intent(inout) :: this
    if(this%get_coulop().eq.'fscs')call this%coul%coulomb_init&
         (this%get_coulop(),this%get_fscsalpha(),this%get_rcutoff(),this%get_pi())
  end subroutine interaction_init

  subroutine set_forcefield(this)
    implicit none
    class(interaction), intent(inout) :: this
    integer                           :: i
    do i=1,this%get_natom()
       this%fax(i)=0.d0
       this%fay(i)=0.d0
       this%faz(i)=0.d0
    end do
    call this%set_virtot(this%get_vircorr())
    call this%set_intraff()
    call this%set_interff()
    call this%set_enpot()
  end subroutine set_forcefield

  subroutine set_intraff(this)
    implicit none
    class(interaction), intent(inout) :: this
    integer                           :: i,j,k,l,m,n,o,ni,nj,nk,nl,nx
    real(8)                           :: xvz,yvz,zvz,dr,virtot,theta,dr1,dr2
    real(8)                           :: eintra
    real(8)                           :: prm(3),drij(3),drik(3),drjk(3)
    real(8)                           :: drkl(3),ri(3),rj(3)
    real(8)                           :: rk(3),rl(3)
    real(8)                           :: vc1x,vc1y,vc1z,vc2x,vc2y,vc2z,phi
    character(7)                      :: ptrm
    nx=0
    eintra=0.d0
    virtot=0.d0
    do i=1,this%get_nmol()
       do j=1,this%ntmol(i)
          do k=1,this%bondscnt(i)
             ni=nx+this%molbond(i,k,1)
             nj=nx+this%molbond(i,k,2)
             call this%mic(ni,nj,xvz,yvz,zvz)
             dr=sqrt(xvz**2+yvz**2+zvz**2)
             do l=1,2
                prm(l)=this%parbnd(i,k,l)
             end do
             ptrm=this%tbonds(i,k)
             call this%bnd%set_bonds(dr,prm,ptrm)
             call this%set_force(ni,nj,xvz,yvz,zvz,this%bnd%get_force())
             call this%bnd%set_virbond(this%bnd%get_force(),dr)
             eintra=eintra+this%bnd%get_enbond()
             virtot=virtot+this%bnd%get_virbond()
          end do
          do k=1,this%bendscnt(i)
             ni=nx+this%molbend(i,k,2)
             nj=nx+this%molbend(i,k,1)
             nk=nx+this%molbend(i,k,3)
             call this%mic(ni,nj,drij(1),drij(2),drij(3))
             call this%mic(ni,nk,drik(1),drik(2),drik(3))
             dr1=sqrt(drij(1)**2+drij(2)**2+drij(3)**2)
             dr2=sqrt(drik(1)**2+drik(2)**2+drik(3)**2)
             theta=acos((drij(1)*drik(1)+drij(2)*drik(2)+drij(3)*drik(3))/(dr1*dr2))
             do l=1,2
                prm(l)=this%parbend(i,k,l)
             end do
             ptrm=this%tbends(i,k)
             call this%tht%set_angles(theta,prm,ptrm)
             call this%set_force&
                  (ni,nj,nk,drij,drik,dr1,dr2,theta,this%tht%get_force())
             call this%tht%set_virbend(this%fbj,this%fbk,drij,drik)
             eintra=eintra+this%tht%get_enbend()
             virtot=virtot+this%tht%get_virbend()
          end do
          do k=1,this%torscnt(i)
             ni=nx+this%moltors(i,k,1)
             nj=nx+this%moltors(i,k,2)
             nk=nx+this%moltors(i,k,3)
             nl=nx+this%moltors(i,k,4)
             call this%mic(ni,nj,drij(1),drij(2),drij(3))
             call this%mic(nj,nk,drjk(1),drjk(2),drjk(3))
             call this%mic(nk,nl,drkl(1),drkl(2),drkl(3))
             vc1x=drij(2)*drjk(3)-drij(3)*drjk(2)
             vc1y=drij(3)*drjk(1)-drij(1)*drjk(3)
             vc1z=drij(1)*drjk(2)-drij(2)*drjk(1)
             vc2x=drjk(2)*drkl(3)-drjk(3)*drkl(2)
             vc2y=drjk(3)*drkl(1)-drjk(1)*drkl(3)
             vc2z=drjk(1)*drkl(2)-drjk(2)*drkl(1)
             dr1=sqrt(vc1x**2+vc1y**2+vc1z**2) !-|rij x rjk|
             dr2=sqrt(vc2x**2+vc2y**2+vc2z**2) !-|rjk x rkn|
             phi=acos((vc1x*vc2x+vc1y*vc2y+vc1z*vc2z)/(dr1*dr2))
             do l=1,3
                prm(l)=this%partors(i,k,l)
             end do
             ptrm=this%ttors(i,k)
             call this%dih%set_dihedrals(phi,prm,ptrm)
             call this%set_force&
                  (ni,nj,nk,nl,drij,drjk,drkl,dr1,dr2,phi,this%dih%get_force())
             ri(1)=this%xa(ni)
             ri(2)=this%ya(ni)
             ri(3)=this%za(ni)
             rj(1)=this%xa(nj)
             rj(2)=this%ya(nj)
             rj(3)=this%za(nj)
             rk(1)=this%xa(nk)
             rk(2)=this%ya(nk)
             rk(3)=this%za(nk)
             rl(1)=this%xa(nl)
             rl(2)=this%ya(nl)
             rl(3)=this%za(nl)
             call this%dih%set_virtors(this%fbi,this%fbj,this%fbk,this%fbl,ri,rj,rk,rl)
             virtot=virtot+this%dih%get_virtors()
          end do
          nx=nx+this%nxmol(i)
       end do
       nx=0
       if(this%nxmol(i).lt.4)goto 2
       do j=1,this%ntmol(i)
          do k=1,this%nxmol(i)
             ni=nx+k
             do l=k+1,this%nxmol(i)
                nj=nx+l
                do m=1,this%bendscnt(i)
                   do n=1,3
                      do o=n+1,3
                         if(ni.eq.this%molbend(i,m,n).and.nj.eq.this%molbend(i,m,o))goto 1
                         if(ni.eq.this%molbend(i,m,o).and.nj.eq.this%molbend(i,m,n))goto 1
                      end do
                   end do
                end do
                call this%mic(ni,nj,xvz,yvz,zvz)
                dr=sqrt(xvz**2+yvz**2+zvz**2)
                if(abs(this%qat(ni)*this%qat(nj)).gt.1.d-8)then
                   call this%coul%set_coulomb(dr,this%qat(ni),this%qat(nj))
                   call this%set_force&
                        (ni,nj,xvz,yvz,zvz,this%coul%get_force()*this%sf_coul(i))
                   call this%coul%set_vircoul(this%coul%get_force()*this%sf_coul(i),dr)
                   eintra=eintra+this%coul%get_encoul()*this%sf_coul(i)
                   virtot=virtot+this%coul%get_vircoul()
                end if
                do m=1,this%get_nvdw()
                   if(this%tpa(ni).eq.this%spcvdw(m,1).and.&
                        this%tpa(nj).eq.this%spcvdw(m,2).or.&
                        this%tpa(ni).eq.this%spcvdw(m,2).and.&
                        this%tpa(nj).eq.this%spcvdw(m,1))then
                      do n=1,2
                         prm(n)=this%parvdw(m,n)
                      end do
                      ptrm='charmm'
                      call this%vdw%set_vanderwaals(dr,prm,ptrm)
                      call this%set_force&
                           (ni,nj,xvz,yvz,zvz,this%vdw%get_force()*this%sf_vdw(i))
                      call this%vdw%set_virvdw(this%vdw%get_force()*this%sf_vdw(i),dr)
                      eintra=eintra+this%vdw%get_envdw()*this%sf_vdw(i)
                      virtot=virtot+this%vdw%get_virvdw()
                   end if
                end do
1               continue
             end do
          end do
          nx=nx+this%nxmol(i)
       end do
2      continue
    end do
    call this%set_eintra(eintra)
    call this%set_virtot(virtot+this%get_virtot())
  end subroutine set_intraff

  subroutine set_interff(this)
    implicit none
    class(interaction), intent(inout) :: this
    integer                           :: i,j,k,l,ni,nj
    real(8)                           :: xvz,yvz,zvz,dr,virtot,einter
    real(8)                           :: prm(3)
    character(7)                      :: ptrm
    einter=0.d0
    virtot=0.d0
    do i=1,this%get_natom()
       do j=1,this%nlist(i)
          ni=i
          nj=this%ilist(i,j)
          call this%mic(ni,nj,xvz,yvz,zvz)
          dr=sqrt(xvz**2+yvz**2+zvz**2)
          if(abs(this%qat(ni)*this%qat(nj)).gt.1.d-8)then
             call this%coul%set_coulomb(dr,this%qat(ni),this%qat(nj))
             call this%set_force(ni,nj,xvz,yvz,zvz,this%coul%get_force())
             call this%coul%set_vircoul(this%coul%get_force(),dr)
             einter=einter+this%coul%get_encoul()
             virtot=virtot+this%coul%get_vircoul()
          end if
          do k=1,this%get_nvdw()
             if(this%tpa(ni).eq.this%spcvdw(k,1).and.this%tpa(nj).eq.this%spcvdw(k,2).or.&
                  this%tpa(ni).eq.this%spcvdw(k,2).and.this%tpa(nj).eq.this%spcvdw(k,1))then
                do l=1,2
                   prm(l)=this%parvdw(k,l)
                end do
                ptrm=this%tvdw(k)
                call this%vdw%set_vanderwaals(dr,prm,ptrm)
                call this%set_force(ni,nj,xvz,yvz,zvz,this%vdw%get_force())
                call this%vdw%set_virvdw(this%vdw%get_force(),dr)
                einter=einter+this%vdw%get_envdw()
                virtot=virtot+this%vdw%get_virvdw()
             end if
          end do
       end do
    end do
    call this%set_einter(einter)
    call this%set_virtot(virtot+this%get_virtot())
  end subroutine set_interff

  subroutine set_force2(this,ni,nj,xvz,yvz,zvz,fr)
    implicit none
    class(interaction), intent(inout) :: this
    integer, intent(in)               :: ni,nj
    real(8), intent(in)               :: fr,xvz,yvz,zvz
    this%fbi(1)=-fr*xvz
    this%fbi(2)=-fr*yvz
    this%fbi(3)=-fr*zvz
    this%fbj(1)=+fr*xvz
    this%fbj(2)=+fr*yvz
    this%fbj(3)=+fr*zvz
    this%fax(ni)=this%fax(ni)+this%fbi(1)
    this%fay(ni)=this%fay(ni)+this%fbi(2)
    this%faz(ni)=this%faz(ni)+this%fbi(3)
    this%fax(nj)=this%fax(nj)+this%fbj(1)
    this%fay(nj)=this%fay(nj)+this%fbj(2)
    this%faz(nj)=this%faz(nj)+this%fbj(3)
  end subroutine set_force2

  subroutine set_force3(this,i1,i2,i3,drij,drik,dr1,dr2,theta,fa)
    implicit none
    class(interaction), intent(inout) :: this
    integer, intent(in)               :: i1,i2,i3
    integer                           :: ix(3),i,j
    real(8), intent(in)               :: fa,dr1,dr2,theta
    real(8), intent(in)               :: drij(3),drik(3)
    real(8)                           :: derij(3,3)
    ix(1)=i1
    ix(2)=i2
    ix(3)=i3
    do j=1,3
       do i=1,3
          derij(i,j)=(kronij(ix(i),ix(2))-kronij(ix(i),ix(1)))*drik(j)/(dr1*dr2) &
               +(kronij(ix(i),ix(3))-kronij(ix(i),ix(1)))*drij(j)/(dr1*dr2) &
               -cos(theta)*((kronij(ix(i),ix(2))-kronij(ix(i),ix(1)))*drij(j)/dr1**2 &
               +(kronij(ix(i),ix(3))-kronij(ix(i),ix(1)))*drik(j)/dr2**2)
       end do
    end do
    this%fbi(1)=fa*derij(1,1)/sin(theta)
    this%fbi(2)=fa*derij(1,2)/sin(theta)
    this%fbi(3)=fa*derij(1,3)/sin(theta)
    this%fbj(1)=fa*derij(2,1)/sin(theta)
    this%fbj(2)=fa*derij(2,2)/sin(theta)
    this%fbj(3)=fa*derij(2,3)/sin(theta)
    this%fbk(1)=fa*derij(3,1)/sin(theta)
    this%fbk(2)=fa*derij(3,2)/sin(theta)
    this%fbk(3)=fa*derij(3,3)/sin(theta)
    this%fax(i1)=this%fax(i1)+this%fbi(1)
    this%fay(i1)=this%fay(i1)+this%fbi(2)
    this%faz(i1)=this%faz(i1)+this%fbi(3)
    this%fax(i2)=this%fax(i2)+this%fbj(1)
    this%fay(i2)=this%fay(i2)+this%fbj(2)
    this%faz(i2)=this%faz(i2)+this%fbj(3)
    this%fax(i3)=this%fax(i3)+this%fbk(1)
    this%fay(i3)=this%fay(i3)+this%fbk(2)
    this%faz(i3)=this%faz(i3)+this%fbk(3)
  end subroutine set_force3

  subroutine set_force4(this,i1,i2,i3,i4,drij,drjk,drkl,dr1,dr2,phi,fd)
    implicit none
    class(interaction), intent(inout) :: this
    integer, intent(in)               :: i1,i2,i3,i4
    integer                           :: ix(4),i,j
    real(8)                           :: dvc(4,3),phi2
    real(8), intent(in)               :: fd,dr1,dr2,phi
    real(8), intent(in)               :: drij(3),drjk(3),drkl(3)
    phi2=max(phi,1.d-8)
    ix(1)=i1 !-i
    ix(2)=i2 !-j
    ix(3)=i3 !-k
    ix(4)=i4 !-n
    do i=1,4
       do j=1,3
          dvc(i,j)=dfunc1(i1,i2,i3,i4,drij,drjk,drkl,ix(i),j)/(dr1*dr2) &
               -0.5d0*cos(phi2)*(dfunc2(i1,i2,i3,drij,drjk,ix(i),j)/dr1**2 &
               +dfunc2(i2,i3,i4,drjk,drkl,ix(i),j)/dr2**2)
       end do
    end do
    this%fbi(1)=fd*dvc(1,1)/sin(phi2)
    this%fbi(2)=fd*dvc(1,2)/sin(phi2)
    this%fbi(3)=fd*dvc(1,3)/sin(phi2)
    this%fbj(1)=fd*dvc(2,1)/sin(phi2)
    this%fbj(2)=fd*dvc(2,2)/sin(phi2)
    this%fbj(3)=fd*dvc(2,3)/sin(phi2)
    this%fbk(1)=fd*dvc(3,1)/sin(phi2)
    this%fbk(2)=fd*dvc(3,2)/sin(phi2)
    this%fbk(3)=fd*dvc(3,3)/sin(phi2)
    this%fbl(1)=fd*dvc(4,1)/sin(phi2)
    this%fbl(2)=fd*dvc(4,2)/sin(phi2)
    this%fbl(3)=fd*dvc(4,3)/sin(phi2)
    this%fax(i1)=this%fax(i1)+this%fbi(1)
    this%fay(i1)=this%fay(i1)+this%fbi(2)
    this%faz(i1)=this%faz(i1)+this%fbi(3)
    this%fax(i2)=this%fax(i2)+this%fbj(1)
    this%fay(i2)=this%fay(i2)+this%fbj(2)
    this%faz(i2)=this%faz(i2)+this%fbj(3)
    this%fax(i3)=this%fax(i3)+this%fbk(1)
    this%fay(i3)=this%fay(i3)+this%fbk(2)
    this%faz(i3)=this%faz(i3)+this%fbk(3)
    this%fax(i4)=this%fax(i4)+this%fbl(1)
    this%fay(i4)=this%fay(i4)+this%fbl(2)
    this%faz(i4)=this%faz(i4)+this%fbl(3)
  end subroutine set_force4

  subroutine set_enpot(this)
    implicit none
    class(interaction), intent(inout) :: this
    this%enpot=this%get_eintra()+this%get_einter()+this%get_encorr()
  end subroutine set_enpot

  double precision function get_enpot(this)
    implicit none
    class(interaction), intent(in) :: this
    get_enpot=this%enpot
  end function get_enpot

  subroutine set_eintra(this,eintra)
    implicit none
    class(interaction), intent(inout) :: this
    real(8), intent(in)               :: eintra
    this%eintra=eintra
  end subroutine set_eintra

  double precision function get_eintra(this)
    implicit none
    class(interaction), intent(in) :: this
    get_eintra=this%eintra
  end function get_eintra

  subroutine set_einter(this,einter)
    implicit none
    class(interaction), intent(inout) :: this
    real(8), intent(in)               :: einter
    this%einter=einter
  end subroutine set_einter

  double precision function get_einter(this)
    implicit none
    class(interaction), intent(in) :: this
    get_einter=this%einter
  end function get_einter

  subroutine set_virtot(this,virtot)
    implicit none
    class(interaction), intent(inout) :: this
    real(8), intent(in)               :: virtot
    this%virtot=virtot
  end subroutine set_virtot

  double precision function get_virtot(this)
    implicit none
    class(interaction), intent(in) :: this
    get_virtot=this%virtot
  end function get_virtot

  subroutine set_vdwcorr(this)
    implicit none
    class(interaction), intent(inout) :: this
    integer                           :: natp1,natp2,i,j
    real(8)                           :: envdw_corr,virvdw_corr,es,vs
    real(8)                           :: prm(2)
    envdw_corr=0.d0
    virvdw_corr=0.d0
    do i=1,this%get_nvdw()
       select case(this%tvdw(i))
       case('charmm')
          do j=1,2
             prm(j)=this%parvdw(i,j)
          end do
          es=prm(1)*&
               (prm(2)**12-6.d0*(this%get_rcutoff()*prm(2))**6)/(9.d0*this%get_rcutoff()**9)
          vs=12.d0*prm(1)*&
               (prm(2)**12-3.d0*(this%get_rcutoff()*prm(2))**6)/(9.d0*this%get_rcutoff()**9)
       case('lj')
          do j=1,2
             prm(j)=this%parvdw(i,j)
          end do
          es=4.d0*prm(1)*&
               (prm(2)**12-3.d0*(this%get_rcutoff()*prm(2))**6)/(9.d0*this%get_rcutoff()**9)
          vs=24.d0*prm(1)*(2.d0*prm(2)**12-3.d0*(this%get_rcutoff()*prm(2))**6)/&
               (9.d0*this%get_rcutoff()**9)
       end select
       natp1=0
       natp2=0
       do j=1,this%get_natom()
          if(this%tpa(j).eq.this%spcvdw(i,1))natp1=natp1+1
          if(this%tpa(j).eq.this%spcvdw(i,2))natp2=natp2+1
       end do
       envdw_corr=envdw_corr+es*natp1*natp2
       virvdw_corr=virvdw_corr+vs*natp1*natp2
    end do
    this%encorr=2.0d0*this%get_pi()*envdw_corr/this%get_volume()
    this%vircorr=-2.0d0*this%get_pi()*virvdw_corr/this%get_volume()
  end subroutine set_vdwcorr

  double precision function get_encorr(this)
    implicit none
    class(interaction), intent(in) :: this
    get_encorr=this%encorr
  end function get_encorr

  double precision function get_vircorr(this)
    implicit none
    class(interaction), intent(in) :: this
    get_vircorr=this%vircorr
  end function get_vircorr

  integer function kronij(i,j)
    integer, intent(in)               :: i,j
    kronij=int((float(i+j)-abs(i-j))/(float(i+j)+abs(i-j)))
  end function kronij

  double precision function dfunc1(i1,i2,i3,i4,drij,drjk,drkn,i,j)
    implicit none
    integer i,j,i1,i2,i3,i4
    real(8) drij(3),drjk(3),drkn(3)
    dfunc1=drij(j)*acomt(drjk,drjk,j)*(kronij(i,i3)-kronij(i,i4)) &
         +drij(j)*acomt(drjk,drkn,j)*(kronij(i,i3)-kronij(i,i2)) &
         +drjk(j)*acomt(drij,drjk,j)*(kronij(i,i4)-kronij(i,i3)) &
         +drjk(j)*acomt(drjk,drkn,j)*(kronij(i,i2)-kronij(i,i1)) &
         +drkn(j)*acomt(drij,drjk,j)*(kronij(i,i3)-kronij(i,i2)) &
         +drkn(j)*acomt(drjk,drjk,j)*(kronij(i,i1)-kronij(i,i2)) &
         +2.d0*drjk(j)*acomt(drij,drkn,j)*(kronij(i,i2)-kronij(i,i3))
  end function dfunc1

  double precision function dfunc2(i1,i2,i3,dri1,dri2,i,j)
    implicit none
    integer i,j,i1,i2,i3
    real(8) dri1(3),dri2(3)
    dfunc2=dri1(j)*acomt(dri2,dri2,j)*(kronij(i,i2)-kronij(i,i1)) &
         +dri1(j)*acomt(dri1,dri2,j)*(kronij(i,i2)-kronij(i,i3)) &
         +dri2(j)*acomt(dri1,dri1,j)*(kronij(i,i3)-kronij(i,i2)) &
         +dri2(j)*acomt(dri1,dri2,j)*(kronij(i,i1)-kronij(i,i2))
    dfunc2=2.d0*dfunc2
  end function dfunc2

  double precision function acomt(x1,x2,k)
    implicit none
    integer i,k
    real(8) sum,x1(3),x2(3)
    sum=0.d0
    do i=1,3
       sum=sum+(1.d0-kronij(i,k))*x1(i)*x2(i)
    end do
    acomt=sum
  end function acomt

end module interaction_module


text:
      program ecg
      integer v, n, i
      parameter(n=32768)
      real signal(n,5)
      common /result/ power(5)
      real power

!hpf$ processors proc(5)
!hpf$ template t(5)
!hpf$ align with t(i):: signal(*,i), power(i)
!hpf$ distribute t(block) onto proc

      external process
      real process
!fcd$ pure process

      do v = 1, 5
         read *, (signal(i, v), i=1, n)
      enddo

!hpf$ independent
      do v = 1, 5
         power(v) = process(signal(1,v))
      enddo

      call printout
      end

      real function process(sig)
!fcd$ pure
!fcd$ pure fft1
      integer f, n
      parameter(n=32768)
      real sig(n)
      real spectre(n)

      call fft1(signal, spectre)

      p = 0.
      do f = 1, n
         p = p + spectre(f)**2
      enddo
      process = sqrt(p)
      end

      subroutine fft1(s, f)
!fcd$ pure
      integer i, n
      parameter(n=32768)
      real s(n), f(n)

      do i = 1, n
         f(i) = s(i)
      enddo

      end

      subroutine printout
      common /result/ power(5)
      real power, p
!hpf$ processors proc(5)
!hpf$ template t(5)
!hpf$ align with t:: power
!hpf$ distribute t(block) onto proc
      integer v

      p = 0.
!hpf$ independent, reduction(p)
      do v = 1, 5
         p = p + power(v)
      enddo

      print *, p

      end


! { dg-do compile }
! { dg-options "-O0" }
! { dg-require-visibility "" }
! Tests the fix for PR24207 and the problems associated
! with the fix for PR21986. In two cases, use associated
! public symbols were taking on the default private access
! attribute of the local namespace. In the third, a private
! symbol was not available to a namelist in contained
! procedure in the same module.
!
! Based on the example in PR24207.
!
module a
  implicit none
  real b
  type :: mytype
    integer :: c
  end type mytype
end module a
module c
  use a
  implicit none
  public d
  private
  real x
  contains
     subroutine d (arg_t)   ! This would cause an error
        type (mytype) :: arg_t
        namelist /e/ b, x   ! .... as would this.
        read(5,e)
	arg_t%c = 42
     end subroutine d
end module c


! Chemfiles, an efficient IO library for chemistry file formats
! Copyright (C) 2015-2019 Guillaume Fraux -- BSD licence
!
! =========================================================================== !
! !!!! AUTO-GENERATED FILE !!!! Do not edit. See bindgen repository for the
! generating code (https://github.com/chemfiles/bindgen).
! This file contains Fortran 2003 ISO C Binding interface to the C API
! =========================================================================== !

interface
! Function "chfl_topology", at topology.h:25
function c_chfl_topology() bind(C, name="chfl_topology")
    use iso_c_binding

    implicit none
    type(c_ptr) :: c_chfl_topology

end function

! Function "chfl_topology_from_frame", at topology.h:38
function c_chfl_topology_from_frame(frame) bind(C, name="chfl_topology_from_frame")
    use iso_c_binding

    implicit none
    type(c_ptr) :: c_chfl_topology_from_frame
    type(c_ptr), value, intent(in) :: frame
end function

! Function "chfl_topology_copy", at topology.h:48
function c_chfl_topology_copy(topology) bind(C, name="chfl_topology_copy")
    use iso_c_binding

    implicit none
    type(c_ptr) :: c_chfl_topology_copy
    type(c_ptr), value, intent(in) :: topology
end function

! Function "chfl_topology_atoms_count", at topology.h:56
function c_chfl_topology_atoms_count(topology, count) bind(C, name="chfl_topology_atoms_count")
    use iso_c_binding
    import chfl_status

    implicit none
    integer(kind=chfl_status) :: c_chfl_topology_atoms_count
    type(c_ptr), value, intent(in) :: topology
    integer(kind=c_int64_t), intent(inout) :: count
end function

! Function "chfl_topology_resize", at topology.h:68
function c_chfl_topology_resize(topology, natoms) bind(C, name="chfl_topology_resize")
    use iso_c_binding
    import chfl_status

    implicit none
    integer(kind=chfl_status) :: c_chfl_topology_resize
    type(c_ptr), value :: topology
    integer(kind=c_int64_t), value :: natoms
end function

! Function "chfl_topology_add_atom", at topology.h:75
function c_chfl_topology_add_atom(topology, atom) bind(C, name="chfl_topology_add_atom")
    use iso_c_binding
    import chfl_status

    implicit none
    integer(kind=chfl_status) :: c_chfl_topology_add_atom
    type(c_ptr), value :: topology
    type(c_ptr), value, intent(in) :: atom
end function

! Function "chfl_topology_remove", at topology.h:86
function c_chfl_topology_remove(topology, i) bind(C, name="chfl_topology_remove")
    use iso_c_binding
    import chfl_status

    implicit none
    integer(kind=chfl_status) :: c_chfl_topology_remove
    type(c_ptr), value :: topology
    integer(kind=c_int64_t), value :: i
end function

! Function "chfl_topology_bonds_count", at topology.h:95
function c_chfl_topology_bonds_count(topology, count) bind(C, name="chfl_topology_bonds_count")
    use iso_c_binding
    import chfl_status

    implicit none
    integer(kind=chfl_status) :: c_chfl_topology_bonds_count
    type(c_ptr), value, intent(in) :: topology
    integer(kind=c_int64_t), intent(inout) :: count
end function

! Function "chfl_topology_angles_count", at topology.h:104
function c_chfl_topology_angles_count(topology, count) bind(C, name="chfl_topology_angles_count")
    use iso_c_binding
    import chfl_status

    implicit none
    integer(kind=chfl_status) :: c_chfl_topology_angles_count
    type(c_ptr), value, intent(in) :: topology
    integer(kind=c_int64_t), intent(inout) :: count
end function

! Function "chfl_topology_dihedrals_count", at topology.h:113
function c_chfl_topology_dihedrals_count(topology, count) bind(C, name="chfl_topology_dihedrals_count")
    use iso_c_binding
    import chfl_status

    implicit none
    integer(kind=chfl_status) :: c_chfl_topology_dihedrals_count
    type(c_ptr), value, intent(in) :: topology
    integer(kind=c_int64_t), intent(inout) :: count
end function

! Function "chfl_topology_impropers_count", at topology.h:122
function c_chfl_topology_impropers_count(topology, count) bind(C, name="chfl_topology_impropers_count")
    use iso_c_binding
    import chfl_status

    implicit none
    integer(kind=chfl_status) :: c_chfl_topology_impropers_count
    type(c_ptr), value, intent(in) :: topology
    integer(kind=c_int64_t), intent(inout) :: count
end function

! Function "chfl_topology_bonds", at topology.h:135
function c_chfl_topology_bonds(topology, data, count) bind(C, name="chfl_topology_bonds")
    use iso_c_binding
    import chfl_status

    implicit none
    integer(kind=chfl_status) :: c_chfl_topology_bonds
    type(c_ptr), value, intent(in) :: topology
    type(c_ptr), value :: data
    integer(kind=c_int64_t), value :: count
end function

! Function "chfl_topology_angles", at topology.h:148
function c_chfl_topology_angles(topology, data, count) bind(C, name="chfl_topology_angles")
    use iso_c_binding
    import chfl_status

    implicit none
    integer(kind=chfl_status) :: c_chfl_topology_angles
    type(c_ptr), value, intent(in) :: topology
    type(c_ptr), value :: data
    integer(kind=c_int64_t), value :: count
end function

! Function "chfl_topology_dihedrals", at topology.h:162
function c_chfl_topology_dihedrals(topology, data, count) bind(C, name="chfl_topology_dihedrals")
    use iso_c_binding
    import chfl_status

    implicit none
    integer(kind=chfl_status) :: c_chfl_topology_dihedrals
    type(c_ptr), value, intent(in) :: topology
    type(c_ptr), value :: data
    integer(kind=c_int64_t), value :: count
end function

! Function "chfl_topology_impropers", at topology.h:176
function c_chfl_topology_impropers(topology, data, count) bind(C, name="chfl_topology_impropers")
    use iso_c_binding
    import chfl_status

    implicit none
    integer(kind=chfl_status) :: c_chfl_topology_impropers
    type(c_ptr), value, intent(in) :: topology
    type(c_ptr), value :: data
    integer(kind=c_int64_t), value :: count
end function

! Function "chfl_topology_add_bond", at topology.h:185
function c_chfl_topology_add_bond(topology, i, j) bind(C, name="chfl_topology_add_bond")
    use iso_c_binding
    import chfl_status

    implicit none
    integer(kind=chfl_status) :: c_chfl_topology_add_bond
    type(c_ptr), value :: topology
    integer(kind=c_int64_t), value :: i
    integer(kind=c_int64_t), value :: j
end function

! Function "chfl_topology_remove_bond", at topology.h:197
function c_chfl_topology_remove_bond(topology, i, j) bind(C, name="chfl_topology_remove_bond")
    use iso_c_binding
    import chfl_status

    implicit none
    integer(kind=chfl_status) :: c_chfl_topology_remove_bond
    type(c_ptr), value :: topology
    integer(kind=c_int64_t), value :: i
    integer(kind=c_int64_t), value :: j
end function

! Function "chfl_topology_clear_bonds", at topology.h:207
function c_chfl_topology_clear_bonds(topology) bind(C, name="chfl_topology_clear_bonds")
    use iso_c_binding
    import chfl_status

    implicit none
    integer(kind=chfl_status) :: c_chfl_topology_clear_bonds
    type(c_ptr), value :: topology
end function

! Function "chfl_topology_residues_count", at topology.h:215
function c_chfl_topology_residues_count(topology, count) bind(C, name="chfl_topology_residues_count")
    use iso_c_binding
    import chfl_status

    implicit none
    integer(kind=chfl_status) :: c_chfl_topology_residues_count
    type(c_ptr), value, intent(in) :: topology
    integer(kind=c_int64_t), intent(inout) :: count
end function

! Function "chfl_topology_add_residue", at topology.h:227
function c_chfl_topology_add_residue(topology, residue) bind(C, name="chfl_topology_add_residue")
    use iso_c_binding
    import chfl_status

    implicit none
    integer(kind=chfl_status) :: c_chfl_topology_add_residue
    type(c_ptr), value :: topology
    type(c_ptr), value, intent(in) :: residue
end function

! Function "chfl_topology_residues_linked", at topology.h:238
function c_chfl_topology_residues_linked(topology, first, second, result) bind(C, name="chfl_topology_residues_linked")
    use iso_c_binding
    import chfl_status

    implicit none
    integer(kind=chfl_status) :: c_chfl_topology_residues_linked
    type(c_ptr), value, intent(in) :: topology
    type(c_ptr), value, intent(in) :: first
    type(c_ptr), value, intent(in) :: second
    logical(kind=c_bool), intent(inout) :: result
end function

! Function "chfl_topology_bond_with_order", at topology.h:251
function c_chfl_topology_bond_with_order(topology, i, j, bond_order) bind(C, name="chfl_topology_bond_with_order")
    use iso_c_binding
    import chfl_status
    import chfl_bond_order

    implicit none
    integer(kind=chfl_status) :: c_chfl_topology_bond_with_order
    type(c_ptr), value :: topology
    integer(kind=c_int64_t), value :: i
    integer(kind=c_int64_t), value :: j
    integer(chfl_bond_order), value :: bond_order
end function

! Function "chfl_topology_bond_orders", at topology.h:265
function c_chfl_topology_bond_orders(topology, orders, nbonds) bind(C, name="chfl_topology_bond_orders")
    use iso_c_binding
    import chfl_status
    import chfl_bond_order

    implicit none
    integer(kind=chfl_status) :: c_chfl_topology_bond_orders
    type(c_ptr), value, intent(in) :: topology
    type(c_ptr), value :: orders
    integer(kind=c_int64_t), value :: nbonds
end function

! Function "chfl_topology_bond_order", at topology.h:278
function c_chfl_topology_bond_order(topology, i, j, order) bind(C, name="chfl_topology_bond_order")
    use iso_c_binding
    import chfl_status
    import chfl_bond_order

    implicit none
    integer(kind=chfl_status) :: c_chfl_topology_bond_order
    type(c_ptr), value, intent(in) :: topology
    integer(kind=c_int64_t), value :: i
    integer(kind=c_int64_t), value :: j
    integer(chfl_bond_order), intent(inout) :: order
end function

end interface