program gridmaker,6
c 9 June 2014. Unaltered version of a program from 1990's previously named gr2.f.
c The progam is simplified version of a code once used for continuous dynamic grid
c adaptation. Here those CDGA routines are used to make a single static grid, placing
c more grid points where a weight function is large.
c Compile and run to make a text file grid.dat, which is read in by boxn.f.
c The file specifies the stretched grid, which places higher resolution in
c the central and lower part of the domain, as is useful for tornado simulations.
include 'gridmaker.size'
parameter (ih=in/2+1,jh=jn/2+1)
common/par/rei,hite,wide,alim,diff,pit,rlx,str,chop,filtn,adjn,
# abc,dtmax,cin,xtra2,xtra3
# ,dx,dtn,time
dimension x(in,jn,kn),y(in,jn,kn),z(in,jn,kn),t(in,jn),
# xs(is,js),ys(is,js),wf(is,js),ug(in,jn),vg(in,jn),
# w(in,jn,kn)
dimension xa(ih,kn),za(ih,kn)
wide=4.
hite=4.
chop=1.0
filtn=0.
str=16.
ampc=7.
dt=0.1 ! bogus time step, any dt works here
xc=wide/2.
yc=hite/2.
stz=.2
c call makebnds
call gridprep
do k=1,kn,kn-1
write (*,*) ' k=',k
rad1=1.0 !original
rad2=1.0 !original
rad3=.3 !original
rad3=.5 !not bad
rad3=.25
rad1=.5
rad2=.5
call gridinit(x(1,1,k),y(1,1,k),wide,hite,in,jn)
call gridinit(xs,ys,wide,hite,is,js)
if (k.eq.1) then
str=64. !original
ampc=15. !original
str=128. !not bad
ampc=15. !not bad
str=128.
ampc=30.
else
str=8.
ampc=4.
end if
do n=1,16
do j=1,jn
do i=1,in
t(i,j)=0.
r=sqrt((x(i,j,k)-xc)**2+0*(y(i,j,k)-yc)**2)
if (r.lt.rad1) t(i,j)=amax1(cos(1.571*r/rad1),.5)
r=sqrt(0*(x(i,j,k)-xc)**2+(y(i,j,k)-yc)**2)
if (r.lt.rad2) t(i,j)=t(i,j)+
# amax1(cos(1.571*r/rad2),.5)
r=sqrt((x(i,j,k)-xc)**2+(y(i,j,k)-yc)**2)
if (r.lt.rad3) t(i,j)=t(i,j)+
# ampc*amax1(cos(1.571*r/rad3),.5)
w(i,j,k)=t(i,j)
end do
end do
call newwf(t,x(1,1,k),y(1,1,k),wf,xs,ys,filtn)
call grid(xs,ys,wf,str,chop,4*js)
call itrpgr(x(1,1,k),y(1,1,k),xs,ys,ug,vg,dt)
do i=1,in
do j=1,jn
x(i,j,k)=x(i,j,k)+dt*ug(i,j)
y(i,j,k)=y(i,j,k)+dt*vg(i,j)
enddo
enddo
write(*,*) x(im/4+1,jm/4+1,k)
end do
end do
do k=1,kn
gz=(k-1.)/(kn-1.)
zz=(stz*gz+(2.-stz)*gz**3)/(1.+gz**2)
do i=1,in
do j=1,jn
z(i,j,k)=zz
stq=.4
zx=(stq*zz+(2.-stq)*zz**3)/(1.+zz**2)
zx=0.+0.5*(1-cos(3.1415*zz))
x(i,j,k)=x(i,j,1)+zx*(x(i,j,kn)-x(i,j,1))
y(i,j,k)=y(i,j,1)+zx*(y(i,j,kn)-y(i,j,1))
end do
end do
end do
do k=1,kn
do i=1,in
do j=1,jn
x(i,j,k)=x(i,j,k)-.5*wide
y(i,j,k)=y(i,j,k)-.5*wide
end do
end do
end do
do k=1,kn
do i=1,ih
xa(i,k)=x(ih-1+i,jh,k)
za(i,k)=z(ih-1+i,jh,k)
end do
end do
open (unit=30,file='grid.dat',
# status='unknown')
write (30,'(f10.6)') x,y,z
open (unit=31,file='grid2d.dat',
# status='unknown')
write (31,'(i4,i4)')ih,kn
write (31,'(f10.6)') xa,za
end
c******************************************************************************
c**** grid adaption subroutines ****
c*
subroutine filtw(wf,dm,is,js) 1
dimension wf(is,js),dm(is,js)
do 7 i=1,is
il=i-1
ir=i+1
qr=1.
if (i.eq.1) il=ir
if (i.eq.is) then
ir=il
qr=0.00000
end if
do 7 j=1,js
jb=j-1
jt=j+1
qt=1.
qb=1.
if (j.eq.1)then
jb=jt
qb=0.000
end if
if (j.eq.js) then
jt=jb
qt=0.0000
end if
7 dm(i,j)=0.5*wf(i,j)+
# .125*(wf(il,j)+qr*wf(ir,j)+qt*wf(i,jt)+qb*wf(i,jb))
c wfmax=0.
do 5 i=1,is
il=i-1
ir=i+1
qr=1.
if (i.eq.1) il=ir
if (i.eq.is) then
ir=il
qr=0.
end if
do 5 j=1,js
jb=j-1
jt=j+1
qt=1.
qb=1.
if (j.eq.1)then
jb=jt
qb=0.
end if
if (j.eq.js) then
jt=jb
qt=0.
end if
5 wf(i,j)=0.5*dm(i,j)+
# .125*(dm(il,j)+qr*dm(ir,j)+qt*dm(i,jt)+qb*dm(i,jb))
return
end
subroutine newwf(v,xg,yg,wf,xs,ys,filtn) 1,1
include 'gridmaker.size'
dimension xg(in,jn),yg(in,jn),xs(is,js),ys(is,js)
dimension wf(is,js),v(in,jn),dum(ng)
common/trng/ir(is),il(is),tr(is),tl(is),
# jt(js),jb(js),tt(js),tb(js),
# ii(ng),jj(ng),iis(ng),jjs(ng)
do 5 n=1,ng
5 dum(n)=0.
do 20 iii=-li,li,2
do 20 jjj=-lj,lj,2
do 10 n=1,(is-2)*(js-2)
ih=(iis(n)-1)*li+1+iii
jh=(jjs(n)-1)*lj+1+jjj
10 dum(n)=dum(n)+v(ih,jh)
20 continue
do 25 n=1,(is-2)*(js-2)
i=iis(n)
j=jjs(n)
25 wf(i,j)=dum(n)
do 30 j=2,js-1
wf(1,j)=wf(2,j)
30 wf(is,j)=wf(is-1,j)
do 40 i=1,is
wf(i,1)=wf(i,2)
40 wf(i,js)=wf(i,js-1)
do 50 a=1.,filtn
50 call filtw(wf,dum,is,js)
return
end
c******************
subroutine itrpgr(xg,yg,xs,ys,dxdt,dydt,dt) 1,2
include 'gridmaker.size'
dimension xg(in,jn),yg(in,jn),xs(is,js),ys(is,js)
dimension dxdt(in,jn),dydt(in,jn)
call bicub(xs,xg,dxdt,dt,1)
call bicub(ys,yg,dydt,dt,2)
return
end
subroutine bicub(f,fg,dfdt,dt,it) 2,1
include 'gridmaker.size'
dimension fg(in,jn),dfdt(in,jn),f(is,js)
dimension f1(is,js),f2(is,js),f12(is,js)
dimension q(4),q1(4),q2(4),q12(4),c(4,4)
common/cub/igv(lsm),jgv(lsm),xxv(lsm),yyv(lsm)
rdt=1./dt
rdx=.5/li
rdy=.5/lj
do 2 i=1,is
do 2 j=1,js
f1(i,j)=0.
f2(i,j)=0.
2 f12(i,j)=0.
if (it.eq.1) then
do 21 j=1,js
f1(1,j)=2.*rdx*(f(2,j)-f(1,j))
21 f1(is,j)=2.*rdx*(f(is,j)-f(is-1,j))
do 23 i=2,is-1
il=i-1
ir=i+1
f1(i,1)=rdx*(f(ir,1)-f(il,1))
23 f1(i,js)=rdx*(f(ir,js)-f(il,js))
end if
if (it.eq.2) then
do 31 i=1,is
f2(i,1)=2.*rdy*(f(i,2)-f(i,1))
31 f2(i,js)=2.*rdy*(f(i,js)-f(i,js-1))
do 33 j=2,js-1
jb=j-1
jt=j+1
f2(1,j)=rdy*(f(1,jt)-f(1,jb))
33 f2(is,j)=rdy*(f(is,jt)-f(is,jb))
end if
do 6 i=2,is-1
il=i-1
ir=i+1
do 6 j=2,js-1
jb=j-1
jt=j+1
f1(i,j)=rdx*(f(ir,j)-f(il,j))
f2(i,j)=rdy*(f(i,jt)-f(i,jb))
6 f12(i,j)=rdx*rdy*(f(ir,jt)+f(il,jb)-f(il,jt)-f(ir,jb))
do 10 i=1,is-1
do 10 j=1,js-1
q(1)=f(i,j)
q(4)=f(i,j+1)
q(3)=f(i+1,j+1)
q(2)=f(i+1,j)
q1(1)=f1(i,j)
q1(4)=f1(i,j+1)
q1(3)=f1(i+1,j+1)
q1(2)=f1(i+1,j)
q2(1)=f2(i,j)
q2(4)=f2(i,j+1)
q2(3)=f2(i+1,j+1)
q2(2)=f2(i+1,j)
q12(1)=f12(i,j)
q12(4)=f12(i,j+1)
q12(3)=f12(i+1,j+1)
q12(2)=f12(i+1,j)
call bcucof(q,q1,q2,q12,1.*li,1.*lj,c)
do 55 nv=1,lsm
ig=(i-1)*li+1+igv(nv)
jg=(j-1)*lj+1+jgv(nv)
xx=xxv(nv)
yy=yyv(nv)
val=
# (((((c(4,4)*yy+c(3,4))*yy+c(2,4))*yy+c(1,4))*xx
# +(((c(4,3)*yy+c(3,3))*yy+c(2,3))*yy+c(1,3)))*xx
# +(((c(4,2)*yy+c(3,2))*yy+c(2,2))*yy+c(1,2)))*xx
# +(((c(4,1)*yy+c(3,1))*yy+c(2,1))*yy+c(1,1))
c val=0.
c do 65 lc=4,1,-1
c65 val=xx*val+
c # ((c(4,lc)*yy+c(3,lc))*yy+c(2,lc))*yy+c(1,lc)
55 dfdt(ig,jg)=(val-fg(ig,jg))*rdt
10 continue
return
end
SUBROUTINE GRIDINIT(X,Y,XMAX,YMAX,In,Jn) 2
DIMENSION X(In,Jn),Y(In,Jn)
DO 10 I=1,In
DO 10 J=1,Jn
X(I,J)=XMAX*(I-1)/(In-1)
10 Y(I,J)=YMAX*(J-1)/(Jn-1)
RETURN
END
function clip(x,c)
clip=alog(1+x/c)*c
return
end
subroutine oldgrid(x,y,wfunc,str,chop,ncrunch)
c** uses eq. 12, page 199 of Thompson et al. book
c** "Numerical Grid Generation". This was used for all research before
c** November 1995
include 'gridmaker.size'
dimension x(is,js),y(is,js),z(is,js),wfunc(is,js),
# w(is,js),dxdt(is,js),dydt(is,js),dzdt(is,js),rate(is,js),
# dxdn(ng),dydn(ng),ratn(ng)
common/trng/ira(is),ila(is),tr(is),tl(is),
# jta(js),jba(js),tt(js),tb(js),
# ii(ng),jj(ng),iis(ng),jjs(ng)
equivalence (dxdn(1),dxdt(1,1))
equivalence (dydn(1),dydt(1,1))
equivalence (ratn(1),rate(1,1))
dummy=0.
do 100 i=1,is
do 100 j=1,js
100 dummy=amax1(wfunc(i,j),dummy)
dummy=chop*dummy
factor=0.
if (dummy.ne.0.) factor=str/dummy
do 101 i=1,is
do 101 j=1,js
101 w(i,j)=1.+factor*amin1(dummy,wfunc(i,j))
cx factor=1/dummy
cx do 91 i=1,is
cx do 91 j=1,js
cx91 w(i,j)=factor*wfunc(i,j)
cx factor=clip(1.,chop)
cx factor=clip(factor,chop)
cx factor=str/factor
cx do 92 i=1,is
cx do 92 j=1,js
cx92 w(i,j)=clip(w(i,j),chop)
cx do 93 i=1,is
cx do 93 j=1,js
cx93 w(i,j)=1.+factor*clip(w(i,j),chop)
do 1001 nn=1,ncrunch
do 10 n=1,ng
i=ii(n)
j=jj(n)
ir=ira(i)
il=ila(i)
jt=jta(j)
jb=jba(j)
x1=0.5*(x(ir,j)-x(i,j))*tr(i)
# +0.5*(x(i,j)-x(il,j))*tl(i)
y1=0.5*(y(ir,j)-y(il,j))
x2=0.5*(x(i,jt)-x(i,jb))
y2=0.5*(y(i,jt)-y(i,j))*tt(j)
# +0.5*(y(i,j)-y(i,jb))*tb(j)
a11=y2
a12=-x2
a21=-y1
a22=x1
g11=a11*a11+a12*a12
g22=a21*a21+a22*a22
g12=a11*a21+a12*a22
p=0.5*(w(ir,j)-w(il,j))/w(i,j)
q=0.5*(w(i,jt)-w(i,jb))/w(i,j)
dxdn(n)=g11*(x(ir,j)-2.*x(i,j)+x(il,j)+
# 0.5*p*(x(ir,j)-x(il,j)))
# +g22*(x(i,jt)-2.*x(i,j)+x(i,jb)+
# 0.5*q*(x(i,jt)-x(i,jb)))
# +.5*g12*(x(ir,jt)-x(ir,jb)-
# x(il,jt)+x(il,jb))
dydn(n)=g11*(y(ir,j)-2.*y(i,j)+y(il,j)+
# 0.5*p*(y(ir,j)-y(il,j)))
# +g22*(y(i,jt)-2.*y(i,j)+y(i,jb)+
# 0.5*q*(y(i,jt)-y(i,jb)))
# +.5*g12*(y(ir,jt)-y(ir,jb)-
# y(il,jt)+y(il,jb))
ratn(n)=0.5/(g22+g11)
10 continue
do 20 i=1,is
do 20 j=2,js-1
20 y(i,j)=y(i,j)+dydt(i,j)*rate(i,j)
do 21 i=2,is-1
do 21 j=1,js
21 x(i,j)=x(i,j)+dxdt(i,j)*rate(i,j)
1001 continue
return
end
subroutine gridprep 1
include 'gridmaker.size'
common/cub/igv(lsm),jgv(lsm),xxv(lsm),yyv(lsm)
common/trng/ir(is),il(is),tr(is),tl(is),
# jt(js),jb(js),tt(js),tb(js),
# ii(ng),jj(ng),iis(ng),jjs(ng)
n=0
do 55 i=0,li
do 55 j=0,lj
n=n+1
igv(n)=i
xxv(n)=1.*i/li
jgv(n)=j
55 yyv(n)=1.*j/lj
n=0
do 10 j=1,js
do 10 i=1,is
n=n+1
ii(n)=i
10 jj(n)=j
n=0
do 15 i=2,is-1
do 15 j=2,js-1
n=n+1
iis(n)=i
jjs(n)=j
15 continue
do 20 i=1,is
ir(i)=i+1
il(i)=i-1
tr(i)=1.
tl(i)=1.
if (i.eq.is) then
ir(i)=i-1
tr(i)=0.
tl(i)=2.
end if
if (i.eq.1) then
il(i)=i+1
tr(i)=2.
tl(i)=0.
end if
20 continue
do 30 j=1,js
jt(j)=j+1
jb(j)=j-1
tt(j)=1.
tb(j)=1.
if (j.eq.js) then
jt(j)=j-1
tt(j)=0.
tb(j)=2.
end if
if (j.eq.1) then
jb(j)=j+1
tt(j)=2.
tb(j)=0.
end if
30 continue
return
end
SUBROUTINE BCUCOF(Y,Y1,Y2,Y12,D1,D2,CL) 1
DIMENSION Y(4),Y1(4),Y2(4),Y12(4),CL(16),X(16),WT(16,16)
DATA WT/1.,0.,-3.,2.,4*0.,-3.,0.,9.,-6.,2.,0.,-6.,
* 4.,8*0.,3.,0.,-9.,6.,-2.,0.,6.,-4.,10*0.,9.,-6.,
* 2*0.,-6.,4.,2*0.,3.,-2.,6*0.,-9.,6.,2*0.,6.,-4.,
* 4*0.,1.,0.,-3.,2.,-2.,0.,6.,-4.,1.,0.,-3.,2.,8*0.,
* -1.,0.,3.,-2.,1.,0.,-3.,2.,10*0.,-3.,2.,2*0.,3.,
* -2.,6*0.,3.,-2.,2*0.,-6.,4.,2*0.,3.,-2.,0.,1.,-2.,
* 1.,5*0.,-3.,6.,-3.,0.,2.,-4.,2.,9*0.,3.,-6.,3.,0.,
* -2.,4.,-2.,10*0.,-3.,3.,2*0.,2.,-2.,2*0.,-1.,1.,
* 6*0.,3.,-3.,2*0.,-2.,2.,5*0.,1.,-2.,1.,0.,-2.,4.,
* -2.,0.,1.,-2.,1.,9*0.,-1.,2.,-1.,0.,1.,-2.,1.,10*0.,
* 1.,-1.,2*0.,-1.,1.,6*0.,-1.,1.,2*0.,2.,-2.,2*0.,-1.,1./
D1D2=D1*D2
DO 11 I=1,4
X(I)=Y(I)
X(I+4)=Y1(I)*D1
X(I+8)=Y2(I)*D2
X(I+12)=Y12(I)*D1D2
11 CONTINUE
DO 13 I=1,16
CL(I)=
# WT(I,1)*X(1)
# +WT(I,2)*X(2)
# +WT(I,3)*X(3)
# +WT(I,4)*X(4)
# +WT(I,5)*X(5)
# +WT(I,6)*X(6)
# +WT(I,7)*X(7)
# +WT(I,8)*X(8)
# +WT(I,9)*X(9)
# +WT(I,10)*X(10)
# +WT(I,11)*X(11)
# +WT(I,12)*X(12)
# +WT(I,13)*X(13)
# +WT(I,14)*X(14)
# +WT(I,15)*X(15)
# +WT(I,16)*X(16)
13 CONTINUE
RETURN
END
subroutine grid(x,y,wfunc,str,chop,ncrunch) 1
c***uses the grid equation derived from Montreal, November 1995
include 'gridmaker.size'
dimension x(is,js),y(is,js),z(is,js),wfunc(is,js),
# w(is,js),dxdt(is,js),dydt(is,js),dzdt(is,js),rate(is,js),
# dxdn(ng),dydn(ng),ratn(ng)
common/trng/ira(is),ila(is),tr(is),tl(is),
# jta(js),jba(js),tt(js),tb(js),
# ii(ng),jj(ng),iis(ng),jjs(ng)
equivalence (dxdn(1),dxdt(1,1))
equivalence (dydn(1),dydt(1,1))
equivalence (ratn(1),rate(1,1))
dummy=0.
do 100 i=1,is
do 100 j=1,js
100 dummy=amax1(wfunc(i,j),dummy)
dummy=chop*dummy
factor=0.
if (dummy.ne.0.) factor=str/dummy
do 101 i=1,is
do 101 j=1,js
101 w(i,j)=1.+factor*amin1(dummy,wfunc(i,j))
do 1001 nn=1,ncrunch
do 10 n=1,ng
i=ii(n)
j=jj(n)
ir=ira(i)
il=ila(i)
jt=jta(j)
jb=jba(j)
x1=0.5*(x(ir,j)-x(i,j))*tr(i)
# +0.5*(x(i,j)-x(il,j))*tl(i)
y1=0.5*(y(ir,j)-y(il,j))
x2=0.5*(x(i,jt)-x(i,jb))
y2=0.5*(y(i,jt)-y(i,j))*tt(j)
# +0.5*(y(i,j)-y(i,jb))*tb(j)
rg=1./(x1*y2-x2*y1)
a11=y2*rg
a12=-x2*rg
a21=-y1*rg
a22=x1*rg
g11=a11*a11+a12*a12
g22=a21*a21+a22*a22
dxdn(n)=a11*(x(ir,j)-2.*x(i,j)+x(il,j))+
# a22*(x(i,jt)-2.*x(i,j)+x(i,jb))+
# (a12+a21)*.25*(x(ir,jt)-x(ir,jb)-
# x(il,jt)+x(il,jb))
# +0.5*(w(ir,j)-w(il,j))/w(i,j)
dydn(n)=a11*(y(ir,j)-2.*y(i,j)+y(il,j))+
# a22*(y(i,jt)-2.*y(i,j)+y(i,jb))+
# (a12+a21)*.25*(y(ir,jt)-y(ir,jb)-
# y(il,jt)+y(il,jb))
# +0.5*(w(i,jt)-w(i,jb))/w(i,j)
c ratn(n)=0.3/sqrt(g22+g11) !original
ratn(n)=0.22/sqrt(g22+g11)
10 continue
do 20 i=1,is
do 20 j=2,js-1
20 y(i,j)=y(i,j)+dydt(i,j)*rate(i,j)
do 21 i=2,is-1
do 21 j=1,js
21 x(i,j)=x(i,j)+dxdt(i,j)*rate(i,j)
1001 continue
return
end