PROGRAM  pltgrid
!
!  Program to plot model grids
!
!  Commands:
!
!    zxncarf77 -o pltgrid pltgrid.f
!    pltgrid < pltgrid.input
!
!  or in ARPS root directory, enter
!
!    makearps pltgrid
!    bin/pltgrid  < input/pltgrid.input
!
! MODIFICATION HISTORY:
!
! 11/27/2001 (K. Brewster)
!
! Added plotting of station locations to plot, correcting some 
! previous test code.
!
  IMPLICIT NONE

  INTEGER :: nx,ny
  INTEGER :: mapproj
  REAL :: dx,dy,ctrlat,ctrlon,trulat1,trulat2,trulon,sclfct,xl,yl

  INTEGER :: nnwgrd, nnwgrdmax
  PARAMETER (nnwgrdmax = 10) ! maximum number of new grids

  INTEGER :: nx1(nnwgrdmax),ny1(nnwgrdmax)
  REAL :: xctr1(nnwgrdmax) , yctr1(nnwgrdmax)
  REAL :: dx1(nnwgrdmax),dy1(nnwgrdmax)

  NAMELIST /grid/ nx,ny,dx,dy,ctrlat,ctrlon

  NAMELIST /projection/ mapproj,trulat1,trulat2,trulon,sclfct

  NAMELIST /newgrid/ nnwgrd,nx1,ny1,xctr1,yctr1,dx1,dy1

  INTEGER :: ovrmap,mapgrid,nmapfile,lmapfile,i, mxmapfile
  PARAMETER (mxmapfile= 10)
  CHARACTER (LEN=132) :: mapfile(mxmapfile)

  REAL :: latgrid,longrid, xorig, yorig
  COMMON /mappar2/ latgrid,longrid
  NAMELIST /map_plot/ ovrmap,mapgrid,latgrid,longrid,nmapfile,          &
         mapfile

  INTEGER :: pltstn
  CHARACTER (LEN=100) :: stnfile
  CHARACTER (LEN=2) :: statelist(60)
  NAMELIST /station_plot/ pltstn, stnfile, statelist

  CHARACTER (LEN=20) :: stnname
  CHARACTER (LEN=2) :: state
  CHARACTER (LEN=4) :: stnsymbol
  REAL :: lat,long, x1,x2,y1,y2
  INTEGER :: elv,TYPE,length
  REAL :: xstn,ystn

  INTEGER :: mxstalo   ! maximum number of observation stations
  PARAMETER (mxstalo=500)

  INTEGER :: ovrstaopt
  INTEGER :: ovrstam,staset,ovrstan,ovrstav,stacol,markprio,wrtstax
  INTEGER :: nsta_typ,sta_typ(10),sta_marktyp(10),sta_markcol(10)
  REAL :: sta_marksz(10)
  REAL :: wrtstad

  INTEGER :: nsta,nstapro(mxstalo),nstatyp(mxstalo)
  REAL :: latsta(mxstalo), lonsta(mxstalo)
  CHARACTER (LEN=2) :: s_state(mxstalo)
  CHARACTER (LEN=5) :: s_name(mxstalo)
  CHARACTER (LEN=20) :: s_site(mxstalo)
  INTEGER :: s_elev(mxstalo)

  CHARACTER (LEN=132) :: stalofl
  INTEGER :: lstalofl

  NAMELIST /plot_sta/ ovrstaopt,ovrstam,ovrstan,ovrstav,wrtstax,        &
      wrtstad, stacol, markprio, nsta_typ, sta_typ, sta_marktyp,        &
      sta_markcol,sta_marksz,stalofl


  REAL :: swx,swy,ctrx,ctry
  LOGICAL :: fexist
  REAL :: truelat(2)
  REAL :: xl1, yl1, ctrlat1, ctrlon1,ysize
  REAL :: cornerlat,cornerlon
  INTEGER :: itype

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!  Input control parameters map plotting
!
!-----------------------------------------------------------------------
!
  WRITE(6,'(/a/)')                                                      &
      ' Please input control parameters for map plotting.'

  READ(5,grid,ERR=100)

  WRITE(6,'(/5x,a,i8)')'Input nx was',nx
  WRITE(6,'(/5x,a,i8)')'Input ny was',ny
  WRITE(6,'(/5x,a,f12.6,a)')'Input dx was',dx,' meters'
  WRITE(6,'(/5x,a,f12.6,a)')'Input dy was',dy,' meters'


  WRITE(6,'(/5x,a,f12.6,a)')                                            &
       'Input ctrlat was ',ctrlat,' degree North'

  WRITE(6,'(/5x,a,f12.6,a)')                                            &
       'Input ctrlon was ',ctrlon,' degree East'

!-----------------------------------------------------------------------
!
!  Input map projection parameters
!
!-----------------------------------------------------------------------
!
  READ (5,projection,ERR=100)

  WRITE(6,'(/5x,a,i4)') 'Input mapproj was ',mapproj

  WRITE(6,'(/5x,a,f12.6,a)')                                            &
       'Input trulat1 was ',trulat1,' degree North'

  WRITE(6,'(/5x,a,f12.6,a)')                                            &
       'Input trulat2 was ',trulat2,' degree North'

  WRITE(6,'(/5x,a,f12.6)')                                              &
      'The latitude of the center of the model domain was ',ctrlat

  WRITE(6,'(/5x,a,f12.6)')                                              &
      'The longitude of the center of the model domain was ',ctrlon

  WRITE(6,'(/5x,a,f12.6,a)')                                            &
       'Input trulon was ',trulon,' degree East'

  IF ( mapproj == 0 ) THEN
    trulat1 = ctrlat
    trulat2 = ctrlat
    trulon  = ctrlon
  END IF
!
!-----------------------------------------------------------------------
!
!  Set the output grid and the variable control parameters
!
!-----------------------------------------------------------------------
!
  READ (5,newgrid)

  PRINT*,' Input nnwgrd=',nnwgrd

  IF( nnwgrd > nnwgrdmax)  THEN
    PRINT*,'Number of new grids more than maximum allowed.'
    PRINT*,'Increase the size of nnwgrdmax in the program.'
    STOP
  END IF

  PRINT*,' Input nx1=',(nx1(i),i=1,nnwgrd)
  PRINT*,' Input ny1=',(ny1(i),i=1,nnwgrd)
  PRINT*,' Input dx1=',(dx1(i),i=1,nnwgrd)
  PRINT*,' Input dy1=',(dy1(i),i=1,nnwgrd)
  PRINT*,' Input xctr1=',(xctr1(i),i=1,nnwgrd)
  PRINT*,' Input yctr1=',(yctr1(i),i=1,nnwgrd)

  READ(5,map_plot,ERR=100)

! READ(5,station_plot,ERR=100)

  GO TO 10
  100   WRITE(6,'(a)') 'Error reading NAMELIST file. Program stopped.'
  STOP

  10    CONTINUE

!
!-----------------------------------------------------------------------
!
!  Initialize ZXPLOT plotting package
!
!-----------------------------------------------------------------------
!
  CALL xdevic
  CALL xafstyl(1)
  CALL xsetclrs(1)
  CALL xcolor(1)
  CALL xartyp(2)

  CALL xdspac(0.9)

  CALL xaxfmt( '(i10)' )

!-----------------------------------------------------------------------
!
!  Set up map projection.
!
!-----------------------------------------------------------------------
!
  xl = (nx-3)*dx * 0.001
  yl = (ny-3)*dy * 0.001
  xorig = 0.0
  yorig = 0.0

  CALL xstpjgrd(mapproj,trulat1,trulat2,trulon,                         &
                ctrlat,ctrlon,xl,yl,xorig,yorig)

  CALL xxytoll(1,1,xorig*1000.0,yorig*1000.0,cornerlat,cornerlon)

  WRITE(6,'(/1x,a,2f12.6)')                                             &
      'Lat/lon at the SW corner of base grid=',cornerlat,cornerlon

  CALL xxytoll(1,1,(xorig+xl)*1000.0,(yorig+yl)*1000.0,                 &
               cornerlat,cornerlon)
  WRITE(6,'(1x,a,2f12.6/)')                                             &
      'Lat/lon at the NE corner of base grid=',cornerlat,cornerlon

  IF( xl >= yl) THEN
    CALL xpspac(0.1,0.9, 0.5-yl/xl*0.4,0.5+yl/xl*0.4)
  ELSE
    CALL xpspac(0.5-xl/yl*0.4,0.5+xl/yl*0.4,0.1, 0.9)
  END IF
!
  CALL xmap(0.0,xl, 0.0, yl)
  CALL xaxsca(0.0,xl, dx/1000.0, 0.0, yl, dy/1000.0 )

  CALL xcolor(9)

  IF( ovrmap /= 0 ) THEN

    DO i=1,MIN(mxmapfile,nmapfile)

      lmapfile=132
      CALL strlnth(mapfile(i), lmapfile)
      WRITE(6,'(1x,a,a)') 'Input was ',mapfile(i)(1:lmapfile)

      INQUIRE(FILE=mapfile(i)(1:lmapfile), EXIST = fexist )
      IF( .NOT.fexist) THEN
        WRITE(6,'(a)') 'Map file '//mapfile(i)(1:lmapfile)//            &
            ' not found. Corresponding map no plotted.'
      ELSE
        CALL xdrawmap(10,mapfile(i)(1:lmapfile),latgrid,longrid)
      END IF

    END DO

  END IF

  CALL xthick(3)
  DO i=1,nnwgrd

    CALL xcolor(5+i)
    xl1 = (nx1(i)-3)*dx1(i) * 0.001
    yl1 = (ny1(i)-3)*dy1(i) * 0.001
    CALL xbox(xctr1(i)/1000.0-xl1/2, xctr1(i)/1000.0+xl1/2,             &
              yctr1(i)/1000.0-yl1/2, yctr1(i)/1000.0+yl1/2)

    CALL xxytoll(1,1,xctr1(i),yctr1(i),ctrlat1,ctrlon1)
    WRITE(6,'(/1x,a,2f15.4,2f12.6)')                                    &
           'xctr1,yctr1,ctrlat1,ctrlon1=',                              &
            xctr1(i),yctr1(i),ctrlat1,ctrlon1

    CALL xxytoll(1,1,xctr1(i)-xl1*500.0,yctr1(i)-yl1*500.0,             &
         cornerlat,cornerlon)
    WRITE(6,'(1x,a,i3,a,2f12.6)')                                       &
           'Lat/lon at the SW corner of new grid no.',i,                &
           '=', cornerlat,cornerlon

    CALL xxytoll(1,1,xctr1(i)+xl1*500.0,yctr1(i)+yl1*500.0,             &
         cornerlat,cornerlon)
    WRITE(6,'(1x,a,i3,a,2f12.6/)')                                      &
          'Lat/lon at the NE corner of new grid no.',i,                 &
          '=', cornerlat,cornerlon

  END DO

  ovrstam=0
  ovrstan=0
  ovrstav=0
  wrtstax=0

  READ(5,plot_sta, ERR=72)

  WRITE(6,'(a)')                                                        &
      'Namelist plot_sta was successfully read.'

  lstalofl=LEN(stalofl)
  CALL strlnth(stalofl,lstalofl)
  WRITE(6,'(1x,a,a)') 'Station file name was ',stalofl(1:lstalofl)

  72 CONTINUE

  nsta=0

  IF(ovrstaopt /= 0) THEN
    INQUIRE(FILE=stalofl(1:lstalofl), EXIST = fexist )
    IF( .NOT.fexist) THEN
      WRITE(6,'(a)') 'WARNING: Surface obs file '                       &
          //stalofl(1:lstalofl)//                                       &
          ' not found. Program continue in ARPSPLT.'
    ELSE
      CALL read_station(stalofl,mxstalo,latsta,lonsta,nstatyp,          &
                  nstapro,nsta,s_name,s_state,s_site,s_elev)
      IF(nsta > 0) staset=1
    END IF

    print *, ' Done reading, nsta= ',nsta
    CALL xcolor(1)

    DO i=1,nsta
      CALL xlltoxy(1,1,latsta(i),lonsta(i), xstn,ystn)
      print *, ' Plotting ',s_name(i),latsta(i),lonsta(i),.001*xstn,.001*ystn
      xstn=xstn*0.001
      ystn=ystn*0.001
      DO itype=1,nsta_typ-1
        IF(nstatyp(i) == sta_typ(itype)) EXIT
      END DO
      print *, ' nstatyp, itype = ',nstatyp(i),itype
      call xcolor(sta_markcol(itype))
      ysize=sta_marksz(itype)*yl
      IF( ovrstam > 0) THEN
        print *, 'plotting marker'
        CALL xmrksz(ysize*100.)
        CALL xmarker(xstn,ystn,sta_marktyp(itype))
      END IF
      IF( ovrstan > 0) THEN
        CALL xchsiz(ysize)
        length = 5
        CALL strlnth(s_name(i),length)
        IF( ovrstam > 0) THEN
          CALL xcharc(xstn, (ystn-ysize), s_name(i)(1:length))
        ELSE
          CALL xcharc(xstn, ystn, s_name(i)(1:length))
        END IF
      END IF
    END DO
  END IF

  CALL xframe
  CALL xgrend

  STOP
END PROGRAM  pltgrid

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


SUBROUTINE strlnth( string, length ) 230
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Return the length of the non-blank part of a character string.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  11/20/91
!
!  MODIFICATION HISTORY:
!
!  5/05/92 (M. Xue)
!  Added full documentation.
!
!  9/10/94 (Weygandt & Y. Lu)
!  Cleaned up documentation.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    string   A character string
!    length   The declared length of the character string 'string'.
!
!  OUTPUT:
!
!    length   The length of the non-blank part of the string.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  CHARACTER (LEN=*     ) :: string ! A character string for the name of
                                   ! this run.
  INTEGER :: length            ! The length of the non-blank part
                               ! of a string.

  INTEGER :: i
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  DO i = length,1,-1

    IF(string(i:i) /= ' ') EXIT

  END DO

  200   CONTINUE

  length = MAX(1,i)

  RETURN
END SUBROUTINE strlnth

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


SUBROUTINE read_station(infile,mxstalo,latsta,lonsta,                   & 2
           nstatyp,nstapro,nsta,sname,state,sitena,nelev)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  This subroutine will read external staion information.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  AUTHOR:
!    Min Zou (6/1/97)
!
!  Modification history:
!
!-----------------------------------------------------------------------
!
!  Variable Declarations
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=*) :: infile
  INTEGER :: nsta,nstapro(mxstalo),nstatyp(mxstalo)
  REAL :: latsta(mxstalo), lonsta(mxstalo)
  CHARACTER (LEN=5) :: sname(mxstalo)
  CHARACTER (LEN=2) :: state(mxstalo)
  CHARACTER (LEN=20) :: sitena(mxstalo)
  INTEGER :: nelev(mxstalo)
  CHARACTER (LEN=132) :: line

  OPEN(1,IOSTAT=ios,FILE=infile,STATUS='old',FORM='formatted')
  IF(ios /= 0) THEN     ! error during read
    istatus = -1
    WRITE(6,650) infile
    650       FORMAT(' +++ ERROR opening: ',a70,' +++')
    WRITE(6,651) ios
    651       FORMAT('     IOS code = ',i5)
    RETURN
  END IF

  nsta = 0

! Read only lines that begin with A-Z, a-z, or 0-9 -- treat the rest as
! comments

  DO i=1,mxstalo
    READ(1,'(A)',END=999) line
    IF ( ( line(:1) >= 'A' .AND. line(:1) <= 'Z' ) .OR.                 &
           ( line(:1) >= 'a' .AND. line(:1) <= 'z' ) .OR.               &
           ( line(:1) >= '0' .AND. line(:1) <= '9' ) ) THEN
      nsta = nsta + 1
      READ(line,101,ERR=999)sname(nsta),state(nsta),sitena(nsta),       &
          latsta(nsta),lonsta(nsta),nelev(nsta),nstatyp(nsta),          &
          nstapro(nsta)
    END IF
    101   FORMAT(a5,2X,a2,1X,a20,4X,f8.3,1X,f8.3,1X,i5,2X,i2,i1)
    102   FORMAT(a5,2X,a2,1X,a20,4X,f8.3,1X,f8.3,1X,i5,2X,i2,1X,i1)
  END DO
  999   CONTINUE
  CLOSE(1)

  RETURN
END SUBROUTINE read_station