program makeplotfiles01 c - Cannibalized from NADCON 5 for VERTCON 3 c - 2018 03 22 update: Still working on getting GMT 5 working, plus now c - have special files kicked out just for CONUS, which use these c - newly invented codes: c "v2" (in place of "cd") for "H88-H29 interpolated from Vertcon 2 grid", CONUS only c "rz" (in place of "cd") for residual between "cd" and "v2", CONUS only c - 2016 08 26 c - Added new code to do reference vectors consistently c - See DRU-12, p. 56-57 c - Also changing the call to "getmapbounds" to give it "olddatum" and "newdatum" c - to aide in filtering out things like the Saint regions in Alaska c - for unsupported transformations. c - 2016 07 29: c Scrapped the code for personalized reference vector location. Just put c all ref vectors outside/below plot. c - 2016 07 21: c - Added code to allow for optional placement of reference vectors, coming from c "map.parameters" as read in subroutine "getmapbounds" c c - Part of the NADCON5 process c c - Program presumes, as input c - the "old datum", "new datum" , "region" , "filter" and "mapflag" c - arguments fed into "doit2.bat" c - Examples: c olddatum = 'ngvd29' c newdatum = 'navd88' c region = 'conus' c filter = 1 (or "001" in file names) c mapflag = 0 c c - Program to take a "work" file and create c - a variety of GMT-ready data files of the following c - 3) Coverage in orthometric height c - 6) Vectors in orthometric height c - And these, for CONUS only: c 8) Vectors in orthometric height from VERTCON 2 grid c 9) Residuals between 6 and 8 c - It furthermore will create batch file to run the c - GMT scripts: c gmtbat01.(olddtm).(newdtm).(region).(filter).(mapflag) c --------------------------------------------------- implicit double precision(a-h,o-z) parameter(maxplots=60) c - Stuff that is in the work file records: character*6 pid character*2 state character*1 rejoht real*8 xlath,xlonh,xohth real*8 dohtm,v2dohtm,residual character*10 olddtm,newdtm,region c - Input work file: character*200 wfname c - Plot/Scale stuff real*8 lorvopc real*8 lorvogohtm,lorvogohtmv2,lorvogohtmrz c - Output GMT-ready files: c - Coverage, All, Coordinate Differences character*200 gfncvacdoht c - Vectors, All, Coordinate Differences, Meters character*200 gfnvmacdoht c - Vectors, All, VERTCON 2 Coordinate Differences, Meters, CONUS only character*200 gfnvmav2oht c - Vectors, All, Residuals from V2 and CD Coordinate Differences, Meters, CONUS only character*200 gfnvmarzoht c - GMT stuff real*8 bw(maxplots),be(maxplots),bs(maxplots),bn(maxplots) real*4 jm(maxplots) real*4 b1(maxplots),b2(maxplots) character*10 fn(maxplots) c - 2016 07 21: logical lrv(maxplots) real*8 rv0x(maxplots),rv0y(maxplots),rl0y(maxplots) c - Output GMT-batch file: character*200 gmtfile character*1 mapflag character*3 filter character*200 suffix1,suffix4 character*1 h1f1,h1f2,h1f3,h1f4,h2f1,h2f2,h2f3,h2f4 logical v2special c ------------------------------------------------------------------ c - BEGIN PROGRAM c ------------------------------------------------------------------ write(6,1001) 1001 format('BEGIN makeplotfiles01.f') c-------------------------------------------- c - Some necessary constants c c - "lorvopc" = Length of Reference Vector on Paper in **CM** c - *** CRITICAL: The units of centimeters align with c - the use, in .gmtdefaults4, of: c - MEASURE_UNIT = cm c - If you change that value, then the units of lorvopc c - will be whatever you change it too, which may prove c - displeasing. c - c - "csm" is the Color Sigma Multiplier...the c - number of Standard Deviations away from c - the average for the color palette to c - span on color plots. c-------------------------------------------- lorvopc = 1.d0 pi = 2.d0*dasin(1.d0) d2r = pi/180.d0 re = 6371000.d0 csm = 3.d0 MultiplierLorvog = 2 c-------------------------------------------- c - User-supplied input c-------------------------------------------- read(5,'(a)')olddtm read(5,'(a)')newdtm read(5,'(a)')region read(5,'(a)')filter read(5,'(a)')mapflag v2special = .false. if(trim(region).eq.'conus' .and. * trim(olddtm).eq.'ngvd29' .and. * trim(newdtm).eq.'navd88') v2special = .true. c ------------------------------------------------------------------ c - Generate the suffixes used in all our files c ------------------------------------------------------------------ suffix1=trim(olddtm)//'.'//trim(newdtm)//'.'//trim(region) *//'.'//trim(filter) suffix4=trim(suffix1)//'.'//trim(mapflag) c ------------------------------------------------------------------ c - Open the work file c ------------------------------------------------------------------ wfname='work.'//trim(suffix1) open(1,file='Work/'//wfname,status='old',form='formatted') write(6,1004)trim(wfname) 1004 format(6x,'makeplotfiles01.f: Opening work file ',a) c ------------------------------------------------------------------ c - Open the GMT batch file for plotting c ------------------------------------------------------------------ gmtfile = 'gmtbat01.'//trim(suffix4) open(99,file=gmtfile,status='new',form='formatted') write(6,1011)trim(gmtfile) 1011 format(6x,'makeplotfiles01.f: Creating GMT batch file ',a) write(99,1030)trim(gmtfile) 1030 format('echo BEGIN batch file ',a) c ------------------------------------------------------------------ c - Open the output GMT-ready files c ------------------------------------------------------------------ gfncvacdoht = 'cvacdoht.'//trim(suffix1) open(13,file=gfncvacdoht,status='new',form='formatted') write(6,1010)trim(gfncvacdoht) c - Vectors, Meters, All, Coordinate Differences gfnvmacdoht = 'vmacdoht.'//trim(suffix1) open(23,file=gfnvmacdoht,status='new',form='formatted') write(6,1010)trim(gfnvmacdoht) c - NGVD29/NAVD88/CONUS only begin ------------------------------------------------------- if(v2special)then c - Vectors, Meters, All, Coordinate Differences from VERTCON 2 grid gfnvmav2oht = 'vmav2oht.'//trim(suffix1) open(25,file=gfnvmav2oht,status='new',form='formatted') write(6,1010)trim(gfnvmav2oht) c - Vectors, Meters, All, Residual Coordinate Differences between "cd" and "v2" gfnvmarzoht = 'vmarzoht.'//trim(suffix1) open(26,file=gfnvmarzoht,status='new',form='formatted') write(6,1010)trim(gfnvmarzoht) endif c - CONUS only end ------------------------------------------------------- 1010 format(6x,'makeplotfiles01.f: Creating file ',a) c - LAZY CODING: Spin through work file, collecting stats c - which we will use to scale the map and vectors, etc. c - Then rewind the work file. c - If I wasn't lazy, this would be a RAM-read, etc etc ndoht = 0 avedohtm = 0.d0 rmsdohtm = 0.d0 stddohtm = 0.d0 nv2doht = 0 avev2dohtm = 0.d0 rmsv2dohtm = 0.d0 stdv2dohtm = 0.d0 nrzdoht = 0 averzdohtm = 0.d0 rmsrzdohtm = 0.d0 stdrzdohtm = 0.d0 891 if(.not.v2special)then read(1,104,end=892)pid,state,rejoht, * xlath,xlonh,xohth,dohtm,olddtm,newdtm, * h1f1,h1f2,h1f3,h1f4,h2f1,h2f2,h2f3,h2f4 else read(1,2104,end=892)pid,state,rejoht, * xlath,xlonh,xohth,dohtm,olddtm,newdtm, * h1f1,h1f2,h1f3,h1f4,h2f1,h2f2,h2f3,h2f4, * v2dohtm,residual endif if(rejoht.eq.' ')then c - Because they can be +/-, we and this is JUST for scaling the map, c - we want the average MAGNITUDE...use ABS. This means that c - I don't need to accept my own latter advice of returning c - here and using RMS. avedohtm = avedohtm + dabs(dohtm) rmsdohtm = rmsdohtm + dabs(dohtm**2) ndoht = ndoht + 1 if(v2special)then avev2dohtm = avev2dohtm + dabs(v2dohtm) rmsv2dohtm = rmsv2dohtm + dabs(v2dohtm**2) nv2doht = nv2doht + 1 averzdohtm = averzdohtm + dabs(residual) rmsrzdohtm = rmsrzdohtm + dabs(residual**2) nrzdoht = nrzdoht + 1 endif endif goto 891 892 continue if(ndoht.gt.0)then avedohtm = avedohtm / dble(ndoht) rmsdohtm = dsqrt(rmsdohtm / dble(ndoht)) fac = dble(ndoht) / dble(ndoht-1) stddohtm = dsqrt(fac*(rmsdohtm**2-avedohtm**2)) if(v2special)then avev2dohtm = avev2dohtm / dble(ndoht) rmsv2dohtm = dsqrt(rmsv2dohtm / dble(nv2doht)) fac = dble(nv2doht) / dble(nv2doht-1) stdv2dohtm = dsqrt(fac*(rmsv2dohtm**2-avev2dohtm**2)) averzdohtm = averzdohtm / dble(ndoht) rmsrzdohtm = dsqrt(rmsrzdohtm / dble(nrzdoht)) fac = dble(nrzdoht) / dble(nrzdoht-1) stdrzdohtm = dsqrt(fac*(rmsrzdohtm**2-averzdohtm**2)) endif else avedohtm = 0.d0 stddohtm = 0.d0 if(v2special)then avev2dohtm = 0.d0 stdv2dohtm = 0.d0 averzdohtm = 0.d0 stdrzdohtm = 0.d0 endif endif if(.not.v2special)then write(6,893)ndoht,avedohtm,stddohtm else write(6,2893)ndoht,avedohtm,stddohtm, * avev2dohtm,stdv2dohtm, * averzdohtm,stdrzdohtm endif 893 format(6x,'makeplotfiles01.f: Vector Stats: ',/, *10x,'Number of Good Ort. Ht. Vectors : ',i10,/, *10x,'Average length (meters) : ',f10.3,/, *10x,'STD length (meters) : ',f10.3) 2893 format(6x,'makeplotfiles01.f: Vector Stats: ',/, *10x,'Number of Good Ort. Ht. Vectors : ',i10,/, *10x,'Average length (meters) : ',f10.3,/, *10x,'STD length (meters) : ',f10.3,/, *10x,'Average length, Vertcon2 (meters) : ',f10.3,/, *10x,'STD length, Vertcon2 (meters) : ',f10.3,/, *10x,'Average length, Residual (meters) : ',f10.3,/, *10x,'STD length, Residual (meters) : ',f10.3) c - END OF LAZY CODING c -------------------------------------------------------------- c - GMT Batch file: Determine Number of areas to plot (=nplots) c - and Boundaries of plots and other stuff c -------------------------------------------------------------- write(6,1006)trim(region) 1006 format(6x,'makeplotfiles01.f: Calling getmapbounds for region ',a) c - 2016 08 29: call getmapbounds(mapflag,maxplots,region,nplots, *olddtm,newdtm, *bw,be,bs,bn,jm,b1,b2,fn,lrv,rv0x,rv0y,rl0y) c - 2016 08 26, DRU-12, p.56-57 c call getmapbounds(mapflag,maxplots,region,nplots, c *bw,be,bs,bn,jm,b1,b2,fn,lrv,rv0x,rv0y,rl0y) c - See DRU-11, p. 126 for choices on grid boundaries for NADCON v5.0 c - 2016 07 21: c call getmapbounds(mapflag,maxplots,region,nplots, c *bw,be,bs,bn,jm,b1,b2,fn,lrv,rv0x,rv0y) c call getmapbounds(mapflag,maxplots,region,nplots, c *bw,be,bs,bn,jm,b1,b2,fn) c --------------------------------------------------------------- c - Compute and report various things about our coverage c - and vector plots. c --------------------------------------------------------------- c - Reference vector set to length of the average c - absolute value of ort ht shift. c - I keep flip-flopping about whether it's more visually c - pleasing to use the average or twice the average. c - As of 2015 10 22, I go with "twice the average" c - Thus: MultiplierLorvog = 2 c - Note "gm2pc" = "ground meters to paper centimeters" c - "gs2pc" = "ground arcseconds to paper centimeters" c - "lorvog*s" = "length of reference vector on ground (element) in arcseconds" c - "lorvog*m" = "length of reference vector on ground (element) in meters" if(ndoht.ne.0)then lorvogohtm = onzd2(MultiplierLorvog*avedohtm) gm2pcoht = lorvopc / lorvogohtm c - VERTCON3, NGVD29/NAVD88/CONUS specials: if(v2special)then lorvogohtmv2 = onzd2(MultiplierLorvog*avev2dohtm) gm2pcohtv2 = lorvopc / lorvogohtmv2 lorvogohtmrz = onzd2(MultiplierLorvog*averzdohtm) gm2pcohtrz = lorvopc / lorvogohtmrz endif endif write(6,894)nplots 894 format(6x,'makeplotfiles01.f: Info about plots:',/, *8x,'Number of sub-region plot sets ', *'being made for this region: ',i2) do 895 i=1,nplots dns = bn(i) - bs(i) dew = be(i) - bw(i) write(6,896)i,region,fn(i), * bs(i),bn(i),bw(i),be(i),dns,dew if(ndoht.ne.0)then write(6,899)lorvogohtm,gm2pcoht if(v2special)then write(6,2899)lorvogohtmv2,gm2pcohtv2 write(6,3899)lorvogohtmrz,gm2pcohtrz endif else write(6,900) endif 895 continue 896 format(50('-'),/, *8x,'Plot # ',i2, *'(',a,1x,a,')',/, *10x,'S/N/W/E/N-S/E-W = ',6f7.1) 899 format( *10x,'Ort. Height plots: ',/, *12x,'Reference Vector ( m) = ',f10.2,/, *12x,'Ground M to Paper cm = ',f20.10) 2899 format( *10x,'Ort. Height plots (vertcon2 grid): ',/, *12x,'Reference Vector ( m) = ',f10.2,/, *12x,'Ground M to Paper cm = ',f20.10) 3899 format( *10x,'Ort. Height plots (residuals): ',/, *12x,'Reference Vector ( m) = ',f10.2,/, *12x,'Ground M to Paper cm = ',f20.10) 900 format( *10x,'Ort. Height plots: ',/, *12x,'No orthometric height data available for plotting') c - pvlon and pvlat are the percentage of the total lon/lat c - span of the plot, from which the reference vector c - begins, starting with the Lower Left corner pvlon = 10.d0 pvlat = 10.d0 pvlat = (pvlat / 100.d0) pvlon = (pvlon / 100.d0) c -------------------------------------------------------------- c -------------------------------------------------------------- c -------------------------------------------------------------- c - MAIN LOOP: TOP c -------------------------------------------------------------- c -------------------------------------------------------------- c -------------------------------------------------------------- c - Read in each record and produce GMT-ready data for c - all appropriate non-rejected records. rewind(1) n=0 ncvoht = 0 nvcoht = 0 1 if(.not.v2special)then read(1,104,end=2)pid,state,rejoht, * xlath,xlonh,xohth,dohtm,olddtm,newdtm, * h1f1,h1f2,h1f3,h1f4,h2f1,h2f2,h2f3,h2f4 else read(1,2104,end=2)pid,state,rejoht, * xlath,xlonh,xohth,dohtm,olddtm,newdtm, * h1f1,h1f2,h1f3,h1f4,h2f1,h2f2,h2f3,h2f4, * v2dohtm,residual endif c 104 format(a6,1x,a2,a1,a1,a1,1x,f14.10,1x,f14.10,1x,f8.3,1x, c *f9.5,1x,f9.5,1x,f9.3,1x,f9.5,1x,f9.5,1x,f9.3,1x,f9.3,1x,f9.3, c *1x,a10,1x,a10) 104 format(a6,1x,a2,a1,1x,f14.10,1x,f14.10,1x,f9.4,1x, * f9.4,1x,a10,1x,a10,1x,8(a1,1x)) 2104 format(a6,1x,a2,a1,1x,f14.10,1x,f14.10,1x,f9.4,1x, * f9.4,1x,a10,1x,a10,1x,8(a1,1x),1x,f9.4,1x,f9.4) c - Experiment 105 format(f16.10,1x,f15.10,1x,f6.2,1x,f12.2,1x,f5.1) 1105 format(f16.10,1x,f15.10,1x,f6.2,1x,a6) c1106 format(f16.10,1x,f15.10,1x,f6.2,1x,f12.2,1x,f5.1,1x,f9.5,1x,a6) c1106 format(f16.10,1x,f15.10,1x,f6.2,1x,f12.2,1x,f9.5,1x,f9.3,1x,a6) c - New for VERTCON 3, not VECTORS but dots of COLOR c1106 format(f16.10,1x,f15.10,1x,6x ,1x,f12.2,1x,9x ,1x,f9.3,1x,a6) 1106 format(f16.10,1x,f15.10,1x,6x ,1x,12x ,1x,9x ,1x,f9.3,1x,a6) c--------------------------------------------------------------------- c - If it's a good ort ht point, write to: c 13) Ort Ht Coverage File c 23) Ort Ht Vector (meters) File c 25) Ort Ht Vector from VERTCON 2 grid (meters) File c 26) Ort Ht Vector residual between true diffs and vertcon 2 diffs, meters, file if(rejoht.eq.' ')then if(dohtm.le.0)then az = 180.d0 else az = 0.d0 endif vcohtm = dabs(dohtm*gm2pcoht) write(13,1105)xlonh,xlath,sngl(1.0),pid c write(23,1106)xlonh,xlath,az,vcohtm,0.d0,dohtm,pid c - NEW for VERTCON 3: c - Not VECTORS, but dots of COLOR c write(23,1106)xlonh,xlath,vcohtm,dohtm,pid write(23,1106)xlonh,xlath,dohtm,pid if(v2special)then write(25,1106)xlonh,xlath,v2dohtm,pid write(26,1106)xlonh,xlath,residual,pid endif ncvoht = ncvoht + 1 nvcoht = nvcoht + 1 endif c - Count total number of points read in the work file n=n+1 goto 1 c -------------------------------------------------------------- c -------------------------------------------------------------- c -------------------------------------------------------------- c - MAIN LOOP: BOTTOM c -------------------------------------------------------------- c -------------------------------------------------------------- c -------------------------------------------------------------- 2 continue write(6,778)n,ncvoht,nvcoht 778 format(6x,'makeplotfiles01.f: Statistics: ',/, *10x,'Number of total records read : ',i10,/, *10x,'Number of oht coverage records prepared: ',i10,/, *10x,'Number of oht vector records prepared: ',i10) c -------------------------------------------------------------- c -------------------------------------------------------------- c -------------------------------------------------------------- c - CREATE GMT BATCH FILE c - There are multiple steps here: c - 1) Determine how many maps to make (for example, for c CONUS its a little easier to look at things when c broken down into a 3x3 tile over CONUS, rather than c one large plot) c - 2) Determine the boundaries of each map to be created c - 3) Determine which shoreline (GMT-provided, or Dru's c personally created detailed shorelines of certain c areas) c - 4) The title of each plot c - 5) The size of each plot c - 6) The file name to use for each plot c - c - As each decision is made, make the right calls to the batch c - file. c -------------------------------------------------------------- c -------------------------------------------------------------- c -------------------------------------------------------------- write(6,800) 800 format(6x,'makeplotfiles01.f: Preparing GMT batch file') c -------------------------------------------------------------- c - GMT Batch file: Make header c -------------------------------------------------------------- write(6,801) 801 format(6x,'makeplotfiles01.f: GMT: Write header') c write(99,901) c 901 format('gmtset GRID_PEN_PRIMARY 0.25p,-') write(99,902) c 902 format('gmtset BASEMAP_TYPE fancy') c 902 format('gmtset BASEMAP_TYPE fancy',/, c *'gmtset HEADER_FONT Helvetica',/, c *'gmtset HEADER_FONT_SIZE 12p',/, c *'gmtset HEADER_OFFSET 0.5c') 902 format('gmtset MAP_GRID_PEN_PRIMARY 0.25p,-',/, *'gmtset MAP_FRAME_TYPE fancy',/, *'gmtset FONT_TITLE 12p,Helvetica',/, *'gmtset MAP_TITLE_OFFSET 0.5c') c *'gmtset GMT_COMPATIBILITY = 4') write(6,802) 802 format(6x,'makeplotfiles01.f: GMT: Num Plots Analysis') write(6,1005)trim(region) 1005 format(6x,'makeplotfiles01.f: REGION = ',a) c -------------------------------------------------------------- c -------------------------------------------------------------- c -------------------------------------------------------------- c - GMT Batch creation loop TOP c -------------------------------------------------------------- c -------------------------------------------------------------- c -------------------------------------------------------------- c - For each "ij" value, we will put GMT calls into batch c - file "gmtbat01***" (unit 99) which will plot the following: c - 3) Coverage in OHT c - 6) Vectors in OHT, meters c - For CONUS specifically, and for NGVD29/NAVD88, do this: c -------------------------------------------------------------- c c - Note on "igridsec": Later in the NADCON5 processing, c - "igridsec" will mean the number of arcseconds at which c - thinning and gridding occur. At this point, though, c - it has no meaning. We set it to "-1" to tell the c - plotting programs to ignore it. igridsec = -1 igridsec2 = -1 do 2001 ij=1,nplots write(99,990)trim(region),trim(fn(ij)), * trim(region),trim(fn(ij)) c-------------------------------------- c - B/W COVERAGE PLOTS c-------------------------------------- c - Make B/W Coverage Plots c - Orthometric Height c - VERTCON3: Added "filter" to call call bwplotcv('oht',gfncvacdoht,bw,be,bs,bn,jm, * b1,b2,maxplots,olddtm,newdtm,region,filter,'OHT',ij, * igridsec,fn) c-------------------------------------- c - VECTORS (not really - color dots...so had to hack some c of the code from "makeplotfiles02.f" from NADCON 5 c into here) c-------------------------------------- csm = 3.0d0 c - Where is the start of the reference vector: c - 2016 07 29: Decided to scrap personalized locations c - and just put all reference vectors outside/below c - the plots c - 2016 08 26 xvlon = rv0x(ij) xvlat = rv0y(ij) xllon = xvlon xllat = rl0y(ij) c - Vectors in OHT c - Not VECTORS anymore, but color DOTS. c call bwplotvc('oht',gfnvmacdoht,bw,be,bs,bn,jm,b1,b2,maxplots, c * olddtm,newdtm,region,'OHT',ij,xvlon,xvlat,xllon,xllat, c * lorvogohtm,lorvopc,igridsec,fn) c - Always do the Hnew-Hold plot... call cpt(avedohtm,stddohtm,csm,cptlocd,cpthicd,cptincd) c - VERTCON3: Added "filter" to call call coplotvc('oht',gfnvmacdoht,bw,be,bs,bn,jm,b1,b2,maxplots, * olddtm,newdtm,region,filter,'OHT',ij,xvlon,xvlat,xllon,xllat, * lorvogohtm,lorvopc,igridsec,igridsec2,fn, * cptlocd,cpthicd,cptincd) c - But for the v2special, we have 2 more to do... if(v2special)then call cpt(avev2dohtm,stdv2dohtm,csm,cptlov2,cpthiv2,cptinv2) c - VERTCON3: Added "filter" to call call coplotvc('oht',gfnvmav2oht,bw,be,bs,bn,jm,b1,b2,maxplots, * olddtm,newdtm,region,filter,'OHT',ij,xvlon,xvlat,xllon,xllat, * lorvogohtmv2,lorvopc,igridsec,igridsec2,fn, * cptlov2,cpthiv2,cptinv2) call cpt(averzdohtm,stdrzdohtm,csm,cptlorz,cpthirz,cptinrz) c - VERTCON3: Added "filter" to call call coplotvc('oht',gfnvmarzoht,bw,be,bs,bn,jm,b1,b2,maxplots, * olddtm,newdtm,region,filter,'OHT',ij,xvlon,xvlat,xllon,xllat, * lorvogohtmrz,lorvopc,igridsec,igridsec2,fn, * cptlorz,cpthirz,cptinrz) endif 2001 continue c -------------------------------------------------------------- c -------------------------------------------------------------- c -------------------------------------------------------------- c - GMT Batch creation loop BOTTOM c -------------------------------------------------------------- c -------------------------------------------------------------- c -------------------------------------------------------------- write(99,1031)trim(gmtfile) 1031 format('echo END batch file ',a) close(99) write(6,9999) 9999 format('END makeplotfiles01.f') 990 format( *'# ------------------------------',/, *'# Plots for region: ',a,', sub-region: ',a,/, *'# ------------------------------',/, *'echo Creating plots for region: ',a,', sub-region: ',a) end c c c include 'Subs/getmapbounds.f' include 'Subs/plotcoast.f' include 'Subs/bwplotcv.f' c include 'Subs/bwplotvc.f' include 'Subs/coplotvc.f' c include '/home/dru/Subs/onzd2.f' include 'Subs/cpt.f'