      program Model
c----------------------------------------------------------------------------------
c LHC Experiment Computing Model Simulation
c
c JJBunn 1996 
c----------------------------------------------------------------------------------
      parameter (lhbook=500000)
      common /PAWC/ h(lhbook)
      parameter (maxinst=200,maxvc=100000,maxtask=100000)
      parameter (maxsteps=6,maxtypes=10,maxmix=10)
      parameter (lget_data_local=1,lput_data_local=2,lproc_data=3)
      parameter (lget_data_CERN=11,lput_data_CERN=12)
      parameter (lget_data_random=21,lput_data_random=32)
c maxinst    : Maximum number of collaborating institutes
c maxvc      : Maximum number of network connections between all institutes
c maxtask    : Maximum number of tasks that can be active at one time
c maxsteps   : Maximum number of steps in a task
c maxtypes   : Maximum number of different types of task
c maxmix     : maximum number of tasks that can be mixed at an institute
      character*200 cline
      character*80 title
      character*30 institute(maxinst)
      character*30 taskname(maxtypes)
      character*12 stepname(50)
      character*5 ifmt
      character*7 ffmt
      integer source(maxvc),sink(maxvc),depend(maxvc)
c source     : source of the data transfer (institute number)
c sink       : sink of the transfer (institute number)
c depend     : the number of the task that initiated the transfer
      integer taskstep(maxtask),tasktype(maxtask),location(maxtask)
      integer tasksmix(maxtask)
c taskstep   : the current step of the task
c tasktype   : the type of the task
c location   : where the task is running
c tasksmix   : the mix number of the task
      integer steps(maxtypes)
      integer steptype(maxtypes,maxsteps)
c steps      : the number of steps in the task type
c steptype   : for each task type, the type of step (PROC,GET,PUT)
      integer nmix(maxinst),mix(maxinst,maxmix)
      integer number(maxinst,maxmix)
      integer TASKSTARTED(maxinst,maxmix),TASKFINISHED(maxinst,maxmix)
      integer ntasks(maxinst),TASKMAX(maxinst)
      integer nincpu(maxinst),connect(maxinst,maxinst),CPUMAX(maxinst)
      integer WANMAX(maxinst)
c nmix       : the number of different tasks in the mix 
c mix        : the task type for each task in the mix
c number     : the number required to run over the simulation period
c TASKSTARTED: the number of tasks started in this mix
c TASKFINISHED the number of tasks that have finished in this mix
c ntasks     : the total number of tasks running at the institute
c TASKMAX    : the maximum number of concurrent tasks seen at the institute
c nincpu     : the number of tasks in a CPU step at the institute
c connect    : the matrix of current connections between institutes
c CPUMAX     : the maximum number of tasks seen in CPU at the institute
c WANMAX     : the maximum number of WAN connections from the institute
      integer seed
      real stepsize(maxtypes,maxsteps)
c stepsize   : the size of the step (seconds if CPU, MBytes if GET,PUT)
      real rate(maxvc)
      real done(maxvc),total(maxvc),startvc(maxvc)
c rate       : the current transfer rate for the Virtual Circuit
c done       : the amount of data already transferred
c total      : the amount of data required to be transferred
c startvc    : the time at which the Virtual Circuit was established
      real starttask(maxtask),startcpu(maxtask)
c starttask  : the time at which the task was started
c startcpu   : the time at which the PROC (CPU) step of the task started
      real frequency(maxinst,maxmix)
      real ratematrix(maxinst,maxinst)
      real diskused(maxinst)
      real costdisk(maxinst)
      real costnet(maxinst)
      real costproc(maxinst)
      real rateLAN(maxinst)
      real TASKTIME(maxinst,maxmix)
c frequency  : the frequency with which tasks of this mix need to be generated
c ratematrix : the ATM connection speed for Virtual Circuits between institutes
c diskused   : the total disk space in use at the Institute
c costdisk   : Cost of disk space in dollars per MByte
c costnet    : Cost of ATM line in dollars per second at specified ATM rate
c costproc   : Cost of 100 MIPS processor in dollars
c rateLAN    : the current LAN aggregate data transfer rate
c TASKTIME   : the accumulated time for all tasks of this mix
      logical debug
      logical filhis
      real*8 time,time_to_snap
c
      data ifmt /'(i 5)'/
      data ffmt /'(g 5.0)'/
      data seed /987654/
c     
c Initialise HBOOK
c
      Call HLIMIT(lhbook)
c
      lun = 6
      ninst = 0
      ntypes = 0
      snapshot = 10.
      endtime = 100.
      atm = 155.
      debug = .false.
      stepname(lget_data_local) = 'GET data local'
      stepname(lput_data_local) = 'PUT data local'
      stepname(lget_data_CERN) = 'GET data from CERN'
      stepname(lput_data_CERN) = 'PUT data to CERN'
      stepname(lget_data_random) = 'GET data from ?'
      stepname(lput_data_random) = 'PUT data to ?'
      stepname(lproc_data)= 'PROCESS data'
c
c We open the data file and read in the model parameters,
c checking for inconsistencies etc.
c
c      open(1,file='e:\model\model.dat',status='old',err=999)
      open(1,file='model.dat',status='old',err=999)
 1000 read(1,'(a)',end=1001,err=900) cline
      lline = lenocc(cline)
      if(lline.eq.0.or.cline(1:1).eq.'#') goto 1000
      if(cline(:6).eq.'Title:') then
         title = cline(7:)
         ltitle = lenocc(title)
         goto 1000
      else if(cline(:5).eq.'debug') then
         debug = .true.
         goto 1000
      else if(cline(:7).eq.'Output:') then
         lun = 10
         open(lun,file=cline(8:lline),status='new')
         Call HOUTPU(lun)
         Call HERMES(lun)
         goto 1000
      else if(cline(:4).eq.'ATM:') then
         write(ffmt(3:4),'(i2)') lline-4
         read(cline(5:lline),ffmt) atm
         goto 1000
      else if(cline(:10).eq.'Institute:') then
         if(ninst.ge.maxinst) stop 'Too many Institutes!'
         ninst = ninst + 1
         institute(ninst) = cline(11:)
         nmix(ninst) = 0
c
c default costs: 
c 50 dollars per GByte
c 100000 dollars per year per external line
c 500 dollars per 1000 MIPS
c
         costdisk(ninst) = 50./1000.               ! cost per MByte
         costnet(ninst) = 100000./365./24./60./60. ! cost per second
         costproc(ninst) = 50.                     ! cost per 100 MIPS
         goto 1000
      else if(cline(:9).eq.'CostDisk:') then
         write(ffmt(3:4),'(i2)') lline-9
         read(cline(10:lline),ffmt) costdisk(ninst)
         costdisk(ninst) = costdisk(ninst)/1000.   ! converts to cost per MByte
         goto 1000
      else if(cline(:8).eq.'CostNet:') then
         write(ffmt(3:4),'(i2)') lline-8
         read(cline(9:lline),ffmt) costnet(ninst)
         costnet(ninst) = costnet(ninst)/365./24./60./60.
         goto 1000
      else if(cline(:8).eq.'CostCPU:') then
         write(ffmt(3:4),'(i2)') lline-8
         read(cline(9:lline),ffmt) costproc(ninst)
         goto 1000
      else if(cline(:13).eq.'Tasks in Mix:') then
         write(ifmt(3:4),'(i2)') lline-13
         read(cline(14:lline),ifmt) nmix(ninst)
         if(nmix(ninst).gt.maxmix) stop 'Too many jobs in mix'
         do imix=1,nmix(ninst)
	      read(1,'(a)',end=900,err=900) cline
            lline = lenocc(cline)
            if(cline(:5).ne.'Type:') stop 'Bad sequence'
            write(ifmt(3:4),'(i2)') lline-5
            read(cline(6:lline),ifmt) mix(ninst,imix)
	      read(1,'(a)',end=900,err=900) cline
            lline = lenocc(cline)
            if(cline(:10).ne.'Frequency:') stop 'Bad sequence'
            write(ffmt(3:4),'(i2)') lline-10
            read(cline(11:lline),ffmt) frequency(ninst,imix)
	    read(1,'(a)',end=900,err=900) cline
            lline = lenocc(cline)
            if(cline(:7).ne.'Number:') stop 'Bad sequence'
            write(ifmt(3:4),'(i2)') lline-7
            read(cline(8:lline),ifmt) number(ninst,imix)
            if(number(ninst,imix).ne.0.and.
     &         frequency(ninst,imix).ne.0.0) then
               Stop 'Bad Frequency/Number combination'
            endif
         end do
         goto 1000
      else if(cline(:5).eq.'Task:') then
         if(ntypes.ge.maxtypes) stop 'Too many task types'
         ntypes = ntypes + 1
         taskname(ntypes) = cline(6:)
         goto 1000   
      else if(cline(:6).eq.'Steps:') then
         write(ifmt(3:4),'(i2)') lline-6
         read(cline(7:lline),ifmt) steps(ntypes)
         if(steps(ntypes).gt.maxsteps) stop 'Too many steps'
         do istep=1,steps(ntypes)
	      read(1,'(a)',end=900,err=900) cline
            lline = lenocc(cline)
            if(cline(:9).ne.'Steptype:') stop 'Bad sequence'
            if(cline(10:).eq.'GET LOCAL') then
               steptype(ntypes,istep) = lget_data_local
            else if(cline(10:).eq.'PUT LOCAL') then
               steptype(ntypes,istep) = lput_data_local
            else if(cline(10:).eq.'GET CERN') then
               steptype(ntypes,istep) = lget_data_CERN
            else if(cline(10:).eq.'PUT CERN') then
               steptype(ntypes,istep) = lput_data_CERN
            else if(cline(10:).eq.'GET RANDOM') then
               steptype(ntypes,istep) = lget_data_random
            else if(cline(10:).eq.'PUT RANDOM') then
               steptype(ntypes,istep) = lput_data_random
            else if(cline(10:).eq.'PROC') then
               steptype(ntypes,istep) = lproc_data
            else
               stop 'Unknown step type'
            endif
	      read(1,'(a)',end=900,err=900) cline
            lline = lenocc(cline)
            if(cline(:9).ne.'Stepsize:') stop 'Bad sequence'
            write(ffmt(3:4),'(i2)') lline-9
            read(cline(10:lline),ffmt) stepsize(ntypes,istep)
         end do
         goto 1000
	else if(cline(:10).eq.'TimeLimit:') then
         write(ffmt(3:4),'(i2)') lline-10
         read(cline(11:lline),ffmt) endtime
         goto 1000
      else if(cline(:9).eq.'SnapShot:') then
         write(ffmt(3:4),'(i2)') lline-9
         read(cline(10:lline),ffmt) snapshot
         goto 1000
      else
         write(lun,*) 'Unknown keyword:'//cline(:lline)
         goto 1000
      endif
 1001 close(1)
c
c We summarise the model parameters read in
c
      write(lun,*) 'Input Summary for "',title(:ltitle),'"'
      write(lun,*) ' Total of ',ntypes,' task types'
      do itype=1,ntypes
         write(lun,60) taskname(itype),steps(itype),
     &               (stepname(steptype(itype,i)),
     &               stepsize(itype,i),i=1,steps(itype))
   60    format(3x,a,i4,' step(s):',10(/,5x,a,' size:',f9.3))      
      end do
      write(lun,*) ' ATM Rate:',atm,' Mbits/second'
      write(lun,*) ' Total of ',ninst,' institutes'
      do inst=1,ninst
         diskused(inst) = 0.0
         WANMAX(inst) = 0
         CPUMAX(inst) = 0
         if(nmix(inst).eq.0) then
            write(lun,50) institute(inst)
         else
            write(lun,50) institute(inst),nmix(inst),
     &        (taskname(mix(inst,i))(:lenocc(taskname(mix(inst,i)))),
     &         frequency(inst,i),number(inst,i),i=1,nmix(inst))
         endif
   50 format(3x,a,i4,' task(s):',10(/,5x,a,' freq:',f9.3,' num:',i10))
         write(lun,*) '    Disk:',costdisk(inst),'$ per GByte'
         write(lun,*) '    WAN: ',costnet(inst),'$ per line per sec'
         write(lun,*) '    CPU: ',costproc(inst),'$ per 100 MIPS'
      end do
c
c Set up rate matrix
c
      do inst1=1,ninst
         do inst2=1,ninst
            ratematrix(inst1,inst2) = atm/8.0  ! ATM
c            if(inst1.eq.inst2) ratematrix(inst1,inst2) = 100.0
         end do
      end do
      write(lun,*) 'Network Point to Point Rate Matrix ...'
      do inst=1,ninst
         write(lun,'(1x,a,10(f6.1,1x))') institute(inst),
     &                           (ratematrix(inst,j),j=1,ninst)
      end do
c
c determine minimum data set size, estimate frequency for tasks that must
c run a given number of times, and zero some arrays
c
      datamin = 9999.
      dtime_tasks = 9999.
      maxrough = 25
      do inst=1,ninst
         do imix=1,nmix(inst)
            TASKTIME(inst,imix) = 0.0
            TASKSTARTED(inst,imix) = 0
            TASKFINISHED(inst,imix) = 0
            itype = mix(inst,imix)
            nrequired = number(inst,imix)
            timetot = 0.0
            datatot = 0.0   
            do istep=1,steps(itype)
               if(steptype(itype,istep).ne.lproc_data) then
                  datamin = min(datamin,stepsize(itype,istep))
                  datatot = datatot + stepsize(itype,istep)
               else 
                  timetot = timetot + stepsize(itype,istep)
               endif
            end do
            if(nrequired.gt.0) then
               timetot = max(timetot,1.0)
               timetot_s = timetot
               nrough = 0
   70          if(timetot.ge.endtime.and.nrough.ne.0) then
                  stop 'Not enough time !'
               endif
               time_needed = real(nrequired)*timetot
               if(time_needed.gt.endtime) then
c
c nrough is approximate number of concurrent tasks of this type
c
                  nrough_new = nint(time_needed/(endtime-timetot))                  
                  if(datatot.gt.0..and.nrough_new.gt.nrough) then
                     nrough = nrough_new
                     timetot = timetot_s + datatot/(ATM/nrough/8.)
                     goto 70
                  endif
                  nrough = nrough_new
                  maxrough = max(nrough,maxrough)
               endif
               frequency(inst,imix)=(endtime-timetot)/nrequired
               dtime_tasks=min(dtime_tasks,frequency(inst,imix)) 
               linst = lenocc(institute(inst))
               ltask = lenocc(taskname(mix(inst,imix)))
               write(lun,700) institute(inst)(:linst),nrough,
     &                        taskname(mix(inst,imix))(:ltask)
  700 format(1x,a15,': need approx ',i5,' concurrent "',a,'" tasks')
            endif
         end do
      end do
c
      write(lun,*) 'Task slice minimum:',dtime_tasks
      write(lun,*) 'Smallest data set size is:',datamin,' MBytes'
c
c determine maximum network bandwidth between institutes
c
      ratemax = 0.0
      do inst1=1,ninst
         do inst2=1,ninst
            ratemax = max(ratemax,ratematrix(inst1,inst2))
         end do
      end do
c
      write(lun,*) 'Maximum network bandwidth is:',ratemax,' MBytes/sec'
      write(lun,'(a,f10.1,a)') ' Simulation over ',endtime,' seconds'
c
c Book some histograms
c
      do i=1,ninst
         if(nmix(i).gt.0) then
           Call HBOOK1(1000+i,'Evolution of tasks at '//institute(i),
     &     50,0.,endtime,0.)
           Call HBOOK1(2000+i,'VC Rate at '//institute(i),
     &                 25,0.,ATM/8.,0.)
           Call HBOOK1(3000+i,'Task Count at '//institute(i),
     &                 25,0.,20.*real(maxrough),0.)
           Call HBOOK1(4000+i,'CPU Count at '//institute(i),
     &                 25,0.,20.*real(maxrough),0.)
         endif
      end do
      Call HMINIM(0,0.0)
c
      nvc = 0
      ntask = 0
      time = 0.0
      tbin = 0.5*endtime/50.
      dtime_nominal = 0.1*datamin/ratemax
      dtime = min(dtime_nominal,dtime_tasks)
      write(lun,*) 'Time slice will be:',dtime,' seconds'
      time_to_snap = snapshot
c
c Main time slice loop ...
c
    1 time = time + dtime
c
c zero some counters
c
      do inst=1,ninst
         rateLAN(inst) = 0.0
         ntasks(inst) = 0
         nincpu(inst) = 0
         do inst2 = 1,ninst
            connect(inst,inst2) = 0
         end do
      end do
c
      if(time.gt.tbin) then
         filhis = .true.
         tbin = tbin + (endtime/50.)
      else
         filhis = .false.
      endif
c
c Calculate connection matrix
c
      do ivc=1,nvc
         isource = source(ivc)
         isink   = sink(ivc)
         connect(isource,isink) = connect(isource,isink) + 1
         connect(isink,isource) = connect(isink,isource) + 1
         if(isource.eq.isink) rateLAN(isink)=rateLAN(isink)+rate(ivc)
      end do
c
c Update distinct WAN count for each institute
c
      do inst=1,ninst
         nlines = 0
         do inst2=1,ninst
            if(inst.ne.inst2.and.connect(inst,inst2).ne.0) then
               nlines = nlines + 1
            endif
         end do
         WANMAX(inst) = max(WANMAX(inst),nlines)
      end do
c
c Output status every snapshot seconds
c
      time_to_snap = time_to_snap - dtime
      if(time_to_snap.lt.0.0) then
         time_to_snap = snapshot
         write(lun,400) sngl(time),dtime,nvc,ntask
  400    format(' Time:',f10.1,' dT:',f10.7,' VCs:',i6,' Tasks:',i6)
         do inst=1,ninst
            linst = lenocc(institute(inst))
            if(rateLAN(inst).ne.0.0) then
             write(lun,500) institute(inst)(:linst),
     &            ': LAN Rate',rateLAN(inst)
	    endif
	end do
  500 format(1x,a,a,f10.5)
      endif
c
c We look at all the Virtual Circuits and check whether the required
c data has been transferred
c
      ivc = 0
    2 ivc = ivc + 1
      if(ivc.le.nvc) then
c The data transferred in this step = rate*time slice
         done(ivc) = done(ivc) + rate(ivc)*dtime
c
c adjust WAN rate according to current number of connections ...
c
         isink = sink(ivc)
         isource = source(ivc)
         if(isink.ne.isource) then
            rate(ivc) = ratematrix(isource,isink)/
     &                  connect(isource,isink)
            Call HF1(2000+isink,rate(ivc),1.)
            Call HF1(2000+isource,rate(ivc),1.)
         endif
         if(done(ivc).ge.total(ivc)) then
c
c Add data transferred to local disk
c
            diskused(isource) = diskused(isource)+total(ivc)
            if(isink.ne.isource) then
               diskused(isink) = diskused(isink)+total(ivc)
            endif
c
c Find task number which depends on this VC
c
            itask = depend(ivc)
c
c move task to next step, or finish it
c
            itype = tasktype(itask)
            istep = taskstep(itask)
            if(istep+1.gt.steps(itype)) then
               taskstep(itask) = 0
            else
c
c A negative step will be used below to trigger the move into that
c step. We could do it here, but it makes the code more compact
c doing it in one place.
c
               taskstep(itask) = -(istep+1)
            endif
            do i=ivc,nvc-1
               source(i) = source(i+1)
               sink(i)   = sink(i+1)
               depend(i) = depend(i+1)
               rate(i)   = rate(i+1)
               done(i)   = done(i+1)
               total(i)  = total(i+1)
               startvc(i)  = startvc(i+1)
            end do
            nvc = nvc - 1
            ivc = ivc - 1
         endif
         goto 2
      endif
c
      itask = 0
    3 itask = itask + 1
      if(itask.le.ntask) then
         itype = tasktype(itask)
         istep = taskstep(itask)
         inst  = location(itask)
c 
c Jump away if the task is marked as finished ...
c
         if(istep.eq.0) goto 4
c
c Update task count at this institute
c
         ntasks(inst) = ntasks(inst) + 1
         TASKMAX(inst) = max(TASKMAX(inst),ntasks(inst))
         if(istep.lt.0) then
c Task needs moving to abs(istep)
            istep = abs(istep)
            taskstep(itask) = istep 
            isteptype = steptype(itype,istep)
            if(isteptype.eq.lput_data_local.or.
     &         isteptype.eq.lput_data_CERN.or.
     &         isteptype.eq.lput_data_random) then
c Step will be PUT of data
               if(nvc.ge.maxvc) stop 'Too many Virtual Circuits!'
               nvc         = nvc + 1
               source(nvc) = location(itask)
               sink(nvc)   = 1
               if(isteptype.eq.lput_data_local) then
                  sink(nvc)= location(itask)
               else if(isteptype.eq.lput_data_random) then
                  sink(nvc)= nint(real(ninst-1)*rndm(dummy))+1
               endif
               depend(nvc) = itask
               nconn       = connect(source(nvc),sink(nvc)) + 1
               rate(nvc)   = ratematrix(sink(nvc),source(nvc))/nconn
               if(sink(nvc).eq.source(nvc)) rate(nvc) = atm/8.
               done(nvc)   = 0.0
               total(nvc)  = stepsize(itype,istep)
               startvc(nvc)= time
            else if(isteptype.eq.lget_data_local.or.
     &              isteptype.eq.lget_data_CERN.or.
     &              isteptype.eq.lget_data_random) then
c Step will be GET of data
               if(nvc.ge.maxvc) stop 'Too many Virtual Circuits!'
               nvc         = nvc + 1
               source(nvc)   = 1
               if(isteptype.eq.lget_data_local) then
                  source(nvc)= location(itask)
               else if(isteptype.eq.lget_data_random) then
                  source(nvc)= nint(real(ninst-1)*rndm(dummy))+1
               endif
               sink(nvc)   = location(itask)
               depend(nvc) = itask
               nconn       = connect(source(nvc),sink(nvc)) + 1
               rate(nvc)   = ratematrix(sink(nvc),source(nvc))/nconn
               if(sink(nvc).eq.source(nvc)) rate(nvc) = atm/8.
               done(nvc)   = 0.0
               total(nvc)  = stepsize(itype,istep)
               startvc(nvc)= time
            else if(steptype(itype,istep).eq.lproc_data) then
c Step will be CPU (processing)
               startcpu(itask) = time
            else
               stop 'Unknown task step !'
            endif
         else if(steptype(itype,istep).eq.lproc_data) then
c Task in CPU mode
            nincpu(inst) = nincpu(inst)+1
            CPUMAX(inst) = max(CPUMAX(inst),nincpu(inst))
            if(time-startcpu(itask).ge.stepsize(itype,istep)) then
c  and has finished ... move to next step or task finish
               if(istep.ge.steps(itype)) goto 4
               taskstep(itask) = -(istep+1)
	    endif            
         endif
c ... and get next task
         goto 3
c
c The task finished, accumulate some statistics and
c remove task from list
c
    4    imix = tasksmix(itask)
         TASKTIME(inst,imix) = TASKTIME(inst,imix) +
     &                         time - starttask(itask)
         TASKFINISHED(inst,imix) = TASKFINISHED(inst,imix) + 1
         if(debug) write(lun,200) time,institute(inst),
     &             ' End Task:',taskname(tasktype(itask))
c
c We shuffle all the tasks down in the stack ...
c
         do i=itask,ntask-1
            tasktype(i) = tasktype(i+1)
            taskstep(i) = taskstep(i+1)
            starttask(i)= starttask(i+1)
            startcpu(i) = startcpu(i+1)
            location(i) = location(i+1)
            tasksmix(i) = tasksmix(i+1)
         end do
         do ivc=1,nvc
            if(depend(ivc).ge.itask) depend(ivc) = depend(ivc)-1
         end do
         ntask = ntask - 1
         itask = itask - 1
         goto 3
      endif 
c
c Determine if a new task starts
c
      dtime_tasks = 9999.
      do inst=1,ninst
         if(nmix(inst).ne.0) then
c
c Fill histograms
c
            if(filhis) Call HF1(1000+inst,sngl(time),real(ntasks(inst)))
            Call HF1(3000+inst,real(ntasks(inst)),1.)
            Call HF1(4000+inst,real(nincpu(inst)),1.)
            do imix=1,nmix(inst)
c
c There are two sorts of task: those which start randomly every 
c "frequency" seconds (on average), and those where the "number"
c is to be reached for the given simulation time. In the latter
c case, we allocate an appropriate frequency ...
c 
               itype = mix(inst,imix)
               ndone = TASKFINISHED(inst,imix)
               val = rndm(seed)
c               call random(val)
               if(number(inst,imix).ne.0.and.ndone.ne.0) then
c We calculate the number of tasks left to do, the current average time needed
c for each task, and then update the frequency at which the tasks need to be
c generated.
                  ntodo = number(inst,imix) - TASKSTARTED(inst,imix)
                  if(ntodo.gt.0) then
                     average_time = TASKTIME(inst,imix)/real(ndone)
                     freq_needed = (endtime-time-average_time)/
     &                             real(ntodo)
                     frequency(inst,imix) = max(dtime,freq_needed)
                     dtime_tasks=min(dtime_tasks,frequency(inst,imix))
                  else
c If no more tasks are required, we ensure another is not generated by
c setting val high ...
                     frequency(inst,imix) = 9999.
                     val = 9999.
                  endif
               endif
               if(dtime/frequency(inst,imix).gt.val) then
                  if(ntask.ge.maxtask) stop 'Too many tasks !'
c
c A new task is to be started .... fill in the required information
c
                  ntask = ntask + 1
                  tasktype(ntask) = itype
                  taskstep(ntask) = 1
                  startcpu(ntask) = 0.0
                  starttask(ntask)= time
                  location(ntask) = inst
                  tasksmix(ntask) = imix
                  TASKSTARTED(inst,imix) = TASKSTARTED(inst,imix)+1
                  if(debug) write(lun,200) time,institute(inst),
     &            ' New Task:',taskname(tasktype(ntask))
  200 format(' Time:',f10.2,' Inst:',3a)
c
c Start the task on the appropriate step
c
                  isteptype = steptype(itype,1)
                  if(isteptype.eq.lput_data_local.or.
     &               isteptype.eq.lput_data_CERN.or.
     &               isteptype.eq.lput_data_random) then
c Step will be PUT of data
                     if(nvc.ge.maxvc) stop 'Too many Virtual Circuits!'
                     nvc         = nvc + 1
                     source(nvc) = inst
                     sink(nvc)   = 1
                     if(isteptype.eq.lput_data_local) then
                        sink(nvc)= inst
                     else if(isteptype.eq.lput_data_random) then
                        sink(nvc)= nint(real(ninst-1)*rndm(dummy))+1
                     endif
                     depend(nvc) = ntask
                     nconn = connect(source(nvc),sink(nvc)) + 1
                     rate(nvc) = ratematrix(sink(nvc),source(nvc))/nconn
                     if(sink(nvc).eq.source(nvc)) rate(nvc) = atm/8.
                     done(nvc)   = 0.0
                     total(nvc)  = stepsize(itype,1)
                     startvc(nvc)= time
                  else if(isteptype.eq.lget_data_local.or.
     &                    isteptype.eq.lget_data_CERN.or.
     &                    isteptype.eq.lget_data_random) then
c Step will be GET of data
                     if(nvc.ge.maxvc) stop 'Too many Virtual Circuits!'
                     nvc         = nvc + 1
                     source(nvc) = 1
                     sink(nvc)   = inst
                     if(isteptype.eq.lget_data_local) then
                        source(nvc)= inst
                     else if(isteptype.eq.lget_data_random) then
                        source(nvc)= nint(real(ninst-1)*rndm(dummy))+1
                     endif
                     depend(nvc) = ntask
                     nconn = connect(source(nvc),sink(nvc)) + 1
                     rate(nvc) = ratematrix(sink(nvc),source(nvc))/nconn
                     if(sink(nvc).eq.source(nvc)) rate(nvc) = atm/8.
                     done(nvc)   = 0.0
                     total(nvc)  = stepsize(itype,1)
                     startvc(nvc)= time
                  else if(steptype(itype,1).eq.lproc_data) then
                     startcpu(ntask) = time
                  else
                     stop 'Unknown task step !'
                  endif
               endif
            end do
         endif
      end do     
c
c We constantly adjust dtime so that iterations are kept to a minimum
c
      if(dtime_tasks.ne.9999.) dtime = min(dtime_tasks,dtime_nominal)
      if(time.lt.endtime) goto 1
      write(lun,*) 'End time:',time,' with ',ntask,' task(s) executing'
c
c Select HBOOK printing options
c
      Call HIDOPT(0,'BLAC')
      Call HSQUEZ('YES')
      Call HPAGSZ(30)
      total_disk_cost = 0.0
      total_CPU_cost = 0.0
      total_WAN_cost = 0.0
      ntask_total = 0
      do inst=1,ninst
         write(lun,'(/,a)') ' Statistics for :'//institute(inst)
         write(lun,*) ' Disk occupancy :',diskused(inst),' MBytes'
         write(lun,*) ' WAN lines used :',WANMAX(inst)
         write(lun,*) ' Processors used:',CPUMAX(inst),' 100MIPS units'
         cd = diskused(inst)*costdisk(inst)
         cn = real(WANMAX(inst))*time*costnet(inst)
         cc = real(CPUMAX(inst))*costproc(inst)
         total_disk_cost = total_disk_cost + cd
         total_CPU_cost = total_CPU_cost + cc
         total_WAN_cost = total_WAN_cost + cn
         ct = cn + cc + cd
         write(lun,'(a,f15.1,a)')   '   Disk cost  : ',cd,'$'
         write(lun,'(a,f15.1,a)')   '   WAN cost   : ',cn,'$'
         write(lun,'(a,f15.1,a)')   '   CPU cost   : ',cc,'$'
         write(lun,'(/,a,f15.1,a,/)') '   Total cost : ',ct,'$'
         do imix=1,nmix(inst)
            itype = mix(inst,imix)
            average_time = TASKTIME(inst,imix)/TASKFINISHED(inst,imix)
            ntask_total = ntask_total+TASKFINISHED(inst,imix)
            proc_tot = 0.0
            data_tot = 0.0
            do istep=1,steps(itype)
               if(steptype(itype,istep).eq.lproc_data) then
                  proc_tot = proc_tot + stepsize(itype,istep)
               else
                  data_tot = data_tot + stepsize(itype,istep)
               endif
            end do
            time_data = average_time - proc_tot
            data_rate = data_tot/time_data
            write(lun,300) taskname(itype)(:lenocc(taskname(itype))),
     &      TASKFINISHED(inst,imix),average_time,data_rate
  300       format(3x,'Task:',a,' Count:',i7,' Time(ave):',f10.5,
     &             ' I/O Rate(ave):',f7.2)
         end do
         if (nmix(inst).gt.0) then
            Call HPRINT(1000+inst)
            Call HPRINT(2000+inst)
            Call HPRINT(3000+inst)
            Call HPRINT(4000+inst)
         endif
      end do
      total_cost = total_CPU_cost + total_disk_cost + total_WAN_cost
      write(lun,'(//)')
      write(lun,*) 'Total CPU Cost this Model:   ',total_CPU_cost,'$'
      write(lun,*) 'Total Disk Cost this Model:  ',total_disk_cost,'$'
      write(lun,*) 'Total WAN Cost this Model:   ',total_WAN_cost,'$'
      write(lun,*) 'Total tasks this Model:      ',ntask_total
      write(lun,*) 'Grand Total Cost this Model: ',total_cost,'$'
      close(lun)
      stop 'Finished Simulation'
  900 stop 'Premature end of file model.dat'
  999 stop 'Cannot open model.dat file'
      end

      integer function lenocc(string)
      character*(*) string
      lenocc = 0
      do i=len(string),1,-1
         if(string(i:i).ne.' ') then
            lenocc = i
            return
         endif
      end do
      end
