         PROGRAM PLOT

C Plotting routine using GKS NCAR Graphics for
C   MAPS/RUC 40-km grids on Unix

C
      INTEGER NX_P,NY_P

      REAL LAT_LL_P, LON_LL_P, LAT_UR_P, LON_UR_P
  
c       PARAMETER (NX_P=151, NY_P=113, NVAR_P=7, NZTN_P=40)
c       PARAMETER (NX_P=301, NY_P=225, NVAR_P=7, NZTN_P=40)
        PARAMETER (NX_P=451, NY_P=337, NVAR_P=7, NZTN_P=50)
        parameter (nvar3 = nvar_p+7)
        parameter (nvar_h = 46)
        parameter (nvar_g = 88)
        parameter (nvar_ms = 44)
C   --- number of surface variables - NCEP
        INTEGER NVARSFC_NCEP_P
        PARAMETER (NVARSFC_NCEP_P=25)


        parameter (nvar_s =  9)
        parameter (maxlev = 50)
        parameter (maxlevp= 37)

      integer p_p,m_p,theta_p,pc_p,u_p,v_p
      parameter (p_p=1,m_p=2,theta_p=3,pc_p=4,u_p=5,v_p=6)

C ----   3-d grids
	REAL G3(NX_P,NY_P,NVAR3,maxlev)
       REAL G3a(NX_P,NY_P,maxlev,nvar3)
       REAL G3p(NX_P,NY_P,maxlevp,6)
        real g4m(nx_p,ny_p,nvar_ms)
        real diffpc(nx_p,ny_p)
        real snratio(nx_p,ny_p)


	equivalence (g3a(1,1,1,1),g3p(1,1,1,1))

cREAL G1(NX_P,NY_P,maxlev-1,NVAR_P)
        real g (nx_p,ny_p,nvar3)
        real g2(nx_p,ny_p,nvar_p)
C ----   2-d grids (MSP files)
        real g4(nx_p,ny_p,nvar_g)
        real gt(nztn_p)
        real gp(maxlev)
        real pv(nx_p,ny_p)
        real pcloud(nx_p,ny_p)
        real sumq (nztn_p)
        real absvor(nx_p,ny_p)
        real spd(nx_p,ny_p)
        real hold (nx_p,ny_p)

	CHARACTER*9 ATIME
	CHARACTER*9 bTIME

        character*80 dir

	REAL CORIO(NX_P,NY_P), SCALE(NX_P,NY_P)
	REAL UWAVE,VWAVE

      CHARACTER*48 LABEL(18) 
      CHARACTER*48 LABELPV
      CHARACTER*51 timelabel(2)
      character*48 LABELV
      character*48 LABELOMG
      character*15 TYPELABEL(3,3)
c      INTEGER*4 NETVARID2(nvar_h)
         INTEGER*4 NETVARID3(NVAR_P)
         INTEGER*4 NETVARIDs(NVAR_s)
 
C  No 2-d fields to be read from NMC data in this program
        DATA NETVARIDx /-99999/
c     DATA NETVARID2 /13,14,18,19,20,22,23,24,27,28,29,30,31,32,
c    1      33/
      data netvarids /7,9,4,5,1,2,3,6,8/

C  $$$  JMB:  Introduced this 9 Mar 95
C  $$$  TBH:  BUGFIX in modelin.f:
C  $$$  TBH:  Re-ordered NETVARID3 to match contents of netCDF file.
C        DATA NETVARID3 /5,3,6,4,1,2/
c       DATA NETVARID3 /4,3,5,6,1,2/
        DATA NETVARID3 /4,3,5,6,1,2,12/
       integer netvarid2p(1)
       integer netvarid3p(6)
       data netvarid2p /-99999/
       data netvarid3p /6,1,2,3,4,5/

       integer ivarx(6)
c      data ivarx /6,1,2,3,4,5,7/
c      data ivarx /1,2,3,4,5,6/
       data ivarx /6,1,2,3,4,5/

        character*60 varname2(nvar_h)
        character*60 varnamex
        character*60 varnames(nvar_s)
        character*60 varname3(NVAR3)
        character*60 varname3_ruc(NVAR3)
        character*60 varname3p(6)
        character*60 varname2_ms(NVAR_ms)
        character*60 varname_sfc_ncep(NVARSFC_ncep_P)
        character*60 varname_sfc_ruc(NVARSFC_ncep_P)




        data varname3_ruc /
     1  'pressure',
     1  'geopotential height',
     1  'virtual potential temperature',
     1  'water vapor mixing ratio',
     1  'u-component of wind',
     1  'v-component of wind',
     1  'vertical velocity',
     1  'cloud water mixing ratio',
     1  'rain water mixing ratio',
     1  'ice mixing ratio',
     1  'snow mixing ratio',
     1  'graupel mixing ratio',
     1  'number concentration for ice particles'
     1  ,'turbulent kinetic energy'
     1  /

        data varname3 /
     1  'pressure - b',
     1  'height - b',
     1  'virtual potential temperature - b',
     1  'water vapor mixing ratio - b',
     1  'u-component of wind - b',
     1  'v-component of wind - b',
     1  'vertical velocity - b',
     1  'cloud water mixing ratio - b',
     1  'rain water mixing ratio - b',
     1  'ice mixing ratio - b',
     1  'snow mixing ratio - b',
     1  'graupel mixing ratio - b',
     1  'number concentration for ice particles - b'
     1  ,'turbulent kinetic energy - b'
     1  /

        data varname3p /
     1  'vertical velocity - p',
     1  'height - p',
     1  'temperature - p',
     1  'relative humidity - p',
     1  'u-component of wind - p',
     1  'v-component of wind - p'
     1  /


        data varname2 /
     1   'MAPS mean sea level pressure '
     1  ,'soil temperature at surface'
     1  ,'sensible heat flux'
     1  ,'latent heat flux'
     1  ,'net longwave radiation at surface '
     1  ,'precipitation rate'
     1  ,'resolvable scale precipitation - 3h'
     1  ,'sub-grid scale precipitation - 3h'
     1  ,'precipitable water'
     1  ,'pressure at tropopause'
     1  ,'potential temperature at tropopause'
     1  ,'u-component of wind at tropopause'
     1  ,'v-component of wind at tropopause'
     1  ,'convective available potential energy'
     1  ,'convective inhibition'
     1  ,'soil temperature at level 1 below surface'
     1  ,'soil temperature at level 2 below surface'
     1  ,'soil temperature at level 3 below surface'
     1  ,'soil temperature at level 4 below surface'
     1  ,'soil temperature at the bottom'
     1  ,'soil volumetric moisture content at surface'
     1  ,'soil volumetric moisture content at level 1'
     1  ,'soil volumetric moisture content at level 2'
     1  ,'soil volumetric moisture content at level 3'
     1  ,'soil volumetric moisture content at level 4'
     1  ,'soil volumetric moisture content at the bottom'
     1  ,'direct evaporation from bare soil'
     1  ,'evaporation of canopy water'
     1  ,'evapotranspiration flux'
     1  ,'water condensation rate near surface'
     1  ,'rate of water dropping from canopy to ground'
     1  ,'net shortwave radiation at surface'
     1  ,'snow accumulation'
     1  ,'snow depth'
     1  ,'surface runoff'
     1  ,'sub-surface runoff'
     1  ,'canopy water'
     1  ,'snow temperature - snow or top soil level'
     1  ,'temperature inside snow or top soil level'
     1  ,'water vapor mixing ratio at surface'
     1  ,'water equivalent of snow depth'
     1  ,'snow density'
     1  ,'temperature, 2m fixed height'
     1  ,'water vapor mixing ratio, 2m fixed height'
     1  ,'u-component of wind, 10m fixed height'
     1  ,'v-component of wind, 10m fixed height'


     1   /

        data varname2_ms /
     1   'MAPS mean sea level pressure '
     1  ,'earth surface pressure'
     1  ,'earth surface height'
     1  ,'soil temperature at surface'
     1  ,'total soil moisture content'
     1  ,'maximum possible total soil moisture content'
     1  ,'water equivalent of snow depth'
     1  ,'2 meter temperature'
     1  ,'2 meter specific humidity'
     1  ,'10 meter component u-wind'
     1  ,'10 meter component v-wind'
     1  ,'total precipitation accumulated'
     1  ,'total convective precipitation accumulated'
     1  ,'snow accumulation'
     1  ,'snow melt accumulated'
     1  ,'surface runoff accumulated'
     1  ,'time averaged latent heat flux'
     1  ,'time averaged sensible heat flux'
     1  ,'time averaged soil heat flux'
     1  ,'time averaged snow phase change heat flux'
     1  ,'surface drag coefficient'
     1  ,'sensible heat flux','latent heat flux'
     1  ,'convective available potential energy'
     1  ,'precipitable water'
     1  ,'low-level cloud cover percentage'
     1  ,'mid-level cloud cover percentage'
     1  ,'high-level cloud cover percentage'
     1  ,'time averaged downward short wave rad. flux'
     1  ,'time averaged downward long wave rad. flux'
     1  ,'time averaged upward short wave rad. flux'
     1  ,'time averaged upward long wave rad. flux'
     1  ,'surface albedo','soil temperature at level 1 below surface'
     1  ,'soil temperature at level 2 below surface'
     1  ,'soil temperature at level 3 below surface'
     1  ,'soil temperature at level 4 below surface'
     1  ,'soil volumetric moisture content at surface'
     1  ,'soil volumetric moisture content at level 1'
     1  ,'soil volumetric moisture content at level 2'
     1  ,'soil volumetric moisture content at level 3'
     1  ,'soil volumetric moisture content at level 4'
     1  ,'potential temperature at 10 m'
     1  ,'specific humidity at 10 m' 
     1    /


        data varnames /
     1   'MAPS mean sea level pressure'
     1  ,'surface pressure at mini topography'
     1  ,'potential temperature'
     1  ,'dew point temperature'
     1  ,'u-component of wind - s'
     1  ,'v-component of wind - s'
     1  ,'3-hr pressure tendency - s'
     1  ,'dew point depression'
     1  ,'altimeter setting'
     1   /

        data varname_sfc_ncep /
     1  'MAPS mean sea level pressure'
     1 ,'altimeter setting'
     1 ,'3-hr pressure tendency - s'
     1 ,'potential temperature'
     1 ,'dew point temperature'
     1 ,'dew point depression'
     1 ,'u-component of wind'
     1 ,'v-component of wind'
     1 ,'temperature'
     1 ,'specific humidity'
     1 ,'equivalent potential temperature'
     1 ,'convective available potential energy'
     1 ,'convective inhibition'
     1 ,'lifted index'
     1 ,'best lifted index'
     1 ,'resolvable scale precipitation - 3h'
     1 ,'sub-grid scale precipitation - 3h'
     1 ,'resolvable scale precipitation - 1h'
     1 ,'sub-grid scale precipitation - 1h'
     1 ,'categorical rain'
     1 ,'categorical freezing rain'
     1 ,'categorical ice pellets'
     1 ,'categorical snow'
     1 ,'snow accumulation'
     1 ,'helicity'
     1   /

        data varname_sfc_ruc /
     1  'MAPS mean sea level pressure'
     1 ,'altimeter setting'
     1 ,'3-hr pressure tendency - s'
     1 ,'potential temperature'
     1 ,'dew point temperature'
     1 ,'dew point depression'
     1 ,'u-component of wind, 10m fixed height'
     1 ,'v-component of wind, 10m fixed height'
     1 ,'temperature, 2m fixed height'
     1 ,'specific humidity'
     1 ,'equivalent potential temperature'
     1 ,'convective available potential energy'
     1 ,'convective inhibition'
     1 ,'lifted index'
     1 ,'best lifted index'
     1 ,'resolvable scale precipitation - 3h'
     1 ,'sub-grid scale precipitation - 3h'
     1 ,'resolvable scale precipitation - 1h'
     1 ,'sub-grid scale precipitation - 1h'
     1 ,'categorical rain'
     1 ,'categorical freezing rain'
     1 ,'categorical ice pellets'
     1 ,'categorical snow'
     1 ,'snow accumulation'
     1 ,'helicity'
     1   /


        data varnamex /
     1   '      '
     1   /


      do i=1,18
        label(i) = '                                                '
      end do
     
      labelpv = '                                                '
      labelv  = '                                                '
      labelomg= '                                                '
      timelabel(2) =
     1       '                                                   '

      LABEL(1)(1:12) = 'PRES    (MB)'
      LABEL(2)(1:12) = 'MONT/G  (DM)'
      LABEL(3)(1:12) = 'THETA-V (K) '
c     LABEL(4)(1:16) = 'COND PRES (MB)  '
      LABEL(4)(1:16) = 'WV MIX RAT -G/KG'

      LABEL(5)(1:12) = 'WIND MAG M/S'
      LABEL(6)(1:12) = 'DIV  E-6  /S'
      LABEL(7)(1:16) = 'STRETCH  DEF /S '
      LABEL(8)(1:16) = 'SHEAR    DEF /S '
      LABEL(9)(1:20) = 'COND PRES DEF (MB)  '
      label(10)(1:20)= 'VERT VEL (MICB/S)   '
      label(11)(1:20)= 'CLD WT MX RAT (G/KG)'
      label(12)(1:20)= 'RN  WT MX RAT (G/KG)'
      label(13)(1:20)= 'ICE    MX RAT (G/KG)'
      label(14)(1:20)= 'SNOW   MX RAT (G/KG)'
      label(15)(1:20)= 'GRPEL  MX RAT (G/KG)'
      label(16)(1:20)= 'ICE PT NO CONC      '
      label(17)(1:20)= 'TURB KIN ENGY (J/KG)'
      label(18)(1:20)= 'CLOUD TOP PRES (MB) '
	
 
      LABELPV(1:20)  = 'POT VORT E-5 K/MB/S '
      LABELV(1:16)   = 'ABS VORT E-5 /S '
      LABELOMG(1:16) = 'OMEGA   MICB/S  '
	DATA TYPELABEL/'MAPS Analysis  ',
     1                 'MAPS Forecast  ',
     1                 'AVN  Forecast  ',
     1                 'RUC  Analysis  ',
     1                 'RUC  Forecast  ',
     1                 'NGM  Forecast  ',
     1                 'MAPSFD Analysis',
     1                 'MAPSFD Forecast',
     1                 'NGM FD Forecast'/

c     polat = 90.
c     polon = -105.
c     lat_ll_p = 22.83730717
c     lon_ll_p = -120.49050192
c     lat_ur_p = 45.99867568
c     lon_ur_p = -60.82944145

      polat = 25.
      polon = -95.
c     lat_ll_p = 22.83730717
      lat_ll_p = 16.2810
      lon_ll_p = -126.1378
      lat_ur_p = 55.4818
      lon_ur_p = -57.3794

 
      IC = 0


      CALL HLEV (GT)

c     call gflas3(1)
c     call frame

C ***  Read in Coriolis parameter and map-scale factor
c     open (unit=9,name='/home/mtn3/benjamin/plot/corsc.asc',
c    1       form='formatted',type='old')
c     read (9,109) corio
c     read (9,109) scale
c09   format (8G10.3)
c     close (9)
C *** READ IN GRDPARM FILE FOR CORIOLIS AND MAP-SCALE FACTOR ARRAYS
         OPEN (UNIT=55, 
     1   FILE='/usr/nfs/common/rapb/13km/prmcmp/parms.mps',
c     1   FILE='/usr/nfs/common/rapb/20km/fslparms/parms.mps',
     1   STATUS='OLD',
     1  FORM='UNFORMATTED')
        READ (55) corio,scale
        CLOSE (55)
        print *,' Read in GRDPARM OK'


      type *,' Soundings and cross-sections'
      type *,'     now available from hybrid-b data.'
      type *,'  See later instructions'

C-----------------
C  Read in data
C-----------------
880   continue

        TYPE *,' ENTER ATIME'

	ACCEPT 10, ATIME
10      FORMAT (A9)
        btime  = atime

        imr = 1


        iaf = 1
        type *,' Enter 1 = hybrid-a (RT GRIB only)'
        type *,'       2 = hybrid-b (RT GRIB or NetCDF)'
        type *,'            Use this for skew-t diagrams from grids!'
        type *,'       3 = 25-mb isobaric (RT GRIB or MAPS DB)'
        type *,'       4 = 2-d fields from 3-d grids (trop,precip.etc.)'
        type *,'                    (RT NetCDF or MAPS DB)'
        type *,'       5 = surface grids'
        type *,'       6 = MORDS sfc grids'
        type *,'       7 = RUC sfc grids'

        accept 102, icord

        if (atime(8:9).ne.'00') iaf=2
        itype = iaf

        type *,' Enter storage format, 2 = GRIB'
        type *,'                       4 = NetCDF'
        accept 102, imap
          type *,' Enter 1/2 for MAPS anx/fcst, 3 for AVN,'
          type *,'         4 for TEST,5 for RUC2,6 - Backup'
          accept 102, ixx
        if (icord.ne.3) then
          if (ixx.eq.1) then
            dir = 'MAPS_DIR '
            if (imap.eq.2) dir = 'MAPS_DIRFG '
            if (icord.eq.7) dir = 'RUC_SFC '
          else if (ixx.eq.2) then
            dir = 'MAPS_DIRF '
            if (imap.eq.2) dir = 'MAPS_DIRFG '
            if (icord.eq.7) dir = 'RUC_SFC '
          else if (ixx.eq.3) then
            dir = 'ETA_HYBB '
            iaf = 3
          else if (ixx.eq.5) then
            dir='RUC_PRES '
            if(icord.eq.7) dir='RRUC_SFC '
            if(icord.eq.2) dir='RUC_DIR '
            if(icord.eq.4) dir='RUC_DIR '
            if(icord.eq.4.and.imap.eq.2) dir='RUC_DIRG '
          else if (ixx.eq.6) then
            dir='BACK_PRES '
            if(icord.eq.2) dir='BACK_DIR '
            if(icord.eq.4) dir='BACK_DIR '
          else 
            dir = 'MAPS_TEST '
          end if
          if (icord.eq.6) dir = 'MAPS_MORDS '
        else
          dir = 'MAPS_PRES '
          if (ixx.eq.3)dir = 'ETA_PRES '
          if (ixx.eq.4)dir = 'MAPS_TEST '
          if (ixx.eq.5)dir = 'RUC_PRES '
          if (ixx.eq.6)dir = 'BACK_PRES '
        end if

        print *,' dir',dir
C *** Set number of levels and number of variables depending on type
C         of file.
        if (icord.eq.1 .or. icord.eq.2) then
          nlev = nztn_p
          nvar = nvar3
        else if (icord.eq.3) then
          nlev = 37
          nvar = 6
          label(1)(1:16) = 'V.VEL (MICB/S)  '
          LABEL(2)(1:12) = 'HEIGHT  (DM)'
          LABEL(3)(1:12) = 'TEMP  (C)   '
          LABEL(4)(1:16) = 'REL HUM (%)     '
        else if (icord.eq.4) then
          nlev = 1
          nvar = nvar_h
        else if (icord.eq.5.or.icord.eq.7) then
          nlev = 1
          nvar = 9
        else if (icord.eq.6) then
          nlev = 1
          nvar = 44
        end if
          
C --- NetCDF data - 3-d
c       else if (imap.eq.4.and.icord.ne.5) then
c         call get_maps_netcdf (atime(1:7)//'00','00'//atime(8:9),
c    1         dir,g3,g4,istatus)

        if (icord.eq.2. or. icord.eq.4) then
         if (imap.eq.4) then 
c         call readnet (dir,atime(1:7)//'00','00'//atime(8:9),
c    1         nx_p,ny_p,nlev,1,nvar_h,nvar_p,netvarid2,netvarid3,
c    1         g4,g3a,istatus)
         if (ixx.ne.5) then
          call readgn (dir,atime(1:7)//'00','00'//atime(8:9),
     1         nx_p,ny_p,nlev,1,nvar_h,nvar3,varname2,varname3,
     1         g4,g3a,istatus)
         else
C -- use different Netcdf names for RUC hybrid-b file
          call readgn (dir,atime(1:7)//'00','00'//atime(8:9),
     1         nx_p,ny_p,nlev,1,nvar_h,nvar3,varname2,varname3_ruc,
     1         g4,g3a,istatus)
         end if

          if (istatus.ne.0 .or. g3a(2,2,2,1).eq.0.) then 
            print *,'istatus',istatus
            print *,' file not found, try again'
            go to 880
          end if


         else if (imap.eq.2) then
           call readgg(dir,atime(1:7)//'00','00'//atime(8:9),
     1         nx_p,ny_p,nztn_p,1,nvar_h,nvar3,varname2,varname3,
     1         g4,g3a,istatus)
          if (istatus.ne.0 .or. g3a(2,2,2,1).eq.0.) then
            print *,'istatus',istatus
            print *,' file not found, try again'
            go to 880
          end if
         end if
C --- Write out SOILSAVE data for soil temp and moisture to
C       assure continuity in soil cycling
        open (23,file='SOILSAVE',form='formatted',status='unknown')
        write (23,877) atime,atime(6:7)
 877     format (a9,2x,a2)
        write (23,878) 
     1      ((g4(i,j,2),i=1,nx_p),j=1,ny_p),
     1      (((g4(i,j,k),i=1,nx_p),j=1,ny_p),k=16,20),
     1      (((g4(i,j,k),i=1,nx_p),j=1,ny_p),k=21,26),
     1      ((g4(i,j,38),i=1,nx_p),j=1,ny_p),
     1      ((g4(i,j,34),i=1,nx_p),j=1,ny_p),
     1      ((g4(i,j,41),i=1,nx_p),j=1,ny_p),
     1      ((g4(i,j,42),i=1,nx_p),j=1,ny_p),
     1      ((g4(i,j,39),i=1,nx_p),j=1,ny_p),
     1      ((g4(i,j,40),i=1,nx_p),j=1,ny_p)
c 878     format (6(22F6.2/),19F6.2)
878     format ((27(11(3PE12.4)/),4(3PE12.4))/)

        close (23)
        open (23,file='SOILSAVE2',form='formatted',status='unknown')
        write (23,877) atime,atime(6:7)
        write (23,878) 
     1      ((g4(i,j,2),i=1,nx_p,2),j=1,ny_p,2),
     1      (((g4(i,j,k),i=1,nx_p,2),j=1,ny_p,2),k=16,20),
     1      (((g4(i,j,k),i=1,nx_p,2),j=1,ny_p,2),k=21,26),
     1      ((g4(i,j,38),i=1,nx_p,2),j=1,ny_p,2),
     1      ((g4(i,j,34),i=1,nx_p,2),j=1,ny_p,2),
     1      ((g4(i,j,41),i=1,nx_p,2),j=1,ny_p,2),
     1      ((g4(i,j,42),i=1,nx_p,2),j=1,ny_p,2),
     1      ((g4(i,j,39),i=1,nx_p,2),j=1,ny_p,2),
     1      ((g4(i,j,40),i=1,nx_p,2),j=1,ny_p,2)
        close (23)
        print *,' Write out new SOILSAVE'

        open (23,file='SFCQDIFF',form='unformatted',status=
     1    'unknown')
          do j=1,ny_p
          do i=1,nx_p
            hold(i,j)=(g4(i,j,44)-g3a(i,j,1,4))/g3a(i,j,1,4)
          end do
          end do
        write (23) hold
        close (23)

C
C     CALC MAX AND MIN OF F
C
      FMAX = -1.E8
      FMIN =  1.E8
      DO J=1,NY_P
        DO I=1,NX_P
          FMAX = AMAX1(FMAX,hold(i,j))
          FMIN = AMIN1(FMIN,hold(i,j))
      end do
      end do

        TYPE *,' qdiff MIN =',FMIN,' MAX =',FMAX

          do 754 k=1,nlev
          do 754 ivar = 1,nvar3
          do 754 j=1,ny_p
          do 754 i=1,nx_p
            g3(i,j,ivar,k) = g3a(i,j,k,ivar)
754       continue

        open (21,file='maxmin',form='formatted',type='unknown')
         do 7541 k=1,nlev
         do 7541 ivar = 1,nvar3
          qmax = -100000.
          qmin = 1000000.
          do j=1,ny_p
          do i=1,nx_p
            qmax = max(qmax,g3(i,j,ivar,k))
            qmin = min(qmin,g3(i,j,ivar,k))
          end do
          end do
          print *,k,qmax,qmin
         write (21,211) k,ivar,qmax,qmin
211      format('k,ivar',2i5,2e12.3)
7541      continue
         k=1
         do 7543 ivar = 1,nvar_h
          qmax = -100000.
          qmin = 1000000.
          do j=1,ny_p
          do i=1,nx_p
            qmax = max(qmax,g4(i,j,ivar))
            qmin = min(qmin,g4(i,j,ivar))
          end do
          end do
          print *,ivar,qmax,qmin
         write (21,211) k,ivar,qmax,qmin
7543      continue
        do 7542 k=1,nlev-1
          qmax = -100000.
          qmin = 1000000.
          do j=1,ny_p
          do i=1,nx_p
            dp = g3(i,j,1,k) - g3(i,j,1,k+1)
            qmax = max(qmax,dp)
            qmin = min(qmin,dp)
          end do
          end do
         write (21,211) k,' dp ',qmax,qmin
7542    continue

        close (21)

          do 755 k=1,nlev
          do 755 j=1,ny_p
          do 755 i=1,nx_p
C --- reduce pressure to mb from Pa
            g3(i,j,1,k) = g3(i,j,1,k)/100.
C --- set water vapor mixing ratio to g/kg
            g3(i,j,4,k) = g3(i,j,4,k)*1000.
C --- reduce ice particle no. conc. by factor of 10**8
            if (g3(i,j,13,k).gt.1.e2) then
               g3(i,j,13,k) = alog10(g3(i,j,13,k))
            else
               g3(i,j,13,k) = 0.
            end if
            do ivar=8,12
              g3(i,j,ivar,k) = 1000.*g3(i,j,ivar,k)-1.e-5
            end do
755       continue
          do 756 j=1,ny_p
          do 756 i=1,nx_p
            g4(i,j,1) = g4(i,j,1)/100.
            g4(i,j,10) = g4(i,j,10)/100.
            g4(i,j,6) = g4(i,j,6)*3600.
            g4(i,j,7) = g4(i,j,7)+g4(i,j,8)
            g4(i,j,28) = g4(i,j,28)
            g4(i,j,29) = g4(i,j,29)
            g4(i,j,33) = g4(i,j,33)*1000.
            g4(i,j,34) = g4(i,j,34)*1000.
            g4(i,j,35) = g4(i,j,35)*1000.
            g4(i,j,21) = -g4(i,j,21)
            g4(i,j,22) = -g4(i,j,22)
            g4(i,j,23) = -g4(i,j,23)
            g4(i,j,24) = -g4(i,j,24)
            g4(i,j,25) = -g4(i,j,25)
            g4(i,j,26) = -g4(i,j,26)
           diffpc(i,j)=g4(i,j,7)-g4(i,j,33)/10.
           if(g4(i,j,7).gt.0.01)
     1     snratio(i,j)=g4(i,j,33)/10./g4(i,j,7) 
          if(g4(i,j,34).le.20.)g4(i,j,34)=-25.
756       continue
        open (17,file='DIFFPCPN',form='unformatted',status='unknown')
           write (17) diffpc
           write(17) snratio
        close (17)
            

          do ivar=1,nvar_h
            print *,ivar,g4(60,60,ivar)
          end do
C-------------------------------------------
C --- 25-mb isobaric data
C-------------------------------------------
        else if (icord.eq.3) then
         if (imap.eq.4) then 
c         call readnet (dir,atime(1:7)//'00','00'//atime(8:9),
c    1         nx_p,ny_p,nlev,1,1,nvar,netvarid2p,varname3p,
c    1         g4,g3p,istatus)
          call readgn (dir,atime(1:7)//'00','00'//atime(8:9),
     1         nx_p,ny_p,nlev,1,1,nvar,varname2,varname3p,
     1         g4,g3p,istatus)
          if (istatus.ne.0) then
            print *,'istatus',istatus
            print *,' file not found, try again'
            go to 880
          end if

          else if (imap.eq.2) then
c          call readgg('MAPS_PRESG ',atime(1:7)//'00','00'//atime(8:9),
           call readgg(dir,atime(1:7)//'00','00'//atime(8:9),
     1         nx_p,ny_p,nlev,1,1,nvar,varnamex,varname3p,
C -- set to nvar3 again when TKE is written to GRIB
c    1         nx_p,ny_p,nlev,1,nvar_h,nvar3,varname2,varname3,
     1         g4,g3p,istatus)
          if (istatus.ne.0) then
            print *,'istatus',istatus
            print *,' file not found, try again'
            go to 880
          end if
         end if

         do ivar=1,nvar
           print *,'isobaric',g3p(15,15,7,ivar)
         end do

         if (imap.eq.4) then
          do 854 k=1,nlev
          do 854 ivar = 1,nvar
          do 854 j=1,ny_p
          do 854 i=1,nx_p
            g3(i,j,ivar,k) = g3p(i,j,k,ivar)
854       continue
         else if (imap.eq.2) then
          do 8541 k=1,nlev
          do 8541 ivar = 1,nvar
          do 8541 j=1,ny_p
          do 8541 i=1,nx_p
            g3(i,j,ivar,k) = g3p(i,j,k,ivarx(ivar))
8541      continue
         end if
          do 855 k=1,nlev
          do 855 j=1,ny_p
          do 855 i=1,nx_p
            g3(i,j,1,k) = g3(i,j,1,k)*10.
855       continue

C --- NetCDF data - surface
        else if (icord.eq.5) then
          call readgn ('MAPS_SFC ',atime(1:7)//'00','00'//atime(8:9),
     1         nx_p,ny_p,nlev,1,nvar_s,1,varnames,varnamex,
     1         g4,g3a,istatus)
c         call readnet ('MAPS_SFC ',atime(1:7)//'00','00'//atime(8:9),
c    1         nx_p,ny_p,1,1,nvar_s,1,netvarids,netvaridx,
c    1         g4,g3a,istatus)
          if (istatus.ne.0.) then
            print *,'istatus',istatus
            print *,' file not found, try again', atime
            go to 880
          end if
          do j=1,ny_p
          do i=1,nx_p
            g4(i,j,1) = g4(i,j,1)/100.
            g4(i,j,2) = g4(i,j,2)/100.
            g4(i,j,9) = g4(i,j,9)/100.
            g4(i,j,4) = g4(i,j,4)-273.15
            g4(i,j,3) = g4(i,j,3)-273.15
            g4(i,j,7) = g4(i,j,7)*10800./100.
          end do
          end do
        do ivar=1,nvar_s
          print *,k,g4(2,2,ivar)
        end do
       do ivar=1,nvar
         print *,' 2-d fields ',ivar,' ',g4(60,60,ivar)
       end do
C --- MORDS sfc data
      else if (icord.eq.6) then
          call readgg (dir,atime(1:7)//'00','99'//atime(8:9),
     1         nx_p,ny_p,nlev,1,nvar_ms,1,varname2_ms,varnamex,
     1         g4m,g3a,istatus)
          if (istatus.ne.0 .or. g4m(2,2,2).eq.0.) then 
            print *,'istatus',istatus
            print *,' file not found, try again'
            go to 880
          end if

       do ivar=1,nvar_ms
         print *,' MORDS-sfc',ivar,' ',g4m(120,120,ivar)
       end do

          do 757 j=1,ny_p
          do 757 i=1,nx_p
            g4m(i,j,1) = g4m(i,j,1)/100.
            g4m(i,j,2) = g4m(i,j,2)/100.
            g4m(i,j,4) = g4m(i,j,4) - 273.15
            g4m(i,j,5) = -g4m(i,j,5)
            g4m(i,j,6) = -g4m(i,j,5)/g4m(i,j,6)
            g4m(i,j,8) = g4m(i,j,8) - 273.15
            g4m(i,j,9) = g4m(i,j,9)*1000.
            g4m(i,j,12) = g4m(i,j,12)
            g4m(i,j,13) = g4m(i,j,13)
            g4m(i,j,14) = g4m(i,j,14)/10.
            g4m(i,j,38) = -g4m(i,j,38)/1000.
            g4m(i,j,39) = -g4m(i,j,39)/1000.
            g4m(i,j,40) = -g4m(i,j,40)/1000.
            g4m(i,j,41) = -g4m(i,j,41)/1000.
            g4m(i,j,42) = -g4m(i,j,42)/1000.
            if(g4m(i,j,7).le.8.)g4m(i,j,7)=-10.
757       continue
       do ivar=1,nvar_ms
         print *,' MORDS-sfc',ivar,' ',g4m(120,120,ivar)
       end do

C --- RUC sfc data
      else if (icord.eq.7) then
          if (imap.eq.4) then
          call readgn (dir,atime(1:7)//'00','00'//atime(8:9),
     1         nx_p,ny_p,nlev,1,nvarsfc_ncep_p,1,varname_sfc_ncep
     1         ,varnamex,
     1         g4m,g3a,istatus)
          if (ixx.eq.5) call readgn 
     1         (dir,atime(1:7)//'00','00'//atime(8:9),
     1         nx_p,ny_p,nlev,1,nvarsfc_ncep_p,1,varname_sfc_ruc
     1         ,varnamex,
     1         g4m,g3a,istatus)
          else if (imap.eq.2) then
          call readgg (dir,atime(1:7)//'00','00'//atime(8:9),
     1         nx_p,ny_p,nlev,1,nvarsfc_ncep_p,1,varname_sfc_ncep,
     1         varnamex,
     1         g4m,g3a,istatus)
          end if

c         if (g4m(2,2,1).eq.0.) then 
c           print *,'istatus',istatus,g4m(2,2,1)
c           print *,' file not found, try again'
c           go to 880
c         end if

          do 7571 j=1,ny_p
          do 7571 i=1,nx_p
            g4m(i,j,1) = g4m(i,j,1)/100.
            g4m(i,j,2) = g4m(i,j,2)/100.
            g4m(i,j,3) = g4m(i,j,3)*10800.
            g4m(i,j,5) = g4m(i,j,5) - 273.15

            g4m(i,j,9) = g4m(i,j,9) - 273.15
            g4m(i,j,10) = g4m(i,j,10)*1000.

            g4m(i,j,16) = g4m(i,j,16)+g4m(i,j,17)
            g4m(i,j,18) = g4m(i,j,18)+g4m(i,j,19)
7571       continue
       do ivar=1,nvarsfc_ncep_p
         print *,' RUC-sfc',ivar,' ',g4m(60,60,ivar)
       end do

        end if

c         call get_maps_sfc60_netcdf (atime,g4,istatus)
c         if (istatus.ne.0 .or. g4(2,2,1).eq.0.) then 
c           print *,' file not found, try again'
c           go to 880
c         end if
c       end if

C     if (imap.eq.2) then
C       print *,' GRIB data'
      if(icord.ne.6) then
        
        do k=1,nlev
          print *,k,(g3(2,2,ivar,k),ivar=1,nvar)
        end do
      endif
C     end if

882   continue

C --------------
C  Now decide what kind of plot
C --------------

      type *,' Horizontal plot/Skew-t/cross-section. Enter 0/1/2'
      accept 102, iskewt
      if (iskewt.eq.1) then
        call skewtdata (atime,g3,nx_p,ny_p,nvar,nlev)
      else if (iskewt.eq.2) then
        type *,
     1 ' x-section of data (1) or actual hybrid-b model levels (2)?'
        accept 102,itypxs
        if (itypxs.eq.1) then
         call hybxsec  (atime,g3,corio,scale,nx_p,ny_p,nvar,nlev)
        else if (itypxs.eq.2) then
         call xsmodel  (atime,g3,corio,scale,nx_p,ny_p,nvar,nlev)
        end if
      else 
C        !!  horizontal plots
      if (icord.eq.1 .or. icord.eq.2) then
        TYPE *,' ENTER LEVEL (1-50) FOR THETA 224K - 500K'
      else if (icord.eq.3) then
        TYPE *,' ENTER LEVEL (1-37) FOR PRES 1000-100mb' 
        type *,' 850mb = 7, 700mb = 13, 500mb = 21, 400mb=25...'
        do ilev = 1,nlev
          gp(ilev) = 1025.-25.*ilev
        end do
      end if

      if (icord.lt.4) then
        ACCEPT 102, K
102     FORMAT (I)
      end if

      if (icord.eq.1 .or. icord.eq.2) then
cencode (48,50,timelabel) atime,gT(k),TYPELABEL(ITYPE)
c  50	format (2x,'TIME= ',a9,'  THETA= ',f5.0,1x,'K ',2X,
c1	A4,10X)
	encode (51,50,timelabel(1)) TYPELABEL(ITYPE,imr),atime,gt(k)
   50	format (2x,a15,2x,'TIME= ',a9,'  THETA= ',f5.0,1x,'K ')
      else if (icord.eq.3.or.icord.eq.6) then
cencode (48,51,timelabel) atime,gp(k),TYPELABEL(ITYPE)
c  51	format (2x,'TIME= ',a9,'   PRES= ',f5.0,1x,'mb',2X,
c1	A4,10X)
	encode (51,51,timelabel(1)) TYPELABEL(ITYPE,imr),btime,gp(k)
   51	format (2x,a15,2x,'TIME= ',a9,'   PRES= ',f5.0,1x,'mb')
      else if (icord.eq.4 .or. icord.eq.7) then
        encode (51,501,timelabel(1)) TYPELABEL(ITYPE,imr),atime
  501   format (2x,a15,5x,'TIME= ',a9,14x)

      end if
      timelabel(2) =
     1       '                                                   '

 
C *** Move data from 3-d arrays into 2-d arrays
      if (icord.eq.1 .or. icord.eq.2) then
        print *,' ivar3=',ivar3
        do ivar=1,nvar3
          print *,' k,ivar,g3',k,ivar,g3(15,15,ivar,k)
        do j=1,ny_p
        do i=1,nx_p
          g(i,j,ivar) = g3(i,j,ivar,k)
          if (ivar.eq.2) g(i,j,ivar) = g3(i,j,ivar,k)/10.
        end do
        end do
        end do
      else if (icord.eq.3) then
        do ivar=1,nvar_p
          print *,' k,ivar,g3',k,ivar,g3(15,15,ivar,k)
        do j=1,ny_p
        do i=1,nx_p
          g(i,j,ivar) = g3(i,j,ivar,k)
          if (ivar.eq.2) g(i,j,ivar) = g3(i,j,ivar,k)/10.
          if (ivar.eq.3) g(i,j,ivar) = g3(i,j,ivar,k)-273.15
          if (ivar.eq.4) g(i,j,ivar) = 100.*
     1        g3(i,j,ivar,k)
        end do
        end do
        end do
      end if

      do ivar=1,nvar_p
        print *,' ivar,g',ivar,g(15,15,ivar)
      end do

887   continue

C  Special plotting routine for 2-d MSP files
      if (icord.eq.4) then
        call plotmsp (g4,spd,nx_p,ny_p,nvar_h,ic,timelabel)
      else if (icord.eq.5) then
C  Special plotting routine for 2-d surface analysis MSX files
        call plotmsx (g4,spd,nx_p,ny_p,nvar_s,ic,atime)
      else if (icord.eq.6) then
        call plotmspm (g4m,spd,nx_p,ny_p,nvar_ms,ic,timelabel)
      else if (icord.eq.7) then
        call plotmspr (g4m,spd,nx_p,ny_p,nvarsfc_ncep_p,ic,timelabel)
      else 

C  Plot for 3-d grids, either hybrid isentropic/sigma or isobaric

      IC = IC + 1

      type *,' Enter variable to be plotted:'
      if (icord.eq.1 .or. icord.eq.2) then
      type *,' 1 - pressure              2 - Montgomery stream function'
      type *,' 3 - potential temp.       4 - water vapor mixing ratio'
      type *,' 5 - winds                 6 - vertical velocity'
      type *,' 7 - potential vorticity   8 - absolute vorticity'
      type *,' 9 - adiabatic omega      10 - flight time (no plot)'
      type *,' 11- divergence           12 - stretching deformation'
      type *,' 13- shearing deformation 14 - condensation pres. deficit'
      type *,' 15- cloud wat mix ratio  16 - rain wtr mixing ratio'
      type *,' 17- ice mixing ratio     18 - snow mixing ratio'
      type *,' 19- graupel mixing ratio 20 - ice particle no. concen'
      type *,' 21- turbulent kin enrgy'
      type *,' 22- cloud top pressure '
      else if (icord.eq.3.or.icord.eq.6) then
      type *,' 1 - vertical velocity     2 - heights '
      type *,' 3 - temperature           4 - relative humidity' 
      type *,' 5 - winds  '
      type *,' 7 - potential vorticity   8 - absolute vorticity'
      type *,'                          10 - flight time (no plot)'
      type *,' 11- divergence           12 - stretching deformation'
      type *,' 13- shearing deformation'
      end if

      ACCEPT 102, IVARTY

      type *,'Do you want max/mins plotted. both-1, hi only-2'
      type *,'                              no-0,   lo only-3'
      accept 102, ihl
	    
       IF (IVARTY.lt.5) THEN
        if (icord.eq.1 .or. icord.eq.2) then
        if (ivarty.eq.3) then
	  print *,' Do you want theta-v/temp/tsfc-tsno/tsfc-t2m,1/2/3/4'
          accept *, itemp
           if (itemp.ge.2) then
             LABEL(3)(1:12) = 'TEMP (C)    '
              if (itemp.eq.3) LABEL(3)(1:12) = 'TSFC-TG (C) '
            do j = 1,ny_p
            do i = 1,nx_p
              temp = g(i,j,3)*((g(i,j,1)/1000.)**.285714)
     1           / (1.+ 0.6078*g(i,j,4)/1000.)
c             call thvpc2tq (g(i,j,3),g(i,j,4),g(i,j,1),
c    1          temp,q)
              g(i,j,3) = temp - 273.15
              if (itemp.eq.3) g(i,j,3) = temp-g4(i,j,38)
              if (itemp.eq.4) g(i,j,3) = temp-g4(i,j,43)
            end do
            end do
           end if
        end if
        end if
         
	type *,' HOW MANY TIMES SMOOTH'
	ACCEPT 102, NSM
	TYPE *,' DO YOU WANT DESMOOTHING?  YES-1, NO-0'
	ACCEPT 102, DSM

	DO ISM = 1,NSM
	  CALL SMOOTH(G(1,1,IVARTY),HOLD,NX_P,NY_P,0.5)
	  IF (DSM.EQ.1) 
	1  CALL SMOOTH(G(1,1,IVARTY),HOLD,NX_P,NY_P,-0.5)
	END DO
C
C     CALC MAX AND MIN OF F
C
      FMAX = -1.E8
      FMIN =  1.E8
      DO 55 J=1,NY_P
        DO 60 I=1,NX_P
	  FMAX = AMAX1(FMAX,G(I,J,IVARTY))
	  FMIN = AMIN1(FMIN,G(I,J,IVARTY))
60      CONTINUE
55    CONTINUE

	TYPE *,' MIN =',FMIN,' MAX =',FMAX
	TYPE *,' ENTER CONTOUR INTERVAL'
	ACCEPT 19, GINTER
19	FORMAT (F10.3)

      CALL PLOT_CONTOUR (1,NX_P,NX_P,NY_P,g(1,1,ivarty),
     1    LABel(ivarty),gINTer,gmin,gmax,
     1             g(1,1,5),g(1,1,6),IVARty,TIMELABEL,IC,1,IHL)

c     CALL PLOT_CONTOUR(1,NX_P,NX_P,NY_P,G(1,1,IVARTY),
c1	LABEL(1,IVARTY),GINTER,GMIN,GMAX,
c1		G(1,1,5),G(1,1,6),IVARTY,TIMELABEL,k,iob)

       ELSE IF (IVARTY.EQ.5) THEN
	      TYPE *,' ISOTACHS? YES-1,NO-0'
	      ACCEPT 102,ISOT
	      TYPE *,' WIND BARBS? YES-1,NO-0'
	      ACCEPT 102,IBARB
              type *,' m/s (1) or kts (2)'
              accept 102, ikts
              IF (Ikts.EQ.2) THEN
                DO J = 1,NY_P
                DO I = 1,NX_P
                  G(I,J,5) = g(i,j,5)*1.94
                  G(I,J,6) = g(i,j,6)*1.94
                END DO
                END DO
                LABEL(5)(1:12) = 'WIND MAG KTS'
              END IF

	      IF (ISOT.GT.0) THEN
	        DO J = 1,NY_P
	        DO I = 1,NX_P
	          spd(I,J)=SQRT(G(I,J,5)*G(I,J,5)+G(I,J,6)*G(I,J,6))
	        END DO
	        END DO
	      ELSE
	        DO J = 1,NY_P
	        DO I = 1,NX_P
	          spd(I,J) = 0.
	        END DO
	        END DO
	      END IF
	      IF (IBARB.EQ.0) THEN
	        DO J = 1,NY_P
	        DO I = 1,NX_P
	          G(I,J,5) = 0.
	          G(I,J,6) = 0.
	        END DO
	        END DO
	      END IF
C
C     CALC MAX AND MIN OF F
C
      FMAX = -1.E8
      FMIN =  1.E8
      DO  J=1,NY_P
        DO  I=1,NX_P
	  FMAX = AMAX1(FMAX,spd(I,J))
	  FMIN = AMIN1(FMIN,spd(I,J))
	end do
	end do

	TYPE *,' MIN =',FMIN,' MAX =',FMAX
	TYPE *,' ENTER CONTOUR INTERVAL'
	ACCEPT 19, GINTER

	print *,' before plot, g(1,1,5)',g(1,1,5)

      CALL PLOT_CONTOUR (1,NX_P,NX_P,NY_P,spd,
     1    LABel(5),gINTer,gmin,gmax,
     1             g(1,1,5),g(1,1,6),5,TIMELABEL,IC,1,IHL)

C -------------------------------
      ELSE IF (IVARTY.EQ.6) THEN
C -------------------------------
        ivart=7 
	type *,' HOW MANY TIMES SMOOTH'
	ACCEPT 102, NSM
	TYPE *,' DO YOU WANT DESMOOTHING?  YES-1, NO-0'
	ACCEPT 102, DSM

	DO ISM = 1,NSM
	  CALL SMOOTH(G(1,1,IVART),HOLD,NX_P,NY_P,0.5)
	  IF (DSM.EQ.1) 
	1  CALL SMOOTH(G(1,1,IVART),HOLD,NX_P,NY_P,-0.5)
	END DO
C
C     CALC MAX AND MIN OF F
C
      FMAX = -1.E8
      FMIN =  1.E8
      DO J=1,NY_P
        DO I=1,NX_P
          if (g(i,j,ivart).gt.fmax) then
            imax=i
            jmax=j
          end if
          if (g(i,j,ivart).lt.fmin) then
            imin=i
            jmin=j
          end if
	  FMAX = AMAX1(FMAX,G(I,J,IVART))
	  FMIN = AMIN1(FMIN,G(I,J,IVART))
        end do
      end do

	TYPE *,' MIN =',FMIN,' MAX =',FMAX
        type *,' min i/j =',imin,jmin
        type *,' max i/j =',imax,jmax
	TYPE *,' ENTER CONTOUR INTERVAL'
	ACCEPT 19, GINTER

      CALL PLOT_CONTOUR (1,NX_P,NX_P,NY_P,g(1,1,ivart),
     1    LABel(10),gINTer,gmin,gmax,
     1             g(1,1,5),g(1,1,6),IVARt,TIMELABEL,IC,1,IHL)

	ELSE IF (IVARTY.EQ.7) THEN
C       --------------------------

        do ivar=1,nvar_p
        do j = 1,ny_p
        do i = 1,nx_p
          g2(i,j,ivar) = g3(i,j,ivar,k+1)
        end do
        end do
        end do


	JJ = NY_P
	II = NX_P
	DO J = 1,JJ
	  J1 = MAX(1,J-1)
	  J2 = MIN(JJ,J+1)
	DO I = 1,II
	  I1 = MAX(1,I-1)
	  I2 = MIN(II,I+1)

	  VORT = 0.25/SCALE(I,J)* (G(I2,J,V_P)-G(I1,J,V_P)
	1	-G(I,J2,U_P)+G(I,J1,U_P))
	1         + 0.25/SCALE(I,J)* (G2(I2,J,V_P)-G2(I1,J,V_P)
	1	-G2(I,J2,U_P)+G2(I,J1,U_P))
          if (icord.eq.1 .or. icord.eq.2) then
	    PDIF = AMAX1(10.,G(I,J,P_P)-G2(I,J,P_P) )
	    STAB = PDIF/(GT(K)-GT(K+1))
          else if (icord.eq.3) then
            pdif = gp(k) - gp(k+1)
            thdif = g(i,j,t_p)*(gp(k)/1000.)**0.286 -
     1             g2(i,j,t_p)*(gp(k+1)/1000.)**0.286
            stab = pdif / max(1,thdif)
          end if

	  IF (STAB.EQ.0.) STAB = 0.01 !
	  ABSVOR(I,J) = 1.E5 * (VORT+CORIO(I,J))
	  PV(I,J) = -ABSVOR(I,J)/STAB
	  if (i.eq.26 .and. j.eq.6 )
	1	print *,' vort,corio,pv,stab',vort,corio(i,j),pv(i,j),
	1	stab, G(I,J,P_P),G2(I,J,P_P)
	IF (ABS(PV(I,J)).GT.100.)
	1	print *,' I/J,vort,corio,pv,stab',I,J,vort,corio(i,j),
	1	pv(i,j),stab, G(I,J,P_P),G2(I,J,P_P)
	END DO
	END DO

	TYPE *,' ENTER CONTOUR INTERVAL'
	ACCEPT 19, GINTER

      CALL PLOT_CONTOUR (1,NX_P,NX_P,NY_P,pv,
     1    LABelpv,gINTer,gmin,gmax,
     1             g(1,1,5),g(1,1,6),4,TIMELABEL,IC,1,IHL)

	ELSE IF (IVARTY.EQ.8) THEN
C       --------------------------

	JJ = NY_P
	II = NX_P
	DO J = 1,JJ
	  J1 = MAX(1,J-1)
	  J2 = MIN(JJ,J+1)
	DO I = 1,II
	  I1 = MAX(1,I-1)
	  I2 = MIN(II,I+1)

	  VORT = 0.5/SCALE(I,J) * (G(I2,J,V_P)-G(I1,J,V_P)
	1	-G(I,J2,U_P)+G(I,J1,U_P))

	  ABSVOR(I,J) = 1.E5 * (VORT+CORIO(I,J))
	END DO
	END DO

      FMAX = -1.E8
      FMIN =  1.E8
      DO  J=1,NY_P
        DO  I=1,NX_P
	  FMAX = AMAX1(FMAX,ABSVOR(I,J))
	  FMIN = AMIN1(FMIN,ABSVOR(I,J))
	end do
	end do

	TYPE *,' MIN =',FMIN,' MAX =',FMAX
	TYPE *,' ENTER CONTOUR INTERVAL'
	ACCEPT 19, GINTER


      CALL PLOT_CONTOUR (1,NX_P,NX_P,NY_P,absvor,
     1    LABelv,gINTer,gmin,gmax,
     1             g(1,1,5),g(1,1,6),4,TIMELABEL,IC,1,IHL)
c     CALL PLOT_CONTOUR(1,NX_P,NX_P,NY_P,ABSVOR,
c	LABELV,2.,GMIN,GMAX,
c1		G(1,1,5),G(1,1,6),4,TIMELABEL,k,iob)
	
	ELSE IF (IVARTY.EQ.9) THEN     ! ADIABATIC OMEGA
C       --------------------------
	
	TYPE *,' ENTER U/V PHASE SPEED FOR WAVE IN M/S'
	TYPE *,' FIRST ENTER U COMPONENT'
	ACCEPT 104, UWAVE
	TYPE *,' NOW ENTER V COMPONENT'
	ACCEPT 104, VWAVE
104	FORMAT (F5.1)


	DO J = 1,NY_P
	  J1 = MAX(J-1,1)
	  J2 = MIN(NY_P,J+1)
	DO I = 1,NX_P
	  I1 = MAX(I-1,1)
	  I2 = MIN(NX_P,I+1)

	  PV(I,J) = 1.E3/SCALE(I,J) *   ! USE PV ARRAY FOR OMEGA
	1	((G(I,J,U_P)-UWAVE)*(G(I2,J,P_P) - G(I1,J,P_P)) 
	1      + (G(I,J,V_P)-VWAVE)*(G(I,J2,P_P) - G(I,J1,P_P)))
	  IF (I.EQ.15.AND.J.EQ.15) THEN
	    PRINT *,' OMEGA,U,V,P1/2/3/4',PV(I,J),
	1	G(I,J,U_P),G(I,J,V_P),G(I2,J,P_P),G(I1,J,P_P),
	1	G(I,J2,P_P),G(I,J1,P_P)
	  END IF

	END DO
	END DO

C
C     CALC MAX AND MIN OF F
C
      FMAX = -1.E8
      FMIN =  1.E8
      DO  J=1,NY_P
        DO  I=1,NX_P
	  FMAX = AMAX1(FMAX,PV(I,J))
	  FMIN = AMIN1(FMIN,PV(I,J))
	end do
	end do

	TYPE *,' MIN =',FMIN,' MAX =',FMAX
	TYPE *,' ENTER CONTOUR INTERVAL'
	ACCEPT 19, GINTER

      CALL PLOT_CONTOUR (1,NX_P,NX_P,NY_P,pv,
     1    LABelomg,gINTer,gmin,gmax,
     1             g(1,1,5),g(1,1,6),3,TIMELABEL,IC,1,IHL)
c      CALL PLOT_CONTOUR(1,nx_p,nx_p,ny_p,PV,
c1	LABELOMG,GINTER,GMIN,GMAX,
c1		G(1,1,5),G(1,1,6),3,TIMELABEL,k,iob)

	ELSE IF (IVARTY.EQ.10) THEN     ! FLIGHT TIME
C       ---------------------------

	  AIRSPEED = 225.    ! 225 M/S

	  TTOT = 0.
	  TYPE *,' ENTER J POINT FOR W->E FLIGHT'
	  ACCEPT *, J
	  TYPE *,' ENTER I W/E LIMITS'
	  ACCEPT *,I1,I2
	  DO I = I1,I2       ! across the country
	    DT = (SCALE(I,J)+SCALE(I+1,J)) /
	1	 (2.*AIRSPEED + G(I,J,5) + G(I+1,J,5))
	    TTOT = TTOT + DT
	  END DO
	  TYPE *,' FLIGHT TIME = ',TTOT

      ELSE IF (IVARTY.EQ.11) THEN	! divergence
C     ---------------------------
	JJ = NY_P
	II = NX_P
	DO J = 1,JJ
	  J1 = MAX(1,J-1)
	  J2 = MIN(JJ,J+1)
	DO I = 1,II
	  I1 = MAX(1,I-1)
	  I2 = MIN(II,I+1)

	  VORT = 0.5/SCALE(I,J) * (G(I2,J,U_P)-G(I1,J,U_P)
	1	+ G(I,J2,V_P)-G(I,J1,V_P))

	  ABSVOR(I,J) = 1.E6 * VORT
	END DO
	END DO

C
C     CALC MAX AND MIN OF F
C
      FMAX = -1.E8
      FMIN =  1.E8
      DO  J=1,NY_P
        DO  I=1,NX_P
	  FMAX = AMAX1(FMAX,ABSVOR(I,J))
	  FMIN = AMIN1(FMIN,ABSVOR(I,J))
	end do
	end do

	TYPE *,' MIN =',FMIN,' MAX =',FMAX
	TYPE *,' ENTER CONTOUR INTERVAL'
	ACCEPT 19, GINTER

      CALL PLOT_CONTOUR (1,NX_P,NX_P,NY_P,absvor,
     1    LABel(6),gINTer,gmin,gmax,
     1             g(1,1,5),g(1,1,6),4,TIMELABEL,IC,1,IHL)
c      CALL PLOT_CONTOUR(1,NX_P,NX_P,NY_P,ABSVOR,
c1	LABEL(1,6),GINTER,GMIN,GMAX,
c1		G(1,1,5),G(1,1,6),4,TIMELABEL,k,iob)

      ELSE IF (IVARTY.EQ.12) THEN	! STRETCHING DEFORMATION
C     ---------------------------

	JJ = NY_P
	II = NX_P
	DO J = 1,JJ
	  J1 = MAX(1,J-1)
	  J2 = MIN(JJ,J+1)
	DO I = 1,II
	  I1 = MAX(1,I-1)
	  I2 = MIN(II,I+1)

	  VORT = 0.5/SCALE(I,J) * (G(I2,J,U_P)-G(I1,J,U_P)
	1	-G(I,J2,V_P)+G(I,J1,V_P))

	  ABSVOR(I,J) = 1.E6 * VORT
	END DO
	END DO

C
C     CALC MAX AND MIN OF F
C
      FMAX = -1.E8
      FMIN =  1.E8
      DO  J=1,NY_P
        DO  I=1,NX_P
	  FMAX = AMAX1(FMAX,ABSVOR(I,J))
	  FMIN = AMIN1(FMIN,ABSVOR(I,J))
	end do
	end do

	TYPE *,' MIN =',FMIN,' MAX =',FMAX
	TYPE *,' ENTER CONTOUR INTERVAL'
	ACCEPT 19, GINTER

      CALL PLOT_CONTOUR (1,NX_P,NX_P,NY_P,absvor,
     1    LABel(7),gINTer,gmin,gmax,
     1             g(1,1,5),g(1,1,6),4,TIMELABEL,IC,1,IHL)
c      CALL PLOT_CONTOUR(1,NX_P,NX_P,NY_P,ABSVOR,
c1	LABEL(1,7),GINTER,GMIN,GMAX,
c1		G(1,1,5),G(1,1,6),4,TIMELABEL,k,iob)

      ELSE IF (IVARTY.EQ.13) THEN	! SHEAR DEFORMATION
C     ---------------------------
	JJ = NY_P
	II = NX_P
	DO J = 1,JJ
	  J1 = MAX(1,J-1)
	  J2 = MIN(JJ,J+1)
	DO I = 1,II
	  I1 = MAX(1,I-1)
	  I2 = MIN(II,I+1)

	  VORT = 0.5/SCALE(I,J) * (G(I2,J,V_P)-G(I1,J,V_P)
	1	+ G(I,J2,U_P)-G(I,J1,U_P))

	  ABSVOR(I,J) = 1.E6 * VORT
	END DO
	END DO

C
C     CALC MAX AND MIN OF F
C
      FMAX = -1.E8
      FMIN =  1.E8
      DO  J=1,NY_P
        DO  I=1,NX_P
	  FMAX = AMAX1(FMAX,ABSVOR(I,J))
	  FMIN = AMIN1(FMIN,ABSVOR(I,J))
	end do
	end do

	TYPE *,' MIN =',FMIN,' MAX =',FMAX
	TYPE *,' ENTER CONTOUR INTERVAL'
	ACCEPT 19, GINTER

      CALL PLOT_CONTOUR (1,NX_P,NX_P,NY_P,absvor,
     1    LABel(8),gINTer,gmin,gmax,
     1             g(1,1,5),g(1,1,6),4,TIMELABEL,IC,1,IHL)
c      CALL PLOT_CONTOUR(1,NX_P,NX_P,NY_P,ABSVOR,
c1	LABEL(1,8),GINTER,GMIN,GMAX,
c1		G(1,1,5),G(1,1,6),4,TIMELABEL,k,iob)

      ELSE IF (IVARTY.EQ.14) THEN	! CONDENSATION PRESSURE DEFICIT
C     ---------------------------

	DO J = 1,NY_P
	DO I = 1,NX_P
	  ABSVOR(I,J) = G(I,J,1) - G(I,J,4)
	END DO
	END DO

C
C     CALC MAX AND MIN OF F
C
      FMAX = -1.E8
      FMIN =  1.E8
      DO  J=1,NY_P
        DO  I=1,NX_P
	  FMAX = AMAX1(FMAX,ABSVOR(I,J))
	  FMIN = AMIN1(FMIN,ABSVOR(I,J))
	end do
	end do

	TYPE *,' MIN =',FMIN,' MAX =',FMAX
	TYPE *,' ENTER CONTOUR INTERVAL'
	ACCEPT 19, GINTER

        type *,' Do you want to set max/min for plotting.'
        type *,' Enter 1 for yes'
        accept *, imxmn
        if (imxmn.eq.1) then
          type *,' Enter max contour value (real)'
          accept *,fmax
          type *,' Enter min contour value (real)'
          accept *,fmin
        end if

      CALL PLOT_CONTOUR (1,NX_P,NX_P,NY_P,absvor,
     1    LABel(9),gINTer,gmin,gmax,
     1             g(1,1,5),g(1,1,6),4,TIMELABEL,IC,1,IHL)
c      CALL PLOT_CONTOUR(1,NX_P,NX_P,NY_P,ABSVOR,
c1	LABEL(1,9),GINTER,GMIN,GMAX,
c1		G(1,1,5),G(1,1,6),4,TIMELABEL,k,iob)

        fmax = 0.0
        fmin = 1.0

       ELSE IF (IVARTY.ge.15.and.ivarty.lt.22) THEN
C     ---------------------------  cloud variables or tke
         
	type *,' HOW MANY TIMES SMOOTH'
	ACCEPT 102, NSM
	TYPE *,' DO YOU WANT DESMOOTHING?  YES-1, NO-0'
	ACCEPT 102, DSM
        ivartz = ivarty - 7
C         cloud water = variable 8 in 3d file
C                     = variable 15 in selection list

	DO ISM = 1,NSM
	  CALL SMOOTH(G(1,1,IVARTz),HOLD,NX_P,NY_P,0.5)
	  IF (DSM.EQ.1) 
	1  CALL SMOOTH(G(1,1,IVARTz),HOLD,NX_P,NY_P,-0.5)
	END DO
C
C     CALC MAX AND MIN OF F
C
      FMAX = -1.E8
      FMIN =  1.E8
      DO J=1,NY_P
        DO I=1,NX_P
	  FMAX = AMAX1(FMAX,G(I,J,IVARTz))
	  FMIN = AMIN1(FMIN,G(I,J,IVARTz))
      end do
      end do

	TYPE *,' MIN =',FMIN,' MAX =',FMAX
	TYPE *,' ENTER CONTOUR INTERVAL'
	ACCEPT 19, GINTER

      CALL PLOT_CONTOUR (1,NX_P,NX_P,NY_P,g(1,1,ivartz),
     1    LABel(ivarty-4),gINTer,gmin,gmax,
     1             g(1,1,5),g(1,1,6),IVARty,TIMELABEL,IC,1,IHL)

       ELSE IF (IVARTY.eq.22) THEN
C     ---------------------------  cloud top pressure

c     cloud_def = 1.e-5
      type *,' Enter cloud mixing ratio threshold in g/kg'
      type *,'This is for integrated mixing ratio -qc/qr/qs/qi/qg'
      accept 19,cloud_def
         
      do j=1,ny_p
      do i=1,nx_p
      do k=nztn_p,1,-1
        sumq(k) = 0.
        do ivar = 8,12
C         sum over all hydrometeor types
          sumq(k) = sumq(k) + g3(i,j,ivar,k)+1.e-5
        end do
      end do

      do k1=nztn_p-1,2,-1
        if (sumq(k1).gt.cloud_def) go to 771
      end do
771   continue
      frac = (cloud_def - sumq(k1)) / (sumq(k1+1)-sumq(k1))
      if (k1.le.2) then
        pcloud(i,j) = 1500.
      else 
        pcloud(i,j) = g3(i,j,p_p,k1) - frac*(g3(i,j,p_p,k1+1)
     1      - g3(i,j,p_p,k1))
      end if
      pcloud(i,j) = max(pcloud(i,j),150.)
      end do
      end do

C
C     CALC MAX AND MIN OF F
C
      FMAX = -1.E8
      FMIN =  1.E8
      DO J=1,NY_P
        DO I=1,NX_P
	  FMAX = AMAX1(FMAX,pcloud(I,J))
	  FMIN = AMIN1(FMIN,pcloud(I,J))
      end do
      end do

	TYPE *,' MIN =',FMIN,' MAX =',FMAX
	TYPE *,' ENTER CONTOUR INTERVAL'
	ACCEPT 19, GINTER

      CALL PLOT_CONTOUR (1,NX_P,NX_P,NY_P,pcloud,
     1    LABel(18),gINTer,gmin,gmax,
     1             g(1,1,5),g(1,1,6),IVARty,TIMELABEL,IC,1,IHL)

      end if

C --- end if for icord.eq.4 vs. other icords
      end if

c     type *,' Do you want a new frame (versus overlay)?'
c     type *,' enter 1-yes, 0-no'
c     accept 102, ixap
c     if (ixap.eq.1) then
c       ic = 0
c       call frame
c     end if

      type *,' Want another plot from the same time and level?'
      type *,' enter 1-yes, 0-no'
      accept 102, ixap
      if (ixap.eq.1) go to 887

      if (icord.lt.4.or.icord.eq.6) then
        type *,' Another plot from same time but different level?'
        type *,' enter 1-yes, 0-no'
        accept 102, ixap
        if (ixap.eq.1) go to 882
      end if

      end if   ! iskewt option

      type *,' Do you desire a new time/file?'
      type *,' enter 1-yes, 0-no'
      accept 102, ixap
      if (ixap.eq.1) go to 880

      call clsgks
  
      STOP
      END


