program vhone C C-------------------------------------------------------------------------- C C C VV VV HH HH 111 C VV VV HH HH 1111 C VV VV HH HH 11 C VV VV HH HH 11 C VV VV HHHHHHHHHHH == 11 C VV VV HH HH 11 C VV VV HH HH 11 C VVVV HH HH 11 C VV HH HH 111111 C C C VIRGINIA HYDRODYNAMICS #1 C C-------------------------------------------------------------------------- C The Virginia Numerical Bull Session ideal hydrodynamics PPMLR C C Version 1.0 October 1991 C This version was removed from the ``loop'' in the continued development C of VH-1 on Sep 1 1990, but has followed a similar evolution until this C release date. C C Output channels: C 4 = aaaaadaa, binary dump files for restarting C 8 = aaaaahst, history file C aaaaamvy, HDF raster movie file C aaaaa.1000, HDF SDS output C GLOBALS C include 'global.h' include 'zone.h' include 'sweepsize.h' C C LOCALS C external sweepx, sweepy, dtcon, grid, init, dump, images C character dmpfile*8, hstfile*8, mvyfile*8, & prefx*5, rstrt*2, sf1*1, sf2*1, char*1 real endtime,timep,tprin,xmin,xmax,ymin,ymax,zmin,zmax integer ncycp, ncycd, ncycm, ncycend, nprin, ndump, & nmovie, tmpcyc C C Namelist for input deck C namelist / hinput / rstrt, prefx & ,gam, imax, jmax, kmax, xmin, xmax, ymin, ymax & ,zmin, zmax, courant & ,ncycend, ndump, nprin, nmovie, tprin, endtime C C Set some floor values for pressure and density C small = 1.e-35 smallp = 1.e-35 smallr = 1.e-35 C C----------------------------------------------------------------------- C C Begin by reading input deck, and defining file names C open (unit=15,file='indat',status='old',form='formatted') read (15,hinput,err=999,end=999) close(15) hstfile = prefx(1:5) // 'hst' dmpfile = prefx(1:5) // 'daa' mvyfile = prefx(1:5) // 'mvy' C C Check that arrays are large enough for desired number of physical zones C if (imax .gt. II) then write(*,*) 'II less than imax = ', imax stop endif if (jmax .gt. JJ) then write(*,*) 'JJ less than jmax = ', jmax stop endif if (kmax .gt. KK) then write(*,*) 'KK less than kmax = ', kmax stop endif if (max(imax,jmax,kmax)+4 .gt. maxsweep) then write(*,*) 'maxsweep too small' stop endif C C Open history file and write out a description of the run C open (unit=8,file=hstfile,form='formatted') write (8,*) 'History File for VH-1 run' write (8,911) prefx write (8,913) mvyfile 911 format(' Filename prefix: ',a5) 913 format(' HDF movie filename: ',a10) write (8,*) 'Courant number = ',courant write (8,*) 'Grid dimensions: imax= ',imax,' jmax= ' & ,jmax,' kmax= ',kmax write (8,*) 'Outer x location = ',xmax write (8,*) 'Inner x location = ',xmin write (8,*) 'Outer y location = ',ymax write (8,*) 'Inner y location = ',ymin write (8,*) 'Outer z location = ',zmax write (8,*) 'Inner z location = ',zmin C if ((rstrt.eq.'NO').or.(rstrt.eq.'no')) then C C Initialize variables for new problem; C timep = 0. ncycle = 0 ncycp = 0 ncycd = 0 ncycm = 0 nfile = 1000 gamm = gam - 1. boltzman = 1.3806e-16 avgmass = 1.21e-24 pi = 3.141592654 G = 6.67e-8 GM = G * 3.7* 1.99e33 C call grid(imax,xmin,xmax,zxa,zdx) call grid(jmax,ymin,ymax,zya,zdy) call grid(kmax,zmin,zmax,zza,zdz) call init call dtcon else C C Restart from old dumpfile... C dmpfile = dmpfile(1:6) // rstrt(1:2) write(8,*) 'Restarting from ',dmpfile open(unit=4,file=dmpfile,status='old',form='unformatted') tmpcyc=ncycend read(4,err=998) ADUMP, BDUMP close(4) ncycend=tmpcyc C C Increment dump filename C if(dmpfile(8:8) .eq. 'z' .or. dmpfile(8:8) .eq. 'Z') then sf1 = dmpfile(7:7) isf1 = ichar(sf1) + 1 sf2 = dmpfile(8:8) isf2 = ichar(sf2) - 25 dmpfile(7:7) = char(isf1) dmpfile(8:8) = char(isf2) else sf2 = dmpfile(8:8) isf2 = ichar(sf2) + 1 dmpfile(8:8) = char(isf2) endif C endif C write (*,*) 'Starting on cycle number: ', ncycle C C############################################################################ C C MAIN COMPUTATIONAL LOOP 1 continue C write(*,*) 'STEP =', ncycle C if ( time + dt .gt. endtime ) then C go exactly to the end dt = endtime - time ncycend = ncycle ncycp = nprin ncycd = ndump endif C C Evolution step: PPM - 3rd order piecewize parabolic. C C Alternate sweeps and radiate to approximate 2nd order operator splitting C if ((ncycle - 2*(ncycle/2)) .eq. 0) then call sweepx if(ndim.gt.1) call sweepy if(ndim.eq.3) call sweepz C call radiate time = time + dt timep = timep + dt else C call radiate if(ndim.eq.3) call sweepz if(ndim.gt.1) call sweepy call sweepx time = time + dt timep = timep + dt call dtcon endif C C Loop Bookkeeping ncycle = ncycle + 1 ncycp = ncycp + 1 ncycd = ncycd + 1 ncycm = ncycm + 1 C call routine to print movie images if ( ncycm.ge.nmovie ) then call images(mvyfile,1) ncycm = 0 endif C call routine to print out numbers if ( ncycp.ge.nprin .or. timep.ge.tprin ) then call prin(prefx) timep = 0. ncycp = 0 endif C dump out a complete dataset if ( ncycd .ge. ndump ) then call dump(dmpfile) ncycd = 0. endif C if ( ncycle.lt.ncycend ) goto 1 C C END OF MAIN LOOP C######################################################################### C close( 8 ) close( 3 ) C stop 998 write(*,*) 'Error reading restart dump' stop 999 write(*,*) 'Error reading namelist file' stop end