Scripts
For each program the relative Example is indicated. Fortran scripts are indicated
E.10.6 Numerical solution to wave equation (fortran)
program gwaves
c parameter nk=49 ,ztop=50000.
dimension t(49),u(49),tz(49),uz(49),uzz(49),qquad(49),f(49),a(49)
dimension b(49),w(49),g(49),wreal(49),wimm(49),z(49),wpro(49)
dimension wmod(49),wfase(49),frea(49),wz(49),wrz(49),wiz(49),m(49)
dimension uprimo(49),fi(49),re(49),imm(49),fm(49),fe(49),wond(49)
real d,ko,kreale,mod,fat, norm, imm, duepi, kappa,zf,zw,m
complex i,c,w,a,b,ti,qquad,f,rapp,vento,wz,uprimo,fi,q,diff
open(10,file='w.dat',status='unknown')
open(20,file='campi.dat',status='unknown')
open(30,file='f.dat',status='unknown')
open(70,file='damp.dat',status='unknown')
open(80,file='wsmorz.dat',status='unknown')
open(90,file='flux.dat',status='unknown')
c
******************************************************************
c this quantity can be changed with nk
ztop=50000.
nk=49
h=1000.
d=ztop/(nk+1)
ztopkm=ztop/1000.
c
c velocity and damping
c
c=(40.,0.)
kappa=6.2831e-5
dam=aimag(c)*kappa
c
c gaussian parameters
c
am=10.
zw=3000.
wid=2*zw
zf=10000.
c
c amplitude at the ground
c
bot=0.
wbot=0.
c
c isothermal atmosphere (iso=1)
c
iso=0.
c
c boundary conditions at the top (lid(0) or radiation (1))
c
cond=0.
c
print*,'these are the paramters used:'
print*,'1)number of levels :',nk
print*,'2)phase velocity:',c
print*,'3)gaussian amplitude:',am
print*,'4) center of the gaussian ',zf,'m'
print*,'5)width of the gaussian:',wid,'m'
print*,'6) vertical wind at the ground :',wbot,'m/s'
c
c define the imaginary unit and pie
c
i=(0.,1.)
duepi=8.*atan(1.)
c
c defined the temperature (k) and wind fields (m/s)
c
tbot=280.
ubot=0.
ttop=270.
utop=-40.
luno=int(12000./d)
ldue=int(15000./d)
do k=1,luno
z(k)=k*d
t(k)=280.-(6.5*z(k)/h)
u(k)=3.3333*z(k)/h
end do
do k=luno+1,ldue
z(k)=k*d
t(k)=202.
u(k)=-2.1052*((z(k)-z(luno))/h)+40.
end do
do k=ldue+1,nk
z(k)=k*d
t(k)=((z(k)-z(ldue))/h)*1.9428+202.0
u(k)=-2.1052*((z(k)-z(luno))/h)+40.0
end do
c
c derivatives of t and u with the appropriate routines
c
if(iso.eq.1) then
call isothermal (t, tz, nk,tbot,ttop)
print*,'7) isothermal atmosphere'
else
call firstderivative (t,nk,d,tz)
tz(1)=-6.5/d
tz(49)=1.9428/d
end if
call firstderivative (u,nk,d,uz)
uz(1)=3.3333/d
uz(49)=-2.1052/d
call seconderivative(u,nk,d,uzz)
uzz(1)=0.
uzz(49)=0.
c
c introduce the quantities q*q e f(z)
c
h=7.4*h
r=287.
r=r/h
gamma=7./5.
ko=(gamma-1)/gamma
am=am/(h*h)
m(1)=1.
do k=1,nk
z(k)=k*d
qquad(k)=r*(tz(k)+ko*t(k)/h)/((c-u(k))**2)
qquad(k)=qquad(k)+(uzz(k)+uz(k)/h)/(c-u(k))
zot=0.25/(h*h)
qquad(k)=qquad(k)-zot
c
c calculation of the modulus k-kf and the exponent of the gaussian (fat)
c
mod=sqrt((z(k)-zf)**2.)
fat=((z(k)-zf)/zw)**2.
c
c forcing term is zero outside twice zw
c
e=z(k)/h
m(k)=exp(-e/2)
if(mod.ge.wid) then
f(k)=0.
g(k)=0.
else
g(k)=am*(exp(-fat)-exp(-4.))
f(k)=m(k)*g(k)
frea(k)=real(f(k))
end if
vvv=z(k)/1000.
write(30,'(4(1x,e11.4))') vvv,frea(k),g(k),m(k)
end do
c
c gaussian algorithm for the solution
c
a(1)=(d*d)*qquad(1)-2.
b(1)=(d*d)*f(1)-wbot
do k=2,nk-1
a(k)=(d*d)*qquad(k)-2.
b(k)=(d*d)*f(k)
end do
c
c top boundary conditions
c
if(cond.eq.0) then
print*,'condition at the top: artificial lid'
a(49)=(d*d)*qquad(49)-2.
b(49)=(d*d)*f(49)
else
print*,'condition at the top: radiation'
q=sqrt(qquad(49)/2.)-(d*q*i)-1.
b(49)=0.
end if
c
c solution of the system
c
do k=2,nk
ti=1./a(k-1)
a(k)=a(k)-ti
b(k)=b(k)-ti*b(k-1)
end do
c
c calcolo di w
c
w(49)=b(49)/a(49)
do j=1,nk-1
k=nk-j
w(k)=(b(k)-w(k+1))/a(k)
end do
c
c print module and phase of w(k)
c
write(10,'(3(1x,e11.4))') bot,wbot,bot
write(20,'((1x,i3),5(1x,f10.4))') k,bot,tbot,ubot,tz(1),uz(1)
do k=1,49
e=z(k)/h
vvv=z(k)/1000.
norm=exp(-e/2.)
wreal(k)=real(w(k))
wimm(k)=aimag(w(k))
wond(k)=sqrt(wreal(k)**2.+wimm(k)**2)
wreal(k)=real(w(k))/norm
wimm(k)=aimag(w(k))/norm
wmod(k)=sqrt(wreal(k)**2+wimm(k)**2)
wpro(k)=wmod(k)*cos((duepi)*(z(k)/5000.))
write(10,'(3(1x,e11.4))') vvv,wmod(k),g(k)
write(20,'(1(1x,i3),5(1x,f10.4))') k,vvv,t(k),u(k),tz(k),uz(k)
end do
write(20,'((1x,i3),5(1x,f10.4))') k,ztop km, ttop,utop,tz(k),uz(k)
c
c calculation of dw/dz and uprimo and fi
c
call firstderivative (wreal,nk,d,wrz)
call firstderivative (wimm,nk,d,wiz)
wrz(1)=0.
wrz(49)=wrz(nk-1)
wiz(1)=0.
wiz(49)=wiz(nk-1)
do k=1,nk
wz(k)=wrz(k)+i*wiz(k)
uprimo(k)=-i*(w(k)/h-wz(k))/kappa
fi(k)=(c-u(k))*uprimo(k)
fi(k)=fi(k)+i*(w(k)*uz(k))/kappa
end do
c
c momentum and energy fluxes
c
do k=1,nk
re(k)=real(uprimo(k))
imm(k)=aimag(uprimo(k))
end do
call flusso(z,re,imm,wreal,wimm,fm,nk,h,d)
do k=1,nk
re(k)=real(fi(k))
imm(k)=aimag(fi(k))
end do
call flusso(z,re,imm,wreal,wimm,fe,nk,h,d)
do k=1,nk
vvv=z(k)/1000.
write(90,'(3(1x,e11.4))') vvv,fm(k),fe(k)
end do
c
c damping term
c
do j=0,48
sec=3600.*j
smorz=exp(-sec*dam)
write(70,'(1x,i2,1x,e11.4)')j,smorz
write(80,'(1x,e11.4)')wbot
do k=1,nk
wsmorz=wmod(k)*smorz
write(80,'(1x,e11.4)') wsmorz
end do
end do
stop
end
c
c Subroutine firstderivative
c
subroutine firstderivative (vett,nk,d,devp)
dimension devp(49),vett(49)
do k=2,nk-1
devp(k)=(vett(k+1)-vett(k-1))/(2.*d)
end do
return
end
c
c Subroutine seconderivative
c
subroutine seconderivative(vett,nk,d,devs)
dimension devs(49),vett(49)
do k=2,nk-1
devs(k)=(vett(k+1)-2*vett(k)+vett(k-1))/(d**2)
end do
return
end
c
c Subroutine flusso
c
subroutine flusso(z,re,imm,wreal,wimm,flux,nk,h,d)
dimension z(49),re(49),imm(49),wreal(49),wimm(49),flux(49)
real med,flux,norm,imm
do k=1,nk
z(k)=k*d
e=z(k)/h
vvv=z(k)/1000.
norm=exp(-e)
med=(re(k)*wreal(k)+imm(k)*wimm(k))/2.
flux(k)=med*norm/(h**2)
end do
return
end
c
c Subroutine isothermal
c
subroutine isothermal (t,tz,nk,tbot,ttop)
dimension t(49),tz(49)
tz(k)=0.
do k=1,nk
t(k)=300.
end do
ttop=t(nk)
tbot=t(nk)
return
end