Monday, December 1, 2008

Vertical Velocity calculation using fortran

c program to find vertical velocity
integer n,p(10)
real lat1,lat2,lon1,lon2,lat,dx,dy,tem,div(10),avgdiv(10),dp
real omega,w
dimension ddd(4),ff(4)
print*,'enter the no of pressure levels'
read(*,*)n
print*,'enter the latitudes of south'
read(*,*)lat1
print*,'enter the latitudes of north'
read(*,*)lat2
print*,'enter the longitude of west'
read(*,*)lon1
print*,'enter the longitude of east'
read(*,*)lon2
lat=(lat1+lat2)/2
dx=(lon2-lon1)*cos(lat)*110*1000
dy=(lat2-lat1)*110*1000
write(*,*)'enter the temp in kelvin'
read(*,*)tem
open(1,file='vert.dat')
do i=1,n
read(1,*)p(i)
write(*,*)p(i)
do j=1,4
read(1,*)ddd(j),ff(j)
write(*,*)ddd(j),ff(j)
enddo
ue=0.5148*(-1)*ff(2)*sin(ddd(2)*3.14/180)
uw=0.5148*(-1)*ff(4)*sin(ddd(4)*3.14/180)
vn=0.5148*(-1)*ff(1)*cos(ddd(1)*3.14/180)
vs=0.5148*(-1)*ff(3)*cos(ddd(3)*3.14/180)
div(i)=(ue-uw)/dx+(vn-vs)/dy

enddo
omega=0.0
do i=1,n-1
avgdiv(i)=(div(i)+div(i+1))/2
dp=(p(i)-p(i+1))
omega=omega+avgdiv(i)*dp
enddo
close(1)

write(*,2)omega
w=((-1)*omega*287*tem)/(9.8*p(n))*100
write(*,3)w

2 format('vertical velocity=',f10.8,' in pressure co-ordinate' )
3 format('vertical velocity=',f10.4,' cm/s' )

stop
end

No comments: