c   thermal performance of water cooled conductor
c
c   input--area of copper, hole dia, total length,incremental length
c   names--acu,dia,x,dx (4f10.3, sq in, in, ft, ft)
c   input--velocity of coolant,inlet temp,ambient copper temp
c   names--v,tw0,t0 (3f10.3,ft/sec,deg c, deg c)
c   input--current, esw time,rep rate, time steps per cycle, number of pulses
c   names--amps,esw,rep,steps,pulses (5f10.3,amps,sec,sec,unit,sec)
      dimension t(500),tw(500),tout(1000),tstart(1000),tmid(1000)
      dimension tend(1000),tcout(1000),tcend(1000)
      dimension risout(1000)
      integer n,nmid,nskip,ncount,ncskip,nccount
100   format(4f10.3)
110   format(3f10.3)
120   format(5f10.3)
130   format(37x,'Thermal Performance of H2O Cooled Conductor')
140   format(37x,'-------------------------------------------')
150   format(5x,'Cross section = ',f10.3,'sq in')
160   format(5x,'Coolant hole dia = ',f10.3,'in')
170   format(5x,'Water velocity = ',f10.3,'fps')
180   format(5x,'Init con water, copper = ',f10.3,',',f10.3,'degC')
190   format(5x,'ESW = ',f10.3,'/',f10.3,'amp-sec every',f10.3,'sec')
200   format(/'      Time = ',f10.3)
210   format(5x,'Resistance =',f10.3,' milliohms')
220   format('      Copper temp '/(1x,20f6.1))
230   format('      Water temp '/(1x,20f6.1))
240   format(5x,'Tau = ',f10.3)
250   format(5x,'Time step dt = ',f10.3)
      open(10,file='cool.d',status='old')
      open(20,file='cool.o',status='old')
      open(30,file='cool.w',status='old')
      open(40,file='cool.r',status='old')
      open(50,file='cool.i2t',status='old')
      read(10,100)acu,dia,x,dx
      read(10,110)v,tw0,t0
      read(10,120)amps,esw,rep,steps,pulses
      read(10,*)jadiabatic
      read(10,*)tauris,risi2t
      write(20,*)'# Pulses = ',pulses
      write(20,130)
      write(20,140)
      write(20,150)acu
      write(20,160)dia
      write(20,170)v
      write(20,180)tw0,t0
      write(20,190)amps,esw,rep
c   copper: resistivity at 20C, res temp coeff,spec heat, density
c   units: ohm-cm, per degC, j/gm-degC, gm/cc
      restvy=1.74e-6
      coeff=0.0041
      sphtcu=0.386
      dens=8.9
c   water: thermal cond, spec heat, viscosity
c   units: (w/sq cm) (degC/cm),j/gm-degC,gm/cm-sec
      tconw=0.00598
      sphtw=4.2
      viscw=0.00994
c   generalize quantities for section dx
      hcapcu=acu*2.54*2.54*dx*12.*2.54*dens*sphtcu
      hcapw=3.141592654*((dia*2.54)**2)/4.*dx*12.*2.54*sphtw
c   hxfrdx=heat transfer conductance from copper to h2o (w/degC)
      a=0.023*((dia*2.54*v*12.*2.54/viscw)**0.8)
      b=((sphtw*viscw/tconw)**0.4)*tconw*3.141592654*dx*12.*2.54
      hxfrdx=a*b
c
c   csa is cross section in m^2
c
      csa=acu*2.54**2./100**2.
c   tau = thermal time constant (approx)
      tau=(hcapcu+hcapw)/hxfrdx
      write(20,240)tau
c   determine time step
      dt=rep/steps
      write(20,250)dt
c   initialize
      res=restvy*dx*12.*2.54/(acu*2.54*2.54)
      r=res*(1.+coeff*(t0-20.))
      pulse=0
      qin=0.
      qxfr=0.
      time=0.
      ctime=0.
      dist=0.
      curr=amps
      ris=0.
      rismax=0.
      tmax=0.
      n=int(x/dx)
      nmid=int(n/2)
      nskip=int(steps/1000)
      ncount=1
      ncskip=int(steps*pulses/1000)
      nccount=1
      jcout=0
      ohms=n*res
      oohms=ohms*1.e3
      do 500 i=1,n
      t(i)=t0
      tw(i)=tw0
500   continue
      jprint=1
c   enter time loop, print at start and end of pulse
505   continue
      if(time.eq.0.)then
c
c   pulse is starting
c
c   reduce current if dt is greater than esw
c
      if(dt.gt.esw)then
		curr=sqrt(amps**2.*esw/dt)
	else
		curr=amps
      end if
      action=0.
      pulse=pulse+1
      if(pulse.gt.pulses)go to 1000
      jout=0
      jprint=-1
      go to 506
      end if
      if(time+dt.lt.esw)then
c
c   pulse is underway
c
      jprint=0
      go to 506
      end if
      if(time+dt.ge.esw.and.curr.ne.0.)then
c
c   pulse is ending
c
c
c   set current to get exact I2T
c
      if(action.lt.amps**2.*esw)then
		curr=sqrt((amps**2.*esw-action)/dt)
	else
      jprint=1
		curr=0.
      end if
      go to 506
      end if
      if(time.gt.esw)then
c
c   pulse is off
c
      curr=0.
      jprint=0
      end if
506   continue
      if(jprint.eq.0)go to 510
      write(20,*)
      if(jprint.eq.-1)write(20,*)'>>>>>>>>>>>>>>>>>>>>>SOP',pulse
      if(jprint.eq.1)write(20,*)'<<<<<<<<<<<<<<<<<<<<<<EOP',pulse
      if(jprint.eq.1)write(20,*)'Action= ',action
      if(jprint.eq.1)write(20,*)'Tmax= ',tmax
      if(jprint.eq.1)write(20,*)'RISmax= ',rismax,rismax/risi2t
      write(20,*)
      write(20,210)oohms
      write(20,220)(t(i),i=1,n)
      write(20,230)(tw(i),i=1,n)
      jprint=0
510   continue
      time=time+dt
      ctime=ctime+dt
      action=action+curr**2.*dt
      ris=ris+(curr**2.-ris/tauris)*dt
      if(time.gt.rep)then
c
c   period complete
c
      time=0.
      go to 505
      end if
      ohms=0.
c   enter section loop
      do 600 i=1,n
      if(curr.eq.0.)then
		qin=0.
		go to 520
      end if
c
c   current is on
c
c   calculate adiabatic temperature rise using g function
c
c
      g0=6.2949e13+2.0934e14*t(i)-2.8709e11*t(i)**2.
      dg=(curr/csa)**2.*dt
      g=g0+dg
      tnew=0.13608+4.5496e-15*g+5.3309e-32*g**2.
      qin=(tnew-t(i))*hcapcu
520   continue
      qxfr=(t(i)-tw(i))*hxfrdx*dt
      t(i)=t(i)+(qin-qxfr)/hcapcu
      if(jadiabatic.eq.1.and.curr.ne.0)then
		t(i)=tnew
		qxfr=0.
      end if
      r=res*(1.+coeff*(t(i)-20.))
      if(t(i).gt.tmax)tmax=t(i)
      if(ris.gt.rismax)rismax=ris
      tw(i)=tw(i)+qxfr/hcapw
      ohms=ohms+r
      if(i.eq.n)then
c
c   record temperature of last section vs. clock time
c
      if(nccount.eq.ncskip)then
		jcout=jcout+1
		tcout(jcout)=ctime
		tcend(jcout)=t(i)
		risout(jcout)=ris
		nccount=1
	else
		nccount=nccount+1
      end if
      end if
      if(i.eq.n.and.pulse.eq.pulses)then
c
c   last pulse, record temperatures vs. time
c
      if(ncount.eq.nskip)then
		jout=jout+1
		tout(jout)=time
		tmid(jout)=t(nmid)
		tend(jout)=t(i)
		tstart(jout)=t(1)
		ncount=1
	else
		ncount=ncount+1
      end if
      end if
600   continue
      oohms=ohms*1.e3
c   determine if water parcel tw(i) should be moved out of x(i)
      dist=dist+v*dt
      if(dist.lt.dx)go to 620
      l=n-1
      do 610 j=1,l
      k=n-j+1
      tw(k)=tw(k-1)
610   continue
      tw(1)=tw0
      dist=0.
620   continue
      go to 505
1000  continue
      do 2000 i=1,jout
      write(30,*)tout(i),tstart(i),tmid(i),tend(i)
2000  continue
      do 3000 i=1,jcout
      write(40,*)tcout(i),tcend(i)
      write(50,*)tcout(i),tcend(i)/tmax,risout(i)/risi2t
3000  continue
      stop
      end