C ********************************************************** C C PLOTS MERIDIONAL HEAT FLUX (INTEGRATED ACROSS ATLANTIC C ON ZIG-ZAG LINES ESTIMATED LATITUDES C (May 93- improvements by glm for MED) C (Dec 93- incl. comparison with ISEMER et al 1989. c ******************** NSATL MODEL ************************** PROGRAM PLXZ C PARAMETER (IM=128,JM=128,KB=16, IMJM=IM+JM) C C Y,Z grid for interpolation PARAMETER ( LY= 112) PARAMETER (AWEST=-100.,AEAST=30.,ASOUTH=-33.,ANORTH=80.) C DIMENSION XX(LY),YY(LY) C DIMENSION OBSX(11),OBSY(11) DIMENSION IZIG(IMJM,LY),JZIG(IMJM,LY),NPTS(LY),EY(LY) DIMENSION VXZ(IM,KB),FLU(IM,JM,KB),FLV(IM,JM,KB) C COMMON/BLK1D/Z(KB),ZZ(KB),DZ(KB),DZZ(KB) C---------------- 2-D ARRAYS ------------------------------------- COMMON/BLK2D/H(IM,JM),DX(IM,JM),DY(IM,JM),FSM(IM,JM), 1 ALON(IM,JM),ALAT(IM,JM),ELB(IM,JM),COR(IM,JM), 3 DUM(IM,JM),DVM(IM,JM),ART(IM,JM),ARU(IM,JM),ARV(IM,JM), 4 UA(IM,JM),VA(IM,JM),UAB(IM,JM),VAB(IM,JM),D(IM,JM) C---------------- 3-D ARRAYS ------------------------------------- COMMON/BLK3D/TB(IM,JM,KB),SB(IM,JM,KB), 1 U(IM,JM,KB),V(IM,JM,KB) 2 ,UM(IM,JM,KB),VM(IM,JM,KB) C-------------------------------------------------------------------- CHARACTER*59 TITLE,TIT(2) CHARACTER*15 XLABEL,YLABEL CHARACTER*10 DATE CHARACTER*2 YEAR CHARACTER*80 fts,fuv C C Observed Heat transport from ISEMER et al. 1989 C DATA OBSX /10., 15., 20., 25., 30., 35., 40., 45., 50., 55., 60./ DATA OBSY /.9, .98, 1.0, .98, .92, .8 , .6 , .5 , .4 , .32, .25 / C C ************** INPUT/OUTPUT FILES ********************************* C grid, bath. and LEV. IC C low res. C OPEN(40,FILE='/archive/o/tne/NATL/initHTS.LEV', C 2 FORM='unformatted') C high res. C OPEN(40,FILE='/archive/tne/NSATL/inHTS.128x', C 2 FORM='unformatted') C C include 'infile' C C OPEN(50,FILE=fts,FORM='unformatted') C OPEN(51,FILE=fuv,FORM='unformatted') C C OPEN(51,FILE='/archive/tne/NSATL/run50-74/ts.1952', C 2 FORM='unformatted') C OPEN(52,FILE='/archive/tne/NSATL/run50-74/uv.1952', C 2 FORM='unformatted') C C ************* OUTPUT FILES *************************************** C OPEN(6,FILE='plot.out',FORM='formatted') C C***************** read geometry & initial cond. ******************** C READ(40) KBB,Z,ZZ,DZ,DZZ,IMM,JMM,ALON,ALAT,DX,DY,H,FSM, 1 DUM,DVM,ART,ARU,ARV,COR,TB,SB DUM=1. DVM=1. C------------------- open GRAPHICS & define B&W CALL OPNGKS CALL GSCR(1,1,0.,0.,0.) CALL GSCR(1,0,1.,1.,1.) c c ************** define axis and titles ********************** c TIME=0. C IY=4 data XL,XR,YB,YT/-40.,70.,-0.5,1.5/ C data XL,XR,YB,YT/-40.,70.,-2.,2./ C IX=11 C IX=ABS(XR-XL)/100. C IY=ABS(YT-YB)/100. X1=0.5*IX X2=X1 Y1=0.5*IY Y2=Y1 XCNT=XL+0.5*IX YCNT=YB+0.5*IY c ------SZ=SIZE(0-1); SQ=DY/DX ; keep SZ < SQ -------- SZ=0.5 SQ=0.6 C C TIT(1)='(A) MERIDIONAL HEAT TRANSPORT - 4680 DAYS' C TIT(2)='(B) MERIDIONAL HEAT TRANSPORT - 4680 DAYS' C XLABEL='LATITUDE (N) ' YLABEL='HEAT FLUX (PW)' C c *************** NREAD= number of pages C IPLTS= 1 or 2 plots/page NREAD=1 IPLTS=1 C IPLTS=2 C NPLOT=0 DO 2000 NPAGE=1,NREAD DO 1000 NY=1,IPLTS DO 1000 NX=1,1 C ISTR=36 C C IF(NY.EQ.2)ISTR=18 C***************** read model simulations (5d interv. ) ************ SB=0. UM=0. VM=0. D=0. DO 29 KDAY=1,ISTR READ(51) TIME,DATE,ELB,TB READ(61) TIME,DATE,UAB,VAB,U,V SB=SB+TB/ISTR D=D+(H+ELB)/ISTR UM=UM+U/ISTR VM=VM+V/ISTR 29 CONTINUE TB=SB U=UM V=VM C DO 104 I=1,IM DO 104 K=1,KB IF(I.LE.66.OR.I.GE.186) V(I,2,K)=0. 104 CONTINUE V(:,1,:)=V(:,2,:) C C WRITE(TITLE,'(''MERIDIONAL HEAT TRANSPORT '',A10)') DATE YEAR=DATE(9:10) WRITE(TITLE,'(''MER. HEAT TRANSPORT (FLSL-C) YEAR= '',A2)')YEAR C WRITE(6,'(''READ '',F10.2,A10)')TIME,DATE C AA=0. VV=0. AX=0. AZ=0. FSZ=1. C C ---- FIND ZIGZAG LINE APPROX. ALONG LATITUDE ------- C DO 106 J=1,LY EY(J)=ASOUTH+(ANORTH-ASOUTH)*(J-1)/LY BLT=EY(J) CALL ZIGZAG(ALON,ALAT,IM,JM,BLT,IZIG,JZIG,IMJM,LY,J,NPTS) 106 CONTINUE C C C ### CALC. HEAT FLUX ACROSS BASIN ON ZIGZAG LINES ######## WRITE(6,'(/,'' LAT(N) HF(PW) '',/)') C C -- factor to change units from V*T to W/M**2 C FAC=4.1876E+6 C DO 115 JJ=1,LY C VTOT=0. C DO 110 K=1,KB-1 DO 110 L=2,NPTS(JJ) C IUV=JZIG(L-1,JJ)-JZIG(L,JJ) I=IZIG(L,JJ) J=JZIG(L,JJ) C IF(IUV.EQ.0)THEN DXH=1.E-15*0.25*((DX(I,J)+DX(I,J-1))*(D(I,J)+D(I,J-1)) ) C DXH=1.E-15*0.25*((DX(I,J)+DX(I,J-1))*(H(I,J)+H(I,J-1)) ) TBTB=0.5*(TB(I,J,K)+TB(I,J-1,K)) VTOT=VTOT+DXH*V(I,J,K)*TBTB*FAC*(Z(K)-Z(K+1)) C VTOT=VTOT+DXH*V(I,J,K)*TB(I,J,K)*FAC*(Z(K)-Z(K+1)) ELSE I=I+1 C*********************** correction NOV96 IF(IUV.EQ.-l)J=J-1 DYH=1.E-15*0.25*((DY(I,J)+DY(I-1,J))*(D(I,J)+D(I-1,J))) C DYH=1.E-15*0.25*((DY(I,J)+DY(I-1,J))*(H(I,J)+H(I-1,J))) TBTB=0.5*(TB(I,J,K)+TB(I-1,J,K)) VTOT=VTOT+DYH*U(I,J,K)*IUV*TBTB*FAC*(Z(K)-Z(K+1)) C VTOT=VTOT+DYH*U(I,J,K)*IUV*TB(I,J,K)*FAC*(Z(K)-Z(K+1)) ENDIF C 110 CONTINUE WRITE(6,'(I5,2X,2F8.2)')JJ,EY(JJ),VTOT C XX(JJ)=EY(JJ) YY(JJ)=VTOT 115 CONTINUE C C IF(2.gt.1) STOP C C C ******* STARTS PLOT SECTION ****************************** C NPLOT=NPLOT+1 TIME=-1. C TITLE=TIT(NPLOT) C c **************** plot frame and titles ************************* c CALL GRD(XL,XR,IX,XLABEL,YB,YT,IY,YLABEL,TITLE,SZ,TIME,CINT, 1 NX,NY,IPLTS,SQ) c C C ----- PLOT MERIDIONAL HEAT FLUX ----- C model CALL CURVED(XX,YY,LY) C C observed C SET LINES WIDER THAN USUAL C CALL GSLWSC(3.) C CALL CURVED(OBSX,OBSY,11) C CALL GSLWSC(1.) C 1000 continue call frame 2000 continue call clsgks stop 'O.K.' end C SUBROUTINE ZIGZAG(ALON,ALAT,IM,JM,BLT,IZIG,JZIG,IMJM,LY,LJ,NPTS) C C FIND ZIGZAG PATH ALONG MODEL GRID POINTS CLOSEST TO BLT C DIMENSION ALON(IM,JM),ALAT(IM,JM) DIMENSION IZIG(IMJM,LY),JZIG(IMJM,LY),NPTS(LY) DIMENSION II(4),JJ(4) DATA II /-1, 0, 1, 0/ DATA JJ / 0,-1, 0, 1/ C NPTS(LJ)=0 C ------------------ find first point on west boundary K=1 DO 10 J=2,JM IF(ALAT(2,J).GT.BLT.AND.ALAT(2,J-1).LE.BLT)THEN IZIG(K,LJ)=2 JZIG(K,LJ)=J-1 ENDIF 10 CONTINUE C 30 CONTINUE I=IZIG(K,LJ) J=JZIG(K,LJ) X=ALON(I,J) Y=ALAT(I,J) C XMAX=-100. LMAX=1 C ---------------------------------- find next points DO 40 L=1,4 IL=I+II(L) JL=J+JJ(L) X1=ALON(IL,JL) Y1=ALAT(IL,JL) C IF(X1.GT.XMAX)THEN XMAX=X1 LMAX=L ENDIF C C WRITE(6,'(I5,4F8.2)')L,X,Y,X1,Y1 D=BLT-Y1 DD=BLT-Y IF(X1.LT.X)GO TO 40 IF(ABS(D).GT.ABS(DD).AND.DD*D.GT.0.)GO TO 40 K=K+1 IZIG(K,LJ)=IL JZIG(K,LJ)=JL 40 CONTINUE C -------------------- if no point found take east one C IF(I.EQ.IZIG(K,LJ).AND.J.EQ.JZIG(K,LJ))THEN K=K+1 IZIG(K,LJ)=I+II(LMAX) JZIG(K,LJ)=J+JJ(LMAX) ENDIF C NPTS(LJ)=K IF(I.GE.IM-1) RETURN GO TO 30 C RETURN END C C ------------------------------------------------------------ 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) TLEV=0. FSM=0. 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.LE.(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. GO TO 200 ELSE K=K+1 GO TO 100 ENDIF ELSE K=K-1 GO TO 100 ENDIF 200 CONTINUE RETURN END C C 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(IM),LINE(IM) 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.4,' 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 SUBROUTINE SMOOTH(A,B,FSM,IM,JM,NITS) C------------------------------------------------------------------- C THIS ROUTINE SMOOTHS DATA WITH A FIVE POINT LAPLACIAN FILTER. C------------------------------------------------------------------- C A is the array to be smoothed. C B is work space that must be provided through the arguement list. C FSM=1 if element is to include in the smoothing operator. Otherwise, FSM=0. C NITS is the number of smoothing operations. C DIMENSION A(IM,JM),B(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 B(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)) 200 CONTINUE DO 210 J=1,JM,JM-1 DO 210 I=2,IM-1 SMFAC=FSM(I+1,J)+FSM(I-1,J)+1.E-5 B(I,J)=A(I,J)+(.5/SMFAC) 1 *(A(I+1,J)*FSM(I+1,J) 2 +A(I-1,J)*FSM(I-1,J) 3 -SMFAC*A(I,J)) 210 CONTINUE DO 220 J=2,JM-1 DO 220 I=1,IM,IM-1 SMFAC=FSM(I,J-1)+FSM(I,J+1)+1.E-5 B(I,J)=A(I,J)+(.5/SMFAC) 1 *(A(I,J-1)*FSM(I,J-1) 2 +A(I,J+1)*FSM(I,J+1) 3 -SMFAC*A(I,J)) 220 CONTINUE A=B*FSM 100 CONTINUE RETURN END C SUBROUTINE SMOOTHOLD(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 GRD(XMIN,XMAX,NXINCS,XLABEL,YMIN,YMAX,NYINCS,YLABEL, 1 TITLE,SZ,TIME,CIR,NX,NY,IPLTS,SQ) C C THIS SUBROUTINE DRAWS AND LABELS A GRID FOR PLOTTING C C XMIN = MINIMUM VALUE ON THE X-AXIS C XMAX = MAXIMUM VALUE ON THE X-AXIS C NXINCS = NUMBER OF DIVISIONS ON THE X-AXIS BETWEEN GRID LINES C XLABEL = LABEL FOR THE X-AXIS C YMIN = MINIMUM VALUE ON THE Y-AXIS C YMAX = MAXIMUM VALUE ON THE Y-AXIS C NYINCS = NUMBER OF DIVISIONS ON THE Y-AXIS BETWEEN GRID LINES C YLABEL = LABEL FOR THE Y-AXIS C TITLE = FRAME TITLE C SZ = SIZE OF PLOT (e.g., .5 IS 50% OF MAXIMUM SIZE) C SQ = RATIO DY/DX [note: SZ 1 C IPLTS = NUMBER OF PLOTS PER PAGE C CHARACTER*59 TITLE CHARACTER*15 XLABEL,YLABEL CHARACTER*4 XCHAR,YCHAR CHARACTER*16 TCHAR C C I-----------------------------------------------------I C I COMPUTE THE GRID LINE SPACING AND GRID CENTER I C I-----------------------------------------------------I XAX=ABS(XMAX-XMIN) YAX=ABS(YMAX-YMIN) DX=XAX/NXINCS DY=YAX/NYINCS XCENT=(XMIN+XMAX)/2. YCENT=(YMIN+YMAX)/2. SINT=0.2 XLEFT=XMIN-XAX*SINT YBOTT=YMIN-YAX*SINT XRIGH=XMAX+XAX*SINT YTTOP=YMAX+YAX*SINT C ********************* SFACT=RATIO OF TOTAL PLOT TO CONTOUR FRAME *** SFACT=SINT/(1+2.*SINT)*SZ bt=(1.-SZ)*.5 tp=1.-(1.-SZ)*.5 dbt=tp-bt ft=bt+0.5*dbt*(sq-1)/sq rt=tp-0.5*dbt*(sq-1)/sq dfr=rt-ft C for multiple plots/page IF(IPLTS.EQ.2)THEN bt=bt-0.5*dbt*(2*ny-3) tp=tp-0.5*dbt*(2*ny-3) ENDIF IF(IPLTS.EQ.4)THEN ft=ft+0.5*dfr*(2*nx-3) rt=rt+0.5*dfr*(2*nx-3) bt=bt-0.5*dbt*(2*ny-3) tp=tp-0.5*dbt*(2*ny-3) ENDIF C I--------------------------------I C I DEFINE AND PLOT THE FRAME I C I--------------------------------I CALL SET(ft,rt,bt,tp,XLEFT,XRIGH,YBOTT,YTTOP,1) CALL FRSTPT(XMIN,YMIN) CALL VECTOR(XMAX,YMIN) CALL VECTOR(XMAX,YMAX) CALL VECTOR(XMIN,YMAX) CALL VECTOR(XMIN,YMIN) C C I---------------------------------------------------I C I WRITE THE FRAME TITLE, CONTOUR INC. AND TIME I C I---------------------------------------------------I C C ****************************** size of titles and numbers ***** ISIZE=12 C IF(SZ.LE.0.5)ISIZE=0 XJ=XMIN YJ=YMAX+.05*YAX CALL PWRITX(XJ,YJ,TITLE,59,ISIZE,0,-1) C XJ=XMIN+0.5*XAX YJ=YMIN+.15*YAX IF(CIR.EQ.0.) GO TO 2011 ICIR=CIR WRITE(TCHAR,'(A4,I2,A3)') 'CI= ',ICIR,' SV' IF(CIR.LT.0.1.OR.CIR.GT.999.) 1 WRITE(TCHAR,'(A4,1PE7.0)') 'CI= ',CIR CALL PWRITX(XJ,YJ,TCHAR,11,ISIZE,0,0) C 2011 XJ=XMAX IF(TIME.LT.0.)GO TO 2002 WRITE(TCHAR,'(A6,F5.1,A5)') 'TIME= ',TIME,' DAYS' CALL PWRITX(XJ,YJ,TCHAR,16,ISIZE,0,1) C C I-----------------------------------------------------I C I COMPUTE THE X-COORDINATES OF THE X-AXIS LABELS I C I-----------------------------------------------------I C C 2002 YJ=YMIN-.04*YAX YYJMIN=YMIN+.02*YAX YYJMAX=YMAX-.02*YAX YYXLBL=YMIN-.1*YAX DO 10 I=0,NXINCS X=XMIN+I*DX IX=ABS(X*1.00001) C IF(MOD(IX,2).EQ.1) GO TO 2010 IF(I.NE.I/2*2) GO TO 2010 WRITE(XCHAR,'(I3,A1)') IX,' ' CALL PWRITX(X,YJ,XCHAR,4,ISIZE,0,0) 2010 CALL LINE(X,YMIN,X,YYJMIN) CALL LINE(X,YMAX,X,YYJMAX) 10 CONTINUE C I-------------------------I C I LABEL THE X-AXIS I C I-------------------------I CALL PWRITX(XCENT,YYXLBL,XLABEL,15,ISIZE,0,0) C C I-----------------------------------------------------I C I COMPUTE THE Y-COORDINATES OF THE Y-AXIS LABELS I C I-----------------------------------------------------I C XJ=XMIN-.02*XAX XXJMIN=XMIN+.02*XAX XXJMAX=XMAX-.02*XAX XXYLBL=XMIN-.15*XAX DO 20 J=0,NYINCS Y= YMIN+J*DY IY=ABS(Y*1.00001) C IF (MOD(IY,2).EQ.1) GO TO 2030 C WRITE(YCHAR,'(I4)') IY WRITE(YCHAR,'(F4.1)') Y CALL PWRITX(XJ,Y,YCHAR,4,ISIZE,0,1) 2030 CALL LINE(XMIN,Y,XXJMIN,Y) CALL LINE(XMAX,Y,XXJMAX,Y) 20 CONTINUE C I-------------------------I C I LABEL THE Y-AXIS I C I-------------------------I C CALL PWRITX(XXYLBL,YCENT,YLABEL,15,ISIZE,90,0) C I-------------------------------I C I SET FRAME FOR CONTOUR AND MAP I C I-------------------------------I ft=ft+SFACT/SQ rt=rt-SFACT/SQ bt=bt+SFACT tp=tp-SFACT CALL SET(ft,rt,bt,tp,XMIN,XMAX,YMIN,YMAX,1) CALL MAPPOS(ft,rt,bt,tp) 9998 RETURN END SUBROUTINE MAKCON(IND,IM,JM,X,Y,T,FSM,CMIN,CINT,CMAX,NU,SZ) 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 = 0; CONTOURS ARE NOT NUMBERED ; NU = 1; LABEL CONT. 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=1000) 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*32 IPAT,IPDS,IPDA 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= lenth/$ ; JSIZE= size of char. in numbered lines C WRITE(6,'('' CO ='',F10.3,'' L,LREM,M ='',3I10)') CO,L,LREM,M JCRT=8 JSIZE=2 C IF(SZ.LE.0.5)JCRT=3 C IF(SZ.LE.0.5)JSIZE=0 ICO=INT(CO) IPAT='$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$' IPDS='$''$''$''$''$''$''$''$''$''$''$''$''$''$''$''$''' C arrow-like C IPDA='$$$$''$''''$$$$''$''''$$$$''$''''$$$$''$''''' C IPDA='$$$$$$$$$$$$<$$$$$$$$$$$$$$$$$$$' IF (CO.LT.0.0.AND.IND.EQ.0)IPAT=IPDS C IF (CO.EQ.10.0)IPAT=IPDA IF(NU.EQ.1.AND.K.EQ.K/2*2)WRITE(IPAT,'(A28,I4)')IPAT,ICO C WRITE(6,'('' IPAT ='',A16)')IPAT C C high quality labeling CALL PCSETI ('QU',0) C CALL DASHDC(IPAT,JCRT,JSIZE) CALL CURVED(XCR,YCR,M) C reset solid line again IPAT='$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$' CALL DASHDC(IPAT,JCRT,JSIZE) C 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