 
*    rcft i=simp,go,lib=(alibcray,fortlib)
c civic o=i,i=simp,lib=(alibcray,stacklib,fortlib),go
Cdc      program h2d (outsi,tape3=outsi)
      program h2d
      common /kls/ kk,ll,debug,version,wher,when,p1d6,pie,igen,p1d2
     x ,dtc,kc,lc,dten,ken,len,ske,hn,siel,cnn,enc,enh,encg,wn
     x    ,ncp
      common /progg/ r0,z0,r1,z1,rp,zp,rr,zz
      common /comn/  r(203,183),z(203,183),u(203,183),rho(203,183)
     x    ,e(203,183),p(203,183),aj(203,183),s(203,183),nbc(203,183)
     x     ,q(203,183),w(203,183),temp(203,183)
     x , a(203,183),b(203,183),cc(203,183),dum(203,183),cbb(203,183)
     x , dbb(203,183),cap(203,183),sig(203,183),ts(203,183)
      common /param/ nycl,tnup,dtnup,dtn,dtnph,dtnmh,edtime,eddt
     x        ,gam,gamz,c0f,c1f,c1,tmax,dtmax,dtmin,tflr,nohyd
     x        ,c2,p2,p3,no,ntty,ned
      common /klspace/ kmn,lmn,kmx,lmx,kmxz,lmxz,kmnp,lmnp,kmxp,lmxp
      common /gencom/ rho0,e0,u0,p0,w0,dr,dz,nbcu,nbcd,nbcl,nbcr
     x ,pb(3),pbb(3),qb(3)
      common /minmax/ xmin,xmax,ymin,ymax,pmin,pmax,qmin,qmax
     x,rmin,rmax,kq,lq,kr,lr,kp,lp
     x ,xminx,xmaxx,yminx,ymaxx
      common /timing/ nbt(40),nct(40),net(40),npt(40),nxt(40),buf(512)
      common /eoscom/ keos,targ1,targ2,targ3,rarg1,rarg2,rarg3,
     x    func1,func2,func3,temps,eps,eps0

      common /com2/ ntsv(2),nrsv(2),msv(2),tes(7),res(9)
     x ,aes(12),bes(12),ces(12),des(12),ees(12),fes(12),ges(12)
     x ,hes(12),pes(12),ites(3),ires(3),izes(3)
       common /tempsr/ v1(203),v2(203),v3(203),v4(203),v5(203),v6(203)
     x ,v7(203),v8(203),v9(203),va(203),vb(203),vc(203)
c ncyl = cycle counter edtime= time ot edit
c tnup = problem time(n+1) eddt = deltat next edit
c dtn = deltat (n) tmax = maximum time
c dtnph= deltat (n+1/2) dtmax = maximum allowed dt
c dtnmh= deltat (n-1/2) dtmin = minimum allowed dt
      dimension kmxprob(7),lmxprob(7),ncylprb(7)
      data kmxprob/ 62, 62, 12,102,202,202,0/
      data lmxprob/ 12, 27, 62, 52 , 5,182,0/
c************************************************************
c     to vary running time change ncylprb, normal is:
c      data ncylprb/800,1000,800,3300,2000,10,0/
c*************************************************************
      data ncylprb/800,1000,800,3300,2000,100,0/
cdc   data version /1./
      version = 1.0
cdc   data ncp/100/
      ncp = 100
CAP   DATA IER/0/
      DATA IER=0
cdc   data ntty/59/
      ntty = 59
cdc   data no/3/
      no = 3
cdc   data debug/0./
      debug =    0.0
cdc   data dtmax/.02/
      dtmax =   0.02
cdc   data dtmin /.0001/
      dtmin =    0.0001
cdc   data tflr/.0001/
       tflr =  0.0001
cdc    data nohyd/0/
      nohyd = 0
cdc   data pie/3.1415926535898/
         pie = 3.1415926535898
cdc   data edtime/0./
      edtime = 0.0
cdc   data eddt/4./
      eddt = 4.0
cdc   data p1d2/.5/
      p1d2 = 0.5
cdc   data tmax/12.001/
      tmax = 12.001
CAP   DATA VCUT/1.E-10/
      DATA VCUT=1.E-10
cdc    data dten/1.e+10/
       dten = 1.e+10
cdc    data dtc/1.e+10/
      dtc = 1.e+10
clll. common dummy(2000)
clll. iq8qdsp=64*.loc.dummy(1)
      time00 = second(0)
c...  call timer(7,'simp0',buf,512,'timsimp',1)
cdc   call link('unit59=terminal,unit3=(outsi,create,text)//')
      open(unit=3,status='NEW')
c     call dropfile(0)
c     call create(3,"hs",3,-1)
c     call fr80id(3hh2d,1)
c     call frame(1)
      sumt = 0
      iprob=1
 1000 if(kmxprob(iprob).eq.0) go to 1001
      do 9 i=1,40
      nbt(i)=0
      nct(i)=0
      net(i)=0
      npt(i)=0
      nxt(i)=0
  9   continue
c
      nycl=0
      tnup=0.
c zero to all arrays
      do 10 l=1,183
      do 10 k=1,203
      r(k,l)=0.
      z(k,l)=0.
      u(k,l)=0.
      rho(k,l)=0.
      q(k,l)=0.
      e(k,l)=0.
      p(k,l)=0.
      aj(k,l)=0.
      s(k,l)=0.
      nbc(k,l)=0
      w(k,l)=0.
      temp(k,l)=0.
      a(k,l)=0.
      b(k,l)=0.
      cc(k,l)=0.
      dum(k,l)=0.
      cbb(k,l)=0.
      dbb(k,l)=0.
      cap(k,l)=0.
      sig(k,l)=0.
      ts(k,l)=0.
  10  continue
c ............ set all arrays to zero
c      call setup
c parameters for test problem 1
      gam=1.4
      p2=6.
      p3=0.
      pb(1)=1.
      pb(2)=0.
      pb(3)=0.
      qb(1)=1.
      qb(2)=0.
      qb(3)=0.
      pbb(1)=0.
      pbb(2)=p2
      pbb(3)=p3
      rho0=1.4
      kmn=2
      lmn=2
      nmax=ncylprb(iprob)
      kmx=kmxprob(iprob)
      lmx=lmxprob(iprob)
      write(no,622) kmx,lmx
c     write(ntty,622) kmx,lmx
  622 format(12h new problem,5h kmx ,i3,5h lmx ,i3)
      iprob=iprob+1
c     read(ntty,5)kmn,kmx,lmn,lmx,eddt,edtime,tmax
c 5   format(4i2,3f5.2)
      dr=1.
      dz=1.
      dtnph=.01
c tmax=10.
      dtn=dtnph
      dtnmh=dtnph
      tnup=0.
      e0=.1
      u0=0.
      p0=.6667*e0*rho0
      w0=0.
      nbcu=1
      nbcr=2
      nbcl=1
      nbcd=1
      c2=1.5
      c0f=c2*.25
      c1=.5
      c1f=.5*c1
      gamz=gam-1.
      kmnp=kmn+1
      lmnp=lmn+1
      kmxp=kmx+1
      kmxz=kmx-1
      lmxp=lmx+1
      lmxz=lmx-1
      dtc=100.
constants
      call gen
call generator
      igen=0
      ned=1
      p1d8=1./8.
      hn=0.
      wn=0.
      i2=necond(i1)
c initialize timer
c
  1   continue
c
c start cycle here
c
      dtc2=1.e+12
      ske=0.
      enc=1.e-20
      enh=1.e-20
      dten=1.e+12
c geometry calculation for boundary zones**********************
      l=lmn
      do 200 k=2,kmxz
      v1(k)=(2.*(z(k+1,l)-z(k,l)))/((r(k+1,l)-r(k,l))**2+
     x     (z(k+1,l)-z(k,l))**2)
      v2(k)=1.-(z(k+1,l)-z(k,l))*v1(k)
      v3(k)=(r(k+1,l)-r(k,l))*v1(k)
      r(k,l-1)=r(k,l)+(r(k,l+1)-r(k,l))*v2(k)+ (z(k,l+1)-z(k,l))*v3(k)
      z(k,l-1)=z(k,l)+(r(k,l+1)-r(k,l))*v3(k)- (z(k,l+1)-z(k,l))*v2(k)
 200  continue
c     1(k-1,l)
c     0(k,l) p(k,l+1)
c     r(k,l-1)
c
      k=kmx
c r0 from above
c z0 from above
      r1=r(k-1,l)
      z1=z(k-1,l)
      rp=r(k,l+1)
      zp=z(k,l+1)
      call projct
      r(k,l-1)=rr
      z(k,l-1)=zz
c     p(k,l+1)
c     1(k-1,l) 0(k,l)
c     r(k,l-1)
c
      l=lmx
c upper side
      do 204 k=2,kmxz
      v1(k)=(2.*(z(k+1,l)-z(k,l)))/((r(k+1,l)-r(k,l))**2+
     x     (z(k+1,l)-z(k,l))**2)
      v2(k)=1.-(z(k+1,l)-z(k,l))*v1(k)
      v3(k)=(r(k+1,l)-r(k,l))*v1(k)
      r(k,l+1)=r(k,l)+(r(k,l-1)-r(k,l))*v2(k)+ (z(k,l-1)-z(k,l))*v3(k)
      z(k,l+1)=z(k,l)+(r(k,l-1)-r(k,l))*v3(k)- (z(k,l-1)-z(k,l))*v2(k)
c     r(k,l+1)
c     0(k,l) 1(k+1,l)
c     p(k,l-1)
c
 204  continue
      k=kmx
      l=lmx
      r0=r(k,l)
      z0=z(k,l)
      r1=r(k-1,l)
      z1=z(k-1,l)
      rp=r(k,l-1)
      zp=z(k,l-1)
      call projct
      r(k,l+1)=rr
      z(k,l+1)=zz
c     r(k,l+1)
c     1(k-1,l) 0(k,l)
c     p(k,l-1)
c
      k=kmn
c left side
      r0=r(k,lmn)
      z0=z(k,lmn)
      do 207 l=lmn,lmxz
      r1=r(k,l+1)
      z1=z(k,l+1)
      rp=r(k+1,l)
      zp=z(k+1,l)
      call projct
      r(k-1,l)=rr
      z(k-1,l)=zz
c     1(k,l+1)
c     r(k-1,l) 0(k,l) p(k+1,l)
c
      r0=r1
      z0=z1
 207  continue
      k=kmn
      l=lmx
      r0=r(k,l)
      z0=z(k,l)
c upper left corner (kmn,lmx)
      r1=r(k,l-1)
      z1=z(k,l-1)
      rp=r(k+1,l)
      zp=z(k+1,l)
      call projct
      r(k-1,l)=rr
c
c     r(k-1,l) 0(k,l) p(k+1,l)
c     1(k,l-1)
      z(k-1,l)=zz
      k=kmx
c right side
      r0=r(k,lmn)
      z0=z(k,lmn)
      do 210 l=lmn,lmxz
      r1=r(k,l+1)
      z1=z(k,l+1)
      rp=r(k-1,l)
      zp=z(k-1,l)
      call projct
      r(k+1,l)=rr
      z(k+1,l)=zz
c
c     1(k,l+1)
c     p(k-1,l) 0(k,l) r(k+1,l)
c
      r0=r1
      z0=z1
 210  continue
      k=kmx
      l=lmx
c upper right corner
      r0=r(k,l)
      z0=z(k,l)
      r1=r(k,l-1)
      z1=z(k,l-1)
      rp=r(k-1,l)
      zp=z(k-1,l)
c
c     p(k-1,l) 0(k,l) r(k+1,l)
c     1(k,l-1)
      call projct
      r(k+1,l)=rr
      z(k+1,l)=zz
      r1=r(k,l+1)
      z1=z(k,l+1)
      rp=r(k-1,l+1)
      zp=z(k-1,l+1)
      call projct
      r(k+1,l+1)=rr
      z(k+1,l+1)=zz
c p(k-1,l+1) 1(k,l+1)r(k+1,l+1)
c     0(k,l)
c upper right corner (kmx,lmx)
      k=kmn
      l=lmn
      r0=r(k,l)
      z0=z(k,l)
      r1=r(k,l-1)
      z1=z(k,l-1)
      rp=r(k+1,l-1)
      zp=z(k+1,l-1)
      call projct
      r(k-1,l-1)=rr
      z(k-1,l-1)=zz
c lower left corner (kmn,lmn)
c     0(k,l)
c     r(k-1,l-1) 1(k,l-1) p(k+1,l-1)
c
      k=kmx
      r0=r(k,l)
      z0=z(k,l)
      r1=r(k+1,l)
      z1=z(k+1,l)
      rp=r(k+1,l+1)
      zp=z(k+1,l+1)
      call projct
      r(k+1,l-1)=rr
      z(k+1,l-1)=zz
c lower right corner (kmx,lmn)
c     r(k+1,l+1)
c     0(k,l) 1(k+1,l)
c     r(k+1,l-1)
      l=lmx
      k=kmn
      r0=r(k,l)
      z0=z(k,l)
      r1=r(k-1,l)
      z1=z(k-1,l)
      rp=r(k-1,l-1)
      zp=z(k-1,l-1)
      call projct
      r(k-1,l+1)=rr
      z(k-1,l+1)=zz
c upper left corner (kmn,lmx)
c
c     i(k-1,l) 0(k,l)
c     p(k-1,l-1) r(k-1,l+1)
c***********************************
c boundary zone set-up
c***********************************
      l=lmn
c lower side
      do 255 k=kmnp,kmx
      rho(k,l)=rho(k,l+1)
      aj(k,l)=aj(k,l+1)
      ip=nbc(k-1,l)
      q(k,l)=qb(ip)*q(k,l+1)
      p(k,l)=pbb(ip)+pb(ip)*p(k,l+1)
 255  continue
c
c     (k,l) = (k,l+1)
c
      k=kmx
c right side
      do 265 l=lmnp,lmx
      rho(k+1,l)=rho(k,l)
      aj(k+1,l)=aj(k,l)
      ip=nbc(k,l)
      q(k+1,l)=qb(ip)*q(k,l)
      p(k+1,l)=pbb(ip)+pb(ip)*p(k,l)
 265  continue
c
c     (k+1,l)=(k,l)
c
      l=lmx
c upper side
      do 275 k=kmnp,kmx
      rho(k,l+1)=rho(k,l)
      aj(k,l+1)=aj(k,l)
      ip=nbc(k-1,l)
      q(k,l+1)=qb(ip)*q(k,l)
      p(k,l+1)=pbb(ip)+pb(ip)*p(k,l)
 275  continue
c
c     (k,l+1)=(k,l)
c
      k=kmn
      do 285 l=lmnp,lmx
      rho(k,l)=rho(k+1,l)
      aj(k,l)=aj(k+1,l)
      ip=nbc(k,l-1)
      q(k,l)=qb(ip)*q(k+1,l)
      p(k,l)=pbb(ip)+pb(ip)*p(k+1,l)
  285 continue
c now set-up corners
c lower left corner
      p(kmn,lmn)=p(kmnp,lmnp)
      q(kmn,lmn)=q(kmnp,lmnp)
      rho(kmn,lmn)=rho(kmnp,lmnp)
      aj(kmn,lmn)=aj(kmnp,lmnp)
c lower right corner
      p(kmxp,lmn)=p(kmxp,lmn+1)
      q(kmxp,lmn)=q(kmxp,lmn+1)
      rho(kmxp,lmn)=rho(kmxp,lmn+1)
      aj(kmxp,lmn)=aj(kmxp,lmn+1)
c upper right corner
      p(kmxp,lmxp)=p(kmxp,lmx)
      q(kmxp,lmxp)=q(kmxp,lmx)
      rho(kmxp,lmxp)=rho(kmxp,lmx)
      aj(kmxp,lmxp)=aj(kmxp,lmx)
c upper left corner
      p(kmn,lmxp)=p(kmnp,lmxp)
      q(kmn,lmxp)=q(kmnp,lmxp)
      rho(kmn,lmxp)=rho(kmnp,lmxp)
      aj(kmn,lmxp)=aj(kmnp,lmxp)
c
c all boundary conditions completed  end of first pass
c
      i2=necond(i1)
      nbt(ned)=nbt(ned)+i2
c bdy condition compute time
      if(debug.eq.0.) go to 442
      igen=1
      write(no,441)
  441 format(10h debug   1)
      call edit
  442 continue
c
c form    p+q in array a and  rho*area in b
c
       do 701 l=lmn,lmxp
       do 701 k=2,kmxp
      a(k,l)=p(k,l)+q(k,l)
      b(k,l)=rho(k,l)*aj(k,l)
  701  continue
c
      dtn2=dtn+dtn
      do 450 l=lmn,lmx
       do 702  k=2,kmx
      v1(k)=a(k,l)*(z(k,l-1)-z(k-1,l))     +
     x        a(k+1,l)*(z(k+1,l)-z(k,l-1)) +
     x        a(k+1,l+1)*(z(k,l+1)-z(k+1,l)) +
     x        a(k,l+1)*(z(k-1,l)-z(k,l+1))
      v2(k)=a(k,l)*(r(k,l-1)-r(k-1,l))     +
     x        a(k+1,l)*(r(k+1,l)-r(k,l-1)) +
     x        a(k+1,l+1)*(r(k,l+1)-r(k+1,l)) +
     x        a(k,l+1)*(r(k-1,l)-r(k,l+1))
      v3(k)=dtn2/(b(k,l)+b(k+1,l)+b(k,l+1)+b(k+1,l+1))
      u(k,l)=u(k,l)-v1(k)*v3(k)
      w(k,l)=w(k,l)+v2(k)*v3(k)
      v4(k)=abs(u(k,l))
      v5(k)=abs(w(k,l))
  702  continue
c
       do 703 k=2,kmx
      if(v4(k).lt.vcut) u(k,l)=0.
      if(v5(k).lt.vcut) w(k,l)=0.
  703  continue
c
  450 continue
      if(nohyd.eq.1) go to 455
c nohyd=1 to skip hydro
      do 452 l=lmn,lmx
      do 451 k=2,kmx
      r(k,l)=r(k,l)+dtnph*u(k,l)
      z(k,l)=z(k,l)+dtnph*w(k,l)
      a(k,l)=u(k,l)*u(k,l)+w(k,l)*w(k,l)
  451 continue
  452 continue
c     acceleration velocity and
c     co-ordinates done           end of second pass
c
c               begin loop 3
c tempry *********
      if(debug.eq.0.) go to 455
      igen=1
      write(no,456)
  456 format(10h debug 2   )
      call edit
c debug print for loop 2
  455 continue
      call hwork
c compute hydro work on the boundary
c
      do 489 k=3,kmx
      vc(k)=0.
  489 continue
c          vc accumulates k.e. for a l row
      do 490 l=lmnp,lmx
c
       do 705 k=3,kmx
      v1(k)=z(k-1,l)-z(k,l-1)
      v2(k)=z(k,l)-z(k-1,l)
      v3(k)=z(k,l)-z(k,l-1)
      v4(k)=z(k-1,l-1)-z(k,l-1)
      v5(k)=z(k-1,l)-z(k-1,l-1)
      v6(k)=v2(k)-v4(k)
      v7(k)=v3(k)+v5(k)
c v6=dzk  v7=dzl
      v9(k)=r(k,l)*v1(k)-r(k-1,l)*v3(k)+r(k,l-1)*v2(k)
      va(k)=r(k-1,l)*v4(k)-r(k-1,l-1)*v1(k)+r(k,l-1)*v5(k)
c v9=aj1 va=aj3
      aj(k,l)=.5*(v9(k)+va(k))
      vb(k)=s(k,l)*rho(k,l)
c v(n)
      s(k,l)=p1d6*((r(k,l)+r(k-1,l)+r(k,l-1))*v9(k)  +
     x           (r(k-1,l)+r(k-1,l-1)+r(k,l-1))*va(k) )
c v4=drk   v5=drl
      v4(k)=r(k,l)-r(k-1,l-1)+r(k,l-1)-r(k-1,l)
      v5(k)=r(k,l)-r(k-1,l-1)-r(k,l-1)+r(k-1,l)
c v8=duk  v9=dul
      v8(k)=u(k,l)-u(k-1,l-1)+u(k,l-1)-u(k-1,l)
      v9(k)=u(k,l)-u(k-1,l-1)-u(k,l-1)+u(k-1,l)
c v3=dwk  v2=dwl
      v3(k)=w(k,l)-w(k-1,l-1)+w(k,l-1)-w(k-1,l)
      v2(k)=w(k,l)-w(k-1,l-1)-w(k,l-1)+w(k-1,l)
      v1(k)=amin1( v4(k)*v2(k)-v6(k)*v9(k),0.)
      va(k)=amin1( v8(k)*v7(k)-v3(k)*v5(k),0.)
c v1=w1     va=w2
      v8(k)=(v1(k)*v1(k))/(v4(k)*v4(k)+v6(k)*v6(k))
      v9(k)=(va(k)*va(k))/(v5(k)*v5(k)+v7(k)*v7(k))
c v8=w3 v9=w4
      v3(k)=1./rho(k,l)
      rho(k,l)=vb(k)/s(k,l)
      v2(k)=(1./rho(k,l))-v3(k)
c v2=delv
      dum(k,l)=rho(k,l)*s(k,l)
  705  continue
c
c    v1=w1      v8=w3      v2=delv   v4=drk  v6=dzk
c    va=w2      v9=w4            v5=drl  v7=dzl
c
c
c
       do 707 k=3,kmx
      v3(k)=sqrt(gam*p(k,l)/rho(k,l))
c                         v3=ca=sound speed
      v1(k)=v8(k)+v9(k)
      q(k,l)=c0f*rho(k,l)*v1(k)+c1f*v3(k)*rho(k,l)*sqrt(v1(k))
      vb(k)=v3(k)+1.e-20
      va(k)=(aj(k,l)*aj(k,l))/(vb(k)*vb(k)*
     x         ( v4(k)*v4(k)+v5(k)*v5(k)+v6(k)*v6(k)+v7(k)*v7(k)))
  707  continue
c
       do 709 k=3,kmx
      if(vb(k).eq.1.e-20) go to 709
      if(dtc2.le.va(k))     go to 709
c  find minimum dt
      dtc2=va(k)
      kc=k
      lc=l
  709  continue
c      v2=delv
c             now - energy update
       do 715 k=3,kmx
      v4(k)=e(k,l)-(p(k,l)+q(k,l))*v2(k)
c e(n+1)
      v7(k)=.6667*v4(k)*rho(k,l)
c p(n+1/2)
      e(k,l)=amax1( e(k,l)-(.5*(v7(k)+p(k,l))+q(k,l))*v2(k),1.e-30)
c
      vc(k)=vc(k)+dum(k,l)*(a(k,l)+a(k-1,l)+a(k,l-1)+a(k-1,l-1))*p1d8
      temp(k,l)=10.*e(k,l)
      p(k,l)=.06667*temp(k,l)*rho(k,l)
c advance theta and pressure in time
  715 continue
c
c
  490 continue
c
      do 491 k=3,kmx
      ske=ske+vc(k)
  491 continue
c                 kinetic energy for problem
c********************************end of loop 3
      if(debug.eq.0.) go to 495
c********** debug print loop-3**********
cdc   write(no,493)
cdc  493 format(10h debug    3  )
      igen=1
      call edit
c debug print loop-3
  495 continue
      igen=0
      i2=necond(i1)
      nxt(ned)=nxt(ned)+i2
c
      call conduct
c             do heat conduction
c
      nycl=nycl+1
c advance cycle counter
      dtnmh=dtnph
      dtc=sqrt(dtc2)
      dtnph=dtc
      dtnph=amin1(dtnph,dten,dtmax)
c limit magnitude of dt
      dtn=.5*(dtnph+dtnmh)
      tnup=tnup+dtnph
      if(dtnph.ge.dtmin) go to 602
c********** dt is below allowed minimum **********
      write(no,601)nycl,tnup,dtnph,dtmin
      write(no,601)nycl,tnup,dtnph,dtmin
  601 format(12h dtstop nycl,i6,3h t ,e12.4,4h dt ,e12.4)
      go to 999
  602 continue
      te=ske+enh
      cn=te-hn-wn
      if(nycl.eq.1) cnold=cn
      cnn=cn-cnold
      encg=encg+cnn
      cnold=cn
      if(mod(nycl,ncp).ne.0) go to 603
      write(no,776)
  776 format(6h cycle,4x,5htime ,7x,2hdt,10x,3hdtc,5x,8h  kc  lc,
     x         4x,3hdte,5x,8h ken len)
      write(no,777)nycl,tnup,dtnph,dtc,kc,lc,dten,ken,len
  777 format(i6,3e12.4,2i4,e12.4,2i4)
      z2=abs(enc-(enh+hn))/enc
      write(no,708)
  708 format(4x,4hetot,8x,4hie    ,8x,4hke  ,8x,4hhn  ,8x,4hwn  ,
     x         8x,5hecons,7x,5hcn(n),7x,4hecng)
      write(no,779)te,enc,ske,hn,wn,z2,cnn,encg
  779 format(8e12.4)
 603  continue
      i2=necond(i1)
      nct(ned)=nct(ned)+i2
c run time for physics
      if(nycl.ne.nmax) go to 605
c time to edit
      call edit
      write(no,604)nycl,tnup,dtnph,rmax,kr,lr
  604 format(12h edit nycl= ,i6,2e12.4,e14.5,2i4)
c message to tty
c advance editme to next value
  605 if(mod(nycl,ncp).eq.0) ned=ned+1
      if(nycl.ne.nmax) go to 610
c********** problem has reached cmax**********
      write(no,607)nycl,tnup,tmax
      write(no,607)nycl,tnup,tmax
  607 format(12h stop cmax  ,i6,2e12.4)
      go to 999
  610 continue
      go to 1
  999 continue
c problem completed get off
c     call plote
      write(no,615)
  615 format(3x,8hboundary,9x,2hec,7x,4hedit,7x,4hplot,6x,5hhydro)
      write(no,616)(nbt(k),nct(k),net(k),npt(k),nxt(k),k=1,ned)
  616 format(5(1x,i10))
      i1=0
      do 618 k=1,ned
      i1=i1+nbt(k)+nct(k)+net(k)+npt(k)+nxt(k)
  618 continue
c  total for various pieces of code
c
      if(i1.eq.0) i1=1
      do 619 k=2,ned
      nbt(1)=nbt(1)+nbt(k)
      nct(1)=nct(1)+nct(k)
      net(1)=net(1)+net(k)
      npt(1)=npt(1)+npt(k)
      nxt(1)=nxt(1)+nxt(k)
 619  continue
c
      write(no,616) nbt(1),nct(1),net(1),npt(1),nxt(1)
      v1(1)=(nbt(1)*100)/i1
      v1(2)=(nct(1)*100)/i1
      v1(3)=(net(1)*100)/i1
      v1(4)=(npt(1)*100)/i1
      v1(5)=(nxt(1)*100)/i1
      write(no,620) (v1(i2),i2=1,5)
  620 format(2x,7hpercent,1x,5e12.4)
      z1=float(i1)*1.e-6
      z2=float( (kmx-kmn)*(lmx-lmn)*nycl)/amax1(z1,.1)
c     write(ntty,621) z1,z2
      write(no,621) z1,z2
  621 format(//,14h total seconds,e12.4,14h zone-cycles/s,e12.4)
      sumt = sumt+z1
      go to 1000
 1001 continue
      write(no,621) sumt
      overhed = second(0)-time00
      write(no,623) overhed
 623  format('total time (includes uncounted overhead) :',f16.6)
c...  call timend
      stop
      end
      subroutine gen
      common /kls/ kk,ll,debug,version,wher,when,p1d6,pie,igen,p1d2
     x ,dtc,kc,lc,dten,ken,len,ske,hn,siel,cnn,enc,enh,encg,wn
     x    ,ncp
      common /progg/ r0,z0,r1,z1,rp,zp,rr,zz
      common /comn/  r(203,183),z(203,183),u(203,183),rho(203,183)
     x    ,e(203,183),p(203,183),aj(203,183),s(203,183),nbc(203,183)
     x     ,q(203,183),w(203,183),temp(203,183)
     x , a(203,183),b(203,183),cc(203,183),dum(203,183),cbb(203,183)
     x , dbb(203,183),cap(203,183),sig(203,183),ts(203,183)
      common /param/ nycl,tnup,dtnup,dtn,dtnph,dtnmh,edtime,eddt
     x        ,gam,gamz,c0f,c1f,c1,tmax,dtmax,dtmin,tflr,nohyd
     x        ,c2,p2,p3,no,ntty,ned
      common /klspace/ kmn,lmn,kmx,lmx,kmxz,lmxz,kmnp,lmnp,kmxp,lmxp
      common /gencom/ rho0,e0,u0,p0,w0,dr,dz,nbcu,nbcd,nbcl,nbcr
     x ,pb(3),pbb(3),qb(3)
      common /minmax/ xmin,xmax,ymin,ymax,pmin,pmax,qmin,qmax
     x,rmin,rmax,kq,lq,kr,lr,kp,lp
     x ,xminx,xmaxx,yminx,ymaxx
      common /timing/ nbt(40),nct(40),net(40),npt(40),nxt(40)
      common /eoscom/ keos,targ1,targ2,targ3,rarg1,rarg2,rarg3,
     x    func1,func2,func3,temps,eps,eps0

      common /com2/ ntsv(2),nrsv(2),msv(2),tes(7),res(9)
     x ,aes(12),bes(12),ces(12),des(12),ees(12),fes(12),ges(12)
     x ,hes(12),pes(12),ites(3),ires(3),izes(3)
cdc   data rs/1./
cdc   data zs/1./
cdc   data igen/1/
      rs = 1.0
      zs = 1.0
      igen = 1
c generate nbc array
      do 52 k=2,kmx
      nbc(k,lmn)=nbcd
c lower boundary
      nbc(k,lmx)=nbcu
c upper boundary
  52  continue
      do 54 l=lmn,lmx
      nbc(kmn,l)=nbcl
c left boundary
      nbc(kmx,l)=nbcr
c right boundary
  54  continue
c generate coordinates and velocities
c igen.ne.0 causes edit routine to print all variables
      xminx=1.e+6
      yminx=1.e+6
      xmaxx=-1.e+6
      ymaxx=-1.e+6
      rp=lmx-lmn
      zp=kmx-kmn
      do 58 k=2,kmx
      z1=10+k-kmn
      do 57 l=lmn,lmx
      rr=l-2
      zz=(-.5+rr/rp)*pie
      r(k,l)=z1*cos(zz)
      z(k,l)=z1*sin(zz)
      ymaxx=amax1(ymaxx,z(k,l))
      xmaxx=amax1(xmaxx,r(k,l))
      xminx=amin1(xminx,r(k,l))
      yminx=amin1(yminx,z(k,l))
  57  continue
  58  continue
c generate zone quantities rho,p, and compute area
      p1d6=1./6.
      do 65 l=lmnp,lmx
      do 63 k=3,kmx
      rho(k,l)=rho0
      p(k,l)=p0
      e(k,l)=e0
compute jacobian
      aj1=r(k,l)*(z(k-1,l)-z(k,l-1))+r(k-1,l)*(z(k,l-1)-z(k,l))
     x      +r(k,l-1)*(z(k,l)-z(k-1,l))
      aj3=r(k-1,l)*(z(k-1,l-1)-z(k,l-1))+r(k-1,l-1)*(z(k,l-1)-z(k-1,l))
     x      +r(k,l-1)*(z(k-1,l)-z(k-1,l-1))
      aj(k,l)=p1d2*(aj1+aj3)
      s(k,l)=p1d6*((r(k,l)+r(k-1,l)+r(k,l-1))*aj1 +
     x           (r(k,l-1)+r(k-1,l)+r(k-1,l-1))*aj3)
  63  continue
  65  continue
      if(debug.eq.0.) go to 70
c debug edit
c write nbc boundary sentinels
      write(no,81) (nbc(k,lmn),k=2,kmx)
      write(no,82) (nbc(k,lmx),k=2,kmx)
      write(no,83) (nbc(kmn,l),l=lmn,lmx)
      write(no,84) (nbc(kmx,l),l=lmn,lmx)
  81  format(3hlmn,80i1)
  82  format(3hlmx,80i1)
  83  format(3hkmn,80i1)
  84  format(3hkmx,80i1)
      call edit
  70  continue
      write(no,75)
  75  format(22h generation completed      )
c     write(ntty,75)
      return
      end
      subroutine edit
      common /kls/ kk,ll,debug,version,wher,when,p1d6,pie,igen,p1d2
     x ,dtc,kc,lc,dten,ken,len,ske,hn,siel,cnn,enc,enh,encg,wn
     x    ,ncp
      common /progg/ r0,z0,r1,z1,rp,zp,rr,zz
      common /comn/  r(203,183),z(203,183),u(203,183),rho(203,183)
     x    ,e(203,183),p(203,183),aj(203,183),s(203,183),nbc(203,183)
     x     ,q(203,183),w(203,183),temp(203,183)
     x , a(203,183),b(203,183),cc(203,183),dum(203,183),cbb(203,183)
     x , dbb(203,183),cap(203,183),sig(203,183),ts(203,183)
      common /param/ nycl,tnup,dtnup,dtn,dtnph,dtnmh,edtime,eddt
     x        ,gam,gamz,c0f,c1f,c1,tmax,dtmax,dtmin,tflr,nohyd
     x        ,c2,p2,p3,no,ntty,ned
      common /klspace/ kmn,lmn,kmx,lmx,kmxz,lmxz,kmnp,lmnp,kmxp,lmxp
      common /gencom/ rho0,e0,u0,p0,w0,dr,dz,nbcu,nbcd,nbcl,nbcr
     x ,pb(3),pbb(3),qb(3)
      common /minmax/ xmin,xmax,ymin,ymax,pmin,pmax,qmin,qmax
     x,rmin,rmax,kq,lq,kr,lr,kp,lp
     x ,xminx,xmaxx,yminx,ymaxx
      common /timing/ nbt(40),nct(40),net(40),npt(40),nxt(40)
      common /eoscom/ keos,targ1,targ2,targ3,rarg1,rarg2,rarg3,
     x    func1,func2,func3,temps,eps,eps0

      common /com2/ ntsv(2),nrsv(2),msv(2),tes(7),res(9)
     x ,aes(12),bes(12),ces(12),des(12),ees(12),fes(12),ges(12)
     x ,hes(12),pes(12),ites(3),ires(3),izes(3)
c tempis subroutine edits all mesh variables
       data n100/100/
      xmin=1.e+6
      ymin=1.e+6
      ymax=1.e-6
      xmax=1.e-6
c find minimum and maximum p,q,rho
      pmin=1.e+6
      rmin=1.e+6
      qmin=1.e+6
      qmax=-1.e+6
      rmax=-1.e+6
      pmax=-1.e+6
      kq=0
      lq=0
      kr=0
      lr=0
      kp=0
      lp=0
      do 715 l=lmnp,lmx
      do 714 k=3,kmx
      if(rho(k,l).le.rmax) go to 701
      rmax=rho(k,l)
      kr=k
      lr=l
  701 if(p(k,l).le.pmax) go to 702
      pmax=p(k,l)
      kp=k
      lp=l
  702 if(q(k,l).le.qmax) go to 703
      qmax=q(k,l)
      kq=k
      lq=l
 703  pmin=amin1(pmin,p(k,l))
      qmin=amin1(qmin,q(k,l))
      rmin=amin1(rmin,rho(k,l))
      xmin=amin1(xmin,r(k,l))
      xmax=amax1(xmax,r(k,l))
      ymin=amin1(ymin,z(k,l))
      ymax=amax1(ymax,z(k,l))
c as used for sound speed
 714  continue
 715  continue
c
c end of loop searching for maximum values
c print problem parameters
      write(no,717) nycl,tnup,dtnph,dtn,version,wher,when
  717 format(6h nycl ,i6,6h time ,e12.4,7h dtnph ,e12.4,5h dtn ,
     x e12.4,9h version ,f4.1,2a10)
      write(no,718) pmax,kp,lp,qmax,kq,lq,rmax,kr,lr
  718 format(14h maximum (k,l),e12.4,2i4,3h p ,e12.4,2i4,3h q ,
     x e12.4,2i4,5h rho )
      uvtest=1.e-5
      kl=kmn
      ll=lmn
      ku=kmx
      lu=lmn+2
      if(igen.eq.0) go to 720
       uvtest=-100.
c print all mesh points
c igen.ne.0 will result in edit of entire mesh,=0 only active zones
      kl=kmn-1
      ll=lmn-1
      ku=kmx+1
      lu=lmx+1
 720  continue
c begin edit
      do 740 l=ll,lu
      write(no,725)
  725 format(8h   l   k,4x,1hr,10x,1hz,10x,1hu,10x,1hw,10x,3hrho,
     x         8x,1he,10x,1hp,10x,1hq,10x,2haj,9x,5htheta)
      do 738 k=kl,ku
      if((abs(u(k,l))+abs(w(k,l))).le.uvtest)go to 738
c dont print variables if no motion
      write(no,726)l,k,r(k,l),z(k,l),u(k,l),w(k,l),rho(k,l),e(k,l)
     x,p(k,l),q(k,l),aj(k,l),temp(k,l)
 726  format(2i4,10e11.3)
 738  continue
 740  continue
      net(ned)=necond(i1)
c     call maps(xminx,xmaxx,yminx,ymaxx)
c     call setch(1.,10.,0,0,1,1,0)
c     write(n100,751)nycl,tnup
c 751 format(5h2-dh ,i6,e12.4)
c     do 760 l=lmn,lmxz
c     do 760 k=2,kmxz
c     call line(r(k,l),z(k,l),r(k+1,l),z(k+1,l))
c     call line(r(k,l),z(k,l),r(k,l+1),z(k,l+1))
c 760  continue
c     l=lmx
c     do 761 k=2,kmxz
c     call line(r(k,l),z(k,l),r(k+1,l),z(k+1,l))
c 761 continue
c     k=kmx
c     do 762 l=lmn,lmxz
c     call line(r(k,l),z(k,l),r(k,l+1),z(k,l+1))
c 762 continue
c     call frame(1)
c     call maps(xminx,xmaxx,yminx,ymaxx)
c      do 765 l=lmn,lmx
c      do 765 k=2,kmx
c     ra=r(k,l)+u(k,l)
c     za=z(k,l)+w(k,l)
c     call plotv(r(k,l),z(k,l),ra,za)
c 765 continue
c     call frame(1)
      npt(ned)=necond(i1)
      if(ned.gt.40) ned=1
      return
      end
      subroutine projct
      common /progg/r0,z0,r1,z1,rp,zp,rr,zz
c (r0,z0) and (r1,z1) are boundary points
c reflect (rp,zp) to (rr,zz)
c
      ww=(2.*(z1-z0))/((r1-r0)**2+(z1-z0)**2)
      alp=1.-(z1-z0)*ww
      bet=(r1-r0)*ww
      rr=r0+(rp-r0)*alp + (zp-z0)*bet
      zz=z0+(rp-r0)*bet - (zp-z0)*alp
      return
      end
      subroutine conduct
      common /kls/ kk,ll,debug,version,wher,when,p1d6,pie,igen,p1d2
     x ,dtc,kc,lc,dten,ken,len,ske,hn,siel,cnn,enc,enh,encg,wn
     x    ,ncp
      common /progg/ r0,z0,r1,z1,rp,zp,rr,zz
      common /comn/  r(203,183),z(203,183),u(203,183),rho(203,183)
     x    ,e(203,183),p(203,183),aj(203,183),s(203,183),nbc(203,183)
     x     ,q(203,183),w(203,183),temp(203,183)
     x , a(203,183),b(203,183),cc(203,183),dum(203,183),cbb(203,183)
     x , dbb(203,183),cap(203,183),sig(203,183),ts(203,183)
      common /param/ nycl,tnup,dtnup,dtn,dtnph,dtnmh,edtime,eddt
     x        ,gam,gamz,c0f,c1f,c1,tmax,dtmax,dtmin,tflr,nohyd
     x        ,c2,p2,p3,no,ntty,ned
      common /klspace/ kmn,lmn,kmx,lmx,kmxz,lmxz,kmnp,lmnp,kmxp,lmxp
      common /gencom/ rho0,e0,u0,p0,w0,dr,dz,nbcu,nbcd,nbcl,nbcr
     x ,pb(3),pbb(3),qb(3)
      common /minmax/ xmin,xmax,ymin,ymax,pmin,pmax,qmin,qmax
     x,rmin,rmax,kq,lq,kr,lr,kp,lp
     x ,xminx,xmaxx,yminx,ymaxx
      common /timing/ nbt(40),nct(40),net(40),npt(40),nxt(40)
      common /eoscom/ keos,targ1,targ2,targ3,rarg1,rarg2,rarg3,
     x    func1,func2,func3,temps,eps,eps0

      common /com2/ ntsv(2),nrsv(2),msv(2),tes(7),res(9)
     x ,aes(12),bes(12),ces(12),des(12),ees(12),fes(12),ges(12)
     x ,hes(12),pes(12),ites(3),ires(3),izes(3)
       common /tempsr/ v1(203),v2(203),v3(203),v4(203),v5(203),v6(203)
     x ,v7(203),v8(203),v9(203),va(203),vb(203),vc(203)
c electron conduction -lu-
c
      do 10 l=lmn,lmx
      do 10 k=2,kmx
      cap(k,l)=.1
      cc(k,l)=(.0001*sqrt(temp(k,l))*temp(k,l)**2)/aj(k,l)
      sig(k,l)=dum(k,l)*cap(k,l)/dtnph
      ts(k,l)=temp(k,l)
 10   continue
c
      do 12 l=lmn,lmx
       do 12 k=2,kmxz
      cbb(k,l)=(2.*cc(k+1,l)*cc(k+1,l+1))/(cc(k+1,l)+cc(k+1,l+1))
     x * (.5*(r(k,l)+r(k+1,l))*((r(k+1,l)-r(k,l))**2
     x      +(z(k+1,l)-z(k,l))**2) )
  12  continue
      do 14 l=lmn,lmxz
      do 14 k=2,kmx
      dbb(k,l)=(2.*cc(k+1,l+1)*cc(k,l+1))/(cc(k+1,l+1)+cc(k,l+1))
     x * (.5*(r(k,l)+r(k,l+1))*((r(k,l+1)-r(k,l))**2
     x      +(z(k,l+1)-z(k,l))**2) )
  14  continue
c
c boundary conditions
      do 17 l=lmn,lmx
      a(kmn,l)=0.
      b(kmn,l)=temp(kmn,l)
      dbb(kmn,l)=0.
  17  continue
c
      do 18 k=2,kmx
  18  cbb(k,lmn)=0.
      do 19 k=2,kmx
      a(k,lmn)=0.
      b(k,lmn)=temp(k,lmn)
      cbb(k,lmx)=0.
  19  continue
      if(lmx.eq.5) go to 45
c if this is a 1-d problem dont do adi
c .................................
c
c
c
c.......... z sweep
c
       do  51 l=lmnp,lmx
       do  51 k=3,kmx
       dum(k,l)=sig(k,l)+cbb(k-1,l)+cbb(k-1,l-1)*(1.-a(k,l-1))
       a(k,l)=cbb(k-1,l)/dum(k,l)
       b(k,l)=(sig(k,l)*temp(k,l)+cbb(k-1,l-1)*b(k,l-1)
     x )/dum(k,l)
   51 continue
c...... alpha,beta forward
      ml=lmx+1
      do 52 l=lmnp,lmx
      ml=ml-1
      do 52 k=3,kmx
      temp(k,ml)=a(k,ml)*temp(k,ml+1)+b(k,ml)
   52 continue
c back substitution
c
c.......... z sweep end
c
c.......... r sweep
c
       do  41 k=3,kmx
       do  41 l=lmnp,lmx
       dum(k,l)=sig(k,l)+dbb(k,l-1)+dbb(k-1,l-1)*(1.-a(k-1,l))
       a(k,l)=dbb(k,l-1)/dum(k,l)
      b(k,l)=(sig(k,l)*temp(k,l)+dbb(k-1,l-1)*b(k-1,l)
     x )/dum(k,l)
   41 continue
c...... alpha beta forward sweep
      ml=kmx+1
      do 42 k=3,kmx
      ml=ml-1
      do 42 l=lmnp,lmx
      temp(ml,l)=a(ml,l)*temp(ml+1,l)+b(ml,l)
  42  continue
c back substitution r direction
c
c.......... r sweep end
c compute dt control for heat conduction
c
  45  ye=0.
      ken=0
      len=0
c
      do 125 k=3,kmx
      va(k)=0.
      vb(k)=0.
  125 continue
c           va and vb are used to accumulate i.e. by l row
      do 130 l=lmnp,lmx
      do 170 k=3,kmx
      a(k,l)=rho(k,l)*s(k,l)
      sig(k,l)=abs((temp(k,l)-ts(k,l))/ts(k,l))
      va(k)=va(k)+e(k,l)*a(k,l)
      e(k,l)=amax1(.1*temp(k,l),1.e-30)
      vb(k)=vb(k)+e(k,l)*a(k,l)
      temp(k,l)=ts(k,l)
  170 continue
      do 111 k=3,kmx
c get new energy
      tempr=sig(k,l)
      if(tempr.le.ye) go to 111
      ye=tempr
      ken=k
      len=l
  111 continue
c
  130 continue
c
      do 132 k=3,kmx
      enc=enc+vb(k)
      enh=enh+va(k)
  132 continue
c                internal energy for total grid
      if(ye.eq.0.) go to 118
      dten=(.1*dtnph)/ye
  118 continue
c energy balance hn
      do 122 k=2,kmx
      v1(k)=dtnph*cbb(k-1,lmn)*(temp(k,lmn)-temp(k,lmn+1))
     x     +dtnph*cbb(k-1,lmx)*(temp(k,lmx+1)-temp(k,lmx))
  122 continue
c
      do 124 l=2,lmx
      v2(l)=dtnph*dbb(kmn,l-1)*(temp(kmn,l)-temp(kmn+1,l))
     x     +dtnph*dbb(kmx,l-1)*(temp(kmx+1,l)-temp(kmx,l))
  124 continue
c
      do 126 k=2,kmx
 126  hn=hn+v1(k)
      do 128 l=2,lmx
  128 hn=hn+v2(l)
      return
      end
      function necond(iarg)
      INTEGER IARG
      aa1=second(iarg)
      necond=(aa1-aa2)*1.e+6
      aa2=aa1
      return
      end
C  Not for Sun
C**********************************************
      FUNCTION   SECOND( OLDSEC)
C***********************************************
C
C     REAL SECOND
C/    IMPLICIT  REAL*8(A-H,O-Z)
C
C     SECOND= Cumulative CPU time for job in seconds.  MKS unit is secon
C             Clock resolution should be less than 2% of Kernel 12 run-t
C
C
C     If your system provides a timing routine that satisfies
C     the definition above; then simply delete this function.
C
C     Else this function must be programmed using some
C     timing routine available in your system.
C     Timing routines with CPU clock resolution are always sufficient.
C     Timing routines with microsec. resolution are usually sufficient.
C
C     Timing routines with much less resolution may require the use
C     of multiple-pass loops around each kernel to make the run time
C     at least 50 times the tick-period of the timing routine.
C     In these cases, you must redefine the DATA statements for IPAS1
C     to increase the value of Loop set in subroutine SIZES.
C
C     If no CPU timer is available, then you can time each kernel by
C     the wall clock using the PAUSE statement at the end of subr. VALUE
C
C     Function TICK measures the overhead time for a call to SECOND.
C     An independent calibration of the running time of at least
C     one kernel may be wise.  See Subr. TEST, label 100+.
C
C
C     The following statement is deliberately incomplete:
C
C     SECOND=
C
C
C     The following statement was used on the DEC  VAX/780  VMS 3.0
C     Also set:   Loop= 10*Loop in Subroutine SIZES to run kernels long
C
C         REAL    SECNDS
C*        SECOND= SECNDS( OLDSEC)
C
C OR:
C
C$      DATA  INITIA /237/
C$      IF( INITIA.NE.237 )  GO TO 1
C$      NSTAT= LIB$INIT_TIMER
C$    1 INITIA= INITIA + 1
C$
C$      NSTAT = LIB$STAT_TIMER(2,ISEC)
C$      SECOND= FLOAT(ISEC)*0.01 - OLDSEC
C
C
C     The following statements were used on the DEC PDP-11/23 RT-11 syst
C     Also set:   Loop=100*Loop in Subroutine SIZES to run kernels long
C
C*       DIMENSION JT(2)
C*       CALL GTIM(JT)
C*       TIME1 = JT(1)
C*       TIME2 = JT(2)
C*       TIME = TIME1 * 65768. + TIME2
C*       SECOND=TIME/60. - OLDSEC
C
C
C     The following statements were used on the Hewlett-Packard HP 9000
C     Also set:   Loop= 10*Loop in Subroutine SIZES to run kernels long
C
C*       INTEGER*4 ITIME(4)
C*       CALL TIMES( ITIME(4))
C*       TIMEX= ITIME(1) + ITIME(2) + ITIME(3) + ITIME(4)
C*       SECOND= TIMEX/60. - OLDSEC
C
C
C     The following statements were used on the  SUN workstation
C     Also set:   Loop= 10*Loop in Subroutine SIZES to run kernels long
C
Cdc      DIMENSION XTIME(2)
Cdc      REAL*4 XTIME
Cdc      XT= ETIME( XTIME)
Cdc      SECOND= (XTIME(1) + XTIME(2)) - OLDSEC
C     The following statements were used on the IBM VM/CMS system
C     Also set:   Loop= 100*Loop in Subroutine SIZES to run kernels long
C
C      REAL*8 XISEC,XSEC,XCNT,XMAX,OLDSEC
C      LOGICAL INIT
C      SAVE INIT
C      SAVE XISEC,XCNT,XMAX
C      DATA INIT/.TRUE./
C      CALL CLOCKX(XSEC,XCNT,XMAX)
C      IF (INIT) THEN
C          SECOND = 0.0
C          XISEC = XSEC
C          INIT = .FALSE.
C      ELSE
C          SECOND = (XSEC - XISEC)/1000000.0
C      ENDIF
Cdc   PRINT *,XISEC,XSEC,SECOND
C
      SECOND = 0.0
      RETURN
      END
      subroutine hwork
      common /kls/ kk,ll,debug,version,wher,when,p1d6,pie,igen,p1d2
     x ,dtc,kc,lc,dten,ken,len,ske,hn,siel,cnn,enc,enh,encg,wn
     x    ,ncp
      common /progg/ r0,z0,r1,z1,rp,zp,rr,zz
      common /comn/  r(203,183),z(203,183),u(203,183),rho(203,183)
     x    ,e(203,183),p(203,183),aj(203,183),s(203,183),nbc(203,183)
     x     ,q(203,183),w(203,183),temp(203,183)
     x , a(203,183),b(203,183),cc(203,183),dum(203,183),cbb(203,183)
     x , dbb(203,183),cap(203,183),sig(203,183),ts(203,183)
      common /param/ nycl,tnup,dtnup,dtn,dtnph,dtnmh,edtime,eddt
     x        ,gam,gamz,c0f,c1f,c1,tmax,dtmax,dtmin,tflr,nohyd
     x        ,c2,p2,p3,no,ntty,ned
      common /klspace/ kmn,lmn,kmx,lmx,kmxz,lmxz,kmnp,lmnp,kmxp,lmxp
      common /gencom/ rho0,e0,u0,p0,w0,dr,dz,nbcu,nbcd,nbcl,nbcr
     x ,pb(3),pbb(3),qb(3)
      common /minmax/ xmin,xmax,ymin,ymax,pmin,pmax,qmin,qmax
     x,rmin,rmax,kq,lq,kr,lr,kp,lp
     x ,xminx,xmaxx,yminx,ymaxx
      common /timing/ nbt(40),nct(40),net(40),npt(40),nxt(40)
      common /eoscom/ keos,targ1,targ2,targ3,rarg1,rarg2,rarg3,
     x    func1,func2,func3,temps,eps,eps0

      common /com2/ ntsv(2),nrsv(2),msv(2),tes(7),res(9)
     x ,aes(12),bes(12),ces(12),des(12),ees(12),fes(12),ges(12)
     x ,hes(12),pes(12),ites(3),ires(3),izes(3)
       common /tempsr/ v1(203),v2(203),v3(203),v4(203),v5(203),v6(203)
     x ,v7(203),v8(203),v9(203),va(203),vb(203),vc(203)
c sum the hydro work on the boundary
c
c inputs: p,q r,z u,w
c outputs: wn
      z1=dtnph/8.
       z1m=-z1
      do 510 k=3,kmx
      v1(k)=z1*(p(k,lmn+1)+p(k,lmn)+q(k,lmn+1)+q(k,lmn))
     x    *( (u(k,lmn)+u(k-1,lmn))*(z(k,lmn)-z(k-1,lmn))
     x      -(w(k,lmn)+w(k-1,lmn))*(r(k,lmn)-r(k-1,lmn))
     x     )*(r(k,lmn)+r(k-1,lmn))
      v2(k)=z1m*(p(k,lmx+1)+p(k,lmx)+q(k,lmx+1)+q(k,lmx))
     x    *( (u(k,lmx)+u(k-1,lmx))*(z(k,lmx)-z(k-1,lmx))
     x      -(w(k,lmx)+w(k-1,lmx))*(r(k,lmx)-r(k-1,lmx))
     x     )*(r(k,lmx)+r(k-1,lmx))
      v5(k)=v1(k)+v2(k)
  510  continue
c
      do 515 l=lmnp,lmx
      v3(l)=z1*(p(kmn+1,l)+p(kmn,l)+q(kmn+1,l)+q(kmn,l))
     x    *( (u(kmn,l)+u(kmn,l-1))*(z(kmn,l)-z(kmn,l-1))
     x      -(w(kmn,l)+w(kmn,l-1))*(r(kmn,l)-r(kmn,l-1))
     x     )*(r(kmn,l)+r(kmn,l-1))
      v4(l)=z1m*(p(kmx+1,l)+p(kmx,l)+q(kmx+1,l)+q(kmx,l))
     x    *( (u(kmx,l)+u(kmx,l-1))*(z(kmx,l)-z(kmx,l-1))
     x      -(w(kmx,l)+w(kmx,l-1))*(r(kmx,l)-r(kmx,l-1))
     x     )*(r(kmx,l)+r(kmx,l-1))
      v6(l)=v3(l)+v4(l)
  515  continue
c
      do 517 k=3,kmx
  517 wn=wn+v5(k)
      do 519 l=lmnp,lmx
  519 wn=wn+v6(l)
      return
      end
