C file PROFILE.FOR
C ----------------------------------------------------------------------!-------
! 22 Sep 2011 - Remove auto-scaling on vertical scale for terrain (bottom) plots
! 25 Aug 2011 - Changes recommended by Larry Baker
C 13 Mar 2008 - Corrected length in common block /PCONST/
C 31 Mar 2006 - Increased max number of kinds of profiles from 14 to 15
C 24 Mar 2006 - Changed COMMON block /RAYDEV/
C 17 Mar 2006 - Increased max number of kinds of profiles from 12 to 14
C 28 Feb 2006 - Added CALL SETDSP
C 14 June 2005 - Added W29 full functionality
C 2 September 2004 - Added a test for when to set TLTS(2)='M'
C 4 August 2004 - Added TLTS(2)='M' FOR M/S
C 2 June 2004 - Changed the CALL to SETOVR to a FUNCTION  usage.
C ----------------------------------------------------------------------!-------
      PROGRAM PROFILE
C ----------------------------------------------------------------------!-------
C
C     PLOTS PROFILES OF TEMPERATURE, SOUND SPEED, ETC. 
C     FOR ATMOSPHERIC/OCEANIC MODELS
C
C***** CHANGED TLTS TO CHAR*40 INSTEAD OF 80.

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

      CHARACTER*40 TLTS(4)
      INTEGER STRIM
      PARAMETER (NPTS=200)
      DOUBLE PRECISION WORK(NPTS,2)
C
      PARAMETER (NWARSZ=1000)
      COMMON /WW/ DUM(10), MAXW, W(NWARSZ)

C.....W(83) - W(86) contain latitude and longitude for two points.

      EQUIVALENCE (EARTHR,W(1)),(RAYFNC,W(29))
C
      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
! 24 Mar 2006     COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVGRP,NDEVBIN,NDEVTXT 
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVTXT,NDEVBIN          ! 24 Mar 2006
C
C     PARAMETER DECK "PGROUPS"
      PARAMETER (NCHPG1=11,NWPV=500,NSPGP=NCHPG1+2*NWPV+1)         ! 25 Aug 2011
      PARAMETER (MNGRP=9,MXGRP=69,MXLIST=MXGRP-MNGRP+2)
C
      LOGICAL START
      INTEGER PLTCNT
      INTEGER MAP(2)
C
C.....FPARMS contains values for computing the REAL distance covered by
C     each increment between graph points as well as factors used in
C     in computing the values to be plotted.  FPARMS(5) = PLTCNT.
C
      DOUBLE PRECISION FPARMS(5)
C
C.....XPLT and YPLT contain the minimum speed, maximum speed and the
C     increment for plotting the X and Y axes.
C
      DOUBLE PRECISION XPLT(3), YPLT(3)
C
C     The following causes XPLT(1) = W(92), XPLT(2) = W(93), XPLT(3) = W(94)
C
      EQUIVALENCE (XPLT(1),W(92))
C
C.....PLTC indicates whether a sound or wind speed profile is to be made
C     or whether a single wind component is to be profiled.  If its
C     value is negative, it suppresses a frame advance between plots, thus
C     resulting in a superimposition of plots.  If its value is greater than
C     40, it indicates that points are to be plotted.
C
      EQUIVALENCE (PLTC,W(90))
C
      LOGICAL PLTMODE(12)

      EXTERNAL READW

C.....The following parameters are offsets into PLTMODE, an array which
C     indicates the Plotting MODE.
C
      INTEGER NFILL, INISCX, INISCY, NSCALE, INPLT, INSFRAM, NORMPLT,
     &        ENDFRAM, ENDPLOT, NFULL, NXSNEG, NYSNEG

      PARAMETER (NFILL=1,INISCX=2,INISCY=3,NSCALE=4,INPLT=5,INSFRAM=6,
     &           NORMPLT=7,ENDFRAM=8,ENDPLOT=9,NFULL=10,NXSNEG=11,
     &           NYSNEG=12)

      DATA PLTMODE /12*.FALSE./

C      NFILL
C           CAUSES WORK ARRAY TO BE FILLED.  NORMAL CASE.  SEE PRFPLT
C           SUBROUTINE.
C      INISCX
C           CAUSES X SCALING ARRAY 'XPLT' TO BE INITIALIZED TO +-INFINI
C      INISCY
C           CAUSES Y SCALING ARRAY 'YPLT' TO BE INITIALIZED TO +-INFINI
C      NSCALE
C           CAUSES EXTREME X-Y VALUES TO BE DETERMINED SO THAT NO POINT
C           WILL BE OUT OF RANGE. MANUALLY SPECIFIED VALUES WILL NOT BE
C           OVERWRITTEN PROVIDED THEY ARE SUFFICIENT TO
C           INCLUDE ALL DATA VALUES.
C      INSFRAM
C           INSERT A FRAME IN THE GRAPHICS OUTPUT
C      NORMPLT
C           NORMAL PLOT. THIS WILL INCLUDE AUTOSCALING IF THE INCREMENT
C           VALUE USED FOR AXIS ANNOTATION IS GIVEN AS ZERO. IF EXTREME
C           VALUES ARE GIVEN BUT INCREMENT IS ZERO, AUTOSCALING WILL BE
C           USED TO DETERMINE THE INCREMENT AND SUBTICK INTERVAL
C           BUT WILL NOT DETERMINE EXTREME VALUES.
C      ENDFRAM
C           CAUSES A FINAL END FRAME TO BE OUTPUT TO GRAPHICS FILE.
C      ENDPLOT
C           COMPLETES OUTPUT TO GRAPHICS FILE, TERMINATING PLOTS.
C

      DATA TLTS/4*' '/
      DATA START/.TRUE./

C
C     PERFORM STANDARD INITIALIZATION
C*****CALL SYSINI   MIGHT WANT THIS IF CRAY IS USED.
C*****              BUT WILL ALSO NEED DFSYS.
c      print*,'Beginning of PROFILE'                             ! 26 Mar 2007
      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.

c      print*,'before setdsp'
      CALL SETDSP                                               ! 28 Feb 2006
c      print*,'before setspd'
      CALL SETSPD
c      print*,'before setwnd'
      CALL SETWND
c      print*,'before setrcv'
      CALL SETRCV
c      print*,'before settop'
      CALL SETTOP
c      print*,'before setset'
      CALL SETSUR
c      print*,'before setabs'
      CALL SETABS

      CALL SETSRD
      CALL STDINT
C
C     ENABLE RETURN OF CONTROL FROM READW AFTER EOF
      NFRMAT=-1
      FPARMS(5)=0.0
C
C.....Initialize counter for number of profiles generated.

      numgend = 0

C********* READ W ARRAY VALUES
C
c      print *, 'before readw'
10    IF(READW(xxx).NE.0.d0) THEN
c      print *, 'after readw'
C        LAST RUN SET HAS BEEN READ
C
C*****   RE-SET PLTMODE TO FALSE, SO THAT ONLY ENDFRAM AND ENDPLOT WILL BE TRUE
C*****   WHICH WILL CAUSE PRFPLT ROUTINE TO TAKE PROPER ACTION.  JAH.
C
C        SET MODE TO TERMINATE PLOTTING
C
	 DO 111 I=1,12
111	     PLTMODE(I) = .FALSE.
         PLTMODE(ENDFRAM)=.TRUE.
         PLTMODE(ENDPLOT)=.TRUE.

         GO TO 1000
      ENDIF
c      print *, 'also after readw'
C
      ICODE=ND2B(IDINT(RAYFNC))                               ! 6/14/05
C
C PROCESS PROFILE CALCULATIONS ONLY IF W(29) is one.  In other words,
C this makes it possible to skip a runset by setting W(29) to zero.
C
      IF ((RAYFNC.NE.0.d0).AND.(IAND(ICODE,1).EQ.0)) then      ! 6/14/05
           PRINT*,' WARNING:  NO PROFILE GENERATED FOR THIS RUNSET.'
           PRINT*,'           because W(29) was zero.'
           GO TO 10
      endif
C
C.....
C
C     OVERRIDE ANY STANDARD VALUES IN W-ARRAY
C      CALL SETOVR(xxx)
      yyy = SETOVR(xxx) ! 2 Jun 2004
C

c      print *, 'before PRFSUP'

      CALL PRFSUP
c      print *, 'after PRFSUP'
      PLTCNT=ABS(PLTC)
      FPARMS(5)=PLTCNT
C
      PLTMODE(NORMPLT)=.TRUE.
      PLTMODE(NFILL)=.TRUE.


      IF(START .OR. PLTC.GT.0.0) THEN

C     INITIALLY AND WITH EACH POSITIVE PLOT PARAMETER, REINITIALIZE
C     SCALING.
C     Note that this is only done for normal plots NOT for point plots.
C     Don't want to reset for point plots since point plots need to use
C     the scaling from the input DINP or from the previous line plot.

	 PLTMODE(NORMPLT)=.TRUE.
         PLTMODE(NFILL)=.TRUE.
         PLTMODE(INSFRAM)=.TRUE.
         PLTMODE(INISCX) = .FALSE.
         PLTMODE(INISCY) = .FALSE.
         PLTMODE(NSCALE) = .FALSE.
         START=.FALSE.
         MARK=0

C*****   THE FOLLOWING ARE READ IN FROM THE "W" ARRAY.  SEE EQUIVALENCES
C*****   ABOVE FOR XPLT(1) AND W(92) etc.  JAH.

         XPLTMN = XPLT(1)
         XPLTMX = XPLT(2)
         XPLTIN = XPLT(3)
C
C     USE auto-SCALING IF INDICATED:
C     1. Minimum value >= Maximum value
C     2. Tick mark interval <= 0
C     3. Tick mark interval > 0 and > Maximum value (i.e. user mistake)
C        Handle the possibility that the user enters a tick mark interval
C        value that is greater than the maximum value.  Here, maximum value
C        means the maximum absolute value of the input max and min.  For example,
C        if we didn't take the absolute value, you could have a an increment
C        value of say .01 and a maximum of 0.0 and minimum of -.02 and you
C        would perform auto-scaling when you didn't want to.
C     4. Topography profiles:  PLTCNT = 6.  Must be auto-scaled.  Otherwise,
C        the max and min and increments are in radians and cause plotting
C        problems.
!        22 Sep 2011 - Actually, only the horizontal axis needs to be auto-scaled.  
!        We have removed autoscaling on the vertical axis now.

c      print *, 'before xmaxval'
         xmaxval = MAX(ABS(xpltmx),ABS(xpltmn))
         IF (XPLTMN.GE.XPLTMX .or. xpltin .le. 0.0 .or.
     &       (xpltin .gt. 0.0 .and. xpltin .gt. ABS(xmaxval)) 
     &       .or.pltcnt .eq. 6
     &       ) THEN
            PLTMODE(NSCALE)=.TRUE.
            PLTMODE(INISCX)=.TRUE.
         ENDIF
c      print *, 'after xmaxval'

! 22 Sep 2011         IF(PLTCNT.EQ.6) THEN
! 22 Sep 2011            YPLTMN=PID2-w(85)
! 22 Sep 2011            YPLTMX=PID2-w(83)
! 22 Sep 2011            YPLTIN= w(87)
! 22 Sep 2011         ELSE
            YPLTMN= w(88)
            YPLTMX= w(89)
            YPLTIN= w(96)
! 22 Sep 2011         ENDIF
C
C     USE auto-SCALING IF INDICATED
C
         ymaxval = MAX(ABS(ypltmx),ABS(ypltmn))
         IF(YPLTMN.GE.YPLTMX.OR.YPLTIN.LE.0.0 .or.
     &       (ypltin .gt. 0.0 .and. ypltin .gt. ABS(ymaxval)) 
! 22 Sep 2011     &       .or.pltcnt .eq. 6
     &       ) THEN
	      PLTMODE(NSCALE)=.TRUE.
            PLTMODE(INISCY)=.TRUE.
	   ENDIF

      ENDIF

C.....For point plots, do not advance frame.
C     Note, however, that all the other elements of PLTMODE are unchanged.

      if (pltc .lt. 0.0) then
	   PLTMODE(INSFRAM) = .FALSE.
      endif
C
      NP=NPTS

C.....If the upper limit of the Y-axis is =< 0, then this is an ocean
C     profile, otherwise it includes (at least part of) the atmosphere.
C     Also, if the entire profile is below the surface, then the max
C     and min values are changed to be positive.

      IF (ypltmx.LE.0.0) THEN
           PLTMODE(NYSNEG)=.TRUE.
           CALL ANNFIL('D#EPTH!',TLTS(3))
      ELSE
           CALL ANNFIL('H#EIGHT ABOVE! MSL',TLTS(3))
      ENDIF
C
C     SET UP MAPPING TO REVERSE X AND Y ORIENTATION
C     SINCE VERTICAL IS INDEPENDENT COORDINATE
C
c      print *, 'before map'
      MAP(1)=2
      MAP(2)=1
      FPARMS(3)=1.0
c      print *, 'after map'
      IF(PLTCNT.EQ.6) THEN
         FPARMS(1)=PID2-w(83)
         FPARMS(2)=PID2-w(85)
C        SET MAPPING FOR STANDARD X-Y ORIENTATION
         MAP(1)=1
         MAP(2)=2
C........If W(300) is 4, then the GLORENZ model was chosen.
C        This model is an east-west ridge and the height of the base of
C        the ridge is given by W(306).   The base can be above or below
C        sea level.
C        If some other bottom/terrain model was chosen, then assume
C        it is above sea level.  The graph is still usable.

         if (w(300) .eq. 4.0) then

              if (w(306) .lt. 0.0) then
                   CALL ANNFIL('D#EPTH!',TLTS(3))
              else
                   CALL ANNFIL('H#EIGHT!',TLTS(3))
              endif
         else
              CALL ANNFIL('H#EIGHT!',TLTS(3))
         endif
c      print *, 'before fparms'
         FPARMS(3)=EARTHR
      ELSEIF(pltcnt .ge. 1 .and. pltcnt .le. 15) THEN            ! 31 Mar 2006
         FPARMS(1)=ypltmn+EARTHR
         FPARMS(2)=ypltmx+EARTHR
      ELSEif (pltcnt .ge. 41 .and. pltcnt .le. 62) then
C
C........FOR POINT PLOTS USE LINEAR INDEXING
C*****   NOTE:  IPTPLT is an entry point in the function PRFSUP.
C
         FPARMS(1)=0.
         PLTCTMP = ABS(PLTC)
         NP=IPTPLT(PLTCTMP,TLTS(1),TLTS(2),TLTS(4))
         FPARMS(2)=NP-1
C
C........Check to see if the vertical tick interval is 0.0 which
C        causes interval to be auto-calc'd AND causes interval to be in meters.
C        If so, need to scale the point plot data by 1000 to convert it
C        to meters (remember that, no matter what units the point plot data
C        uses in the DINP file, the internal units of the program are KMs).

         if (w(96) .eq. 0.0) fparms(3) = fparms(3)*1000.0

      else
C.....Warn the user that W(90) needs to be set and go to next runset.
           PRINT*,' WARNING:  NO PROFILE GENERATED FOR THIS RUNSET.'
           PRINT*,' Invalid W(90) value.',                       ! 31 Mar 2006
     &            '  Valid (absolute) values are 1 to 15 and 41 to 62.'
           PRINT*,' Absolute value of W(90) = ',pltcnt
           go to 10
      ENDIF
c      print *, 'after endif'
C
C.....TLTS(3) is the Y-axis label.
C     IF DEPTH(HEIGHT) RANGE IS LARGE ENOUGH USE KILOMETERS in the Y-axis label,
C     otherwise, use meters.

      IF(START .OR. PLTC.GT.0.0) THEN
         IF( (YPLTMX-YPLTMN)*FPARMS(3).GT.10.
     &      .OR. YPLTIN.NE.0.) then
            CALL ANNFIL(TLTS(3)(:STRIM(TLTS(3)))//' (#KM!)',TLTS(3))
         ELSE
            FPARMS(3)=FPARMS(3)*1000.
            CALL ANNFIL(TLTS(3)(:STRIM(TLTS(3)))//' (#M!)',TLTS(3))
         ENDIF
      ENDIF
C
      FPARMS(4)=1.

C.....TLTS(2) is the X-axis label.  This label is later extended with other
C     characters in PRFSUP.

      TLTS(2)='KM'
      IF(XPLT(1).EQ.XPLT(2).OR.XPLT(3).EQ.0.D0.OR.XPLT(2).GT.1.D0)
     1 TLTS(2)='M'      ! 2 Sept 2004
C
C.....Set MARK to negative for point plots.  Indicates that symbols
C     will be drawn at each point.

      MRK=MARK
      IF(PLTCNT .gt. 40) MRK=-MARK
C
C.....Increment number of profiles generated.  Note that it is incremented
C     here, before the "1000 continue" statement.  Do not want to increment
C     for the last call to PLOTAL.

      numgend = numgend + 1
c      print *, '1000 continue'

1000  CONTINUE
c      print *, 'after 1000 continue'
C.....Make sure any changes are in the arrays before being passed out.
C
      XPLT(1) = XPLTMN
      XPLT(2) = XPLTMX
      XPLT(3) = XPLTIN
      YPLT(1) = YPLTMN
      YPLT(2) = YPLTMX
      YPLT(3) = YPLTIN

c      print *, 'before PLOTAL'

      CALL PLOTAL
     &           (PLTMODE,PLTCNT,TLTS,NP,MAP,WORK,FPARMS,XPLT,YPLT,MRK)
      MARK=MARK+1
c      print *, 'after PLOTAL'
C
C.....If the plotting mode has not been set to end frame and end plot
C     then continue to read data sets.
C
      IF (.NOT.PLTMODE(ENDFRAM) .AND. .NOT. PLTMODE(ENDPLOT)) THEN
           GO TO 10
      ENDIF

      PRINT*,' '
      if (numgend .gt. 0) then
           PRINT*,' PROFILE Completed Normally.'
           PRINT*,' '
           PRINT 222, numgend
222        FORMAT('  Number of Profiles Generated: ',I2)
      else
           PRINT*,' WARNING:  No Profiles were Generated.'
           PRINT*,'           Check input data file for invalid values.'
      endif
      PRINT*,' '
      PRINT*,' Output File:  TXTOUT.'
C
C ----------------------------------------------------------------------!-------
      END PROGRAM PROFILE
C ----------------------------------------------------------------------!-------
