!-----------------------------------------------------------------------
!     Copyright (c) 2002 GFD Dennou Club. All Rights Reserved.
!-----------------------------------------------------------------------
program sample_3d06

  use dcl
  integer,parameter :: nx=37, ny=37
  real,parameter ::  xmin=   0,  xmax=360,  ymin= -90,  ymax= 90
  real,parameter :: vxmin=-0.4, vxmax=0.4, vymin=-0.3, vymax=0.3
  real,parameter ::  zmin=-1.9,  zmax=1.9
  real,parameter :: vzmin=-0.2, vzmax=0.2
  real,parameter :: pi=3.14159, drad=pi/180

  real :: alon(nx), alat(ny), p(nx,ny)
  real :: xp(3), yp(3), zp(3)

    do i=1,nx
      alon(i) = xmin + (xmax-xmin)*(i-1)/(nx-1)
    end do
    do j=1,ny
      alat(j) = ymin + (ymax-ymin)*(j-1)/(ny-1)
    end do

    do j=1,ny
      slat = sin(alat(j)*drad)
      do i=1,nx
        p(i,j) = 3*sqrt(1-slat**2)*slat*cos(alon(i)*drad)-0.5*(3*slat**2-1)
      end do
    end do

    call DclOpenGraphics()
    call DclNewFrame

!-- x-y plane ----

    call DclSetWindow(  xmin,  xmax,  ymin,  ymax )
    call DclSetViewPort( vxmin, vxmax, vymin, vymax )
    call DclSetTransNumber( 1 )
    call DclSetTransFunction

    call DclSet2DPlane( 1, 2, vzmin)
    call DclSet3DEyePoint( -1.1, -1.1, 2.5 )
    call DclSet3DObjectPoint(  0.0,  0.0, 0.0 )
    call DclSet3DProjection

    call DclDrawAxis( 'bt', 10., 60. )
    call DclDrawTitle( 'b', 'longitude', 0. )
    call DclDrawAxis( 'lr', 10., 30. )
    call DclDrawTitle( 'l', 'latitude', 0. )

!-- x-z plane ----

    call DclSetWindow(  xmin,  xmax,  zmin,  zmax )
    call DclSetViewPort( vxmin, vxmax, vzmin, vzmax )
    call DclSetTransNumber( 1 )
    call DclSetTransFunction

    call DclSet2DPlane( 1, 3, vymax)
    call DclSet3DProjection

    call uzinit
    call DclDrawAxis( 't', 10., 60. )
    call DclSetParm( 'AXIS:labelxb', .false. )
    call DclDrawAxis( 'b', 10., 60. )
    call DclDrawAxis( 'lr', 0.2, 1.0 )
    call DclDrawTitle( 'l', 'amplitude', 0. )

!-- y-z plane ----

    call DclSetWindow(  ymin,  ymax,  zmin,  zmax )
    call DclSetViewPort( vymin, vymax, vzmin, vzmax )
    call DclSetTransNumber( 1 )
    call DclSetTransFunction

    call DclSet2DPlane( 2, 3, vxmax)
    call DclSet3DProjection

    call uzinit
    call DclDrawAxis( 'tb', 10., 30. )
    call DclSetParm( 'AXIS:labelyl', .false. )
    call DclDrawAxis( 'lr', 0.2, 1.0 )

!---------------- 3-d ------------------

    call DclSet3DViewPort(vxmin, vxmax, vymin, vymax, vzmin, vzmax)
    call DclSet3DWindow( xmin,  xmax,  ymin,  ymax,  zmin,  zmax)
    call DclSet3DLogAxis(.false., .false., .false.)
    call DclSet3DTransNumber(1)
    call DclSet3DTransFunction

!    call szl3op(1)
!      call szt3op(2999,4999)

    do j=ny-1, 1, -1
     do i=nx-1, 1, -1
      xp(1) = alon(i)
      yp(1) = alat(j)
      zp(1) = p(i,j)
      xp(2) = alon(i)
      yp(2) = alat(j+1)
      zp(2) = p(i,j+1)
      xp(3) = alon(i+1)
      yp(3) = alat(j+1)
      zp(3) = p(i+1,j+1)

!        call szt3zu(xp, yp, zp)
!      call szl3zu(3, xp, yp, zp)
      call DclDraw3DLine( xp, yp, zp )

      xp(1) = alon(i+1)
      yp(1) = alat(j+1)
      zp(1) = p(i+1,j+1)
      xp(2) = alon(i+1)
      yp(2) = alat(j)
      zp(2) = p(i+1,j)
      xp(3) = alon(i)
      yp(3) = alat(j)
      zp(3) = p(i,j)

!        call szt3zu(xp, yp, zp)
!      call szl3zu(3, xp, yp, zp)
      call DclDraw3DLine( xp, yp, zp )
     end do
    end do

!    call szl3cl
!      call szt3cl

    call DclCloseGraphics

end program
