*deck wenos
c-----------------------------------------------------------
csssssssssssssssssssssssssssssssssssssssssssssssssssssssssss
c
subroutine wenos(idirex)
implicit none
c
c------------------------------------------------------------
include 'com_main.h'
include 'params.h'
c------------------------------------------------------------
c
integer idirex
integer j , i
integer ip , I am , jp , jm
c
real*8 ucon , vcon , ftx , fty , eno ,
& tc , tr , tl , tf , ta ,
& dtr , dtl , d2tcx, d2tlx, d2trx ,
& dtrr , dtll , dtff , dtaa ,
& dtf , dta , d2tcy, d2tfy, d2tay ,
& ur , ul , vf , va , rmnmd,
& trr , tll , tff , taa
real*8 S1,S2,S3,v1,v2,v3,v4,v5,phix1,phix2,phix3
real*8 a11,a22,a33,a1,a2,a3,epss,tlll,taaa,trrr
real*8 tfff,dtrrr,dtlll ,phixr,phixl,phixt,phixb
real*8 dtfff,dtaaa,tempp,tempq,tempr,temps,tempt,tempu
c
c-----all parameters used here are defined----------------------------------------------------------------
c
eno = one
c
do 1000 j = 2, jbp1
c
jp = j + 1
jm = j - 1
c
do 1000 i = 2, ibp1
c
ip = i + 1
I am = i - 1
ur = u(i ,j )
ul = u(im,j )
vf = v(i ,j )
va = v(i ,jm)
c
tc = lsfn(i ,j )
tr = lsfn(ip ,j )
tl = lsfn(im ,j )
tf = lsfn(i ,jp )
ta = lsfn(i ,jm )
trr = lsfn(ip+1,j )
tll = lsfn(im-1,j )
tff = lsfn(i ,jp+1)
taa = lsfn(i ,jm-1)
tlll = lsfn(im-2,j)
taaa = lsfn(i,jm-2)
trrr = lsfn(ip+2,j)
tfff = lsfn(i,jp+2)
c
c------------------------------------------------------------
c >>> x-sweep
c------------------------------------------------------------
c
if( idirex .eq. 1 ) then
c
ucon = half*( ul + ur )
c
dtr = rdxy * (tr - tc)
dtl = rdxy * (tc - tl)
c
if( i .lt. ibr ) then
dtrr = (trr - tr) * rdxy
else
dtrr = dtr
endif
c
if( i .gt. 2 ) then
dtll = (tl - tll) * rdxy
else
dtll = dtl
endif
c
if( i .gt. 3 ) then
dtlll = (tll - tlll) * rdxy
else
dtlll = dtll
endif
c
if( i .lt. ibr-1 ) then
dtrrr = (trrr - trr) * rdxy
else
dtrrr = dtrr
endif
c
v1 = dtlll
v2 = dtll
v3 = dtl
v4 = dtr
v5 = dtrr
tempp = v1 − 2.0d0*v2 + v3
tempq = v1 − 4.0d0*v2 + 3.0d0*v3
tempr = v2 − 2.0d0*v3 + v4
temps = v2 − v4
tempt = v3 − 2.0d0*v4 + v5
tempu = 3.0d0*v3 − 4.0d0*v4 + v5
S1=(13./12.)*(tempp**2.)+(1./4.)*(tempq**2.)
S2=(13./12.)*(tempr**2.)+(1./4.)*(temps**2.)
S3=(13./12.)*(tempt**2.)+(1./4.)*(tempu**2.)
phix1=(v1/6.)-(7.*v2/6.)+(11.*v3/6.)
phix2=(-v2/6.)-(5.*v3/6.)+(v4/3.)
phix3=(-v3/3.)-(5.*v4/6.)+(v5/6.)
tempu=max(v1**2,v2**2,v3**2,v4**2,v5**2)
epss=(10.0d-6)*tempu + (10.0d-99)
a11=0·1/((S1+epss)**2)
a22=0·6/((S2+epss)**2)
a33=0·3/((S3+epss)**2)
a1=a11/(a11+a22+a33)
a2=a22/(a11+a22+a33)
a3=a33/(a11+a22+a33)
phixl=a1*phix1+a2*phix2+a3*phix3
c
v1=dtll
v2=dtl
v3=dtr
v4=dtrr
v5=dtrrr
tempp=v1 − 2.0d0*v2 + v3
tempq=v1 − 4.0d0*v2 + 3.0d0*v3
tempr=v2 − 2.0d0*v3 + v4
temps=v2 − v4
tempt=v3 − 2.0d0*v4 + v5
tempu=3.0d0*v3 − 4.0d0*v4 + v5
S1=(13./12.)*(tempp**2.)+(1./4.)*(tempq**2.)
S2=(13./12.)*(tempr**2.)+(1./4.)*(temps**2.)
S3=(13./12.)*(tempt**2.)+(1./4.)*(tempu**2.)
phix1=(v1/6.)-(7.*v2/6.)+(11.*v3/6.)
phix2=(-v2/6.)-(5.*v3/6.)+(v4/3.)
phix3=(-v3/3.)-(5.*v4/6.)+(v5/6.)
tempu=max(v1**2,v2**2,v3**2,v4**2,v5**2)
epss=(10.0d-6)*tempu + (10.0d-99)
a11=0·1/((S1+epss)**2)
a22=0·6/((S2+epss)**2)
a33=0·3/((S3+epss)**2)
a1=a11/(a11+a22+a33)
a2=a22/(a11+a22+a33)
a3=a33/(a11+a22+a33)
phixr=a1*phix1+a2*phix2+a3*phix3
c
if( ucon .gt. zero ) then
ftx = ucon * phixl
else
ftx = ucon * phixr
endif
c
lsf(i,j) = lsfn(i,j) - dt * ftx
c
endif
c------------------------------------------------------------
c >>> y-sweep
c------------------------------------------------------------
c
if ( idirex .eq. 2 ) then
c
vcon = half*( va + vf )
c
dtf = rdxy * (tf - tc)
dta = rdxy * (tc - ta)
c
if( j .lt. jbp1 ) then
dtff = (tff - tf) * rdxy
else
dtff = dtf
endif
c
if( j .gt. 2 ) then
dtaa = (ta - taa) * rdxy
else
dtaa = dta
endif
c
if( j .gt. 3 ) then
dtaaa = (taa - taaa) * rdxy
else
dtaaa = dtaa
endif
c
if( j .lt. jbr-1 ) then
dtfff = (tfff - tff) * rdxy
else
dtfff = dtff
endif
c
v1=dtaaa
v2=dtaa
v3=dta
v4=dtf
v5=dtff
tempp=v1 − 2.0d0*v2 + v3
tempq=v1 − 4.0d0*v2 + 3.0d0*v3
tempr=v2 − 2.0d0*v3 + v4
temps=v2 − v4
tempt=v3 − 2.0d0*v4 + v5
tempu=3.0d0*v3 − 4.0d0*v4 + v5
S1=(13./12.)*(tempp**2.)+(1./4.)*(tempq**2.)
S2=(13./12.)*(tempr**2.)+(1./4.)*(temps**2.)
S3=(13./12.)*(tempt**2.)+(1./4.)*(tempu**2.)
phix1=(v1/6.)-(7.*v2/6.)+(11.*v3/6.)
phix2=(-v2/6.)-(5.*v3/6.)+(v4/3.)
phix3=(-v3/3.)-(5.*v4/6.)+(v5/6.)
tempu=max(v1**2,v2**2,v3**2,v4**2,v5**2)
epss=(10.0d-6)*tempu + (10.0d-99)
a11=0·1/((S1+epss)**2)
a22=0·6/((S2+epss)**2)
a33=0·3/((S3+epss)**2)
a1=a11/(a11+a22+a33)
a2=a22/(a11+a22+a33)
a3=a33/(a11+a22+a33)
phixb=a1*phix1+a2*phix2+a3*phix3
c
v1=dtaa
v2=dta
v3=dtf
v4=dtff
v5=dtfff
tempp=v1 − 2.0d0*v2 + v3
tempq=v1 − 4.0d0*v2 + 3.0d0*v3
tempr=v2 − 2.0d0*v3 + v4
temps=v2 − v4
tempt=v3 − 2.0d0*v4 + v5
tempu=3.0d0*v3 − 4.0d0*v4 + v5
S1=(13./12.)*(tempp**2.)+(1./4.)*(tempq**2.)
S2=(13./12.)*(tempr**2.)+(1./4.)*(temps**2.)
S3=(13./12.)*(tempt**2.)+(1./4.)*(tempu**2.)
phix1=(v1/6.)-(7.*v2/6.)+(11.*v3/6.)
phix2=(-v2/6.)-(5.*v3/6.)+(v4/3.)
phix3=(-v3/3.)-(5.*v4/6.)+(v5/6.)
tempu=max(v1**2,v2**2,v3**2,v4**2,v5**2)
epss=(10.0d-6)*tempu + (10.0d-99)
a11=0·1/((S1+epss)**2)
a22=0·6/((S2+epss)**2)
a33=0·3/((S3+epss)**2)
a1=a11/(a11+a22+a33)
a2=a22/(a11+a22+a33)
a3=a33/(a11+a22+a33)
phixt=a1*phix1+a2*phix2+a3*phix3
c
c d2tcy = rdxy * (dtf - dta)
c d2tfy = rdxy * (dtff - dtf)
c d2tay = rdxy * (dta - dtaa)
c
if( vcon .gt. zero ) then
fty = vcon * phixb
else
fty = vcon * phixt
endif
c
lsf(i,j) = lsfn(i,j) - dt * fty
c
endif
c____________________________________________________________
c
1000 continue
c
return