      PROGRAM PLOT_SKEWT
C
C---------------------------------------------------------------------------+
C                         SKEWT Plotting Program                            +
C                                                                           +
C                            David Blanchard                                +
C                         NOAA/NSSL & UCAR/COMET                            +
C                                                                           +
C This Fortran code plots a Skew-T Log-P sounding diagram.  It requires     +
C NCAR GKS V3.0 or higher.  Subroutine SKEWTB expects arrays of pressure,   +
C temperature, dewpoint, wind direction, and wind speed.  The thermodynamic +
C arrays use the NC counter; the wind arrays use the NN counter.            +
C                                                                           +
C---------------------------------------------------------------------------+
C
C Parameter statements set to 7 levels.  Resize accordingly.
C
      PARAMETER (NC=7,NN=7)
      DIMENSION PRES(NC),DWPT(NC),TEMP(NC)
      DIMENSION PW(NN),DIR(NN),SPD(NN)
      CHARACTER ISTAT*20, IDATE*20
C
C Canned data for demonstration purposes.
C
      DATA PRES / 1000., 850., 700., 500.,  300.,  200.,  100./
      DATA TEMP /   30.,  20.,  10., -10.,  -40.,  -60.,  -70./
      DATA DWPT /   15.,  10.,   4., -20.,  -60., -999., -999./
C
      DATA PW   / 1000., 850., 700., 500.,  300.,  200.,  100./
      DATA DIR  /   90., 180., 200., 240.,  270.,  270.,  200./
      DATA SPD  /    8.,  15.,  20.,  25.,   30.,   35.,   20./
C
      DATA ISTAT  /'DEN'/
      DATA BADVAL /-999./
      DATA IDATE  /'01-JAN-2001 00'/
C
C  Open up the GKS graphics package.
C
      IWKID1=1
      CALL GOPKS(6,0)
      CALL GOPWK(IWKID1,2,1)
      CALL GACWK(IWKID1)
C
C Define the color table.
C
      CALL GSCR(IWKID1, 0, 0.0, 0.0, 0.0)		!black
      CALL GSCR(IWKID1, 1, 1.0, 1.0, 1.0)		!white
      CALL GSCR(IWKID1, 2, 1.0, 0.0, 0.0)		!red
      CALL GSCR(IWKID1, 3, 0.0, 1.0, 0.0)		!green
      CALL GSCR(IWKID1, 4, 0.0, 0.0, 1.0)		!blue
      CALL GSCR(IWKID1, 5, 0.0, 1.0, 1.0)		!cyan
      CALL GSCR(IWKID1, 6, 0.3, 0.3, 0.3)		!dark grey
      CALL GSCR(IWKID1, 7, 0.7, 0.7, 0.7)		!medium grey
C
C Set up basic needs: font; clipping, etc.
C
      CALL GSCLIP(0)		! Clipping off
      CALL GSTXFP(12,2)		! Helvetica
C
C-----------------------------------------------------------------
C Insert your code here to read in data and prepare it for plotting
C instead of using the canned data from the DATA statements.
C-----------------------------------------------------------------
C
C Draw the sounding and advance the frame.
C
      CALL SKEWTPL(PRES,TEMP,DWPT,NC,PW,SPD,DIR,NN,IDATE,ISTAT,BADVAL)
      CALL FRAME
C
C  Close GKS
C
      CALL GDAWK(IWKID1)
      CALL GCLWK(IWKID1)
      CALL GCLKS
C
C That's all, folks!
C
      CALL EXIT(0)
      END
C***********************************************************************
      SUBROUTINE SKEWTPL(PRES, TEMP, DWPT, NC, PW, SPD, DIR, NN,
     1                  IDATE, ISTAT, BADVAL)
C
C Carl Mohr NCAR
C    Original version.  Date unknown (early 70's)
C
C David Blanchard NOAA/WRP   Fall 1985
C    Fixed labels for Theta.  Added wind barbs (instead of vectors)
C    Numerous cosmetic changes.
C
C David Blanchard NOAA/NSSL & UCAR/COMET  11/07/93
C    Converted from old NCARG to NCARG GKS
C
C----------------------------------------------------
C PRES   - pressure array for thermodynamic data (mb)
C TEMP   - temperature array (celsius)
C DWPT   - dew point array (celsius)
C NC     - number of thermodynamic points
C PW     - pressure array for wind data (mb)
C SPD    - wind speed array (m/s)
C DIR    - wind direction array (degrees-0 is north)
C NN     - number of wind data points
C IDATE  - char*20 of date/time
C ISTAT  - char*20 word containing station name
C BADVAL - value assigned to missing data  --temp,dwpt traces are
C         terminated as soon as this number is encountered.
C
C Output parameters.....all input parameters remain unchanged by
C this routine.
C
      DIMENSION PRES(1),TEMP(1),DWPT(1),SPD(1),DIR(1),PW(1)
C
C  Text buffer
C
      CHARACTER TEXT*80, ISTAT*3, IDATE*14
C
C  Skewt border
C
      DIMENSION XB(7),YB(7)
      DATA XB/-19.,27.1,27.1,18.6,18.6,-19.,-19./
      DATA YB/-.9346217,-.9346217,9.,17.53,44.061,44.061,-.9346217/
C
C  Pressure line specs
C
      DIMENSION PLV(11),PLN(11,2)
      DATA PLV/100.,200.,300.,400.,500.,600.,700.,800.,900.,1000.,1050./
      DATA PLN/11*-19.,4*18.6,22.83,26.306,5*27.1/
C
C  Temperature line specs
C
      DIMENSION TP(15,2)
      DATA TP/8*1050.,855.,625.,459.,337.,247.,181.,132.,730.,580.,500.,
     &      430.,342.,251.,185.,135.,7*100./
C
C  Mixing ratio specs
C
      DIMENSION RAT(8)
      CHARACTER LRAT(8)*2
      DATA RAT/20.,12.,8.,5.,3.,2.,1.,0.4/
      DATA LRAT/'20','12',' 8',' 5',' 3',' 2',' 1','.4'/
C
C Dry/saturated adiabat buffers
C
      DIMENSION SX(162),SY(162),Y45(162)
C
C Degrees to radians, absolute zero
C
      DATA DTR,ABZ/ 0.0174532925, 273.16 /
C
C A wind vector of speed VSCALE m/s will be 2.5 cm in length on skewt
C XM is x-location on skewt where wind vectors are drawn
C
      DATA VSCALE /20./
      DATA VLONG/ 2.5 /
      DATA XM / 22.5 /
C
C Logical definitions
C
      LOGICAL FIRST
      DATA FIRST /.TRUE./
C
C Mappings from (P,T) to cm on skewt
C
      FY(P)=132.182-44.061*ALOG10(P)
      FX(T,Y)=0.54*T+0.90692*Y
C
C    Draw skewt border
C
      WSCALE = VLONG/VSCALE
      XM1 = XM + VLONG
      CALL SET(.05,.95,.05,.95,-19.0,27.1,-.9346217,44.061,1)
      CALL GSPLCI(7)
      CALL CURVE(XB,YB,7)
C
C  Draw the pressure lines
C
      CALL GSPLCI(7)
      CALL GSTXCI(7)
      DO 40 K=1,11
         Y1=FY(PLV(K))
         IF(K.NE.1.AND.K.NE.11) CALL LINE(PLN(K,1),Y1,PLN(K,2),Y1)
         ITS=NINT(PLV(K))
         WRITE(TEXT,20) ITS
20       FORMAT(I4)
         CALL PLCHLQ(-20.5,Y1,TEXT(1:4),12.,0.,0.)
40    CONTINUE
C
C Draw temperature lines
C
      CALL GSPLCI(7)
      T=40.
      DO 120 I=1,15
         Y1=FY(TP(I,1))
         Y2=FY(TP(I,2))
         X1=FX(T,Y1)
         X2=FX(T,Y2)
         CALL LINE(X1,Y1,X2,Y2)
         ITS=NINT(T)
         X2=X2+0.4
         Y2=Y2+0.441
         WRITE(TEXT,20) ITS
         CALL PLCHLQ(X2,Y2,TEXT(1:4),10.,45.,0.)
100      T=T-10.
120   CONTINUE
C
C Tick marks at 500 mb
C
      CALL GSPLCI(7)
      Y1=13.2627
      Y2=13.75
      T=-52.
      DO 140 I=1,31
         T=T+2.
         IF(AMOD(T,10.).EQ.0.)GO TO 140
         X1=FX(T,Y1)
         X2=FX(T,Y2)
         CALL LINE(X1,Y1,X2,Y2)
140   CONTINUE
C
C Draw mixing ratio lines
C
      CALL SFLUSH
      CALL GSPLCI(7)
      CALL DASHDB('7417'O)
      Y1=FY(1050.)
      Y2=FY(700.)
      DO 160 I=1,8
         X1=FX(TMR1(RAT(I),1050.)-ABZ,Y1)
         X2=FX(TMR1(RAT(I), 700.)-ABZ,Y2)
         CALL LINED(X1,Y1,X2,Y2)
         CALL PLCHLQ(X2,Y2+0.6,LRAT(I),10.,0.,1.)
160   CONTINUE
C
C Draw saturated adiabats
C
      CALL SFLUSH
      CALL GSPLCI(7)
      CALL DASHDB('70707'O)
      TS=32.
      DO 220 I=1,7
         P=1060.
         TK=TS+ABZ
         AOS=OS1(TK,1000.)
         DO 180 J=1,86
            P=P-10.
            ATSA=TSA1(AOS,P)-ABZ
            SY(J)=FY(P)
            SX(J)=FX(ATSA,SY(J))
180      CONTINUE
         CALL CURVED(SX,SY,86)
         ITS=NINT(TS)
         IF(ITS .GE. 20) THEN
           WRITE(TEXT,200) ITS
200        FORMAT(I2)
           CALL PLCHLQ(SX(86),SY(86)+0.6,TEXT(1:2),10.,0.,1.)
         END IF
         TS=TS-4.0
220   CONTINUE
C
C Draw dry adiabat lines
C
      CALL SFLUSH
      CALL GSPLCI(7)
      CALL DASHDB('4444'O)
      T=51.
      DO 240 I=1,162
         Y45(I)=66.67*(5.7625544-ALOG(T+ABZ))
         T=T-1.0
240   CONTINUE
      T=450.
      TD=52.
      DO 320 I=1,20
         T=T-10.
         K=0
         YD=66.67*(ALOG(T)-5.7625544)
         DO 260 J=1,162
            YPD=Y45(J)+YD
            TX=TD-FLOAT(J)
            IF(YPD.GT.44.061) GO TO 280
            IF(YPD.LT.-.9346217) GO TO 260
            XPD=FX(TX,YPD)
            IF(XPD.LT.-19.0)GO TO 280
            IF(XPD.GT.27.1)GO TO 260
            IF(XPD.GT.18.6.AND.T.GT.350.0)GO TO 260
            K=K+1
            SX(K)=XPD
            SY(K)=YPD
260      CONTINUE
C
280      CALL CURVED(SX,SY,K)
         ITS=NINT(T)
         WRITE(TEXT,300) ITS
300      FORMAT(I3)
C
         IF(ITS .GE. 320) THEN
           CALL PLCHLQ(SX(K-3),43.0,TEXT(1:3),10.,0.,1.)
         ELSE
           CALL PLCHLQ(-17.5,SY(K-3),TEXT(1:3),10.,0.,1.)
         END IF
320   CONTINUE
C
C Skewt background has been generated-- plot the sounding
C
      CALL SFLUSH
      CALL GSTXCI(3)
      CALL GSPLCI(3)
      CALL PLCHLQ(-18.5,-2.4,ISTAT,20.,0.,-1.)
      CALL PLCHLQ(-01.0,-2.4,IDATE,20.,0.,-1.)
C
C  Solid dash pattern, increased spot size
C
      CALL SFLUSH
      CALL DASHDB('177777'O)
      CALL GSLWSC(2.)
C
C Draw temperature curve
C
      CALL GSPLCI(3)
      DO 340 I=1,NC
         IF(TEMP(I).EQ.BADVAL)GO TO 360
            Y=FY(PRES(I))
            X=FX(TEMP(I),Y)
            IF(I.EQ.1)CALL FRSTPT(X,Y)
            CALL VECTOR(X,Y)
340   CONTINUE
360   CONTINUE
C
C Draw dewpoint curve
C
      CALL GSPLCI(3)
      DO 380 I=1,NC
         IF(DWPT(I).EQ.BADVAL)GO TO 400
         IF(PRES(I).LE. 250.) GO TO 400
         Y=FY(PRES(I))
         X=FX(DWPT(I),Y)
         IF(I.EQ.1)CALL FRSTPT(X,Y)
         CALL VECTOR(X,Y)
380   CONTINUE
400   CONTINUE
C
C Plot wind vectors
C
      CALL SFLUSH
      CALL GSPLCI(7)
      IF(NN.GT.0) THEN
         CALL LINE(XM,-.9346217,XM,44.061)
         CALL SFLUSH
         CALL GSPLCI(3)
         CALL GSLWSC(2.)
         PWOLD = 1050.	!Set first level to "sea-level"
         SPACE=0.03
         DO 420 I=1,NN
            PDIFF = LOG(PWOLD)-LOG(PW(I)) !Avoid overcrowding of winds
            IF(DIR(I) .LE. 360. .AND. PDIFF .GT. SPACE) THEN
               PWOLD = PW(I)
               ANG=DIR(I)*DTR
               U = -SPD(I)*SIN(ANG)
               V = -SPD(I)*COS(ANG)
               Y1=FY(PW(I))
               CALL WNDBARB(XM,Y1,U,V)
            END IF
420      CONTINUE
      ENDIF
C
C  Reset to normal spot size and return
C
      CALL SFLUSH
      CALL GSPLCI(7)
      CALL GSLWSC(1.)
C
      RETURN
      END
C...........................................................
      FUNCTION  ESAT1(T)
C  ESAT(MILLIBARS),T(KELVIN)
      DATA ABZ/273.16/
      TC=T-ABZ
      ESAT1=6.1078*EXP((17.2693882*TC)/(TC+237.3))
      RETURN
      END
C...........................................................
      FUNCTION   OS1(T,P)
C  OS AND T (KELVIN) , P (MILLIBARS )
      OS1=T*((1000./P)**.286)/(EXP(-2.6518986*W1(T,P)/T))
      RETURN
      END
C...........................................................
      FUNCTION  TMR1(W,P)
C  TMR(KELVIN),W(GRAMS WATER VAPOR/KILOGRAM DRY AIR),P(MILLIBAR)
C  ALOG10  15   LOG TO THE BASE  TEN.
      X =  ALOG10(   W*P/(622.+ W)  )
      TMR1=10.**(.0498646455*X+2.4082965)-7.07475+38.9114*
     X                   ((10.**( .0915*X ) - 1.2035 )**2 )
      RETURN
      END
C...............................................................
      FUNCTION TSA1(OS,P)
C  TSA AND OS(KELVIN),P(MILLIBARS)
C  SIGN(A,B) REPLACES THE ALGEBRAIC SIGN OF A WITH THE SIGN OF B
          A  =  OS
          TQ  =  253.16
          D =  120
          DO 20  I = 1,12
             D = D/2.
C  IF THE TEMPERATURE DIFFERENCE,X, IS SMALL,EXIT THIS LOOP
             X=A*EXP(-2.6518986*W1(TQ,P)/TQ)-TQ*((1000./P)**.286)
             IF(ABS(X).LT.0.01)GO TO 40
             TQ = TQ + SIGN(D,X)
20        CONTINUE
40    TSA1=TQ
      RETURN
      END
C...............................................................
      FUNCTION  W1(T,P)
C  W(GRAMS WATER VAPOR/KILOGRAM DRY AIR ), P(MILLIBAR )
      IF(T.GE.999.)GO TO 20
         X  =       ESAT1(T)
         W1 =  621.97*X/(P-X)
      RETURN
20    W1=0.0
      RETURN
      END
C****************************************************************
      SUBROUTINE WNDBARB(X,Y,U,V)
C
C J. Meitin      Original version (c. 1976)
C I. Watson      Various improvements
C B. Stankov     Added flag fill code
C D. Blanchard   Converted to NCAR GKS code
C
C Windbarbs are plotted in knots, yet the U and V components are
C entered in meters per seconds (as they should be).  This
C subroutine generates the classic windbarbs plots with a shaft
C and the appropriate feathers to represent 5 and 10 knot
C increments.
C
C  X,Y - Coordinates of origin of windbarbs
C  U,V - Wind components
C  SC  - Length of windshaft in crt plotting units
C  PF  - Full length of feather as a fraction of the shaft length
C  PT  - Distance between feathers as a fraction of the shaft length
C  BA  - Value at which triangle replaces feathers
C  BB  - Value for which a full feather is drawn
C      - half a feather corresponds to half of BB
C
      DATA SC/75./, PF/.4/, PT/.2/
      DATA BA/50./, BB/10./
C
C First point is along vertical staff
C
      CALL FRSTPT(X,Y)
      IF(U.EQ.0..AND.V.EQ.0.)RETURN
      CALL MXMY(M,N)
      XN=SQRT(U*U+V*V)*1.94
      MX=-SC*U/XN
      NY=-SC*V/XN
      IA=PF*NY
      JA=-PF*MX
      MMX=MX+M
      NNY=NY+N
      CALL VECTOR(CPUX(MMX),CPUY(NNY))
C
C Shaft is drawn.  Draw feathers/pennant
C
      PC=1.
      XT=XN+.25*BB
  5   IF(XT .LT. BA) GO TO 15
      I=PC*MX+M+.5
      J=PC*NY+N+.5
      CALL FRSTPT(CPUX(I),CPUY(J))
      II=I+IA
      JJ=J+JA
      CALL VECTOR(CPUX(II),CPUY(JJ))
      PC=PC-PT
      I=PC*MX+M+.5
      J=PC*NY+N+.5
      CALL VECTOR(CPUX(I),CPUY(J))
C******************************** Fill Barbs starts here
      PC=PC+PT*7./8.
      I=PC*MX+M+.5
      J=PC*NY+N+.5
      CALL FRSTPT(CPUX(I),CPUY(J))
      CALL VECTOR(CPUX(II),CPUY(JJ))
      PC=PC-PT*6./8.
      I=PC*MX+M+.5
      J=PC*NY+N+.5
      CALL VECTOR(CPUX(I),CPUY(J))
      PC=PC+PT*5./8.
      I=PC*MX+M+.5
      J=PC*NY+N+.5
      CALL FRSTPT(CPUX(I),CPUY(J))
      CALL VECTOR(CPUX(II),CPUY(JJ))
      PC=PC-PT*4./8.
      I=PC*MX+M+.5
      J=PC*NY+N+.5
      CALL VECTOR(CPUX(I),CPUY(J))
      PC=PC+PT*3./8.
      I=PC*MX+M+.5
      J=PC*NY+N+.5
      CALL FRSTPT(CPUX(I),CPUY(J))
      CALL VECTOR(CPUX(II),CPUY(JJ))
      PC=PC-PT*2./8.
      I=PC*MX+M+.5
      J=PC*NY+N+.5
      CALL VECTOR(CPUX(I),CPUY(J))
      PC=PC+PT*1./8.
      I=PC*MX+M+.5
      J=PC*NY+N+.5
      CALL FRSTPT(CPUX(I),CPUY(J))
      CALL VECTOR(CPUX(II),CPUY(JJ))
      PC=PC-PT/2.
      I=PC*MX+M+.5
      J=PC*NY+N+.5
      CALL VECTOR(CPUX(I),CPUY(J))
      PC=PC-PT/2.
C*********************************  Fill Barbs ends here.
C     PC=PC-PT
      XT=XT-BA
      GO TO 5
  15  CONTINUE
      IF(XT .LT. BB) GO TO 25
      I=PC*MX+M+.5
      J=PC*NY+N+.5
      CALL FRSTPT(CPUX(I),CPUY(J))
      II=I+IA
      JJ=J+JA
      CALL VECTOR(CPUX(II),CPUY(JJ))
      PC=PC-PT
      XT=XT-BB
      GO TO 15
  25  IF(XT .LT. .5*BB) GO TO 30
      I=PC*MX+M
      J=PC*NY+N
      CALL FRSTPT(CPUX(I),CPUY(J))
      II=I+.5*IA
      JJ=J+.5*JA
      CALL VECTOR(CPUX(II),CPUY(JJ))
  30  CONTINUE
      RETURN
      END
