PROGRAM figpedav C Also need: bplot.f C Compile: ncargf77 -o figpedav.run figpedav.f bplot.f C Run as: figpedav.run C Program to plot Figure 4-like figure, but with Pedersen-averaged winds C ************************************************************************** c PARAMETER (NUMK=200,NUMAZI=1,NUMNP=300) PARAMETER (NUMK=200,NUMAZI=1,NUMNP=1) C ************************************************************************** DIMENSION DAT16(NUMK,NUMAZI,NUMNP),RLAT16(NUMK,NUMAZI,NUMNP), | RMLT16(NUMK,NUMAZI,NUMNP),CD16(NUMK,NUMAZI,NUMNP), | SD16(NUMK,NUMAZI,NUMNP),N16(NUMAZI,NUMNP), | DAT2P(NUMK,NUMAZI,NUMNP), | ERDAT16(NUMK,NUMAZI,NUMNP),ERDAT26(NUMK,NUMAZI,NUMNP) CHARACTER*3 NAM3 CHARACTER*13 INFILE CHARACTER*50 CHAR50 CHARACTER*14 DATLBL CHARACTER*8 APLABEL C BLATPL is the bottom latitude plotted. Is usually 50 degrees. COMMON/BOTLAT/BLATPL c CHARACTER*8 SPMESS,LINE4,TDATE CHARACTER*8 SPMESS,TDATE CHARACTER*26 LINE4 CHARACTER*30 MESS COMMON /DATTIM/ MESS,SPMESS,LINE4 COMMON /HEM/ IHSOLN CHARACTER*26 LABL3,LABL4 COMMON/FXFY/TWOPI,XLONMX,ROT C Suitable for 244 basis funcs PARAMETER (ITHMX=54,LONMX=36,ITHPLT=24) COMMON/ERR/ERAR(0:LONMX,0:ITHPLT),ERRMX,ISOLF,IERRD COMMON /LBLS/ LABL3,LABL4,TDATE common/edges/xvpl,xvpr,yvpb,yvpt common/vecscale/ivecscale,ilabl3,iline4,ipass INTEGER TDMY(3) CHARACTER*1 FLGB CHARACTER*24 DESC1 CHARACTER*16 DESC2 CHARACTER*4 SCWRD1 c DATA DESC1/' Average Wind Velocity'/ DATA DESC1/' '/ C DATA DESC2/'(rotated -90deg)'/ C DATA IROT/-90/ DATA IROT/0/ C iscale1 = number of m/s in a vector of length 10 degrees latitude. DATA ISCALE1/200/ DATA SCWRD1/' m/s'/ DATA IERRD/0/ DATA BLATPL/50./ DATA SPMESS/' '/ DATA LINE4/' '/ DATA IHSOLN/1/ DATA LABL3/' '/ DATA LABL4/' '/ CALL IDATE (TDMY) WRITE (TDATE,"(2(i2.2,'/'),i2.2)") TDMY(3)-1900,TDMY(2),TDMY(1) WRITE (6,"(1x,'Running fig3 on ',a8)") TDATE TDATE = ' ' PII = 4.*ATAN(1.) RAD = 180./PII TWOPI=2.*PII XLONMX=24 ROT = 0. QCIRCLE = PII/2. * (90.-BLATPL)/90. open(11,file='wnpedavmy',status='old') open(12,file='wnpedavy',status='old') open(13,file='wnpedavmz',status='old') open(14,file='wnpedavz',status='old') open(15,file='y-3.2z0',status='old') open(16,file='y+3.2z0',status='old') open(17,file='y0z-2.0',status='old') open(18,file='y0z+2.0',status='old') CALL OPNGKS CALL GSCLIP (0) CALL GSLWSC (4.0) Do 2000 iun=11,18 ipass = iun - 10 iyplace = 1 LABL3(9:26) = ' ' c if (iyplace.eq.1.and.ipass.eq.1) c 1 LABL3(9:26) = ' 140 km ' if (iun.ge.15) then ipass = iun - 14 iyplace = 2 endif xvpl = .005 + (ipass-1)*.25 xvpr = .245 + (ipass-1)*.25 yvpb = .33 - (iyplace-2)*.3 yvpt = .57 - (iyplace-2)*.3 ivecscale = 0 if (ipass.eq.4) ivecscale=1 iline4 = 0 if (iyplace.eq.1) iline4 = 1 LINE4 = ' ' if (iyplace.eq.1) then if (ipass.eq.1) LINE4 = ' B:B:y:N: = -3.2 :L:n:U:T' if (ipass.eq.2) LINE4 = ' B:B:y:N: = 3.2 :L:n:U:T' if (ipass.eq.3) LINE4 = ' B:B:z:N: = -2.0 :L:n:U:T' if (ipass.eq.4) LINE4 = ' B:B:z:N: = 2.0 :L:n:U:T' endif c X = QCIRCLE c CALL SET(xvpl,xvpr,yvpb,yvpt,-X,X,-X,X,1) DO 50 I=1,NUMAZI DO 50 J=1,NUMNP 50 N16(I,J) = 0 C iflsc = 1 means plot each scan separately between uts1 and uts2 C iflsc = 2 means plot all the points together C ilos = 1,2 for l-o-s or vector data C nhd = number of lines in the header IFLSC = 1 ILOS = 2 NHD = 3 NAM3 = 'win' NUMRNG = 1 NUMAZ = 1 SCANSC = 3000 NPL = 1 NUM = 145 UTSPAC = SCANSC*FLOAT(NPL) UTPL = UTSPAC/FLOAT(NPL) C Read input file READ (iun,"(6x,a14,2x,a8)") DATLBL,APLABEL c WRITE (LABL3,"(a14,12x)") DATLBL c WRITE (DESC2,"(8x,a8)") APLABEL c WRITE (LINE4,"(a8)") APLABEL IF (APLABEL.EQ.'ap25-ap5') DESC1 = 'Difference Wind Velocity' READ (iun,"(a50)") CHAR50 X = 8. CALL SET(xvpl,xvpr,yvpb,yvpt,-X,X,-X,X,1) CALL GSELNT(1) CALL FULCIR(8.,LABL3,LABL4) X = QCIRCLE CALL SET(xvpl,xvpr,yvpb,yvpt,-X,X,-X,X,1) IDAP = -1 KLOCP = 100 NSCN = 0 NP = 0 AVIDM = 0. KLOC = 0 DO 1000 NRD=1,NUM C Read vector file READ (iun,"(i4,f6.2,f6.1,f7.1,f5.0,f6.1,f7.1,f6.2,9f7.1)", 1 END=1001) | IDAYN,UTHR,GLAT,GLON,ALT,ALAT,ALON,RMLT,UNEW,ERRUN,VNNS,ERRVN, | VIEW,ERREW,VINS,ERRNS,EMISS IF (EMISS .LT. 0.001) GO TO 1000 DAT = VINS ERDAT = ERRNS DAT2 = VIEW ERDAT2 = ERREW IDA = 0 KLOC = KLOC + 1 ANG1N = 0. UTSC = 0. FLGB = ' ' KAZ = 0 NP = 1 N16(KAZ+1,NP) = N16(KAZ+1,NP) + 1 K = N16(KAZ+1,NP) N100 = NRD/1000 NRD100 = N100*1000 IF (NRD100 .EQ. NRD) | WRITE (6,"(1x,'nrd n16(1)=',17i5)") NRD,(N16(I,1),I=1,NUMAZ) IF (NP .GT. NUMNP .OR. K .GT. NUMK .OR. KAZ+1 .GT. NUMAZI) THEN WRITE (6,"(1x,'Increase parameters: nrd np numnp k numk' | ' kaz+1 numazi=',7i5)") NRD,NP,NUMNP,K,NUMK,KAZ+1,NUMAZI STOP ENDIF RLAT16(K,KAZ+1,NP) = ALAT RMLT16(K,KAZ+1,NP) = RMLT CD16(K,KAZ+1,NP) = COS(ANG1N/RAD) SD16(K,KAZ+1,NP) = SIN(ANG1N/RAD) DAT16(K,KAZ+1,NP) = DAT ERDAT16(K,KAZ+1,NP) = ERDAT IF (ILOS .EQ. 1) GO TO 900 DAT2P(K,KAZ+1,NP) = DAT2 ERDAT26(K,KAZ+1,NP) = ERDAT2 900 CONTINUE 950 IDAP = IDA 1000 CONTINUE 1001 CONTINUE IALT = ALT c WRITE(LINE4,"(i3,' km',2x)") IALT AVIDM = AVIDM/FLOAT(NRD) WRITE (6,"(1X,'NRD AVIDM =',I3,F7.2)") NRD,AVIDM C Plot at every azimuth write(6,*),NPL,NUMAZ DO 9020 NP=1,NPL DO 9020 NAZ=1,NUMAZ IF (N16(NAZ,NP) .LE. 0) GO TO 9020 c MESS = 'WINDII 557.7 nm ' MESS = ' ' IHSOLN = -1 IF (RLAT16(1,NAZ,NP) .GT. 0.) IHSOLN = 1 WRITE (6,"(1x,'np npts =',i3,i6,2x,a30,i2)") NP,N16(NAZ,NP), | MESS,IHSOLN CALL PLTV (2,N16(NAZ,NP),RLAT16(1,NAZ,NP),RMLT16(1,NAZ,NP), | CD16(1,NAZ,NP),SD16(1,NAZ,NP),DAT16(1,NAZ,NP),DAT2P(1,NAZ,NP), | DAT16(1,NAZ,NP),ERDAT16(1,NAZ,NP),ERDAT26(1,NAZ,NP), | ERDAT16(1,NAZ,NP),IROT,ISCALE1,SCWRD1,ISCALE1,SCWRD1,DESC1, | DESC2) 9020 CONTINUE 2000 CONTINUE CALL GSCLIP (1) CALL CLSGKS STOP END SUBROUTINE PLTV (IFLG,NUM,RLAT,RMLT,COSD,SIND,DAT1,DAT2,DAT3, | ERDAT1,ERDAT2,ERDAT3,IROT,ISCALE1,SCWRD1,ISCALE3,SCWRD3,DESC1, | DESC2) C Need bplot.f and fulcir.f C Subroutine to plot vectors on a polar plot C iflg = 1 if only have dat1 filled (i.e., single valued) C iflg = 2 if have both dat1 and dat2 (i.e., a vector) C iflg = 3 if have both dat1 and dat2 and dat3 (i.e. deltaB) C num = number of points in rest of arrays C rlat = magnetic latitude in degrees C rmlt = magnetic local time in hours C cosd,sind = cos and sin of angle between measurement and magnetic north C dat1 = data in first array, could be toward mag north C dat2 = data in second array, could be toward mag east C dat3 = data in third array, could be upward or downward C erdat1 = error of data in first array, could be toward mag north C erdat2 = error of data in second array, could be toward mag east C erdat3 = error of data in third array, could be upward or downward C irot = rotation of vectors for plots. Only values are: -90,0,90 C irot=-90, if dat1,2 are E fields plotted in dir of ion drift C irot=90, if dat1,2 are deltaBs rotated in dir of horiz current C irot=0, if dat1,2 are ion or neutral velocities. C irot=-90 or 90 if dat1 is a conductance or other scalar to be C plotted perpendicular to the satellite track. C iscale1 = scale for vector (dat1,2) or line (dat1) (i.e., 50 if mV/m) C scwrd1 = 4 character descriptor for scale1 ('mV/m','m/s ','mho ',etc) C iscale3 = scale for 3rd component (dat3) (i.e., 200 nT Z component) C scwrd3 = 4 character descriptor for scale3 ('nT ', for Z component) C desc1 = 24 character description in upper right hand side of plot C desc2 = 16 character description in 2cnd row of upper rhs of plot SAVE DIMENSION RLAT(NUM),RMLT(NUM),COSD(NUM),SIND(NUM),DAT1(NUM), | DAT2(NUM),DAT3(NUM),ERDAT1(NUM),ERDAT2(NUM),ERDAT3(NUM) C BLATPL is the bottom latitude plotted. Is usually 50 degrees. COMMON/BOTLAT/BLATPL c CHARACTER*8 SPMESS,LINE4,TDATE CHARACTER*8 SPMESS,TDATE CHARACTER*26 LINE4 CHARACTER*30 MESS COMMON /DATTIM/ MESS,SPMESS,LINE4 COMMON /HEM/ IHSOLN CHARACTER*26 LABL3,LABL4 CHARACTER*24 DESC1 CHARACTER*16 DESC2 CHARACTER*4 SCWRD1,SCWRD3 CHARACTER*8 SCLAB1,SCLAB3 COMMON /LBLS/ LABL3,LABL4,TDATE common/edges/xvpl,xvpr,yvpb,yvpt common/vecscale/ivecscale,ilabl3,iline4,ipass DATA IFRST/0/ IF (IROT .EQ. -90 .OR. IROT .EQ. 0 .OR. IROT .EQ. 90) GO TO 100 WRITE (6,"(1x,'irot .ne. -90, 0, or 90 -- stop. irot=',i4)") | IROT STOP 100 IFRST = IFRST + 1 IF (IFRST .EQ. 1) THEN c SZLBL = 0.016 c SZLBL = 0.02 SZLBL = 0.007 SZDAT = 0.010 SZSML = 0.006 XLN1 = 0.923 c YLNE = 0.96 YLNE = 0.86 PII=4.*ATAN(1.) RAD = 180./PII QCIRCLE = PII/2. * (90.-BLATPL)/90. HTR = PII/12. ENDIF C SCL=.1/rad implies 50 mV/m in 10 degrees on plot down to 50 degrees C SCL = 5/rad implies 1 in 10 degrees SCALE1 = ISCALE1 SCALE3 = ISCALE3 SCL = (500./SCALE1)*(0.01/RAD) SCLZ = (500./SCALE3)*(0.01/RAD) WRITE (SCLAB1,"(i4,a4)") ISCALE1,SCWRD1 Cc write (sclab1,"(i3,1x,a4)") iscale1,scwrd1 WRITE (SCLAB3,"(i3,1x,a4)") ISCALE3,SCWRD3 c CALL SET (0.1,0.9,0.1,0.9,-QCIRCLE,QCIRCLE,-QCIRCLE,QCIRCLE,1) X = QCIRCLE CALL SET(xvpl,xvpr,yvpb,yvpt,-X,X,-X,X,1) IF (NUM .LT. 1) RETURN DO 1000 N=1,NUM DASH = 0. COLAT = 90. - ABS(RLAT(N)) R = COLAT/RAD C Do not plot if beyond outside latitude IF (R .GT. QCIRCLE) GO TO 1000 IF (IFLG .EQ. 1) THEN C Plot single value DTRUN = DAT1(N)*COSD(N) DTRUE = DAT1(N)*SIND(N) ERLG = 4.*(ERDAT1(N)**2)/(DAT1(N)**2+3.) ELSE C Plot vector value (often dat1=magN and dat2=magE) DTRUN = DAT1(N)*COSD(N) - DAT2(N)*SIND(N) DTRUE = DAT2(N)*COSD(N) + DAT1(N)*SIND(N) ERLG =4.*(ERDAT1(N)**2+ERDAT2(N)**2)/(DAT1(N)**2+DAT2(N)**2+3.) ENDIF C Note that erlg was determined for E fields in mV/m, so could be C different for different quantities. IF (ERLG .GT. 1.) DASH = -1. X = R*SIN(RMLT(N)*HTR) Y = -R*COS(RMLT(N)*HTR) C Set values according to rotation BZ = 0. IF (IFLG .EQ. 3) BZ = DAT3(N) IF (IROT .EQ. -90) THEN BN = -DTRUN BE = -DTRUE ENDIF IF (IROT .EQ. 0) THEN BN = DTRUE BE = DTRUN ENDIF IF (IROT .EQ. 90) THEN BN = DTRUN BE = DTRUE ENDIF CALL GSLWSC (2.0) CALL BPLOT (X,Y,R,SCL,SCLZ,BN,BE,BZ,1.,DASH) CALL GSLWSC (1.0) 1000 CONTINUE 4506 CALL GSELNT (0) c CALL PLCHHQ (YLNE,0.17,SCLAB1,SZLBL,0.,1.) c CALL PLCHHQ (.94,0.17,SCLAB1,SZLBL,0.,1.) c CALL PLCHHQ (.89,0.15,SCLAB1,SZLBL,0.,1.) CALL PLCHHQ (.995,0.34,SCLAB1,SZLBL,0.,1.) CALL PLCHHQ (.995,0.64,SCLAB1,SZLBL,0.,1.) c CALL PLCHHQ (.5,0.92,'Neutral Wind at 140 km',.012,0.,0.) CALL PLCHHQ (.5,0.92,'Pedersen-Weighted Neutral Wind',.012,0.,0.) CALL PLCHHQ (.5,0.595,'Convection Velocity',.012,0.,0.) CALL GSELNT(1) c X = 0.925 * QCIRCLE X = 0.8 * QCIRCLE Y = -X R = SQRT(X*X + Y*Y) C 97/9/13 ADR C 353.6 = 250.*sqrt(2.). Next two lines could be simplified; same result C is obtained as long as SCL*HSCFAC=50./(RAD*sqrt(2.)) and C SCLZ*HSCFAC=50./(RAD*sqrt(2.)), i.e, scale vector = 10 deg. latitude. HSCFAC = 353.6 * (SCALE1/500.) CALL GSLWSC (2.0) if (ivecscale.eq.1) 1 CALL BPLOT(X,Y,R,SCL,SCLZ,HSCFAC,HSCFAC,0.,1.,0.) CALL GSLWSC (1.0) IF (IFLG .EQ. 3) THEN IF (ISCALE1 .NE. ISCALE3) THEN CALL GSELNT(0) CALL PLCHHQ (YLNE,.08,SCLAB3,SZLBL,0.,1.) CALL GSELNT(1) ENDIF X = 0.837 * QCIRCLE Y = -1.037 * QCIRCLE R = SQRT(X*X + Y*Y) ZFAC = 353.6 * (SCALE3/500.) CALL BPLOT (X,Y,R,SCL,SCLZ,0.,0.,ZFAC,0.,1.) X = 1.012 * QCIRCLE R = SQRT(X*X + Y*Y) CALL BPLOT (X,Y,R,SCL,SCLZ,0.,0.,-ZFAC,0.,1.) ENDIF CALL GSELNT (0) c CALL PLCHHQ(YLNE,XLN1,DESC1,SZLBL,0.0,1.0) c CALL PLCHHQ(YLNE,XLN1-0.04,DESC2,SZLBL,0.0,1.0) c CALL PLCHHQ(.03,.03,TDATE,SZDAT,0.,-1.) CALL GSELNT(1) C c CALL FULCIR(QCIRCLE,LABL3,LABL4) c CALL FRAME 4510 CONTINUE RETURN END SUBROUTINE FULCIR (Q,LABL3,LABL4) C C BLATPL is the bottom latitude plotted in MGPLT and FULCIR to show C the ground magnetometer data. Is usually 50 degrees. COMMON/BOTLAT/BLATPL C include "Headers/dattim.h" C include "Headers/hem.h" c CHARACTER*8 SPMESS,LINE4,LINE5 CHARACTER*8 SPMESS,TDATE CHARACTER*26 LINE4 CHARACTER*30 MESS COMMON /DATTIM/ MESS,SPMESS,LINE4 COMMON /HEM/ IHSOLN common/vecscale/ivecscale,ilabl3,iline4,ipass CHARACTER*14 LATBND CHARACTER*26 LABL3,LABL4 DIMENSION XMLT(4),YMLT(4) DIMENSION X(145),Y(145),X1(145),Y1(145) SAVE IFIRST,X,Y,SZLBL,XLN1 DATA IFIRST,PII/1,3.1415926535898/ C J10 = (90.1-BLATPL)/10 IF (IFIRST.NE.1) GO TO 100 IFIRST = 0 SZLBL = 0.016 SZLBL = 0.023 SZLBL = 0.020 SZLBL = 0.006 SZLBL = 0.007 SZLBL = 0.008 XLN1 = 0.923 DO 10 I=1,145 T = (I+52)*PII/72. X(I) = SIN(T) Y(I) = -COS(T) 10 CONTINUE 100 CALL GSELNT(0) CALL GSELNT(1) CALL GSLWSC (1.0) C Change default function character from colon to semi-colon CALL PCSETC ('FC - function code delimiter', ';') C MESS has a colon in it! c CALL PLCHHQ (0.04, XLN1, MESS, SZLBL, 0.0, -1.0) C Reset default function character from semi-colon to colon CALL PCSETC ('FC - function code delimiter', ':') C TEMP!! c CALL PLCHHQ (0.22, 0.07, SPMESS, SZLBL, 0., -1.0) c CALL PLCHHQ (0.04, XLN1-0.04, LABL3, SZLBL, 0., -1.0) c CALL PLCHHQ (0.04, XLN1-0.08, LABL4, SZLBL, 0., -1.0) c CALL PLCHHQ (0.04, XLN1-0.12, LINE4, SZLBL, 0., -1.0) c CALL PLCHHQ (0.04, XLN1-0.16, LINE5, SZLBL, 0., -1.0) CALL GSELNT(1) DO 11 J=1,J10-1 R = Q * FLOAT(J) / FLOAT(J10) DO 5 I=1,145 X1(I) = R*X(I) 5 Y1(I) = R*Y(I) C 89/1/17 MAKE COORDINATE GRID LIGHT CALL GSLWSC (1.0) Cc CALL SETUSV ('IN',7000) CALL CURVE(X1,Y1,145) 11 CONTINUE DO 6 I=1,145 X1(I) = Q*X(I) 6 Y1(I) = Q*Y(I) CALL GSLWSC (4.0) CALL CURVE(X1,Y1,141) CALL GSLWSC (1.0) C LINE does not recognize the command GSLWSC, but CURVE does. C CALL LINE (-Q,0., Q,0.) C CALL LINE (0.,-Q, 0.,Q) QQ = .88*Q XMLT(1) = -Q*.85 XMLT(2) = Q*.85 YMLT(1) = 0. YMLT(2) = 0. XMLT(3) = 0. XMLT(4) = 0. YMLT(3) = -QQ YMLT(4) = QQ CALL CURVE (XMLT,YMLT,2) CALL CURVE (XMLT(3),YMLT(3),2) c CALL GSLWSC (4.0) Cc CALL SETUSV ('IN',9000) c QQ=Q*1.06 QQ=Q*.94 QQM=-QQ c Q06=Q*1.07 Q06=Q*.93 Q18=-Q06 CALL GSELNT(1) CALL PLCHHQ (-.75*Q, .935*Q, LABL3, SZLBL, 0., 0.0) c CALL PLCHHQ ( .75*Q, .935*Q, LINE4, SZLBL, 0., 0.0) c CALL PLCHHQ ( 0.,1.2*Q, LINE4, .012, 0., 0.0) CALL PLCHHQ ( 0.,1.15*Q, LINE4, .010, 0., 0.0) c xofst = -.57*Q c if (ipass.eq.2.or.ipass.eq.4) xofst = -.45*Q xofst = -.475*Q if (ipass.eq.2.or.ipass.eq.4) xofst = -.375*Q if (iline4.eq.1) c 1 CALL PLCHHQ ( xofst,1.3*Q, '-', .012, 0., 0.0) 1 CALL PLCHHQ ( xofst,1.25*Q, '-', .012, 0., 0.0) CALL PLCHHQ (0., QQM, '00', SZLBL, 0., 0.0) CALL PLCHHQ (Q06, 0., '06', SZLBL, 0., 0.0) CALL PLCHHQ (0., QQ, '12', SZLBL, 0., 0.0) CALL PLCHHQ (Q18, 0., '18', SZLBL, 0., 0.0) IFIFTY = BLATPL*IHSOLN WRITE (LATBND,"(I3,':KRL:.:PRU:')") IFIFTY CALL PLCHHQ (X1(143),Y1(143),LATBND,SZLBL, 0., 0.) C CALL PLCHHQ (X1(143),Y1(143),'50:KRL:.:PRU:',SZLBL, 0., 0.) CALL GSELNT(0) RETURN END