program map04
use dcl
integer,parameter :: nx=36, ny=36
real,parameter :: xmin= 0, xmax=360, ymin=-90, ymax=+90
real,parameter :: pi=3.141592, drad=pi/180, dz=0.05
integer,parameter :: fact=10
real,dimension(0:nx,0:ny) :: p,u,v
real :: alon(0:nx), alat(0:ny)
rmiss = DclGetReal( 'GLOBAL:RMISS' )
call DclSetParm( 'GLOBAL:lmiss', .true. )
alon = xmin + (xmax-xmin) * (/(i, i=0,nx)/) /nx
alat = ymin + (ymax-ymin) * (/(j, j=0,ny)/) /ny
do j = 0, ny
do i = 0, nx
slat = sin(alat(j)*drad)
p(i,j) = cos(alon(i)*drad) * (1-slat**2) * sin(2*pi*slat) + dz
end do
end do
do j = 1, ny-1
do i = 0, nx
u(i,j) = ( p(i,j-1) - p(i,j+1) ) * fact
v(i,j) = ( p(modulo(i,nx)+1,j) - p(modulo(i-1,nx),j) )* fact
end do
end do
u(:,1) = rmiss; u(:,ny) = rmiss
v(:,1) = rmiss; v(:,ny) = rmiss
call DclOpenGraphics()
call DclNewFrame
call DclSetWindow( xmin, xmax, ymin, ymax )
call DclSetViewPort( 0.1, 0.9, 0.1, 0.9 )
call DclSetSimilarity( 0.4, 0.0, 0.0 )
call DclSetMapProjectionAngle( 165.0, 60.0, 0.0 )
call DclSetMapProjectionWindow( -180.0, 180.0, 0.0, 90.0 )
call DclSetTransNumber( 30 )
call DclSetTransFunction
call DclSetParm( 'GRAPH:lclip', .true. )
call DclDrawMap( 'coast_world' )
call DclDrawGlobe()
call DclDrawContour( p )
do j=0,ny
do i=0,nx
if (.not.(u(i,j).eq.rmiss .or. v(i,j).eq.rmiss)) then
call DclDrawArrow( &
alon(i),alat(j),alon(i)+u(i,j),alat(j)+v(i,j), 1, 3 )
end if
end do
end do
call DclCloseGraphics
end program
|
