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.
!
! 03/03/2003 (K. Brewster)
!
! Modified scaling of station markers and name sizes.
!
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,chsize,yoffset
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))
chsize=sta_marksz(itype)*yl*7.
yoffset=sta_marksz(itype)*yl*10.
IF( ovrstam > 0) THEN
print *, 'plotting marker'
CALL xmrksz(sta_marksz(itype))
CALL xmarker(xstn,ystn,sta_marktyp(itype))
END IF
IF( ovrstan > 0) THEN
CALL xchsiz(chsize)
length = 5
CALL strlnth(s_name(i),length)
IF( ovrstam > 0) THEN
CALL xcharc(xstn, (ystn-yoffset), 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 ) 235
!
!-----------------------------------------------------------------------
!
! 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)
length = LEN_TRIM(string)
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