c c PLOTTING PROGRAM FOR HORIZONTAL PROPERTY CONTOURS C AND VELOCITY VECTORS C PROGRAM MAIN c============================================================ c c nread = total read records c c zlev = the interpolated z-level to be plotted. c If Zlev < -5000., the bottom sigma c level is plotted. c c The above variables are redefined in runpl through 'params' c This code created the two .ps files in the main directory. C Without modification NCAR graphics are required. c============================================================ c PARAMETER (IM=65,JM=49,KB=11,IMM1=IM-1,JMM1=JM-1,KBM1=KB-1) PARAMETER (IMM2=IM-2,JMM2=JM-2,KBM2=KB-2) PARAMETER (NXTIC=8,NYTIC=5) c----------------- 1-D ARRAYS ------------------------------- COMMON/BLK1D/ 1 Z(KB),ZZ(KB),DZ(KB),DZZ(KB),X(IM),Y(JM) C---------------- 2-D ARRAYS -------------------------------------- COMMON/BLK2D/H(IM,JM),DX(IM,JM),DY(IM,JM), 1 ART(IM,JM),ARU(IM,JM),ARV(IM,JM), 2 DUM(IM,JM),DVM(IM,JM),FSM(IM,JM),COR(IM,JM), 3 TPS(IM,JM), 4 UA(IM,JM), 5 VA(IM,JM),EL(IM,JM),TLEV(IM,JM),PSI(IM,JM), 6 WX1(IM,JM,2),WX2(IM,JM,2), 7 X2D(IM,JM),Y2D(IM,JM) C---------------- 3-D ARRAYS -------------------------------------- COMMON/BLK3D/ 1 U(IM,JM,KB), 2 V(IM,JM,KB),W(IM,JM,KB), 3 T(IM,JM,KB),S(IM,JM,KB),RHO(IM,JM,KB),RMEAN(IM,JM,KB) C-------------------------------------------------------------------- dimension fx(im,jm),fy(im,jm) DIMENSION XTIC(NXTIC),YTIC(NYTIC) DIMENSION PROP(IM,JM,KB) DIMENSION ZPROP(IM,JM),ZU(IM,JM),ZV(IM,JM) DIMENSION PI(20),PJ(20),XP(20),YP(20) CHARACTER*59 TITLE CHARACTER*20 TCI CHARACTER*20 TIM CHARACTER*20 TTIME c CHARACTER*3 XLABEL,YLABEL CHARACTER*3 XLAB(NXTIC),YLAB(NYTIC) DATA ISKP/1/ c-----------------------activate graphics routines-------------------- CALL OPNGKS c-----------------------turn off clipping indicator------------------- c-----------------------to plot labels outside map-------------------- CALL GSCLIP(0) C------------------------define axis and titles----------------------- TIME=0. DATA XL,XR,YB,YT/0.,350.,0.,260./ DATA XTIC/0.,50.,100.,150.,200.,250.,300.,350./ DATA YTIC/0.,50.,100.,150.,200./ DATA XLAB/'0','50','100','150','200','250','300','350'/ c DATA ZLAT/'30','35','40','45'/ DATA XLABEL/'X'/ DATA YLABEL/'Y'/ c------------------------plots in x and y--------------------------- JSIZE=12 C----- plotting parameters ---------------------------------------- DAYS=2. ZLEV=0. NITS=4 NSM=10 LMAX=1 INCLUDE 'params' C-------------------------------------------------------------- C*** READ GEOMETRY AND INITIAL CONDITIONS *************************** REWIND 40 DO 1000 LOOP=1,LMAX READ(40) TIME,Z,ZZ,X,Y,DX,DY,H,FSM,U,V,UA,VA,T DO I=1,IM DO J=1,JM X2D(I,J)=1.E-3*X(I) Y2D(I,J)=1.E-3*Y(J) ENDDO ENDDO CALL PRXY('H',0.,H,IM,2,JM,2,0.) CALL PRXY('X2D',0.,X2D,IM,2,JM,2,0.) CALL PRXY('Y2D',0.,Y2D,IM,2,JM,2,0.) C CALL PRXY ('FSM ',TIME,FSM ,IM,2,JM,2,0.0 ) ISK=2 JSK=2 CALL PRXZ('TB' ,TIME,T,IM,ISK,JM,10,25,KB,0.,DT,ZZ) C---------------------------------------------- DO I=1,IM DO J=1,JM DO K=1,KB PROP(I,J,K)=T(I,J,K) ENDDO ENDDO ENDDO c-------------------- Set up and frame plot -------------------------- X1=0.1 X2=0.95 Y2=0.95 Y1=Y2-(X2-X1)*(YT-YB)/(XR-XL) CALL SET(X1,X2,Y1,Y2,XL,XR,YB,YT,1) call GSLWSC(5.) call frstpt(XL,YB) call vector(XR,YB) call frstpt(XL,YT) call vector(XR,YT) call GSLWSC(1.) c-------------------- Insert tic marks ------------------------------- TCSZ=(YT-YB)/50. DO 50 I =1,NYTIC call frstpt(XL,YTIC(I)) call vector(XL+TCSZ,YTIC(I)) CALL PWRITX(XL-2.*TCSZ,YTIC(I),YLAB(I),3,JSIZE,0,0) 50 CONTINUE DO 55 I =1, NXTIC call frstpt(XTIC(I),YB) call vector(XTIC(I),YB+TCSZ) CALL PWRITX(XTIC(I),YB-2.*TCSZ,XLAB(I),3,JSIZE,0,0) 55 CONTINUE c-------------------------------------------------------------------- CMIN=6.0 CINT=.05 CMAX=8.0 ICINT=CINT WRITE(TCI,'(''CI = '',f4.2,'' DEG '')') CINT CALL PWRITX(XR-80,YB-30 ,TCI,20,JSIZE,0,-1) WRITE(TIM,'(''TIME = '',f4.0,'' DAYS '')') TIME CALL PWRITX(XL+80,YB-30 ,TIM,20,JSIZE,0,-1) c WRITE(TTIME,'(''TIME = '',F4.0 ,'' days '')')TIME c CALL PWRITX(XL,YB-30 ,TTIME,20,JSIZE,0,-1) c---------plot contours---------------------------------------------- ZLEV=-2000. CALL SIGTOZ(ZZ,H,PROP,ZPROP,ZLEV,TPS,IM,JM,KB) CALL SMOOTH(ZPROP,TPS,IM,JM,NITS) CALL PRXY(' T ',TIME,ZPROP,IM,2,JM,2,0.0) CALL MAKCON(0,IM,JM,X2D,Y2D,ZPROP,TPS,CMIN,CINT,CMAX,1) c---------plot vectors----------------------------------------------- CALL SIGTOZ(ZZ,H,U,ZU,ZLEV,DUM,IM,JM,KB) CALL SIGTOZ(ZZ,H,V,ZV,ZLEV,DUM,IM,JM,KB) IF(nits.gt.0) then CALL SMOOTH(ZU,DUM,TPS,IM,JM,nits) CALL SMOOTH(ZV,DVM,TPS,IM,JM,nits) ENDIF C------------------------------------------------------------------- c CALL PRXY ('ZU ' , 0., UZ ,IM,2 ,JM,2 ,0.001 ) c CALL PRXY ('ZV ' , 0., VZ ,IM,2 ,JM,2 ,0.001 ) C DO 300 I=1,IM-1 DO 300 J=2,JM-1 IF(FSM(I,J).LT.0.1) GOTO 300 C Select launch points IF(MOD(J,5).EQ.0.and.MOD(I+4,5).EQ.0) 2 THEN C CALL TRAJ(I,J,DAYS,ZU,ZV,DX,DY,IM,JM,PI,PJ,NS,NSM) write(6,'('' NSM,NS ='', 2I6)') NSM,NS C Solve for XP, YP and convert to plot coordinates. DO N=1,NS IP=PI(N) JP=PJ(N) XP(N)=X2D(IP,JP) 1 +(X2D(IP+1,JP)-X2D(IP,JP))*(PI(N)-FLOAT(IP)) 2 +(X2D(IP,JP+1)-X2D(IP,JP))*(PJ(N)-FLOAT(JP)) YP(N)=Y2D(IP,JP) 1 +(Y2D(IP+1,JP)-Y2D(IP,JP))*(PI(N)-FLOAT(IP)) 2 +(Y2D(IP,JP+1)-Y2D(IP,JP))*(PJ(N)-FLOAT(JP)) ENDDO write(6,'(/,5x,10F10.4)') (XP(n),n=1,NS) write(6,'(5x,10F10.4)') (YP(n),n=1,NS) CALL CURVED(XP,YP,NS) SS=0. DO N=2,NS SS=SS+SQRT((PI(N)-PI(N-1))**2+(PJ(N)-PJ(N-1))**2) ENDDO c write(6,'(5x,10F10.4)') (PI(n),n=1,NS) c write(6,'(5x,10F10.4)') (PJ(n),n=1,NS) c write(6,'('' SS ='',F10.4)') SS IF(SS.GE.0.006.AND.NS.EQ.NSM) THEN CALL ARRWHD(XP(NS-2),XP(NS),YP(NS-2),YP(NS),TCSZ,0.4) ENDIF C ENDIF 300 CONTINUE C CALL FRAME 1000 CONTINUE CALL CLSGKS STOP END C SUBROUTINE SIGTOZ(ZZ,H,T,TLEV,ZLEV,FSM,IM,JM,KB) C------------------------------------------------------------------- C THIS ROUTINE LINERLY INTERPOLATES TLEV AT THE LEVEL, C ZLEV, FROM T LOCATED ON SIGMA LEVELS, ZZ. C NOTE THAT A NEW MASK ,FSM, IS CREATED. C------------------------------------------------------------------- DIMENSION ZZ(KB),H(IM,JM),T(IM,JM,KB),TLEV(IM,JM),FSM(IM,JM) do 99 j= 1, jm do 99 i=1, im TLEV(i,j)=0.E0 FSM(i,j)=0.E0 99 continue K=1 DO 200 J=1,JM DO 200 I=1,IM 100 CONTINUE IF(K.LT.1) GO TO 200 IF(K.GT.(KB-1)) THEN K=KB-2 GO TO 200 ENDIF C C FIND K AND K+1 INTERVAL THAT BRACKETS ZLEV, THEN INTERPOLATE C IF(ZLEV.LT.(ZZ(K)*H(I,J))) THEN IF(ZLEV.GE.(ZZ(K+1)*H(I,J))) THEN TLEV(I,J)=T(I,J,K)+(ZLEV-ZZ(K)*H(I,J)) 1 *(T(I,J,K+1)-T(I,J,K))/((ZZ(K+1)-ZZ(K))*H(I,J)) FSM(I,J)=1.E0 GO TO 200 ELSE K=K+1 GO TO 100 ENDIF ELSE K=K-1 GO TO 100 ENDIF 200 CONTINUE RETURN END SUBROUTINE PRXY (LABEL,TIME,A,IM,ISKP,JM,JSKP,SCALA) C >>> C THIS WRITES A 2-D FIELD C TIME=TIME IN DAYS C A = ARRAY(IM,JM) TO BE PRINTED C ISKP=PRINT SKIP FOR I C JSKP=PRINT SKIP FOR J C SCALE=DIVISOR FOR VALUES OF A C DIMENSION A(IM,JM),NUM(150),LINE(150) CHARACTER LABEL*(*) DATA ZERO /1.E-10/ C SCALE=SCALA IF (SCALE.GT.ZERO) GO TO 160 AMX=ZERO DO 150 J=1,JM,JSKP DO 150 I=1,IM,ISKP AMX=MAX1(ABS(A(I,J)),AMX) 150 CONTINUE SCALE=10.**(INT(LOG10(AMX)+1.E2)-103) 160 CONTINUE SCALEI=1./SCALE WRITE(6,170) LABEL 170 FORMAT(1X,A40) WRITE(6,180) TIME,SCALE 180 FORMAT(' TIME =',F9.2,' DAYS MULTIPLY ALL VALUES BY',1PE10.3) DO 190 I=1,IM 190 NUM(I)=I IB=1 C 200 CONTINUE IE=IB+23*ISKP IF(IE.GT.IM) IE=IM WRITE(6,210) (NUM(I),I=IB,IE,ISKP) 210 FORMAT(/,2X,24I5,/) DO 260 J=1,JM,JSKP JWR=JM+1-J DO 220 I=IB,IE,ISKP 220 LINE(I)=INT(SCALEI*A(I,JWR)) WRITE(6,240) JWR,(LINE(I),I=IB,IE,ISKP) 240 FORMAT(1X,I3,24I5) 260 CONTINUE WRITE(6,280) 280 FORMAT(//) IF(IE.GE.IM) RETURN IB=IB+24*ISKP GO TO 200 END C SUBROUTINE PRXZ(LABEL,TIME,A,IM,ISKP,JM,J1,J2,KB,SCALA,DT,ZZ) DIMENSION A(IM,JM,KB),NUM(100),LINE(100),IDT(100), 1 ZZ(KB),DT(IM,JM) CHARACTER LABEL*(*) DATA ZERO/1.E-10/ C C THIS SUBROUTINE WRITES VERTICAL SECTIONS OF A 3-D FIELD C C TIME=TIME IN DAYS C A= ARRAY(IM,JM,KB) TO BE PRINTED C ISPL=PRINT SKIP FOR I C JSPL=PRINT SKIP FOR J C SCALE=DIVISOR FOR VALUES OF A C SCALE=SCALA IF (SCALE.GT.ZERO) GO TO 150 AMX=ZERO DO 140 K=1,KB DO 140 J=1,JM DO 140 I=1,IM,ISKP AMX=MAX(ABS(A(I,J,K)),AMX) 140 CONTINUE IF(AMX.EQ.0.) THEN SCALEI=0. GOTO 165 ENDIF SCALE=10.E0**(INT(LOG10(AMX)+1.E2)-103) 150 CONTINUE SCALEI=1./SCALE 165 CONTINUE WRITE(6,9980) LABEL 9980 FORMAT(' ',A30) WRITE(6,9981) TIME,SCALE 9981 FORMAT(' TIME = ',F9.3,' MULTIPLY ALL VALUES BY',E8.2) DO 10 JPP=1,2 J=J1 IF(JPP.EQ.2) J=J2 IF(J.EQ.0) RETURN WRITE(6,30) J 30 FORMAT(3X,/' SECTION J =',I3) IB=1 IE=IB+23*ISKP 50 CONTINUE IF(IE.GT.IM) IE=IM DO 100 I=IB,IE,ISKP IDT(I)=DT(I,J) 100 NUM(I)=I WRITE(6,9990) (NUM(I),I=IB,IE,ISKP) 9990 FORMAT(/,' I = ',24I5,/) WRITE(6,9991) (IDT(I),I=IB,IE,ISKP) 9991 FORMAT(8X,'D =',24I5.0,/,' ZZ ') DO 200 K=1,KB DO 190 I=IB,IE,ISKP 190 LINE(I)=SCALEI*A(I,J,K) WRITE(6,9966) K,ZZ(K),(LINE(I),I=IB,IE,ISKP) 9966 FORMAT(1X,I2,2X,F6.3,24I5) 200 CONTINUE WRITE(6,9984) IF(IE.EQ.IM) GO TO 10 IB=IB+24*ISKP IE=IB+23*ISKP GO TO 50 9984 FORMAT(//) 10 CONTINUE RETURN END C SUBROUTINE SMOOTH(A,FSM,IM,JM,NITS) C------------------------------------------------------------------- C THIS ROUTINE SMOOTHS DATA WITH A FIVE POINT LAPLACIAN FILTER. C------------------------------------------------------------------- DIMENSION A(IM,JM),FSM(IM,JM) DO 100 N=1,NITS DO 200 J=2,JM-1 DO 200 I=2,IM-1 SMFAC=FSM(I+1,J)+FSM(I,J-1)+FSM(I-1,J)+FSM(I,J+1)+1.E-5 A(I,J)=A(I,J)+(.5/SMFAC) 1 *(A(I+1,J)*FSM(I+1,J)+A(I,J-1)*FSM(I,J-1) 2 +A(I-1,J)*FSM(I-1,J)+A(I,J+1)*FSM(I,J+1) 3 -SMFAC*A(I,J)) A(I,J)=A(I,J)*FSM(I,J) 200 CONTINUE 100 CONTINUE RETURN END C----------------------------------------------------------------------- SUBROUTINE MAKCON(IND,IM,JM,X,Y,T,FSM,CMIN,CINT,CMAX,NU) C---------------------------------------------------------------------- C CREATES CONTOURS IN X,Y SPACE FROM T FIELD. X,Y AND T ARE 2-D C ARRAYS AND THEREFORE CAN REPRESENT ANY COORDINATE SYSTEM. C A SPECIAL NON-LINEAR INTERPOLATING FEATURE INVOLVING THE C FUNCTION COR HAS BEEN REMOVED FROM THIS VERSION. C C IND = 0; NEGATIVE CONTOURS ARE DASHED. IND = 1; ALL CONTOURS ARE SOLID. C NU = 1; CONTOURS ARE NUMBERED C IM,JM = ARRAY SIZE C X,Y = PLOTTING COORDINATES C T = VARIABLE TO BE CONTOURED C FSM(I,J) = 0; GRID IS MASKED OUT. FSM(I,J)=1; GRID IS PLOTTED. C CMIN,CINT,CMAX = CONTOUR MINIMUM, INTERVAL AND MAXIMUM C---------------------------------------------------------------------- PARAMETER (LMAX=10000) DIMENSION X(IM,JM),Y(IM,JM),T(IM,JM),FSM(IM,JM),C(500) DIMENSION XF(LMAX),YF(LMAX),XL(LMAX),YL(LMAX),XC(LMAX),YC(LMAX) CHARACTER*16 IPAT CHARACTER*5 JCO REAL XCR(LMAX),YCR(LMAX) DATA X1,X2,X3,X4,X5,X6,X7,X8/8*0./ DATA Y1,Y2,Y3,Y4,Y5,Y6,Y7,Y8/8*0./ EC=.00010*ABS(CMAX-CMIN) EC=0.0 DO 220 K=1,500 C(K)=CMIN+CINT*(K-1)+EC IF(C(K).GE.CMAX) GO TO 221 220 KCM=K 221 CONTINUE C----------------SOLVE FOR C(K) CONTOUR -------------------------------- DO 9999 K=1,KCM CO=C(K) C**** THIS SECTION COLLECTS ALL CELL VECTORS FOR EACH C(K) CONTOUR. **** C**** S=ANY NUMBER NOT IN (X,Y) RANGE. **** C**** REINSTATE ABS'S EX=.00001*ABS(X(1,JM)-X(1,1))+.00001*ABS(X(IM,JM)-X(IM,1)) EY=.00001*ABS(Y(1,JM)-Y(1,1))+.00001*ABS(Y(IM,JM)-Y(IM,1)) S=-1000. IMM1=IM-1 JMM1=JM-1 L=0 DO 960 J=1,JMM1 DO 960 I=1,IMM1 IF(L.GT.(LMAX-5)) THEN WRITE(6,'('' INCREASE LMAX IN MAKCON TO COMPLETE CONTOUR'')') GO TO 9990 ENDIF TSW=T(I,J) TSE=T(I+1,J) TNE=T(I+1,J+1) TNW=T(I,J+1) IF(CO.LE.TSW.AND.CO.LE.TNW.AND.CO.LE.TNE. 1 AND.CO.LE.TSE) GO TO 960 IF(CO.GT.TSW.AND.CO.GT.TNW.AND.CO.GT.TNE. 1 AND.CO.GT.TSE) GO TO 960 X1=S X2=S X3=S X4=S X5=S X6=S X7=S X8=S XSE=X(I+1,J) XSW=X(I,J) XNE=X(I+1,J+1) XNW=X(I,J+1) YSE=Y(I+1,J) YSW=Y(I,J) YNE=Y(I+1,J+1) YNW=Y(I,J+1) C****** FIND CELL INTERSECTIONS WITH CONTOUR *************************** IF(CO.GT.TNW.AND.CO.LT.TSW) GO TO 430 IF(CO.LT.TNW.AND.CO.GT.TSW) GO TO 430 GO TO 470 430 RA=(CO-TSW)/(TNW-TSW+.01*EC) RACOR=RA Y1=YSW+RACOR*(YNW-YSW) X1=XSW+RACOR*(XNW-XSW) 470 IF(CO.GT.TSE.AND.CO.LT.TSW) GO TO 480 IF(CO.LT.TSE.AND.CO.GT.TSW) GO TO 480 GO TO 520 480 RA=(CO-TSW)/(TSE-TSW+.01*EC) RACOR=RA X3=XSW+RACOR*(XSE-XSW) Y3=YSW+RACOR*(YSE-YSW) 520 IF(CO.GT.TSE.AND.CO.LT.TNE) GO TO 530 IF(CO.LT.TSE.AND.CO.GT.TNE) GO TO 530 GO TO 560 530 RA=(CO-TSE)/(TNE-TSE+.01*EC) RACOR=RA X5=XSE+RACOR*(XNE-XSE) Y5=YSE+RACOR*(YNE-YSE) 560 IF(CO.GT.TNW.AND.CO.LT.TNE) GO TO 570 IF(CO.LT.TNW.AND.CO.GT.TNE) GO TO 570 GO TO 600 570 RA=(CO-TNW)/(TNE-TNW+.01*EC) RACOR=RA X7=XNW+RACOR*(XNE-XNW) Y7=YNW+RACOR*(YNE-YNW) 600 CONTINUE C****** CONNECT HEAD AND TAILS OF CONTOUR VECTORS ****************** C ONE=1.0 IF(FSM(I,J).LT.ONE.OR.FSM(I,J+1).LT.ONE) X1=S IF(FSM(I,J+1).LT.ONE.OR.FSM(I+1,J+1).LT.ONE) X3=S IF(FSM(I+1,J).LT.ONE.OR.FSM(I+1,J+1).LT.ONE) X5=S IF(FSM(I,J).LT.ONE.OR.FSM(I+1,J).LT.ONE) X7=S IF(X3.EQ.S) GO TO 930 IF(X1.EQ.S) GO TO 910 L=L+1 XF(L)=X1 YF(L)=Y1 XL(L)=X3 YL(L)=Y3 IF(X5.NE.S) GO TO 950 910 IF(X5.EQ.S) GO TO 920 L=L+1 XF(L)=X5 YF(L)=Y5 XL(L)=X3 YL(L)=Y3 920 IF(X7.EQ.S.OR.X5.NE.S.OR.X1.NE.S) GO TO 930 L=L+1 XF(L)=X7 YF(L)=Y7 XL(L)=X3 YL(L)=Y3 930 IF(X1.EQ.S) GO TO 950 IF(X5.EQ.S.OR.X3.NE.S.OR.X7.NE.S) GO TO 940 L=L+1 XF(L)=X1 YF(L)=Y1 XL(L)=X5 YL(L)=Y5 940 IF(X7.EQ.S) GO TO 950 L=L+1 XF(L)=X1 YF(L)=Y1 XL(L)=X7 YL(L)=Y7 950 IF(X5.EQ.S.OR.X7.EQ.S) GO TO 960 L=L+1 XF(L)=X5 YF(L)=Y5 XL(L)=X7 YL(L)=Y7 960 CONTINUE LREM=L IF(LREM.LE.1) GO TO 9999 C****** SORT ASSEMBLAGE OF VECTORS TO FORM CONTOURS *********** 9990 CONTINUE IPASS=1 L=1 XC(1)=XF(L) YC(1)=YF(L) XC(2)=XL(L) YC(2)=YL(L) M=2 C 20 L=L+1 DO 50 LLL=L,LREM XF(LLL-1)=XF(LLL) YF(LLL-1)=YF(LLL) XL(LLL-1)=XL(LLL) YL(LLL-1)=YL(LLL) 50 CONTINUE LREM=LREM-1 60 CONTINUE IF(LREM.LT.1) GO TO 106 DO 105 L=1,LREM IF(XF(L).GT.(XC(M)+EX).OR.XF(L).LT.(XC(M)-EX)) GO TO 100 IF(YF(L).GT.(YC(M)+EY).OR.YF(L).LT.(YC(M)-EY)) GO TO 100 M=M+1 XC(M)=XL(L) YC(M)=YL(L) GO TO 20 100 CONTINUE IF(XL(L).GT.(XC(M)+EX).OR.XL(L).LT.(XC(M)-EX)) GO TO 105 IF(YL(L).GT.(YC(M)+EY).OR.YL(L).LT.(YC(M)-EY)) GO TO 105 M=M+1 XC(M)=XF(L) YC(M)=YF(L) GO TO 20 105 CONTINUE C***** CHECK TO SEE IF CLOSED CURVE ******************************** IF(XC(1).LT.(XC(M)+EX).AND.XC(1).GT.(XC(M)-EX).AND. 1 YC(1).LT.(YC(M)+EY).AND.YC(1).GT.(YC(M)-EY)) GO TO 106 C**** OTHERWISE CURVE TERMINATES ON EDGE; THEN FIND SECOND PIECE *** IF(LREM.GT.2.AND.IPASS.EQ.1) THEN CALL TRANPOS(M,XC,YC) IPASS=2 GO TO 60 ENDIF 106 CONTINUE C**** WEED OUT CLOSELY SPACED POINTS ********************************** ME=M M=1 EXY=100.0*(EX+EY) XCR(1)=XC(1) YCR(1)=YC(1) DO 110 I=2,ME IF(ABS(XC(I)-XC(M)).LT.EXY.AND.ABS(YC(I)-YC(M)).LT.EXY) GO TO 110 M=M+1 XCR(M)=XC(I) YCR(M)=YC(I) 110 CONTINUE C************ CONTOUR COMPLETED **************************************** C DASH: JCRT= length/$ ; JSIZE= size of char. in numbered lines C COSKP = # of contour lines where labelling is skipped WRITE(6,'('' CO ='',F10.3,'' L,LREM,M ='',3I10)') CO,L,LREM,M JCRT=10 JSIZE=10 COSKP=2 IPAT='$$$$$$$$$$$$$$$$' IF (CO.LT.0.0.AND.IND.EQ.0)IPAT='$$''$$''$$''$$''' c ICO=INT(CO) C WRITE(JCO,'(I4 )')ICO WRITE(JCO,'(F4.1)')CO CALL DELZER(JCO,5) WRITE(6,'('' JCO= '',A5)')JCO c COINT=COSKP*CINT kCO=(ABS(CO)+.000001)/COINT IF(NU.EQ.1.AND.ABS(CO)-kCO*COINT.LT.0.00001) 1 WRITE(IPAT,'(A11,A5)')IPAT,JCO c IF(NU.EQ.1.AND.CO.EQ.0.)WRITE(IPAT,'(A12,A4)')IPAT,' O ' IF(NU.EQ.1.AND.CO.EQ.0.)JCRT=10 c CALL DASHDC(IPAT,JCRT,JSIZE) CALL CURVED(XCR,YCR,M) c WRITE(6,'('' IPAT ='',A16)')IPAT 9988 IF(LREM.GT.2) GO TO 9990 9999 CONTINUE RETURN END C---------------------------------------------------------------------- SUBROUTINE TRANPOS(M,XC,YC) DIMENSION XC(1000),YC(1000),XT(1000),YT(1000) DO 100 I=1,M XT(M-I+1)=XC(I) 100 YT(M-I+1)=YC(I) DO 200 I=1,M XC(I)=XT(I) 200 YC(I)=YT(I) RETURN END subroutine julia(time,iyear,month,iday) c julian=int(time) c ----------calculate year---------- iy=julian/360 iyear=iy+1 c ----------calculate day----------- my=360*iy iday=julian-my if(iday.eq.0.or.iday.eq.360) then iyear=iy iday=360 endif c ----------calculate month--------- month=(iday-1)/30+1 return end c SUBROUTINE DELZER (ICH,N) C C Replace ZERO in character ICH with Os for better contour labeling C CHARACTER*5 ICH CHARACTER*1 IC(10),IOH,IZR IOH='O' IZR='0' DO 10 I=1,N IC(I)=ICH(I:I) IF(IC(I).EQ.IZR) IC(I)=IOH 10 CONTINUE WRITE(ICH,'(5A1)')(IC(I),I=1,N) RETURN END c SUBROUTINE FINDPSI PARAMETER (IM=65,JM=49,KB=11,IMM1=IM-1,JMM1=JM-1,KBM1=KB-1) COMMON/BLK2D/H(IM,JM),DX(IM,JM),DY(IM,JM), 1 ART(IM,JM),ARU(IM,JM),ARV(IM,JM),ALON(IM,JM),ALAT(IM,JM), 2 DUM(IM,JM),DVM(IM,JM),FSM(IM,JM),COR4(IM,JM),ANG(IM,JM), 3 TPS(IM,JM), 4 UAB(IM,JM), 5 VAB(IM,JM),ELB(IM,JM),TLEV(IM,JM),PSI(IM,JM), 6 WX1(IM,JM,2),WX2(IM,JM,2) DIMENSION D(IM,JM) C---------------- 3-D ARRAYS -------------------------------------- DO 9004 J=1,JM DO 9004 I=1,IM D(I,J)=H(I,J)+ELB(I,J) 9004 PSI(I,J)=0. DO 9005 J=2,JM DO 9005 I=1,IMM1 9005 PSI(I+1,J)=(PSI(I,J)-VAB(I,J)*.25E0*(D(I,J)+D(I,J-1)) 1 *(DX(I,J)+DX(I,J-1)))*FSM(I+1,J) CALL PRXY(' PSI FROM V ',TIME,PSI,IM,3,JM,3,1.E4) DO 9006 J=1,JM DO 9006 I=1,IM 9006 PSI(I,J)=0. DO 9007 J=2,JMM1 DO 9007 I=1,IM 9007 PSI(I,J+1)=(PSI(I,J)+.25E0*UAB(I,J)*(D(I,J)+D(I-1,J)) 1 *(DY(I,J)+DY(I-1,J)))*FSM(I,J+1) CALL PRXY(' PSI FROM U ',TIME,PSI,IM,3,JM,3,1.E4) RETURN END C SUBROUTINE TRAJ(IL,JL,DAYS,U,V,DX,DY,IM,JM,PI,PJ,NS,NSM) c c This code creates "Eulerian trajectories" in a c continuous I, J coordinate system called PI, PJ. c c IL,JL = launch point c U, V = horizontal velocity fields c DX, DY = Size of grid elements. c DAYS = # of days for the trajectory. c NSM = specified # of trajectory segments c NS = actual # of trajectory segments c C---------------- 2-D ARRAYS -------------------------------------- DIMENSION U(IM,JM),V(IM,JM),DX(IM,JM),DY(IM,JM) DIMENSION PI(NSM),PJ(NSM) TTRAJ=DAYS*86400. DTSTEP=TTRAJ/NSM PI(1)=IL PJ(1)=JL IP=PI(1)+0.5 JP=PJ(1)+0.5 NS=1 DO 100 L=2,NSM c Take half step to evaluate interpolated velocities. UP=U(IP,JP)+(U(IP+1,JP)-U(IP,JP))*(PI(L-1)-FLOAT(IP)+0.5) VP=V(IP,JP)+(V(IP,JP+1)-V(IP,JP))*(PJ(L-1)-FLOAT(JP)+0.5) c write(6,'('' PI(L-1),UP,DX,DTSTEP'',4F13.3)') c 1 PI(L-1),UP,DX(I,J),DTSTEP C23456 | | | | | | | xxxxxxx| PI(L)=PI(L-1)+0.5*DTSTEP*UP/DX(IP,JP) PJ(L)=PJ(L-1)+0.5*DTSTEP*VP/DY(IP,JP) IP=PI(L) JP=PJ(L) c write(6,'('' PI(L)'',F10.3)') PI(L) c c Take full step UP=U(IP,JP)+(U(IP+1,JP)-U(IP,JP))*(PI(L)-FLOAT(IP)+0.5) VP=V(IP,JP)+(V(IP,JP+1)-V(IP,JP))*(PJ(L)-FLOAT(JP)+0.5) c write(6,'('' PI(L-1),UP,DX,DTSTEP'',4F13.3)') c 1 PI(L-1),UP,DX(I,J),DTSTEP IF(ABS(UP).GT.1.E-10) PI(L)=PI(L-1)+ DTSTEP*UP/DX(IP,JP) IF(ABS(VP).GT.1.E-10) PJ(L)=PJ(L-1)+ DTSTEP*VP/DY(IP,JP) IP=PI(L)+0.5 JP=PJ(L)+0.5 c write(6,'('' PI(L)'',F10.3)') PI(L) IF(PI(L).LT.1.or.PI(L).GT.(IM)) RETURN IF(PJ(L).LT.1.or.PJ(L).GT.(JM)) RETURN NS=NS+1 100 CONTINUE RETURN END c SUBROUTINE ARRWHD(xf,xl,yf,yl,L,f) C This draws an arrowhead. The length is L. C The half width is f*L. (xf, yf) and (xl, yl) C defines the arrow direction. real L dimension x(4),y(4) x(1)=xl y(1)=yl x(4)=xl y(4)=yl XX=xl-xf YY=yl-yf SS=SQRT(XX**2+YY**2) IF(SS.LT.0.000001) RETURN sinang=YY/SS cosang=XX/SS W=f*L x(2)=xl-L*cosang-W*sinang y(2)=yl-L*sinang+W*cosang x(3)=xl-L*cosang+W*sinang y(3)=yl-L*sinang-W*cosang CALL GSLWSC(2.) CALL CURVED(x,y,4) CALL GSLWSC(1.) c call GSCR(1,5,0.,0.,0.) c call SFSETI('TYPE OF FILL',0) c call SFSGFA(X,Y,3,DST,200,IND,350,5) RETURN END