MODULE domain
        implicit none
        save
        real :: XW,XE,YS,YN
      END MODULE domain
      PROGRAM contour90
      USE domain
      implicit none
!
!         READ ME
!
!      your source code to generate the 2-d arrays goes between the 
!        first pair of lines that contain c repeated many times. To
!        change the calls to the plot generation, make mods between
!        the next pair of lines containing c repeated many times.
!
!      run the program by:
!         ncargf90 contour90.f
!         ./a.out
!       
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!      For this example, arrays PSIP and PSI are created.
!
      integer, PARAMETER :: NY=21,Nx=41
!      NY=GRID PTS IN Y, NX GRID PTS IN X, 
      real :: Y(NY),X(NX),PSIP(NX,Ny),PSI(NX,Ny)
      real :: pi,g,cint
      integer :: i,j
!
      CHARACTER (len=8) :: CHRS
      PI=acos(-1.)
      eye : DO I=1,NX
      X(I)=FLOAT(I-1)*0.1-2.
      jay : DO J=1,NY
      Y(J)=FLOAT(J-1)*0.1-1.
      PSIP(I,J)=17.*EXP(-X(I)**2-Y(J)**2)
      PSI(I,J)=27.2*SIN(0.5*PI*X(I))*COS(PI*Y(J)/2.-X(I)/3.)
      end do jay
      end do eye
!
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!
!      NOW PLOT PSIP and PSI
!
!      start plot commands
      CALL OPNGKS
      g=1.
!      input ranges, neg values may be needed due to constraints of the
!        call to subroutine set in subroutine ploty (see note below)
!      these ranges are also used to label axes.
!      left & right edges of "x" range
      XW=X(1)
      XE=X(NX)
!      bottom & top edges of "y" range
      YS=Y(1)
      YN=Y(NY)
!      note: the set call requires: xwd<xed, ysd<ynd.  
!       If you have a vertical x-section with pressure
!        as a vertical coordinate, you may want to use negative values 
!        for ysd and ynd: i.e. ysd=-1000  ynd=-200 will put "-200" label
!        at the top of the axis and "-1000" at the bottom.  Of course...
!      you could always modify subroutine ploty to suit specific needs...
!
!      create a plot title 8 characters long
      WRITE(CHRs,'(8hPSI SB P)')
!      cint is the contour interval
      cint=1.5
!      create a plot and advance to the next page when done
      call ploty(PSIP,nx,nx,ny,g,chrs,cint)
!      start creation of second plot, etc.
      WRITE(CHRs,'(8h   PSI   )')
!      using cint=0. tells ncar graphics to choose the contour interval
      cint=0.
      call ploty(PSI,nx,nx,ny,g,chrs,cint)
!      end plot commands
      CALL CLSGKS
!
      stop
      end program contour90
!
!------------------------------------------------------------------
!
      SUBROUTINE PLOTY(PLTA,NDIM,NX,NY,G,CHRS,CINT)
      USE domain
!      While this sub can be used as is, you may wish to improve the
!       appearance of your plots by modifying it as you become 
!       familiar with ncar graphics.
!      Sub plots array plta.  
!      Plot the first nx columns of the array plta which
!       is dimensioned to have ndim columns.  plot the first ny rows
!       which must be less than or equal to the dimensioned size of
!       plta in the main program.
!      chrs contains a short label that is passed through
!      cint is the contour interval which is also passed through
!      g is an optional variable numeric label which could be encoded
!        using a write statement and plotted using a plchhq call.
!        (the latter is included, but commented out).
      implicit none
      integer :: ndim,nx,ny
      integer :: IX,JY,ISIZ,JSIZ,IOR,IG,JG,IGTYP,nx3,ny3
      real :: xdif,ydif,g,cint
      real :: PLTA(NDIM,NY)
      CHARACTER (len=8) :: CHRS,chr
!      CHR: EXAMPLE FOR ENCODING A VARIABLE USED IN A GRAPH LABEL (OPT.)
      WRITE(CHR,'(5htime=,F3.1)')G
!      SET SPECIAL VALUE FOR CONTOURING TO AVOID GRID PTS WITH SPV VALUE
      CALL CPSETR('SPV',1.E36)
!      GSCLIP(0) CALL SO THAT WHAT IS OUTSIDE SET CALL RANGE IS ALLOWED
      CALL GSCLIP(0)
!      SET CALL USED TO DEFINE THE LOCATION OF PLOT ON MONITOR SCREEN
!      CALL SET(.2,.8,.25,.75,0.,1.,0.,1.,1)
!      NEXT 6 LINES NEEDED TO GET CNTR PLT TO FIT INSIDE SET RANGES
!      CALL SET(.2,.7,.25,.75,XW,XE,YS,YN,1)
      CALL SET(.1,.9,.35,.75,XW,XE,YS,YN,1)
      CALL CPSETI('SET',0) !1=CONPACK CALLS SET, =0 NO SET CALL
      CALL CPSETR('XC1',XW)
      CALL CPSETR('XCM',XE)
      CALL CPSETR('YC1',YS)
      CALL CPSETR('YCN',YN)
!
      CALL CPSETR('SPV',1.E36) ! SPECIAL VALS TO AVOID WITH CNTR LINE
      CALL CPSETR('CWM',1.8)  ! SCALE FACTR TO MULT ALL DEFAULT LABL SIZES
!
!      CALL CPSETI('LLP',-2)  !QUALITY OF LABL PLACEMNT 1=WORST 3=BEST
      CALL CPSETI('LIS',4)  !INTRVL BETWEEN LABLED LINES
!      SECOND TEST
      CALL CPSETI ('LLP - LINE LABEL POSITIONING',3)! USE FORTRAN FORMAT 
!     +   SCIENTIFIC NOTATION IF VALUE =3
!
      CALL CPSETI ('HLB - HIGH/LOW LABEL BOX FLAG',1) !ASK 4 BOX AROUND
!     +   HI/LO LABEL
!
      CALL CPSETI ('HLO - HIGH/LOW LABEL OVERLAP',4) !NO LABL OUTSIDE SET
!     +   RANGE
!
!      CALL CPSETR('HLS - HI/LO LABEL SIZ',.012)  ! DEFAULT IS .012
!     CALL CPSETI ('NCL - CONTOUR LEVEL SELECTION FLAG', NCL)
!
      CALL CPSETR ('PC1 - PENALTY SCHEME CONSTANT 1',1.5)
      CALL CPSETR ('T2D - TENSION ON THE 2D SMOOTHER',1.)
!      NEXT LINE FINALLY DRAWS THE CONTOUR LINES
      CALL CPCNRC(PLTA,NDIM,NX,NY,0.,0.,CINT,-1,0,-1)
!      CALL CPCNRC(PLTA,NDIM,NX,NY,0.,0.,0.,-1,0,-1)
!
!      NOW ADD VARIOUS LABELS AND TICK MARKS
!      AXES/PERIMETER LABELS
!      IF IX OR JY =1 USE RT EDGE OR TOP EDGE FOR LABELS
      IX=0
      JY=1
!      BORDER LABELS LETTER SIZE
      ISIZ=12
      JSIZ=12
!      IOR=0 FOR HORZ LABLS,  =1 FOR VERTICAL ORIENTED LABELS,
      IOR=0
!      LABMOD ALSO SPECIFIES FORMAT OF LABELS, SIZE, ETC.
      CALL LABMOD('(F5.1)','(F5.1)',0,0,ISIZ,JSIZ,IX,JY,IOR)
!      CALL LABMOD('(F5.3)','(F5.3)',0,0,ISIZ,JSIZ,IX,JY,IOR)
!      NOW DETERMINE TYPE OF BORDER AROUND CONTOURS, E.G. A GRID, A BOX
!       WITH TICK MARKS, ETC.  SPECIFY IF TICK MARKS HAVE LABELS, ETC.
!      PLOT BORDER/LABELS,  1ST 4 ARGS OF GRIDAL LIKE OLD SUB PERIM
!      IG,JG =0 THEN LIKE PERIM, =1 THEN LIKE PERIML
      IG=0
      JG=0
      IF(G.GE.1.) THEN
      IG=1
      JG=1
      END IF
!      IGTYP =0 FOR GRID, =5 PERIMETR, =10 AXES, OTHER COMBS POSS.
      IGTYP=5
!      LAST 2 ARGS GRIDAL ARE INTERSECTN OF AXES IN X,Y VALUES
!      CALL GRIDAL(NX-1,0,NY-1,0,IG,JG,IGTYP,0.,0.)
      NX3=(NX-1)/3
      NY3=(NY-1)/3
!      CALL GRIDAL(NX3,3,NY3,3,IG,JG,IGTYP,0.,0.)
      XDIF=(XE-XW)/3.
!      YDIF=(YN-YS)/5.
      YDIF=(YN-YS)/2.
      NX3=IFIX(XDIF)+1
      NY3=IFIX(YDIF)+1
!      CALL GRIDAL(NX3,3,NY3,5,IG,JG,IGTYP,0.,0.)
      CALL GRIDAL(NX3,3,NY3,2,IG,JG,IGTYP,0.,0.)
!      NOW ADD GRAPH TITLE AND OTHER LABELS
!      THE SET CALL IS FOR PLACEMENT OF THESE LABELS
      CALL SET(.05,.9,.05,.9,0.,1.,0.,1.,1)
!      GRAPH TITLE CHRS PASSED THROUGH 
      CALL PLCHHQ(.5,0.05,CHRS,12.,0.,0.)
!      PLOT ANOTHER GRAPH LABEL CHR (OPT.)
      CALL PLCHHQ(.5,0.1,CHR,12.,0.,0.)
      CALL FRAME
      RETURN
      END subroutine ploty