USE mg3d.cpl
FILE ALTEZZE
ALTEZZE=CREATE('sto')
nx=2∗8
ny=nx∗6
nz=2∗nx
hmax=1.
INTEGER ny0
ny0=FLOOR(2∗nx∗hmax) +5
pi=ATAN(1)∗4
pfo=0.5
pbc=2
protra = 2
prolo = 3
inco=1
REAL fx0
IF protra=1 THEN
fx0 =0
ELSE IF protra=2 THEN
fx0= 0.05
ELSE IF protra=3 THEN
fx0= 0.25
ELSE IF protra=4 THEN
fx0= pi/2
END IF
INTEGER it,js,wp,btm(-1..nx+1,-1..nz+1),btf(0..nx,0..nz)
REAL ad(0..nx,0..ny0,0..nz),al(0..nx,0..ny0,0..nz),ar(0..nx,0..ny0,0..nz), \
ab(0..nx,0..ny0,0..nz),af(0..nx,0..ny0,0..nz),rhs(0..nx,0..ny,0..nz)
REAL up(0..nx,0..ny0,0..nz),right(0..nx,0..ny0,0..nz),down(0..nx,0..ny0,0..nz),\
left(0..nx,0..ny0,0..nz),back(0..nx,0..ny0,0..nz),forw(0..nx,0..ny0,0..nz)
REAL bcu(0..nx,0..ny0,0..nz),bcr(0..nx,0..ny0,0..nz),bcd(0..nx,0..ny0,0..nz), \
bcl(0..nx,0..ny0,0..nz),bcb(0..nx,0..ny0,0..nz),bcf(0..nx,0..ny0,0..nz)
REAL coenoto(0..nx,0..ny0,0..nz)
REAL m1(-1..nx+1,0..ny,-1..nz+1),m2(-1..nx+1,0..ny,-1..nz+1), \
m3(-1..nx+1,0..ny,-1..nz+1),fi(-1..nx+1,0..ny,-1..nz+1)
REAL vzn(1..1500),vsalfa(1..1500),vcalfa(1..1500),vdel(1..1500),vk1(1..1500),vk2(1..1500)
REAL zn,salfa,calfa,del,k1,k2,c2,c4,c6,c2l
INTEGER nid(0..nx,0..ny0,0..nz)
REAL fill,fillm,fils,deno,rapp,hlon,htra,dh,prec,ka
REAL denx,deny,denz,ris(1..50,1..50)
REAL ax0,axl,axr,ay0,ayd,ayu,az0,azf,azb
INTEGER datof
jsta=9
numh=1
numk=1
numf=16
k0=1
delth=0.200
deltf=0.050
fil0=0.10-deltf
rap0=0.300- delth
lx=1.
h=lx/nx/2
ly=ny∗h
lz=nz∗h∗4.
eps=h/100
dcor=2.05∗h
REAL FUNCTION tanh( REAL xx)
RESULT=(EXP(xx)-EXP(-xx))/(EXP(xx)+EXP(-xx))
END tanh
REAL FUNCTION fx( REAL xx)
REAL xs
xs=0.
IF protra = 1 THEN
RESULT=1
IF xx > xs THEN RESULT=0
ELSE IF protra = 2 THEN
RESULT=(1-xx∗4)**5
IF xx>0.25 THEN RESULT= 0
ELSE IF protra = 3 THEN
RESULT=1-xx∗4
IF xx>0.25 THEN RESULT= 0
ELSE IF protra = 4 THEN
RESULT=COS(2∗pi∗xx)
IF xx>0.25 THEN RESULT= 0
END IF
IF RESULT <0 THEN RESULT=0
END fx
REAL FUNCTION gz(REAL zz)
REAL l1,l2,l3,l4,zr,a1,a2,a3,b1,b2,y1
nr=4
IF prolo =1 THEN
RESULT=rapp
ELSE IF prolo=2 THEN
RESULT=rapp∗(1.+tanh(ka∗2./rapp∗(fill∗lz/2.-zz)))/2
ELSE IF prolo=3 THEN
tal=1
cal=1/(1+tal**2)**0.5
sal=tal∗cal
r0=0.347826
l0=fill
l2=(rapp-r0+tal∗l0+r0∗(1+tal**2)**0.5)/tal
l1=l2-(rapp-r0+r0/(1+tal**2)**0.5)/tal
y1=rapp-r0+r0∗cal
l1=r0∗sal+l0
l2=l1+y1/tal
l3=l1+(rapp-2∗r0∗(1-cal))/tal
l4=l3+r0∗sal
b1=(l2-fill+0.125-(nr+1)∗rapp)/(l2-fill+0.125)**nr
b2=(nr∗rapp-l2+fill-0.125)/(l2-fill+0.125)**(nr+1)
IF zz <l0 THEN
RESULT =rapp
ELSE IF zz <l1 THEN
RESULT=rapp-r0+(r0**2-(zz-l0)**2)**0.5
ELSE IF zz<l2 THEN
RESULT = y1+tal∗(l1-zz)
ELSE IF zz<lz-l2 THEN
RESULT=0
ELSE IF zz<lz-l1 THEN
RESULT = y1+tal∗(zz-lz+l1)
ELSE IF zz<lz-l0 THEN
RESULT=rapp-r0+(r0**2-(zz-lz+l0)**2)**0.5
ELSE IF zz <lz THEN
RESULT =rapp
END IF
IF RESULT<0 THEN RESULT=0
END IF
END gz
REAL FUNCTION fondo(REAL xx,zz)
RESULT=fx(xx)∗gz(zz)+fx(0.5-xx)∗gz(lz/2.-zz)
END fondo
REAL FUNCTION vert(REAL yy,x1,x2,zz)
REAL f1,f2,xm,fm,x1c,x2c
f1=fondo(x1,zz)-yy
f2=fondo(x2,zz)-yy
IF f1∗f2 > 0 THEN
IF ABS(f1) < ABS(f2) THEN
RESULT=x1
ELSE
RESULT=x2
END IF
EXIT vert
END IF
x1c=x1
x2c=x2
DO
xm=(f2∗x1c-f1∗x2c)/(f2-f1)
fm=fondo(xm,zz)-yy
IF fm∗f1 > 0 THEN
x1c=xm
f1=fm
ELSE
x2c=xm
f2=fm
END IF
UNTIL (ABS(fm) < 0.0001) OR (ABS(x2c-x1c)) < 0.00005
RESULT=xm
END vert
REAL FUNCTION vertz(REAL yy,z1,z2,xx)
REAL f1,f2,zm,fm,z1c,z2c
f1=fondo(xx,z1)-yy
f2=fondo(xx,z2)-yy
IF(f1∗f2 > 0) THEN
IF(ABS(f1) < ABS(f2)) THEN
RESULT =z1
ELSE
RESULT =z2
END IF
EXIT vertz
END IF
z1c=z1
z2c=z2
DO
zm=(f2∗z1c-f1∗z2c)/(f2-f1)
fm=fondo(xx,zm)-yy
IF (fm∗f1 > 0) THEN
z1c=zm
f1=fm
ELSE
z2c=zm
f2=fm
END IF
UNTIL (ABS(fm) < 0.0001) OR (ABS(z2c-z1c)) <0.00005
RESULT=zm
END vertz
REAL FUNCTION delta(REAL xx,zz,calf,salf)
REAL dz,dzp,xv,talfa,fp,gp,pv
IF xx <0.25 THEN
xv=0
pv=1
ELSE
xv=0.5
pv=-1
END IF
IF calf = 0 THEN
talfa=100000000000000.
ELSE
talfa=salf/calf
END IF
fp=(fondo(xv,zz)-fondo(xv+pv∗eps,zz))/eps
IF xx=0 THEN
gp=0
ELSE
gp=(fondo(xv,zz+eps)-fondo(xv,zz))/eps
END IF
IF rapp=0 THEN
RESULT=pi/2
ELSE
RESULT=ATAN(fx0/rapp∗(1-gp∗talfa)**0.5)
END IF
END delta
REAL FUNCTION kappa(REAL xx, zz,calf,salf)
REAL cot,s,e,kv,kn
IF ABS(delta(xx,zz,calf,salf)) = 0 THEN
RESULT=0.5
ELSE
kn=0.5
s=SIN(2.∗delta(xx,zz,calf,salf))
DO
kv=kn
IF wp = 1 THEN
cot=((1.-s**2∗kv**2)**0.5-1.)/s/kv
kn=(ATAN(1./cot)+pi)/(pi-delta(xx,zz,calf,salf))
ELSE
cot=(1-(1.-s**2∗kv**2)**0.5)/s/kv
kn=ATAN(1./cot)/(pi-delta(xx,zz,calf,salf))
END IF
WHILE ABS(kv-kn) > 0.0001
RESULT=kn
END IF
END kappa
SUBROUTINE catalogo
REAL xx,yy,zz,dis,zve,znu,zno,z1,z2,d1,d2,dm,xv,yv,calf,salf,dpm,z3,norma
INTEGER kv,kvm,kid
kid=0
LOOP FOR ica=0 TO nid.HI
LOOP FOR jca=0 TO nid(0).HI
LOOP FOR kca=0 TO nid(0,0).HI
nid(ica,jca,kca)=0
REPEAT LOOP
REPEAT LOOP
REPEAT LOOP
LOOP FOR ica=0 TO nx
LOOP FOR kca=0 TO nz
IF ica < 2 OR ica > nx-2 THEN
LOOP FOR jca=btm(ica,kca) TO btm(0,0)+2
xx=ica∗h
yy=jca∗h
zz=kca∗h
IF ica∗h < 0.25 THEN
xv=0.
ELSE
xv=0.5
END IF
IF kca > 0 THEN
z1=0.
z2=2.∗nz∗h
d1=(xx-xv)**2+(yy-fondo(xv,z1))**2+(zz-z1)**2
d2=(xx-xv)**2+(yy-fondo(xv,z2))**2+(zz-z2)**2
DO
zno=0.5∗z1+0.5∗z2
dm=(xx-xv)**2+(yy-fondo(xv,zno))**2+(zz-zno)**2
dpm=(xx-xv)**2+(yy-fondo(xv,zno+eps))**2+(zz-zno-eps)**2-
((xx-xv)**2+(yy-fondo(xv,zno-eps))**2+(zz-zno+eps)**2)
IF dpm > 0 THEN
z2=zno
d2=dm
ELSE
z1=zno
d1=dm
END IF
WHILE ABS(z1-z2) > h/1000000
ELSE
zno=0.
calf=1.0
salf=0.
dis=((xx-xv)**2+(yy-fondo(xv,0.))**2)**0.5
END IF
yv=fondo(xv,zno)
dis=((xx-xv)**2+(yy-yv)**2+(zz-zno)**2)**0.5
IF yv > 1.05∗h AND jca>1 THEN
IF dis < dcor THEN
norma=((yy-yv)**2+(zz-zno)**2)**0.5
IF yy < yv THEN norma=-norma
IF ABS(norma) < 0.001 THEN
calf=1.
salf=0.
ELSE
calf=ABS((yy-yv)/norma)
salf=(zz-zno)/norma
END IF
z3=(zz-zno)∗calf-(yy-yv)∗salf
kid=kid+1
vzn(kid)=zno
vsalfa(kid)=salf
vcalfa(kid)=calf
nid(ica,jca,kca)=kid
vdel(kid)=delta(xx,zno,calf,salf)
vk1(kid)=kappa(xx,zno,calf,salf)
vk2(kid)=pi/(pi-vdel(kid))∗(5.+wp)/4.
END IF
END IF
REPEAT LOOP
END IF
REPEAT LOOP
REPEAT LOOP
WRITE 'punti da correggere',kid
END catalogo
REAL FUNCTION derc(REAL fc,fs,fd,ds,dd,den)
RESULT=(fd∗ds**2-fs∗dd**2+fc∗(dd**2-ds**2))/den/h
END derc
REAL FUNCTION ders(REAL f0,f1,f2,d1,d2,den)
RESULT=(f1∗d2**2-f2∗d1**2-f0∗(d2**2-d1**2))/den/h
END ders
REAL FUNCTION des(REAL fl,fc,fr,frr,le,ri,rr)
REAL den,alr,alrr
alr=le+ri
alrr=alr+rr
IF le< pfo AND rr>0 THEN
den =alr∗alrr∗rr
RESULT=(fr∗alrr**2-frr∗alr**2-fl∗(alrr**2-alr**2))/h/den
ELSE
den=le∗ri∗alr
RESULT=(fc∗alr**2-fr∗le**2-fl∗(alr**2-le**2))/h/den
END IF
END des
REAL FUNCTION teta(REAL xx,yy)
REAL tet,rq
INTEGER si
IF xx = 0 THEN si=0
IF xx < 0 THEN si=-1
IF xx > 0 THEN si=1
rq=xx∗xx+yy∗yy
IF yy = 0 THEN
tet=pi/2∗si
ELSE
tet=ATAN(xx/yy)
IF yy < 0 THEN tet=tet+pi∗si
END IF
RESULT=tet
END teta
SUBROUTINE fun1(REAL fi1^,m11^,m21^,xx,yy)
REAL tet,e,a,d,b,fil,fis,k3,es,cot,s,c,r
r=(xx∗xx+yy∗yy)**0.5
tet=teta(xx,yy)
IF del > 0 THEN
s=SIN(2∗del)
c=COS(2∗del)
k3=(pi-del)∗k1
IF ABS(pi/2-k3) < 0.005 THEN
cot=0
ELSE
cot=1/TAN(k3)
END IF
IF wp = 1 THEN
e=1-2∗cot/(cot+TAN(del))
d=cot/(cot∗(1+e∗(1+k1)+k1∗c)-s∗k1)
b=1
a=4∗d-1
ELSE
e=2.∗TAN(del)∗cot/(1.-cot∗TAN(del))+1.
d=cot/(cot∗(1-k1∗c+e∗(1+k1))+k1∗s)
a=1
b=1-4∗d
END IF
es=1.+k1
IF (wp = 1) THEN
fis=r**es∗COS((1.-k1)∗tet)
fil=r**es∗COS(es∗tet)
fi1=d∗(fis+e∗fil)
m11=a∗r**k1∗SIN(k1∗tet)
m21=b∗r**k1∗COS(k1∗tet)
ELSE
fis=r**es∗SIN((1.-k1)∗tet)
fil=r**es∗SIN(es∗tet)
fi1=d∗(fis+e∗fil)
m11=a∗r**k1∗COS(k1∗tet)
m21=b∗r**k1∗SIN(k1∗tet)
END IF
ELSE
IF wp = -1 THEN
m11=r**0.5∗COS(tet/2.)
m21=0.
fis=r**1.5∗SIN(tet∗0.5)
fil=r**1.5∗SIN(tet∗1.5)
fi1=0.25∗(fis+fil)
ELSE
m11=-r**0.5∗SIN(0.5∗tet)/2
m21=r**0.5∗COS(0.5∗tet)
fis=r**1.5∗COS(tet∗0.5)
fil=r**1.5∗COS(tet∗1.5)
fi1=(fis+fil∗3)/8
END IF
END IF
END fun1
SUBROUTINE fun2(REAL fi2^,m12^,m22^,m32^,xx,yy)
REAL tet,es,k23,r
k23=pi/2./(pi-del)
r=(xx∗xx+yy∗yy)**0.5
tet=teta(xx,yy)
IF del > 0 THEN
es=k2-1.
IF (wp = 1) THEN
fi2=-r**k2∗COS(k2∗tet)/k2
m12=r**es∗SIN(es∗tet)
m22=-r**es∗COS(es∗tet)
m32=r**k23∗COS(k23∗tet)
ELSE
fi2=r**k2∗SIN(k2∗tet)/k2
m12=r**es∗COS(es∗tet)
m22=r**es∗SIN(es∗tet)
m32=r**k2∗SIN(k2∗tet)
END IF
ELSE
IF wp = -1 THEN
fi2=xx
m12=1.
m22=xx
m32=xx
ELSE
fi2=-2./3.∗r**1.5∗COS(1.5∗tet)
m12=r**0.5∗SIN(0.5∗tet)
m22=-r**0.5∗COS(0.5∗tet)
m32=r**0.5∗COS(0.5∗tet)
END IF
END IF
END fun2
SUBROUTINE costanti( INTEGER icor,jcor,kcor)
INTEGER ne,ive,kve,ke,n6,k6
neq=120
REAL xv,yv,den1,den2,denc6,dc6
REAL m1a,m2a,m3a,den,xc,yc,zc,lun,er
REAL fia,fib,fic,m11,m12,m21,m22,m32,fi1,fi2
REAL b1,b2,b3,ai(1..neq,1..2),a6(1..neq),bi(1..neq)
REAL b6(1..neq),bl(1..neq),ail(1..neq),a11,a12,a21,a22
REAL f1l,f2l,fl,fip,fip1,fip2,xp,yp,erv
IF icor∗h < 0.25 THEN
xv=0
ive=0
ELSE
xv=0.5
ive=nx
END IF
yv=fondo(xv,zn)
jve=btm(ive,kcor)
LOOP FOR iv=1 TO neq
bi(iv)=0.
LOOP FOR jv=1 TO 2
ai(iv,jv)=0.
REPEAT LOOP
REPEAT LOOP
LOOP FOR iv=1 TO neq
a6(iv)=0.
ail(iv)=0.
b6(iv)=0.
bl(iv)=0.
REPEAT LOOP
ke=0
k6=0
LOOP FOR iv=-1 TO 1
LOOP FOR kv=-1 TO 1
LOOP FOR jv=-1 TO 1
IF icor+iv > -1 AND icor+iv < nx+1 THEN
IF kcor+kv < nz+1 AND kcor+kv > -1 THEN
IF jcor+jv > btm(icor+iv,kcor+kv) -1 THEN
xc=(icor+iv)∗h-xv
yc=((jcor+jv)∗h-yv)∗calfa+(h∗kcor+kv∗h-zn)∗salfa
zc=-((jcor+jv)∗h-yv)∗salfa+(h∗kcor+kv∗h-zn)∗calfa
lun=1
m1a=m1(icor+iv,jcor+jv,kcor+kv)
m2a=m2(icor+iv,jcor+jv,kcor+kv)
m3a=m3(icor+iv,jcor+jv,kcor+kv)
fic=fi(icor+iv,jcor+jv,kcor+kv)
fun1(fi1,m11,m21,xc,yc)
fun2(fi2,m12,m22,m32,xc,yc)
ke=ke+1
k6=k6+1
ai(ke,1)=m11
ai(ke,2)=m12
bi(ke)=m1a
ail(k6)=m22
a6(k6)=m32
b6(k6)=m3a∗calfa-m2a∗salfa
bl(k6)=m2a∗calfa+m3a∗salfa
ke=ke+1
ai(ke,1)=fi1/lun
ai(ke,2)=fi2/lun
bi(ke)=fic/lun
ke=ke+1
ai(ke,1)=m21
ai(ke,2)=m22
bi(ke)=m2a∗calfa+m3a∗salfa
END IF
END IF
END IF
REPEAT LOOP
REPEAT LOOP
REPEAT LOOP
ne=ke
n6=k6
a11=0.
a12=0.
a21=0.
a22=0.
b1=0.
b2=0.
c6=0.
denc6=0.
dc6=0.
LOOP FOR ke=1 TO n6
denc6=denc6+a6(ke)**2
dc6=dc6+a6(ke)∗b6(ke)
REPEAT LOOP
c6=dc6/denc6
denc6=0.
dc6=0.
LOOP FOR ke=1 TO n6
denc6=denc6+ail(ke)**2
dc6=dc6+ail(ke)∗bl(ke)
REPEAT LOOP
c2l=dc6/denc6
LOOP FOR iv=1 TO ne
a11=ai(iv,1)∗ai(iv,1)+a11
a12=ai(iv,2)∗ai(iv,1)+a12
a21=ai(iv,1)∗ai(iv,2)+a21
a22=ai(iv,2)∗ai(iv,2)+a22
b1=bi(iv)∗ai(iv,1)+b1
b2=bi(iv)∗ai(iv,2)+b2
REPEAT LOOP
den=a11∗a22-a12∗a21
den1=b1∗a22-b2∗a12
den2=a11∗b2-b1∗a21
c2=den1/den
c4=den2/den
END costanti
REAL FUNCTION m1loc(REAL xx,yy)
REAL fi1,fi2,m11,m12,m21,m22,m32
fun1(fi1,m11,m21,xx,yy)
fun2(fi2,m12,m22,m32,xx,yy)
RESULT=c2∗m11+c4∗m12
END m1loc
REAL FUNCTION m2loc(REAL xx,yy)
REAL fi1,fi2,m11,m12,m21,m22,m32
fun1(fi1,m11,m21,xx,yy)
fun2(fi2,m12,m22,m32,xx,yy)
RESULT=c2∗m21+c4∗m22
END m2loc
REAL FUNCTION m3loc(REAL xx,yy)
REAL fi2,m12,m22,m32
fun2(fi2,m12,m22,m32,xx,yy)
RESULT=c6∗m32
END m3loc
REAL FUNCTION filoc(REAL xx,yy)
REAL fi1,fi2,m11,m12,m21,m22,m32
fun1(fi1,m11,m21,xx,yy)
fun2(fi2,m12,m22,m32,xx,yy)
RESULT=c2∗fi1+c4∗fi2
END filoc
REAL FUNCTION fico(INTEGER ic,jc,kc)
REAL xx,yy,xv
IF ic∗h < 0.25 THEN
xv=0
ELSE
xv=0.5
END IF
xx=ic∗h-xv
yy=(jc∗h-fondo(xv,zn))∗calfa+(kc∗h-zn)∗salfa
RESULT=filoc(xx,yy)
END fico
REAL FUNCTION m1co(INTEGER ic,jc,kc)
REAL xx,yy,xv
IF ic∗h < 0.25 THEN
xv=0
ELSE
xv=0.5
END IF
xx=ic∗h-xv
yy=(jc∗h-fondo(xv,zn))∗calfa+(kc∗h-zn)∗salfa
RESULT=m1loc(xx,yy)
END m1co
REAL FUNCTION m2co(INTEGER ic,jc,kc)
REAL xx,yy,xv
IF ic∗h < 0.25 THEN
xv=0
ELSE
xv=0.5
END IF
xx=ic∗h-xv
yy=(jc∗h-fondo(xv,zn))∗calfa+(kc∗h-zn)∗salfa
RESULT=m2loc(xx,yy)∗calfa-m3loc(xx,yy)∗salfa
END m2co
REAL FUNCTION m3co( INTEGER ic,jc,kc)
REAL xx,yy,xv
IF ic∗h < 0.25 THEN
xv=0
ELSE
xv=0.5
END IF
xx=ic∗h-xv
yy=(jc∗h-fondo(xv,zn))∗calfa+(kc∗h-zn)∗salfa
RESULT=m2loc(xx,yy)∗salfa+m3loc(xx,yy)∗calfa
END m3co
REAL FUNCTION fxd(REAL tang; INTEGER i,j,k,tra,co)
INTEGER ile,ifu,kfu
REAL fr,fd,fu,fuu,ld
IF tang>0 THEN
ile=1
ELSE
ile=-1
END IF
ifu=i+tra∗ile
kfu=k+(1-tra)∗ile
ld=ad(ifu,j,kfu)-ad(i,j,k)
IF co=0 THEN
fr=fi(ifu,j,kfu)
IF btm(ifu,kfu)> j THEN fr=0
fd=fi(ifu,j-1,kfu)
IF btm(ifu,kfu)> j-1 THEN fd=0
fu=fi(ifu,j+1,kfu)
fuu=fi(ifu,j+2,kfu)
ELSE
fr=fico(ifu,j,kfu)
IF btm(ifu,kfu)> j THEN fr=0
fd=fico(ifu,j-1,kfu)
IF btm(ifu,kfu)> j-1 THEN fd=0
fu=fico(ifu,j+1,kfu)
fuu=fico(ifu,j+2,kfu)
END IF
RESULT =(fd/h+des(fd,fr,fu,fuu,ad(ifu,j,kfu),1,1)∗ld)∗ile
END fxd
REAL FUNCTION fyl( INTEGER i,j,k,tra,co,leftside)
INTEGER ile,ifu,kfu
REAL fl,fc,fr,frr,ld,lrd
ile=leftside
IF ile=1 THEN
ld=al(i,j+1,k)
lrd=ld-al(i,j,k)
ELSE
ld=ar(i,j+1,k)
lrd=ld-ar(i,j,k)
END IF
ifu=i-tra∗ile
kfu=k-(1-tra)∗ile
IF co=0 THEN
fr=fi(i+tra∗ile,j+1,k+(1-tra)∗ile)
IF btm(i+tra∗ile,k+(1-tra∗ile))> j+1 THEN fr=0
frr=fi(i+2∗tra∗ile,j+1,k+2∗(1-tra)∗ile)
IF btm(i+2∗tra∗ile,k+2∗(1-tra∗ile))> j+1 THEN frr=0
fl=fi(ifu,j+1,kfu)
IF btm(ifu,kfu)> j+1 THEN fl=0
fc=fi(i,j+1,k)
ELSE
fr=fico(i+tra∗ile,j+1,k+(1-tra)∗ile)
IF btm(i+tra∗ile,k+(1-tra∗ile))> j+1 THEN fr=0
frr=fico(i+2∗tra∗ile,j+1,k+2∗(1-tra)∗ile)
IF btm(i+2∗tra∗ile,k+2∗(1-tra∗ile))> j+1 THEN frr=0
fl=fico(ifu,j+1,kfu)
IF btm(ifu,kfu)> j+1 THEN fl=0
fc=fico(i,j+1,k)
END IF
RESULT =fl/h+des(fl,fc,fr,frr,ld,1,1)∗lrd
END fyl
SUBROUTINE notom1
REAL den,cl,cr,cd,cf,cb,ml,mr,md,mf,mb,fix,alr,rr,ll,
tang,ficen,fl,fll,fr,frr,fb,fff,fd,fu,flad
LOOP FOR ino=0 TO nx
LOOP FOR kno=0 TO nz
LOOP FOR jno=btm(ino,kno) TO btf(ino,kno)-1
alr=al(ino,jno,kno)+ar(ino,jno,kno)
den=al(ino,jno,kno)∗ar(ino,jno,kno)∗alr
cl=0.
cd=0.
cr=0.
cb=0.
cf=0.
mr=0.
md=0.
ml=0.
mf=0.
mb=0.
fu=fi(ino,jno+1,kno)
ficen=fi(ino,jno,kno)
IF btm(ino-1,kno) > jno THEN
fl=0
rr=0
ELSE
fl=fi(ino-1,jno,kno)
rr=al(ino-1,jno,kno)
END IF
IF ino=0 THEN
fll=wp∗fi(ino+2,jno,kno)
IF btm(ino+2,kno)> jno THEN fll=0
ELSE
fll=fi(ino-2,jno,kno)
IF btm(ino-2,kno)> jno THEN fll=0
END IF
IF btm(ino+1,kno) > jno THEN
fr=0
ll=0
ELSE
fr=fi(ino+1,jno,kno)
ll=ar(ino+1,jno,kno)
END IF
IF ino=nx THEN
frr=wp∗fi(ino-2,jno,kno)
IF btm(ino-2,kno)> jno THEN frr=0
ELSE
frr=fi(ino+2,jno,kno)
IF btm(ino+2,kno)> jno THEN frr=0
END IF
IF btm(ino,kno-1) > jno THEN
fb=0
ELSE
fb=fi(ino,jno,kno-1)
END IF
IF btm(ino,kno+1) > jno THEN
fff=0
ELSE
fff=fi(ino,jno,kno+1)
END IF
IF al(ino,jno,kno) < pfo THEN
fix=ders(ficen,fr,frr,1,2,2)
ELSE IF ar(ino,jno,kno) < pfo THEN
fix=-ders(ficen,fl,fll,1,2,2)
ELSE
fix=derc(ficen,fl,fr,al(ino,jno,kno),ar(ino,jno,kno),den)
END IF
IF btm(ino-1,kno) > jno THEN
cl=bcl(ino,jno,kno)
IF al(ino,jno,kno) < pfo THEN
ml=ders(0,fr,frr,alr,alr+1,alr∗(alr+1))
ELSE
ml=ders(0,ficen,fr,al(ino,jno,kno),alr,den)
END IF
END IF
IF btm(ino+1,kno) > jno THEN
cr=bcr(ino,jno,kno)
IF ar(ino,jno,kno) < pfo THEN
mr=-ders(0,fl,fll,alr,alr+1,alr∗(alr+1))
ELSE
mr=-ders(0,ficen,fl,ar(ino,jno,kno),alr,den)
END IF
END IF
IF jno = btm(ino,kno) THEN
cd=bcd(ino,jno,kno)
tang=(fondo(ino∗h-eps,kno∗h)-fondo(ino∗h+eps,h∗kno))/2/eps
IF pbc=1 THEN
md=fix
ELSE IF pbc=2 THEN
IF ino=0 OR ino= nx THEN
md=fix
ELSE
IF tang=0 THEN
md=0
ELSE
md=fxd(tang,ino,jno,kno,1,0)
END IF
END IF
END IF
END IF
IF btm(ino,kno-1) > jno THEN
cb=bcb(ino,jno,kno)
mb=fix
END IF
IF btm(ino,kno+1) > jno THEN
cf=bcf(ino,jno,kno)
mf=fix
END IF
rhs(ino,jno,kno)=cl∗ml+cr∗mr+cd∗md+cf∗mf+cb∗mb
REPEAT LOOP
LOOP FOR jno=btf(ino,kno) TO ny-1
rhs(ino,jno,kno)=0.
REPEAT LOOP
rhs(ino,ny,kno)=h/3∗(1-wp)/2
REPEAT LOOP
REPEAT LOOP
END notom1
SUBROUTINE notom2
REAL cl,cr,cf,cd,cb,ml,mr,md,mb,mf,den,fiy,fd,ad1,tang
REAL fir,ficen,fu,fl,fr,fb,fff
LOOP FOR ino=0 TO nx
LOOP FOR kno=0 TO nz
LOOP FOR jno=btm(ino,kno) TO btf(ino,kno)-1
cl=0.
cd=0.
cr=0.
cf=0.
cb=0.
ml=0.
md=0.
mr=0.
mf=0.
mb=0.
fu=fi(ino,jno+1,kno)
ficen=fi(ino,jno,kno)
IF btm(ino-1,kno) > jno THEN
fl=0
ELSE
fl=fi(ino-1,jno,kno)
END IF
IF btm(ino+1,kno) > jno THEN
fr=0
ELSE
fr=fi(ino+1,jno,kno)
END IF
IF btm(ino,kno-1) > jno THEN
fb=0
ELSE
fb=fi(ino,jno,kno-1)
END IF
IF btm(ino,kno+1) > jno THEN
fff=0
ELSE
fff=fi(ino,jno,kno+1)
END IF
ad1=ad(ino,jno,kno)+1.
den=ad(ino,jno,kno)∗ad1
IF ad(ino,jno,kno) < pfo THEN
fiy=ders(ficen,fu,fi(ino,jno+2,kno),1,2,2)
ELSE
fd=fi(ino,jno-1,kno)
IF jno = btm(ino,kno) THEN fd=0.
fiy=derc(ficen,fd,fu,ad(ino,jno,kno),1,den)
END IF
tang=(fondo(ino∗h-eps,kno∗h)-fondo(ino∗h+eps,h∗kno))/2/eps
IF btm(ino-1,kno) > jno THEN
cl=bcl(ino,jno,kno)
IF pbc=1 THEN
ml =fiy
ELSE IF pbc=2 THEN
IF ino =1 AND jno = btm(ino-1,kno) THEN
ml=0
ELSE
ml =fyl(ino,jno,kno,1,0,1)
END IF
END IF
END IF
IF btm(ino+1,kno) > jno THEN
cr=bcr(ino,jno,kno)
IF pbc=1 THEN
mr =fiy
ELSE IF pbc=2 THEN
IF ino =nx-1 AND jno = btm(ino+1,kno) THEN
mr=0
ELSE
mr =fyl(ino,jno,kno,1,0,-1)
END IF
END IF
END IF
IF jno = btm(ino,kno) THEN
cd=bcd(ino,jno,kno)
IF( ino=0 OR ino=nx) AND jno >1 THEN
md=0
ELSE
IF ad(ino,jno,kno) < pfo THEN
md=ders(0,fu,fi(ino,jno+2,kno),ad1,ad1+1,ad1∗(ad1+1))
ELSE
md=ders(0,ficen,fu,ad(ino,jno,kno),ad1,den)
END IF
END IF
END IF
IF btm(ino,kno-1) > jno THEN
cb=bcb(ino,jno,kno)
IF pbc=1 THEN
mb=fiy
ELSE IF pbc=2 THEN
mb =fyl(ino,jno,kno,0,0,1)
END IF
END IF
IF btm(ino,kno+1) > jno THEN
cf=bcf(ino,jno,kno)
IF pbc=1 THEN
mf=fiy
ELSE IF pbc=2 THEN
mf =fyl(ino,jno,kno,0,0,-1)
END IF
END IF
rhs(ino,jno,kno)=cl∗ml+cr∗mr+cd∗md+mb∗cb+mf∗cf
REPEAT LOOP
LOOP FOR jno=btf(ino,kno) TO ny-1
rhs(ino,jno,kno)=0.
REPEAT LOOP
rhs(ino,ny,kno)=0
REPEAT LOOP
REPEAT LOOP
END notom2
SUBROUTINE notom3
REAL den,fiz,cf,cb,cl,cr,cd,ml,mr,md,mf,mb,abf,fiq,fit,
tangz,fu,ficen,fl,fr,fb,fff,ffff,fbb
LOOP FOR ino=0 TO nx
LOOP FOR kno=0 TO nz
LOOP FOR jno=btm(ino,kno) TO btf(ino,kno)-1
abf=ab(ino,jno,kno)+af(ino,jno,kno)
den=ab(ino,jno,kno)∗af(ino,jno,kno)∗abf
cb=0.
cf=0.
cl=0.
cr=0.
mf=0.
mb=0.
cd=0.
md=0.
ml=0.
mr=0.
fu=fi(ino,jno+1,kno)
ficen=fi(ino,jno,kno)
IF btm(ino-1,kno) > jno THEN
fl=0
ELSE
fl=fi(ino-1,jno,kno)
END IF
IF btm(ino+1,kno) > jno THEN
fr=0
ELSE
fr=fi(ino+1,jno,kno)
END IF
IF btm(ino,kno-1) > jno THEN
fb=0
ELSE
fb=fi(ino,jno,kno-1)
END IF
IF btm(ino,kno+1) > jno THEN
fff=0
ELSE
fff=fi(ino,jno,kno+1)
END IF
ad1=ad(ino,jno,kno)+1.
IF ab(ino,jno,kno) < pfo AND kno < nz THEN
fiz=ders(ficen,fff,fi(ino,jno,kno+2),1,2,2)
ELSE IF af(ino,jno,kno) < pfo AND kno > 0 THEN
fiz=-ders(ficen,fb,fi(ino,jno,kno-2),1,2,2)
ELSE
fiz=derc(ficen,fb,fff,ab(ino,jno,kno),af(ino,jno,kno),den)
END IF
IF (ab(ino,jno,kno) < pfo AND kno=nz) OR (af(ino,jno,kno) < pfo AND kno =0) THEN \
fiz=(fff- fb)/abf/h
IF btm(ino,kno-1) > jno THEN
cb=bcb(ino,jno,kno)
IF ab(ino,jno,kno) < pfo AND kno < nz THEN
mb=ders(0,fff,fi(ino,jno,kno+2),abf,abf+1,abf∗(abf+1))
ELSE
mb=ders(0,ficen,fff,ab(ino,jno,kno),abf,den)
END IF
IF ab(ino,jno,kno) < pfo AND kno=nz THEN mb=fiz
IF ino=0 OR ino=nx THEN mb=0
END IF
IF btm(ino,kno+1) > jno THEN
cf=bcf(ino,jno,kno)
IF af(ino,jno,kno) < pfo AND kno > 0 THEN
mf=-ders(0,fb,fi(ino,jno,kno-2),abf,abf+1,abf∗(abf+1))
ELSE
mf=-ders(0,ficen,fb,af(ino,jno,kno),abf,den)
END IF
IF ino=0 OR ino=nx THEN mf=0
END IF
IF af(ino,jno,kno) < pfo AND kno=0 THEN mf=fiz
IF jno = btm(ino,kno) THEN
cd=bcd(ino,jno,kno)
tangz=(fondo(ino∗h,kno∗h-eps)-fondo(ino∗h,h∗kno+eps))/2/eps
IF pbc =1 THEN
md=fiz
ELSE IF pbc =2 THEN
md=fxd(tangz,ino,jno,kno,0,0)
IF ino=0 OR ino=nx THEN md=0
END IF
IF tangz=0 THEN md=0
IF kno=0 AND wp=-1 THEN md=0
END IF
IF btm(ino-1,kno) > jno THEN
cl=bcl(ino,jno,kno)
IF pbc =1 THEN
ml=fiz
ELSE IF pbc =2 THEN
ml=0
END IF
END IF
IF btm(ino+1,kno) > jno THEN
cr=bcr(ino,jno,kno)
IF pbc =1 THEN
mr=fiz
ELSE IF pbc =2 THEN
mr=0
END IF
END IF
rhs(ino,jno,kno)=cd∗md+cb∗mb+cf∗mf+cl∗ml+cr∗mr
REPEAT LOOP
LOOP FOR jno=btf(ino,kno) TO ny-1
rhs(ino,jno,kno)=0.
REPEAT LOOP
rhs(ino,ny,kno)=h/3∗(wp+1)/2
REPEAT LOOP
REPEAT LOOP
END notom3
SUBROUTINE notofia
REAL den1,den2,den3,m1x,m2y,m3z,alr,abf,fff,fb,md,ml,mr,mb,mf,ad1
REAL ficen,fu,fl,fr,fll,frr
LOOP FOR ino=0 TO nx
LOOP FOR kno=0 TO nz
LOOP FOR jno=btm(ino,kno) TO btf(ino,kno)-1
alr=al(ino,jno,kno)+ar(ino,jno,kno)
abf=ab(ino,jno,kno)+af(ino,jno,kno)
ad1=ad(ino,jno,kno)+1.
den1=al(ino,jno,kno)∗ar(ino,jno,kno)∗alr
den2=ad(ino,jno,kno)∗ad1
den3=af(ino,jno,kno)∗ab(ino,jno,kno)∗abf
ficen=fi(ino,jno,kno)
fu=fi(ino,jno+1,kno)
ml=m1(ino-1,jno,kno)
mr=m1(ino+1,jno,kno)
md=m2(ino,jno-1,kno)
mb=m3(ino,jno,kno-1)
mf=m3(ino,jno,kno+1)
fb=fi(ino,jno,kno-1)
fff=fi(ino,jno,kno+1)
IF btm(ino,kno-1) > jno THEN fb=0.
IF btm(ino,kno+1) > jno THEN fff=0.
IF btm(ino-1,kno) > jno THEN
fl=0
ELSE
fl=fi(ino-1,jno,kno)
END IF
IF ino=0 THEN
fll=wp∗fi(ino+2,jno,kno)
IF btm(ino+2,kno)> jno THEN fll=0
ELSE
fll=fi(ino-2,jno,kno)
IF btm(ino-2,kno)> jno THEN fll=0
END IF
IF btm(ino+1,kno) > jno THEN
fr=0
ELSE
fr=fi(ino+1,jno,kno)
END IF
IF ino=nx THEN
frr=wp∗fi(ino-2,jno,kno)
IF btm(ino-2,kno)> jno THEN frr=0
ELSE
frr=fi(ino+2,jno,kno)
IF btm(ino+2,kno)> jno THEN frr=0
END IF
IF jno = btm(ino,kno) THEN
IF ad(ino,jno,kno) < pfo THEN
md=ders(0,fu,fi(ino,jno+2,kno),ad1,ad1+1,ad1∗(ad1+1))
ELSE
md=ders(0,ficen,fu,ad(ino,jno,kno),ad1,den2)
END IF
IF (ino=0 OR ino= nx) AND jno >1 THEN md=0
END IF
IF btm(ino+1,kno) > jno THEN
IF ar(ino,jno,kno) < pfo THEN
mr=-ders(0,fl,fll,alr,alr+1,alr∗(alr+1))
ELSE
mr=-ders(0,ficen,fl,ar(ino,jno,kno),alr,den1)
END IF
END IF
IF btm(ino-1,kno) > jno THEN
IF al(ino,jno,kno) < pfo THEN
ml=ders(0,fr,frr,alr,alr+1,alr∗(alr+1))
ELSE
ml=ders(0,ficen,fr,al(ino,jno,kno),alr,den1)
END IF
END IF
IF btm(ino,kno+1) > jno THEN
IF af(ino,jno,kno) < pfo AND kno > 0 THEN
mf=-ders(0,fi(ino,jno,kno-1),fi(ino,jno,kno-2),abf,abf+1,abf∗(abf+1))
ELSE
mf=-ders(0,fi(ino,jno,kno),fi(ino,jno,kno-1),af(ino,jno,kno),abf,den3)
IF af(ino,jno,kno) < pfo AND kno=0 THEN mf=(fff-fb)/abf/h
END IF
IF ino=0 OR ino= nx THEN mf=0
END IF
IF btm(ino,kno-1) > jno THEN
IF ab(ino,jno,kno) < pfo AND kno < nz THEN
mb=ders(0,fi(ino,jno,kno+1),fi(ino,jno,kno+2),abf,abf+1,abf∗(abf+1))
ELSE
mb=ders(0,fi(ino,jno,kno),fi(ino,jno,kno+1),ab(ino,jno,kno),abf,den3)
IF ab(ino,jno,kno) < pfo AND kno =nz THEN mb=(fff-fb)/abf/h
END IF
IF ino=0 OR ino= nx THEN mb=0
END IF
IF al(ino,jno,kno) < pfo AND ino < nx THEN
m1x=ders(m1(ino,jno,kno),mr,m1(ino+2,jno,kno),1,2,2)
ELSE IF ar(ino,jno,kno) < pfo AND ino > 0 THEN
m1x=-ders(m1(ino,jno,kno),ml,m1(ino-2,jno,kno),1,2,2)
ELSE
m1x=derc(m1(ino,jno,kno),ml,mr,al(ino,jno,kno),ar(ino,jno,kno),den1)
END IF
IF ad(ino,jno,kno) < pfo THEN
m2y=ders(m2(ino,jno,kno),m2(ino,jno+1,kno),m2(ino,jno+2,kno),1,2,2)
ELSE
m2y=derc(m2(ino,jno,kno),md,m2(ino,jno+1,kno),ad(ino,jno,kno),1,den2)
END IF
IF ab(ino,jno,kno) < pfo AND kno < nz THEN
m3z=ders(m3(ino,jno,kno),mf,m3(ino,jno,kno+2),1,2,2)
ELSE IF af(ino,jno,kno) < pfo AND kno > 0 THEN
m3z=-ders(m3(ino,jno,kno),mb,m3(ino,jno,kno-2),1,2,2)
ELSE
m3z=derc(m3(ino,jno,kno),mb,mf,ab(ino,jno,kno),af(ino,jno,kno),den3)
IF(af(ino,jno,kno) < pfo AND kno=0) OR (ab(ino,jno,kno) < pfo AND kno=nz ) THEN m3z=(mf-mb)/abf/h
END IF
rhs(ino,jno,kno)=-(m1x+m2y+m3z)∗coenoto(ino,jno,kno)∗h∗h
REPEAT LOOP
LOOP FOR jno=btf(ino,kno) TO ny-1
rhs(ino,jno,kno)=-(m1(ino+1,jno,kno)-m1(ino-1,jno,kno)+m2(ino,jno+1,kno)-m2(ino,jno-1,kno)+m3(ino,jno,kno+1)-m3(ino,jno,kno-1))∗h/12
REPEAT LOOP
rhs(ino,ny,kno)=-(m1(ino+1,ny,kno)-m1(ino-1,ny,kno)-2.∗m2(ino,ny-1,kno)+m3(ino,ny,kno+1)-m3(ino,ny,kno-1))∗h/12
REPEAT LOOP
REPEAT LOOP
END notofia
SUBROUTINE corm1
REAL ficen,ml,mr,md,mf,mb,fix,alr,rhs1,cen,lef,rig,u,d,ba,fo,xv
REAL den,cl,cr,cd,cf,cb,fl,fll,fr,fb,fff,fu,fd,frr,tang,ll,rr
INTEGER ndis
LOOP FOR k=0 TO nz
LOOP FOR i=0 TO nx
IF i < 2 OR i > nx-2 THEN
IF i∗h < 0.25 THEN
xv=0.
ELSE
xv=0.5
END IF
IF fondo(xv,k∗h)>2∗h THEN
LOOP FOR j=btm(i,k) TO btm(0,0)+2
ndis=nid(i,j,k)
IF ndis > 0 THEN
del=vdel(ndis)
k1=vk1(ndis)
k2=vk2(ndis)
zn=vzn(ndis)
calfa=vcalfa(ndis)
salfa=vsalfa(ndis)
costanti(i,j,k)
alr=al(i,j,k)+ar(i,j,k)
den=al(i,j,k)∗ar(i,j,k)∗alr
u=m1co(i,j+1,k)
cen=m1co(i,j,k)
ba=0.
fo=0.
lef=0.
rig=0.
mr=0.
md=0.
ml=0.
mf=0.
mb=0.
cr=0.
cl=0.
cd=0.
cb=0.
cf=0.
d=0
ficen=fico(i,j,k)
fu=fico(i,j+1,k)
IF btm(i-1,k) > j THEN
fl=0
rr=0
ELSE
fl=fico(i-1,j,k)
rr=al(i-1,j,k)
END IF
IF i=0 THEN
fll=wp∗fico(i+2,j,k)
IF btm(i+2,k)> j THEN fll=0
ELSE
fll=fico(i-2,j,k)
IF btm(i-2,k)> j THEN fll=0
END IF
IF btm(i+1,k) > j THEN
fr=0
ll=0
ELSE
fr=fico(i+1,j,k)
ll=al(i+1,j,k)
END IF
IF i=nx THEN
frr=wp∗fico(i-2,j,k)
IF btm(i-2,k)> j THEN frr=0
ELSE
frr=fico(i+2,j,k)
IF btm(i+2,k)> j THEN frr=0
END IF
IF btm(i,k-1) > j THEN
fb=0
ELSE
fb=fico(i,j,k-1)
END IF
IF btm(i,k+1) > j THEN
fff=0
ELSE
fff=fico(i,j,k+1)
END IF
IF btm(i,k)=j THEN
fd=0
ELSE
fd=fico(i,j-1,k)
END IF
IF al(i,j,k) < pfo THEN
fix=ders(ficen,fr,frr,1,2,2)
ELSE IF ar(i,j,k) < pfo THEN
fix=-ders(ficen,fl,fll,1,2,2)
ELSE
fix=derc(ficen,fl,fr,al(i,j,k),ar(i,j,k),den)
END IF
IF btm(i-1,k) > j THEN
cl=bcl(i,j,k)
lef=0.
IF al(i,j,k) < pfo THEN
ml=ders(0,fr,frr,alr,alr+1,alr∗(alr+1))
ELSE
ml=ders(0,ficen,fr,al(i,j,k),alr,den)
END IF
ELSE
lef=m1co(i-1,j,k)
END IF
IF btm(i+1,k) > j THEN
cr=bcr(i,j,k)
rig=0.
IF ar(i,j,k) < pfo THEN
mr=-ders(0,fl,fll,alr,alr+1.,alr∗(alr+1.))
ELSE
mr=-ders(0,ficen,fl,ar(i,j,k),alr,den)
END IF
ELSE
rig=m1co(i+1,j,k)
END IF
IF j > btm(i,k) THEN
d=m1co(i,j-1,k)
ELSE
cd=bcd(i,j,k)
d=0.
tang=(fondo(i∗h-eps,k∗h)-fondo(i∗h+eps,h∗k))/2/eps
IF pbc=1 THEN
md=fix
ELSE IF pbc=2 THEN
IF i=0 OR i= nx THEN
md=fix
ELSE
IF tang=0 THEN
md=0
ELSE
md=fxd(tang,i,j,k,1,1)
END IF
END IF
END IF
END IF
IF btm(i,k-1) > j THEN
cb=bcb(i,j,k)
ba=0.
IF pbc =1 THEN
mb=fix
ELSE IF pbc =2 THEN
mb=fix
END IF
ELSE
ba=m1co(i,j,k-1)
END IF
IF btm(i,k+1) > j THEN
cf=bcf(i,j,k)
fo=0
IF pbc =1 THEN
mf=fix
ELSE IF pbc =2 THEN
mf=fix
END IF
ELSE
fo=m1co(i,j,k+1)
END IF
rhs1=cl∗ml+cd∗md+cb∗mb+cf∗mf+cr∗mr
rhs(i,j,k)=rhs(i,j,k)-left(i,j,k)∗lef-right(i,j,k)∗rig-
down(i,j,k)∗d-up(i,j,k)∗u-back(i,j,k)∗ba-forw(i,j,k)∗fo+cen-rhs1
END IF
REPEAT LOOP
END IF
END IF
REPEAT LOOP
REPEAT LOOP
END corm1
SUBROUTINE corm2
REAL ficen,ml,mr,md,mf,mb,fiy,fd,rhs1,cen,lef,rig,u,d,ba,fo,xv
REAL den,cl,cr,cd,cf,cb,tang,fiy0,fir,fu,fl,fr,fb,fff
INTEGER ndis
LOOP FOR k=0 TO nz
LOOP FOR i=0 TO nx
IF i < 2 OR i > nx-2 THEN
IF i∗h < 0.25 THEN
xv=0.
ELSE
xv=0.5
END IF
IF fondo(xv,k∗h)>2∗h THEN
LOOP FOR j=btm(i,k) TO btm(0,0)+2
ndis=nid(i,j,k)
IF ndis > 0 THEN
del=vdel(ndis)
k1=vk1(ndis)
k2=vk2(ndis)
zn=vzn(ndis)
calfa=vcalfa(ndis)
salfa=vsalfa(ndis)
costanti(i,j,k)
ad1=ad(i,j,k)+1.
den=ad(i,j,k)∗ad1
cl=0.
cd=0.
cr=0.
cf=0.
cb=0.
ml=0.
md=0.
mr=0.
mf=0.
mb=0.
ba=0
fo=0
lef=0
rig=0
d=0
u=m2co(i,j+1,k)
cen=m2co(i,j,k)
ficen=fico(i,j,k)
fu=fico(i,j+1,k)
IF btm(i-1,k) > j THEN
fl=0
ELSE
fl=fico(i-1,j,k)
END IF
IF btm(i+1,k) > j THEN
fr=0
ELSE
fr=fico(i+1,j,k)
END IF
IF btm(i,k-1) > j THEN
fb=0
ELSE
fb=fico(i,j,k-1)
END IF
IF btm(i,k+1) > j THEN
fff=0
ELSE
fff=fico(i,j,k+1)
END IF
IF btm(i,k)=j THEN
fd=0
ELSE
fd=fico(i,j-1,k)
END IF
IF ad(i,j,k) < pfo THEN
fiy=ders(ficen,fu,fico(i,j+2,k),1,2,2)
ELSE
fiy=derc(ficen,fd,fu,ad(i,j,k),1,den)
END IF
IF j=btm(i,k) THEN
IF ad(i,j,k) < pfo THEN
fiy0=ders(0,fu,fico(i,j+2,k),ad1,ad1+1,ad1∗(ad1+1))
ELSE
fiy0=ders(0,ficen,fu,ad(i,j,k),ad1,den)
END IF
END IF
tang=(fondo(i∗h-eps,k∗h)-fondo(i∗h+eps,h∗k))/2/eps
IF btm(i-1,k) > j THEN
cl=bcl(i,j,k)
lef=0.
IF pbc=1 THEN
ml=fiy
ELSE IF pbc=2 THEN
IF i=1 AND j= btm(i-1,k) THEN
ml=0
ELSE
ml=fyl(i,j,k,1,1,1)
END IF
END IF
ELSE
lef=m2co(i-1,j,k)
END IF
IF btm(i+1,k) > j THEN
cr=bcr(i,j,k)
rig=0.
IF pbc=1 THEN
mr=fiy
ELSE IF pbc=2 THEN
IF i=nx-1 AND j= btm(i+1,k) THEN
mr=0
ELSE
mr=fyl(i,j,k,1,1,-1)
END IF
END IF
ELSE
rig=m2co(i+1,j,k)
END IF
IF j > btm(i,k) THEN
d=m2co(i,j-1,k)
ELSE
cd=bcd(i,j,k)
d=0.
IF (i=0 OR i=nx) AND j >1 THEN
md=0
ELSE
md=fiy0
END IF
END IF
IF btm(i,k-1) > j THEN
cb=bcb(i,j,k)
ba=0
IF pbc=1 THEN
mb=fiy
ELSE IF pbc=2 THEN
mb=fyl(i,j,k,0,1,1)
END IF
ELSE
ba=m2co(i,j,k-1)
END IF
IF btm(i,k+1) > j THEN
cf=bcf(i,j,k)
fo=0
IF pbc=1 THEN
mf=fiy
ELSE IF pbc=2 THEN
mf=fyl(i,j,k,0,1,-1)
END IF
ELSE
fo=m2co(i,j,k+1)
END IF
rhs1=cl∗ml+cd∗md+cb∗mb+cf∗mf+cr∗mr
rhs(i,j,k)=rhs(i,j,k)-left(i,j,k)∗lef-right(i,j,k)∗rig-down(i,j,k) \
∗d-up(i,j,k)∗u-back(i,j,k)∗ba-forw(i,j,k)∗fo+cen-rhs1
END IF
REPEAT LOOP
END IF
END IF
REPEAT LOOP
REPEAT LOOP
END corm2
SUBROUTINE corm3
REAL ficen,ml,mr,md,mf,mb,fiz,rhs1,cen,lef,rig,u,d,ba,fo,abf,xv
REAL den,cl,cr,cd,cf,cb,fiq,fit,tangz,fl,fr,fb,fff,fu,ffff,fbb
INTEGER ndis
LOOP FOR k=0 TO nz
LOOP FOR i=0 TO nx
IF i < 2 OR i > nx-2 THEN
IF i∗h < 0.25 THEN
xv=0.
ELSE
xv=0.5
END IF
IF fondo(xv,h∗k)>2∗h THEN
LOOP FOR j=btm(i,k) TO btm(0,0)+2
ndis=nid(i,j,k)
IF ndis > 0 THEN
del=vdel(ndis)
k1=vk1(ndis)
k2=vk2(ndis)
zn=vzn(ndis)
calfa=vcalfa(ndis)
salfa=vsalfa(ndis)
costanti(i,j,k)
abf=ab(i,j,k)+af(i,j,k)
den=ab(i,j,k)∗af(i,j,k)∗abf
ad1=ad(i,j,k)+1
cb=0.
cf=0.
cl=0.
cr=0.
cd=0.
mf=0.
mb=0.
md=0.
ml=0.
mr=0.
u=m3co(i,j+1,k)
cen=m3co(i,j,k)
ficen=fico(i,j,k)
fu=fico(i,j+1,k)
IF btm(i-1,k) > j THEN
fl=0
ELSE
fl=fico(i-1,j,k)
END IF
IF btm(i+1,k) > j THEN
fr=0
ELSE
fr=fico(i+1,j,k)
END IF
IF btm(i,k-1) > j THEN
fb=0
ELSE
fb=fico(i,j,k-1)
END IF
IF btm(i,k+1) > j THEN
fff=0
ELSE
fff=fico(i,j,k+1)
END IF
IF ab(i,j,k) < pfo AND k < nz THEN
fiz=ders(ficen,fff,fico(i,j,k+2),1,2,2)
ELSE IF af(i,j,k) < pfo AND k > 0 THEN
fiz=-ders(ficen,fb,fico(i,j,k-2),1,2,2)
ELSE
fiz=derc(ficen,fb,fff,ab(i,j,k),af(i,j,k),den)
END IF
IF(ab(i,j,k) < pfo AND k=nz) OR (af(i,j,k)<pfo AND k=0) THEN fiz=(fff-fb)/abf/h
IF btm(i,k-1) > j THEN
ba=0.
cb=bcb(i,j,k)
IF ab(i,j,k) < pfo AND k < nz THEN
mb=ders(0,fff,fico(i,j,k+2),abf,abf+1,abf∗(abf+1))
ELSE
mb=ders(0,ficen,fff,ab(i,j,k),abf,den)
END IF
IF i=0 OR i=nx THEN mb=0
ELSE
ba=m3co(i,j,k-1)
END IF
IF btm(i,k+1) > j THEN
fo=0.
cf=bcf(i,j,k)
IF af(i,j,k) < pfo AND k > 0 THEN
mf=-ders(0,fb,fico(i,j,k-2),abf,abf+1,abf∗(abf+1))
ELSE
mf=-ders(0,ficen,fb,af(i,j,k),abf,den)
END IF
IF i=0 OR i=nx THEN mf=0
ELSE
fo=m3co(i,j,k+1)
END IF
IF j > btm(i,k) THEN
d=m3co(i,j-1,k)
ELSE
d=0
cd=bcd(i,j,k)
tangz=(fondo(i∗h,k∗h-eps)-fondo(i∗h,h∗k+eps))/2/eps
IF pbc=1 THEN
md=fiz
ELSE IF pbc=2 THEN
md=fxd(tangz,i,j,k,0,1)
IF i=0 OR i=nx THEN md=0
END IF
IF tangz=0 THEN md=0
IF k=0 AND wp=-1 THEN md=0
END IF
IF btm(i-1,k) > j THEN
cl=bcl(i,j,k)
lef=0
IF pbc=1 THEN
ml=fiz
ELSE IF pbc=2 THEN
ml=0
END IF
ELSE
lef=m3co(i-1,j,k)
END IF
IF btm(i+1,k) > j THEN
cr=bcr(i,j,k)
rig=0
IF pbc=1 THEN
mr=fiz
ELSE IF pbc=2 THEN
mr=0
END IF
ELSE
rig=m3co(i+1,j,k)
END IF
rhs1=cl∗ml+cd∗md+cb∗mb+cf∗mf+mr∗cr
rhs(i,j,k)=rhs(i,j,k)-left(i,j,k)∗lef-right(i,j,k)∗rig-down(i,j,k) \
∗d-up(i,j,k)∗u-back(i,j,k)∗ba-forw(i,j,k)∗fo+cen-rhs1
END IF
REPEAT LOOP
END IF
END IF
REPEAT LOOP
REPEAT LOOP
END corm3
SUBROUTINE corfi
REAL den1,den2,den3
REAL m1x,m2y,m3z,alr,abf,ml,mr,mb,mf,ad1,md,ficen,lef,rig,d,fu,ba,xv,fo
REAL fl,fr,fb,fff,fll,frr
INTEGER ndis
LOOP FOR k=0 TO nz
LOOP FOR i=0 TO nx
IF i < 2 OR i > nx-2 THEN
IF i∗h < 0.25 THEN
xv=0.
ELSE
xv=0.5
END IF
IF fondo(xv,h∗k)>2∗h THEN
LOOP FOR j=btm(i,k) TO btm(0,0)+2
ndis=nid(i,j,k)
IF ndis > 0 THEN
del=vdel(ndis)
k1=vk1(ndis)
k2=vk2(ndis)
zn=vzn(ndis)
calfa=vcalfa(ndis)
salfa=vsalfa(ndis)
costanti(i,j,k)
alr=al(i,j,k)+ar(i,j,k)
abf=ab(i,j,k)+af(i,j,k)
ad1=ad(i,j,k)+1.
den1=al(i,j,k)∗ar(i,j,k)∗alr
den2=ad(i,j,k)∗ad1
den3=af(i,j,k)∗ab(i,j,k)∗abf
mr=m1co(i+1,j,k)
rig=fico(i+1,j,k)
fu=fico(i,j+1,k)
fff=fico(i,j,k+1)
IF btm(i,k+1) > j THEN fff=0.
fb=fico(i,j,k-1)
IF btm(i,k-1) > j THEN fb=0.
ficen=fico(i,j,k)
IF btm(i-1,k) > j THEN
fl=0
ELSE
fl=fico(i-1,j,k)
END IF
IF i=0 THEN
fll=wp∗fico(i+2,j,k)
IF btm(i+2,k)> j THEN fll=0
ELSE
fll=fico(i-2,j,k)
IF btm(i-2,k)> j THEN fll=0
END IF
IF btm(i+1,k) > j THEN
fr=0
ELSE
fr=fico(i+1,j,k)
END IF
IF i=nx THEN
frr=wp∗fico(i-2,j,k)
IF btm(i-2,k)> j THEN frr=0
ELSE
frr=fico(i+2,j,k)
IF btm(i+2,k)> j THEN frr=0
END IF
IF j > btm(i,k) THEN
d=fico(i,j-1,k)
md=m2co(i,j-1,k)
ELSE
d=0.
IF ad(i,j,k) < pfo THEN
md=ders(0,fu,fico(i,j+2,k),ad1,ad1+1,ad1∗(ad1+1))
ELSE
md=ders(0,ficen,fu,ad(i,j,k),ad1,den2)
END IF
IF (i=0 OR i=nx) AND j>1 THEN md=0
END IF
IF btm(i+1,k) > j THEN
rig=0.
IF ar(i,j,k) < pfo THEN
mr=-ders(0,fl,fll,alr,alr+1,alr∗(alr+1))
ELSE
mr=-ders(0,ficen,fl,ar(i,j,k),alr,den1)
END IF
ELSE
rig=fico(i+1,j,k)
mr=m1co(i+1,j,k)
END IF
IF btm(i-1,k) > j THEN
lef=0.
IF al(i,j,k) < pfo THEN
ml=ders(0,fr,frr,alr,alr+1,alr∗(alr+1))
ELSE
ml=ders(0,ficen,fr,al(i,j,k),alr,den1)
END IF
ELSE
lef=fico(i-1,j,k)
ml=m1co(i-1,j,k)
END IF
IF btm(i,k+1) > j THEN
fo=0.
IF af(i,j,k) < pfo AND k > 0 THEN
mf=-ders(0.,fico(i,j,k-1),fico(i,j,k-2),abf,abf+1,abf∗(abf+1))
ELSE
mf=-ders(0.,fico(i,j,k),fico(i,j,k-1),af(i,j,k),abf,den3)
END IF
IF (af(i,j,k) < pfo OR ab(i,j,k) < pfo) AND (k=nz OR k=0) THEN mf=(fff-fb)/abf/h
IF i=0 OR i=nx THEN mf=0
ELSE
fo=fico(i,j,k+1)
mf=m3co(i,j,k+1)
END IF
IF btm(i,k-1) > j THEN
ba=0.
IF ab(i,j,k) < pfo AND k < nz THEN
mb=ders(0.,fico(i,j,k+1),fico(i,j,k+2),abf,abf+1,abf∗(abf+1))
ELSE
mb=ders(0.,fico(i,j,k),fico(i,j,k+1),ab(i,j,k),abf,den3)
END IF
IF(af(i,j,k) < pfo OR ab(i,j,k) < pfo) AND (k=nz OR k=0 ) THEN mb=(fff-fb)/abf/h
IF i=0 OR i=nx THEN mb=0
ELSE
ba=fico(i,j,k-1)
mb=m3co(i,j,k-1)
END IF
IF al(i,j,k) < pfo AND i < nx THEN
m1x=ders(m1co(i,j,k),mr,m1co(i+2,j,k),1.,2.,2.)
ELSE IF ar(i,j,k) < pfo AND i > 0 THEN
m1x=-ders(m1co(i,j,k),ml,m1co(i-2,j,k),1.,2.,2.)
ELSE
m1x=derc(m1co(i,j,k),ml,mr,al(i,j,k),ar(i,j,k),den1)
END IF
IF ad(i,j,k) < pfo THEN
m2y=ders(m2co(i,j,k),m2co(i,j+1,k),m2co(i,j+2,k),1,2,2)
ELSE
m2y=derc(m2co(i,j,k),md,m2co(i,j+1,k),ad(i,j,k),1,den2)
END IF
IF ab(i,j,k) < pfo AND k < nz THEN
m3z=ders(m3co(i,j,k),mf,m3co(i,j,k+2),1.,2.,2.)
ELSE IF af(i,j,k) < pfo AND k > 0 THEN
m3z=-ders(m3co(i,j,k),mb,m3co(i,j,k-2),1.,2.,2.)
ELSE
m3z=derc(m3co(i,j,k),mb,mf,ab(i,j,k),af(i,j,k),den3)
IF ( af(i,j,k) < pfo OR ab(i,j,k) < pfo) AND (k = nz OR k=0) THEN m3z=(mf-mb)/abf/h
END IF
rhs(i,j,k)=rhs(i,j,k)+(m1x+m2y+m3z)∗coenoto(i,j,k)∗h∗h+ficen-lef∗ \
left(i,j,k)-rig∗right(i,j,k)-d∗down(i,j,k)-fu∗up(i,j,k)-fo∗ \
forw(i,j,k)-ba∗back(i,j,k)
END IF
REPEAT LOOP
END IF
END IF
REPEAT LOOP
REPEAT LOOP
END corfi
LOOP kapp FOR fk=k0 TO numk
ka=fk
rapp=rap0
LOOP FOR fr=1 TO numh
rapp=rapp+delth
LOOP FOR ff=1 TO numf
fill=fil0+deltf∗ff
WRITE 'ly,lz.fill,rapp,rapp/h',ly,lz,fill,rapp,rapp/h
WRITE'kappa',ka
LOOP FOR i=-1 TO nx+1
LOOP FOR k=-1 TO nz+1
btm(i,k)=0
REPEAT LOOP
REPEAT LOOP
LOOP FOR i=0 TO nx
LOOP FOR k=0 TO nz
btf(i,k)=0
REPEAT LOOP
REPEAT LOOP
LOOP FOR i=-1 TO nx+1
LOOP FOR j=0 TO ny
LOOP FOR k=-1 TO nz+1
fi(i,j,k)=0.
m1(i,j,k)=0.
m2(i,j,k)=0.
m3(i,j,k)=0.
REPEAT LOOP
REPEAT LOOP
REPEAT LOOP
LOOP FOR i=0 TO nx
LOOP FOR j=0 TO ny0
LOOP FOR k=0 TO nz
ad(i,j,k)=1.
al(i,j,k)=1.
ar(i,j,k)=1.
af(i,j,k)=1.
ab(i,j,k)=1.
left(i,j,k)=0.
right(i,j,k)=0.
up(i,j,k)=0.
down(i,j,k)=0.
back(i,j,k)=0.
forw(i,j,k)=0.
bcl(i,j,k)=0.
bcr(i,j,k)=0.
bcd(i,j,k)=0.
bcu(i,j,k)=0.
bcb(i,j,k)=0.
bcf(i,j,k)=0.
coenoto(i,j,k)=0.
REPEAT LOOP
REPEAT LOOP
REPEAT LOOP
LOOP FOR i=0 TO nx
LOOP FOR j=0 TO ny
LOOP FOR k=0 TO nz
rhs(i,j,k)=0
REPEAT LOOP
REPEAT LOOP
REPEAT LOOP
LOOP FOR i=0 TO nx
x=i∗h
LOOP FOR k=0 TO nz
z=k∗h
btm(i,k)=1+FLOOR(fondo(x,z)/h)
IF ABS(1+fondo(x,z)/h- btm(i,k))> 1-h/100 THEN btm(i,k)=btm(i,k)+1
REPEAT LOOP
REPEAT LOOP
LOOP FOR k=-1 TO nz+1
btm(nx+1,k)=btm(nx-1,k)
btm(-1,k)=btm(1,k)
REPEAT LOOP
LOOP FOR i=-1 TO nx+1
btm(i,-1)=btm(i,1)
btm(i,nz+1)=btm(nx-i,nz-1)
REPEAT LOOP
LOOP FOR k=0 TO nz
LOOP FOR i=0 TO nx
btf(i,k)=MAX(btm(i,k)+1,btm(i-1,k),btm(i+1,k))
btf(i,k)=MAX(btf(i,k),btm(i,k-1),btm(i,k+1))
REPEAT LOOP
REPEAT LOOP
WRITE 'fine calcolo btf '
LOOP FOR k=0 TO nz
z=k∗h
LOOP FOR i=0 TO nx
x=h∗i
LOOP FOR j=btm(i,k) TO btf(i,k)-1
y=h∗j
IF j =btm(i,k) THEN ad(i,j,k)=j-fondo(x,z)/h
IF btm(i-1,k) > j THEN al(i,j,k)=(x-vert(y,x,x-h,z))/h
IF btm(i+1,k) > j THEN ar(i,j,k)=(vert(y,x,x+h,z)-x)/h
IF i = nx THEN ar(i,j,k)=al(i,j,k)
IF btm(i,k-1) > j THEN ab(i,j,k)=(z-vertz(y,z,z-h,x))/h
IF btm(i,k+1) > j THEN af(i,j,k)=(vertz(y,z,z+h,x)-z)/h
IF z=0 THEN ab(i,j,k)=af(i,j,k)
deny=ad(i,j,k)∗(1.+ad(i,j,k))
denx=al(i,j,k)∗ar(i,j,k)∗(al(i,j,k)+ar(i,j,k))
denz=af(i,j,k)∗ab(i,j,k)∗(af(i,j,k)+ab(i,j,k))
dud=deny/2.
dlr=denx/2.
dbf=denz/2.
axr=al(i,j,k)
axl=ar(i,j,k)
ax0=axl+axr
ayu=ad(i,j,k)
ayd=1.
ay0=ayu+ayd
azf=ab(i,j,k)
azb=af(i,j,k)
az0=azf+azb
deno=ax0∗dbf∗dud+ay0∗dbf∗dlr+az0∗dlr∗dud
bcl(i,j,k)=axl∗dud∗dbf/deno
bcd(i,j,k)=ayd∗dlr∗dbf/deno
bcr(i,j,k)=axr∗dud∗dbf/deno
bcu(i,j,k)=ayu∗dlr∗dbf/deno
bcb(i,j,k)=azb∗dud∗dlr/deno
bcf(i,j,k)=azf∗dlr∗dud/deno
coenoto(i,j,k)=dud∗dlr∗dbf/deno
REPEAT LOOP
LOOP FOR j= btf(i,k) TO ny0
bcl(i,j,k)=1./6.
bcd(i,j,k)=1./6.
bcr(i,j,k)=1./6.
bcu(i,j,k)=1./6.
bcf(i,j,k)=1./6.
bcb(i,j,k)=1./6.
coenoto(i,j,k)=1./6.
REPEAT LOOP
REPEAT LOOP
REPEAT LOOP
LOOP FOR i=0 TO nx
LOOP FOR k=0 TO nz
LOOP FOR j=btm(i,k) TO ny0
left(i,j,k)=bcl(i,j,k)
IF btm(i-1,k) > j THEN left(i,j,k)=0.
right(i,j,k)=bcr(i,j,k)
IF btm(i+1,k) > j THEN right(i,j,k)=0.
back(i,j,k)=bcb(i,j,k)
IF btm(i,k-1) > j THEN back(i,j,k)=0.
forw(i,j,k)=bcf(i,j,k)
IF btm(i,k+1) > j THEN forw(i,j,k)=0.
down(i,j,k)=bcd(i,j,k)
IF btm(i,k)= j THEN down(i,j,k)=0.
up(i,j,k)=bcu(i,j,k)
REPEAT LOOP
REPEAT LOOP
REPEAT LOOP
LOOP FOR fw=0 TO 1
wp=-1+fw∗2
LOOP FOR i=-1 TO nx+1
LOOP FOR j=0 TO ny
LOOP FOR k=-1 TO nz+1
fi(i,j,k)=0.
m1(i,j,k)=0.
m2(i,j,k)=0.
m3(i,j,k)=0.
REPEAT LOOP
REPEAT LOOP
REPEAT LOOP
IF inco=1 THEN catalogo
WRITE'inizio ciclo,btm(0,0),btm(0,1),btm(0,-1)',btm(0,0),btm(0,1),btm(0,-1)
js=0
it=1
DO
IF js=0 THEN prec=((wp+1)∗m3(0,ny,0)+(1-wp)∗m1(0,ny,0))/2
IF it<jsta THEN prec= 1000
notom1
IF it > 10 AND inco=1 THEN corm1
vbmg(m1,btm,btf,left,down,right,up,back,forw,rhs,nx,ny,nz,-wp,-wp)
notofia
IF it > 10 AND inco=1 THEN corfi
vbmg(fi,btm,btf,left,down,right,up,back,forw,rhs,nx,ny,nz,wp,-wp)
notom2
IF it > 10 AND inco=1 THEN corm2
vbmg(m2,btm,btf,left,down,right,up,back,forw,rhs,nx,ny,nz,wp,-wp)
notom3
IF it > 10 AND inco=1 THEN corm3
vbmg(m3,btm,btf,left,down,right,up,back,forw,rhs,nx,ny,nz,wp,wp)
IF js =jsta THEN
WRITE'm1c,m2c,m3c,fi',m1(0,ny,0):9.9,m2(0,ny,0),m3(0,ny,0):9.11,fi(0,ny,0),it
js=0
ELSE
js=js+1
END IF
it=it+1
IF wp = 1 THEN hlon=m3(0,ny,0)-ny∗h+fondo(0.,0.)
IF wp = -1 THEN htra=m1(0,ny,0)-ny∗h+fondo(0.,0.)
UNTIL ABS(prec-((wp+1)∗m3(0,ny,0)+(1-wp)∗m1(0,ny,0))/2) < 0.000001
REPEAT LOOP
dh=hlon-htra
WRITE 'hlon,htra,dh',hlon,htra,dh
fillm=fill
fils=fill
IF prolo =3 THEN
fillm=(fill +rapp/2 +0.144074)∗2/lz
fils=(fill+0.144074)∗2/lz
END IF
ALTEZZE rapp:10.6,fils:10.6,fillm:10.6,dh:15.10,hlon:15.10,htra:15.10
<∗ fflush(ALTEZZE_); ∗>
REPEAT LOOP
REPEAT LOOP
REPEAT kapp