c c------------------------------------------------------------------ c Begin file /home/sting/foster/libplt/velvct.f c------------------------------------------------------------------ c SUBROUTINE VELVCT (U,LU,V,LV,M,N,FLO,HI,NSET,LENGTH,ISPV,SPV, + MM,INC,incyv) c call velvct(udat,imx,vdat,imx,imx,nlats, c + 0.,vgcmx,1,length,ispvalue,spval,mproj,inc(j0,mj),incyv) C C --------------------------------------------------------------------- C C DECLARATIONS - C COMMON /VEC1/ ASH ,EXT ,ICTRFG ,ILAB , + IOFFD ,IOFFM ,ISX ,ISY , + RMN ,RMX ,SIDE ,SIZE , + XLT ,YBT ,ZMN ,ZMX C COMMON /VEC2/ BIG ,INCX ,INCY DIMENSION INC(*) C C FORCE THE BLOCK DATA ROUTINE, WHICH SETS DEFAULT VARIABLES, TO LOAD. C EXTERNAL VELDAT C C ARGUMENT DIMENSIONS. C DIMENSION U(LU,N) ,V(LV,N) ,SPV(2) CHARACTER*10 LABEL character*2 chproj CHARACTER*12 LABV REAL WIND(4), VIEW(4), IAR(4) C C --------------------------------------------------------------------- C SAVE C FX(XX,YY) = XX C FY(XX,YY) = YY DIST(XX,YY) = SQRT(XX*XX+YY*YY) C MXF(XX,YY,UU,VV,SFXX,SFYY,MXX,MYY) = MXX+IFIX(SFXX*UU) C MYF(XX,YY,UU,VV,SFXX,SFYY,MXX,MYY) = MYY+IFIX(SFYY*VV) SCALEX(MM,NN,INCXX,INCYY,HAA,XX1,XX2,YY1,YY2,XX3,XX4,YY3,YY4, 1 LENN) = LENN/HAA SCALEY(MM,NN,INCXX,INCYY,HAA,XX1,XX2,YY1,YY2,XX3,XX4,YY3,YY4, 1 LENN) = SCALEX(MM,NN,INCXX,INCYY,HAA,XX1,XX2,YY1,YY2,XX3, 2 XX4,YY3,YY4,LENN) VLAB(UU,VV,II,JJ) = DIST(UU,VV) C C --------------------------------------------------------------------- C C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR. C CALL Q8QST4 ('NSSL','VELVCT','VELVCT','VERSION 6') C C INITIALIZE AND TRANSFER SOME ARGUMENTS TO LOCAL VARIABLES. C CALL MAPGTC('PR',chproj) incy = incyv LEN = LENGTH BIG = -R1MACH(2) MX = LU MY = LV NX = M NY = N GL = FLO HA = HI ISP = ISPV NC = 0 C C COMPUTE CONSTANTS BASED ON THE ADDRESSABILITY OF THE PLOTTER. C CALL GETUSV('XF',ISX) CALL GETUSV('YF',ISY) ISX = 2**(15-ISX) ISY = 2**(15-ISY) LEN = LENGTH*ISX C C SET UP THE SCALING OF THE PLOT. C CALL GQCNTN(IERR,IOLDNT) CALL GQNT(IOLDNT,IERR,WIND,VIEW) X1 = VIEW(1) X2 = VIEW(2) Y1 = VIEW(3) Y2 = VIEW(4) X3 = WIND(1) X4 = WIND(2) Y3 = WIND(3) Y4 = WIND(4) CALL GETUSV('LS',IOLLS) C C SAVE NORMALIZATION TRANSFORMATION 1 C CALL GQNT(1,IERR,WIND,VIEW) C IF (NSET) 101,102,106 C 101 X3 = 1. X4 = FLOAT(NX) Y3 = 1. Y4 = FLOAT(NY) GO TO 105 C 102 X1 = XLT X2 = XLT+SIDE Y1 = YBT Y2 = YBT+SIDE X3 = 1. Y3 = 1. X4 = FLOAT(NX) Y4 = FLOAT(NY) IF (AMIN1(X4,Y4)/AMAX1(X4,Y4) .LT. EXT) GO TO 105 C IF (NX-NY) 103,105,104 103 X2 = XLT+SIDE*X4/Y4 GO TO 105 104 Y2 = YBT+SIDE*Y4/X4 C 105 CALL SET(X1,X2,Y1,Y2,X3,X4,Y3,Y4,1) IF (NSET .EQ. 0) CALL PERIM (1,0,1,0) C C CALCULATE A LENGTH IF NONE PROVIDED. C 106 IF (LEN .NE. 0) GO TO 107 CALL FL2INT(FX(1.,1.),FY(1.,1.),MX,MY) CALL FL2INT(FX(FLOAT(1+INCX),FLOAT(1+INCY)), + FY(FLOAT(1+INCX),FLOAT(1+INCY)),LX,LY) LEN = SQRT((FLOAT(MX-LX)**2+FLOAT(MY-LY)**2)/2.) C C SET UP SPECIAL VALUES. C 107 IF (ISP .EQ. 0) GO TO 108 SPV1 = SPV(1) SPV2 = SPV(2) IF (ISP .EQ. 4) SPV2 = SPV(1) C C FIND THE MAXIMUM VECTOR LENGTH. C 108 IF (HA .GT. 0.) GO TO 118 C HA = BIG IF (ISP .EQ. 0) GO TO 115 C DO 114 J=1,NY,INCY DO 113 I=1,NX,INCX IF (ISP-2) 109,111,110 109 IF (U(I,J) .EQ. SPV1) GO TO 113 GO TO 112 110 IF (U(I,J) .EQ. SPV1) GO TO 113 111 IF (V(I,J) .EQ. SPV2) GO TO 113 112 HA = AMAX1(HA,DIST(U(I,J),V(I,J))) 113 CONTINUE 114 CONTINUE GO TO 126 C 115 DO 117 J=1,NY,INCY DO 116 I=1,NX,INCX HA = AMAX1(HA,DIST(U(I,J),V(I,J))) 116 CONTINUE 117 CONTINUE C C BRANCH IF NULL VECTOR SIZE. C 126 IF (HA .LE. 0.) GO TO 125 C C COMPUTE SCALE FACTORS. C 118 SFX = SCALEX(M,N,INCX,INCY,HA,X1,X2,Y1,Y2,X3,X4,Y3,Y4,LEN) SFY = SCALEY(M,N,INCX,INCY,HA,X1,X2,Y1,Y2,X3,X4,Y3,Y4,LEN) IOFFDT = IOFFD IF (GL.NE.0.0 .AND. (ABS(GL).LT.0.1 .OR. ABS(GL).GE.1.E5)) 1 IOFFDT = 1 IF (HA.NE.0.0 .AND. (ABS(HA).LT.0.1 .OR. ABS(HA).GE.1.E5)) 1 IOFFDT = 1 ASH = 1.0 IF (IOFFDT .NE. 0) 1 ASH = 10.**(3-IFIX(ALOG10(AMAX1(ABS(GL),ABS(HA)))-500.)-500) IZFLG = 0 C C COMPUTE ZMN AND ZMX, WHICH ARE USED IN DRWVEC. C ZMN = LEN*(GL/HA) ZMX = FLOAT(LEN)+.01 c c (cx,cy) = center of projection c r = radius of projection c if (chproj.eq.'SV'.or.chproj.eq.'ST') then cx = x1+0.5*(x2-x1) cy = y1+0.5*(y2-y1) r = x2-cx endif C C DRAW THE VECTORS. c c write(6,"('velvct: ny=',i2,' incy=',i2,' nx=',i2,' inc=', c + /(12i3))") ny,incy,nx,(inc(j),j=1,ny) DO 123 J=2,NY,INCY INCX=INC(J) DO 122 I=2,NX,INCX UI = U(I,J) VI = V(I,J) IF (ISP-1) 121,119,120 119 IF (UI-SPV1) 121,122,121 120 IF (VI .EQ. SPV2) GO TO 122 IF (ISP .GE. 3) GO TO 119 121 X = I Y = J fxx = fx(x,y) fyy = fy(x,y) if ((fxx.eq.1.e12.or.fyy.eq.1.e12).or. + (fxx.eq.0..and.fyy.eq.0.)) goto 122 CALL FL2INT(FXX,FYY,MX,MY) c c Warning: ui and vi cannot both be 0 going into MXF (see gcmdif) c write(6,"('velvct before mxf: x y=',2e12.4,' ui vi=', c + 2e12.4,/' sfx y=',2e12.4,' mx,y=',2i8,' i j=',2i4)") c + x,y,ui,vi,sfx,sfy,mx,my,i,j LX = MAX0(1,MXF(X,Y,UI,VI,SFX,SFY,MX,MY)) LY = MAX0(1,MYF(X,Y,UI,VI,SFX,SFY,MX,MY)) if ((chproj.eq.'SV'.or.chproj.eq.'ST').and. + dist(cmfx(lx)-cx,cmfy(ly)-cy).gt.r) goto 122 IZFLG = 1 IF (ILAB .NE. 0) CALL ENCD(VLAB(UI,VI,I,J),ASH,LABEL,NC, + IOFFDT) CALL DRWVEC (MX,MY,LX,LY,LABEL,NC) c write(6,"('velvct: i=',i2,' j=',i2,' fxx=',f6.1,' fyy=', c + f5.1,' lx,ly=',2i6,' mx,my=',2i6,' ui,vi=',2f7.2)") c + i,j,fxx,fyy,lx,ly,mx,my,ui,vi 122 CONTINUE 123 CONTINUE C IF (IZFLG .EQ. 0) GO TO 125 IF (IOFFM .NE. 0) GO TO 200 IHA = HA WRITE(LABV,"(I4,' M/S')") IHA C C TURN OFF CLIPPING SO ARROW CAN BE DRAWN C CALL GQCLIP(IER,ICLP,IAR) CALL GSCLIP(0) CALL GSPLCI(1) CALL DRWVEC(26000,700,26000+LEN,700,0,0) c CALL DRWVEC(28000,10000,28000+LEN,10000,0,0) C C RESTORE CLIPPING C CALL GSCLIP(ICLP) IX = 1+(28768+LEN/2)/ISX IY = 1+(608-(5*ISX*MAX0(256/ISX,8))/4)/ISY CALL GQCNTN(IER,ICN) CALL GSELNT(0) XC = CPUX(IX) YC = CPUY(IY) call plchhq(cpux(800),cpuy(50),labv,.01,0.,-1.) c call plchhq(cpux(850),cpuy(285),labv,.01,0.,-1.) CALL GSELNT(ICN) C C DONE. C GOTO 200 C C ZERO-FIELD ACTION. C 125 IX = 1+16384/ISX IY = 1+16384/ISY CALL GQCNTN(IER,ICN) CALL GSELNT(0) XC = CPUX(IX) YC = CPUY(IY) call plchhq(xc,yc,'ZERO FIELD',.01,0.,0.) CALL GSELNT(ICN) C C RESTORE TRANS 1 AND LOG SCALING AND ORIGINAL TRANS NUMBER C 200 CONTINUE IF (NSET .LE. 0) THEN CALL SET(VIEW(1),VIEW(2),VIEW(3),VIEW(4), - WIND(1),WIND(2),WIND(3),WIND(4),IOLLS) ENDIF CALL GSELNT(IOLDNT) RETURN END c FUNCTION FX(X,Y) C **** TRANSFORMATION FOR CONREC COMMON/TRANS/M,N,XLONMI,XLONMA,YLATMI,YLATMA,IF,IQZQZ save xans,yans IF(IQZQZ.EQ.0)GO TO 5 ASSIGN 2 TO IXY 4 CONTINUE IFF=0 IF(IF.EQ.1)GO TO 1 XLON=XLONMI+(XLONMA-XLONMI)*(X-1.)/(FLOAT(M)-1.) IF(XLON.GE.180.)XLON=XLON-360. YLAT=YLATMI+(YLATMA-YLATMI)*(Y-1.)/(FLOAT(N)-1.) if (iqzqz.eq.1) then CALL MAPTRN(YLAT,XLON,XANS,YANS) else xans = xlon yans = ylat endif IFF=1 1 CONTINUE IF=IFF GO TO IXY(2,3) 2 CONTINUE FX=XANS RETURN C ENTRY FY(X,Y) IF(IQZQZ.EQ.0)GO TO 6 ASSIGN 3 TO IXY GO TO 4 3 CONTINUE FY=YANS RETURN 5 CONTINUE FX=X RETURN 6 CONTINUE FY=Y RETURN END