!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE IMG3D                      ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE img3d(b,x,y,z,                                               &,16
           nx,ibgn,iend,isk, ny,jbgn,jend,jsk, nz,kbgn,kend,ksk,        &
           bmax,bmin, label,fnkey,time,mode,kslice,jslice,islice,       &
           n,xp,yp,b1,b2,zs2, mgrid, nestgrd, tem1)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Produce 2-D HDF images for specified slices of 3-d array b.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!
!  MODIFICATION HISTORY:
!
!  9/1/94 (Y. Lu)
!  Cleaned up documentation
!
!  2/1/1995 (J. Zhang)
!    Added constant p-level plots (mode.eq.6)
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    b        3-dimensional array of data
!
!    x        x coordinate of grid points in physical/comp. space (m)
!    y        y coordinate of grid points in physical/comp. space (m)
!    z        z coordinate of grid points in computational space (m)
!
!    nx       first dimension of b
!    ibgn     index of first i grid point to be used.
!    iend     index of last  i grid point to be used.
!    isk      skip or interpolation parameter in i-dir.
!             If isk >0, every isk'th data is used,
!             If isk <0, each grid zone is sub-divided into isk zones.
!
!    ny       second dimension of b
!    jbgn     index of first j grid point to be used.
!    jend     index of last  j grid point to be used.
!    jsk      skip or interpolation parameter in j-dir.
!             If jsk >0, every jsk'th data is used,
!             If jsk <0, each grid zone is sub-divided into jsk zones.
!
!    nz       third dimension of b
!    kbgn     index of first k grid point to be used.
!    kend     index of last  k grid point to be used.
!    ksk      skip or interpolation parameter in k-dir.
!             If ksk >0, every ksk'th data is used,
!             If ksk <0, each grid zone is sub-divided into ksk zones.
!
!    label    character string describing the contents of b
!    fnkey    file name key
!    time     time of data in seconds
!
!    mode     slice orientation indicator
!             mode = 1, x-y slice at z index kslice is plotted.
!             mode = 2, x-z slice at y index jslice is plotted.
!             mode = 3, y-z slice at x index islice is plotted.
!             mode = 0, all of the three slices above are plotted.
!
!    kslice   k index of plane for mode=1 x-y slice
!    jslice   j index of plane for mode=2 x-z slice
!    islice   i index of plane for mode=1 y-z slice
!
!  OUTPUT :
!
!    n        horizontal dimension of arbitary vertical cross section
!    xp       x-coordinate of grid points on arbitary
!             vertical cross-section
!    yp       y-coordinate of grid points on arbitary
!             vertical cross-section
!    b1       2-D field interpolated to a horizontal cross-section
!    b2       2-D field interpolated to a vertical corss-section
!    zs2      z-coordinate of grid points on arbitary vertical
!             cross-section
!    mgrid    the grid number
!    nestgrd  flag for nested grid run.
!
!  WORK ARRAYS:
!
!    tem1     Temporary work array.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz       ! dimemsions of 3-D input data arrray b
  INTEGER :: n              ! horizontal dimension of arbitary
                            ! vertical cross section.

  REAL :: b(nx,ny,nz)       ! 3-dimensional array of data
  REAL :: x(nx,ny,nz)       ! x coordinate of grid points in physical
                            ! /comp. space (m)
  REAL :: y(ny,ny,nz)       ! y coordinate of grid points in physical
                            ! /comp. space (m)
  REAL :: z(nz,ny,nz)       ! z coordinate of grid points in physical
                            ! /comp. space (m)

  REAL :: bmax,bmin         ! upper and lower bounds of b used in IMG2D

  INTEGER :: kslice         ! k index of plane for mode=1 x-y slice
  INTEGER :: jslice         ! j index of plane for mode=1 x-z slice
  INTEGER :: islice         ! i index of plane for mode=1 y-z slice

  REAL :: b1(nx,ny)         ! 2-D field interpolated to a horizontal
                            ! cross-section
  REAL :: b2(n,nz)          ! 2-D field interpolated to a vertical
                            ! corss-section.
  REAL :: zs2(n,nz)         ! z-coordinate of grid points on arbitary
                            ! vertical cross-section
  REAL :: xp(n)             ! x-coordinate of grid points on arbitary
                            ! vertical cross-section
  REAL :: yp(n)             ! y-coordinate of grid points on arbitary
                            ! vertical cross-section

  INTEGER :: ibgn,iend,isk  ! index of first and last i grid point;
                            ! skip/interpolation parameter in i-dir
  INTEGER :: jbgn,jend,jsk  ! index of first and last j grid point;
                            ! skip/interpolation parameter in j-dir
  INTEGER :: kbgn,kend,ksk  ! index of first and last k grid point;
                            ! skip/interpolation parameter in k-dir
  INTEGER :: mode           ! slice orientation indicator
  INTEGER :: mgrid          ! the grid number
  INTEGER :: nestgrd        ! flag for nested grid run.
  REAL :: tem1(*)           ! temporary work array

  INTEGER :: length         ! the length of string

  CHARACTER (LEN=*   ) :: label ! character string describing the contents
                                ! of b
  CHARACTER (LEN=80  ) :: timsnd  ! character string describing the time of b
  CHARACTER (LEN=*   ) :: fnkey ! file name key
  INTEGER :: tmstrln
  REAL :: time              ! time of data (s)
!
!-----------------------------------------------------------------------
!
!  Common blocks for plotting control parameters
!
!-----------------------------------------------------------------------
!
  REAL :: x01,y01            ! the first  point of interpolation
  REAL :: x02,y02            ! the second point of interpolation
  REAL :: z01                ! the given height of the slice
  REAL :: sinaf,cosaf,dist,sqrtdxy
  COMMON /slicev/x01,y01,x02,y02,sinaf,cosaf,dist,sqrtdxy
  COMMON /sliceh/z01
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k,ij,ik,jk,idata,jdata,kdata
  CHARACTER (LEN=120) :: imgfn
  CHARACTER (LEN=35) :: gridnum
  DATA gridnum /'123456789abcdefghijklmnopqrstuvwxyz'/
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!  mode=1   Plot x-y cross-section
!
!-----------------------------------------------------------------------
!
  CALL cvttsnd( time, timsnd, tmstrln )

  IF(mode == 1 .OR. mode == 0 ) THEN

    k = kslice
    DO j=jbgn,jend
      DO i=ibgn,iend
        ij = i-ibgn+1 + (j-jbgn)*(iend-ibgn+1)
        tem1(ij)=b(i,j,k)
      END DO
    END DO

    imgfn = fnkey
    length = LEN( fnkey )
    WRITE(imgfn(length+1:length+5),'(a,i3.3)') '.k',k
    imgfn(length+6:120) = label//timsnd(1:tmstrln)
    length = 120
    CALL strlnth( imgfn, length)

    IF( nestgrd == 1 ) THEN
      WRITE(imgfn(length+1:length+3),'(''.G'',A)')                      &
            gridnum(mgrid:mgrid)
      length = length + 3
    END IF

    idata = iend-ibgn+1
    jdata = jend-jbgn+1

    CALL img2d(tem1,idata,1,idata,isk,jdata,1,jdata,isk,                &
               bmax,bmin,                                               &
               imgfn(1:length))

!
!-----------------------------------------------------------------------
!
!  mode=2   Plot x-z cross-section
!
!-----------------------------------------------------------------------
!
  ELSE IF (mode == 2 .OR. mode == 0 ) THEN

    j = jslice
    DO k=kbgn,kend
      DO i=ibgn,iend
        ik = i-ibgn+1 + (k-kbgn)*(iend-ibgn+1)
        tem1(ik)=b(i,j,k)
      END DO
    END DO

    imgfn = fnkey
    length = LEN( fnkey )
    WRITE(imgfn(length+1:length+5),'(a,i3.3)') '.j',j
    imgfn(length+6:120) = label//timsnd(1:tmstrln)
    length = 120
    CALL strlnth( imgfn, length)

    IF( nestgrd == 1 ) THEN
      WRITE(imgfn(length+1:length+3),'(''.G'',A)')                      &
            gridnum(mgrid:mgrid)
      length = length + 3
    END IF

    idata = iend-ibgn+1
    kdata = kend-kbgn+1

    CALL img2d(tem1,idata,1,idata,isk,kdata,1,kdata,ksk,                &
               bmax,bmin,                                               &
               imgfn(1:length))

!
!-----------------------------------------------------------------------
!
!  mode=3   Plot y-z cross-section
!
!-----------------------------------------------------------------------
!
  ELSE IF ( mode == 3 .OR. mode == 0) THEN

    i = islice
    DO k=kbgn,kend
      DO j=jbgn,jend
        jk = j-jbgn+1 + (k-kbgn)*(jend-jbgn+1)
        tem1(jk)=b(i,j,k)
      END DO
    END DO

    imgfn = fnkey
    length = LEN( fnkey )
    WRITE(imgfn(length+1:length+5),'(a,i3.3)') '.i',i
    imgfn(length+6:120) = label//timsnd(1:tmstrln)
    length = 120
    CALL strlnth( imgfn, length)

    IF( nestgrd == 1 ) THEN
      WRITE(imgfn(length+1:length+3),'(''.G'',A)')                      &
            gridnum(mgrid:mgrid)

      length = length + 3
    END IF

    jdata = jend-jbgn+1
    kdata = kend-kbgn+1

    CALL img2d(tem1,jdata,1,jdata,jsk,kdata,1,kdata,ksk,                &
               bmax,bmin,                                               &
               imgfn(1:length))
!
!-----------------------------------------------------------------------
!
!  mode=4   Plot horizontal slice at given height
!
!-----------------------------------------------------------------------
!
  ELSE IF( mode == 4 ) THEN

    CALL secthrz(nx,ny,nz,b,z,b1)

    DO j=jbgn,jend
      DO i=ibgn,iend
        ij = i-ibgn+1 + (j-jbgn)*(iend-ibgn+1)
        tem1(ij)=b1(i,j)
      END DO
    END DO

    imgfn = fnkey
    length = LEN( fnkey )
    WRITE(imgfn(length+1:length+5),'(a,i3.3)') '.z',nint(z01)
    imgfn(length+6:120) = label//timsnd(1:tmstrln)
    length = 120
    CALL strlnth( imgfn, length)

    IF( nestgrd == 1 ) THEN
      WRITE(imgfn(length+1:length+3),'(''.G'',A)')                      &
            gridnum(mgrid:mgrid)

      length = length + 3
    END IF

    idata = iend-ibgn+1
    jdata = jend-jbgn+1

    CALL img2d(tem1,idata,1,idata,isk,jdata,1,jdata,jsk,                &
               bmax,bmin,                                               &
               imgfn(1:length))
!
!-----------------------------------------------------------------------
!
!  mode=5   Plot vectical slice through two given points
!
!-----------------------------------------------------------------------
!
  ELSE IF( mode == 5 ) THEN

    CALL sectvrt(nx,ny,nz,b,x,y,z,b2,zs2,n,xp,yp)

    DO k=kbgn,kend
      DO i=ibgn,iend
        ik = i-ibgn+1 + (k-kbgn)*(iend-ibgn+1)
        tem1(ik)=b2(i,k)
      END DO
    END DO

    imgfn = fnkey
    length = LEN( fnkey )
    WRITE(imgfn(length+1:length+5),'(a,i3.3)') '.v',nint(x01)
    imgfn(length+6:120) = label//timsnd(1:tmstrln)
    length = 120
    CALL strlnth( imgfn, length)

    IF( nestgrd == 1 ) THEN
      WRITE(imgfn(length+1:length+3),'(''.G'',A)')                      &
            gridnum(mgrid:mgrid)

      length = length + 3
    END IF

    idata = iend-ibgn+1
    kdata = kend-kbgn+1

    CALL img2d(tem1,idata,1,idata,isk,kdata,1,kdata,ksk,                &
               bmax,bmin,                                               &
               imgfn(1:length+4))

!
!-----------------------------------------------------------------------
!
!  mode=6   Plot on constant pressure-level
!
!-----------------------------------------------------------------------
!
  ELSE IF( mode == 6 ) THEN

    CALL secthrz(nx,ny,nz,b,z,b1)

    DO j=jbgn,jend
      DO i=ibgn,iend
        ij = i-ibgn+1 + (j-jbgn)*(iend-ibgn+1)
        tem1(ij)=b1(i,j)
      END DO
    END DO

    imgfn = fnkey
    length = LEN( fnkey )
    WRITE(imgfn(length+1:length+5),'(a,i3.3)') '.z',nint(z01)
    imgfn(length+6:120) = label//timsnd(1:tmstrln)
    length = 120
    CALL strlnth( imgfn, length)

    IF( nestgrd == 1 ) THEN
      WRITE(imgfn(length+1:length+3),'(''.G'',A)')                      &
            gridnum(mgrid:mgrid)

      length = length + 3
    END IF

    idata = iend-ibgn+1
    jdata = jend-jbgn+1

    CALL img2d(tem1,idata,1,idata,isk,jdata,1,jdata,jsk,                &
               bmax,bmin,                                               &
               imgfn(1:length))
!

  END IF

  RETURN
END SUBROUTINE img3d
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE IMG3D0                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE img3d0(b,nx,ibgn,iend,isk, ny,jbgn,jend,jsk,                 & 4,7
           nz,kbgn,kend,ksk,                                            &
           bmax,bmin, label,fnkey,time,mode,kslice,jslice,islice,       &
           mgrid,nestgrd, tem1)

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Produce 2-D HDF images for specified slices of 3-d array b.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!
!  MODIFICATION HISTORY:
!
!  9/1/94 (Y. Lu)
!  Cleaned up documentation
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    b        3-dimensional array of data
!
!    x        x coordinate of grid points in physical/comp. space (m)
!    y        y coordinate of grid points in physical/comp. space (m)
!    z        z coordinate of grid points in computational space (m)
!
!    nx       first dimension of b
!    ibgn     index of first i grid point to plot
!    iend     index of last  i grid point to plot
!             if isk >0, every isk'th data is used,
!             if isk <0, each grid zone is sub-divided into isk zone.
!
!    ny       second dimension of b
!    jbgn     index of first j grid point to plot
!    jend     index of last  j grid point to plot
!             if jsk >0, every jsk'th data is used,
!             if jsk <0, each grid zone is sub-divided into jsk zone.
!
!    nz       third dimension of b
!    kbgn     index  f first k grid point to plot
!    kend     index of last  k grid point to plot
!             if ksk >0, every ksk'th data is used,
!             if ksk <0, each grid zone is sub-divided into ksk zone.
!
!    label    character string describing the contents of b
!    fnkey    File name key
!    time     time of data in seconds
!
!    mode     slice orientation indicator
!             mode = 1, x-y slice at z index kslice is plotted.
!             mode = 2, x-z slice at y index jslice is plotted.
!             mode = 3, y-z slice at x index islice is plotted.
!             mode = 0, all of the three slices above are plotted.
!
!    kslice   k index of plane for mode=1 x-y slice
!    jslice   j index of plane for mode=2 x-z slice
!    islice   i index of plane for mode=1 y-z slice
!
!    mgrid    The grid number
!    nestgrd  Flag for nested grid run.
!
!  WORK ARRAYS:
!
!    tem1     Temporary work array.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz       ! 3 dimemsions of b
!  integer n              ! dimension on horizontal dirction

  REAL :: b(nx,ny,nz)       ! 3-dimensional array of data

  REAL :: bmax,bmin         ! the maximum and minimum value of b

  INTEGER :: ibgn,iend,isk  ! index of first and last i grid point;
                            ! skip/interpolation parameter in i-dir
  INTEGER :: jbgn,jend,jsk  ! index of first and last j grid point;
                            ! skip/interpolation parameter in j-dir
  INTEGER :: kbgn,kend,ksk  ! index of first and last k grid point;
                            ! skip/interpolation parameter in k-dir
  INTEGER :: length         ! the length of string

  INTEGER :: idata          ! total number of grid points in i-dir
  INTEGER :: jdata          ! total number of grid points in j-dir
  INTEGER :: kdata          ! total number of grid points in k-dir

  CHARACTER (LEN=*   ) :: label ! character string describing the contents
                                ! of b
  CHARACTER (LEN=80  ) :: timsnd  ! character string describing thetime of b
  CHARACTER (LEN=*   ) :: fnkey ! file name key
  INTEGER :: tmstrln        ! slice orientation indicator
  REAL :: time              ! time of data (s)

  INTEGER :: mode           ! slice orientation indicator
  INTEGER :: kslice         ! k index of plane for mode=1 x-y slice
  INTEGER :: jslice         ! j index of plane for mode=1 x-z slice
  INTEGER :: islice         ! i index of plane for mode=1 y-z slice

  INTEGER :: mgrid          ! the grid number
  INTEGER :: nestgrd        ! flag for nested grid run.

  REAL :: tem1(*)           ! temporary work array
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k,ij,ik,jk
  CHARACTER (LEN=120) :: imgfn
  CHARACTER (LEN=35) :: gridnum
  DATA gridnum /'123456789abcdefghijklmnopqrstuvwxyz'/
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!  mode=1   Plot x-y cross-section
!
!-----------------------------------------------------------------------
!
  CALL cvttsnd( time, timsnd, tmstrln )

  IF(mode == 1 .OR. mode == 0 ) THEN

    k = kslice
    DO j=jbgn,jend
      DO i=ibgn,iend
        ij = i-ibgn+1 + (j-jbgn)*(iend-ibgn+1)
        tem1(ij)=b(i,j,k)
      END DO
    END DO

    imgfn = fnkey
    length = LEN( fnkey )
    WRITE(imgfn(length+1:length+5),'(a,i3.3)') '.k',k
    imgfn(length+6:120) = label//timsnd(1:tmstrln)
    length = 120
    CALL strlnth( imgfn, length)

    IF( nestgrd == 1 ) THEN
      WRITE(imgfn(length+1:length+3),'(''.G'',A)')                      &
            gridnum(mgrid:mgrid)

      length = length + 3
    END IF

    idata = iend-ibgn+1
    jdata = jend-jbgn+1

    CALL img2d(tem1,idata,1,idata,isk,jdata,1,jdata,isk ,bmax,bmin,     &
               imgfn(1:length+4))
!
!-----------------------------------------------------------------------
!
!  mode=2   Plot x-z cross-section
!
!-----------------------------------------------------------------------
!
  ELSE IF (mode == 2 .OR. mode == 0 ) THEN

    j = jslice
    DO k=kbgn,kend
      DO i=ibgn,iend
        ik = i-ibgn+1 + (k-kbgn)*(iend-ibgn+1)
        tem1(ik)=b(i,j,k)
      END DO
    END DO

    imgfn = fnkey
    length = LEN( fnkey )
    WRITE(imgfn(length+1:length+5),'(a,i3.3)') '.j',j
    imgfn(length+6:120) = label//timsnd(1:tmstrln)
    length = 120
    CALL strlnth( imgfn, length)

    IF( nestgrd == 1 ) THEN
      WRITE(imgfn(length+1:length+3),'(''.G'',A)')                      &
            gridnum(mgrid:mgrid)

      length = length + 3
    END IF


    idata = iend-ibgn+1
    kdata = kend-kbgn+1

    CALL img2d(tem1,idata,1,idata,isk,kdata,1,kdata,ksk,                &
               bmax,bmin,                                               &
               imgfn(1:length+4))
!
!-----------------------------------------------------------------------
!
!  mode=3   Plot y-z cross-section
!
!-----------------------------------------------------------------------
!
  ELSE IF ( mode == 3 .OR. mode == 0) THEN

    i = islice
    DO k=kbgn,kend
      DO j=jbgn,jend
        jk = j-jbgn+1 + (k-kbgn)*(jend-jbgn+1)
        tem1(jk)=b(i,j,k)
      END DO
    END DO

    imgfn = fnkey
    length = LEN( fnkey )
    WRITE(imgfn(length+1:length+5),'(a,i3.3)') '.i',i
    imgfn(length+6:120) = label//timsnd(1:tmstrln)
    length = 120
    CALL strlnth( imgfn, length)

    IF( nestgrd == 1 ) THEN
      WRITE(imgfn(length+1:length+3),'(''.G'',A)')                      &
            gridnum(mgrid:mgrid)

      length = length + 3
    END IF

    jdata = jend-jbgn+1
    kdata = kend-kbgn+1

    CALL img2d(tem1,jdata,1,jdata,jsk,kdata,1,kdata,ksk,                &
               bmax,bmin,                                               &
               imgfn(1:length+4))

  ELSE IF( mode == 4 .OR. mode == 5 .OR. mode == 6 ) THEN

    WRITE(6,'(a)')                                                      &
        'Mode 4 or 5 is not available with IMG3D0, call IGM3D instead.'

  END IF

  RETURN
END SUBROUTINE img3d0
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE IMG2D                      ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE img2d(z,nx,ibgn,iend,isk,ny,jbgn,jend,jsk,                   & 9,3
           zmax,zmin,imgfn)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Create an HDF image file 'imgfn' for a 2-D field stored in z.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  4/14/93
!
!  MODIFICATION HISTORY:
!
!  9/1/94 (Y. Lu)
!  Cleaned up documentation
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    z        a two array containing the field for which a image file
!             is to be created.
!    nx, ny   the dimensions of array z.
!    ibgn     index of first i grid point to be used.
!    iend     index of last  i grid point to be used.
!    isk      skip or interpolation parameter in i-dir.
!             If isk >0, every isk'th data is used,
!             If isk <0, each grid zone is sub-divided into isk zones.
!
!    jbgn     index of first j grid point to be used.
!    jend     index of last  j grid point to be used.
!    jsk      skip or interpolation parameter in j-dir.
!             If jsk >0, every jsk'th data is used,
!             If jsk <0, each grid zone is sub-divided into jsk zones.
!
!    zmax,zmin
!             prespecified upper and lower bounds of the array value
!             to be used to scale the image.
!
!    imgfn    the name of the output image file.
!
!  OUTPUT:
!
!    None.
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!

  IMPLICIT NONE

  INTEGER :: nx, ny         ! dimensions of array z
  REAL :: z(nx,ny)       ! A two array containing the field for
                         ! which a image file is to be created.

  INTEGER :: ibgn,iend,isk  ! index of first and last i grid point;
                            ! skip/interpolation parameter in i-dir
  INTEGER :: jbgn,jend,jsk  ! index of first and last j grid point;
                            ! skip/interpolation parameter in j-dir

  REAL :: zmax,zmin      ! prespecified upper and lower bounds
                         ! of the array value to be used to scale
                         ! the image.

  CHARACTER (LEN=*) :: imgfn    ! name of the output image file


  INTEGER :: maxdim         ! maximum number of gird points in work
                            ! array image
  INTEGER :: icompres
  INTEGER :: idata          ! total number of grid points in i-dir
  INTEGER :: jdata          ! total number of grid points in j-dir

  INTEGER :: d8pimg, iret, imax,imin, nimage

  PARAMETER (maxdim= 512 , icompres=11)
  CHARACTER (LEN=1) :: image(maxdim*maxdim)

  REAL :: zx, zn
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  IF( isk > 0 )THEN
    idata = (iend-ibgn)/isk + 1
  ELSE
    idata = (iend-ibgn)*ABS(isk) + 1
  END IF

  IF( jsk > 0 )THEN
    jdata = (jend-jbgn)/jsk + 1
  ELSE
    jdata = (jend-jbgn)*ABS(jsk) + 1
  END IF


  IF( idata*jdata > maxdim*maxdim ) THEN
    WRITE(6,'(1x,a,a/)') 'Character work array image ',                 &
        'defined in IMG2D not big enough. Job stopped in IMG2D.'
    CALL arpsstop('arpsstop called from IMG2D due to character length',1)
  END IF
!
!-----------------------------------------------------------------------
!
!  RASTERIZE and store as HDF image
!
!-----------------------------------------------------------------------
!
  CALL rasteriz(z, nx,ibgn,iend,isk, ny,jbgn,jend,jsk,zmax,zmin,        &
                image,nimage,imax,imin)

!mx
  print*,'D8pimg is not called. No HDF image produced.'

  RETURN  ! D8pimg not found in HDF library for some reason

! iret = D8pimg (imgfn , image, idata,jdata, icompres)
!
!  IF (iret /= 0) THEN
!    WRITE(6,'(/1x,a,a,a,i3)')                                           &
!        'Error writing HDF file ',imgfn,'. Error flag was ',iret
!  END IF

!
!-----------------------------------------------------------------------
!
!  Find Max & Min (for output to screen only)
!
!-----------------------------------------------------------------------
!
  CALL a3dmax0(z, 1,nx,ibgn,iend, 1,ny,jbgn,jend,1,1,1,1, zx, zn)

  WRITE(6,'(/1x,3a)') 'Image file ',imgfn,' created.'
  WRITE(6,'(/1x,a,e15.5,a,e15.5,a,i4,a,i4/)')                           &
        'zmax=',zx,' zmin=',zn,' imin=',imin,' imax=',imax

!  100  FORMAT (1X,(a,2X),1P,2G16.7,0P, 2X,2I4,1X,2F5.2)

END SUBROUTINE img2d
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE RASTERIZ                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE rasteriz (z, nx,ibgn,iend,isk0, ny,jbgn,jend,jsk0,           & 1,1
           zmax,zmin, image, nimage,imax,imin)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Convert a 2-D real array into a raster form.
!
!  This routine should be equivalent to the C routine  floattor8.c,
!  p. 2.10 of NCSA HDF manual.
!
!  Values are scaled to be in the range 1 to 254 (Icharmin to Icharmax).
!  Image is a 1-D array and must be stored upside-down.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  4/14/93
!
!  MODIFICATION HISTORY:
!
!  9/1/94 (Y. Lu)
!  Cleaned up documentation
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    z          a two array containing the field for which a image
!               file is to be created.
!    nx, ny     the dimensions of array z.
!    zmax,zmin  the maximum or minimum value of z
!
!  OUTPUT:
!
!    image      1-D array
!    imax,imin  the maximum and minimum value of image
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx, ny         ! dimensions of array z
  REAL :: z(nx,ny)       ! A two array containing the field for
                         ! which a image file is to be created.

  INTEGER :: ibgn,iend,isk  ! index of first and last i grid point;
                            ! skip/interpolation parameter in i-dir
  INTEGER :: jbgn,jend,jsk  ! index of first and last j grid point;
                            ! skip/interpolation parameter in j-dir

  INTEGER :: isk0,jsk0
  REAL :: zmax, zmin     ! the maximum or minimum value of z

  INTEGER :: imax, imin, nimage
                             ! the maximum and minimum value and
! dimension of imag
  CHARACTER (LEN=1) :: image(nimage)
!
!-----------------------------------------------------------------------
!
!  Misc. variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: icharmax,icharmin,icp1,iblack
  PARAMETER (icharmax=254, icharmin=2, icp1=icharmax+1, iblack=255)

  INTEGER :: i,j,k, ivalue, rtoi, isub, jsub
  REAL :: step, rvalue ,arg, zj1,zj2
!
!-----------------------------------------------------------------------
!
!  Inline funtion
!
!-----------------------------------------------------------------------
!
  rtoi( arg )  = nint ((arg-zmin)/step) + icharmin
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  step      = (zmax - zmin) / icp1

  IF (step <= 0.) THEN
    PRINT *, '.....Error in RASTERIZ:  Zmax <= Zmin'
    RETURN
  END IF


  imax = -32000
  imin =  32000

  isk = isk0
  jsk = jsk0

  IF( isk*jsk < 0 ) THEN
    IF( isk < 0 .AND. jsk == 1 ) THEN
      jsk = -1
    ELSE IF( jsk < 0 .AND. isk == 1 ) THEN
      isk = -1
    ELSE
      WRITE(6,'(/1x,a,/1x,a,i2,a,i2,/1x,a/)')                           &
          'Error in RASTERIZ: isk and jsk must have the same sign.',    &
          'isk = ', isk, ' jsk= ', jsk,                                 &
          'Job stopped in the subroutine.'
      CALL arpsstop('arpsstop called from RASTERIZ due to sign error',1)
    END IF
  END IF


  IF( isk > 0 .AND. jsk > 0 ) THEN

    k = 0
    DO j=jend,jbgn,-jsk
      DO i=ibgn,iend,isk
        rvalue = z(i,j)
        ivalue = MIN( MAX( rtoi(rvalue),icharmin ), icharmax )

        imax   = MAX (imax, ivalue)
        imin   = MIN (imin, ivalue)

        k      = k + 1
        image(k)= CHAR (ivalue)
      END DO
    END DO

    nimage = k

  ELSE IF( isk < 0 .AND. jsk < 0 ) THEN
                                  ! Case with binlinear interpolation.

    k = 0
    DO j=jend,jbgn+1,-1
      DO jsub=1,ABS(jsk)

        DO i=ibgn,iend-1
          DO isub=1,ABS(isk)
            zj1=z(i,j  )+(isub-1.0)/ABS(isk)*(z(i+1,j  )-z(i,j  ))
            zj2=z(i,j-1)+(isub-1.0)/ABS(isk)*(z(i+1,j-1)-z(i,j-1))
            rvalue = zj1+(jsub-1.0)/ABS(jsk)*(zj2-zj1)

            ivalue = MIN( MAX( rtoi(rvalue),icharmin ), icharmax )
            imax   = MAX (imax, ivalue)
            imin   = MIN (imin, ivalue)
            k      = k + 1
            image(k)= CHAR (ivalue)
          END DO
        END DO
        rvalue = z(iend,j)+(jsub-1.0)/ABS(jsk) *                        &
                (z(iend,j-1)-z(iend,j))
        ivalue = MIN( MAX( rtoi(rvalue),icharmin ), icharmax )
        imax   = MAX (imax, ivalue)
        imin   = MIN (imin, ivalue)
        k      = k + 1
        image(k)= CHAR (ivalue)
      END DO
    END DO

    DO i=ibgn,iend-1
      DO isub=1,ABS(isk)
        zj1=z(i,jbgn)+(isub-1.0)/ABS(isk)*(z(i+1,jbgn)-z(i,jbgn))
        rvalue = zj1

        ivalue = MIN( MAX( rtoi(rvalue),icharmin ), icharmax )
        imax   = MAX (imax, ivalue)
        imin   = MIN (imin, ivalue)
        k      = k + 1
        image(k)= CHAR (ivalue)
      END DO
    END DO

    rvalue = z(iend,jbgn)
    ivalue = MIN( MAX( rtoi(rvalue),icharmin ), icharmax )
    imax   = MAX (imax, ivalue)
    imin   = MIN (imin, ivalue)
    k      = k + 1
    image(k)= CHAR (ivalue)

    nimage = k

  END IF

  RETURN
END SUBROUTINE rasteriz