
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
