c	this is the latest version edited for f90 and lblbr: 17 Sep 2015
c	Track colour order: Black, White, Dark Blue, Red, Orange

c	modified on 22 Jun 2006 for five tracks
c	14 Oct 2014: output ps/pdf file can be generated. Lenth of output file character fixed. monthly plots possible. change mnthplt in ngrd sub
c       19 sep 2015: for no label bar for marking plot. search lblbr
c	modify the b/w routine for many of the changes ref y/n commands
c	keywords: 
c	for keywords search for conditions
c	mhtb top/bottom heading; 
c	mrk for markings
c	lblbr for label bar
c	spl label sizes are fixed for 6 chars if more are required change the character *
c	clrt for colour labels

c	special commands:
c	for superscripts: watts/m:S:2:N: (for wattrs/m**2)
c	for ':'  Fig. 1:::N:Figure caption
      	real zreg(1000,500),xreg(1000),yreg(1000),gx(1000),gy(1000),sumw
     1(1000,1000),sumz(1000,1000)
	print '(/a,$)',' No. of grids for interpolation in  x & y direction
     1s: '
	read *,mreg,nreg
	call total(mreg,nreg,zreg,xreg,yreg,gx,gy,sumw,sumz)
	print *, 'Track colour order: Black, White, Dark Blue, 
	1Red, Orange'
	stop	
	end

	subroutine total(mreg,nreg,zreg,xreg,yreg,gx,gy,sumw,sumz)
c	to have superscript (eg. square) put :S:2:N:
c	column numbers of subroutines:
c	colcon-89 ccpscm-316, fill-581, sfill-621, msk-641, getdat-671, color-731, clr-862, 
c	grd-878, ngrd-900, landfill-1022, colram-1041, colrln-1058
c	info: if mx and min values are zero, defaults obtained from
c	the data set are used; otherwise the max and min specified values
c	are used. if a value is greater than mx, same clr is put for >mx
c	but better to consider the actual mx value while giving mx/mn
c	values.
c	does not plot x/y axes at all places if more than one plot is chosen
c	in case it is required,comment xaxt, yaxt 
c	imp: if max value of the contour is given more than the available value
c	white colour appears in the figure. hence max value should 
c	be .le. available max value. similarly for min value!
c	if landfill is no and boundary is y, plots only land boundary

c	modified on 6 nov 97	
c	modified on 17 Feb 98 for landfill
c	modified on 20 Aug 99 for land boundary only
c	modified on 21 Aug 99 for markings
c	modified on 29 Aug 99 for f90
c	modified on 10 Sep 99 for spl labels
c	modified on 5 May for no gaps 
c	colour codes used in gsplci
C	1-Black, 2-White, 3-Blue Dark, 4-Greenish Blue (Light), 5-Greenish Blue (Dark), 6-Light Green,
C	7-Dark Green, 8-light red, 9-red
C	12-Dark Redish Yellow, 13-Rose, 14-Red.

c define error file, fortran unit number, and workstation type,
c and workstation id.
      parameter (ierrf=6, lunit=2, iwkid=1,nran=1000)
c      parameter (mreg=160,nreg=120)
c	mreg: for x-axis;  nreg: for y-axis
      real zreg(mreg,nreg),xreg(mreg),yreg(nreg),gx(mreg),gy(nreg)
      real sumw(mreg,nreg),sumz(mreg,nreg)
      integer rlvl(25),bw
      character *1 splct,shdng,onlbl
      character (20) off
c	character (:) off !, allocatable :: off(14)
      character *25 gm*1

      external color
	common lndfl


	open (unit=1,file='ncrinp',status='old')
cc	constants defined
	onlbl='y'	!y to plot only one label
	read (1,'(a)')onlbl


c     
c
c get data array.
c
c
c open gks, and turn clipping off.
c
	print '(/a,$)',' Output type (1/8/11/20): '
c	 1 for ngcm, 8 for display, 11 for pdf, 20 for postscript
	read *,iwtype
	print '(/a,$)',' Output file: '
	read '(a)',  off
	noff=len_trim(off)
	open (unit=555,status='unknown',file=off(1:noff))
	print '(/a,$)',' Land overlay requied (y/n) ? '
	read '(a)',lndfl
        call gopks (ierrf, iszdm)
	call ngsetc('ME',off)
        call gopwk (iwkid, lunit, iwtype)
        call gacwk (iwkid)
        call gsclip(0)
	
	if (onlbl.eq.'y') then
		print '(/a,$)',' No. of contours tbp: '
		read *, nc 
	end if
	print '(/a,$)',' Value ntbc: '
	read *, thv 
	read (1,*)ipwr
	read (1,*)wf
	read (1,'(a)')gm
c
c call conpack color fill routine.

	print '(/a,$)',' Colour, black & white or line drawing (1/2/3) ?
     1   '
	read *,bw 
	if (bw.eq.3) then
		print '(/a,$)',' Shading required (y/n) ? '
		read '(a)',shdng
		if (shdng.eq.'y') then
			print '(/a,$)',' min/max contour valuess for 
     1                  shading: '
			read *,smn,smx
		end if
		print '(/a,$)',' Contour label size (in fractions of 1)
     1          : '
		read *,rls
	end if
	  print '(/a,$)',' Special contours requied (y/n) ? '
	  read '(a)',splct
	  if (splct.eq.'y'.and.onlbl.eq.'y') then
		print '(/a,$)',' Feed values:  '
		do kk=1,nc
			read*,rlvl(kk)
		end do
	  else if (onlbl.eq.'y') then
		print '(/a,$)',' Contour min/max (0,0 for dflt): '
		read *,cmnn,cmxx
	  end if
	print '(/a,$)',' No. of sets tbp: '
	read *, ns
	npts=ns
	do kkk = 1,ns
	print'(/a,$)','xmin,xmax,ymin,ymax:'
	read*,xmn,xmx,ymn,ymx
	ymd=(ymn+ymx)*.5
	xmd=(xmn+xmx)*.5
     	call ngrd(xreg,yreg,zreg,mreg,nreg,xmn,xmx,ymn,ymx,izmn,izmx,
     1	sumw,gy,gx,sumz,thv,ipwr,wf,gm)
        if (bw.eq.1.or.bw.eq.2) call colcon(onlbl,bw,zreg,mreg,nreg,-15,
     1	color,'ce',ymn,ymx,xmn,xmx,ymd,xmd,iwkid,kkk,npts
     1	,izmn,izmx,nc,cmnn,cmxx,rlvl,splct)
cc     1	,izmn,izmx,nc,cmnn,cmxx,rlvl,splct,nc,rls)
        if (bw.eq.3) call ccpscm(zreg,mreg,nreg,-15,color,
     1       'ce',ymn,ymx,xmn,xmx,ymd,xmd,iwkid,kkk,npts
     1	,izmn,izmx,nc,cmnn,cmxx,smn,smx,shdng,rlvl,splct,rls)
	end do
c
c close frame and close gks.
c
        call frame
        call gdawk (iwkid)
        call gclwk (iwkid)
        call gclks
	return
c
	end

      	subroutine colcon(onlbl,bw,zreg,mreg,nreg,ncl,color,proj,rlatmn,
     1	rlatmx,rlonmn,rlonmx,plat,plon,iwkid,kkk,npts,izmn,izmx,
     1	nc,cmnn,cmxx,rlvl,splct)
cc   1	nc,cmnn,cmxx,rlvl,splct,splvl)

c	xmn,xmx,ymn,ymx are changed while calling this routine

	parameter (lrwk=5000,liwk=5000,lmap=2000000,nwrk=150000,
     1  nogrps=5)
      	integer mreg,nreg,iwrk(liwk), lfin(75)
      	real zreg(mreg,nreg),rwrk(lrwk), xwrk(nwrk), ywrk(nwrk)
      	integer map(lmap),iarea(nogrps),igrp(nogrps)
      	character*4 proj
      	character*10 lbls(40)
	character *40 mrkfl
	character *1 lndfl,splct,xaxt,yaxt,splay,splax,mrk,onlbl,pl,
     1  cntrs,
     1	bndry,clrt,side,lblbr,mrkonplt,lndstd,clrl
	dimension ifil(1),px(2000),py(2000),px1(2000),py1(2000),px2
     1  (2000)
     1	,py2(2000),px3(2000),py3(2000),px4(2000),py4(2000)
	character*6 laby(50),labx(50)
	common lndfl
	integer rlvl(25),bw
	character (90) hd,xax,yax,hdm,hdm1
cc	character *60 xax,yax,hdm,hdm1
	character *1 hdng(90),xaxs(60),yaxs(60),hdngm(90),splv
        character *1 hdngm1(90)
cc	equivalence (hd,hdng),(xaxs,xax),(yaxs,yax),(hdm,hdngm),(hdm1,
cc     1  hdngm1)
	character *1 mhtb,rlin,grdonmp,stmp,rlaxs,lrg
      	external fill
      	external mask
      	external color
	data ifil /1/
c	conditions:
c	bndry='y'	!plts only land boundary if bndry is y and landfill through term is n
	if (kkk.eq.1) read (1,'(a)')bndry
c	lndstd ='y'	!lnd bndry w/o fill - for land studies
	if (kkk.eq.1) read (1,'(a)')lndstd
	if (lndstd.eq.'y') then
	  if (lndfl.eq.'n') bndry='y' !cmmnt this for land based studies
	end if
c	cntrs='n'	!y to have cntrs on clrs. for cntr lbls/size look for llp/lls
	if (kkk.eq.1) read (1,'(a)')cntrs
	splay='n'    		!y-axis lbl through input 
	splax='n'    		!x-axis lbl through input 
c	mnhdng=2	!for two main headings
	if (kkk.eq.1) read (1,*)mnhdng
c	mrk='n'
	if (kkk.eq.1) read (1,'(a)')mrk
	if (kkk.eq.1) read (1,*)nutr
c	mrkonplt='n'	!marks on the existing plot
	if (kkk.eq.1) read (1,'(a)')mrkonplt
c	clrt='y'	!y for colour titles
c	clrl='y'	!y for color trakcs
	if (kkk.eq.1) read (1,'(a)')clrl
	if (kkk.eq.1) read (1,'(a)')clrt
c	lblbr='y'
	if (kkk.eq.1) read (1,'(a)')lblbr
	if (kkk.eq.1.and.mrk.eq.'y') lblbr='n' !doest pritn label bar
	if (bw.eq.2) clrt='n'
	if (mrk.eq.'y') then
		print '(a,$)',' File for marking: '  !after search rad
		read '(a)',mrkfl
		open (unit=73,status='old',file=mrkfl)
	print *,'entd'
	end if
c	test for x & y axes titles
	if (npts.ne.1) then
		xaxt='n'
		yaxt='n'
	else 
		xaxt='y'
		yaxt='y'
	end if
	pl='l'		!p/l for portrait/landscape
	if (pl.eq.'p'.and.npts.eq.6) then
		if (kkk.eq.1.or.kkk.eq.3.or.kkk.eq.5) yaxt='y'
		if (kkk.eq.5.or.kkk.eq.6) xaxt='y'
	else if (pl.eq.'l'.and.npts.eq.6) then
		if (kkk.eq.4.or.kkk.eq.5.or.kkk.eq.6) xaxt='y'
		if (kkk.eq.1.or.kkk.eq.4) yaxt='y'
	end if
	if (npts.eq.5.and.(kkk.eq.1.or.kkk.eq.2.or.kkk.eq.4)) yaxt='y'	
	if (npts.eq.5.and.(kkk.eq.1.or.kkk.eq.4.or.kkk.eq.5)) xaxt='y'	
	if (npts.eq.3) yaxt='y'	! plts yaxix ttl
	if (npts.eq.3.and.kkk.eq.3) xaxt='y'	! plts xaxix ttl
	if (npts.eq.4.and.(kkk.eq.1.or.kkk.eq.3)) yaxt='y'	! plts yaxix ttl
	if (npts.eq.4.and.(kkk.eq.3.or.kkk.eq.4)) xaxt='y'	! plts xaxix ttl
c	side='y'
	if (kkk.eq.1) read (1,'(a)')side
c	following two statements for printing two figs side by side
	if (side.eq.'y') then
		if (npts.eq.2.and.(kkk.eq.1.or.kkk.eq.2)) xaxt='y'	! plts xaxix ttl
		if (npts.eq.2.and.kkk.eq.1) 		  yaxt='y'	! plts yaxix ttl
c	following two statements for writing one below the other 
	else if (side.eq.'n') then
		if (npts.eq.2.and.(kkk.eq.1.or.kkk.eq.2)) yaxt='y'	! plts xaxix ttl
		if (npts.eq.2.and.kkk.eq.2) 		  xaxt='y'	! plts yaxix ttl
	end if

c	write x and y axes if lable bar is to be put after all plots
	if (onlbl.eq.'n') then
		xaxt='y'
		yaxt='y'
	end if
c
c set color fill to solid.
c
      call gsfais (1)
c
c initialize areas.
	if (kkk.eq.1) then
		tlxmn=999
		tlxmx=-999
		tlymn=999
		tlymx=-999
	end if
c
c
c initialize ezmap and add to area map.
c
      	call mapstr ('gr',0.)
      	call mapstc ('ou','co')

	print'(/a,$)',' No. of divisions on x & y axes: '
	read*,nx,ny
	print '(/a,$)', ' Special x-axis required (y/n)? '
	read '(a)',splax
	if (splax.eq.'y') then
		print '(/a,$)', ' Special lable for x-axis (xdiv+1): '
		do kk = 1,nx+1
			read '(a)',labx(kk)
		end do
	end if
	print '(/a,$)', ' Special y-axis required (y/n)? '
	read '(a)',splay
	if (splay.eq.'y') then
		print '(/a,$)', ' Special lable for y-axis (ydiv+1): '
		do k = 1,ny+1
		read '(a)',laby(k)
		end do
	end if
	print '(/a,$)',' Data gaps to be skipped (y/n): '
	read '(a)', splv
	print '(/a,$)',' 4 coordinate positions for figure (in fractions o
     cf 1) '
	read *,vpxmn,vpxmx,vpymn,vpymx
	call mappos(vpxmn,vpxmx,vpymn,vpymx)
      	call maproj(proj,0.,0.,0.0)
        mpMinLonF = 0. 
        mpMaxLonF = 200. 
        mpCenterLonF = (mpMinLonF+mpMaxLonF)/2. 
c      	if (proj.eq.'sv') call mapstr ('sa',10.)
      	call mapset('co',rlatmn,rlonmn,rlatmx,rlonmx)
      	call mapint
      	call arinam(map, lmap)
      	call mapbla(map)
c
c initialize conpack and add to area map.
c
	if (onlbl.eq.'n') then
		print '(/a,$)',' No. of contours tbp: '
		read *, nc 
	end if 
	  if (splct.eq.'y'.and.onlbl.eq.'n') then
		print '(/a,$)',' Feed values:  '
		do kk=1,nc
			read*,rlvl(kk)
		end do
	  else if (onlbl.eq.'n') then
		print '(/a,$)',' Contour min/max (0,0 for dflt): '
		read *,cmnn,cmxx
	  end if
	if (cmnn.ne.0.or.cmxx.ne.0) then
		rmx=cmxx
		rmn=cmnn
	else
		rmx=izmx
		rmn=izmn
	end if
c	nc=(rmx-rmn)/cis+1
	cis=(rmx-rmn)/float(nc)
	ncll=nc+1
	if (splct.eq.'y') ncll=nc	!splvl
     	call cpseti('set - do-set-call flag',0)
      	call cpseti('map - mapping flag',1)
      	call cpsetr('orv - out of range flag',1.e12)
      	call cpsetr('xc1 - x coordinate at index 1',rlonmn)
    	call cpsetr('xcm - x coordinate at index m',rlonmx)
      	call cpsetr('yc1 - y coordinate at index 1',rlatmn)
      	call cpsetr('ycn - y coordinate at index n',rlatmx)
      	if (splct.eq.'y') then
		call cpseti('cls - contour level slctn flg',0) !to use cntrs frm rlvl
      	else
        	call cpseti('cls - contour level slctn flg',+nc) ! 
      	end if
	call cpseti('ncl - no of cntr lvls', ncll)

	if (splct.eq.'y') then
      	   do 10, i=1,ncll
         	call cpseti('pai - parameter array index',i)
         	call cpsetr('clv - contour level values',float(rlvl(i)))
         	call cpseti('clu - contour level use flag',3) 	!0/1/2/3 no cntr/only cntr/only lbl/lbl&cntr
		mxcntr=clv		!mx cntr value
 10   	   continue
        end if
	if (kkk.eq.1) read (1,*)smthfc
      	call cpsetr('t2d - smthng fctr', smthfc)
      	if (splct.eq.'n') then
      		call cpsetr('cis - contour interval specifier',cis)
      		call cpsetr('lis - label interval specifier',2.)
      		call cpsetr('cmx - contour max  specifier',rmx)	
      		call cpsetr('cmn - contour min  specifier',rmn)
      	end if
      	if(splv.eq.'y')call cpsetr('spv - spl vl for n/a data',-888.8)
	if (kkk.eq.1) read (1,*)lblpst
	if (kkk.eq.1) read (1,*)lblorn
	if (kkk.eq.1) read (1,*)szlbl
      	if (cntrs.eq.'y')
     1	call cpseti('llp - line label positioning flag',lblpst)	!2/3 for cntr lbls; 1 for no lbls. also clu
      	call cpseti('llo - line label orientation flag',lblorn)	
      	call cpsetr('lls - line label size ',szlbl)
      	call cpseti('llb - line label box flag',0)
      	call cpseti('hlb - high/low label box flag',0)
      	call cpseti('hlc - high/low label color index',2)
	call cpseti('ilp bottom info',0)		!c to have btm info
	call cpsetc('ilt bottom info'," ")		!  -- do --
	call cpsetc('lot high/low'," ")			!to sprs h/l info
	call cpsetc('hit  high/low'," ")			!  -- do --
      	call cprect(zreg, mreg, mreg, nreg, rwrk, lrwk, iwrk, liwk)
      	call cpclam(zreg, rwrk, iwrk, map)
c      	call cplbam(zreg, rwrk, iwrk, map)  !for high/low label

c
c choose a color for every contour level.
c
      	if (onlbl.eq.'n'.and.kkk.ge.1) call color (0,iwkid,kkk,bw)
      	call color (ncll+3,iwkid,kkk,bw)
      	if (lndfl.eq.'y') call landfill(rlatmn,rlonmn,rlatmx,rlonmx)
c
c fill contours and areas over land.
c
c	supress this statement if only dot plots are required
      	if (mrkonplt.eq.'y') ! if 'n' no contours come
     1	call arscam(map, xwrk, ywrk, nwrk, iarea, igrp, nogrps, fill)
c
c draw continental outlines, labels, and masked contours.
c
	call gsplci(1)   !this sets the black colour to the  border
      	if (lndfl.eq.'y'.or.bndry.eq.'y') call maplot	!draw continental lines
      	call cplbdr(zreg,rwrk,iwrk)
      	if (cntrs.eq. 'y')
     1	call cpcldm(zreg,rwrk,iwrk,map,mask)	! to have cntrs ovr clrs
c
c draw and fill a label bar.
c
c	call gsplci(1)   !this sets the black colour to the  border
      	call getset(xminn,xmaxx,yminn,ymaxx,dum1,dum2,dum3,dum4,idum)
	tlxmn=min(tlxmn,xminn)
	tlxmx=max(tlxmx,xmaxx)
	tlymn=min(tlymn,yminn)
	tlymx=max(tlymx,ymaxx)
	ybot=yminn-.06		
c	read *,yymm
	if (kkk.eq.1) read (1,*) ybtad
	if (kkk.eq.1) read (1,*) ybtad1
	if (kkk.eq.1) read (1,*) ytpad
	if (kkk.eq.1) read (1,*) ytpad1
	ybot=tlymn -ybtad
	ybot1=yminn-ybtad1
	ytop=ybot+ytpad
	ytop1=ybot1+ytpad1
	if (kkk.eq.1) read (1,'(a)')rlin
c	if (onlbl.eq.'n'.or.(kkk.eq.npts.and.lblbr.eq.'y')) then
	if (onlbl.eq.'n'.or.(kkk.eq.npts.and.lblbr.eq.'y')) then !chang
      	  do i=1,ncll
		lfin(i)=i+2
         	call cpseti('pai - parameter array index',i)
         	call cpgetr('clv - contour level values',clv)
		iclv=clv
		if (rlin.eq.'i') then
         		write (lbls(i), '(i4)') iclv
		else if (rlin.eq.'r') then
         		write (lbls(i), '(f5.1)') clv
		end if
		mxcntr=clv		!mx cntr value
		rxcntr=clv		!mx cntr value
          end do
      	  lfin(ncll+1)=ncll+3
c      	  write (lbls(ncll+2),'(a4)') 'land'
      	  if (splv.eq.'y') then
		write (lbls(ncll+2),'(a3)') 'gap'
		if (rlin.eq.'i') then
      			write (lbls(ncll+1),'(a1,i4)') '>',mxcntr
		else if (rlin.eq.'r') then
      			write (lbls(ncll+1),'(a1,f4.1)') '>',rxcntr
		end if
  	  else if (splv.eq.'n'.and.onlbl.eq.'y') then
		if (rlin.eq.'i') then
      			write (lbls(ncll+1),'(a1,i4)') '>',mxcntr
		else if (rlin.eq.'r') then
      			write (lbls(ncll+1),'(a1,f4.1)') '>',rxcntr
		end if
	  else if (splv.eq.'n'.and.onlbl.eq.'n') then
		if (rlin.eq.'i') then
      			write (lbls(ncll+0),'(a1,i4)') '>',mxcntr
		else if (rlin.eq.'r') then
      			write (lbls(ncll+0),'(a1,f4.1)') '>',rxcntr
		end if
	  end if 
      	  lfin(ncll+2)=2
	  call gsplci(1)   !this sets the black colour to the label border
      	  if (splv.eq.'y'.and.onlbl.eq.'y') call lblbar(0,tlxmn,tlxmx,
     1	  ybot,ytop,ncll+2,1.,.5,lfin,1,lbls,ncll+2,1)	!for gap
      	  if (splv.eq.'y'.and.onlbl.eq.'n') call lblbar(0,xminn,xmaxx,
     1	  ybot1,ytop1,ncll+2,1.,.5,lfin,1,lbls,ncll+2,1)	!for gap
      	  if (splv.eq.'n'.and.onlbl.eq.'y') call lblbar(0,tlxmn,tlxmx,
     1	  ybot,ytop,ncll+1,1.,.5,lfin,1,lbls,ncll+1,1)	!for no gap
      	  if (splv.eq.'n'.and.onlbl.eq.'n') call lblbar(0,xminn,xmaxx,
     1	  ybot1,ytop1,ncll+1,1.,.5,lfin,1,lbls,ncll+2,1)	!for no gap
	end if
	
	if (kkk.eq.1) read (1,*)yxvl
	if (kkk.eq.1) read (1,*)xxvl
	if (kkk.eq.1) read (1,'(a)')rlaxs

	if (splay.eq.'y') then
		ixdc=9999
	else
		ixdc=yxvl		!for y axis values pstn like lats
	end if
	if (splax.eq.'y') then
		iydc=9999
	else
		iydc=xxvl		!for x axis
	end if
	if (xaxt.eq.'n') iydc=999 
	if (yaxt.eq.'n') ixdc=999 
	if (rlaxs.eq.'i') then
		call labmod('(i4)','(i4)',5,5,14,14,ixdc,iydc,0)
	else if (rlaxs.eq.'r') then
		call labmod('(f5.1)','(f4.1)',5,5,14,14,ixdc,iydc,0)
	end if
	call gridal(nx,0,ny,0,1,1,5,0.,0.)
	if (kkk.eq.1) read (1,'(a)') grdonmp
	if (kkk.eq.1) read (1,*)i1,i2,i3,i4
	if (grdonmp.eq.'y') then
		call gridl(i1,i2,i3,i4)	!to put the grids on the map
	end if
	print '(/a,$)',' Sub-heading: '
cc	read '(q,"ncr",a)',ncr,(hdng(i),i=1,ncr)
c	read '(q,a)',ncr,hd
	read '(a)',hd
	ncr=len_trim(hd)
	print '(/a,$)',' Title for x-axis: '
c	read '(q,a)',ncrx,xax
	read '(a)',xax
	ncrx=len_trim(xax)
	print '(/a,$)',' Title for y-axis: '
	read '(a)',yax
	ncry=len_trim(yax)
c	read '(q,a)',ncry,yax
c
c	plots markers on the map
c	if (kkk.eq.1) read (1,*)mrktp
	if (kkk.eq.1) read (1,*)mrkcn
	if (kkk.eq.1) read (1,*)mrkgspl
	if (kkk.eq.1) read (1,*)iaxgspl
	if (kkk.eq.1) read (1,*)ihdgspl
	if (kkk.eq.1) read (1,*)ibtgspl
	if (kkk.eq.1) read (1,*)istpgspl
	if (mrk.eq.'y') then
		mi=0
c		call set(tlxmn,tlxmx,tlymn,tlymx,rlonmn,rlonmx,rlatmn,rlatmx,0)
		call set(xminn,xmaxx,yminn,ymaxx,rlonmn,rlonmx,rlatmn,r
     1          latmx,0)
555		mi=mi+1	
		if (nutr.eq.1) then  ! 2 for plotting 2 tracks cmk study
		    read (73,*,end=99)px(mi),py(mi)
		else if (nutr.eq.2) then
		    read (73,*,end=99)py(mi),px(mi),py1(mi),px1(mi)
		else if (nutr.eq.3) then
		    read (73,*,end=99)py(mi),px(mi),py1(mi),px1(mi),
     1	            py2(mi),px2(mi)
		else if (nutr.eq.4) then
		    read (73,*,end=99)py(mi),px(mi),py1(mi),px1(mi),
     1	            py2(mi),px2(mi),py3(mi),px3(mi)
		else if (nutr.eq.5) then
		    read (73,*,end=99)py(mi),px(mi),py1(mi),px1(mi),
     1	            py2(mi),px2(mi),py3(mi),px3(mi),py4(mi),px4(mi)
		end if
		go to 555
99	        continue
		close (73)
		mi=mi-1
		call gsplci(1)
     		call points(px,py,mi,-3,mrkcn)	!-1 for .; -2 +; -3 *; -4 0; -5 x
     		if (nutr.eq.2) then
			if (clrl.eq.'y') call gsplci(2)
			call points(px1,py1,mi,-4,mrkcn)	
		end if
     		if (nutr.eq.3) then 
			if (clrl.eq.'y') call gsplci(2)
			call points(px1,py1,mi,-4,mrkcn)	
			if (clrl.eq.'y') call gsplci(3)
			call points(px2,py2,mi,-2,mrkcn)
		end if
     		if (nutr.eq.4) then
			if (clrl.eq.'y') call gsplci(2)
			call points(px1,py1,mi,-4,mrkcn)	
			if (clrl.eq.'y') call gsplci(3)
			call points(px2,py2,mi,-2,mrkcn)
			if (clrl.eq.'y') call gsplci(14)
			call points(px3,py3,mi,-1,mrkcn)
		end if
     		if (nutr.eq.5) then
			if (clrl.eq.'y') call gsplci(2)
			call points(px1,py1,mi,-4,mrkcn)	
			if (clrl.eq.'y') call gsplci(3)
			call points(px2,py2,mi,-2,mrkcn)
			if (clrl.eq.'y') call gsplci(14)
			call points(px3,py3,mi,-1,mrkcn)
			if (clrl.eq.'y') call gsplci(11)
			call points(px4,py4,mi,-5,mrkcn)
		end if
	end if
	call gsplci(mrkgspl)	!sets the colour of the mark

	if (npts.eq.kkk) then
		print '(/a,$)',' Main heading: '
c		read '(q,a)',ncrh,hdm
		read '(a)',hdm
		ncrh=len_trim(hdm)
		if (mnhdng.eq.2) then
		print '(/a,$)',' Main heading2: '
c		read '(q,a)',ncrh1,hdm1
		read '(a)',hdm1
		ncrh1=len_trim(hdm1)
		end if
		print '(/a,$)',' Main heading at top/bottom? (t/b)'
		read '(a)',mhtb
	end if
	xstp=(xminn+xmaxx)/2.	!stng pstn for xax lbl
	ystp=(yminn+ymaxx)/2.	!stng pstn for yax lbl
        call gselnt(0)
c	wdy=0.025			!to takecare of overflow of labels
c	wds=0.05			!width of letters for lables
	if (kkk.eq.1) read (1,*) wdy
	if (kkk.eq.1) read (1,*) wds
	if (kkk.eq.1) read (1,*)sylblpst
	if (kkk.eq.1) read (1,*)sxlblpst
	call gsplci(1)   !sets black colour
	if (splay.eq.'y') then
		rnyd=(ymaxx-yminn)/(ny)
		rnydd=yminn-rnyd-.002
		do ik=1,ny+1
			rnydd=rnydd+rnyd
 			call pchiqu(xminn-sylblpst,rnydd,laby(ik),0.009
     1                  ,0.,1.)
		end do
	end if

	if (splax.eq.'y') then
		rnxd=(xmaxx-xminn)/(nx)
		rnxdd=xminn-rnxd-.002
		do ik=1,nx+1
			rnxdd=rnxdd+rnxd
 			call pchiqu(rnxdd,yminn-sxlblpst,labx(ik),0.009
     1                  ,0.,0.)
		end do
	end if
c11111111111111111111111111111111111111111111111111111111111111111111111
 	call pchiqu((xminn+xmaxx)/2.,ymaxx+.02,hd(1:ncr),0.012,0.,0.)	
	if (xaxt.eq.'n') xstp=999 
	if (yaxt.eq.'n') ystp=999 
	if (clrt.eq.'y') call gsplci(iaxgspl)
	if (kkk.eq.1) read (1,*) hdcr
	if (kkk.eq.1) read (1,*) xaxcr
	if (kkk.eq.1) read (1,*) sxaxcr
	if (kkk.eq.1) read (1,*) yaxcr

c 	call pchiqu((xminn+xmaxx)/2.,ymaxx+hdcr,hd(1:ncr),0.012,0.,0.)	
	if (clrt.eq.'y') call gsplci(iaxgspl)
 	if (splax.eq.'n')call pchiqu(xstp,yminn-xaxcr,xax(1:ncrx),
     1  0.012,0.,0.) 
 	if (splax.eq.'y')call pchiqu(xstp,yminn-sxaxcr,xax(1:ncrx),0.012
     1  ,0.,0.) 
 	call pchiqu(xminn-yaxcr,ystp,yax(1:ncry),0.012,90.,0.)

	if (kkk.eq.1) read (1,*)thdcr1
	if (kkk.eq.1) read (1,*)thdcr2
	if (kkk.eq.1) read (1,*)bhdcr1
	if (kkk.eq.1) read (1,*)bhdcr2
	if (kkk.eq.1) read (1,'(a)')lrg
	if (kkk.eq.1) read (1,*)hdmps
	if (kkk.eq.1) read (1,*)hdmpsl
     	if (npts.eq.kkk) then
 		if (mhtb.eq.'t') then
		if (clrt.eq.'y') call gsplci(ihdgspl)
		if (mnhdng.eq.2) then
			rmmhh1=thdcr1	!.12
			rmmhh2=thdcr2	!.08
		else
			rmmhh1=thdcr2	!.08
		end if
		call pchiqu((tlxmn+tlxmx)/2.,tlymx+rmmhh1,hdm(1:ncrh),
     1          0.02,0.,0.)
		call pchiqu((tlxmn+tlxmx)/2.,tlymx+rmmhh2,hdm1(1:ncrh1)
     1          ,0.02,0.,0.)
		end if
		call gsplci(1)   !sets black colour
 		if (mhtb.eq.'b') then
			if (clrt.eq.'y') call gsplci(ibtgspl)
			if (mnhdng.eq.2) then
				rmmhh1=bhdcr1	!-0.16
				rmmhh2=bhdcr2	!-0.19
			else
				rmmhh1=bhdcr2	!-0.18
			end if

			if (lrg.eq.'n') then
c              		    call pchiqu(hdmps,tlymn+rmmhh1,hdm(1:ncrh)
c     1			    ,0.014,0.,-1.)              		    
c			    call pchiqu(hdmps,tlymn+rmmhh2,
c     1			    hdm1(1:ncrh1),0.014,0.,-1.)
               		    call pchiqu(tlxmn-0.054,tlymn+rmmhh1,
     1			    hdm(1:ncrh),0.014,0.,-1.)
               		    call pchiqu(tlxmn-0.054,tlymn+rmmhh2,
     1			    hdm1(1:ncrh1),0.014,0.,-1.)
c               use the following for one fig and large text
               		 else if (lrg.eq.'y') then
               		    call pchiqu(tlxmn-hdmpsl,tlymn+rmmhh1,
     1			    hdm(1:ncrh),0.014,0.,-1.)
               		   call pchiqu(tlxmn-hdmpsl,tlymn+rmmhh2,
     1			   hdm1(1:ncrh1),0.014,0.,-1.)
               		 end if
		end if
	end if
	call gsplci(istpgspl)
	if (kkk.eq.1) read (1,'(a)')stmp
	print *,'stmp', stmp
	if (kkk.eq.npts.and.stmp.eq.'y') then
	 	call pchiqu(tlxmx-.2,tlymn+.02,"NRSA/DOS",0.014,0.,-1.)
	end if
      	return
c
      end



      	subroutine ccpscm(zreg,mreg,nreg,ncl,color,
     1  proj,rlatmn,rlatmx,rlonmn,rlonmx,plat,plon,iwkid,kkk,npts
     1	,izmn,izmx,nc,cmnn,cmxx,smn,smx,shdng,rlvl,splct,rls)

      	parameter (lrwk=5000,liwk=10000,lmap=50000,nwrk=1500,nogrps=5)
      	integer mreg,nreg,iwrk(liwk)
      	real zreg(mreg,nreg),rwrk(lrwk), xwrk(nwrk), ywrk(nwrk)
      	integer map(lmap),iarea(nogrps),igrp(nogrps),rlvl(25)
      	character*4 proj

	dimension ifil(1)
	character*14 laby(50),labx(50)
	data ifil /1/
	character (60) hd,xax,yax,hdm,hdm1
	character *1 hdng(60),xaxs(60),yaxs(60),hdngm(60),splv,hdngm1
     1  (60),
     1  xaxt,yaxt,pl,splay,splax,splct,bndry,onlbl,shdng,lndfl,clrt,mhtb
	equivalence (hd,hdng),(xaxs,xax),(yaxs,yax),(hdm,hdngm),(hdm1,
     1  hdngm1)
	common lndfl
      	external sfill
      	external cpdrpl
      	external mask
      	external clr
	mnhdng=2
c	test for x & y axes titles
	if (npts.ne.1) xaxt='n'
	if (npts.ne.1) yaxt='n'
	pl='p'		!p/l for portrait/landscape
	if (pl.eq.'p'.and.npts.eq.6) then
		if (kkk.eq.1.or.kkk.eq.3.or.kkk.eq.5) yaxt='y'
		if (kkk.eq.5.or.kkk.eq.6) xaxt='y'
	else if (pl.eq.'l'.and.npts.eq.6) then
		if (kkk.eq.4.or.kkk.eq.5.or.kkk.eq.6) xaxt='y'
		if (kkk.eq.1.or.kkk.eq.4) yaxt='y'
	end if
	if (npts.eq.4.and.(kkk.eq.1.or.kkk.eq.3)) yaxt='y'	! plts yaxix ttl
	if (npts.eq.4.and.(kkk.eq.3.or.kkk.eq.4)) xaxt='y'	! plts xaxix ttl
	if (npts.eq.2.and.(kkk.eq.1.or.kkk.eq.2)) yaxt='y'	! plts yaxix ttl
	if (npts.eq.2.and.kkk.eq.2) 		  xaxt='y'	! plts xaxix ttl
c
c set color fill to solid.
c
      	call clr (2,iwkid,kkk)
	call gsfais(1)
c
c initialize areas.
	if (kkk.eq.1) then
		tlxmn=999
		tlxmx=-999
		tlymn=999
		tlymx=-999
	end if
c
c
c initialize ezmap and add to area map.
c
      call mapstr ('gr',0.)
      call mapstc ('ou','co')

	splay='n'    		!y-axis lbl through input 
	splax='n'    		!x-axis lbl through input 
	print'(/a,$)','No. of divisions on x & y axes: '
	read*,nx,ny
	print '(/a,$)', ' Special x-axis label required (y/n)? '
	read '(a)',splax
	if (splax.eq.'y') then
		print '(/a,$)', ' Special label for x-axis (xdiv+1): '
		do kk = 1,nx+1
		read '(a)',labx(kk)
		end do
	end if
	print '(/a,$)', ' Special y-axis label required (y/n)? '
	read '(a)',splay
	if (splay.eq.'y') then
		print '(/a,$)', ' Special label for y-axis (ydiv+1): '
		do k = 1,ny+1
		read '(a)',laby(k)
		end do
	end if
	print '(/a,$)',' Data gaps to be skipped (y/n): '
	read '(a)', splv
	print '(/a,$)',' Four coordinate positions for figure 
     1	(in fractions of 1): '
	read *,vpxmn,vpxmx,vpymn,vpymx
c	compues cntr lbl sz
	vpxd=vpxmx-vpxmn
	vpyd=vpymx-vpymn
	vpmx=max(vpxd,vpyd)
	rls=0.012/vpmx

c
	call mappos(vpxmn,vpxmx,vpymn,vpymx)
      call maproj(proj,0.,0.,0.0)
c      if (proj.eq.'sv') call mapstr ('sa',10.)
      call mapset('co',rlatmn,rlonmn,rlatmx,rlonmx)
      call mapint
      call arinam(map, lmap)
      call mapbla(map)
c
c initialize conpack and add to area map.
c
	if (cmnn.ne.0.and.cmxx.ne.0) then
		rmx=cmxx
		rmn=cmnn
	else
		rmx=izmx
		rmn=izmn
	end if
	cis=(rmx-rmn)/float(nc)
	ncll=nc+1
	if (splct.eq.'y') ncll=nc		!splvl
      	call cpseti('nof',3)  !to omit decimal vals
      	call cpseti('set - do-set-call flag',0)
      	call cpseti('map - mapping flag',1)
c      	call cpsetr('orv - out of range flag',1.e12)
      	call cpsetr('xc1 - x coordinate at index 1',rlonmn)
      	call cpsetr('xcm - x coordinate at index m',rlonmx)
      	call cpsetr('yc1 - y coordinate at index 1',rlatmn)
      	call cpsetr('ycn - y coordinate at index n',rlatmx)
      	call cpsetr('t2d - smthng fctr', 25.)
	if (splct.eq.'y') then
		call cpseti('cls - cntr lvl slctn flg',0)
      		call cpsetr('lis - label interval specifier',0.)
	else
		call cpseti('cls - cntr lvl slctn flg',+nc)
      		call cpsetr('cis - contour interval specifier',cis)
      		call cpsetr('cmx - contour max  specifier',rmx)	
      		call cpsetr('cmn - contour min  specifier',rmn)
	end if
      	call cpseti('ncl - number of contour levels',ncll)
      	do i=1,ncll
        	call cpseti('pai - parameter array index',i)
		if (splct.eq.'y') then
		call cpsetr('clv - cntr lvl vl',float(rlvl(i)))
		else
		call cpsetr('clv - cntr lvl vl',clv)
		end if
       		call cpseti('clu - contour level use flag',3) 	!0/1/2/3 no cntr/only cntr/only lbl/lbl&cntr
         	call cpseti('aia - area identifier above',0)
         	call cpseti('aib - area identifier below',0)
	end do 
c	 call cpsetr('cll - cntr ln thcknss', 2)
      	if(splv.eq.'y')call cpsetr('spv - spl value for n/a data',
     1  -888.8)
      	call cpseti('llp - line label positioning flag',2)	!to write cntr vls
	rc1=0.01
	call cpsetr('rc1 - regular sch const 1',rc1)
c	call cpsetr('rc2 - regular sch const 2',rc2)
c	call cpsetr('pc4- penalty const sch 4',.1)
c	call cpsetr('pc5- penalty const sch 5',.2)
c	call cpsetr('pc6- penalty const sch 6',.05)
c	call cpsetr('pw4- penalty weight sch 6',5.)
      	call cpsetr('lls - line lebel size', rls)
      	call cpsetc('hlt - hi/lo lbl text',' '' ')
      	call cpseti('llo - line lbl pstn',1)
c      	call cpseti('llb - line label box flag',0)	!for hi/lo vls
      	call cpseti('hlb - high/low label box flag',0)
      	call cpseti('hlc - high/low label color index',2)
	call cpseti('ilp bottom info',0)		!c to have btm info
	call cpsetc('ilt bottom info'," ")		!  -- do --
	call cpsetc('lot high/low'," ")			!to sprs h/l info
	call cpsetc('hit  high/low'," ")			!  -- do --
      	call cprect(zreg, mreg, mreg, nreg, rwrk, lrwk, iwrk, liwk)
      	call cpclam(zreg, rwrk, iwrk, map)
      	call cplbam(zreg, rwrk, iwrk, map)  !for high/low label

c set number of contour levels and initialize conpack
c
      	call cprect (zreg, mreg, mreg, nreg, rwrk, lrwk, iwrk, liwk)	 !repeated ?
      	call cppkcl (zreg, rwrk, iwrk)
c
c turn on line labeling and turn off area identifiers for all levels
c
c?      	call cpgeti('ncl - number of contour levels',ncll)
      	do i=1,ncll
        	call cpseti('pai - parameter array index',i)
		call cpgetr('clv - cntr lvl vl',clv)
c        	call cpseti('clu - contour level use flag',3) 	!0/1/2/3 no cntr/only cntr/only lbl/lbl&cntr
         	call cpseti('aia - area identifier above',0)
         	call cpseti('aib - area identifier below',0)
	end do 
c
c add contour levels at min and max, and set area ids so that 
c you can fill between them
c
	if (shdng.eq.'y') then	
      	call cpseti('ncll - number of contour levels',ncll+2)
      	call cpseti('pai - parameter array index',ncll+1)
      	call cpsetr('clv - contour level value',smn)  !shdng mn
      	call cpseti('clu - contour level use flag',3)
      	call cpseti('aia - area identifier above',1)
      	call cpseti('aib - area identifier below',2)
      	call cpseti('pai - parameter array index',ncll+2)
	call sfseti('do - 0 fr dsh 1 fr dts and chr',1)
c	call sfsetc('ch- chr input','*')	! fr any chrtr
      	call cpsetr('clv - contour level value',smx)   !shdng mx
      	call cpseti('clu - contour level use flag',3)
      	call cpseti('aia - area identifier above',3)
      	call cpseti('aib - area identifier below',1)
	end if
c
c draw perimeter
c
      	call cpback(zreg, rwrk, iwrk)
c
c initialize areas
c
      	call arinam(map, lmap)
c
c add contours to area map
c
      	call cpclam(zreg, rwrk, iwrk, map)
c
c add label boxes to area map
c
      	call cplbam(zreg, rwrk, iwrk, map)
c
c fill contours
c
      	call arscam(map, xwrk, ywrk, nwrk, iarea, igrp, nogrps, sfill)
c
c draw contours, masking label boxes
c
      	call cpcldm(zreg, rwrk, iwrk, map, cpdrpl)
c
c draw labels
c
	call cpsetr('ilx',.99)
	call cpsetr('ily',-.1)
      	call cplbdr(zreg, rwrk, iwrk)	!hi/lo bxs & cntr inf
c
c choose a color for every contour level.
c
c      	call clr (ncll+3,iwkid,kkk)
     	if (lndfl.eq.'y') call landfill(rlatmn,rlonmn,rlatmx,rlonmx)
c
c fill contours and areas over land.
c
c
c draw continental outlines, labels, and masked contours.
c
      	if (lndfl.eq.'y') call maplot	!draw continental lines
      	call cplbdr(zreg,rwrk,iwrk)
      	call cpcldm(zreg,rwrk,iwrk,map,mask)
c
c draw and fill a label bar.
c
      	call getset(xminn,xmaxx,yminn,ymaxx,dum1,dum2,dum3,dum4,idum)
	tlxmn=min(tlxmn,xminn)
	tlxmx=max(tlxmx,xmaxx)
	tlymn=min(tlymn,yminn)
	tlymx=max(tlymx,ymaxx)
	ybot=yminn-.16		
	ybot=tlymn-.14		
	ytop=ybot+.06
	if (kkk.eq.npts) then
	end if

	if (splay.eq.'y') then
		ixdc=9999
	else
		ixdc=2		!for y axis values
	end if
	if (splax.eq.'y') then
		iydc=9999
	else
		iydc=2		!for x axis
	end if
	if (xaxt.eq.'n') iydc=999 
	if (yaxt.eq.'n') ixdc=999 
	call labmod('(i3)','(i3)',5,5,14,14,ixdc,iydc,0)
	call gridal(nx,0,ny,0,1,1,5,0.,0.)
	print '(/a,$)',' Sub-heading: '
c	read '(q,a)',ncr,hd
	print '(/a,$)',' Title for x-axis: '
c	read '(q,a)',ncrx,xax
	print '(/a,$)',' Title for y-axis: '
c	read '(q,a)',ncry,yax

	if (npts.eq.kkk) then
		print '(/a,$)',' Main heading: '
c		read '(q,a)',ncrh,hdm
		read '(a)',hdm
		ncrh=len_trim(hdm)
		if (mnhdng.eq.2) then
		print '(/a,$)',' Main heading2: '
		read '(a)',hdm1
		ncrh1=len_trim(hdm1)
c		read '(q,a)',ncrh1,hdm1
		end if
		print '(/a,$)',' Main heading at top/bottom? (t/b)'
		read '(a)',mhtb
	end if
	xstp=(xminn+xmaxx)/2.	!stng pstn for xax lbl
	ystp=(yminn+ymaxx)/2.	!stng pstn for yax lbl
        call gselnt(0)
	wdy=0.025			!to takecare of overflow of labels
	wds=0.05			!width of letters for lables
	if (splay.eq.'y') then
		rnyd=(ymaxx-yminn)/(ny)
		rnydd=yminn-rnyd-.002
		do ik=1,ny+1
			rnydd=rnydd+rnyd
 			call pchiqu(xminn-0.04,rnydd,laby(ik),0.009,0.
     2                  ,1.)
		end do
	end if

	if (splax.eq.'y') then
		rnxd=(xmaxx-xminn)/(nx)
		rnxdd=xminn-rnxd-.002
		do ik=1,nx+1
			rnxdd=rnxdd+rnxd
 			call pchiqu(rnxdd,yminn-.015,labx(ik),0.009,0.
     1                  ,0.)
		end do
	end if

	if (xaxt.eq.'n') xstp=999 
	if (yaxt.eq.'n') ystp=999 
 	call pchiqu((xminn+xmaxx)/2.,ymaxx+0.02,hd(1:ncr),0.012,0.,0.)
 	if (splax.eq.'n')call pchiqu(xstp,yminn-0.04,xax(1:ncrx),0.012
     1  ,0.,0.) 
 	if (splax.eq.'y')call pchiqu(xstp,yminn-0.05,xax(1:ncrx),0.012
     1  ,0.,0.) 
 	call pchiqu(xminn-0.064,ystp,yax(1:ncry),0.012,90.,0.)

     	if (npts.eq.kkk) then
 		if (mhtb.eq.'t') then
		if (mnhdng.eq.2) then
			rmmhh1=.12
			rmmhh2=.08
		else
			rmmhh1=.08
		end if
		call pchiqu((tlxmn+tlxmx)/2.,tlymx+rmmhh1,hdm(1:ncrh),
     1          0.02,0.,0.)
		call pchiqu((tlxmn+tlxmx)/2.,tlymx+rmmhh2,hdm1(1:ncrh1)
     1          ,0.02,0.,0.)
		end if
 		if (mhtb.eq.'b') then
		if (mnhdng.eq.2) then
			rmmhh1=-0.16
			rmmhh2=-0.18
		else
			rmmhh1=-0.18
		end if
		call pchiqu(tlxmn-0.01,tlymn+rmmhh1,hdm(1:ncrh),0.014,
     1          0.,-1.)
		call pchiqu(tlxmn-0.01,tlymn+rmmhh2,hdm1(1:ncrh1),0.014
     1          ,0.,-1.)
		end if
	end if
c	call gsplci(1)
      	return
      	end

      	subroutine fill (xwrk,ywrk,nwrk,iarea,igrp,ngrps)
c
      	dimension xwrk(*),ywrk(*),iarea(*),igrp(*)
	character *1 lndfl,lbwf
	common lndfl
c	conditions:
	lbwf='n'	! y for land boundary w/o filling the land
      	idmap=-1
      	idcont=-1

      	do 10, i=1,ngrps
         if (igrp(i).eq.1) idmap=iarea(i)
         if (igrp(i).eq.3) idcont=iarea(i)
 10   	continue
c
c if the area is defined by 2 or fewer points, return to arscam.
c
      	if (nwrk .le. 3) return
c
c check if the area is over the map 
c
      	if ((idmap .gt. 0) .and. (idcont .gt. 0)) then
c
c if the area is over water, fill the contours with colors depending
c on their level.
c
        if (mapaci(idmap).eq.1.or.lndfl.eq.'n') then
            call gsfaci(idcont+2)
            call gfa(nwrk-1,xwrk,ywrk)
c
c if the area is over land, fill with gray.
c
        else if (lndfl.eq.'y') then 
            if (lbwf.eq.'n') call gsfaci(30)	!0/30 white/brown land  
            call gfa(nwrk-1,xwrk,ywrk)
         endif
      	endif        
c
      	return
c
      	end

      	subroutine sfill (xwrk,ywrk,nwrk,iarea,igrp,ngrps)
c
      	real xwrk(*),ywrk(*),iscr(50000)
      	integer iarea(*),igrp(*),rscr(50000)
      	do 10, i=1,ngrps
         if (igrp(i).eq.3) iarea3=iarea(i)
 10   	continue

      	if (iarea3 .eq. 1) then
c
c if the area is defined by 3 or more points, fill it
c
         call sfsetr('spacing',.009)
         call sfnorm(xwrk,ywrk,nwrk,rscr,50000,iscr,50000)
      	endif
c   
c otherwise, do nothing
c
      	return
      	end
      	subroutine mask(xwrk,ywrk,nwrk,iarea,igrp,ngrps)
c
      	integer iarea(ngrps),igrp(ngrps)
      	real xwrk(nwrk),ywrk(nwrk)
c
      	idmap=-1
      	idcont=-1
c
      	do 10, i=1,ngrps
         if (igrp(i).eq.1) idmap=iarea(i)
         if (igrp(i).eq.3) idcont=iarea(i)
 10   	continue
c
c if the line is defined by 1 or fewer points, return to cpcldm
c
      	if (nwrk .lt. 2) return
c
c draw the line if the area is over the map, and not over a label, or
c over land.
c
      	if ((idmap.gt.0).and.(idcont.gt.0).and.(mapaci(idmap).eq.1))
     1     then
         call curve(xwrk,ywrk,nwrk)
      	endif
c
      	return
c
      	end


c
c
      	subroutine color (n,iwkid,kkk,bw)

	integer bw
	dimension rgb25(3,27),mtrx(27,25)
	data mtrx/
     1  	05,09,25,00,00,00,00,00,00,00,00,00,00,00,
     1     	00,00,00,00,00,00,00,00,00,00,00,00,00,
     1	     	05,09,12,25,00,00,00,00,00,00,00,00,00,00,
     1		00,00,00,00,00,00,00,00,00,00,00,00,00,
     1		05,06,09,12,25,00,00,00,00,00,00,00,00,00,
     1		00,00,00,00,00,00,00,00,00,00,00,00,00,
     1		05,06,09,12,13,25,00,00,00,00,00,00,00,00,
     1		00,00,00,00,00,00,00,00,00,00,00,00,00,
     1		05,06,09,12,13,16,25,00,00,00,00,00,00,00,
     1		00,00,00,00,00,00,00,00,00,00,00,00,00,
     1		05,06,08,09,12,13,16,25,00,00,00,00,00,00,
     1		00,00,00,00,00,00,00,00,00,00,00,00,00,
     1		05,06,08,09,12,13,16,25,26,00,00,00,00,00,
     1		00,00,00,00,00,00,00,00,00,00,00,00,00,
     1		05,06,08,09,11,12,13,16,25,26,00,00,00,00,
     1		00,00,00,00,00,00,00,00,00,00,00,00,00,
     1		03,05,06,08,09,11,12,13,16,25,26,00,00,00,
     1		00,00,00,00,00,00,00,00,00,00,00,00,00,
     1		03,05,06,08,09,11,12,13,14,16,25,26,00,00,
     1		00,00,00,00,00,00,00,00,00,00,00,00,00,
     1		03,04,05,06,08,09,11,12,13,14,16,25,26,00,
     1		00,00,00,00,00,00,00,00,00,00,00,00,00,
     1		03,04,05,06,08,09,11,12,13,14,16,22,25,26,
     1		00,00,00,00,00,00,00,00,00,00,00,00,00,
     1		02,03,04,05,06,08,09,11,12,13,14,16,22,25,
     1		26,00,00,00,00,00,00,00,00,00,00,00,00,
     1		02,03,04,05,06,08,09,11,12,13,14,16,22,24,
     1		25,26,00,00,00,00,00,00,00,00,00,00,00,
     1		02,03,04,05,06,07,08,09,11,12,13,14,16,22,
     1		24,25,26,00,00,00,00,00,00,00,00,00,00,
     1		02,03,04,05,06,07,08,09,11,12,13,14,16,22,
     1		23,24,25,26,00,00,00,00,00,00,00,00,00,
     1		01,02,03,04,05,06,07,08,09,11,12,13,14,16,
     1		22,23,24,25,26,00,00,00,00,00,00,00,00,
     1		01,02,03,04,05,06,07,08,09,11,12,13,14,15,
     1		16,22,23,24,25,26,00,00,00,00,00,00,00,
     1		01,02,03,04,05,06,07,08,09,10,11,12,13,14,
     1		15,16,22,23,24,25,26,00,00,00,00,00,00,
     1		01,02,03,04,05,06,07,08,09,10,11,12,13,14,
     1		15,16,17,22,23,24,25,26,00,00,00,00,00,
     1		01,02,03,04,05,06,07,08,09,10,11,12,13,14,
     1		15,16,17,22,23,24,25,26,27,00,00,00,00,
     1		01,02,03,04,05,06,07,08,09,10,11,12,13,14,
     1		15,16,17,21,22,23,24,25,26,27,00,00,00,
     1		01,02,03,04,05,06,07,08,09,10,11,12,13,14,
     1		15,16,17,20,21,22,23,24,25,26,27,00,00,
     1		01,02,03,04,05,06,07,08,09,10,11,12,13,14,
     1		15,16,17,19,20,21,22,23,24,25,26,27,00,
     1		01,02,03,04,05,06,07,08,09,10,11,12,13,14,
     1		15,16,17,18,19,20,21,22,23,24,25,26,27/
     	data rgb25/ 
     1		50,    0, 170,
     1		50, 0, 255,
     1		0, 0, 255, 
     1		0, 50, 200,
     1	 	0, 100, 255,
     1	 	0, 191,255,
     1	 	50, 205, 100,
     1	 	5, 255, 80,
     1	 	5, 220, 0,
     1	 	0,   190, 0,
     1	 	124, 252, 0,
     1	 	255, 255,   0,
     1	 	255, 215, 0,
     1	 	255, 165,   0,
     1	 	250,   150,   0,
     1	 	250,   120,   0,
     1	 	255, 105, 0,
c    1	 	240, 90, 155,
     1	 	191, 128, 255,
     1	 	191, 92, 255,
     1	 	195, 48, 255,
     1	 	180,  32, 240,
     1	 	210, 1, 255,
     1	 	255,   0, 255,
     1	 	240,   0, 150,
     1	 	255,   10,   100,
     1	 	230,   0,   0,
     1	 	200,   0,   0/
c    	read(100,*)rgb25
	nocl=n-3
	nn=nocl+1
c
c	first background color is white
	if (kkk.eq.1) call gscr(iwkid,0,1.,1.,1.)
c	first foreground color is black
	call gscr(iwkid,1,0.,0.,0.)
c	second foreground color is gray
	call gscr(iwkid,2,1.,1.,1.)
	if (bw.eq.1) then
	if (nocl.gt.26) then
		print *
		print *
		print *,'colours >25'
		stop
	else 
        	do i=1,nn
			kk=mtrx(i,nn-2)
        		r=rgb25(1,kk)
			g=rgb25(2,kk)
        		b=rgb25(3,kk)
        		call gscr(iwkid,i+2,r/255.,g/255.,b/255.)
        	end do
	end if
	end if

      if (bw.eq.2) then
	do i=1,nocl+1
      	 fi=1-(i/float(nocl+5))
    	 r=fi
	 b=fi
	 g=fi
         call gscr(iwkid,i+2,r,g,b)
	end do 
      end if
      call gscr(iwkid,30,.5,.3,.3) ! for land 
      close(10)
c
      return
      end
      subroutine clr (n,iwkid,kkk)

	nocl=n-3
c
c	first background color is white
	if (kkk.eq.1) call gscr(iwkid,0,1.,1.,1.)
c	first foreground color is black
	call gscr(iwkid,1,0.,0.,0.)
c	second foreground color is gray
	call gscr(iwkid,1,0.,0.,0.)
c	call gscr(iwkid,2,1.,1.,1.)
      call gscr(iwkid,30,0.5,0.5,0.5) ! for land in ccpcm
c
      return
c
      end


	subroutine ngrd(xreg,yreg,zreg,mreg,nreg,xmn,xmx,ymn,ymx,izmn
     1	,izmx,sumw,gy,gx,sumz,thv,ipwr,wf,gm)
c	weighted average gridding
        integer mreg, nreg
        real xreg(mreg), yreg(nreg), zreg(mreg,nreg)
        character *24 iff,gm*1,std*1
	character *1 mnthplt
        dimension sumw(mreg,nreg),vrbl(20)
	dimension gy(nreg),gx(mreg)
        integer xmnlmt,ymnlmt,xmxlmt,ymxlmt,tprm
	real sumz(mreg,nreg)
	mnthplt='n' !n for all data together
	izmn=999999.
	izmx=-999999.
	rmnv=99999
	rmxv=-99999
	n=0
	nn=0
	vna=-888.8	!vl tbw whn nt avlbl
	print '(/a,$)',' Input file: '
	read '(a)',  iff
	if (mnthplt.eq.'y') then
       		print '(/a,$)',' Month tbc: '
        	read *,mntbc
	end if
	open (unit=55,status='old',file=iff)
	print *,'entd'
	print '(/a,$)',' Total No. of parameters tbr: '
	read *,tprm
	print '(/a,$)', 'x-axis/y-axis/value order:'
	read *, ix,iy,iv
cc	gm='w'  !b for box average; w for weighted mean
	print '(/a,$)',' Search radius: '
	read *,srad
	nx=mreg
	ny=nreg
	dx=(xmx-xmn)/(nx-1)
	dy=(ymx-ymn)/(ny-1)
c	initialise x & y grid mode:
	do i=1,ny
		gy(i)=float(i-1)*dy+ymn-dy/2
	end do
	do i=1,nx
		gx(i)=float(i-1)*dx+xmn-dx/2
	end do
	do i=1,nx
	do j=1,ny
		sumz(i,j)=0.
		sumw(i,j)=0.
	end do
	end do
	read (55,*)
	print *,
	print *,' First record skipped'
	print *,
  5	read(55,*,end=90)(vrbl(kk),kk=1,tprm)
	x=vrbl(ix)
	if (mnthplt.eq.'y') then
		mon=vrbl(3)
		if (mntbc.ne.mon) go to 5
	end if
	if (x.gt.180.) x=x-360.
	y=vrbl(iy)
	z=vrbl(iv)
        write (100,*)x,y,z
	n=n+1
	nn=nn+1
	if (x.lt.xmn.or.x.gt.xmx.or.y.lt.ymn.or.y.gt.ymx.or.z.eq.thv) 
     1	go to 5
	if (z.ne.thv) rmnv=min(rmnv,z)
	if (z.ne.thv) rmxv=max(rmxv,z)
	ncd=ncd+1
	if (gm.eq.'w') then
		xmnlmt=int((x-xmn-srad)/dx)-5
		ymnlmt=int((y-ymn-srad)/dy)-5
		xmxlmt=int((x-xmn+srad)/dx)+5
		ymxlmt=int((y-ymn+srad)/dy)+5
	end if
	if (gm.eq.'b') then	
		xmnlmt=int((x-xmn-dx)/dx)-1
		ymnlmt=int((y-ymn-dy)/dy)-1
		xmxlmt=int((x-xmn+dx)/dx)+1
		ymxlmt=int((y-ymn+dy)/dy)+1
	end if
	if (xmnlmt.lt.0) xmnlmt=1  !0
	if (ymnlmt.lt.0) ymnlmt=1  !0
	if (xmxlmt.gt.nx) xmxlmt=nx
	if (ymxlmt.gt.ny) ymxlmt=ny
	do i=xmnlmt,xmxlmt
	do j=ymnlmt,ymxlmt
		dx1=abs(x-gx(i))
		dy1=abs(y-gy(j))
		dist=sqrt(dx1*dx1+dy1*dy1)
		if (dist.lt.srad.and.gm.eq.'w') then


c			weighting functions:
c                       1: gaussian  2: Cressman   3: inverse square   4: Barnes

                        if (wf.eq.1) then
                          wfnc=exp(-1*(dist**ipwr)/(2*srad**ipwr))
                        else if (wf.eq.2) then
                          wfnc=(-dist**ipwr+srad**ipwr)/(dist**ipwr+
     1			  srad**ipwr)
                        else if (wf.eq.3) then
                          wfnc=1./(1.+dist**ipwr)
                        else if (wf.eq.4) then
                          wfnc=exp(-4*(dist**ipwr)/srad**ipwr)
                        end if
                        sumz(i,j)=sumz(i,j)+(z*wfnc)
			sumw(i,j)=sumw(i,j)+wfnc

		else if (dx1.lt.dx.and.x.lt.gx(i)+dx/2..and. 
     1			dy1.lt.dy.and.y.lt.gy(j)+dy/2..and.gm.eq.'b') 
     1                  then
			sumz(i,j)=sumz(i,j)+z ! (z/(dist**ipwr+1.))
			sumw(i,j)=sumw(i,j)+1 ! (1./(dist**ipwr+1.))
		end if
	end do
	end do
	go to 5
90	continue
	print '(//)'
	print *,' total/considered values: ',nn,ncd
	do i=1,nx
	rmdln=gx(i) !-dx/2 	!remove dy/2 if cntr pnt not rqrd
	do j=1,ny
		if (sumw(i,j).gt.0) then
			sumz(i,j)=sumz(i,j)/sumw(i,j)
		else
			sumz(i,j)=vna
		end if
		if (sumz(i,j).gt.izmx.and.sumw(i,j).gt.0)izmx=sumz(i,j)
		if (sumz(i,j).lt.izmn.and.sumw(i,j).gt.0)izmn=sumz(i,j)
		rmdlt=gy(j) !-dy/2 
		yn=ymn-dy/2.
		xn=xmn-dx/2.	
		xreg(i)=rmdln
		yreg(j)=rmdlt
		zreg(i,j)=sumz(i,j)
	end do
	end do 
	print *
	print *,' actl  min/max: ',rmnv,rmxv
	print *,' fnl   min/max: ',izmn,izmx
	close (55)
	return
	end



      subroutine landfill(rlat1,rlon1,rlat2,rlon2)
        common /clrcmn/ iam(250000)
        dimension xcs(10000),ycs(10000)
        dimension iai(10),iag(10)
        dimension ioc(14)
        external colram,colrln
        data ioc / 6,2,5,12,10,11,1,3,4,8,9,7,13,14 /
        call gsfais (1)
        call mapstc ('ou','co')
        call maproj ('ce',0.,0.,0.)
        call mapset ('co',rlat1,rlon1,rlat2,rlon2)
        call mapint
        call arinam (iam,250000)
        call mapbla (iam)
        call arpram (iam,0,0,0)
        call arscam (iam,xcs,ycs,10000,iai,iag,10,colram)
      return
      end
c
      subroutine colram (xcs,ycs,ncs,iai,iag,nai)
      dimension xcs(*),ycs(*),iai(*),iag(*)
      idmap=-1
      idcont=-1

      do 10, i=1,nai
         if (iag(i).eq.1) idmap=iai(i)
         if (iag(i).eq.3) idcont=iai(i)
 10   continue
      if (iai(1).ge.0.and.iai(2).ge.0) then
       if (mapaci(idmap).ne.1) then
         call gsfaci (30)
         call gfa (ncs-1,xcs,ycs)
       end if
      end if
      return
      end
      subroutine colrln (xcs,ycs,ncs,iai,iag,nai)
      dimension xcs(*),ycs(*),iai(*),iag(*)
        if (iai(1).ge.0.and.iai(2).ge.0) then
          itm=max0(iai(1),iai(2))
          if (mapaci(itm).eq.1) then
            if (ncs.gt.150) print * , 'colrln - ncs too big - ',ncs
            call gpl (ncs,xcs,ycs)
          end if
        end if
        return
      end