C**** VECNP4.FOR Northern Polar Vector PostScript program 99/02/25 C**** C**** Compile: FCE90 VECNP4.FOR C**** C**** Compilation notes: The directory /u/cmrun/ should contain C**** the files: Z72X46N , Z144X90N , Z8XT70 , and VECNP.I . C**** This directory can be changed in this Fortran source listing. C**** The executable module a.out should be renamed to VECNP . C**** The PP command sends the PostScript file to the 6-th floor C**** network printer. The printer's correct address should be part C**** of the argument to the CALL SYSTEM Fortran statment. C**** PARAMETER (MAXREC=128, MAXIJR=72*46*MAXREC) C**** C**** MAXREC = dimension of TITLE array C**** MAXIJR = dimension of XYZD array C**** LOGICAL*4 QOPEN,QFOCEN,QFOPEN,QPOPEN COMMON /PARMCB/ IM,JM,NRECM(4),QOPEN(3),QFOCEN DATA QFOPEN/.FALSE./, QPOPEN/.FALSE./ C**** C**** IM,JM = longitude and latitude horizontal resolution C**** NRECM = number of records read from data file so far C**** NRECM(4) = number of difference records calculated so far C**** QOPEN = whether data file is currently open C**** QFOCEN = whether ocean fraction record has been read in C**** QFOPEN = whether VECNP.O and VECNP.IF have been opened C**** CHARACTER TITLE*80, JGRID*1,JLAND*1,JCONT*1, LABEL*16 COMMON /DATAFI/ UXYZD(MAXIJR,4),VXYZD(MAXIJR,4),TITLE(MAXREC,4), * JGRID(MAXREC,4),JLAND(MAXREC,4),JCONT(MAXREC,4), * JROTA(MAXREC,4),VALUE(MAXREC,4),LABEL(MAXREC,4) C**** C**** UXYZD = U component of longitude by latitude data read from files C**** VXYZD = V component of longitude by latitude data read from files C**** TITLE = title describing each record C**** JGRID = location of components on grid: A, B or C C**** JLAND = boxes shown: A = All, L = land, O = ocean C**** JCONT = continents: B = black outline, G = land is grey C**** JROTA = plot is rotated JROTA boxes to the left C**** VALUE = value of medium size vector C**** LABEL = label describing medium size vector in lower left corner C**** INTEGER*4 NREC(4), KMIN(3),KMAX(3) CHARACTER FILEIN*80, ARG*80, ERRARG*80, * UNNN*4,UNN2*4,UNN3*4, CDSCL*6 COMMON /FIXDCB/ FOCEAN(144,90) EQUIVALENCE (NRX,NREC(1)),(NRY,NREC(2)),(NRZ,NREC(3)) DATA NREC/4*1/ C**** C**** FOCEAN = ocean fraction for each horizontal grid box C**** NREC = most recently displayed record of file C**** IARGS = IARGC() IF(IARGS.le.0) GO TO 800 C**** Read in VECNP.I CALL MATCH0 C**** C**** Process command line arguments C**** IARG = 1 CALL GETARG (IARG,ARG) C**** Determine resolution LCS = LOCACS ('X',ARG,80) IF(LCS.le.0) LCS = LOCACS ('x',ARG,80) IF(LCS.le.0) GO TO 120 READ (ARG(1:LCS-1) ,*,IOSTAT=IOSQ) IM IF(IOSQ.ne.0) GO TO 120 READ (ARG(LCS+1:80),*,IOSTAT=IOSQ) JM IF(IOSQ.ne.0) GO TO 110 IARG = 2 GO TO 120 110 IM = 72 C**** Open X file 120 CALL GETARG (IARG,ARG) OPEN (1,FILE=ARG,FORM='UNFORMATTED',STATUS='OLD',ERR=810) QOPEN(1) = .TRUE. NRX = 1 CALL READDF (1,NRX,*820) CALL MATCHI (IM*JM,UXYZD(1,1),VXYZD(1,1),TITLE(1,1),JGRID(1,1), * JLAND(1,1),JCONT(1,1),JROTA(1,1),VALUE(1,1),LABEL(1,1)) IARG = IARG + 1 IF(IARG.gt.IARGS) GO TO 200 C**** Open Y file CALL GETARG (IARG,ARG) OPEN (2,FILE=ARG,FORM='UNFORMATTED',STATUS='OLD',ERR=810) QOPEN(2) = .TRUE. NRY = 1 CALL READDF (2,NRY,*820) CALL MATCHI (IM*JM,UXYZD(1,2),VXYZD(1,2),TITLE(1,2),JGRID(1,2), * JLAND(1,2),JCONT(1,2),JROTA(1,2),VALUE(1,2),LABEL(1,2)) IARG = IARG + 1 IF(IARG.gt.IARGS) GO TO 200 C**** Open Z file CALL GETARG (IARG,ARG) OPEN (3,FILE=ARG,FORM='UNFORMATTED',STATUS='OLD',ERR=810) QOPEN(3) = .TRUE. NRZ = 1 CALL READDF (3,NRZ,*820) CALL MATCHI (IM*JM,UXYZD(1,3),VXYZD(1,3),TITLE(1,3),JGRID(1,3), * JLAND(1,3),JCONT(1,3),JROTA(1,3),VALUE(1,3),LABEL(1,3)) C**** C**** Read ocean fraction file C**** C**** Read in horizontal ocean coverage 200 QFOCEN = .FALSE. IF(IM.ne.72 .or. JM.ne.46) GO TO 210 ARG = '/u/cmrun/Z72X46N' GO TO 230 210 IF(IM.ne.144 .or. JM.ne.90) GO TO 220 ARG = '/u/cmrun/Z144X90N' GO TO 230 220 IF(IM.ne.36 .or. JM.ne.24) GO TO 240 ARG = '/u/cmrun/Z8XT70' 230 QFOCEN = .TRUE. OPEN (26,FILE=ARG,FORM='UNFORMATTED',STATUS='OLD',ERR=830) READ (26) ARG,(FOCEAN(I,1),I=1,IM*JM) CLOSE (26) C**** 240 NR = NRX IU = 1 WRITE (UNNN,923) 1000+NR UNNN(1:1) = CHAR (Z'57'+IU) C**** C**** Display color plot C**** 300 IMJMNR = 1 + IM*JM*(NR-1) CALL CPSIJ (IM,JM,FOCEAN,UXYZD(IMJMNR,IU),VXYZD(IMJMNR,IU), * TITLE(NR,IU),JGRID(NR,IU),JLAND(NR,IU),VALUE(NR,IU),UNNN) C**** C**** Receive keyboard input C**** 400 READ (5,940,END=750) ARG ERRARG = ARG IF(ARG(1:1).EQ.'Q' .OR. ARG(1:1).EQ.'q') GO TO 750 IF(ARG(1:1).EQ.'A' .OR. ARG(1:1).EQ.'a') GO TO 410 IF(ARG(1:1).EQ.'L' .OR. ARG(1:1).EQ.'l') GO TO 420 IF(ARG(1:1).EQ.'O' .OR. ARG(1:1).EQ.'o') GO TO 430 IF(ARG(1:1).EQ.'W' .OR. ARG(1:1).EQ.'w') GO TO 440 IF(ARG(1:1).EQ.'B' .OR. ARG(1:1).EQ.'b') GO TO 450 IF(ARG(1:1).EQ.'G' .OR. ARG(1:1).EQ.'g') GO TO 460 IF(ARG(1:1).EQ.'R' .OR. ARG(1:1).EQ.'r') GO TO 470 IF(ARG(1:1).EQ.'V' .OR. ARG(1:1).EQ.'v') GO TO 480 IF(ARG(1:1).EQ.'S' .OR. ARG(1:1).EQ.'s') GO TO 490 IF(ARG(1:1).EQ.'M' .OR. ARG(1:1).EQ.'m') GO TO 500 IF(ARG(1:1).EQ.'X' .OR. ARG(1:1).EQ.'x') GO TO 510 IF(ARG(1:1).EQ.'Y' .OR. ARG(1:1).EQ.'y') GO TO 510 IF(ARG(1:1).EQ.'Z' .OR. ARG(1:1).EQ.'z') GO TO 510 IF(ARG(1:1).EQ.'D' .OR. ARG(1:1).EQ.'d') GO TO 550 IF(ARG(1:1).EQ.'T' .OR. ARG(1:1).EQ.'t') GO TO 590 IF(ARG(1:1).EQ.'N' .OR. ARG(1:1).EQ.'n') GO TO 600 IF(ARG(1:1).EQ.'F' .OR. ARG(1:1).EQ.'f') GO TO 640 IF(ARG(1:1).EQ.'P' .OR. ARG(1:1).EQ.'p') GO TO 650 GO TO 400 C**** A: display All grid boxes 410 IF(QFOCEN) JLAND(NR,IU) = 'A' GO TO 300 C**** L: display only Land grid boxes 420 IF(QFOCEN) JLAND(NR,IU) = 'L' GO TO 300 C**** O: display only Ocean grid boxes 430 JLAND(NR,IU) = 'O' GO TO 300 C**** W: display grey continental outline over land boxes 440 JCONT(NR,IU) = 'G' GO TO 300 C**** B: display continental outline or unused boxes in Black 450 JCONT(NR,IU) = 'B' GO TO 300 C**** G: display data on different Grid 460 IF(ARG(2:2).EQ.'A' .OR. ARG(2:2).EQ.'a') JGRID(NR,IU) = 'A' IF(ARG(2:2).EQ.'B' .OR. ARG(2:2).EQ.'b') JGRID(NR,IU) = 'B' IF(ARG(2:2).EQ.'C' .OR. ARG(2:2).EQ.'c') JGRID(NR,IU) = 'C' GO TO 300 C**** R: Rotate color plot I grid boxes to the left 470 READ (ARG(2:10),*,IOSTAT=IOSQ) I IF(IOSQ.eq.0) JROTA(NR,IU) = MOD(I+JROTA(NR,IU)+1000*IM,IM) GO TO 300 C**** V: use new Value for medium size arrow 480 READ (ARG(2:20),*,IOSTAT=IOSQ) V IF(IOSQ.eq.0) VALUE(NR,IU) = V GO TO 300 C**** S: use new label or Scale in lower left corner of plot 490 LABEL(NR,IU) = ARG(2:17) GO TO 400 C**** M: Multiply current data 500 RMULK = 0. READ (ARG(2:80),*,IOSTAT=IOSQ) RMULK IF(RMULK.eq.0.) GO TO 400 VALUE(NR,IU) = VALUE(NR,IU)*RMULK WRITE (LABEL(NR,IU)(1:6),950) VALUE(NR,IU) DO 501 I=1+IM*JM*(NR-1),IM*JM*NR UXYZD(I,IU) = UXYZD(I,IU)*RMULK 501 VXYZD(I,IU) = VXYZD(I,IU)*RMULK GO TO 300 C**** Display another record 510 CALL NEWREC (ARG,NREC,IU,NR,UNNN,*400) NREC(IU) = NR GO TO 300 C**** C**** D: display a difference record C**** 550 CALL PARSE (ARG,80,3,KMIN,KMAX) IF(KMAX(3).gt.0) GO TO 551 C**** Display previously calculated differences CALL NEWREC (ARG,NREC,IU,NR,UNNN,*400) NREC(4) = NR GO TO 300 C**** Calculate and display difference C Check that new difference does not exceed limits of MAXREC, MAXIJR 551 CALL NEWREC (ARG(KMIN(2):80),NREC,KU2,KN2,UNN2,*400) CALL NEWREC (ARG(KMIN(3):80),NREC,KU3,KN3,UNN3,*400) IU = 4 NR = NRECM(4) + 1 NRECM(4) = NR NREC(4) = NR DO 552 K=1,72 552 IF(TITLE(KN2,KU2)(K:K) .ne. TITLE(KN3,KU3)(K:K)) GO TO 553 TITLE(NR,4) = TITLE(KN2,KU2)(1:76) // ' dif' GO TO 554 553 TITLE(NR,4) = TITLE(KN2,KU2)(1:72) //' - '// TITLE(KN3,KU3)(K:K+3) 554 DO 555 I=1,IM*JM UXYZD(I+IM*JM*(NR-1),4) = UXYZD(I+IM*JM*(KN2-1),KU2) * - UXYZD(I+IM*JM*(KN3-1),KU3) 555 VXYZD(I+IM*JM*(NR-1),4) = VXYZD(I+IM*JM*(KN2-1),KU2) * - VXYZD(I+IM*JM*(KN3-1),KU3) JGRID(NR,4) = JGRID(KN2,KU2) JLAND(NR,4) = JLAND(KN2,KU2) JCONT(NR,4) = JCONT(KN2,KU2) JROTA(NR,4) = JROTA(KN2,KU2) VALUE(NR,4) = VALUE(KN2,KU2) LABEL(NR,4) = LABEL(KN2,KU2) WRITE (UNNN,923) 1000+NR UNNN(1:1) = 'D' WRITE (6,930) UNNN // ' = ' // UNN2 // ' - ' // UNN3 GO TO 300 C**** C**** T: edit current title C**** 590 READ (5,940,END=750) ARG DO 591 K=1,80 591 IF(ARG(K:K).ne.'"') TITLE(NR,IU)(K:K) = ARG(K:K) GO TO 300 C**** C**** N: load New data file C**** 600 KU = 0 IF(ARG(2:2).EQ.'X' .OR. ARG(2:2).EQ.'x') KU = 1 IF(ARG(2:2).EQ.'Y' .OR. ARG(2:2).EQ.'y') KU = 2 IF(ARG(2:2).EQ.'Z' .OR. ARG(2:2).EQ.'z') KU = 3 IF(KU.le.0) GO TO 400 IF(QOPEN(KU)) CLOSE (KU) QOPEN(KU) = .FALSE. CALL PARSE (ARG(3:80),78,2,KMIN,KMAX) CALL FILEEX (ARG(2+KMIN(1):2+KMAX(1)),KMAX(1)-KMIN(1)+1,FILEIN) OPEN (KU,FILE=FILEIN,FORM='UNFORMATTED',STATUS='OLD',ERR=603) QOPEN(KU) = .TRUE. DO 601 N=1,NRECM(KU) 601 VALUE(N,KU) = 0. NRECM(KU) = 0 KN = 1 IF(KMAX(2).gt.0) READ (ARG(2+KMIN(2):2+KMAX(2)),*,IOSTAT=IOSQ) KN 602 CALL READDF (KU,KN,*400) IMJMKN = 1 + IM*JM*(KN-1) CALL MATCHI (IM*JM,UXYZD(IMJMKN,KU),VXYZD(IMJMKN,KU),TITLE(KN,KU), * JGRID(KN,KU),JLAND(KN,KU),JCONT(KN,KU), * JROTA(KN,KU),VALUE(KN,KU),LABEL(KN,KU)) IU = KU NR = KN NREC(IU) = NR WRITE (UNNN,923) 1000+NR UNNN(1:1) = CHAR (Z'57'+IU) GO TO 300 603 WRITE (6,930) 'Unable to open file: ' // FILEIN(1:56) GO TO 400 C**** C**** F: write Files VECNP.O and VECNP.IF of title, data and parameters C**** 640 IF(QFOPEN) GO TO 641 OPEN (7,FILE='VECNP.O',FORM='UNFORMATTED') OPEN (8,FILE='VECNP.IF') 641 WRITE (7) TITLE(NR,IU),(UXYZD(I,IU),I=1+IM*JM*(NR-1),IM*JM*NR), * (VXYZD(I,IU),I=1+IM*JM*(NR-1),IM*JM*NR) WRITE (8,964) TITLE(NR,IU)(1:32),JGRID(NR,IU),JLAND(NR,IU), * JCONT(NR,IU),JROTA(NR,IU),VALUE(NR,IU),LABEL(NR,IU) WRITE (6,930) * 'Files VECNP.O and VECNP.IF written for record: ' // UNNN GO TO 400 C**** C**** P: Print record on color printer C**** 650 CALL PARSE (ARG(3:80),78,2,KMIN,KMAX) CALL NEWREC (ARG(2+KMIN(1):80),NREC,KU,KN,UNN2,*400) IF(.not.QPOPEN) OPEN (9,FILE='VECNP.PS') QPOPEN = .TRUE. WRITE (9,965) IMJMKN = 1 + IM*JM*(KN-1) CALL WRITPS (IM,JM,FOCEAN,UXYZD(IMJMKN,KU),VXYZD(IMJMKN,KU), * TITLE(KN,KU),JGRID(KN,KU),JLAND(KN,KU),JCONT(KN,KU), * JROTA(KN,KU),VALUE(KN,KU),LABEL(KN,KU),*653) WRITE (9,968) IF(ARG(2:2).eq.'P' .or. ARG(2:2).eq.'p') GO TO 652 WRITE (6,930) * 'PostScript file VECNP.PS written containing: ' // UNN2 GO TO 400 652 CLOSE (9) QPOPEN = .FALSE. CALL SYSTEM ('lp -d pub6 VECNP.PS') WRITE (6,930) 'PostScript file VECNP.PS written and printed ' // * 'containing: ' // UNN2 GO TO 400 653 WRITE (9,968) GO TO 400 C**** C**** Q: Quit VECNP program C**** 750 IF(.NOT.QFOPEN) GO TO 751 CLOSE (7) CLOSE (8) 751 GO TO 999 C**** 800 WRITE (0,*) 'Usage: VECNP [144x90] filex [filey filey] 99/02/25' GO TO 999 810 WRITE (6,*) ' Unable to open file: ',ARG(1:56) STOP 820 STOP 830 WRITE (0,*) ' File not found: ',ARG(1:16) STOP C**** 923 FORMAT (I4) 930 FORMAT (3X,A) 940 FORMAT (A) 950 FORMAT (F6.2) 964 FORMAT (A32,2X,3A2,I4,F10.2,2X,A16) 965 FORMAT ( * '%!PS VECNP.PS PostScript output file from VECNP 99/02/25'/ * ' gsave 306 396 translate') 968 FORMAT (/'grestore showpage'/) 999 END SUBROUTINE READDF (IU,NR,*) C**** C**** READDF reads data from unit IU up to record NR C**** C**** NRECM = number of records read from file so far C**** PARAMETER (MAXREC=128, MAXIJR=72*46*MAXREC) LOGICAL*4 QOPEN,QFOCEN CHARACTER TITLE*80 COMMON /PARMCB/ IM,JM,NRECM(4),QOPEN(3),QFOCEN COMMON /DATAFI/ UXYZD(MAXIJR,4),VXYZD(MAXIJR,4),TITLE(MAXREC,4) C**** IF(NR.le.NRECM(IU)) RETURN IF(.not.QOPEN(IU)) GO TO 800 IF(NR.gt.MAXREC) GO TO 810 IF(IM*JM*NR.gt.MAXIJR) GO TO 820 DO 10 N=NRECM(IU)+1,NR IMIN = IM*JM*(N-1) + 1 IMAX = IM*JM*N READ (IU,IOSTAT=IOSQ) TITLE(N,IU),(UXYZD(I,IU),I=IMIN,IMAX), * (VXYZD(I,IU),I=IMIN,IMAX) 10 IF(IOSQ.ne.0) GO TO 830 NRECM(IU) = NR RETURN C**** 800 WRITE (6,980) CHAR(IU+87),NRECM(IU) RETURN 1 810 WRITE (6,981) 'Requested record exceeds internal dimens' // * 'ion of the parameter MAXREC = 128.' RETURN 1 820 WRITE (6,981) 'Requested record exceeds internal dimens' // * 'ion of the parameter MAXIJR.' RETURN 1 830 WRITE (6,983) N,CHAR(IU+87) CLOSE(IU) NRECM(IU) = N-1 QOPEN(IU) = .FALSE. RETURN 1 C**** 980 FORMAT (' File',A2,' was closed.',I6,' records were read in.') 981 FORMAT (3X,A) 983 FORMAT (' Error encountered reading record',I4,' of file',A2, * '. File is closed.') END SUBROUTINE MATCH0 C**** C**** MATCH0 loads data from the file VECNP.I into CPIJoI common block C**** PARAMETER (MAXoI=1024) REAL*4 UXYZD(*),VXYZD(*) CHARACTER TITLE*80, JGRID*1,JLAND*1,JCONT*1, LABEL*16, * TITLEI*32, IGRID*1,ILAND*1,ICONT*1, LABEI*16 COMMON /CPIJoI/ TITLEI(MAXoI),IGRID(MAXoI),ILAND(MAXoI), * ICONT(MAXoI),IROTA(MAXoI),VALUI(MAXoI),LABEI(MAXoI) C**** C**** TITLEI= matches first 32 characters of data record TITLE C**** IGRID = location of components on grid: A, B or C C**** ILAND = boxes shown: A = All, L = land, O = ocean C**** ICONT = continents: B = black outline, G = land is grey C**** IROTA = plot is rotated IROTA boxes to the left C**** VALUI = value of medium size vector C**** LABEI = label describing medium size vector in lower left corner C**** DATA IIMAX /0/ C**** C**** IIMAX = number of lines read from the file VECNP.I C**** C**** Read in VECNP.I C**** OPEN (4,FILE='VECNP.I',STATUS='OLD',ERR=10) WRITE (0,*) 'File VECNP.I read in.' GO TO 20 10 OPEN (4,FILE='/u/cmrun/VECNP.I',STATUS='OLD',ERR=60) WRITE (0,*) 'File /u/cmrun/VECNP.I read in.' 20 DO 30 I=1,4 30 READ (4,904) DO 40 I=1,MAXoI READ (4,904,IOSTAT=IOSQ) TITLEI(I),IGRID(I),ILAND(I),ICONT(I), * IROTA(I),VALUI(I),LABEI(I) 40 IF(IOSQ.ne.0) GO TO 50 50 IIMAX = I-1 CLOSE (4) WRITE (0,*) 'File VECNP.I contains',IIMAX,' lines of data.' IF(IIMAX.le.0) STOP 'MATCH0 50' RETURN 60 WRITE (0,*) 'Neither VECNP.I nor /u/cmrun/VECNP.I was found.' RETURN C**** C**** ENTRY MATCHI (IMJM,UXYZD,VXYZD,TITLE,JGRID,JLAND,JCONT, * JROTA,VALUE,LABEL) C**** C**** MATCHI matches the data record TITLE with TITLEIs in the file C**** VECNP.I to get initial parameters for the vector plot. C**** If the data record TITLE does not match any VECNP.I TITLEI, C**** then the initial parameters are estimated. C**** IF(VALUE.ne.0.) RETURN DO 130 L=1,IIMAX DO 110 K=1,32 110 IF(TITLEI(L)(K:K).ne.TITLE(K:K) .and. TITLEI(L)(K:K).ne.'.') * GO TO 130 C**** TITLEI(L) matches TITLE JGRID = IGRID(L) JLAND = ILAND(L) JCONT = ICONT(L) JROTA = IROTA(L) VALUE = VALUI(L) LABEL = LABEI(L) RETURN 130 CONTINUE C**** C**** Estimate vector plot parameters C**** 200 SQMAX = 0. DO 210 IJ=1,IMJM IF(SQMAX.lt.UXYZD(IJ)**2+VXYZD(IJ)**2) * SQMAX = UXYZD(IJ)**2+VXYZD(IJ)**2 210 CONTINUE JGRID = 'A' JLAND = 'A' JCONT = 'G' JROTA = 0 VALUE = ROUND (.5*SQRT(SQMAX)) WRITE (LABEL,921) VALUE RETURN C**** 904 FORMAT (A32,2X,3A2,I4,F10.0,2X,A16) 921 FORMAT (F6.2) END SUBROUTINE CPSIJ (IM,JM,FOCEAN,UXYZD,VXYZD,TITLE, * JGRID,JLAND,VALUE,UNNN) C**** C**** CPSIJ writes maximum magnitude to unit 6 C**** REAL*4 FOCEAN(IM,JM),UXYZD(IM,JM),VXYZD(IM,JM) CHARACTER TITLE*80, JGRID*1,JLAND*1, UNNN*4 C**** Calculate maximum magnitude SQMAX = 0. DO 110 I=1,IM*JM IF((JLAND.eq.'O' .and. FOCEAN(I,1).le..5) .or. * (JLAND.eq.'L' .and. FOCEAN(I,1).gt..5)) GO TO 110 IF(SQMAX.lt.UXYZD(I,1)**2+VXYZD(I,1)**2) * SQMAX = UXYZD(I,1)**2+VXYZD(I,1)**2 110 CONTINUE C**** Display TITLE above the data WRITE (6,920) TITLE(1:79) 320 WRITE (6,932) SQRT(SQMAX),VALUE,UNNN RETURN C**** 920 FORMAT (A) 932 FORMAT (' Maximum magnitude:',F10.2, * ' Value of medium arrow: ',F8.2,4X,A4) END SUBROUTINE WRITPS (IM,JM,FOCEAN,UXYZD,VXYZD,TITLE, * JGRID,JLAND,JCONT,JROTA,VALUE,LABEL) C**** C**** Write PostScript file for record KN from unit KU. C**** C**** IM,JM = horizontal resolution C**** FOCEAN= ocean fration array C**** UXYZD = U component of velocity C**** VXYZD = V component of velocity C**** TITLE = title at top of plot C**** JGRID = location of velocity points on grid: A, B or C C**** JLAND = boxes shown: A = all, L = land, O = ocean C**** JCONT = continents: B = black outline, G = land boxes are grey C**** JROTA = plot is rotated JROTA boxes to the left C**** VALUE = value of medium sized arrow C**** LABEL = label of arrow in lower left corner of plot C**** PARAMETER (MAXC=63, MAXBLK=64, MAXCB=64, IBACKC=0, IFOREC=1, * TWOPI=6.2831853) REAL*4 FOCEAN(IM,JM),UXYZD(IM,JM),VXYZD(IM,JM), * UA(360,180),VA(360,180) INTEGER*4 ILAND(360), IANG(360),IAMP(360) CHARACTER TITLE*80, JGRID*1,JLAND*1,JCONT*1, LABEL*16, ARG*80 C**** C**** Copy velocity components to UA and VA. Interpolate if necessary. C**** IF(JGRID.eq.'A') GO TO 10 IF(JGRID.eq.'B') GO TO 20 IF(JGRID.eq.'C') GO TO 30 C**** Velocity points are located on A grid 10 DO 11 J=2,JM-1 DO 11 I=1,IM UA(I,J) = UXYZD(I,J) 11 VA(I,J) = VXYZD(I,J) UNP = 0. VNP = 0. DO 12 I=1,IM ANGLE = TWOPI*(I-.5)/IM UNP = UNP + UXYZD(I,JM)*COS(ANGLE) - VXYZD(I,JM)*SIN(ANGLE) 12 VNP = VNP + VXYZD(I,JM)*COS(ANGLE) + UXYZD(I,JM)*SIN(ANGLE) UNP = UNP/IM VNP = VNP/IM GO TO 100 C**** Velocity points are located on B grid 20 STOP 'WRITPS: B grid not programmed.' C**** Velocity points on C grid are interpolated to A grid 30 IM1=IM DO 31 J=2,JM-1 DO 31 I=1,IM UA(I,J) = .5*(UXYZD(IM1,J)+UXYZD(I,J)) VA(I,J) = .5*(VXYZD(I,J-1)+VXYZD(I,J)) 31 IM1=I C USP = 0. C VSP = 0. UNP = 0. VNP = 0. DO 32 I=1,IM ANGLE = TWOPI*(I-.5)/IM C USP = USP + VXYZD(I,1 )*SIN(ANGLE) C VSP = VSP + VXYZD(I,1 )*COS(ANGLE) UNP = UNP - VXYZD(I,JM-1)*SIN(ANGLE) 32 VNP = VNP + VXYZD(I,JM-1)*COS(ANGLE) C USP = USP*2./IM C VSP = VSP*2./IM UNP = UNP*2./IM VNP = VNP*2./IM C DO 33 I=1,IM C ANGLE = TWOPI*(I-.5)/IM C UA(I,1 ) = USP*COS(ANGLE) - VSP*SIN(ANGLE) C VA(I,1 ) = VSP*COS(ANGLE) + USP*SIN(ANGLE) C UA(I,JM) = UNP*COS(ANGLE) + VNP*SIN(ANGLE) C 33 VA(I,JM) = VNP*COS(ANGLE) - UNP*SIN(ANGLE) C**** C**** Define spatial parameters C**** 100 DLON = 360./IM DLAT = NINT(360./(JM-1)) / 2. FJPOLE = 90./DLAT - (JM/2-1) WRITE (9,900) IM,JM, DLON J40N = 1 + NINT(40./DLAT) + JM/2 SCGCM = 72.*3.75 / (JM-J40N+FJPOLE) C**** C**** Write title, labels and rectangle to PostScript file C**** C**** Display title at top of plot LENGTH = LEN_TRIM (TITLE) WRITE (9,910) TITLE(1:LENGTH) C**** Display NASA/GISS logo in lower right corner WRITE (9,911) C**** Display surrounding circle WRITE (9,912) C**** Display label for standard arrow in lower left corner WRITE (9,913) '= '//LABEL C**** Display standard arrow in lower left corner WRITE (9,914) SCGCM,SCGCM C**** C**** Continental outline C**** JROTB = JROTA + IM/2 IF(JROTB.ge.IM) JROTB = JROTB - IM C IF(JCONT.ne.'G') GO TO ??? C**** Fill in each continental grid box with light grey WRITE (9,930) SCGCM,SCGCM DO 320 J=J40N,JM-1 NLAND = 0 I = 1+JROTB DO 310 IR=1,IM IP1 = IR+1+JROTB IF(IP1.gt.IM) IP1 = IP1-IM IF(FOCEAN(I,J).gt..5) GO TO 310 NLAND = NLAND + 1 ILAND(NLAND) = IR 310 I = IP1 IF(NLAND.le.0) GO TO 320 FJ = JM-J+FJPOLE WRITE (9,931) J,FJ, (ILAND(N),N=1,NLAND) WRITE (9,932) NLAND 320 CONTINUE C**** C**** Write longitude by latitude vector plot to PostScript file C**** WRITE (9,940) C**** Latitudes from 40N to JM-3 DO 420 J=J40N,JM-3 DO 410 IR=1,IM I = IR+JROTB IF(I.gt.IM) I = I-IM AMP = SQRT(UA(I,J)**2+VA(I,J)**2) / VALUE ANG = 0. IF(AMP.gt.0.) ANG = ATAN2(VA(I,J),UA(I,J)) IAMP(IR) = NINT(100.*SQRT(AMP)) IF(IAMP(IR).gt.999) IAMP(IR) = 999 IANG(IR) = NINT(ANG*360./TWOPI) IF(IANG(IR).lt.0.) IANG(IR) = IANG(IR) + 360 410 CONTINUE FJ = JM-J-.5+FJPOLE RNSEW = SQRT(IM/(TWOPI*FJ)) WRITE (9,941) J,FJ,RNSEW, (IANG(IR),IAMP(IR),IR=1,IM) 420 WRITE (9,942) IM WRITE (9,*) ' grestore' C**** Latitudes JM-2 WRITE (9,943) DLON*2. J=JM-2 DO 430 IR=1,IM,2 I =IR+JROTB IP1=IR+JROTB+1 IF(I .gt.IM) I =I -IM IF(IP1.gt.IM) IP1=IP1-IM UAAVE = .5*(UA(I,J)+UA(IP1,J)) VAAVE = .5*(VA(I,J)+VA(IP1,J)) AMP = SQRT(UAAVE**2+VAAVE**2) / VALUE ANG = 0. IF(AMP.gt.0.) ANG = ATAN2(VAAVE,UAAVE) IAMP(IR) = NINT(100.*SQRT(AMP)) IF(IAMP(IR).gt.999) IAMP(IR) = 999 IANG(IR) = NINT(ANG*360./TWOPI) IF(IANG(IR).lt.0.) IANG(IR) = IANG(IR) + 360 430 CONTINUE FJ = JM-J-.5+FJPOLE RNSEW = SQRT(IM/(2.*TWOPI*FJ)) WRITE (9,941) J,FJ,RNSEW, (IANG(IR),IAMP(IR),IR=1,IM,2) WRITE (9,944) IM/2 C**** Latitude JM-1 WRITE (9,943) DLON*4. J=JM-1 DO 450 IR=1,IM,4 I =IR+JROTB IP1=IR+JROTB+1 IP2=IR+JROTB+2 IP3=IR+JROTB+3 IF(I .gt.IM) I =I -IM IF(IP1.gt.IM) IP1=IP1-IM IF(IP2.gt.IM) IP2=IP2-IM IF(IP3.gt.IM) IP3=IP3-IM UAAVE = .25*(UA(I,J)+UA(IP1,J)+UA(IP2,J)+UA(IP3,J)) VAAVE = .25*(VA(I,J)+VA(IP1,J)+VA(IP2,J)+VA(IP3,J)) AMP = SQRT(UAAVE**2+VAAVE**2) / VALUE ANG = 0. IF(AMP.gt.0.) ANG = ATAN2(VAAVE,UAAVE) IAMP(IR) = NINT(100.*SQRT(AMP)) IF(IAMP(IR).gt.999) IAMP(IR) = 999 IANG(IR) = NINT(ANG*360./TWOPI) IF(IANG(IR).lt.0.) IANG(IR) = IANG(IR) + 360 450 CONTINUE FJ = JM-J-.5+FJPOLE RNSEW = SQRT(IM/(4.*TWOPI*FJ)) WRITE (9,941) J,FJ,RNSEW, (IANG(IR),IAMP(IR),IR=1,IM,4) WRITE (9,944) IM/4 C**** North Pole arrow AMP = SQRT(UNP**2+VNP**2) / VALUE ANG = 0. IF(AMP.gt.0.) ANG = ATAN2(VNP,UNP) NPAMP = NINT(100.*SQRT(AMP)) IF(NPAMP.gt.999) NPAMP = 999 NPANG = NINT(ANG*360./TWOPI) IF(NPANG.lt.0.) NPANG = NPANG + 360 WRITE (9,947) JROTB*DLON,NPANG,NPAMP WRITE (9,*) ' grestore' RETURN C**** 900 FORMAT (/'%%% Parameters:' / * ' /IM',I4,' def /JM',I4,' def /DLON',F9.3,' def') 910 FORMAT (/'%%% Produce TITLE at top of plot' / * ' /Helvetica-Bold findfont 20 scalefont setfont' / * '(',A,')' / * ' gsave -270 298 moveto dup stringwidth pop' / * ' 540 exch div 1 scale show grestore') 911 FORMAT (/'%%% Draw NASA/GISS logo in lower right corner' / * ' /Helvetica-Bold findfont 12 scalefont setfont' / * ' 204 -284 moveto (NASA/GISS) show') 912 FORMAT (/'%%% Draw surrounding circle' / ' gsave ', * '2 setlinewidth 1 setlinecap newpath 0 272 moveto' / *' IM 1 sub {DLON rotate 0 272 lineto} repeat' / *' closepath stroke grestore') 913 FORMAT (/'%%% Draw label for standard arrow in lower left corner'/ * ' /Helvetica-Bold findfont 16 scalefont setfont' / * ' -248 -284 moveto (',A,') show') 914 FORMAT (/'%%% Standard arrow in unit box [-.5:.5, -.5:.5]' / *' /ARROW {newpath 0 .1875 moveto -.1875 0 lineto' / *13X,'-.0625 0 lineto -.0625 -.1875 lineto .0625 -.1875 lineto' / *13X,'.0625 0 lineto .1875 0 lineto closepath fill} def' / *' gsave -262 -280 translate ',2F9.3,' scale -90 rotate' / *' ARROW grestore') 930 FORMAT (/3X,2F9.3,' scale %%% Scale is now 1 unit per box' / *'%%% Fill in continental grid boxes with light gray' / *' .8 setgray' / *' /LAND {gsave DLON mul rotate newpath FJ 0 moveto' / *' -1 0 rlineto DLON neg rotate FJ 1 sub 0 lineto' / *' 1 0 rlineto closepath fill grestore} def') 931 FORMAT ( '/J',I4,' def /FJ',F9.3,' def' / (18I4)) 932 FORMAT ( I8,' {LAND} repeat') 940 FORMAT (/'%%% Draw arrows' / *' /AI {DLON neg rotate gsave FJ 0 translate .01 mul' / *' dup RNSEW mul exch RNSEW div scale rotate' / *' ARROW grestore} def' / *' 0 setgray gsave DLON 2 div rotate') 941 FORMAT ('/J',I4,' def /FJ',F9.3,' def /RNSEW',F9.3,' def' / * (18I4)) 942 FORMAT (I8,' {AI} repeat') 943 FORMAT (' gsave /DLON',F9.3,' def DLON 2 div rotate') 944 FORMAT (I8,' {AI} repeat grestore') 947 FORMAT ('/J JM def' / *' gsave',F9.2,' neg',I6,' add rotate'/ *I7,' .01 mul dup scale ARROW') END SUBROUTINE PARSE (S,LENS,NTOK,KMIN,KMAX) C**** C**** PARSE parses a character string S of length LENS into its first C**** NTOK tokens. KMIN(n) and KMAX(n) are the first and last C**** characters of S for the n-th token. C**** INTEGER*4 KMIN(*),KMAX(*) CHARACTER S*(*) N = 0 K = 0 10 IF(N.GE.NTOK) RETURN N = N + 1 C**** Skip over blank characters between tokens 20 IF(K.GE.LENS) GO TO 50 K = K + 1 IF(S(K:K).EQ.' ') GO TO 20 KMIN(N) = K C**** Contiguous nonblank characters are part of the token 30 IF(K.GE.LENS) GO TO 40 K = K + 1 IF(S(K:K).NE.' ') GO TO 30 KMAX(N) = K-1 GO TO 10 40 KMAX(N) = K N = N + 1 C**** Remaining tokens were not included in the string 50 DO 60 N0=N,NTOK KMIN(N0) = 0 60 KMAX(N0) = 0 RETURN END SUBROUTINE NEWREC (S,NREC,IU,NR,UNNN,*) C**** C**** NEWREC decodes the token S into a file and a record number, C**** and checks that the record is available C**** C**** Input: S = string containing X,Y,Z,D and/or record number C**** NREC = most recently dispaled records of files X, Y or Z C**** Output: IU = unit number 1,2,3,4 corresponding to X,Y,Z,D C**** NR = record numbr of unit IU determined from S and/or NREC C**** UNNN = 4 byte character string showing file and record numbr C**** PARAMETER (MAXREC=128, MAXIJR=72*46*MAXREC) INTEGER*4 NREC(4) CHARACTER S*(*), UNNN*4 CHARACTER TITLE*80, JGRID*1,JLAND*1,JCONT*1, LABEL*16 COMMON /PARMCB/ IM,JM,NRECM(4),QOPEN(3),QFOCEN COMMON /DATAFI/ UXYZD(MAXIJR,4),VXYZD(MAXIJR,4),TITLE(MAXREC,4), * JGRID(MAXREC,4),JLAND(MAXREC,4),JCONT(MAXREC,4), * JROTA(MAXREC,4),VALUE(MAXREC,4),LABEL(MAXREC,4) C**** IF(S(1:1).EQ.'D' .OR. S(1:1).EQ.'d') GO TO 10 C**** File is X, Y or Z KU = 0 IF(S(1:1).EQ.'X' .OR. S(1:1).EQ.'x') KU = 1 IF(S(1:1).EQ.'Y' .OR. S(1:1).EQ.'y') KU = 2 IF(S(1:1).EQ.'Z' .OR. S(1:1).EQ.'z') KU = 3 IF(KU.EQ.0) GO TO 810 IF(NRECM(KU).LE.0) GO TO 820 KN = NREC(KU) IOSQ = 0 IF(S(2:2).NE.' ') READ (S(2:80),*,IOSTAT=IOSQ) KN IF(IOSQ.NE.0) GO TO 830 CALL READDF (KU,KN,*800) IMJMKN = 1 + IM*JM*(KN-1) CALL MATCHI (IM*JM,UXYZD(IMJMKN,KU),VXYZD(IMJMKN,KU),TITLE(KN,KU), * JGRID(KN,KU),JLAND(KN,KU),JCONT(KN,KU), * JROTA(KN,KU),VALUE(KN,KU),LABEL(KN,KU)) IU = KU NR = KN WRITE (UNNN,900) 1000+NR UNNN(1:1) = CHAR(Z'57'+IU) RETURN C**** File is D 10 KN = NREC(4) IOSQ = 0 IF(S(2:2).NE.' ') READ (S(2:80),*,IOSTAT=IOSQ) KN IF(IOSQ.NE.0) GO TO 830 IF(KN.GT.NRECM(4)) GO TO 840 IU = 4 NR = KN WRITE (UNNN,900) 1000+NR UNNN(1:1) = 'D' RETURN C**** 800 RETURN 1 810 WRITE (6,981) 'File should be X, Y, Z or D. Not: ' // S(1:12) RETURN 1 820 WRITE (6,981) S(1:1) // ' file was never opened.' RETURN 1 830 WRITE (6,981) 'Integer record number should follow file: ' // * S(1:12) RETURN 1 840 WRITE (6,981) 'Requested difference record is not defined: ' // * S(1:12) RETURN 1 C**** 900 FORMAT (I4) 981 FORMAT (3X,A) END SUBROUTINE FILEEX (FIN,LENF,FEX) C**** C**** FILEEX expands a filename FIN by environment variables defined C**** in the user's .profile file. C**** CHARACTER FIN*80, FEX*80 IF(FIN(1:1).eq.'$') GO TO 10 FEX = FIN(1:LENF) RETURN C**** Leading $ was found, determine / or . 10 DO 20 K=3,LENF 20 IF(FIN(K:K).eq.'/' .or. FIN(K:K).eq.'.') GO TO 30 FEX = FIN(1:LENF) RETURN C**** Internal / or . was found 30 KSLASH = K CALL GETENV (FIN(2:KSLASH-1),FEX) IF(FEX(1:1).gt.' ') GO TO 40 FEX = FIN(1:LENF) RETURN C**** Insert the rest of FIN into the output filename FEX 40 LENE = LEN_TRIM (FEX) FEX(LENE+1:80) = FIN(KSLASH:LENF) RETURN END BLOCK DATA C**** C**** Default values for PARMCB C**** LOGICAL*4 QOPEN,QFOCEN COMMON /PARMCB/ IM,JM,NRECM(4),QOPEN(3),QFOCEN C**** C**** IM,JM = longitude and latitude horizontal resolution C**** NRECM = number of records read from data file so far C**** NRECM(4) = number of difference records calculated so far C**** QOPEN = whether data file is currently open C**** QFOCEN = whether ocean fraction record has been read in C**** C**** DATA IM,JM/72,46/, NRECM/4*0/, QOPEN/3*.FALSE./, QFOCEN/.FALSE./ END FUNCTION LOCACS (C,S,LMAX) C**** C**** Find the first LOCAtion of the Character C in the String S C**** CHARACTER S*256, C, CHAR0 CHAR0 = CHAR(0) LOCACS = 0 DO 10 L=1,LMAX IF(S(L:L).EQ.CHAR0) RETURN 10 IF(S(L:L).EQ.C) GO TO 20 RETURN 20 LOCACS = L RETURN END C FUNCTION LENOFS (S,KMAX) replaced by LEN_TRIM (S) C**** C**** LENOFS calculates the LENgth OF the character String S, or C**** the location of the last nonblank character in the string S C**** C CHARACTER S*(*) C DO 10 K=KMAX,1,-1 C 10 IF(S(K:K).NE.' ' .AND. ICHAR(S(K:K)).NE.0) GO TO 20 C 20 LENOFS = K C RETURN C END FUNCTION ROUND (X) C**** C**** ROUND returns a rounded value of the positive number X C**** IF(X.LT.9000.) GO TO 10 ROUND = 80000. IF(X.GE.65000.) RETURN ROUND = NINT(X*.0001)*10000 RETURN C**** 10 IF(X.LT.900.) GO TO 20 ROUND = 8000. IF(X.GE.6500.) RETURN ROUND = NINT(X*.001)*1000 RETURN C**** 20 IF(X.LT.90.) GO TO 30 ROUND = 800. IF(X.GE.650.) RETURN ROUND = NINT(X*.01)*100 RETURN C**** 30 IF(X.LT.9.) GO TO 40 ROUND = 80. IF(X.GE.65.) RETURN ROUND = NINT(X*.1)*10 RETURN C**** 40 IF(X.LT..9) GO TO 50 ROUND = 8. IF(X.GE.6.5) RETURN ROUND = NINT(X) RETURN C**** 50 IF(X.LT..09) GO TO 60 ROUND = .8 IF(X.GE..65) RETURN ROUND = NINT(X*10.)*.1 RETURN C**** 60 IF(X.LT..015) GO TO 70 ROUND = .08 IF(X.GE..065) RETURN ROUND = NINT(X*100.)*.01 RETURN C**** 70 ROUND = .01 RETURN END