C C $Id: cssex02.f,v 1.1 2007/09/24 23:25:40 haley Exp $ C PROGRAM CSSEX02 C C Define random points and functional values on the globe, C triangulate, interpolate to a uniform grid, then use Ezmap C and Conpack to draw a contour plot. C C C External for CPCLDM that draws the contours. C EXTERNAL DRAWCL C C Specify GKS output data. C PARAMETER (IERRF=6, LUNIT=2, IWTYPE=1, IWKID=1) C C Dimension the arrays for holding the points for drawing C circular arcs. C PARAMETER (NPARC=50) DIMENSION ARCLAT(NPARC),ARCLON(NPARC) C C Size of the interpolated output grid. C PARAMETER (NI=73, NJ=145) C C Size of workspaces, etc. for triangulation and interpolation. C PARAMETER (NMAX=500, NTMX=2*NMAX, NT6=6*NMAX, LWK=27*NMAX) C C Array storage for the triangulation, work space, and nodal C coordinates. Note that the real workspace must be double C precision. C DOUBLE PRECISION RWK(13*NMAX) REAL RLAT(NMAX), RLON(NMAX), FVAL(NMAX), + PLAT(NI), PLON(NJ) REAL FF(NI,NJ), ZDAT(NJ,NI) C C Storage for the triangulation, work space, and Voronoi vertices. C INTEGER LTRI(3,NTMX), IWK(LWK), NV(NMAX) C C Storage for circumcenters and circumcircle radii. C REAL RLATO(NTMX), RLONO(NTMX), RC(NTMX) C C Sizes of work arrays for Conpack and Ezmap. C PARAMETER (IRWRK=1000, IIWRK=2000, IAMASZ=10000) DIMENSION RWRK(IRWRK), IWRK(IIWRK), IAMA(IAMASZ) C C Generate a default set of nodes as latitudinal and longitudinal C coordinates (latitudes in the range -90. to 90. and longitudes C in the range -180. to 180). The input function is generated C using the local subroutines CSGENRS, CSGENI, and CSGENPNT. C N = NMAX CALL CSGENRS(N,RLAT,RLON) CALL CSGENI(5,10,-200.,500.) DO 10 I=1,N FVAL(I) = CSGENPNT(RLAT(I),RLON(I)) 10 CONTINUE C C Create the triangulation. C CALL CSSTRI(N,RLAT,RLON, NT,LTRI, IWK,RWK,IER) C C Draw the triangular spherical patches using Ezmap with a C satellite view projection. C C Open GKS. C CALL GOPKS (IERRF, ISZDM) CALL GOPWK (IWKID, LUNIT, IWTYPE) CALL GACWK (IWKID) C C Color table. C CALL GSCR(IWKID,0,1.,1.,1.) CALL GSCR(IWKID,1,0.,0.,0.) CALL GSCR(IWKID,2,1.,0.,0.) CALL GSCR(IWKID,3,0.,0.,1.) C C Plot title. C CALL PLCHHQ(0.50,0.95,':F26:Triangulation',0.030,0.,0.) C C Define the map projection (satellite view). C CALL MAPSTR ('SA - SATELLITE HEIGHT',4.) CALL MAPROJ ('SV - SATELLITE-VIEW',40.,-105.,0.) CALL MAPSET ('MA - MAXIMAL AREA',0.,0.,0.,0.) CALL MAPPOS (0.06, 0.94, 0.02, 0.90) C C Draw the triangles. C CALL GSPLCI(2) DO 130 NP=1,NT DO 140 NS=1,3 C C Calculate points along the arc of each side of the triangle. C CALL MAPGCI(RLAT(LTRI(NS,NP)), RLON(LTRI(NS,NP)), + RLAT(LTRI(MOD(NS,3)+1,NP)), + RLON(LTRI(MOD(NS,3)+1,NP)), + NPARC, ARCLAT, ARCLON) C C Draw the arcs. C CALL MAPIT(RLAT(LTRI(NS,NP)), RLON(LTRI(NS,NP)),0) DO 120 I=1,NPARC CALL MAPIT(ARCLAT(I),ARCLON(I),1) 120 CONTINUE CALL MAPIT(RLAT(LTRI(MOD(NS,3)+1,NP)), + RLON(LTRI(MOD(NS,3)+1,NP)),1) 140 CONTINUE 130 CONTINUE C C Mark the original data points (create a circular C marker by drawing concentric circles at the data points). C CALL GSPLCI(1) DO 420 K=1,8 DO 400 I=1,N CALL NGGCOG(RLAT(I),RLON(I),0.1*REAL(K),ARCLAT,ARCLON,NPARC) CALL MAPIT(ARCLAT(1),ARCLON(1),0) DO 410 J=2,NPARC-1 CALL MAPIT(ARCLAT(J),ARCLON(J),1) 410 CONTINUE CALL MAPIT(ARCLAT(NPARC),ARCLON(NPARC),1) CALL MAPIQ 400 CONTINUE 420 CONTINUE C C Draw the Voronoi polygons. C DO 150 I=1,N C C Get the Voronoi polygon containing the original data point C (RLAT(I),RLON(I)). C IF (I .EQ. 1) THEN CALL CSVORO(N,RLAT,RLON,I,1,IWK,RWK,NTMX,RLATO,RLONO,RC, + NCA,NUMV,NV,IER) ELSE CALL CSVORO(N,RLAT,RLON,I,0,IWK,RWK,NTMX,RLATO,RLONO,RC, + NCA,NUMV,NV,IER) ENDIF C C Draw the polygons. C DO 855 NN=2,NUMV C C Get a polygonal segment. C RLAT1 = RLATO(NV(NN-1)) RLON1 = RLONO(NV(NN-1)) RLAT2 = RLATO(NV(NN)) RLON2 = RLONO(NV(NN)) C C Plot it. C CALL GSPLCI(3) CALL MAPGCI(RLAT1,RLON1,RLAT2,RLON2,NPARC,ARCLAT,ARCLON) CALL MAPIT(RLAT1,RLON1,0) DO 920 J=1,NPARC CALL MAPIT(ARCLAT(J),ARCLON(J),1) 920 CONTINUE CALL MAPIT(RLAT2,RLON2,1) CALL MAPIQ 855 CONTINUE 150 CONTINUE C C Flush mapping buffer and end the picture. C CALL MAPIQ CALL NGEZLOGO CALL FRAME C C Grid the data using NI longitudes and NJ latitudes. C DO 200 I=1,NI PLAT(I) = -90.+(I-1)*2.5 200 CONTINUE DO 210 J=1,NJ PLON(J) = -180.+(J-1)*2.5 210 CONTINUE C C Do the interpolation to the above uniform grid. The call C to the triangulation routine CSSTRI above is not necessary C before calling CSSGRID - it was done in order to illustrate C its usage. C C CALL CSSETR('SIG',2.) C CALL CSSETI('IGR',0) CALL CSSGRID(N,RLAT,RLON,FVAL,NI,NJ,PLAT,PLON,FF,IWK,RWK,IER) C C Draw a contour map of the gridded data on the globe. C CALL MAPSTI ('PE - PERIMETER',0) CALL MAPSTI ('GR - GRID',0) CALL MAPSTC ('OU - OUTLINE DATASET','PS') CALL MAPSTI ('DO - DOTTING OF OUTLINES',1) CALL MAPDRW C CALL GSPLCI(3) CALL GSTXCI(3) CALL GSTXFP(-13,2) CALL CPSETI ('SET - DO SET-CALL FLAG',0) CALL CPSETR ('DPS - LINE LABEL SIZE',0.02) CALL CPSETR ('T2D - TENSION ON THE 2D SPLINES',1.) CALL CPSETI ('CLS - CONTOUR LEVEL SELECTION FLAG',16) CALL CPSETC ('HLT - TURN OFF HIGH/LOW LABELS',' ') CALL CPSETC ('ILT - TURN OFF INFORMATIONAL LABEL',' ') C CALL CPSETR ('XC1 - X COORDINATE AT I = 1',-180.) CALL CPSETR ('XCM - X COORDINATE AT I = M', 180.) CALL CPSETR ('YC1 - Y COORDINATE AT J = 1', -90.) CALL CPSETR ('YCN - Y COORDINATE AT J = N', 90.) CALL CPSETI ('MAP - MAPPING FLAG',1) CALL CPSETR ('ORV - OUT-OF-RANGE VALUE',1.E12) C C Reverse the indices, since CSSGRID returns FF as a function C of latitude and longitude, whereas Conpack wants longitude as C the first dimension. C DO 350 I=1,NI DO 360 J=1,NJ ZDAT(J,I) = FF(I,J) 360 CONTINUE 350 CONTINUE C CALL CPRECT (ZDAT,NJ,NJ,NI,RWRK,IRWRK,IWRK,IIWRK) CALL ARINAM (IAMA,IAMASZ) CALL CPLBAM (ZDAT,RWRK,IWRK,IAMA) CALL CPCLDM (ZDAT,RWRK,IWRK,IAMA,DRAWCL) C C Plot picture title. C CALL SET(0.,1.,0.,1.,0.,1.,0.,1.,1) CALL PLCHHQ(0.50,0.95,':F26:Contour Plot of Gridded Data', + 0.025,0.,0.) C CALL NGEZLOGO CALL FRAME C CALL GDAWK (IWKID) CALL GCLWK (IWKID) CALL GCLKS C STOP END SUBROUTINE CSGENRS(N,RLAT,RLON) C C Generate random positions on a sphere in latitudes and longitudes C (latitude values between -90. and 90. and longitude values C between -180. and 180. C C First select a random longitude, and then select a random value C on the Z axis -R to R (-1 to 1 in the case of the unit sphere). C From the random Z value, calculate the latitude. C DIMENSION RLAT(N),RLON(N) PARAMETER (R=1., R2D=57.29578) C DO 10 I=1,N RLON(I) = -180.+360.*CSRAND() RZ = R*(2.*CSRAND()-1.) RLAT(I) = R2D*ASIN(RZ) 10 CONTINUE C RETURN END REAL FUNCTION CSRAND() C C Portable random number generator. C PARAMETER (MPLIER=16807,MODLUS=2147483647,MOBYMP=127773, + MOMDMP=2836) C INTEGER HVLUE, LVLUE, TESTV, NEXTN SAVE NEXTN DATA JSEED,IFRST/9,0/ C IF (IFRST .EQ. 0) THEN NEXTN = JSEED IFRST = 1 ENDIF C HVLUE = NEXTN / MOBYMP LVLUE = MOD(NEXTN, MOBYMP) TESTV = MPLIER*LVLUE - MOMDMP*HVLUE IF (TESTV .GT. 0) THEN NEXTN = TESTV ELSE NEXTN = TESTV + MODLUS ENDIF CSRAND = REAL(NEXTN)/REAL(MODLUS) C RETURN END SUBROUTINE DRAWCL (XCS,YCS,NCS,IAI,IAG,NAI) DIMENSION XCS(*),YCS(*),IAI(*),IAG(*) C IDR=1 DO 101 I=1,NAI IF (IAI(I) .LT. 0) IDR = 0 101 CONTINUE IF (IDR .NE. 0) CALL CURVED (XCS,YCS,NCS) C RETURN END SUBROUTINE CSGENI (MLOW,MHGH,DLOW,DHGH) C C The routines CSGENI and CSGENPT are used to generate test C data on the globe. The data generated will have approximately C "MLOW" lows and "MHGH" highs within each of the twenty areas C defined by an inscribed icosahedron, a minimum value of C approximately "DLOW" and a maximum value of approximately C "DHGH". CSGENI is called to initialize the process and the C function CSGENPNT returns a computed functional value at the C point (RLAT,RLON). C C The function used is a sum of exponentials. C C This code was originally written by David Kennison at NCAR. C COMMON /CSSBD1/R0,R1,R2 COMMON /CSSCMN/ CMUL(160),NCNT,NLOW,NHGH,RLOW,RHGH,RMIN,RMAX, + XCOC(160),YCOC(160),ZCOC(160) C C Define the vertices and faces of an icosahedron. C DIMENSION XCVI(12),YCVI(12),ZCVI(12),IVOF(3,20) DIMENSION R0(220),R1(220),R2(220) DATA XCVI / .9510565162952 , -.9510565162951 , .4253254041760 , + -.4253254041760 , .4253254041760 , -.4253254041760 , + .4253254041760 , -.4253254041760 , .4253254041760 , + -.4253254041760 , .4253254041760 , -.4253254041760 / DATA YCVI / .0000000000000 , .0000000000000 , .8506508083520 , + -.8506508083520 , .2628655560596 , -.2628655560596 , + -.6881909602356 , .6881909602356 , -.6881909602356 , + .6881909602356 , .2628655560595 , -.2628655560596 / DATA ZCVI / .0000000000000 , .0000000000000 , .0000000000000 , + .0000000000000 , .8090169943749 , -.8090169943749 , + .5000000000000 , -.5000000000000 , -.5000000000000 , + .5000000000000 , -.8090169943749 , .8090169943749 / DATA IVOF / 1, 3, 5 , 1, 5, 7 , 1, 7, 9 , + 1, 9,11 , 1, 3,11 , 3, 8,11 , + 3, 8,10 , 3, 5,10 , 5,10,12 , + 5, 7,12 , 7, 4,12 , 4, 7, 9 , + 4, 6, 9 , 6, 9,11 , 6, 8,11 , + 6, 2, 8 , 2, 4, 6 , 2, 4,12 , + 2,10,12 , 2, 8,10 / DATA DTOR / .017453292519943 / DATA RTOD / 57.2957795130823 / C NLOW = MAX0(0,MIN0(4,MLOW)) NHGH = MAX0(0,MIN0(4,MHGH)) IF (NLOW+NHGH.GE.12 .AND. NLOW.GE.6) NLOW = NLOW-1 IF (NLOW+NHGH.GE.12 .AND. NHGH.GE.6) NHGH = NHGH-1 NCNT = 20*(NLOW+NHGH) C C Compute the center positions for the exponentials which are C summed to create the function. C DO 100 K=1,NCNT INDX = MOD(K-1,20)+1 ILOH = (K+19)/20 WGT1 = R0(K) WGT2 = R1(K) WGT3 = R2(K) WSUM = WGT1+WGT2+WGT3 WGT1 = WGT1/WSUM WGT2 = WGT2/WSUM WGT3 = WGT3/WSUM XTMP = WGT1*XCVI(IVOF(1,INDX))+ + WGT2*XCVI(IVOF(2,INDX))+ + WGT3*XCVI(IVOF(3,INDX)) YTMP = WGT1*YCVI(IVOF(1,INDX))+ + WGT2*YCVI(IVOF(2,INDX))+ + WGT3*YCVI(IVOF(3,INDX)) ZTMP = WGT1*ZCVI(IVOF(1,INDX))+ + WGT2*ZCVI(IVOF(2,INDX))+ + WGT3*ZCVI(IVOF(3,INDX)) TEMP = SQRT(XTMP*XTMP+YTMP*YTMP+ZTMP*ZTMP) XCOC(K) = XTMP/TEMP YCOC(K) = YTMP/TEMP ZCOC(K) = ZTMP/TEMP IF (.NOT.(ILOH .LE. NLOW)) GO TO 110 CMUL(K) = -1. GO TO 120 110 CONTINUE CMUL(K) = +1. 120 CONTINUE 100 CONTINUE C RLOW = DLOW RHGH = DHGH RMIN = +1.E36 RMAX = -1.E36 DO 130 I=-180,175,5 RLON = DTOR*REAL(I) CRLN = COS(RLON) SRLN = SIN(RLON) DO 140 J=-85,85,5 RLAT = DTOR*REAL(J) CRLT = COS(RLAT) SRLT = SIN(RLAT) RVAL = 0.5*(RLOW+RHGH) DO 150 K=1,NCNT DIST = SQRT((CRLT*CRLN-XCOC(K))**2 + +(CRLT*SRLN-YCOC(K))**2 + +(SRLT -ZCOC(K))**2) ANGL = 2.*RTOD*ASIN(DIST/2.) DIST = ANGL/18. RVAL = RVAL+.5*(RHGH-RLOW)*CMUL(K)*EXP(-DIST*DIST) + *(2.-SIN(6.283*DIST)/2.) 150 CONTINUE RMIN = AMIN1(RMIN,RVAL) RMAX = AMAX1(RMAX,RVAL) 140 CONTINUE 130 CONTINUE RETURN END FUNCTION CSGENPNT (RLAT,RLON) COMMON /CSSCMN/ CMUL(160),NCNT,NLOW,NHGH,RLOW,RHGH,RMIN,RMAX, + XCOC(160),YCOC(160),ZCOC(160) DATA DTOR / .017453292519943 / DATA RTOD / 57.2957795130823 / C TLAT = DTOR*RLAT TLON = DTOR*RLON CRLT=COS(TLAT) CRLN=COS(TLON) SRLT=SIN(TLAT) SRLN=SIN(TLON) RVAL=.5*(RLOW+RHGH) DO 100 K=1,NCNT DIST = SQRT((CRLT*CRLN-XCOC(K))**2 + +(CRLT*SRLN-YCOC(K))**2 + +(SRLT -ZCOC(K))**2) ANGL = 2.*RTOD*ASIN(DIST/2.) DIST = ANGL/18. RVAL = RVAL+.5*(RHGH-RLOW)*CMUL(K)*EXP(-DIST*DIST) + *(2.-SIN(6.283*DIST)/2.) 100 CONTINUE CSGENPNT = RLOW+(RHGH-RLOW)*(RVAL-RMIN)/(RMAX-RMIN) C RETURN END BLOCKDATA CSSCCBD C COMMON /CSSBD1/R0,R1,R2 DIMENSION R0(220),R1(220),R2(220) C DATA R0/ + .968,.067,.478,.910,.352,.933,.654,.021,.512,.202,.940,.204, + .379,.793,.288,.267,.357,.128,.703,.737,.803,.915,.511,.762, + .456,.471,.300,.613,.073,.498,.220,.041,.565,.698,.951,.917, + .630,.605,.938,.143,.807,.878,.347,.186,.671,.635,.453,.028, + .763,.157,.765,.566,.072,.276,.328,.528,.747,.627,.141,.821, + .126,.360,.862,.690,.058,.813,.607,.689,.419,.545,.831,.226, + .422,.178,.412,.093,.813,.866,.121,.576,.023,.886,.142,.095, + .162,.470,.623,.910,.097,.764,.730,.223,.124,.593,.913,.183, + .406,.520,.871,.825,.065,.703,.051,.488,.881,.463,.581,.694, + .329,.702,.270,.352,.587,.412,.446,.750,.882,.069,.659,.979, + .833,.390,.202,.957,.982,.115,.140,.389,.635,.011,.213,.701, + .714,.264,.188,.594,.727,.769,.288,.056,.471,.558,.408,.058, + .970,.854,.808,.851,.923,.467,.830,.756,.857,.032,.713,.839, + .147,.852,.228,.783,.863,.441,.483,.577,.705,.671,.171,.432, + .441,.459,.489,.911,.017,.897,.969,.987,.751,.777,.838,.674, + .244,.669,.430,.101,.701,.143,.940,.848,.995,.168,.631,.858, + .608,.114,.435,.313,.785,.606,.746,.226,.065,.234,.137,.082, + .131,.106,.069,.882,.883,.907,.556,.127,.576,.986,.228,.276, + .128,.168,.124,.123/ C DATA R1/ + .023,.003,.880,.088,.237,.017,.170,.368,.123,.239,.250,.006, + .146,.806,.134,.722,.791,.361,.998,.920,.529,.122,.043,.864, + .877,.025,.808,.746,.442,.065,.400,.464,.068,.280,.552,.305, + .297,.722,.673,.420,.961,.923,.426,.107,.729,.560,.829,.520, + .921,.827,.440,.451,.950,.483,.315,.827,.508,.123,.573,.949, + .188,.973,.414,.256,.253,.966,.561,.550,.689,.234,.970,.650, + .157,.396,.757,.885,.956,.587,.405,.877,.414,.845,.328,.363, + .328,.643,.190,.836,.766,.763,.786,.954,.737,.199,.211,.990, + .165,.772,.540,.854,.006,.510,.504,.163,.906,.261,.048,.862, + .848,.454,.740,.262,.299,.068,.625,.627,.711,.815,.464,.477, + .579,.249,.431,.315,.449,.642,.305,.614,.414,.845,.468,.420, + .355,.972,.582,.261,.234,.631,.123,.082,.084,.863,.343,.383, + .930,.968,.011,.641,.784,.474,.118,.362,.723,.549,.678,.172, + .191,.983,.786,.605,.828,.254,.024,.183,.226,.607,.444,.460, + .237,.567,.541,.322,.430,.885,.705,.361,.853,.715,.002,.637, + .190,.119,.999,.913,.668,.677,.085,.859,.660,.871,.464,.488, + .124,.489,.671,.350,.095,.115,.810,.333,.683,.351,.654,.113, + .236,.359,.473,.089,.075,.475,.726,.264,.594,.725,.177,.263, + .402,.262,.122,.062/ C DATA R2/ + .337,.417,.503,.020,.769,.158,.133,.005,.517,.606,.094,.591, + .081,.820,.855,.675,.545,.033,.938,.947,.294,.060,.009,.427, + .646,.559,.684,.721,.781,.291,.892,.118,.708,.395,.138,.476, + .552,.270,.481,.069,.877,.575,.660,.957,.395,.516,.633,.939, + .548,.570,.886,.843,.630,.895,.270,.276,.455,.953,.998,.236, + .244,.889,.354,.952,.284,.492,.428,.837,.762,.909,.906,.639, + .484,.566,.596,.879,.082,.229,.818,.631,.799,.704,.473,.430, + .600,.743,.706,.055,.696,.704,.291,.940,.593,.645,.892,.877, + .137,.321,.714,.899,.230,.620,.538,.714,.186,.134,.593,.268, + .364,.411,.899,.163,.116,.372,.593,.716,.115,.298,.770,.812, + .002,.061,.752,.595,.706,.645,.472,.843,.965,.186,.742,.195, + .806,.280,.910,.992,.414,.503,.260,.778,.915,.159,.941,.030, + .531,.533,.746,.647,.832,.516,.458,.834,.577,.211,.429,.283, + .855,.901,.126,.821,.087,.868,.016,.893,.148,.926,.885,.562, + .429,.145,.340,.343,.304,.281,.374,.835,.814,.120,.482,.646, + .636,.940,.479,.213,.151,.908,.497,.006,.809,.623,.827,.895, + .490,.843,.788,.638,.769,.673,.200,.198,.817,.540,.541,.121, + .821,.915,.956,.635,.035,.438,.280,.671,.377,.760,.884,.528, + .668,.381,.534,.477/ C END