program tariff

! Tariff Program with HPF Directives (XL HPF)
!
! This program solves the linear asymmetric spatial 
! price equilibrium problem with tariffs using the 
! modified projection method.  About the Program
!
! Read data in by subroutine 'getdata'

! The following parameter controls number of iterations in short runs:

        parameter (niter = 5)

	Include 'params.inc'

	real*8 rtc,ti1,ti2,ti3,t1,t2,t3,to1,to2,to3
        

!HPF$ DISTRIBUTE APATHO (*,CYCLIC)

!HPF$ ALIGN (:) WITH apatho(*,:)  ::  demand,demp

! The following array holds temporary apath values:
	real, dimension(maxexten,maxexten):: apathn

!HPF$ ALIGN (i,j) WITH apatho(i,j):: apathn

      INTERFACE

      EXTRINSIC(HPF_LOCAL) &
          SUBROUTINE sayit (t1,t2,message)
             real*8, intent(in)   :: t1,t2
             character*(*), intent(in) :: message
          END SUBROUTINE sayit

      EXTRINSIC(HPF_LOCAL) &
	  SUBROUTINE getdata(sslope,shxlin,dslope,dhxlin,thxlin,ind1,ind2,tar)
             real,intent(out)::  sslope(:, :), shxlin(:)
             real,intent(out)::  dslope(:, :), dhxlin(:)
             real,intent(out)::  thxlin(:, :),tar(:, :)
             integer,intent(out):: ind1(:, :), ind2(:, :)
	  END SUBROUTINE getdata

       END INTERFACE

! Read input data:
! Read in number of supply markets, number of demand markets,
! number of commodities, and number of cross-terms in the price functions

        ti1 = rtc()
        call getdata(sslope,shxlin,dslope,dhxlin,thxlin,ind1,ind2,tar)
        ti2 = rtc()

! Computation

        t1 = rtc()
        npath = nsupm * ndemm
        itcnt = 1
        rho = .0001

! Construct an initial feasible flow assignment
! Most likely can be parallelized (parallelize on j, block on j)

! Initializations

        do 65 j = 1, ndemm
          do 70 l = 1, nsupm
            apatho(l, j) = 0.
70        continue
65      continue
        do 651 i=1, ndemm
           demand(i) = 0.0
 651    continue
        do 652 i=1, nsupm
            supply(i) = 0.0
 652    continue

! (Outer loop starts here)

! Construct new supply price functions
! Compute supply prices, demand prices, and transportation costs

!  Exchange commenting on the following 2 statements for limited run:
!       do while (iict .lt. npath)     ! production run
        do while (itcnt .lt. NITER)    ! limited run

        do 1175 j = 1, nsupm
          supp(j) = shxlin(j)
          do 1176 k = 1, ncr
            supp(j) = supp(j) + sslope(j, k) * supply(ind1(j, k))
1176      continue
1175    continue
        do 1185 j = 1, ndemm
          demp(j) = dhxlin(j)
          do 1187 k = 1, ncr
            demp(j) = demp(j) + dslope(j, k) * demand(ind2(j, k))
1187      continue
1185    continue
        do j = 1, nsupm
          supply(j) = 0.
        enddo

! This is one of the two main iterative steps and can be parallelized

	demand(1:ndemm)=0.

!HPF$ INDEPENDENT,NEW(k,test,apath)

        do 2225 j = 1, ndemm
          do 2245, k = 1, nsupm
            test = (-supp(k) - thxlin(k,j)) * (1. + tar(k,j)) + demp(j)

            apath = max(0., apatho(k,j) + (rho * test))
	    apathn(k, j) = apath

            demand(j) = demand(j) + apath
2245      continue
2225    continue
        supply(1:nsupm) = sum( apathn(1:ndemm, 1:nsupm), 2)

        do 11175 j = 1, nsupm
          supp(j) = shxlin(j)
          do 11176 k = 1, ncr
            supp(j) = supp(j) + sslope(j, k) * supply(ind1(j, k))
11176     continue
11175   continue
        do 11185 j = 1, ndemm
          demp(j) = dhxlin(j)
          do 11187 k = 1, ncr
            demp(j) = demp(j) + dslope(j, k) * demand(ind2(j, k))
11187     continue
11185   continue
        do j = 1, nsupm
          supply(j) = 0.
        enddo

! This iterative step can also be parallelized

        iict = 0
	demand(:ndemm) = 0.

!HPF$ INDEPENDENT,NEW(k,test,apath)
        do 12225 j = 1, ndemm
          do 12230 k = 1, nsupm
            test =(-supp(k) - thxlin(k, j))*(1. + tar(k, j)) + demp(j)

            apath = max( 0., apatho(k, j) + rho * test)
	    apathn(k, j) = apath

            demand(j) = demand(j) + apath
12230     continue
12225   continue

        iict = count(abs(apathn(:nsupm,:ndemm)-apatho(:nsupm,:ndemm)).le. .01)
	supply(1:nsupm) = sum ( apathn(:nsupm,:ndemm), 2)
	apatho(:nsupm,:ndemm)=apathn(:nsupm,:ndemm)

        itcnt = itcnt + 1

        end do
! (End of outer loop)

        t2 = rtc()

! Output section

        to1 = rtc()
        write (*, 10) nsupm, ndemm, ncomm, ncr
10      format(' nsupm=',I5,'   ndemm=',I5,'   ncomm=',I5,'   ncr=',I5/)

        ti3 = ti2 - ti1
        t3 = t2 - t1

        write (*, 890) ti3, t3, itcnt, iict
890     format ('    DATA INPUT WALLCLOCK TIME IN SECONDS=',f12.4,/        &
     & '   COMPUTATION WALLCLOCK TIME IN SECONDS=',f12.4,//                &
     & ' itcnt=',i7,' iict=',i9)
        
	itc  =  count (apatho(:nsupm, :ndemm) .gt. 0.)
        
        perc = (float(itc) / float(npath)) * 100.
        write (*, 138) perc, itc
138     format (' perc=', f10.4, 2x, 'itc = ', i5/)

        to2 = rtc()

        to3 = to2 - to1

        write (*, 895) to3
 895    format (' REPORT OUTPUT WALLCLOCK TIME IN SECONDS=',f12.4)

        write (*, 899) (to2 - ti1)
 899    format (' OVERALL WALLCLOCK TIME IN SECONDS=',f12.4)

        stop

        end