c some minor revision in 18 March 2025
c but mostly from ~1996
c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
subroutine abcoef(p,q) 1
common/adams/ab(3,3)
ab(3,1)=0.
ab(2,1)=0.
ab(1,1)=1.
ab(3,2)=0.
ab(2,2)=-0.5/p
ab(1,2)=1-ab(2,2)
ab(3,3)=(2+3*p)/(6*q*(q-p))
ab(2,3)=-(1+2*q*ab(3,3))/(2*p)
ab(1,3)=1.-ab(2,3)-ab(3,3)
return
end
subroutine march(ko,na,nb,nc,dudt,dvdt,dwdt) 1
include 'box.size'
include 'box.common'
dimension dudt(in,jn,kn,3),
& dvdt(in,jn,kn,3),dwdt(in,jn,kn,3)
do k=slip,kn
do j=1,jn
do i=2,in-1
u(i,j,k)=u(i,j,k)+dt*(ab(1,ko)*dudt(i,j,k,na)
# +ab(2,ko)*dudt(i,j,k,nb)
# +ab(3,ko)*dudt(i,j,k,nc))
enddo
u(1,j,k)=0.
u(in,j,k)=0.
enddo
enddo
do k=slip,kn
do i=1,in
do j=2,jn-1
v(i,j,k)=v(i,j,k)+dt*(ab(1,ko)*dvdt(i,j,k,na)
# +ab(2,ko)*dvdt(i,j,k,nb)
# +ab(3,ko)*dvdt(i,j,k,nc))
enddo
v(i,1,k)=0.
v(i,jn,k)=0.
enddo
enddo
do k=2,kn-1
do j=1,jn
do i=1,in
w(i,j,k)=w(i,j,k)+dt*(ab(1,ko)*dwdt(i,j,k,na)
# +ab(2,ko)*dwdt(i,j,k,nb)
# +ab(3,ko)*dwdt(i,j,k,nc))
w(i,j,1)=0.
w(i,j,kn)=0.
enddo
enddo
enddo
return
end
subroutine pstuff(q,qq) 1
c* compute average p on u-grid; compute density (actually, 1/density) at
c* u-grid nodes
include 'box.size'
include 'box.common'
dimension q(im,jm,km),qq(in,jn,kn)
do k=1,kn
do j=1,jn
do i=1,in
qq(i,j,k)=0.125*(q(ilu(i),jbu(j),kuu(k))+q(ilu(i),jtu(j),kuu(k))
& +q(iru(i),jtu(j),kuu(k))+q(iru(i),jbu(j),kuu(k))
& +q(ilu(i),jbu(j),kdu(k))+q(ilu(i),jtu(j),kdu(k))
& +q(iru(i),jtu(j),kdu(k))+q(iru(i),jbu(j),kdu(k)))
enddo
enddo
enddo
return
end
subroutine maketrns(dtdif) 1,1
include 'box.size'
include 'box.common'
common/spare/dummy(in,jn,kn)
c---- transformation relations on u-grid----
do k=1,kn
do j=1,jn
do i=1,in
c---- compute jacobian for u-grid ----
x1=0.25*(x(ir(i),j,k)-x(i,j,k))*tr(i)
& +0.25*(x(i,j,k)-x(il(i),j,k))*tl(i)
y1=0.25*(y(ir(i),j,k)-y(il(i),j,k))
z1=0.25*(z(ir(i),j,k)-z(il(i),j,k))
x2=0.25*(x(i,jt(j),k)-x(i,jb(j),k))
y2=0.25*(y(i,jt(j),k)-y(i,j,k))*tt(j)
& +0.25*(y(i,j,k)-y(i,jb(j),k))*tb(j)
z2=0.25*(z(i,jt(j),k)-z(i,jb(j),k))
x3=0.25*(x(i,j,ku(k))-x(i,j,kd(k)))
y3=0.25*(y(i,j,ku(k))-y(i,j,kd(k)))
z3=0.25*(z(i,j,ku(k))-z(i,j,k))*tu(k)
& +0.25*(z(i,j,k)-z(i,j,kd(k)))*td(k)
g=x1*(y2*z3-y3*z2)+y1*(x3*z2-x2*z3)+
& z1*(x2*y3-x3*y2)
dummy(i,j,k)=g**.666/vnu(k)
rg=1./g
c---- contravariant base vectors ----
a11(i,j,k)=rg*(y2*z3-y3*z2)
a12(i,j,k)=rg*(x3*z2-x2*z3)
a13(i,j,k)=rg*(x2*y3-x3*y2)
a21(i,j,k)=rg*(y3*z1-y1*z3)
a22(i,j,k)=rg*(x1*z3-x3*z1)
a23(i,j,k)=rg*(x3*y1-x1*y3)
a31(i,j,k)=rg*(y1*z2-y2*z1)
a32(i,j,k)=rg*(z1*x2-z2*x1)
a33(i,j,k)=rg*(x1*y2-x2*y1)
enddo
enddo
enddo
c---- basis vectors on p-grid----
do k=1,km
do j=1,jm
do i=1,im
x1=0.125*(
& (x(irc(i),jtc(j),kuc(k))-x(ilc(i),jtc(j),kuc(k)))
& +(x(irc(i),jbc(j),kuc(k))-x(ilc(i),jbc(j),kuc(k)))
& +(x(irc(i),jtc(j),kdc(k))-x(ilc(i),jtc(j),kdc(k)))
& +(x(irc(i),jbc(j),kdc(k))-x(ilc(i),jbc(j),kdc(k))))
y1=0.125*(
& (y(irc(i),jtc(j),kuc(k))-y(ilc(i),jtc(j),kuc(k)))
& +(y(irc(i),jbc(j),kuc(k))-y(ilc(i),jbc(j),kuc(k)))
& +(y(irc(i),jtc(j),kdc(k))-y(ilc(i),jtc(j),kdc(k)))
& +(y(irc(i),jbc(j),kdc(k))-y(ilc(i),jbc(j),kdc(k))))
z1=0.125*(
& (z(irc(i),jtc(j),kuc(k))-z(ilc(i),jtc(j),kuc(k)))
& +(z(irc(i),jbc(j),kuc(k))-z(ilc(i),jbc(j),kuc(k)))
& +(z(irc(i),jtc(j),kdc(k))-z(ilc(i),jtc(j),kdc(k)))
& +(z(irc(i),jbc(j),kdc(k))-z(ilc(i),jbc(j),kdc(k))))
x2=0.125*(
& (x(irc(i),jtc(j),kuc(k))-x(irc(i),jbc(j),kuc(k)))
& +(x(ilc(i),jtc(j),kuc(k))-x(ilc(i),jbc(j),kuc(k)))
& +(x(irc(i),jtc(j),kdc(k))-x(irc(i),jbc(j),kdc(k)))
& +(x(ilc(i),jtc(j),kdc(k))-x(ilc(i),jbc(j),kdc(k))))
y2=0.125*(
& (y(irc(i),jtc(j),kuc(k))-y(irc(i),jbc(j),kuc(k)))
& +(y(ilc(i),jtc(j),kuc(k))-y(ilc(i),jbc(j),kuc(k)))
& +(y(irc(i),jtc(j),kdc(k))-y(irc(i),jbc(j),kdc(k)))
& +(y(ilc(i),jtc(j),kdc(k))-y(ilc(i),jbc(j),kdc(k))))
z2=0.125*(
& (z(irc(i),jtc(j),kuc(k))-z(irc(i),jbc(j),kuc(k)))
& +(z(ilc(i),jtc(j),kuc(k))-z(ilc(i),jbc(j),kuc(k)))
& +(z(irc(i),jtc(j),kdc(k))-z(irc(i),jbc(j),kdc(k)))
& +(z(ilc(i),jtc(j),kdc(k))-z(ilc(i),jbc(j),kdc(k))))
x3=0.125*(
& (x(irc(i),jtc(j),kuc(k))-x(irc(i),jtc(j),kdc(k)))
& +(x(ilc(i),jtc(j),kuc(k))-x(ilc(i),jtc(j),kdc(k)))
& +(x(irc(i),jbc(j),kuc(k))-x(irc(i),jbc(j),kdc(k)))
& +(x(ilc(i),jbc(j),kuc(k))-x(ilc(i),jbc(j),kdc(k))))
y3=0.125*(
& (y(irc(i),jtc(j),kuc(k))-y(irc(i),jtc(j),kdc(k)))
& +(y(ilc(i),jtc(j),kuc(k))-y(ilc(i),jtc(j),kdc(k)))
& +(y(irc(i),jbc(j),kuc(k))-y(irc(i),jbc(j),kdc(k)))
& +(y(ilc(i),jbc(j),kuc(k))-y(ilc(i),jbc(j),kdc(k))))
z3=0.125*(
& (z(irc(i),jtc(j),kuc(k))-z(irc(i),jtc(j),kdc(k)))
& +(z(ilc(i),jtc(j),kuc(k))-z(ilc(i),jtc(j),kdc(k)))
& +(z(irc(i),jbc(j),kuc(k))-z(irc(i),jbc(j),kdc(k)))
& +(z(ilc(i),jbc(j),kuc(k))-z(ilc(i),jbc(j),kdc(k))))
g=x1*(y2*z3-y3*z2)+y1*(x3*z2-x2*z3)+
& z1*(x2*y3-x3*y2)
rg=1./g
c---- contravariant base vectors on p-grid----
c11(i,j,k)=rg*(y2*z3-y3*z2)
c12(i,j,k)=rg*(x3*z2-x2*z3)
c13(i,j,k)=rg*(x2*y3-x3*y2)
c21(i,j,k)=rg*(y3*z1-y1*z3)
c22(i,j,k)=rg*(x1*z3-x3*z1)
c23(i,j,k)=rg*(x3*y1-x1*y3)
c31(i,j,k)=rg*(y1*z2-y2*z1)
c32(i,j,k)=rg*(z1*x2-z2*x1)
c33(i,j,k)=rg*(x1*y2-x2*y1)
c---alpha used in incmp ---
alpha(i,j,k)=rlx/(
& c11(i,j,k)*c11(i,j,k)+c12(i,j,k)*c12(i,j,k)
& +c13(i,j,k)*c13(i,j,k)+
& c21(i,j,k)*c21(i,j,k)+c22(i,j,k)*c22(i,j,k)
& +c23(i,j,k)*c23(i,j,k)+
& c31(i,j,k)*c31(i,j,k)+c32(i,j,k)*c32(i,j,k)
& +c33(i,j,k)*c33(i,j,k))
enddo
enddo
enddo
c---for diffusion terms
do i=1,in
do j=1,jn
do k=1,kn
g11(i,j,k)=a11(i,j,k)*a11(i,j,k)+a12(i,j,k)*a12(i,j,k)
& +a13(i,j,k)*a13(i,j,k)
g12(i,j,k)=a21(i,j,k)*a11(i,j,k)+a22(i,j,k)*a12(i,j,k)
& +a23(i,j,k)*a13(i,j,k)
g13(i,j,k)=a31(i,j,k)*a11(i,j,k)+a32(i,j,k)*a12(i,j,k)
& +a33(i,j,k)*a13(i,j,k)
g22(i,j,k)=a21(i,j,k)*a21(i,j,k)+a22(i,j,k)*a22(i,j,k)
& +a23(i,j,k)*a23(i,j,k)
g23(i,j,k)=a31(i,j,k)*a21(i,j,k)+a32(i,j,k)*a22(i,j,k)
& +a33(i,j,k)*a23(i,j,k)
g33(i,j,k)=a31(i,j,k)*a31(i,j,k)+a32(i,j,k)*a32(i,j,k)
& +a33(i,j,k)*a33(i,j,k)
end do
end do
end do
do i=1,in
do j=1,jn
do k=1,kn
c---- laplacian of logical coordinate, for diffusion term transformation ----
x11=0.25*(x(ild(i),j,k)-2.*x(icd(i),j,k)+x(ird(i),j,k))
y11=0.25*(y(ild(i),j,k)-2.*y(icd(i),j,k)+y(ird(i),j,k))
z11=0.25*(z(ild(i),j,k)-2.*z(icd(i),j,k)+z(ird(i),j,k))
x22=0.25*(x(i,jbd(j),k)-2.*x(i,jcd(j),k)+x(i,jtd(j),k))
y22=0.25*(y(i,jbd(j),k)-2.*y(i,jcd(j),k)+y(i,jtd(j),k))
z22=0.25*(z(i,jbd(j),k)-2.*z(i,jcd(j),k)+z(i,jtd(j),k))
x33=0.25*(x(i,j,kdd(k))-2.*x(i,j,kcd(k))+x(i,j,kud(k)))
y33=0.25*(y(i,j,kdd(k))-2.*y(i,j,kcd(k))+y(i,j,kud(k)))
z33=0.25*(z(i,j,kdd(k))-2.*z(i,j,kcd(k))+z(i,j,kud(k)))
c---- cross derivatives ----
x12=0.0625*(
& +(x(ir(i),jt(j),k)-x(ir(i),jb(j),k))
& -(x(il(i),jt(j),k)-x(il(i),jb(j),k)))
y12=0.0625*(
& +(y(ir(i),jt(j),k)-y(ir(i),jb(j),k))
& -(y(il(i),jt(j),k)-y(il(i),jb(j),k)))
z12=0.0625*(
& +(z(ir(i),jt(j),k)-z(ir(i),jb(j),k))
& -(z(il(i),jt(j),k)-z(il(i),jb(j),k)))
x13=0.0625*(
& +(x(ir(i),j,ku(k))-x(ir(i),j,kd(k)))
& -(x(il(i),j,ku(k))-x(il(i),j,kd(k))))
y13=0.0625*(
& +(y(ir(i),j,ku(k))-y(ir(i),j,kd(k)))
& -(y(il(i),j,ku(k))-y(il(i),j,kd(k))))
z13=0.0625*(
& +(z(ir(i),j,ku(k))-z(ir(i),j,kd(k)))
& -(z(il(i),j,ku(k))-z(il(i),j,kd(k))))
x23=0.0625*(
& +(x(i,jt(j),ku(k))-x(i,jb(j),ku(k)))
& -(x(i,jt(j),kd(k))-x(i,jb(j),kd(k))))
y23=0.0625*(
& +(y(i,jt(j),ku(k))-y(i,jb(j),ku(k)))
& -(y(i,jt(j),kd(k))-y(i,jb(j),kd(k))))
z23=0.0625*(
& +(z(i,jt(j),ku(k))-z(i,jb(j),ku(k)))
& -(z(i,jt(j),kd(k))-z(i,jb(j),kd(k))))
d2z1(i,j,k)=-a11(i,j,k)*
& (g11(i,j,k)*x11+2*g12(i,j,k)*x12+g22(i,j,k)*x22
& +2*g13(i,j,k)*x13+2*g23(i,j,k)*x23+g33(i,j,k)*x33)
& -a12(i,j,k)*
& (g11(i,j,k)*y11+2*g12(i,j,k)*y12+g22(i,j,k)*y22
& +2*g13(i,j,k)*y13+2*g23(i,j,k)*y23+g33(i,j,k)*y33)
& -a13(i,j,k)*
& (g11(i,j,k)*z11+2*g12(i,j,k)*z12+g22(i,j,k)*z22
& +2*g13(i,j,k)*z13+2*g23(i,j,k)*z23+g33(i,j,k)*z33)
d2z2(i,j,k)=-a21(i,j,k)*
& (g11(i,j,k)*x11+2*g12(i,j,k)*x12+g22(i,j,k)*x22
& +2*g13(i,j,k)*x13+2*g23(i,j,k)*x23+g33(i,j,k)*x33)
& -a22(i,j,k)*
& (g11(i,j,k)*y11+2*g12(i,j,k)*y12+g22(i,j,k)*y22
& +2*g13(i,j,k)*y13+2*g23(i,j,k)*y23+g33(i,j,k)*y33)
& -a23(i,j,k)*
& (g11(i,j,k)*z11+2*g12(i,j,k)*z12+g22(i,j,k)*z22
& +2*g13(i,j,k)*z13+2*g23(i,j,k)*z23+g33(i,j,k)*z33)
d2z3(i,j,k)=-a31(i,j,k)*
& (g11(i,j,k)*x11+2*g12(i,j,k)*x12+g22(i,j,k)*x22
& +2*g13(i,j,k)*x13+2*g23(i,j,k)*x23+g33(i,j,k)*x33)
& -a32(i,j,k)*
& (g11(i,j,k)*y11+2*g12(i,j,k)*y12+g22(i,j,k)*y22
& +2*g13(i,j,k)*y13+2*g23(i,j,k)*y23+g33(i,j,k)*y33)
& -a33(i,j,k)*
& (g11(i,j,k)*z11+2*g12(i,j,k)*z12+g22(i,j,k)*z22
& +2*g13(i,j,k)*z13+2*g23(i,j,k)*z23+g33(i,j,k)*z33)
end do
end do
end do
c do n=1,in
c write (*,*) n,dummy(n,n,1)
c end do
dtdif=rmin(dummy)
return
end
subroutine diffus(dudt,dvdt,dwdt) 1
c*
c** updated 11 August to account for variable viscosity
c
*
include 'box.size'
include 'box.common'
dimension dudt(in,jn,kn),dvdt(in,jn,kn),
& dwdt(in,jn,kn)
do 10 k=slip,kn
do 10 j=1,jn
do i=1,in
c---- calculate diffusion term
dnud3=a33(i,j,k)*0.25*(vnu(ku(k))-vnu(kd(k)))
c* u
dud1=0.25*(qr(i)*u(ir(i),j,k)
& -ql(i)*u(il(i),j,k))
dud2=0.25*(u(i,jt(j),k)
& -u(i,jb(j),k))
dud3=0.25*(u(i,j,ku(k))
& -u(i,j,kd(k)))
diffu=
& g11(i,j,k)*0.25*(qr(i)*u(ir(i),j,k)
& -2.*u(i,j,k)
& +ql(i)*u(il(i),j,k))
& +g22(i,j,k)*0.25*(u(i,jt(j),k)
& -2.*u(i,j,k)
& +u(i,jb(j),k))
& +g33(i,j,k)*0.25*(u(i,j,ku(k))
& -2.*u(i,j,k)
& +u(i,j,kd(k)))
& +2.*g12(i,j,k)*0.0625*(
& +(qr(i)*u(ir(i),jt(j),k)
& -ql(i)*u(il(i),jt(j),k))
& -(qr(i)*u(ir(i),jb(j),k)
& -ql(i)*u(il(i),jb(j),k)))
& +2.*g13(i,j,k)*0.0625*(
& +(qr(i)*u(ir(i),j,ku(k))
& -ql(i)*u(il(i),j,ku(k)))
& -(qr(i)*u(ir(i),j,kd(k))
& -ql(i)*u(il(i),j,kd(k))))
& +2.*g23(i,j,k)*0.0625*(
& +(u(i,jt(j),ku(k))
& -u(i,jb(j),ku(k)))
& -(u(i,jt(j),kd(k))
& -u(i,jb(j),kd(k))))
& +d2z1(i,j,k)*dud1
& +d2z2(i,j,k)*dud2
& +d2z3(i,j,k)*dud3
c* v
dvd1=0.25*(v(ir(i),j,k)
& -v(il(i),j,k))
dvd2=0.25*(qt(j)*v(i,jt(j),k)
& -qb(j)*v(i,jb(j),k))
dvd3=0.25*(v(i,j,ku(k))
& -v(i,j,kd(k)))
diffv=
& g11(i,j,k)*0.25*(v(ir(i),j,k)
& -2.*v(i,j,k)
& +v(il(i),j,k))
& +g22(i,j,k)*0.25*(qt(j)*v(i,jt(j),k)
& -2.*v(i,j,k)
& +qb(j)*v(i,jb(j),k))
& +g33(i,j,k)*0.25*(v(i,j,ku(k))
& -2.*v(i,j,k)
& +v(i,j,kd(k)))
& +2.*g12(i,j,k)*0.0625*(
& +(qt(j)*v(ir(i),jt(j),k)
& -qt(j)*v(il(i),jt(j),k))
& -(qb(j)*v(ir(i),jb(j),k)
& -qb(j)*v(il(i),jb(j),k)))
& +2.*g13(i,j,k)*0.0625*(
& +(v(ir(i),j,ku(k))
& -v(il(i),j,ku(k)))
& -(v(ir(i),j,kd(k))
& -v(il(i),j,kd(k))))
& +2.*g23(i,j,k)*0.0625*(
& +(qt(j)*v(i,jt(j),ku(k))
& -qb(j)*v(i,jb(j),ku(k)))
& -(qt(j)*v(i,jt(j),kd(k))
& -qb(j)*v(i,jb(j),kd(k))))
& +d2z1(i,j,k)*dvd1
& +d2z2(i,j,k)*dvd2
& +d2z3(i,j,k)*dvd3
c* w
dwd1=0.25*(w(ir(i),j,k)
& -w(il(i),j,k))
dwd2=0.25*(w(i,jt(j),k)
& -w(i,jb(j),k))
dwd3=0.25*(qu(k)*w(i,j,ku(k))
& -qd(k)*w(i,j,kd(k)))
diffw=
& g11(i,j,k)*0.25*(w(ir(i),j,k)
& -2.*w(i,j,k)
& +w(il(i),j,k))
& +g22(i,j,k)*0.25*(w(i,jt(j),k)
& -2.*w(i,j,k)
& +w(i,jb(j),k))
& +g33(i,j,k)*0.25*(qu(k)*w(i,j,ku(k))
& -2.*w(i,j,k)
& +qd(k)*w(i,j,kd(k)))
& +2.*g12(i,j,k)*0.0625*(
& +(w(ir(i),jt(j),k)
& -w(ir(i),jb(j),k))
& -(w(il(i),jt(j),k)
& -w(il(i),jb(j),k)))
& +2.*g13(i,j,k)*0.0625*(
& +(qu(k)*w(ir(i),j,ku(k))
& -qu(k)*w(il(i),j,ku(k)))
& -(qd(k)*w(ir(i),j,kd(k))
& -qd(k)*w(il(i),j,kd(k))))
& +2.*g23(i,j,k)*0.0625*(
& +(qu(k)*w(i,jt(j),ku(k))
& -qu(k)*w(i,jb(j),ku(k)))
& -(qd(k)*w(i,jt(j),kd(k))
& -qd(k)*w(i,jb(j),kd(k))))
& +d2z1(i,j,k)*dwd1
& +d2z2(i,j,k)*dwd2
& +d2z3(i,j,k)*dwd3
c---- terms due to variable viscosity:
unuterm=(a13(i,j,k)*dud1+a23(i,j,k)*dud2+a33(i,j,k)*dud3
& +a11(i,j,k)*dwd1+a21(i,j,k)*dwd2+a31(i,j,k)*dwd3)*dnud3
vnuterm=(a13(i,j,k)*dvd1+a23(i,j,k)*dvd2+a33(i,j,k)*dvd3
& +a12(i,j,k)*dwd1+a22(i,j,k)*dwd2+a32(i,j,k)*dwd3)*dnud3
wnuterm=(a13(i,j,k)*dwd1+a23(i,j,k)*dwd2+a33(i,j,k)*dwd3)
& *dnud3*2.
c---- add diffusion to other terms ----
dudt(i,j,k)=dudt(i,j,k)+vnu(k)*diffu+unuterm
dvdt(i,j,k)=dvdt(i,j,k)+vnu(k)*diffv+vnuterm
dwdt(i,j,k)=dwdt(i,j,k)+vnu(k)*diffw+wnuterm
enddo
10 continue
return
end
subroutine advect(dudt,dvdt,dwdt,vcmax) 1
c*
c** compute advection and buoyancy terms
c*
include 'box.size'
include 'box.common'
dimension dudt(in,jn,kn),dvdt(in,jn,kn),
& dwdt(in,jn,kn)
common/spare/uc(in,jn,kn),vc(in,jn,kn),wc(in,jn,kn),
# dwd1(in),dwd2(in),dwd3(in),
# dud1(in),dud2(in),dud3(in),dvd1(in),dvd2(in),dvd3(in)
denom=1./120.
factor=1.-exp(-2.*time)
dx=1.e20
do k=1,kn
do j=1,jn
do i=1,in
c---- compute contravariant velocity components
d1=u(i,j,k)
d2=v(i,j,k)
d3=w(i,j,k)
uc(i,j,k)=a11(i,j,k)*d1+a12(i,j,k)*d2+a13(i,j,k)*d3
vc(i,j,k)=a21(i,j,k)*d1+a22(i,j,k)*d2+a23(i,j,k)*d3
wc(i,j,k)=a31(i,j,k)*d1+a32(i,j,k)*d2+a33(i,j,k)*d3
c---- compute horizontal advection terms
c**** 5th order upwind (horizontal and vertical) ****
if(uc(i,j,k).gt.0.) then
dwd1(i)=(-3.*w(ir2(i),j,k)+30.*w(ir(i),j,k)+20.*w(i,j,k)
& -60.*w(il(i),j,k)+15.*w(il2(i),j,k)
& -2.*w(il3(i),j,k))*denom
dvd1(i)=(-3.*v(ir2(i),j,k)+30.*v(ir(i),j,k)+20.*v(i,j,k)
& -60.*v(il(i),j,k)+15.*v(il2(i),j,k)
& -2.*v(il3(i),j,k))*denom
dud1(i)=(-3.*qr2(i)*u(ir2(i),j,k)+30.*qr(i)*u(ir(i),j,k)
& +20.*u(i,j,k)
& -60.*ql(i)*u(il(i),j,k)+15.*ql2(i)*u(il2(i),j,k)
& -2.*ql3(i)*u(il3(i),j,k))*denom
else
dwd1(i)=(2.*w(ir3(i),j,k)-15.*w(ir2(i),j,k)+60.*w(ir(i),j,k)
& -20.*w(i,j,k)-30.*w(il(i),j,k)+3.*w(il2(i),j,k))*denom
dvd1(i)=(2.*v(ir3(i),j,k)-15.*v(ir2(i),j,k)+60.*v(ir(i),j,k)
& -20.*v(i,j,k)-30.*v(il(i),j,k)+3.*v(il2(i),j,k))*denom
dud1(i)=(2.*qr3(i)*u(ir3(i),j,k)-15.*qr2(i)*u(ir2(i),j,k)
& +60.*qr(i)*u(ir(i),j,k)-20.*u(i,j,k)
& -30.*ql(i)*u(il(i),j,k)+3.*ql2(i)*u(il2(i),j,k))*denom
endif
if(vc(i,j,k).gt.0.) then
dwd2(i)=(-3.*w(i,jt2(j),k)+30.*w(i,jt(j),k)+20.*w(i,j,k)
& -60.*w(i,jb(j),k)+15.*w(i,jb2(j),k)
& -2.*w(i,jb3(j),k))*denom
dvd2(i)=(-3.*qt2(j)*v(i,jt2(j),k)+30.*qt(j)*v(i,jt(j),k)
& +20.*v(i,j,k)
& -60.*qb(j)*v(i,jb(j),k)+15.*qb2(j)*v(i,jb2(j),k)
& -2.*qb3(j)*v(i,jb3(j),k))*denom
dud2(i)=(-3.*u(i,jt2(j),k)+30.*u(i,jt(j),k)+20.*u(i,j,k)
& -60.*u(i,jb(j),k)+15.*u(i,jb2(j),k)
& -2.*u(i,jb3(j),k))*denom
else
dwd2(i)=(2.*w(i,jt3(j),k)-15.*w(i,jt2(j),k)+60.*w(i,jt(j),k)
& -20.*w(i,j,k)-30.*w(i,jb(j),k)+3.*w(i,jb2(j),k))*denom
dvd2(i)=(2.*qt3(j)*v(i,jt3(j),k)-15.*qt2(j)*v(i,jt2(j),k)
& +60.*qt(j)*v(i,jt(j),k)-20.*v(i,j,k)
& -30.*qb(j)*v(i,jb(j),k)+3.*qb2(j)*v(i,jb2(j),k))*denom
dud2(i)=(2.*u(i,jt3(j),k)-15.*u(i,jt2(j),k)+60.*u(i,jt(j),k)
& -20.*u(i,j,k)-30.*u(i,jb(j),k)+3.*u(i,jb2(j),k))*denom
endif
c---- combine terms ----
dudt(i,j,k)=-uc(i,j,k)*dud1(i)-vc(i,j,k)*dud2(i)
& +swirl*v(i,j,k)
dvdt(i,j,k)=-uc(i,j,k)*dvd1(i)-vc(i,j,k)*dvd2(i)
& -swirl*u(i,j,k)
dwdt(i,j,k)=-uc(i,j,k)*dwd1(i)-vc(i,j,k)*dwd2(i)
& +b(i,j,k)*factor
enddo
enddo
enddo
c---- FOR VERTICAL LEVELS ABOVE K=3:
do k=4,kn
do j=1,jn
do i=1,in
if(wc(i,j,k).gt.0.) then
dwd3(i)=(-3.*qu2(k)*w(i,j,ku2(k))+30.*qu(k)*w(i,j,ku(k))
& +20.*w(i,j,k)-60.*qd(k)*w(i,j,kd(k))
& +15.*qd2(k)*w(i,j,kd2(k))
& -2.*qd3(k)*w(i,j,kd3(k)))*denom
dvd3(i)=(-3.*v(i,j,ku2(k))+30.*v(i,j,ku(k))+20.*v(i,j,k)
& -60.*v(i,j,kd(k))+15.*v(i,j,kd2(k))
& -2.*v(i,j,kd3(k)))*denom
dud3(i)=(-3.*u(i,j,ku2(k))+30.*u(i,j,ku(k))+20.*u(i,j,k)
& -60.*u(i,j,kd(k))+15.*u(i,j,kd2(k))
& -2.*u(i,j,kd3(k)))*denom
else
dwd3(i)=(2.*qu3(k)*w(i,j,ku3(k))-15.*qu2(k)*w(i,j,ku2(k))
& +60.*qu(k)*w(i,j,ku(k))-20.*w(i,j,k)
& -30.*qd(k)*w(i,j,kd(k))+3.*qd2(k)*w(i,j,kd2(k)))*denom
dvd3(i)=(2.*v(i,j,ku3(k))-15.*v(i,j,ku2(k))+60.*v(i,j,ku(k))
& -20.*v(i,j,k)-30.*v(i,j,kd(k))+3.*v(i,j,kd2(k)))*denom
dud3(i)=(2.*u(i,j,ku3(k))-15.*u(i,j,ku2(k))+60.*u(i,j,ku(k))
& -20.*u(i,j,k)-30.*u(i,j,kd(k))+3.*u(i,j,kd2(k)))*denom
endif
dudt(i,j,k)=dudt(i,j,k)-wc(i,j,k)*dud3(i)
dvdt(i,j,k)=dvdt(i,j,k)-wc(i,j,k)*dvd3(i)
dwdt(i,j,k)=dwdt(i,j,k)-wc(i,j,k)*dwd3(i)
enddo
enddo
enddo
c---- FOR 2ND VERTICAL GRID LEVEL ABOVE SURFACE ----
k=3
do j=1,jn
do i=1,in
if(wc(i,j,k).gt.0.) then
dwd3(i)=(2.*qu(k)*w(i,j,ku(k))+3.*w(i,j,k)
& -6.*qd(k)*w(i,j,kd(k))+qd2(k)*w(i,j,kd2(k)))/12.
dvd3(i)=(2.*v(i,j,ku(k))+3.*v(i,j,k)-6.*v(i,j,kd(k))
& +v(i,j,kd2(k)))/12.
dud3(i)=(2.*u(i,j,ku(k))+3.*u(i,j,k)-6.*u(i,j,kd(k))
& +u(i,j,kd2(k)))/12.
else
dwd3(i)=(-qu2(k)*w(i,j,ku2(k))+6.*qu(k)*w(i,j,ku(k))
& -3.*w(i,j,k)-2.*qd(k)*w(i,j,kd(k)))/12.
dvd3(i)=(-v(i,j,ku2(k))+6.*v(i,j,ku(k))-3.*v(i,j,k)
& -2.*v(i,j,kd(k)))/12.
dud3(i)=(-u(i,j,ku2(k))+6.*u(i,j,ku(k))-3.*u(i,j,k)
& -2.*u(i,j,kd(k)))/12.
endif
dudt(i,j,k)=dudt(i,j,k)-wc(i,j,k)*dud3(i)
dvdt(i,j,k)=dvdt(i,j,k)-wc(i,j,k)*dvd3(i)
dwdt(i,j,k)=dwdt(i,j,k)-wc(i,j,k)*dwd3(i)
enddo
enddo
c---- FOR 1ST VERTICAL GRID LEVEL ABOVE SURFACE:
k=2
do j=1,jn
do i=1,in
if(wc(i,j,k).gt.0.) then
dwd3(i)=(w(i,j,k)-qd(k)*w(i,j,kd(k)))*0.5
dvd3(i)=(v(i,j,k)-v(i,j,kd(k)))*0.5
dud3(i)=(u(i,j,k)-u(i,j,kd(k)))*0.5
else
dwd3(i)=(qu(k)*w(i,j,ku(k))-w(i,j,k))*0.5
dvd3(i)=(v(i,j,ku(k))-v(i,j,k))*0.5
dud3(i)=(u(i,j,ku(k))-u(i,j,k))*0.5
endif
dudt(i,j,k)=dudt(i,j,k)-wc(i,j,k)*dud3(i)
dvdt(i,j,k)=dvdt(i,j,k)-wc(i,j,k)*dvd3(i)
dwdt(i,j,k)=dwdt(i,j,k)-wc(i,j,k)*dwd3(i)
enddo
enddo
vcmax=amax1(rm(uc),rm(vc),rm(wc))
return
end
subroutine advect3(dudt,dvdt,dwdt,vcmax) 1
c*
c** compute advection and buoyancy terms
c*
include 'box.size'
include 'box.common'
dimension dudt(in,jn,kn),dvdt(in,jn,kn),
& dwdt(in,jn,kn)
common/spare/uc(in,jn,kn),vc(in,jn,kn),wc(in,jn,kn),
# dwd1(in),dwd2(in),dwd3(in),
# dud1(in),dud2(in),dud3(in),dvd1(in),dvd2(in),dvd3(in)
denom=1./120.
twl=1./12.
factor=1.-exp(-2.*time)
dx=1.e20
do k=1,kn
do j=1,jn
do i=1,in
c---- compute contravariant velocity components
d1=u(i,j,k)
d2=v(i,j,k)
d3=w(i,j,k)
uc(i,j,k)=a11(i,j,k)*d1+a12(i,j,k)*d2+a13(i,j,k)*d3
vc(i,j,k)=a21(i,j,k)*d1+a22(i,j,k)*d2+a23(i,j,k)*d3
wc(i,j,k)=a31(i,j,k)*d1+a32(i,j,k)*d2+a33(i,j,k)*d3
c---- compute horizontal advection terms
c**** 5th order upwind (horizontal and vertical) ****
if(uc(i,j,k).gt.0.) then
dwd1(i)=(2.*w(ir(i),j,k)+3.*w(i,j,k)
& -6.*w(il(i),j,k)+w(il2(i),j,k))*twl
dvd1(i)=(2.*v(ir(i),j,k)+3.*v(i,j,k)
& -6.*v(il(i),j,k)+v(il2(i),j,k))*twl
dud1(i)=(2.*qr(i)*u(ir(i),j,k)
& +3.*u(i,j,k)
& -6.*ql(i)*u(il(i),j,k)+ql2(i)*u(il2(i),j,k)
& )*twl
else
dwd1(i)=(-w(ir2(i),j,k)+6.*w(ir(i),j,k)
& -3.*w(i,j,k)-2.*w(il(i),j,k))*twl
dvd1(i)=(-v(ir2(i),j,k)+6.*v(ir(i),j,k)
& -3.*v(i,j,k)-2.*v(il(i),j,k))*twl
dud1(i)=(-qr2(i)*u(ir2(i),j,k)
& +6.*qr(i)*u(ir(i),j,k)-3.*u(i,j,k)
& -2.*ql(i)*u(il(i),j,k))*twl
endif
if(vc(i,j,k).gt.0.) then
dwd2(i)=(2.*w(i,jt(j),k)+3.*w(i,j,k)
& -6.*w(i,jb(j),k)+1.*w(i,jb2(j),k)
& )*twl
dvd2(i)=(2.*qt(j)*v(i,jt(j),k)
& +3.*v(i,j,k)
& -6.*qb(j)*v(i,jb(j),k)+1.*qb2(j)*v(i,jb2(j),k)
& )*twl
dud2(i)=(2.*u(i,jt(j),k)+3.*u(i,j,k)
& -6.*u(i,jb(j),k)+1.*u(i,jb2(j),k)
& )*twl
else
dwd2(i)=(-w(i,jt2(j),k)+6.*w(i,jt(j),k)
& -3.*w(i,j,k)-2.*w(i,jb(j),k))*twl
dvd2(i)=(-1.*qt2(j)*v(i,jt2(j),k)
& +6.*qt(j)*v(i,jt(j),k)-3.*v(i,j,k)
& -2.*qb(j)*v(i,jb(j),k))*twl
dud2(i)=(-u(i,jt2(j),k)+6.*u(i,jt(j),k)
& -3.*u(i,j,k)-2.*u(i,jb(j),k))*twl
endif
c---- combine terms ----
dudt(i,j,k)=-uc(i,j,k)*dud1(i)-vc(i,j,k)*dud2(i)
& +swirl*v(i,j,k)
dvdt(i,j,k)=-uc(i,j,k)*dvd1(i)-vc(i,j,k)*dvd2(i)
& -swirl*u(i,j,k)
dwdt(i,j,k)=-uc(i,j,k)*dwd1(i)-vc(i,j,k)*dwd2(i)
& +b(i,j,k)*factor
enddo
enddo
enddo
c---- FOR VERTICAL LEVELS ABOVE K=2:
do k=3,kn
do j=1,jn
do i=1,in
if(wc(i,j,k).gt.0.) then
dwd3(i)=(2.*qu(k)*w(i,j,ku(k))+3.*w(i,j,k)
& -6.*qd(k)*w(i,j,kd(k))+qd2(k)*w(i,j,kd2(k)))/12.
dvd3(i)=(2.*v(i,j,ku(k))+3.*v(i,j,k)-6.*v(i,j,kd(k))
& +v(i,j,kd2(k)))/12.
dud3(i)=(2.*u(i,j,ku(k))+3.*u(i,j,k)-6.*u(i,j,kd(k))
& +u(i,j,kd2(k)))/12.
else
dwd3(i)=(-qu2(k)*w(i,j,ku2(k))+6.*qu(k)*w(i,j,ku(k))
& -3.*w(i,j,k)-2.*qd(k)*w(i,j,kd(k)))/12.
dvd3(i)=(-v(i,j,ku2(k))+6.*v(i,j,ku(k))-3.*v(i,j,k)
& -2.*v(i,j,kd(k)))/12.
dud3(i)=(-u(i,j,ku2(k))+6.*u(i,j,ku(k))-3.*u(i,j,k)
& -2.*u(i,j,kd(k)))/12.
endif
dudt(i,j,k)=dudt(i,j,k)-wc(i,j,k)*dud3(i)
dvdt(i,j,k)=dvdt(i,j,k)-wc(i,j,k)*dvd3(i)
dwdt(i,j,k)=dwdt(i,j,k)-wc(i,j,k)*dwd3(i)
enddo
enddo
enddo
c---- FOR 1ST VERTICAL GRID LEVEL ABOVE SURFACE:
k=2
do j=1,jn
do i=1,in
if(wc(i,j,k).gt.0.) then
dwd3(i)=(w(i,j,k)-qd(k)*w(i,j,kd(k)))*0.5
dvd3(i)=(v(i,j,k)-v(i,j,kd(k)))*0.5
dud3(i)=(u(i,j,k)-u(i,j,kd(k)))*0.5
else
dwd3(i)=(qu(k)*w(i,j,ku(k))-w(i,j,k))*0.5
dvd3(i)=(v(i,j,ku(k))-v(i,j,k))*0.5
dud3(i)=(u(i,j,ku(k))-u(i,j,k))*0.5
endif
dudt(i,j,k)=dudt(i,j,k)-wc(i,j,k)*dud3(i)
dvdt(i,j,k)=dvdt(i,j,k)-wc(i,j,k)*dvd3(i)
dwdt(i,j,k)=dwdt(i,j,k)-wc(i,j,k)*dwd3(i)
enddo
enddo
vcmax=amax1(rm(uc),rm(vc),rm(wc))
return
end
subroutine incmp(dudt,dvdt,dwdt) 1,2
c*
c** solve for pressure of incompressible flow -- invert laplcian using jacobi
c** iteration
c*
include 'box.size'
include 'box.common'
dimension dudt(in,jn,kn),dvdt(in,jn,kn),dwdt(in,jn,kn)
common/spare/dqdt(in,jn,kn),drdt(in,jn,kn),dsdt(in,jn,kn)
# ,dum(im,jm,km)
c rdt=0.5/dt
rdt=gamma/dt
c---- continuity equation ----
do k=1,km
do j=1,jm
do i=1,im
divu1=0.125*(
& u(irc(i),jtc(j),kuc(k))-u(ilc(i),jtc(j),kuc(k))
& +u(irc(i),jbc(j),kuc(k))-u(ilc(i),jbc(j),kuc(k))
& +u(irc(i),jtc(j),kdc(k))-u(ilc(i),jtc(j),kdc(k))
& +u(irc(i),jbc(j),kdc(k))-u(ilc(i),jbc(j),kdc(k)))
divu2=0.125*(
& u(irc(i),jtc(j),kuc(k))-u(irc(i),jbc(j),kuc(k))
& +u(ilc(i),jtc(j),kuc(k))-u(ilc(i),jbc(j),kuc(k))
& +u(irc(i),jtc(j),kdc(k))-u(irc(i),jbc(j),kdc(k))
& +u(ilc(i),jtc(j),kdc(k))-u(ilc(i),jbc(j),kdc(k)))
divu3=0.125*(
& u(irc(i),jtc(j),kuc(k))-u(irc(i),jtc(j),kdc(k))
& +u(ilc(i),jtc(j),kuc(k))-u(ilc(i),jtc(j),kdc(k))
& +u(irc(i),jbc(j),kuc(k))-u(irc(i),jbc(j),kdc(k))
& +u(ilc(i),jbc(j),kuc(k))-u(ilc(i),jbc(j),kdc(k)))
c*
divv1=0.125*(
& v(irc(i),jtc(j),kuc(k))-v(ilc(i),jtc(j),kuc(k))
& +v(irc(i),jbc(j),kuc(k))-v(ilc(i),jbc(j),kuc(k))
& +v(irc(i),jtc(j),kdc(k))-v(ilc(i),jtc(j),kdc(k))
& +v(irc(i),jbc(j),kdc(k))-v(ilc(i),jbc(j),kdc(k)))
divv2=0.125*(
& v(irc(i),jtc(j),kuc(k))-v(irc(i),jbc(j),kuc(k))
& +v(ilc(i),jtc(j),kuc(k))-v(ilc(i),jbc(j),kuc(k))
& +v(irc(i),jtc(j),kdc(k))-v(irc(i),jbc(j),kdc(k))
& +v(ilc(i),jtc(j),kdc(k))-v(ilc(i),jbc(j),kdc(k)))
divv3=0.125*(
& v(irc(i),jtc(j),kuc(k))-v(irc(i),jtc(j),kdc(k))
& +v(ilc(i),jtc(j),kuc(k))-v(ilc(i),jtc(j),kdc(k))
& +v(irc(i),jbc(j),kuc(k))-v(irc(i),jbc(j),kdc(k))
& +v(ilc(i),jbc(j),kuc(k))-v(ilc(i),jbc(j),kdc(k)))
c*
divw1=0.125*(
& w(irc(i),jtc(j),kuc(k))-w(ilc(i),jtc(j),kuc(k))
& +w(irc(i),jbc(j),kuc(k))-w(ilc(i),jbc(j),kuc(k))
& +w(irc(i),jtc(j),kdc(k))-w(ilc(i),jtc(j),kdc(k))
& +w(irc(i),jbc(j),kdc(k))-w(ilc(i),jbc(j),kdc(k)))
divw2=0.125*(
& w(irc(i),jtc(j),kuc(k))-w(irc(i),jbc(j),kuc(k))
& +w(ilc(i),jtc(j),kuc(k))-w(ilc(i),jbc(j),kuc(k))
& +w(irc(i),jtc(j),kdc(k))-w(irc(i),jbc(j),kdc(k))
& +w(ilc(i),jtc(j),kdc(k))-w(ilc(i),jbc(j),kdc(k)))
divw3=0.125*(
& w(irc(i),jtc(j),kuc(k))-w(irc(i),jtc(j),kdc(k))
& +w(ilc(i),jtc(j),kuc(k))-w(ilc(i),jtc(j),kdc(k))
& +w(irc(i),jbc(j),kuc(k))-w(irc(i),jbc(j),kdc(k))
& +w(ilc(i),jbc(j),kuc(k))-w(ilc(i),jbc(j),kdc(k)))
dum(i,j,k)=rdt*
& (c11(i,j,k)*divu1+c21(i,j,k)*divu2+c31(i,j,k)*divu3
& +c12(i,j,k)*divv1+c22(i,j,k)*divv2+c32(i,j,k)*divv3
& +c13(i,j,k)*divw1+c23(i,j,k)*divw2+c33(i,j,k)*divw3)
enddo
enddo
enddo
do 100 sn=1.,pit
call pforce(dqdt,dudt,drdt,dvdt,dsdt,dwdt)
do k=1,km
do j=1,jm
do i=1,im
divu1=0.125*(
& dqdt(irc(i),jtc(j),kuc(k))-dqdt(ilc(i),jtc(j),kuc(k))
& +dqdt(irc(i),jbc(j),kuc(k))-dqdt(ilc(i),jbc(j),kuc(k))
& +dqdt(irc(i),jtc(j),kdc(k))-dqdt(ilc(i),jtc(j),kdc(k))
& +dqdt(irc(i),jbc(j),kdc(k))-dqdt(ilc(i),jbc(j),kdc(k)))
divu2=0.125*(
& dqdt(irc(i),jtc(j),kuc(k))-dqdt(irc(i),jbc(j),kuc(k))
& +dqdt(ilc(i),jtc(j),kuc(k))-dqdt(ilc(i),jbc(j),kuc(k))
& +dqdt(irc(i),jtc(j),kdc(k))-dqdt(irc(i),jbc(j),kdc(k))
& +dqdt(ilc(i),jtc(j),kdc(k))-dqdt(ilc(i),jbc(j),kdc(k)))
divu3=0.125*(
& dqdt(irc(i),jtc(j),kuc(k))-dqdt(irc(i),jtc(j),kdc(k))
& +dqdt(ilc(i),jtc(j),kuc(k))-dqdt(ilc(i),jtc(j),kdc(k))
& +dqdt(irc(i),jbc(j),kuc(k))-dqdt(irc(i),jbc(j),kdc(k))
& +dqdt(ilc(i),jbc(j),kuc(k))-dqdt(ilc(i),jbc(j),kdc(k)))
c*
divv1=0.125*(
& drdt(irc(i),jtc(j),kuc(k))-drdt(ilc(i),jtc(j),kuc(k))
& +drdt(irc(i),jbc(j),kuc(k))-drdt(ilc(i),jbc(j),kuc(k))
& +drdt(irc(i),jtc(j),kdc(k))-drdt(ilc(i),jtc(j),kdc(k))
& +drdt(irc(i),jbc(j),kdc(k))-drdt(ilc(i),jbc(j),kdc(k)))
divv2=0.125*(
& drdt(irc(i),jtc(j),kuc(k))-drdt(irc(i),jbc(j),kuc(k))
& +drdt(ilc(i),jtc(j),kuc(k))-drdt(ilc(i),jbc(j),kuc(k))
& +drdt(irc(i),jtc(j),kdc(k))-drdt(irc(i),jbc(j),kdc(k))
& +drdt(ilc(i),jtc(j),kdc(k))-drdt(ilc(i),jbc(j),kdc(k)))
divv3=0.125*(
& drdt(irc(i),jtc(j),kuc(k))-drdt(irc(i),jtc(j),kdc(k))
& +drdt(ilc(i),jtc(j),kuc(k))-drdt(ilc(i),jtc(j),kdc(k))
& +drdt(irc(i),jbc(j),kuc(k))-drdt(irc(i),jbc(j),kdc(k))
& +drdt(ilc(i),jbc(j),kuc(k))-drdt(ilc(i),jbc(j),kdc(k)))
c*
divw1=0.125*(
& dsdt(irc(i),jtc(j),kuc(k))-dsdt(ilc(i),jtc(j),kuc(k))
& +dsdt(irc(i),jbc(j),kuc(k))-dsdt(ilc(i),jbc(j),kuc(k))
& +dsdt(irc(i),jtc(j),kdc(k))-dsdt(ilc(i),jtc(j),kdc(k))
& +dsdt(irc(i),jbc(j),kdc(k))-dsdt(ilc(i),jbc(j),kdc(k)))
divw2=0.125*(
& dsdt(irc(i),jtc(j),kuc(k))-dsdt(irc(i),jbc(j),kuc(k))
& +dsdt(ilc(i),jtc(j),kuc(k))-dsdt(ilc(i),jbc(j),kuc(k))
& +dsdt(irc(i),jtc(j),kdc(k))-dsdt(irc(i),jbc(j),kdc(k))
& +dsdt(ilc(i),jtc(j),kdc(k))-dsdt(ilc(i),jbc(j),kdc(k)))
divw3=0.125*(
& dsdt(irc(i),jtc(j),kuc(k))-dsdt(irc(i),jtc(j),kdc(k))
& +dsdt(ilc(i),jtc(j),kuc(k))-dsdt(ilc(i),jtc(j),kdc(k))
& +dsdt(irc(i),jbc(j),kuc(k))-dsdt(irc(i),jbc(j),kdc(k))
& +dsdt(ilc(i),jbc(j),kuc(k))-dsdt(ilc(i),jbc(j),kdc(k)))
p(i,j,k)=p(i,j,k)-alpha(i,j,k)*(dum(i,j,k)+
& (c11(i,j,k)*divu1+c21(i,j,k)*divu2+c31(i,j,k)*divu3
& +c12(i,j,k)*divv1+c22(i,j,k)*divv2+c32(i,j,k)*divv3
& +c13(i,j,k)*divw1+c23(i,j,k)*divw2+c33(i,j,k)*divw3))
enddo
enddo
enddo
100 continue
call pforce(dudt,dudt,dvdt,dvdt,dwdt,dwdt)
return
end
subroutine pforce(dqdt,dudt,drdt,dvdt,dsdt,dwdt) 2
include 'box.size'
include 'box.common'
dimension dudt(in,jn,kn),dvdt(in,jn,kn),
& dwdt(in,jn,kn),dqdt(in,jn,kn),drdt(in,jn,kn),
& dsdt(in,jn,kn)
do k=slip,kn
do j=1,jn
do i=1,in
c---- compute pressure gradient terms ----
dpdxi1=0.125*(
& +(p(iru(i),jtu(j),kuu(k))-p(ilu(i),jtu(j),kuu(k)))
& +(p(iru(i),jbu(j),kuu(k))-p(ilu(i),jbu(j),kuu(k)))
& +(p(iru(i),jtu(j),kdu(k))-p(ilu(i),jtu(j),kdu(k)))
& +(p(iru(i),jbu(j),kdu(k))-p(ilu(i),jbu(j),kdu(k))))
dpdxi2=0.125*(
& +(p(iru(i),jtu(j),kuu(k))-p(iru(i),jbu(j),kuu(k)))
& +(p(ilu(i),jtu(j),kuu(k))-p(ilu(i),jbu(j),kuu(k)))
& +(p(iru(i),jtu(j),kdu(k))-p(iru(i),jbu(j),kdu(k)))
& +(p(ilu(i),jtu(j),kdu(k))-p(ilu(i),jbu(j),kdu(k))))
dpdxi3=0.125*(
& +(p(iru(i),jtu(j),kuu(k))-p(iru(i),jtu(j),kdu(k)))
& +(p(ilu(i),jtu(j),kuu(k))-p(ilu(i),jtu(j),kdu(k)))
& +(p(iru(i),jbu(j),kuu(k))-p(iru(i),jbu(j),kdu(k)))
& +(p(ilu(i),jbu(j),kuu(k))-p(ilu(i),jbu(j),kdu(k))))
c---- combine terms ----
dqdt(i,j,k)=dudt(i,j,k)
& -(a11(i,j,k)*dpdxi1+a21(i,j,k)*dpdxi2
& +a31(i,j,k)*dpdxi3)
drdt(i,j,k)=dvdt(i,j,k)
& -(a12(i,j,k)*dpdxi1+a22(i,j,k)*dpdxi2
& +a32(i,j,k)*dpdxi3)
dsdt(i,j,k)=dwdt(i,j,k)
& -(a13(i,j,k)*dpdxi1+a23(i,j,k)*dpdxi2
& +a33(i,j,k)*dpdxi3)
enddo
enddo
enddo
return
end
function rmin(t) 1
include 'box.size'
dimension t(in,jn,kn)
rmin=t(1,1,1)
do 10 k=1,kn
do 10 j=1,jn
do 10 i=1,in
10 rmin=amin1(rmin,t(i,j,k))
return
end
function rm(t)
c***vectorizable maximum finder
include 'box.size'
dimension t(in,jn,kn)
dimension ts(jn,kn),tt(kn)
do 10 k=1,kn
do 20 j=1,jn
20 ts(j,k)=abs(t(1,j,k))
do 40 i=2,in
do 40 j=1,jn
40 ts(j,k)=amax1(ts(j,k),abs(t(i,j,k)))
10 continue
do 50 k=1,kn
50 tt(k)=ts(1,k)
do 60 j=2,jn
do 60 k=1,kn
60 tt(k)=amax1(ts(j,k),tt(k))
rm=tt(1)
do 70 k=2,kn
70 rm=amax1(rm,tt(k))
return
end
function rmll(t,kk)
c---- diagnose maximum for "low levels"
include 'box.size'
dimension t(in,jn,kn)
rmll=-9999.9
do 10 k=1,kk
c do 10 j=1,jn
c do 10 i=1,in
c changed July 8, 2008
do 10 j=jn/4,3*jn/4
do 10 i=in/4,3*in/4
10 rmll=amax1(rmll,t(i,j,k))
return
end
function rmpll(t,kk)
c---- diagnose minimum pressure for "low levels"
include 'box.size'
dimension t(im,jm,km)
rmpll=+9999.9
do 10 k=1,kk
c do 10 j=1,jm
c do 10 i=1,im
c changed July 8, 2008
do 10 j=jm/4,3*jm/4
do 10 i=im/4,3*im/4
10 rmpll=amin1(rmpll,t(i,j,k))
return
end
function rmp(t)
include 'box.size'
dimension t(im,jm,km)
rmp=t(1,1,1)
do 10 k=1,km
do 10 j=1,jm
do 10 i=1,im
10 rmp=amax1(rmp,abs(t(i,j,k)))
return
end
subroutine makebnds 1
include 'box.size'
include 'box.common'
do i=1,in
ir(i)=i+1
il(i)=i-1
ir2(i)=i+2
il2(i)=i-2
ir3(i)=i+3
il3(i)=i-3
qr(i)=1.
ql(i)=1.
qr2(i)=1.
ql2(i)=1.
qr3(i)=1.
ql3(i)=1.
tr(i)=1.
tl(i)=1.
if(i.eq.1) then
il(i)=ir(i)
ql(i)=-1.
il2(i)=ir2(i)
ql2(i)=-1.
il3(i)=ir3(i)
ql3(i)=-1.
tr(i)=2.
tl(i)=0.
else if(i.eq.2) then
il2(i)=i
ql2(i)=-1.
il3(i)=i+1
ql3(i)=-1.
else if(i.eq.3) then
il3(i)=i-1
ql3(i)=-1.
else if(i.eq.in-2) then
ir3(i)=i+1
qr3(i)=-1
else if(i.eq.in-1) then
ir2(i)=i
qr2(i)=-1.
ir3(i)=i-1
qr3(i)=-1.
else if(i.eq.in) then
ir(i)=il(i)
qr(i)=-1.
ir2(i)=il2(i)
qr2(i)=-1.
ir3(i)=il3(i)
qr3(i)=-1.
tr(i)=0.
tl(i)=2.
endif
enddo
do j=1,jn
jt(j)=j+1
jb(j)=j-1
jt2(j)=j+2
jb2(j)=j-2
jt3(j)=j+3
jb3(j)=j-3
qt(j)=1.
qb(j)=1.
qt2(j)=1.
qb2(j)=1.
qt3(j)=1.
qb3(j)=1.
tt(j)=1.
tb(j)=1.
if(j.eq.1) then
jb(j)=jt(j)
qb(j)=-1.
jb2(j)=jt2(j)
qb2(j)=-1.
jb3(j)=jt3(j)
qb3(j)=-1.
tt(j)=2.
tb(j)=0.
else if(j.eq.2) then
jb2(j)=j
qb2(j)=-1.
jb3(j)=j+1
qb3(j)=-1.
else if(j.eq.3) then
jb3(j)=j-1
qb3(j)=-1.
else if(j.eq.jn-2) then
jt3(j)=j+1
qt3(j)=-1.
else if(j.eq.jn-1) then
jt2(j)=j
qt2(j)=-1.
jt3(j)=j-1
qt3(j)=-1.
else if(j.eq.jn) then
jt(j)=jb(j)
qt(j)=-1.
jt2(j)=jb2(j)
qt2(j)=-1.
jt3(j)=jb3(j)
qt3(j)=-1.
tt(j)=0.
tb(j)=2.
endif
enddo
do k=1,kn
ku(k)=k+1
kd(k)=k-1
ku2(k)=k+2
kd2(k)=k-2
ku3(k)=k+3
kd3(k)=k-3
qu(k)=1.
qd(k)=1.
qu2(k)=1.
qd2(k)=1.
qu3(k)=1.
qd3(k)=1.
tu(k)=1.
td(k)=1.
if(k.eq.1) then
kd(k)=ku(k)
qd(k)=-1.
kd2(k)=ku2(k)
qd2(k)=-1.
kd3(k)=ku3(k)
qd3(k)=-1.
tu(k)=2.
td(k)=0.
else if(k.eq.2) then
kd2(k)=k
qd2(k)=-1.
kd3(k)=k+1
qd3(k)=-1.
else if(k.eq.3) then
kd3(k)=k-1
qd3(k)=-1.
else if(k.eq.kn-2) then
ku3(k)=k+1
qu3(k)=-1.
else if(k.eq.kn-1) then
ku2(k)=k
qu2(k)=-1.
ku3(k)=k-1
qu3(k)=-1.
else if(k.eq.kn) then
ku(k)=kd(k)
qu(k)=-1.
ku2(k)=kd2(k)
qu2(k)=-1.
ku3(k)=kd3(k)
qu3(k)=-1.
tu(k)=0.
td(k)=2.
endif
enddo
do i=1,in
iru(i)=i
ilu(i)=i-1
if(i.eq.in) iru(i)=ilu(i)
if(i.eq.1) ilu(i)=iru(i)
enddo
do j=1,jn
jtu(j)=j
jbu(j)=j-1
if(j.eq.jn) jtu(j)=jbu(j)
if(j.eq.1) jbu(j)=jtu(j)
enddo
do k=1,kn
kuu(k)=k
kdu(k)=k-1
if(k.eq.kn) kuu(k)=kdu(k)
if(k.eq.1) kdu(k)=kuu(k)
enddo
do i=1,im
ilc(i)=i
irc(i)=i+1
enddo
do j=1,jm
jbc(j)=j
jtc(j)=j+1
enddo
do k=1,km
kdc(k)=k
kuc(k)=k+1
enddo
do i=1,in
ild(i)=i-1
icd(i)=i
ird(i)=i+1
if(i.eq.1) then
ild(i)=i
icd(i)=i+1
ird(i)=i+2
endif
if(i.eq.in) then
ild(i)=i
icd(i)=i-1
ird(i)=i-2
endif
enddo
do j=1,jn
jbd(j)=j-1
jcd(j)=j
jtd(j)=j+1
if(j.eq.1) then
jbd(j)=j
jcd(j)=j+1
jtd(j)=j+2
endif
if(j.eq.jn) then
jbd(j)=j
jcd(j)=j-1
jtd(j)=j-2
endif
enddo
do k=1,kn
kdd(k)=k-1
kcd(k)=k
kud(k)=k+1
if(k.eq.1) then
kdd(k)=k
kcd(k)=k+1
kud(k)=k+2
endif
if(k.eq.kn) then
kdd(k)=k
kcd(k)=k-1
kud(k)=k-2
endif
enddo
do i=1,im
irp(i)=i+1
ilp(i)=i-1
if(i.eq.1) then
ilp(i)=irp(i)
else if(i.eq.im) then
irp(i)=ilp(i)
endif
enddo
do j=1,jm
jtp(j)=j+1
jbp(j)=j-1
if(j.eq.1) then
jbp(j)=jtp(j)
else if(j.eq.jm) then
jtp(j)=jbp(j)
endif
enddo
do k=1,km
kup(k)=k+1
kdp(k)=k-1
if(k.eq.1) then
kdp(k)=kup(k)
else if(k.eq.km) then
kup(k)=kdp(k)
endif
enddo
return
end