* * * * ncargf77 skewt.f rdsnd.f -o asc_to_skewt * program skewtest c c From gill@ncar.UCAR.EDU Tue Nov 17 17:28:03 1992 c c ... how many values, here we make the assumption that c winds and temp are at the same level c PARAMETER(IERRF=6, LUNIT=2, IWKID=1, IWTYPE=20) parameter (nn=20000) dimension pres(nn),temp(nn),dwpt(nn),spd(nn),dir(nn) COMMON /SKWDRW/ ISKDRW character *80 title c c--------------------------------------------------------------- c pres pressure (mb) c temp temperature (C) c dwpt dew point temperature (C) c spd wind speed (m/s) c dir wind direction (degrees) c c you can choose the value of a flag to mark suspect c data *BADVAL* c c wind speeds on the plot have the following interpretation c full pennant 50 kts c full barb 10 kts c half barb 5 kts c--------------------------------------------------------------- c c ... flag that says activate the background, modified internally c iskdrw=0 c c c open gks, shut off clipping, and allow for flashing c * CALL OPNGKS * CALL GSCLIP (0) * CALL GOPWK (2,9,3) C C OPEN GKS, OPEN WORKSTATION OF TYPE 1, ACTIVATE WORKSTATION C CALL GOPKS (IERRF, ISZDM) CALL GOPWK (IWKID, LUNIT, IWTYPE) CALL GACWK (IWKID) * call gsclip (0) * write(*,*)' ' * write(*,*)' Where is the input file?' * read(*,200)infile * 200 format(a60) * open(unit=9,file='sounding.txt') iunit = 9 * call rdsnd(iunit,pres,temp,dwpt,spd,dir,nn,kk,title,badval) * c c ... pass this routine the data, it does it all c write(*,*)' ' write(*,*)' Getting ready to create the plot...' call skewt (pres,temp,dwpt,spd,dir,nn,kk,title,badval) write(*,*)' Plot completed!' write(*,*)' ' * call frame * C C CLOSE GKS C CALL GDAWK (IWKID) CALL GCLWK (IWKID) CALL GCLKS * c c ... shut down GKS c * CALL GCLWK (2) ! CLOSE GFLASH FOR MAP AND SKEW-T BACKGROUND * CALL CLSGKS c c ... that's all folks c stop 99999 end SUBROUTINE SKEWT (PRES,TEMP,DWPT,SPD,DIR,NN,KK,TITLINE,BADVAL) C C SKEWT- PLOTS SOUNDINGS ON A SKEWT, LOG P THERMODYNAMIC DIAGRAM C C PRES- PRESSURE ARRAY FOR THERMODYNAMIC DATA (MB) C TEMP- TEMPERATURE ARRAY (CELSIUS) C DWPT- DEW POINT ARRAY (CELSIUS) C SPD- WIND SPEED ARRAY (M/S) C DIR- WIND DIRECTION ARRAY (DEGREES-0 IS NORTH) C NN- MAX NUMBER OF DATA LEVELS C KK- ACTUAL NUMBER OF DATA LEVELS C BADVAL- VALUE ASSIGNED TO MISSING DATA --TEMP,DWPT TRACES ARE C TERMINATED AS SOON AS THIS NUMBER IS ENCOUNTERED. C C OUTPUT PARAMETERS..... C ALL INPUT PARAMETERS REMAIN UNCHANGED BY THIS ROUTINE. C DIMENSION PRES(1),TEMP(1),DWPT(1),SPD(1),DIR(1) CHARACTER LAB*120, ISTAT*80 C CHARACTER *80 TITLINE INTEGER ITSTRT, * ITITLEN, * LCNTR COMMON /SKWDRW/ ISKDRW C C DEGREES TO RADIANS PARAMETER (DTR = 0.0174532925) C WIND BARB DATA PARAMETER (XM = 24.2) C C LINE WIDTH VARIABLES C * DATA ISPOPT / 'LW' / DATA ISWIDE / 2000 / 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 TEST TO SEE IF A BACKGROUND HAS BEEN DRAWN, IF NOT CALL SKWTBKG C IF (ISKDRW .EQ. 0) THEN CALL SKWTBKG ISKDRW = 1 END IF C C SKEWT BACKGROUND HAS BEEN GENERATED-- PLOT THE SOUNDING C * CALL GFLAS3 (1) CALL SET(.05,.95,.05,.95,-19.0,27.1,-.9346217,44.061,1) C C PUT ON TITLE ITITLEN=LEN(TITLINE) CALL WTSTR (4.05,-2.4,TITLINE(1:ITITLEN),12,0,0) IF(KK.GT.0) CALL LINE(XM,-.9346217,XM,44.061) C C SOLID DASH PATTERN, INCREASED SPOT SIZE (DEFAULT=8) C CALL DASHDB (65535) * CALL GETUSV (ISPOPT,ISNORM) CALL GETUSV ('LW',ISNORM) ISWIDE = ISWIDE + 3000 * CALL SETUSV (ISPOPT,ISWIDE) CALL SETUSV ('LW',ISWIDE) ISWIDE = ISWIDE - 3000 C CALL SETUSV('LW',6000) * DO 60 I=1,NN DO 60 I=1,KK IF(TEMP(I).EQ.BADVAL)GO TO 61 Y=FY(PRES(I)) X=FX(TEMP(I),Y) IF (I.EQ.1) CALL FRSTPT(X,Y) CALL VECTOR(X,Y) 60 CONTINUE 61 CONTINUE CALL SETUSV('LW',4000) * DO 70 I=1,NN DO 70 I=1,KK IF(DWPT(I).EQ.BADVAL)GO TO 71 Y=FY(PRES(I)) X=FX(DWPT(I),Y) IF(I.EQ.1)CALL FRSTPT(X,Y) CALL VECTOR(X,Y) 70 CONTINUE 71 CONTINUE CALL SETUSV('LW',2000) C C PLOT WIND VECTORS C * IF (NN.LE.0) GO TO 76 IF (KK.LE.0) GO TO 76 C * CALL SETUSV (ISPOPT,ISWIDE) CALL SETUSV ('LW',ISWIDE) ** DO 75 I=1,NN * DO 75 I=1,KK * NTHIN = KK/25 NTHIN = 50 ! Digicora RS-80 NTHIN = 80 ! IMetOS version 1 NTHIN = 10 ! IMetOS-II write(*,*)' NTHIN=',NTHIN DO 75 I=2,KK,NTHIN IF (DIR(I) .GT. 360.) GO TO 75 ANG=DIR(I)*DTR U = -SPD(I)*SIN(ANG) V = -SPD(I)*COS(ANG) Y1=FY(PRES(I)) CALL WNDBARB (XM,Y1,U,V) 75 CONTINUE 76 CONTINUE C C RESET TO NORMAL SPOT SIZE AND EXIT C * CALL SETUSV (ISPOPT,ISNORM) CALL SETUSV ('LW',ISNORM) C NUMFRAME=NUMFRAME+1 CALL FRAME C RETURN END SUBROUTINE SKWTBKG C C SKWTBKG- PLOTS BACKGROUND FOR A SKEWT, LOG P THERMODYNAMIC DIAGRAM C COMMON /SKWDRW/ ISKDRW C C ENCODE BUFFER C CHARACTER*4 ITIT 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 REAL RAT(8) CHARACTER*2 LRAT(8) 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 ADAIBAT BUFFERS C DIMENSION SX(162),SY(162),Y45(162) C C DEGREES TO RADIANS, ABSOLUTE ZERO C PARAMETER (ABZ = 273.16) 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 * CALL GFLAS1 (1) CALL SET(.05,.95,.05,.95,-19.0,27.1,-.9346217,44.061,1) CALL CURVE(XB,YB,7) C C DRAW THE PRESSURE LINES C DO 10 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 (ITIT,101) ITS 101 FORMAT(I4) C CALL PWRY(-20.9,Y1,ITIT,4,1.9,0,0) CALL WTSTR (-19.2,Y1,ITIT,11,0,1) 10 CONTINUE C C DRAW TEMPERATURE LINES C T=40. DO 20 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) IF(ITS.EQ.20) GO TO 19 X2=X2+0.4 Y2=Y2+0.441 WRITE (ITIT,101) ITS C CALL PWRY(X2,Y2,ITIT,4,12,.83422,0) CALL WTSTR (X2,Y2,ITIT,12,47,-1) 19 T=T-10. 20 CONTINUE C C TICK MARKS AT 500 MB C Y1=13.2627 Y2=13.75 T=-52. DO 25 I=1,31 T=T+2. IF(AMOD(T,10.).EQ.0.)GO TO 25 X1=FX(T,Y1) X2=FX(T,Y2) CALL LINE(X1,Y1,X2,Y2) 25 CONTINUE C C DRAW MIXING RATIO LINES C CALL DASHDB (3855) ! PATTERN = 0000111100001111 Y1=FY(1050.) Y2=FY(700.) DO 30 I=1,8 X1=FX(TMR(RAT(I),1050.)-ABZ,Y1) X2=FX(TMR(RAT(I), 700.)-ABZ,Y2) CALL LINED(X1,Y1,X2,Y2) C CALL PWRY(X2,Y2+0.6,LRAT(I),2,10,0,1) CALL WTSTR (X2,Y2+0.6,LRAT(I),10,0,0) 30 CONTINUE C C DRAW SATURATED ADIABATS C CALL DASHDB (31710) ! PATTERN = 0111101111011110 TS=32. DO 40 I=1,7 P=1060. TK=TS+ABZ AOS=OS(TK,1000.) DO 35 J=1,86 P=P-10. ATSA=TSA(AOS,P)-ABZ SY(J)=FY(P) SX(J)=FX(ATSA,SY(J)) 35 CONTINUE CALL CURVED(SX,SY,86) ITS=NINT(TS) WRITE (ITIT,102) ITS 102 FORMAT(I2) C CALL PWRY(SX(86),SY(86)+0.6,ITIT,2,10,0,1) CALL WTSTR (SX(86),SY(86)+0.6,ITIT(1:2),10,0,0) TS=TS-4.0 40 CONTINUE C C DRAW DRY ADIABAT LINES C CALL DASHDB (21845) ! PATTERN = 0101010101010101 C CALL DASHD(4444B) T=51. DO 45 I=1,162 Y45(I)=66.67*(5.7625544-ALOG(T+ABZ)) T=T-1.0 45 CONTINUE T=450. TD=52. DO 55 I=1,20 T=T-10. K=0 YD=66.67*(ALOG(T)-5.7625544) DO 50 J=1,162 YPD=Y45(J)+YD TX=TD-FLOAT(J) IF(YPD.GT.44.061) GO TO 54 IF(YPD.LT.-.9346217) GO TO 50 XPD=FX(TX,YPD) IF(XPD.LT.-19.0)GO TO 54 IF(XPD.GT.27.1)GO TO 50 IF(XPD.GT.18.6.AND.T.GT.350.0)GO TO 50 K=K+1 SX(K)=XPD SY(K)=YPD 50 CONTINUE C 54 CALL CURVED(SX,SY,K) ITS=NINT(T) WRITE (ITIT,103) ITS 103 FORMAT(I3) C IF(ITS .GE. 320) THEN C CALL PWRY(SX(K-3),43.0,ITIT,3,10,0,1) C ELSE C CALL PWRY(-18.0,SY(K-3),ITIT,3,10,0,1) C END IF X=SX(K-3) Y=SY(K-3) IF(X.LT.-15.0) X = -17.95 IF(Y.GT.40.0) Y = 42.9 CALL WTSTR (X,Y,ITIT(1:3),10,0,0) 55 CONTINUE C * CALL GFLAS2 ISKDRW = 1 RETURN END FUNCTION ESAT(T) C ESAT(MILLIBARS),T(KELVIN) PARAMETER (ABZ=273.16) TC=T-ABZ ESAT=6.1078*EXP((17.2693882*TC)/(TC+237.3)) RETURN END FUNCTION OS(T, P) C OS AND T (KELVIN) , P (MILLIBARS ) OS = T*((1000./P)**.286)/(EXP(-2.6518986*W(T,P)/T)) RETURN END FUNCTION TMR(W, P) C TMR(KELVIN),W(GRAMS WATER VAPOR/KILOGRAM DRY AIR),P(MILLIBAR) X = ALOG10 (W * P / (622.+ W)) TMR = 10.**(.0498646455*X+2.4082965)-7.07475+38.9114* & ((10.**( .0915*X ) - 1.2035 )**2 ) RETURN END FUNCTION TSA(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 1 I = 1,12 D = D/2. C IF THE TEMPERATURE DIFFERENCE,X, IS SMALL,EXIT THIS LOOP X=A*EXP(-2.6518986*W(TQ,P)/TQ)-TQ*((1000./P)**.286) IF(ABS(X).LT.0.01)GO TO 2 TQ = TQ + SIGN(D,X) 1 CONTINUE 2 TSA=TQ RETURN END FUNCTION W(T, P) C W(GRAMS WATER VAPOR/KILOGRAM DRY AIR ), P(MILLIBAR ) IF (T .GE. 999.) GO TO 10 X = ESAT(T) W = 621.97 * X / (P - X) RETURN 10 W = 0.0 RETURN END SUBROUTINE WNDBARB (XBASE,YBASE,U,V) C*********************************************************************** C WNDBARB - FOR THE GRAPH PORTION OF GRIN C THIS ROUTINE DRAWS A SINGLE WIND BARB PER CALL. C C ON INPUT - FOUR VARIABLES COME IN. XBASE CONTAINS THE HORIZONTAL COOR C DINATE OF THE BASE OF THE WIND BARB. YBASE CONTAINS THE C VERTICAL COORDINATE OF THE BASE OF THE WIND BARB. U CONTAI C THE EAST-WEST WIND COMPONENT IN METERS PER SECOND. V CONTA C THE NORTH-SOUTH WIND COMPONENT IN METERS PER SECOND. C C ON OUTPUT - ONE BARB HAS BEEN DRAWN TO UNIT NUMBER 2, WHICH CORRESPOND C TO GMETA.CGM, THE GKS OUTPUT META CODE FILE. C C ASSUMPTIONS - THIS ROUTINE ASSUMES THAT GKS HAS BEEN OPENED AND A WORK C STATION HAS BEEN SET. C C REVISED BY - JEREMY ASBILL ON MAY 3, 1990. WNDBARB.20 C*********************************************************************** WNDBARB.21 WNDBARB.22 C INPUT VARIABLE DECLARATIONS ... WNDBARB.23 WNDBARB.24 REAL XBASE,YBASE, WNDBARB.25 * U,V WNDBARB.26 WNDBARB.27 C LOCAL PARAMETER ... C SC SPECIFIES IN THE NORMALIZED (FRACTIONAL) GRAPHICS COORDINATE C SYSTEM HOW LONG THE BARB SHAFT IS TO BE C COORDINATE SYSTEMS ARE EXPLAINED IN NCAR GRAPHICS USER'S GUIDE VERS C 2.00 ON PAGE 46. PARAMETER (SC = 0.05493) C LOCAL VARIABLE DECLARATIONS ... INTEGER LLSV ! SAVE VARIABLE, SCALING FOR SET LOGICAL DONE ! T => SUBROUTINE ENDS, F => LOOP A REAL WINDVCT, ! WIND VECTOR MAGNITUDE * FLSV,FRSV,FBSV,FTSV, ! SAVE VARIABLES, FRACTIONAL COORDI * ULSV,URSV,UBSV,UTSV, ! SAVE VARIABLES, INCOMING USER COO * NEWXBASE, ! FRACTIONAL X COORD. FOR BARB BASE * NEWYBASE, ! FRACTIONAL Y COORD. FOR BARB BASE * XCOMP, ! X COMPONENT OF GRAPHICAL VECTOR ( * YCOMP, ! Y COMPONENT OF GRAPHICAL VECTOR ( * PK, ! PLACE KEEPER * FETHLENX, ! X COMPONENT OF GRAPHICAL VECT. (F * FETHLENY ! Y COMPONENT OF GRAPHICAL VECT. (F C LOCAL ARRAY DECLARATIONS ... INTEGER IJUNK(5) ! CALCULATION ARRAY FOR SFSGFA REAL POINTX(3), ! USED TO SPECIFY POINTS TO DRAW BE * POINTY(3), ! USED TO SPECIFY POINTS TO DRAW BE * JUNK(7) ! CALCULATION ARRAY FOR SFSGFA C***************************** SUBROUTINE BEGIN ************************ C INITIALIZE LOOP, BOOLEAN INDICATOR DONE = .FALSE. C CALCULATE THE WIND VECTOR MAGNITUDE IN KNOTS IF ((U .EQ. 0) .AND. (V .EQ. 0)) THEN WINDVCT = 1.0 ELSE WINDVCT = SQRT(U**2 + V**2) * 1.94 END IF C SAVE INCOMING USER COORDINATES AND CHANGE BACK TO NORMALIZED COORDINA C DOCUMENTATION FOR SET AND GETSET CAN BE FOUND IN NCAR GRAPHICS USER'S C GUIDE VERSION 2.00 ON PAGES 49 (GETSET) AND 53 (SET). CALL GETSET (FLSV,FRSV,FBSV,FTSV,ULSV,URSV,UBSV,UTSV,LLSV) CALL SET (FLSV,FRSV,FBSV,FTSV, 0.0, 1.0, 0.0, 1.0, 1) C DETERMINE WHERE THE BASE OF THE BARB IS IN THE NORMALIZED COORDINATES NEWXBASE = (XBASE - ULSV)/(URSV - ULSV) NEWYBASE = (YBASE - UBSV)/(UTSV - UBSV) C CALCULATE THE X DISTANCE AND Y DISTANCE FROM THE BASE OF THE BARB THA C DEFINES THE BARBS TIP (NORMALIZED COORD'S) XCOMP = -SC * U * 1.94/WINDVCT YCOMP = -SC * V * 1.94/WINDVCT C DETERMINE THE ACTUAL LOCATION IN NORMALIZED COORDINATES OF THE BARB'S POINTX(1) = NEWXBASE + XCOMP POINTY(1) = NEWYBASE + YCOMP C DRAW THE BARB SHAFT, DOCUMENTATION FOR THE LINE SUBROUTINE CAN BE FOU C IN NCAR GRAPHICS USER'S GUIDE VERSION 2.00 ON PAGE 50 CALL LINE (NEWXBASE,NEWYBASE,POINTX(1),POINTY(1)) C DETERMINE THE FEATHER LENGTH FETHLENX = 0.3 * YCOMP FETHLENY = -0.3 * XCOMP C SET THE PLACE KEEPER AND BOOST THE WIND MAGNITUDE PK = 0.9 WINDVCT = WINDVCT + 2.5 C BEGIN MAKING FEATHERS 10 CONTINUE C DRAW A FLAG FOR EVERY 50 KNOTS WIND MAGNITUDE IF (WINDVCT .GE. 50.0) THEN C DETERMINE THE POSITION OF THE FLAG TIP, POINT_(2) C AND DETERMINE POSITION WHERE FLAG BOTTOM MEETS THE SHAFT, POINT_( POINTX(2) = POINTX(1) + FETHLENX + 0.0005 POINTY(2) = POINTY(1) + FETHLENY + 0.0005 POINTX(3) = PK * XCOMP + NEWXBASE POINTY(3) = PK * YCOMP + NEWYBASE C DRAW FLAG CALL LINE (POINTX(1),POINTY(1),POINTX(2),POINTY(2)) CALL LINE (POINTX(3),POINTY(3),POINTX(2),POINTY(2)) C FILL IN FLAG, DOCUMENTATION FOR SFSGFA CAN BE FOUND IN NCAR C GRAPHICS GUIDE TO NEW UTILITIES VERSION 3.00 ON PAGE 4-8 CALL SFSETR ('SP',0.000001) CALL SFSGFA (POINTX,POINTY,3,JUNK,5,IJUNK,7,2) C REMOVE 50 KNOTS FROM WIND MAGNITUDE (ALREADY DRAWN IN) WINDVCT = WINDVCT - 50.0 C DETERMINE NEW BEGIN POINT FOR NEXT FLAG OR FEATHER PK = PK - 0.05 POINTX(1) = PK * XCOMP + NEWXBASE POINTY(1) = PK * YCOMP + NEWYBASE PK = PK - 0.1 C DRAW A FULL FEATHER FOR WIND MAGNITUDE OF EVERY 10 KNOTS ELSE IF (WINDVCT .GE. 10.0) THEN C CALCULATE POSITION OF FEATHER END POINTX(2) = POINTX(1) + FETHLENX + 0.0005 POINTY(2) = POINTY(1) + FETHLENY + 0.0005 C DRAW FEATHER CALL LINE (POINTX(1),POINTY(1),POINTX(2),POINTY(2)) C REMOVE 10 KNOTS FROM WIND MAGNITUDE (ALREADY DRAWN IN) WINDVCT = WINDVCT - 10.0 C DETERMINE NEW START POINT FOR NEXT FEATHER OR FLAG POINTX(1) = PK * XCOMP + NEWXBASE POINTY(1) = PK * YCOMP + NEWYBASE PK = PK - 0.1 C DRAW A HALF FEATHER FOR EVERY 5 KNOTS OF WIND MAGNITUDE ELSE IF (WINDVCT .GE. 5.0) THEN C CALCULATE POSITION OF TIP OF HALF FEATHER POINTX(2) = POINTX(1) + 0.5 * FETHLENX + 0.0005 POINTY(2) = POINTY(1) + 0.5 * FETHLENY + 0.0005 C DRAW IN FEATHER CALL LINE (POINTX(1),POINTY(1),POINTX(2),POINTY(2)) C TELL LOOP TO QUIT DONE = .TRUE. ELSE DONE = .TRUE. END IF C IF THERE IS STILL MORE WIND MAGNITUDE (>= 5 KNOTS) LOOP AGAIN IF (.NOT. DONE) GOTO 10 C RESET USER COORDINATES TO THE INCOMING VALUES CALL SET (FLSV,FRSV,FBSV,FTSV,ULSV,URSV,UBSV,UTSV,LLSV) C****************************** SUBROUTINE END ************************* RETURN END