c some revision in 18 March 2025
c but mostly from ~1996
program box,13
include 'box.size'
include 'box.common'
dimension dwdt(in,jn,kn,3),dudt(in,jn,kn,3),
# dvdt(in,jn,kn,3)
common/spare/pp(in,jn,kn),others(in,jn,kn,4)
character case4*4,fname*15,hname*9,outpath*60
LOGICAL :: fileExists
equivalence (ap(95),case4)
c--- define constants, initialize variables ----
c--- the characters 'case' will be replaced by the 4 characters
c--- in the variable named case4. (Unfortunately, case is a reserved
c--- word in fortran)
outpath='../boxout/'
fname='cased000000.sav'
hname='case.hist'
open (unit=22,file='filetimes.txt',status='old')
write(*,*) 'continue from previous case?'
read(*,*) icont
c---- start fresh ----
if (icont.eq.0) then
open (unit=10,file='grid.dat',status='unknown')
read(10,*) x,y,z
write (*,*) ' enter four characters: ''case'''
read (*,*) case4
call change(1)
ap(98)=in
ap(99)=jn
ap(100)=kn
call init(dudt,dvdt,dwdt)
time=0.
read (22,*) tdisk
cput=0.
dt1=1.
dt2=1.
na=1
nb=2
nc=3
end if
c--- else restart ----
if(icont.eq.1) then
open (unit=30,file=trim(outpath)//'restart.sav',
# status='unknown',form='unformatted')
read(30) ap,x,y,z,b,u,v,w,pp,p,vnu,dudt,dvdt,dwdt
close (unit=30)
write (*,*) ' continuing with ',case4
c time is from ap(15)
write (*,*) ' time now: ',time
c--- apparently there was a problem with vnu not being updated, so dammit
do 1200 k=1,kn
dammit=vnu(k)
vnu(k)=amax1( rei,reit*2*(z(1,1,k)-.5) )
c print *,k,z(1,1,k),vnu(k),dammit
1200 continue
39 read (22,*) tdisk
if (tdisk.le.time) goto 39
write (*,*) 'will write .sav to disk at time:',tdisk
call change(0)
na=ap(88)
nb=ap(87)
nc=ap(86)
endif
c--- for either ---
write(hname(1:4),'(a4)') case4
INQUIRE(FILE=trim(outpath)//hname, EXIST=fileExists)
open (unit=55,file=trim(outpath)//hname,
# status='unknown',form='formatted')
# if (icont.eq.0) write (55,'(es15.5)') ap
if (.not.fileExists) write (55,'(es15.5)') ap
close (unit=55)
call makebnds
call maketrns(dtdif)
kl=kn
do while (z(1,1,kl).gt.0.5)
kl=kl-1
end do
c write (*,*) 'dtdif=',dtdif
write(*,*) 'time now =',time,' enter tstop,tdiag'
read(*,*) tstop,tdiag
timel=time+1.e-5
c npict=-1
c299 read (20,*) tpict
c npict=npict+1
c if (tpict.lt.time) goto 299
c write (*,*) ' tpict=',tpict
nloc=0
elapsed=timef()
write (*,*)
# 'steps time dt ',
# ' s w p'
c--- start main loop ----
300 continue
c---- calculate advection and buoyancy terms ----
if (order.ne.3.) then
call advect(dudt(1,1,1,na),dvdt(1,1,1,na),
# dwdt(1,1,1,na),vcmax)
else
call advect3(dudt(1,1,1,na),dvdt(1,1,1,na),
# dwdt(1,1,1,na),vcmax)
end if
call diffus(dudt(1,1,1,na),dvdt(1,1,1,na),
# dwdt(1,1,1,na),vnu)
c---- set time step ----
dt=amin1(alim/(1.e-10+vcmax),diff*dtdif,dtmax)
c---- solve diagnostic equation for pressure, calc pgf terms ---
call incmp(dudt(1,1,1,na),dvdt(1,1,1,na),dwdt(1,1,1,na))
c---- recompute coefficients for A-B scheme ----
call abcoef(dt1/dt,(dt1+dt2)/dt)
c---step forward---
ko=amin1(steps+1.,3.)
call march(ko,na,nb,nc,dudt,dvdt,dwdt)
timel=time
time=time+dt
nloc=nloc+1
steps=steps+1
dt2=dt1
dt1=dt
ndum=nc
nc=nb
nb=na
na=ndum
c---write diagnostics
if (amod(time,tdiag).lt.amod(timel,tdiag).or.
# time.gt.tstop.or.nloc.eq.1) then
call speed(pp,u,v,w)
write (*,90)steps,time,dt,rmll(pp,kl),rmll(w,kl),rmpll(p,kl)
open (unit=55,file=trim(outpath)//hname,position='append',
# status='unknown',form='formatted')
c write (55,90)steps,time,dt,rm(b),rm(u),rm(v),rm(w),rmp(p)
write (55,90)steps,time,dt,rmll(pp,kl),rmll(w,kl),rmpll(p,kl)
close (unit=55)
90 format(f9.0,f9.3,f9.6,5(f10.5))
endif
if (time.ge.tdisk.or.time.ge.tstop) then
!itime=1000*time
itime=1000.*tdisk+.5
if (time.ge.tdisk) read (22,*) tdisk
elapsed=timef()
cput=cput+elapsed
ap(88)=na
ap(87)=nb
ap(86)=nc
call pstuff(p,pp)
write(fname(1:4),'(a4)') case4
if (in.eq.181) write(fname(5:5),'(a1)') 'D'
write(fname(6:11),'(i6.6)') itime
open (unit=69,file=trim(outpath)//fname,STATUS= 'unknown',
# form='unformatted')
write (69) ap,x,y,z,b,u,v,w,pp,p,vnu,dudt,dvdt,dwdt
close (unit=69)
write (*,*) 'write file: ',fname
end if
if (time.lt.tstop) goto 300
c --- end main loop ---
write (*,*) ' FINISHED!'
write (*,*) 'time=',time,' steps=',steps
c write (*,*) 'nloc=',nloc,' elapsed/nloc', elapsed/nloc,
c # ' millisecs'
c write (*,*) 'time=',time,' cput=',cput,' steps=',steps
end
subroutine init(dudt,dvdt,dwdt) 1
include 'box.size'
include 'box.common'
dimension dwdt(in,jn,kn,3),dudt(in,jn,kn,3),
# dvdt(in,jn,kn,3)
cput=0.
time=0.
steps=0.
do 120 k=1,kn
! vnu(k)=rei+(reit-rei)*.5*(1-cos(3.1415*z(1,1,k)))
vnu(k)=amax1( rei,reit*2*(z(1,1,k)-.5) )
do 120 j=1,jn
do 120 i=1,in
u(i,j,k)=0.
v(i,j,k)=0.
w(i,j,k)=0.
rad2=x(i,j,k)**2+y(i,j,k)**2
b(i,j,k)=1.2626*exp(-10.*rad2)*
# exp(-20*(z(i,j,k)-.5)**2)
120 continue
do 121 n=1,3
do 121 k=1,kn
do 121 j=1,jn
do 121 i=1,in
dudt(i,j,k,n)=0.
dvdt(i,j,k,n)=0.
121 dwdt(i,j,k,n)=0.
iseed=1
do 122 k=1,kn
do 122 j=1,jn
do 122 i=1,in
iseed=mod(iseed*7141+54773,259200)
122 b(i,j,k)=(1.+(iseed/259199. -0.5)*.001)*b(i,j,k)
return
end
subroutine speed(s,u,v,w) 1
include 'box.size'
dimension u(in,jn,kn),v(in,jn,kn),w(in,jn,kn),s(in,jn,kn)
do i=1,in
do j=1,jn
do k=1,kn
s(i,j,k)=sqrt(u(i,j,k)**2+v(i,j,k)**2+w(i,j,k)**2)
end do
end do
end do
return
end
c--- following subroutine is not used
subroutine slice(npict,a,q,x,y,case4,time)
include 'box.size'
dimension q(in,jn),x(in,jn),y(in,jn)
character fname*16,case4*4,a*1
fname='sub/casew000.dat'
write(fname(5:8),'(a4)') case4
write(fname(10:12),'(i3.3)') npict
write(fname(9:9),'(a1)') a
open (unit=69,file=fname,STATUS= 'unknown')
write (69,301) q,x,y,time
301 format(4(e15.7))
close (unit=69)
write (*,*) ' write movie file= ',fname,' at',time
return
end
subroutine change(iread) 2
include 'box.size'
include 'box.common'
if (iread.gt.0) open (unit=62,file='box.input',status='unknown')
if (iread.eq.1) read (62,f)
321 write (*,*) ' Current parameters are:'
write (*,f)
write (*,*)' Do you want to change any parameters?'
write (*,*) ' enter yes, and store=2, yes, no store=1, no=0'
read (*,*) ichn
if (ichn.ne.0) then
write (*,*)
# ' enter changes with: (space)$f example=1.0, etc=2. $end'
read(*,f)
rewind 62
if (ichn.eq.2) write(62,f)
goto 321
end if
close (unit=62)
return
end