C conplt.for
C ----------------------------------------------------------------------!-------
! 26 Sep 2011 - Allow terrain contours for vertical plot projections
!               Needed only for testing terrain (bottom) models
C 16 Mar 2010 - Increased dimension of ..GP.. arrays to 611
C 15 Mar 2010 - Added common blocks B21,B22,B23,B24 
C 15 Mar 2010 - Added calls to PLTNOTE to label contours
C 10 Sep 2008 - Made data types consistent.
C 13 Mar 2008 - Corrected length in common block /MM/
C 13 Mar 2008 - Correction to common block /HDR/
C 13 Mar 2008 - Corrected length in common block /PCONST/
C      1 Nov 2006 - Added printing of date and time
C      1 Nov 2006 - Fixed modmtxt and modttxt problem
C     27 Mar 2006 - Fixed type for common block ddscale
C      6 Mar 2006 - Added call to setdsp
C      6 Mar 2006 - Moved call to idisper
C      3 Mar 2006 - Changed units of Brunt-Vaisala frequency to Hz
C      3 Mar 2006 - Changed units of acoustic cutoff frequency to Hz
C      3 Mar 2006 - Corrections
C      3 Mar 2006 - Added square of Brunt-Vaisala frequency 
C      3 Mar 2006 - Added Brunt-Vaisala frequency 
C      3 Mar 2006 - Added acoustic cutoff frequency 
C      2 Mar 2006 - Added an argument to READW  and SETOVR
C ----------------------------------------------------------------------!-------
Contains CONPLT
C ----------------------------------------------------------------------!-------
C     28 June 2005 - Modified "Models Used" column
C     28 June 2005 - Changed km/s printing to K for temperature contours
C     23 June 2005 - Added temperature contour plotting with labeling
c     14 June 2005 - Added W29 full functionality
c     2 May 2005 - changed the arguments to the first call to CONHEAD,
c     so that arguments 6 and 7 are double precision
C ----------------------------------------------------------------------!-------
      PROGRAM CONPLT
C     PLOTS CONTOURS OF TEMPERATURE, SOUND SPEED, WIND/CURRENT VELOCITY, ETC. 
C     MAIN PROGRAM FOR THE RAY TRACING PACKAGE TO BE RUN AS A
C     PERIPHERAL PROGRAM.
C
C          CONTOURS OF SPEED FUNCTION IN ANY VERTICAL PLANE
C.....This program calculates contours of 1 of 6 different parameters. The
C     parameter to be contoured is indicated by the value of W(90) = PLTC.
C     The possibilities are:
C
C     1 = Sound Speed
C     2 = Current/Wind Speed
C     3,4,5 = Current/Wind Speed Component (Vertical,South,East)
C     6 = Terrain Height
C     7 = Temperature
C     8 = Wind direction
C     9 = Acoustic cutoff frequency
C    10 = Brunt-Vaisala frequency
C    11 = Square of Brunt-Vaisala frequency
C
C     These contours will be made in 1 of 4 different planes/projections which
C     is indicated by the value of W(81) = PROJ.  The possibilities are:
C
C     1 = Vertical Plane, Polar Plot, Rectangular Expansion
C     2 = Horizontal Plane, lateral expansion
C     3 = Vertical Plane, Polar Plot, Radial Expansion
C     4 = Vertical Plane, Rectangular Plot
C
C     NOTE:  Projections 1,3,4 are essentially treated the same throughout
C            the program, but projection 2 of course requires completely
C            different processing for the frame and labeling.
C
C.....HISTORY
C
C     This program was written over the course of almost 20 years.  Therefore,
C     much of it is unstructured and difficult to follow.  Because the
C     massive task of restructuring it would be too time-consuming and because
C     restructuring this program would make it incompatible with the HARPO pro-
C     gram, most of the program has been left untouched.
C
C     The major changes to this program were made to make the program
C     PC-compatible and to make its output compatible with the PSGRAPH program.
C     The PSGRAPH-compatibility changes are located in the main routine, the
C     PLOTXY subroutine and in routines which interact with it.
C
C     In a few places, GOTO's have been replaced with IF-THEN logic.  Also,
C     the error handling for the W-array in the main routine has been changed.
C     Previously, there was usually no message to the user that an invalid
C     value was causing the program to skip the runset entirely.
C     Indentation has also been added.
C     Most of the new changes are in lower case.

      implicit double precision (a-h,o-z)

      TYPE :: value_and_derivs                                   ! 3 Mar 2006
         REAL(KIND=8) :: value          
         REAL(KIND=8) :: pt,pr,pth,pph
      END TYPE 

      TYPE (value_and_derivs) :: kA2                             ! 3 Mar 2006
      TYPE (value_and_derivs) :: Ns                              ! 3 Mar 2006

      CHARACTER CTEMP,U*4,LABEL*80
      DIMENSION AM(46,46),LONG(46),LATT(46),CC(150)

C.....FYI:  In the following COMMONs, the variables beginning with "P"
C           are usually a partial derivative e.g., PTR is the partial of
C           temperature with respect to radius and PTTH is the partial of
C           temperature with respect to theta (latitude).

C     COMMON DECK "FLAG" INSERTED HERE
      LOGICAL NEWWR,NEWWP,NEWTRC,PENET
      COMMON /FLG/ NTYP,NEWWR,NEWWP,NEWTRC,PENET,LINES,IHOP,HPUNCH
      COMMON/FLGP/NSET

C     COMMON DECK "TT" INSERTED HERE:  Temperature Model Common
      double precision MODT
      COMMON/TT/MODT(4), T,PTT,PTR,PTTH,PTPH

C     COMMON DECK "MM" INSERTED HERE:  Molecular Weight Model Common
      double precision M,MODM
      COMMON/MM/MODM(4),M,PMT,PMR,PMTH,PMPH
     1 ,XO2,XN2,XO,XN,XH2O,XO3                        ! 5 Mar 2007 ! 13 Mar 2008

C     COMMON DECK "UU" INSERTED HERE:  Ocean Current Model Common
      double precision MODU
      COMMON/UU/MODU(4)
     1  ,V  ,PVT  ,PVR  ,PVTH  ,PVPH
     2  ,VR ,PVRT ,PVRR ,PVRTH ,PVRPH
     3  ,VTH,PVTHT,PVTHR,PVTHTH,PVTHPH
     4  ,VPH,PVPHT,PVPHR,PVPHTH,PVPHPH

C     COMMON DECK "HDR" INSERTED HERE
      CHARACTER*10 INITID*80,DAT,TOD
      COMMON/HDR/SEC,stsec                            !Rev 7/13/90 ! 13 Mar 2008
      COMMON/HDRC/INITID,DAT,TOD

C     COMMON DECK "GG" INSERTED HERE: Terrain Model Common
      double precision MODG
      COMMON/GG/MODG(4)
      COMMON/GG/G,PGR,PGRR,PGRTH,PGRPH
      COMMON/GG/PGTH,PGPH,PGTHTH,PGPHPH,PGTHPH,GSELECT,GTIME

C     COMMON DECK "CONST" INSERTED HERE
      COMMON/PCONST/CREF,RGAS,GAMMA
     1 ,MO2,MN2,MO,MN,MH2O,MO3                        ! 9 Mar 2007 ! 13 Mar 2008
      DOUBLE PRECISION MO2,MN2,MO,MN,MH2O,MO3                      ! 13 Mar 2008
      COMMON/MCONST/PI,PIT2,PID2,DEGS,RAD,ALN10

C     COMMON DECK "CC" INSERTED HERE: Sound Speed Model Common
      double precision MODC
      COMMON/CC/MODC(4),CS,PCST,PCSR,PCSTH,PCSPH

C     COMMON DECK "RKAM" INSERTED HERE
      double precision KR,KTH,KPH
      COMMON//R,TH,PH,KR,KTH,KPH,RKVARS(14),TPULSE,CSTEP,DRDT(20)

C     COMMON DECK "WWR" INSERTED HERE
      character*80 id
      common /ident/ id
      PARAMETER (NWARSZ=1000)
      COMMON/WW/ dum(10),MAXW,W(NWARSZ)
      double precision LLAT,LLON
      EQUIVALENCE (EARTHR,W(1)), (RAYFNC,W(29)), (PROJ,W(81)),
     &            (PFACTR,W(82)), (LLAT,W(83)), (LLON,W(84)),
     &            (RLAT,W(85)), (RLON,W(86)), (TIC,W(87)),
     &            (HB,W(88)), (HT,W(89)), (TICV,W(96)),
     &            (PLTC,W(90)), (CBEG,W(92)), (CEND,W(93)),
     &            (CNSTEP,W(94))
C
      character*8 modmtxt(4)
      equivalence (modm(1),modmtxt(1))
      character*8 modutxt(4)
      equivalence (modu(1),modutxt(1))
      character*8 modttxt(4)
      equivalence (modt(1),modttxt(1))

C                             
C.....TXTSPEED is filled with 'WIND' or 'CURRENT' depending on whether
C     the top of the vertical plot W(89) is above or below sea level.

      character*7 txtspeed

C.....STRIM is an integer function that returns the length of a string
c     minus the trailing blanks.

      INTEGER STRIM

c     real xmin,xmax,ymin,ymax,scx,scy
      REAL(KIND=8) :: xmin,xmax,ymin,ymax,scx,scy            ! 27 Mar 2006

      COMMON /DDSCALE/ XMIN,XMAX,YMIN,YMAX,MINX,MAXX,MINY,MAXY,SCX,SCY,
     1 NSCX,NXCY,MSCX,MSCY,ISCX,ISCY

      COMMON/DD/IN,IOR,IT,IS,IC,ICC,IX,IY


      double precision LONG,LATT

C     COMMON DECK "RAYDEV" INSERTED HERE
C     DEVICE ASSIGNED TO RAYTRC INPUT FILE
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVTXT,NDEVBIN
C
      COMMON/RAYCON/MCONP

C.....Flag for deciding whether the "moveto" for the terrain profile
C     has been made yet.

C.....OCEAN is set true when top of vertical slice is at or below sea level.

      LOGICAL movedto,
     &        ocean
C
      EXTERNAL PLOTXY,PLOTLB,PLOTHB
      external readw
C
C     COMMON DECK "CB21" INSERTED HERE                            ! 15 Mar 2010
      INTEGER MX21,NTBL21,ITBL21,FRMTB21,IDS21(20) 
      COMMON/B21/MX21,NTBL21(10),ITBL21(10),FRMTB21(10),GP21(611) 
      EQUIVALENCE (GP21,IDS21),(AN21,GP21(11)) 
      EQUIVALENCE (CHAR21,ids21)
      character*80 char21

C     COMMON DECK "CB22" INSERTED HERE                            ! 15 Mar 2010
      INTEGER MX22,NTBL22,ITBL22,FRMTB22,IDS22(20) 
      COMMON/B22/MX22,NTBL22(10),ITBL22(10),FRMTB22(10),GP22(611) 
      EQUIVALENCE (GP22,IDS22),(AN22,GP22(11)) 

C     COMMON DECK "CB23" INSERTED HERE                            ! 15 Mar 2010
      INTEGER MX23,NTBL23,ITBL23,FRMTB23,IDS23(20) 
      COMMON/B23/MX23,NTBL23(10),ITBL23(10),FRMTB23(10),GP23(611) 
      EQUIVALENCE (GP23,IDS23),(AN23,GP23(11)) 

C     COMMON DECK "CB24" INSERTED HERE                            ! 15 Mar 2010
      INTEGER MX24,NTBL24,ITBL24,FRMTB24,IDS24(20) 
      COMMON/B24/MX24,NTBL24(10),ITBL24(10),FRMTB24(10),GP24(611) 
      EQUIVALENCE (GP24,IDS24),(AN24,GP24(11)) 

      INTEGER PX21    ,QTBL21(10),LTBL21(10),IRMTBL21(10)       
      INTEGER PX22    ,QTBL22(10),LTBL22(10),IRMTBL22(10)       
      INTEGER PX23    ,QTBL23(10),LTBL23(10),IRMTBL23(10)       
      INTEGER PX24    ,QTBL24(10),LTBL24(10),IRMTBL24(10)       

      DATA PX21/2/                                           
      DATA QTBL21/1,11,612,7*0/    
      DATA LTBL21/1,100,8*0/       
      DATA IRMTBL21/1,2,8*0/     
      DATA AQ21/0.d0/

      DATA PX22/2/                                           
      DATA QTBL22/1,11,612,7*0/    
      DATA LTBL22/1,100,8*0/       
      DATA IRMTBL22/1,2,8*0/     
      DATA AQ22/0.d0/

      DATA PX23/2/                                           
      DATA QTBL23/1,11,612,7*0/    
      DATA LTBL23/1,100,8*0/       
      DATA IRMTBL23/1,2,8*0/     
      DATA AQ23/0.d0/

      DATA PX24/2/                                           
      DATA QTBL24/1,11,612,7*0/    
      DATA LTBL24/1,100,8*0/       
      DATA IRMTBL24/1,2,8*0/     
      DATA AQ24/0.d0/

      DATA CTEMP/' '/
      DATA NLON,NLAT,NTRRES,IOFF/46,46,10,0/
C
      AN21=AQ21       
      MX21=PX21                 
      CALL IMOVE(NTBL21,QTBL21,10)     
      CALL IMOVE(ITBL21,LTBL21,10)     
      CALL IMOVE(FRMTB21,IRMTBL21,10) 
                              
      AN22=AQ22       
      MX22=PX22                 
      CALL IMOVE(NTBL22,QTBL22,10)     
      CALL IMOVE(ITBL22,LTBL22,10)     
      CALL IMOVE(FRMTB22,IRMTBL22,10) 
                              
      AN23=AQ23       
      MX23=PX23                 
      CALL IMOVE(NTBL23,QTBL23,10)     
      CALL IMOVE(ITBL23,LTBL23,10)     
      CALL IMOVE(FRMTB23,IRMTBL23,10) 
                              
      AN24=AQ24       
      MX24=PX24                 
      CALL IMOVE(NTBL24,QTBL24,10)     
      CALL IMOVE(ITBL24,LTBL24,10)     
      CALL IMOVE(FRMTB24,IRMTBL24,10) 

      CALL SYSDAT(DAT)       !  1 Nov 2006
      CALL SYSTIM(TOD)       !  1 Nov 2006

      CALL CONBLK

C.....The following 6 calls replace the call to SETDSP.  SETDSP is an entry
C     point in the dispersion relation routines.  At that entry point, the
C     following calls are made.  Note, however, that SETWND is for use with
C     current models and SETABS is for use with absorption models.  Since
C     these calls simply move some initial values into common blocks, there
C     is no harm in calling them here even when no current or absorption
C     model is being used.  Thus, the program does not have to perform
C     differently (at least at this early point) according to which models
C     were chosen in the DINP file.
C     The FINDMODL program (which is used for the PC version) always links
C     with one current model and one absorption model.

      CALL SETDSP                       ! 6 Mar 2006
      CALL SETSPD
      CALL SETWND
      CALL SETRCV
      CALL SETTOP
      CALL SETSUR
      CALL SETABS

      CALL SETSRD
      CALL stdint
      NEWWR=.TRUE.
C
      MCONP=1
      AMNCON=1.E30
      AMXCON=-AMNCON
      CALL SETUP
      NRUNSET=0
1     NRUNSET=NRUNSET+1

C     READ W ARRAY
c      print *, 'before readw'               ! 6 Mar 2006
c      pause       ! 6 Mar 2006
      if (READW(xx) .ne. 0.d0) then         ! 2 Mar 2005 - added argument xx
c      print *, 'after readw'               ! 6 Mar 2006
c      pause       ! 6 Mar 2006
           CALL DDFR
           CALL DDEND
           print*,'CONPLT Completed Normally.'
           stop   'Output file:  TXTOUT'
      endif
c      print *, 'after readw and endif'               ! 6 Mar 2006
c      pause       ! 6 Mar 2006
c      CALL IDISPER                                            ! 28 Mar 2006

      ICODE=ND2B(IDINT(RAYFNC))   ! 6/14/05

C PROCESS CONTOUR CALCULATIONS ONLY IF W(29) IS BEING USED
      IF((RAYFNC.NE.0.d0).AND.(IAND(ICODE,2).EQ.0)) then      ! 6/14/05
C      IF (RAYFNC .eq. 0.0) then
           print 77, nrunset
77         FORMAT(' WARNING:  NO CONTOURS GENERATED FOR RUNSET #',i2)
           print*,'Reason:   W(29) is zero.'
           print*,' '
           GO TO 1   ! rayfnc = W(29)
      endif
C
      if (setovr(xx) .ne. 1.0) then         ! 2 Mar 2005 - added argument xx
          STOP 'ERROR CALLING SETOVR'
      endif
C
C
      proj=ABS(proj)
      PLTCNT=ABS(PLTC)

! 26 Sep 2011      if (pltcnt .eq. 6 .and. proj .ne. 2) then
! 26 Sep 2011           PRINT 77, nrunset
! 26 Sep 2011           print*,'Reason:   For terrain contours [ W(90) = 6 ],'
! 26 Sep 2011           print*,'          plot projection [ W(81) ] must be 2.'
! 26 Sep 2011           print*,' '
! 26 Sep 2011           GO TO 1
! 26 Sep 2011      endif
      print 22 ,NINT(proj)
22    format (' Plot projection = ',I1)

      IF(pltcnt .gt. 11.) then                            !  3 Mar 2006
           PRINT 77, nrunset
           print*,'Reason:   W(90) has an invalid value.'
           print*,'          Valid values are 1 to 11.'         !  3 Mar 2006
           print*,' '
           GO TO 1
      endif

C     CUMULATIVE TITLE OFFSET IS BASED ON NEGATIVE VALUE OF PLOT TYPE

      IF(IOFF.NE.0.AND.PLTC.GT.0.0) CALL DDFR
      IF(PLTC.GT.0.0) IOFF=0
C
      IF(CNSTEP.LE.0.0) then
           PRINT 77, nrunset
           print*,'Reason:   Horizontal interval [ W(94) ] has an'
     &           ,'          invalid value.'
           print*,' '
           GO TO 1
      endif

      NCNLVS=1.5+(CEND-CBEG)/CNSTEP
      IF(NCNLVS.gt.150) then
           PRINT 77, nrunset
           PRINT *,'Reason:  Number of contour levels = ',NCNLVS
     &            ,'         exceeds current limit of 150.'
           print*,'Limit is a constant in main routine.  See variable'
     &           ,' NCNLVS.'
           print*,' '
           CNSTEP=0.
           GO TO 1
      endif

      CALL ITOPOG
      IF(PLTCNT.EQ.1.0.OR.PLTCNT.EQ.7.0) CALL ISPEED
      IF(PLTCNT.GE.2.0.AND.PLTCNT.LE.8.0.AND.PLTCNT.NE.6.0) CALL IWINDR ! 3Mar06
      IF(PLTCNT.EQ.7.0) CALL ITEMP                              ! 23 June 05
      IF(PLTCNT.EQ.1.0.OR.PLTCNT.EQ.7.0) CALL IMOLWT            ! 28 June 05
      IF(PLTCNT.GE.9.0.AND.PLTCNT.LE.11.0) CALL IDISPER         ! 3 Mar 2006
      THMIN=0.0
      THMAX=GCIRCLE(LLAT,LLON,RLAT,RLON)
C     PREVENT UNECESSARY ABORTS
      IF(THMAX.LE.0.0) then
           PRINT 77, nrunset
           print*,'Reason:   Invalid latitude and/or longitude values'
     &           ,'          in input data file.'
           print*,' '
           GO TO 1
      endif
C
      PRINT 11,nrunset
11    format (' Processing runset #',I2)

      IF(proj.EQ.2.) then
           YMAX=THMAX/(2.0*PFACTR)
           YMIN=-YMAX
           FX=EARTHR
           FY=EARTHR
      else
           YMIN=EARTHR+HB
           YMAX=EARTHR+HT
           FX=EARTHR
           FY=1.
      endif
C
C     DEFINE NUMBER OF STEPS NEEDED FOR TERRAIN PROFILE
      NGRTCR=(NLON-1)*NTRRES+1
      XMIN=THMIN
      XMAX=THMAX
      DPHI=(THMAX-THMIN)/(NGRTCR-1)
      RMIN=YMIN
      RMAX=YMAX
      CALL SETDDST
C
C
      IF(proj.eq.2.0) then
           MINX=165   !378
           MINY=180 !140   !180
           MAXX=888   !1023
           MAXY=903 !863   !825
      else if (proj.eq.4) then
           MINX=165   !250
           MINY=120   !100   !80    !250
           MAXX=978   !1000
           MAXY=1290  !1300  !1250  !1000  !893   !1000
      else  ! projections 1 and 3 are square.  Both sides length=813.
           MINX=165   !250
           MINY=140   !120   !100   !80    !250
           MAXX=978   !1000
           MAXY=953
      endif
C
C     DO NOT REMOVE THIS CALL SINCE IT IS NEEDED TO RESET CHARACTER SET
      ymintmp = ymin
      ymaxtmp = ymax
ccccc      CALL PLOTXY(THMIN,YMINtmp,1)

C.....DDINIT will write the ID containing user identification to the
C     output file, TXTOUT.  The graphing program, PSGRAPH, which reads
C     TXTOUT will also use this as an indication to translate the plot
C     away from the absolute lower left corner for PostScript files.

      call ddinit (8,id)

C.....The following indicates to PSGRAPH that the TXTOUT file contains
C     rayplot-like plot commands, i.e. command "1".

      write (ndevtxt,111) 1,9999,9999
      write (ndevtxt,222) 0
111   format (i2,2i5)
222   format (i5)

      CALL SETXY(proj,THMIN,RMIN,THMAX,RMAX)
C
CCCCC      IF(PLTC.LE.0.0) GO TO 70
C
ccccc WRITE(3,1200) ID
ccccc 1200 FORMAT('1',A//'    SPEED, TERRAIN HEIGHT IN METERS,   '
ccccc     1    ,'GRND RNG(KM)     LATITUDE AND LONGITUDE IN DEGS')
C
C          CREATE ARRAY TO BE CONTOURED
70    DH=(RMAX-RMIN)/(NLAT-1)
C
C     INITIALIZE GEOGRAPHICAL ROTATION PLACING CONTOUR NET ON A
C     GLOBAL GRID HAVING ITS EQUATOR THRU HORIZONTAL CENTER OF PLOT REG
C
      CALL IZPZ(PID2-LLAT,LLON,PID2-RLAT,RLON)
C
      print*,'Computing maximum and minimum contours.'
      print*,'This may take a few minutes...'
      print*,' '

      movedto = .false.

      DO 100 I=1,NGRTCR
           NSET=(I-1)/NTRRES+1
           SCX=THMIN+(I-1)*DPHI
           PHP=SCX

           IF(proj.EQ.2.) GO TO 20

           GRRANG=EARTHR*PHP
           R=EARTHR
C
C     TRANSFORM POINT ON LOCAL EQUATOR TO GEOGRAPHICAL GRID
C
           CALL ZPZ(PID2,PHP,TH,PH)
C
           CALL TOPOG
           R=R-G/PGR
C
C     PLOT TERRAIN PROFILE
C     Don't plot the terrain profile if it does not extend above the
C     bottom of the plot or if the bottom model is the horizontal model.

           scxtmp = scx

C..........For the first terrain point, perform a "moveto".

           if (.not. movedto .and. r .ge. ymin .and. w(300) .ne. 1) then
                call PLOTXY (SCXtmp,R,1)
                movedto = .true.
           endif

C..........For subsequent points, perform a "lineto".

           IF (movedto .and. r .ge. ymin .and. w(300) .ne. 1)
     &          CALL PLOTXY(SCXtmp,R,0)
C
C..........Calculate the value of the function to be contoured at each point.

      SELECT CASE (INT(PLTCNT))                                   ! 3 Mar 2006

      CASE (1) ! Sound speed
           CALL SPEED
           XX=SQRT(CS)

      CASE (2) ! Wind speed
           CALL WINDR
           XX=SQRT(VR*VR+VTH*VTH+VPH*VPH)

      CASE (3) ! Vertical wind component 
           CALL WINDR
           XX=VR

      CASE (4) ! Southerly wind component
           CALL WINDR
           XX=VTH

      CASE (5) ! Easterly wind component
           CALL WINDR
           XX=VPH

      CASE (6) ! Terrain 
           CALL TOPOG                                              ! 26 Sep 2011
           XX=G                                                    ! 26 Sep 2011

      CASE (7) ! Temperature 23 JUNE 05
           CALL TEMP                
           XX=T                 

      CASE (8) ! Wind direction 3 Mar 2006
           CALL WINDR                      
           XX=DEGS*ATAN2(-VPH,VTH)     
           IF(XX.LT.0.D0) XX = XX + 360.0D0 

      CASE (9) ! Acoustic cutoff frequency 3 Mar 2006
           CALL DISPER
           CALL speed
           CALL cutoff(kA2)
           XX=DSQRT(Cs*kA2%value)/pit2                              ! 3 Mar 2006

      CASE (10) ! Brunt-Vaisala frequency 3 Mar 2006
           CALL DISPER 
           CALL Brunt(Ns) 
           IF(Ns%value.GE.0.d0) XX=DSQRT(Ns%value)/pit2             ! 3 Mar 2006
           IF(Ns%value.LT.0.d0) XX=-DSQRT(-Ns%value)/pit2           ! 3 Mar 2006

      CASE (11) ! Square of Brunt-Vaisala frequency  3 Mar 2006
           CALL DISPER 
           CALL Brunt(Ns) 
           XX=Ns%value/(pit2*pit2)                                  ! 3 Mar 2006

      END SELECT


ccccc      WRITE(3,2200) XX,(R-EARTHR)*1000.,GRRANG
ccccc     1    ,EARTHR*(PID2-TH),EARTHR*PH,FLOAT(I),FLOAT(NSET)
C
C     USE ONLY LOW RESOLUTION POINTS FOR CONTOURING

20         IF((NSET-1)*NTRRES.NE.(I-1)) GO TO 100

C     PUT CURRENT LONGITUDE ANGLE IN U-ARRAY
           LONG(NSET)=PHP
C
           R=EARTHR
           IF(PLTCNT.NE.6.) R=EARTHR+HB

           DO 4 J=1,NLAT

C     PUT CURRENT RANGE(OR LATITUDE) VALUE IN AM-ARRAY

                IF(I.EQ.1) LATT(J)=RMIN+(J-1)*DH
                IF(proj.EQ.2.) GO TO 120
                R=LATT(J)
                GO TO 150
120             CALL ZPZ(PID2-LATT(J),PHP,TH,PH)
150             CALL TOPOG

      SELECT CASE (INT(PLTCNT))                                     ! 3 Mar 2006

      CASE (1) ! Sound speed
                CALL SPEED
                XX=SQRT(CS)

      CASE (2) ! Wind speed
                CALL WINDR
                XX=SQRT(VR*VR+VTH*VTH+VPH*VPH)

      CASE (3) ! Vertical wind component 
                CALL WINDR
                XX=VR

      CASE (4) ! Southerly wind component
                CALL WINDR
                XX=VTH

      CASE (5) ! Easterly wind component
                CALL WINDR
                XX=VPH

      CASE (6) ! Terrain 
           IF(proj.EQ.2.) then                                     ! 26 Sep 2011
                XX=-G/PGR
           else                                                    ! 26 Sep 2011
                XX=G                                               ! 26 Sep 2011
           endif                                                   ! 26 Sep 2011

      CASE (7) ! Temperature 23 JUNE 05
                CALL TEMP               
c                print*,'temperature = ',t
c                print*,'h temp = ',r
                XX=T

      CASE (8) ! Wind direction 3 Mar 2006
                CALL WINDR                      
                XX=DEGS*ATAN2(-VPH,VTH)     
                IF(XX.LT.0.D0) XX = XX + 360.0D0 

      CASE (9) ! Acoustic cutoff frequency 3 Mar 2006
                CALL DISPER
                CALL speed
                CALL cutoff(kA2)
c                print*,'Cs conplt = ',Cs
c                print*,'kA2 conplt = ',kA2%value
c                print*,'h conplt = ',r
                XX=DSQRT(Cs*kA2%value)/pit2                         ! 3 Mar 2006

      CASE (10) ! Brunt-Vaisala frequency 3 Mar 2006
                CALL DISPER 
                CALL Brunt(Ns) 
                IF(Ns%value.GE.0.d0) XX=DSQRT(Ns%value)/pit2        ! 3 Mar 2006
                IF(Ns%value.LT.0.d0) XX=-DSQRT(-Ns%value)/pit2      ! 3 Mar 2006

      CASE (11) ! Square of Brunt-Vaisala frequency  3 Mar 2006
                CALL DISPER 
                CALL Brunt(Ns) 
                XX=Ns%value/(pit2*pit2)                             ! 3 Mar 2006

      END SELECT

                AMNCON=DMIN1(AMNCON,XX)                            ! 10 Sep 2008
                AMXCON=DMAX1(AMXCON,XX)                            ! 10 Sep 2008
4               AM(NSET,J)=XX
100   CONTINUE
C
ccccc      PRINT 77,AMXCON
ccccc77    format (' Maximum Function Value = ',F14.8)
ccccc      PRINT 88,AMNCON
ccccc88    format (' Minimum Function Value = ',F14.8)
C
      AMNCON=DMAX1(AMNCON,CBEG)                                    ! 10 Sep 2008
      AMXCON=DMIN1(AMXCON,CEND)                                    ! 10 Sep 2008
      RND=(AMNCON-CBEG)/CNSTEP
      AMNCON=CBEG+AINT(RND)*CNSTEP
      IF(RND.GT.AINT(RND)) AMNCON=AMNCON+CNSTEP
      AMXCON=CBEG+AINT((AMXCON-CBEG)/CNSTEP)*CNSTEP
C
      U='km/s'
      IF(PLTCNT.EQ.6.) U='km'
      IF(PLTCNT.EQ.7.) U='K'                    ! 28 June 2005
      WRITE(LABEL,140) AMNCON,AMXCON,CNSTEP,U
140   FORMAT('Min =',F13.6,', Max =',F13.6,', Step =',F11.6,' ',A)
C
      PRINT *,'For Plotting Region and Contour Intervals Specified:'
      PRINT *,LABEL
C
C.....MODMTXT and MODTTXT are filled in by ITEMP.  If they are
C     blank, then the default "stub" routines were used e.g. NOTEMP.
      if (modttxt(1) .lt. 'A') modttxt(1) = 'NOTEMP  '
      if (modttxt(3) .lt. 'A') modttxt(3) = 'NPTEMP  '

      if (w(89) .le. 0.0) then
           txtspeed = 'CURRENT'
           ocean = .true.
      else
           txtspeed = 'WIND'
           ocean = .false.
      endif
      lenspeed = STRIM(txtspeed)

c      IF(pltc.EQ.1.0) CALL CONHEAD('SOUND SPEED CONTOURS',CTEMP,IOFF,
c     &                               4,MODC,MODTtxt,MODMtxt,MODG,
c     &                               label, STRIM(label), ocean)
c 2 May 2005
      IF(pltc.EQ.1.0) CALL CONHEAD('SOUND SPEED CONTOURS',CTEMP,IOFF,
     &                               4,MODC,MODT,MODM,MODG,
     &                               label, STRIM(label), ocean)

      IF(pltc.EQ.2.0) CALL CONHEAD (
     &                         txtspeed(1:lenspeed)//' SPEED CONTOURS',
     &                         CTEMP,IOFF,2,MODU,MODG, DUM,DUM,label,
     &                         STRIM(label), ocean)
      IF(pltc.EQ.3.0) CALL CONHEAD (
     1   txtspeed(1:lenspeed)//' Z COMPONENT CONTOURS',CTEMP,IOFF,2,
     &   MODU,MODG,DUM,DUM,label, STRIM(label), ocean)
      IF(pltc.EQ.4.0) CALL CONHEAD (
     1   txtspeed(1:lenspeed)//' S. COMPONENT CONTOURS',CTEMP,IOFF,2,
     &   MODU,MODG,DUM,DUM,label, STRIM(label), ocean)
      IF(pltc.EQ.5.0) CALL CONHEAD (
     1   txtspeed(1:lenspeed)//' E. COMPONENT CONTOURS',CTEMP,IOFF,2,
     &   MODU,MODG,DUM,DUM,label, STRIM(label), ocean)
      IF(pltc.EQ.6.0) CALL CONHEAD (
     1   'TERRAIN CONTOURS',CTEMP,IOFF,1,MODG,DUM,DUM,DUM,label,
     &   STRIM(label), ocean)
      IF(pltc.EQ.7.0) CALL CONHEAD('TEMPERATURE CONTOURS',CTEMP,IOFF,
     &                               4,MODC,MODT,MODM,MODU,
     &                               label, STRIM(label), ocean)
      IF(pltc.EQ.8.0) CALL CONHEAD (                                ! 3 Mar 2006
     1   txtspeed(1:lenspeed)//' DIRECTION CONTOURS',CTEMP,IOFF,2,
     &   MODU,MODG,DUM,DUM,label, STRIM(label), ocean)
      IF(pltc.EQ.9.0) CALL CONHEAD('ACOUSTIC CUTOFF FREQ. CONTOURS',! 3 Mar 2006
     & CTEMP,IOFF,
     &                               4,MODC,MODT,MODM,MODU,
     &                               label, STRIM(label), ocean)
      IF(pltc.EQ.10.0) CALL CONHEAD('BRUNT-VAISALA FREQ. CONTOURS', ! 3 Mar 2006
     & CTEMP,IOFF,
     &                               4,MODC,MODT,MODM,MODU,
     &                               label, STRIM(label), ocean)
      IF(pltc.EQ.11.0) CALL CONHEAD(                                ! 3 Mar 2006
     & 'SQUARE OF BRUNT-VAISALA FREQ. CONTOURS',CTEMP,IOFF,
     &                               4,MODC,MODT,MODM,MODU,
     &                               label, STRIM(label), ocean)
      IOFF=IY
C
      IF(PLTC.lt.0.) GO TO 160
C
      IF(proj.EQ.2.) then

C..........DRAWTKS is for horizontal plot projections only: W(81) = 2.

           CALL DRAWTKS (TIC,TICV/EARTHR,THMIN,THMAX,YMINtmp,YMAXtmp,
     &                   PLOTXY)

      else

C..........Vertical Plane Projections.

           TICY=.01*(RMAX-RMIN)

C..........First call to ARCTIC does the bottom X-axis.

           CALL ARCTIC(THMIN,THMAX,RMIN,TICY,PLOTHB)

C..........Second call to ARCTIC does the top X-axis.

           CALL ARCTIC(THMIN,THMAX,RMAX,-TICY,PLOTHB)
           TICKX=0.01*(THMAX-THMIN)

C..........First call to TIKLINE does the left Y-axis.

           CALL TIKLINE(THMIN,RMIN,THMIN,RMAX,TICV,TICKX,PLOTLB)

C..........Second call does the right Y-axis.

           CALL TIKLINE(THMAX,RMIN,THMAX,RMAX,TICV,-TICKX,PLOTXY)
      endif
C
      CALL PLTANOT(ID,0.d0,XFQMDL,YFQMDL,THMIN,RMIN,THMAX,RMAX,DEGS,
C    1  LLAT,LLON,RLAT,RLON,HB,HT,dABS(PLT),TIC*EARTHR,TICV,PLOTXY)
     1  LLAT,LLON,RLAT,RLON,HB,HT,dABS(proj),tic,earthr,PLOTXY)
C
C          WRITE ARRAY
C
160   continue  !!!!!DO 35 IBEG=1,NLON,7
ccccc           IEND=MIN0(NLON,IBEG+6)
ccccc           WRITE(3,2000) ID,PLT,THMIN*FX,RMIN*FY,THMAX*FX,RMAX*FY
ccccc     1       ,(LONG(I)*FX,I=IBEG,IEND)
ccccc 2000      FORMAT ('1PROGRAM CONPLT',10X,A/
ccccc     1   '  PLT=',F3.0,' XMIN=',G13.6,',YMIN=',G13.6,
ccccc     2   ',XMAX=',G13.6,',YMAX=',G13.6/(16X,7G16.6))
ccccc      DO 35 K=1,NLAT
ccccc   35      WRITE(3,2200) LATT(K)*FY,(AM(I,K),I=IBEG,IEND)
ccccc 2200      FORMAT (8G16.6)
C
      I=1
  40  CC(I)=CBEG+(I-1)*CNSTEP
      IF (CC(I).GT.CEND-.5*CNSTEP) GO TO 45
      I=I+1
      GO TO 40
  45  CONTINUE

C     PLOT THE CONTOURS

      print*,'number of contour intervals = ',i
      numpts = 0

      CALL DDCON  (LONG,LATT,AM,CC,46,46,46,I)
C
      print*,' '

C.....Process the next runset...

c      do kk = 1,4                           ! 1 Nov 2006
c           modmtxt(kk) = ' '                ! 1 Nov 2006
c           modttxt(kk) = ' '                ! 1 Nov 2006
c      enddo                                 ! 1 Nov 2006

c      print *, ' IDS21 = ',char21
      CALL PLTNOTE(W(401),IDS21)                                   ! 15 Mar 2010
      CALL PLTNOTE(W(426),IDS22)                                   ! 15 Mar 2010
      CALL PLTNOTE(W(451),IDS23)                                   ! 15 Mar 2010
      CALL PLTNOTE(W(476),IDS24)                                   ! 15 Mar 2010

      GO TO 1

C ----------------------------------------------------------------------!-------
      END PROGRAM CONPLT
C ----------------------------------------------------------------------!-------

