PROGRAM pltradcol,7
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 PROGRAM PLTRADCOL                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!  PURPOSE:
!
!  Plots data written by wtradcol.
!
!
!  AUTHOR:
!
!  Keith Brewster, CAPS, April, 1996
!
!  LINKING:
!
!  ncargf77 -o pltradcol pltradcol.f pltmap.f maproj3d.f timelib3d.f
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'
!
  INTEGER :: pltref,pltvel,pltnyq,plttim
  PARAMETER (pltref=1,                                                  &
             pltvel=1,                                                  &
             pltnyq=0,                                                  &
             plttim=0)
  CHARACTER (LEN=4) :: radid
  CHARACTER (LEN=5) :: lvlid
  INTEGER :: ireftim,itime
  INTEGER :: vcpnum
  INTEGER :: mxradvr
  PARAMETER (mxradvr=10)
  INTEGER :: iradvr
  INTEGER :: nradvr
!-----------------------------------------------------------------------
!
!  Dimension in x, y, and z direction
!
!-----------------------------------------------------------------------
  INTEGER :: nx, ny, nz
!-----------------------------------------------------------------------
!
!  Read-in variables
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: readk(:)
  REAL, ALLOCATABLE :: readhgt(:)
  REAL, ALLOCATABLE :: readref(:)
  REAL, ALLOCATABLE :: readvel(:)
  REAL, ALLOCATABLE :: readnyq(:)
  REAL, ALLOCATABLE :: readtim(:)
!
!-----------------------------------------------------------------------
!
!  Plotting variables
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: gridref(:,:,:)
  REAL, ALLOCATABLE :: gridvel(:,:,:)
  REAL, ALLOCATABLE :: gridnyq(:,:,:)
  REAL, ALLOCATABLE :: gridtim(:,:,:)
!
!-----------------------------------------------------------------------
!
!  Statistics variables
!
!-----------------------------------------------------------------------
!
  INTEGER, ALLOCATABLE :: kntref(:)
  INTEGER, ALLOCATABLE :: kntvel(:)
!
  REAL :: refint,velint,nyqint,timint
  PARAMETER (refint=0.,                                                 &
             velint=0.,                                                 &
             nyqint=0.,                                                 &
             timint=0.)
!
!-----------------------------------------------------------------------
!
!  Map orientation variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: maxpts
  PARAMETER (maxpts = 100)
  REAL :: latmap(maxpts),lonmap(maxpts)
  REAL :: xmap(maxpts),ymap(maxpts)
!
!-----------------------------------------------------------------------
!
!  Misc local variables
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=3) :: chplt
  CHARACTER (LEN=18) :: timplt
  CHARACTER (LEN=72) :: fname
  CHARACTER (LEN=72) :: mapfile
  PARAMETER                                                             &
      (mapfile='spcounty.mapdata')
  REAL :: rdummy
  REAL :: x,y,xrd,yrd,gridlat,gridlon,elev
  REAL :: latnot(2)
  REAL :: ctrx,ctry,swx,swy,nex,ney
  INTEGER :: i,j,k,kk,klev,ipt, munit, mapgrid
  INTEGER :: idummy
  INTEGER :: iyr, imon, idy, ihr, imin, isec
  INTEGER :: istatus
!-----------------------------------------------------------------------
!
!  NAMELIST declaration
!
!-----------------------------------------------------------------------
  NAMELIST /grid_dims/ nx, ny, nz
  NAMELIST /file_name/ fname

!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Execute code begin
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!-----------------------------------------------------------------------
!
!  Read in the dimensions
!
!-----------------------------------------------------------------------
  READ(5,grid_dims,END=100)
  WRITE(6,'(/a)') 'Namelist block grid_dims successfully read.'
  
!-----------------------------------------------------------------------
!
!  Allocate and Initializations
!
!-----------------------------------------------------------------------
!
  ALLOCATE(readk(nz) ,STAT=istatus)
  readk=0
  ALLOCATE(readhgt(nz),STAT=istatus)
  readhgt=0
  ALLOCATE(readref(nz),STAT=istatus)
  readref=0
  ALLOCATE(readvel(nz),STAT=istatus)
  readvel=0
  ALLOCATE(readnyq(nz),STAT=istatus)
  readnyq=0
  ALLOCATE(readtim(nz),STAT=istatus)
  readtim=0
  
  ALLOCATE(gridref(nx,ny,nz),STAT=istatus)
  gridref=-999999.
  ALLOCATE(gridvel(nx,ny,nz),STAT=istatus)
  gridvel=-999999.
  ALLOCATE(gridnyq(nx,ny,nz),STAT=istatus)
  gridnyq=-999999.
  ALLOCATE(gridtim(nx,ny,nz),STAT=istatus)
  gridtim=-999999.
  
  ALLOCATE(kntref(nz),STAT=istatus)
  kntref=0
  ALLOCATE(kntvel(nz),STAT=istatus)
  kntvel=0
  
!
!-----------------------------------------------------------------------
!
!  Get user input
!
!-----------------------------------------------------------------------
!
  READ(5,file_name,END=100)
  WRITE(6,'(/a)') 'Namelist block file_name successfully read.'

  OPEN(11,FILE=fname,FORM='unformatted',STATUS='old')
!
!-----------------------------------------------------------------------
!
!  Read radar columns
!
!-----------------------------------------------------------------------
!
  READ(11) radid
  READ(11) ireftim,itime,vcpnum,idummy,idummy,                          &
           idummy,idummy,idummy,idummy,idummy
!
  READ(11) runname
  READ(11) hdmpfmt,strhopt,mapproj,idummy,idummy,                       &
           idummy,idummy,idummy,idummy,idummy
  READ(11) dx,dy,dz,dzmin,ctrlat,                                       &
           ctrlon,trulat1,trulat2,trulon,sclfct,                        &
           rdummy,rdummy,rdummy,rdummy,rdummy
  READ(11) nradvr,iradvr
!
  DO ipt=1,(nx*ny)
    READ(11,END=51) i,j,xrd,yrd,                                        &
                   gridlat,gridlon,elev,klev
    READ(11,END=52) (readk(kk),kk=1,klev)
    READ(11,END=52) (readhgt(kk),kk=1,klev)
    READ(11,END=52) (readref(kk),kk=1,klev)
    READ(11,END=52) (readvel(kk),kk=1,klev)
    READ(11,END=52) (readnyq(kk),kk=1,klev)
    READ(11,END=52) (readtim(kk),kk=1,klev)
    DO kk=1,klev
      k=nint(readk(kk))
      gridref(i,j,k)=readref(kk)
      gridvel(i,j,k)=readvel(kk)
      gridnyq(i,j,k)=readnyq(kk)
      gridtim(i,j,k)=readtim(kk)
      IF(gridref(i,j,k) > -200. .AND.                                   &
          gridref(i,j,k) < 200.) kntref(k)=kntref(k)+1
      IF(gridvel(i,j,k) > -200. .AND.                                   &
          gridvel(i,j,k) < 200.) kntvel(k)=kntvel(k)+1
    END DO
  END DO
  51 CONTINUE
  ipt=ipt-1
  WRITE(6,'(a,i6,a)') ' End of file reached after reading',             &
                     ipt,' columns'
  GO TO 55
  52 CONTINUE
  WRITE(6,'(a,i6,a)') ' End of file reached while reading',             &
                     ipt,' column'
  55 CONTINUE
  CLOSE(11)
!
!-----------------------------------------------------------------------
!
!  Write statistics
!
!-----------------------------------------------------------------------
!
  WRITE(6,'(a)') '  k   n vel    n vel'
  DO k=1,nz
    WRITE(6,'(i3,2i10)') k,kntref(k),kntvel(k)
  END DO
!
  CALL abss2ctim(itime, iyr, imon, idy, ihr, imin, isec )
  iyr=MOD(iyr,100)
  WRITE(timplt,815) imon,idy,iyr,ihr,imin
  815 FORMAT(i2.2,'/',i2.2,'/',i2.2,1X,i2.2,':',i2.2,' UTC')
!
  CALL opngks
!
!-----------------------------------------------------------------------
!
!  Set up map
!
!-----------------------------------------------------------------------
!
  IF(ctrlon > 180.) ctrlon=ctrlon-360.
  IF(trulon > 180.) trulon=trulon-360.
  latnot(1)=trulat1
  latnot(2)=trulat2
  CALL setmapr(mapproj,sclfct,latnot,trulon)
  CALL lltoxy(1,1,ctrlat,ctrlon,ctrx,ctry)
  swx=ctrx-((nx-3)/2)*dx
  swy=ctry-((ny-3)/2)*dy
  nex=swx+(nx-3)*dx
  ney=swy+(ny-3)*dy
!
!-----------------------------------------------------------------------
!
!  Vertical level loop
!
!-----------------------------------------------------------------------
!
  DO k=10,nz,2
    PRINT *, ' plotting k= ',k
!    CALL CONREC(gridref(1,1,k),
!    +              nx,nx,ny,0.,0.,refint,-1,-1,-127)
!
    WRITE(lvlid,'(a,i3)') 'k=',k
!
!-----------------------------------------------------------------------
!
!  Plot reflectivity
!
!-----------------------------------------------------------------------
!
    IF( pltref > 0 .AND. kntref(k) > 0 ) THEN
      CALL set( 0.,1.0,0.0,1.0,                                         &
             -.1,1.1,-.1,1.1,1)
      CALL wtstr(0.05,1.05,radid,16,0,0)
      CALL wtstr(0.2,1.05,lvlid,16,0,0)
      CALL wtstr(0.6,1.05,timplt,16,0,0)
      CALL wtstr(0.9,1.05,'Reflectivity',16,0,0)
      CALL set( 0.,1.0,0.0,1.0,                                         &
             swx,nex,swy,ney,1)
      DO j=1,ny
        DO i=1,nx
          IF(gridref(i,j,k) > -200. .AND. gridref(i,j,k) < 200.) THEN
            WRITE(chplt,820) nint(gridref(i,j,k))
            820         FORMAT(i3)
            x=swx+FLOAT(i-1)*dx
            y=swy+FLOAT(j-1)*dy
            CALL wtstr(x,y,chplt,8,0,0)
          END IF
        END DO
      END DO
      CALL pltmap(maxpts,mapfile,latmap,lonmap,xmap,ymap)
      CALL frame
    END IF
!
!-----------------------------------------------------------------------
!
!  Plot velocity
!
!-----------------------------------------------------------------------
!
    IF( pltvel > 0 .AND. kntvel(k) > 0 ) THEN
      CALL set( 0.,1.0,0.0,1.0,                                         &
           -.1,1.1,-.1,1.1,1)
      CALL wtstr(0.1,1.05,radid,16,0,0)
      CALL wtstr(0.2,1.05,lvlid,16,0,0)
      CALL wtstr(0.4,1.05,timplt,16,0,0)
      CALL wtstr(0.8,1.05,'Radial Velocity',16,0,0)
      CALL set( 0.,1.0,0.0,1.0,                                         &
               swx,nex,swy,ney,1)
      DO j=1,ny
        DO i=1,nx
          IF(gridvel(i,j,k) > -200. .AND. gridvel(i,j,k) < 200.) THEN
            WRITE(chplt,820) nint(gridvel(i,j,k))
            x=swx+FLOAT(i-1)*dx
            y=swy+FLOAT(j-1)*dy
            CALL wtstr(x,y,chplt,8,0,0)
          END IF
        END DO
      END DO
!      CALL CONREC(gridvel(1,1,k),
!    +                nx,nx,ny,0.,0.,velint,-1,-1,-127)
      CALL pltmap(maxpts,mapfile,latmap,lonmap,xmap,ymap)
      CALL frame
    END IF
!
!-----------------------------------------------------------------------
!
!  Plot Nyquist velocity
!
!-----------------------------------------------------------------------
!
    IF(pltnyq > 0) THEN
      CALL set( 0.,1.0,0.0,1.0,                                         &
             -.1,1.1,-.1,1.1,1)
      CALL wtstr(0.1,1.05,radid,16,0,0)
      CALL wtstr(0.2,1.05,lvlid,16,0,0)
      CALL wtstr(0.4,1.05,timplt,16,0,0)
      CALL wtstr(0.8,1.05,'Nyquist',16,0,0)
      CALL set( 0.,1.0,0.0,1.0,                                         &
               swx,nex,swy,ney,1)
      DO j=1,ny
        DO i=1,nx
          IF(gridnyq(i,j,k) > 0. .AND. gridnyq(i,j,k) < 1000.) THEN
            WRITE(chplt,820) nint(10.*gridnyq(i,j,k))
            x=swx+FLOAT(i-1)*dx
            y=swy+FLOAT(j-1)*dy
            CALL wtstr(x,y,chplt,8,0,0)
          END IF
        END DO
      END DO
!      CALL CONREC(gridnyq(1,1,k),
!    +                nx,nx,ny,0.,0.,nyqint,-1,-1,-127)
      CALL pltmap(maxpts,mapfile,latmap,lonmap,xmap,ymap)
      CALL frame
    END IF
!
!-----------------------------------------------------------------------
!
!  Plot data time
!
!-----------------------------------------------------------------------
!
    IF(plttim > 0) THEN
      CALL set( 0.,1.0,0.0,1.0,                                         &
             -.1,1.1,-.1,1.1,1)
      CALL wtstr(0.1,1.05,radid,16,0,0)
      CALL wtstr(0.2,1.05,lvlid,16,0,0)
      CALL wtstr(0.4,1.05,timplt,16,0,0)
      CALL wtstr(0.8,1.05,'Time ',16,0,0)
      CALL set( 0.,1.0,0.0,1.0,                                         &
               swx,nex,swy,ney,1)
      DO j=1,ny
        DO i=1,nx
          IF(gridtim(i,j,k) > -1000. .AND. gridtim(i,j,k) < 1000.) THEN
            WRITE(chplt,820) nint(0.1*gridtim(i,j,k))
            x=swx+FLOAT(i-1)*dx
            y=swy+FLOAT(j-1)*dy
            CALL wtstr(x,y,chplt,8,0,0)
          END IF
        END DO
      END DO
!      CALL CONREC(gridtim(1,1,k),
!    +                nx,nx,ny,0.,0.,timint,-1,-1,-127)
      CALL pltmap(maxpts,mapfile,latmap,lonmap,xmap,ymap)
      CALL frame
    END IF
!
  END DO
!
  CALL clsgks
!
  GOTO 101
  
  100 CONTINUE
  
  WRITE(6,'(/a,a)') 'Error reading NAMELIST file. The program will abort.'
  
  101 CONTINUE  
 
  STOP
END PROGRAM pltradcol