c c From kennison@redcloud.scd.ucar.EDU Thu May 14 11:35:02 1992 c Return-Path: c Received: from redcloud.scd.ucar.EDU by sting.hao.ucar.edu.hao (4.1/SMI-4.1) c id AA02939; Thu, 14 May 92 11:34:57 MDT c Date: Thu, 14 May 92 11:34:52 MDT c From: kennison@redcloud.scd.ucar.EDU (Dave Kennison) c Message-Id: <9205141734.AA04298@redcloud.scd.ucar.EDU> c Received: by redcloud.scd.ucar.EDU (5.65/ NCAR Mail Server 04/10/90) c id AA04298; Thu, 14 May 92 11:34:52 MDT c To: foster@sting.hao.ucar.edu c Subject: CPCICA c Status: R c c c I have printed a copy of your message and will look at it this afternoon. c Meanwhile, I am enclosing a test program for CPCICA. It is in IFTRAN. c If you put it in a file called "TestCpcica.ift", you can create the c FORTRAN source, compile, run, and look at the metafile using the commands c c iftran < TestCpcica.ift > TestCpcica.f c ncargf77 -o TestCpcica TestCpcica.f c TestCpcica c ctrans gmeta c c or something similar. It asks you for a projection; type CE, MO, OR, etc. c It then asks you for the pole latitude and longitude; give it any values c that make sense. It then asks you how many spin-up calls to a random c number generator are to be done; give it an integer between 1 and 1000. c At that point, it generates a fractal field of global data and uses c CPCICA to construct a plot of the data. c c I don't know whether this is the best example for CPCICA; still, it may c help you. A word of warning: As part of the testing process, I set this c program up to put a particular color anywhere that isn't associated with c some latitude and longitude. If you use an orthographic projection, for c example, the area around the globe is likely to be red. This looks like c a mistake, but it isn't; it's just the way I left the program set up. c c Program follows ............................................................. PROGRAM GLCOPL C C Parameterize the number of latitudes, the number of longitudes, the C sizes of the real and integer workspaces, and the dimensions of the C cell array. C PARAMETER (NLON=361,NLAT=181,LRWK=5000,LIWK=5000, + ICAM=512,ICAN=512) C C Declare required data arrays and workspace arrays. C DIMENSION ZDAT(NLON,NLAT),RWRK(LRWK),IWRK(LIWK),ICRA(ICAM,ICAN) DIMENSION IASF(13) C C Define a character variable for the projection type. C CHARACTER*2 PTYPE C C Define an array for GKS aspect source flags. C DATA IASF / 13*1 / C C Read the desired projection type. C PRINT * , 'ENTER DESIRED PROJECTION TYPE' READ '(A2)' , PTYPE PRINT * , 'PROJECTION TYPE = ',PTYPE C C Read the desired latitude and longitude. C PRINT * , 'ENTER DESIRED LATITUDE AND LONGITUDE' READ * , PLAT,PLON PRINT * , 'LATITUDE = ',PLAT,' LONGITUDE = ',PLON C C Read the desired number of spin-up calls to FRAN. C PRINT * , 'ENTER DESIRED NUMBER OF SPIN-UP CALLS TO FRAN' READ * , NCTFR PRINT * , 'NUMBER OF CALLS WILL BE ',NCTFR C C Open GKS. C CALL OPNGKS C C Turn off the clipping indicator. C CALL GSCLIP (0) C C Set all aspect source flags to "individual". C CALL GSASF (IASF) C C Define color indices. C CALL DFCLRS C C Generate an array of test data. C CALL GGDINI (0.,1.,NCTFR,.9) C ZMIN= 1.E36 ZMAX=-1.E36 C DO (I=1,NLON) RLON=.017453292519943*(-180.+360.*REAL(I-1)/REAL(NLON-1)) DO (J=1,NLAT) RLAT=.017453292519943*(-90.+180.*REAL(J-1)/REAL(NLAT-1)) ZDAT(I,J)=GGDPNT(RLAT,RLON) + +.5*COS(4.*RLAT) C ZDAT(I,J)=.25*(1.+COS(RLAT))+ C + .25*(1.+SIN(RLON))*COS(RLAT) ZMIN=MIN(ZMIN,ZDAT(I,J)) ZMAX=MAX(ZMAX,ZDAT(I,J)) END DO END DO C DO (I=1,NLON) DO (J=1,NLAT) ZDAT(I,J)=((ZDAT(I,J)-ZMIN)/(ZMAX-ZMIN))*130.-10. END DO END DO C C Put in some special values. C DO (I=175,185) DO (J=85,95) ZDAT(I,J)=1.E36 END DO END DO C C Initialize EZMAP/EZMAPA. C CALL MAPPOS (.01,.99,.01,.99) CALL MAPROJ (PTYPE,PLAT,PLON,0.) CALL MAPSET ('MA - MAXIMAL AREA',0.,0.,0.,0.) CALL MAPINT C C Tell CONPACK not to do the SET call (since it's already been done) and C to use mapping function 1 (EZMAP background), and specify the range of C X and Y coordinates to send into the mapping function. C CALL CPSETI ('SET - DO-SET-CALL FLAG',0) CALL CPSETI ('MAP - MAPPING FLAG',1) 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.) C C Tell CONPACK exactly what contour levels to use. C CALL CPSETI ('CLS - CONTOUR LEVEL SELECTOR',1) CALL CPSETR ('CMN - CONTOUR LEVEL MINIMUM',0.) CALL CPSETR ('CMX - CONTOUR LEVEL MAXIMUM',110.) CALL CPSETR ('CIS - CONTOUR INTERVAL SPECIFIER',10.) C C Tell CONPACK what to use as the out-of-range flag. C CALL CPSETR ('ORV - OUT-OF-RANGE VALUE',1.E12) C C Tell CONPACK what the special value is and what area identifier to C use for special-value areas. C CALL CPSETR ('SPV - SPECIAL VALUE',1.E36) C C Tell CONPACK to use the area identifier plus 2 for the color index C and tell it what area identifiers to use for special-value areas and C outside the grid. C C CALL CPSETI ('CAF - CELL ARRAY FLAG',2) C CALL CPSETI ('PAI - PARAMETER ARRAY INDEX',-1) C CALL CPSETI ('AIA - AREA IDENTIFIER OUTSIDE THE GRID',-100) C CALL CPSETI ('PAID- PARAMETER ARRAY INDEX',-2) C CALL CPSETI ('AIA - AREA IDENTIFIER - SPECIAL-VALUE AREAS',-100) C C Tell CONPACK to call the routine CPSCAE to set cell array elements C and tell it what area identifiers to use for special-value areas and C outside the grid. C CALL CPSETI ('CAF - CELL ARRAY FLAG',-1) CALL CPSETI ('PAI - PARAMETER ARRAY INDEX',-1) CALL CPSETI ('AIA - AREA IDENTIFIER OUTSIDE THE GRID',14) CALL CPSETI ('PAID- PARAMETER ARRAY INDEX',-2) CALL CPSETI ('AIA - AREA IDENTIFIER - SPECIAL-VALUE AREAS',15) C C Initialize the drawing of the contour plot. C CALL CPRECT (ZDAT,NLON,NLON,NLAT,RWRK,LRWK,IWRK,LIWK) C C Fill the cell array. C DO (I=1,ICAM) DO (J=1,ICAN) ICRA(I,J)=0 END DO END DO C CALL CPCICA (ZDAT,RWRK,IWRK,ICRA,ICAM,ICAM,ICAN,0.,0.,1.,1.) C C Draw the cell array. C CALL GCA (CFUX(0.),CFUY(0.),CFUX(1.),CFUY(1.),ICAM,ICAN,1,1, + ICAM,ICAN,ICRA) C C Draw contour lines over the cell array in white. C CALL SFLUSH CALL GSPLCI (1) CALL CPCLDR (ZDAT,RWRK,IWRK) C C Draw a map over the cell array in black. C CALL SFLUSH CALL GSPLCI (0) CALL MAPDRW C C Advance the frame. C CALL FRAME C C Close GKS. C CALL CLSGKS C C Done. C STOP C END SUBROUTINE DFCLRS C C Define 16 different color indices, for indices 0 through 15. The C color corresponding to index 0 is black, the color corresponding to C index 1 is white, and the color corresponding to index 2 is cyan. C CALL GSCR (1, 0,0.000,0.000,0.000) CALL GSCR (1, 1,1.000,1.000,1.000) CALL GSCR (1, 2,0.500,1.000,1.000) C C Bailey's colors: C CALL GSCR (1, 3,0.700,0.700,0.700) CALL GSCR (1, 4,0.750,0.500,1.000) CALL GSCR (1, 5,0.500,0.000,1.000) CALL GSCR (1, 6,0.000,0.000,1.000) CALL GSCR (1, 7,0.000,0.500,1.000) CALL GSCR (1, 8,0.000,1.000,0.600) CALL GSCR (1, 9,0.000,1.000,0.000) CALL GSCR (1,10,0.700,1.000,0.000) CALL GSCR (1,11,1.000,1.000,0.000) CALL GSCR (1,12,1.000,0.750,0.000) CALL GSCR (1,13,1.000,0.380,0.380) CALL GSCR (1,14,1.000,0.000,0.380) CALL GSCR (1,15,1.000,0.000,0.000) C C Roy's colors, with a multiplicative factor: C C CALL GSCR (1, 3,.85*0.000,.85*0.000,.85*1.000) C CALL GSCR (1, 4,.85*0.000,.85*0.381,.85*1.000) C CALL GSCR (1, 5,.85*0.000,.85*0.762,.85*1.000) C CALL GSCR (1, 6,.85*0.000,.85*1.000,.85*0.476) C CALL GSCR (1, 7,.85*0.000,.85*1.000,.85*0.095) C CALL GSCR (1, 8,.85*0.286,.85*1.000,.85*0.000) C CALL GSCR (1, 9,.85*0.667,.85*1.000,.85*0.000) C CALL GSCR (1,10,.85*1.000,.85*0.952,.85*0.000) C CALL GSCR (1,11,.85*1.000,.85*0.571,.85*0.000) C CALL GSCR (1,12,.85*1.000,.85*0.190,.85*0.000) C CALL GSCR (1,13,.85*1.000,.85*0.000,.85*0.190) C CALL GSCR (1,14,.85*1.000,.85*0.000,.85*0.571) C CALL GSCR (1,15,.85*1.000,.85*0.000,.85*0.952) C C Mixes of nothing but red and blue: C C DO (I=3,15) C CALL GSCR (1,I,MAX(0.,MIN(1.,REAL( I-3)/12.)), C + 0., C + MAX(0.,MIN(1.,REAL( 15-I)/12.))) C END DO C C Exponentials for the RGB curves, each centered at a different point; C didn't work too well (I think the constant 10 needs to be increased): C C DO (I=3,15) C CALL GSCR (1,I,MAX(0.,MIN(1.,EXP(-( I-3)*( I-3)/10.))), C + MAX(0.,MIN(1.,EXP(-( I-9)*( I-9)/10.))), C + MAX(0.,MIN(1.,EXP(-(15-I)*(15-I)/10.)))) C END DO C C Upside-down V's for the RGB curves; gives a nice pastel effect: C C DO (I=3,15) C CALL GSCR (1,I,MAX(0.,MIN(1.,1.-REAL(ABS(I- 3)/16.))), C + MAX(0.,MIN(1.,1.-REAL(ABS(I- 9)/16.))), C + MAX(0.,MIN(1.,1.-REAL(ABS(I-15)/16.)))) C END DO C C Same thing, but with steeper slopes, so the colors aren't quite so C pastel: C C DO (I=3,15) C CALL GSCR (1,I,MAX(0.,MIN(1.,1.-REAL(ABS(I- 3)/10.))), C + MAX(0.,MIN(1.,1.-REAL(ABS(I- 9)/10.))), C + MAX(0.,MIN(1.,1.-REAL(ABS(I-15)/10.)))) C END DO C C Same thing, but with even steeper slopes. C C DO (I=3,15) C CALL GSCR (1,I,MAX(0.,MIN(1.,1.-REAL(ABS(I- 3)/8.))), C + MAX(0.,MIN(1.,1.-REAL(ABS(I- 9)/8.))), C + MAX(0.,MIN(1.,1.-REAL(ABS(I-15)/8.)))) C END DO C C Done. C RETURN C END SUBROUTINE GGDINI (DLOW,DHGH,NRND,FRCT) C C The routines GGDINI and GGDPNT are used to generate test data on the C globe, using a method which is essentially fractal and is based on C successive elaborations of an icosahedron. C C To initialize the generation of a particular data field, execute the C statement C C CALL GGDINI (DLOW,DHGH,NRND,FRCT) C C DLOW is the desired minimum value and DHGH the desired maximum value. C NRND is the number of random numbers (from the function FRAN) to be C discarded before beginning to generate the field. FRCT is a fraction C between -1. and 1. which affects the "roughness" of the generated C data; values near 0. will give smoother values than values further C from 0. C C After a call to GGDINI, the function GGDPNT returns the value of the C generated data field at the point (RLAT,RLON); the angles RLAT and C RLON are given in radians. C C Parameters determine how many elaboration steps will be used (NELA) C and therefore how many faces, edges, and vertices will be used (MFCE, C MEDG, and MVTX, respectively). Choose values from the table which C follows: C C NELA MFCE MEDG MVTX C ---- ---- ---- ---- C 0 20 30 12 C 1 100 150 42 C 2 420 630 162 C 3 1700 2550 642 C 4 6820 10230 2562 C 5 27300 40950 10242 C 6 109220 163830 40962 C PARAMETER (NELA=5,MFCE=27300,MEDG=40950,MVTX=10242,NELP=NELA+1) C C NOTE: A COPY OF THE ABOVE PARAMETER STATEMENT MUST BE PLACED IN THE C FUNCTION GGDPNT, AS WELL. C C Information about the structure of the elaborated icosahedron is C stored in variables in the labelled common block GGDCMN. Each of C NFCE entries in IFCE defines a triangular face; it consists of three C pointers to elements of IEDG plus another element which, if non-zero, C points to the faces into which that face was broken by the next C elaboration step. Each of NEDG entries in IEDG defines an edge; C the first two elements are pointers to elements of VRTX, and the C second two elements, if non-zero, are pointers back into IEDG, to C the two edges into which the edge was broken by the next elaboration C step. Each of NVTX entries in VRTX defines a vertex, including the C position and the data value associated with it. C C The arrays LFCE, LEDG, and LVTX point to the beginning and ending of C the faces, edges, and vertices, respectively, added by the step whose C number is given by the second index. C C RLOW and RHGH are just copies of the user's DLOW and DHGH; RMIN and C RMAX are the actual minimum and maximum value of the data values C stored in the array VRTX. C COMMON /GGDCMN/ NFCE,IFCE(4,MFCE),NEDG,IEDG(4,MEDG), + NVTX,VRTX(4,MVTX), + LFCE(2,NELP),LEDG(2,NELP),LVTX(2,NELP), + RLOW,RHGH,RMIN,RMAX SAVE /GGDCMN/ C C Declare local arrays which are used for initialization. C DIMENSION JFCE(3,20),JEDG(2,30),XCVI(12),YCVI(12),ZCVI(12) C C Define the twenty faces of an icosahedron (pointers to edges). C DATA JFCE(1, 1),JFCE(2, 1),JFCE(3, 1) / 1,11, 2 / DATA JFCE(1, 2),JFCE(2, 2),JFCE(3, 2) / 2,19, 3 / DATA JFCE(1, 3),JFCE(2, 3),JFCE(3, 3) / 3,25, 4 / DATA JFCE(1, 4),JFCE(2, 4),JFCE(3, 4) / 4,29, 5 / DATA JFCE(1, 5),JFCE(2, 5),JFCE(3, 5) / 1,14, 5 / DATA JFCE(1, 6),JFCE(2, 6),JFCE(3, 6) / 12,28,14 / DATA JFCE(1, 7),JFCE(2, 7),JFCE(3, 7) / 12,27,13 / DATA JFCE(1, 8),JFCE(2, 8),JFCE(3, 8) / 11,20,13 / DATA JFCE(1, 9),JFCE(2, 9),JFCE(3, 9) / 20,30,21 / DATA JFCE(1,10),JFCE(2,10),JFCE(3,10) / 19,26,21 / DATA JFCE(1,11),JFCE(2,11),JFCE(3,11) / 16,18,26 / DATA JFCE(1,12),JFCE(2,12),JFCE(3,12) / 16,25,17 / DATA JFCE(1,13),JFCE(2,13),JFCE(3,13) / 15,23,17 / DATA JFCE(1,14),JFCE(2,14),JFCE(3,14) / 23,29,24 / DATA JFCE(1,15),JFCE(2,15),JFCE(3,15) / 22,28,24 / DATA JFCE(1,16),JFCE(2,16),JFCE(3,16) / 7, 8,22 / DATA JFCE(1,17),JFCE(2,17),JFCE(3,17) / 6,15, 7 / DATA JFCE(1,18),JFCE(2,18),JFCE(3,18) / 6,18,10 / DATA JFCE(1,19),JFCE(2,19),JFCE(3,19) / 9,30,10 / DATA JFCE(1,20),JFCE(2,20),JFCE(3,20) / 8,27, 9 / C C Define the thirty edges of the icosahedron (pointers to vertices). C DATA JEDG(1, 1),JEDG(2, 1) / 1, 3 / DATA JEDG(1, 2),JEDG(2, 2) / 1, 5 / DATA JEDG(1, 3),JEDG(2, 3) / 1, 7 / DATA JEDG(1, 4),JEDG(2, 4) / 1, 9 / DATA JEDG(1, 5),JEDG(2, 5) / 1,11 / DATA JEDG(1, 6),JEDG(2, 6) / 2, 4 / DATA JEDG(1, 7),JEDG(2, 7) / 2, 6 / DATA JEDG(1, 8),JEDG(2, 8) / 2, 8 / DATA JEDG(1, 9),JEDG(2, 9) / 2,10 / DATA JEDG(1,10),JEDG(2,10) / 2,12 / DATA JEDG(1,11),JEDG(2,11) / 3, 5 / DATA JEDG(1,12),JEDG(2,12) / 3, 8 / DATA JEDG(1,13),JEDG(2,13) / 3,10 / DATA JEDG(1,14),JEDG(2,14) / 3,11 / DATA JEDG(1,15),JEDG(2,15) / 4, 6 / DATA JEDG(1,16),JEDG(2,16) / 4, 7 / DATA JEDG(1,17),JEDG(2,17) / 4, 9 / DATA JEDG(1,18),JEDG(2,18) / 4,12 / DATA JEDG(1,19),JEDG(2,19) / 5, 7 / DATA JEDG(1,20),JEDG(2,20) / 5,10 / DATA JEDG(1,21),JEDG(2,21) / 5,12 / DATA JEDG(1,22),JEDG(2,22) / 6, 8 / DATA JEDG(1,23),JEDG(2,23) / 6, 9 / DATA JEDG(1,24),JEDG(2,24) / 6,11 / DATA JEDG(1,25),JEDG(2,25) / 7, 9 / DATA JEDG(1,26),JEDG(2,26) / 7,12 / DATA JEDG(1,27),JEDG(2,27) / 8,10 / DATA JEDG(1,28),JEDG(2,28) / 8,11 / DATA JEDG(1,29),JEDG(2,29) / 9,11 / DATA JEDG(1,30),JEDG(2,30) / 10,12 / C C Define the vertices of the icosahedron (note radius less than one). C 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 / C C Define constants for converting from degrees to radians and back. C DATA DTOR / .017453292519943 / DATA RTOD / 57.2957795130823 / C C Spin up the random number generator. C DO (I=1,NRND) TEMP=FRAN() END DO C C Initialize the list of faces. C NFCE=20 C DO (I=1,NFCE) IFCE(1,I)=JFCE(1,I) IFCE(2,I)=JFCE(2,I) IFCE(3,I)=JFCE(3,I) IFCE(4,I)=0 END DO C C Initialize the list of edges. C NEDG=30 C DO (I=1,NEDG) IEDG(1,I)=JEDG(1,I) IEDG(2,I)=JEDG(2,I) IEDG(3,I)=0 IEDG(4,I)=0 END DO C C Initialize the list of vertices. C NVTX=12 C DO (I=1,NVTX) TEMP=SQRT(XCVI(I)**2+YCVI(I)**2+ZCVI(I)**2) VRTX(1,I)=XCVI(I)/TEMP VRTX(2,I)=YCVI(I)/TEMP VRTX(3,I)=ZCVI(I)/TEMP VRTX(4,I)=FRAN() END DO C C Define pointers to the existing faces, edges, and vertices. C LFCE(1,1)=1 LFCE(2,1)=NFCE C LEDG(1,1)=1 LEDG(2,1)=NEDG C LVTX(1,1)=1 LVTX(2,1)=NVTX C C Define the initial multiplier for the random number generator. C RMUL=1. C C Do a number of successive elaboration steps. C DO (IELA=2,NELP) C C Update pointers which will tell us what faces, edges, and vertices C were added during this step. C LFCE(1,IELA)=NFCE+1 LEDG(1,IELA)=NEDG+1 LVTX(1,IELA)=NVTX+1 C C Change the multiplier for the random number generator. C RMUL=FRCT*RMUL C C Add a vertex at the midpoint of each edge added during the last step. C DO (I=LEDG(1,IELA-1),LEDG(2,IELA-1)) NVTX=NVTX+1 IF (NVTX.GT.MVTX) GO TO 901 VRTX(1,NVTX)=.5*(VRTX(1,IEDG(1,I))+VRTX(1,IEDG(2,I))) VRTX(2,NVTX)=.5*(VRTX(2,IEDG(1,I))+VRTX(2,IEDG(2,I))) VRTX(3,NVTX)=.5*(VRTX(3,IEDG(1,I))+VRTX(3,IEDG(2,I))) VRTX(4,NVTX)=.5*(VRTX(4,IEDG(1,I))+VRTX(4,IEDG(2,I))) TEMP=SQRT(VRTX(1,NVTX)**2+VRTX(2,NVTX)**2+VRTX(3,NVTX)**2) VRTX(1,NVTX)=VRTX(1,NVTX)/TEMP VRTX(2,NVTX)=VRTX(2,NVTX)/TEMP VRTX(3,NVTX)=VRTX(3,NVTX)/TEMP NEDG=NEDG+2 IF (NEDG.GT.MEDG) GO TO 902 IEDG(3,I)=NEDG-1 IEDG(4,I)=NEDG IEDG(1,NEDG-1)=IEDG(1,I) IEDG(2,NEDG-1)=NVTX IEDG(3,NEDG-1)=0 IEDG(4,NEDG-1)=0 IEDG(1,NEDG )=NVTX IEDG(2,NEDG )=IEDG(2,I) IEDG(3,NEDG )=0 IEDG(4,NEDG )=0 END DO C C Modify the data values at each of the vertices of the new object. C DO (I=1,NVTX) VRTX(4,I)=VRTX(4,I)+RMUL*FRAN() END DO C C Break each face from the last step into four new ones and add them C back into the list. C DO (I=LFCE(1,IELA-1),LFCE(2,IELA-1)) IFCE(4,I)=NFCE+1 NEDG=NEDG+3 IF (NEDG.GT.MEDG) GO TO 902 IEDG(1,NEDG-2)=IEDG(2,IEDG(3,IFCE(1,I))) IEDG(2,NEDG-2)=IEDG(2,IEDG(3,IFCE(2,I))) IEDG(3,NEDG-2)=0 IEDG(4,NEDG-2)=0 IEDG(1,NEDG-1)=IEDG(2,IEDG(3,IFCE(2,I))) IEDG(2,NEDG-1)=IEDG(2,IEDG(3,IFCE(3,I))) IEDG(3,NEDG-1)=0 IEDG(4,NEDG-1)=0 IEDG(1,NEDG )=IEDG(2,IEDG(3,IFCE(3,I))) IEDG(2,NEDG )=IEDG(2,IEDG(3,IFCE(1,I))) IEDG(3,NEDG )=0 IEDG(4,NEDG )=0 NFCE=NFCE+4 IF (NFCE.GT.MFCE) GO TO 903 IF (IEDG(1,IFCE(1,I)).EQ.IEDG(1,IFCE(2,I))) IF (IEDG(2,IFCE(1,I)).EQ.IEDG(1,IFCE(3,I))) IFCE(1,NFCE-3)=IEDG(3,IFCE(1,I)) IFCE(2,NFCE-3)=IEDG(3,IFCE(2,I)) IFCE(3,NFCE-3)=NEDG-2 IFCE(4,NFCE-3)=0 IFCE(1,NFCE-2)=IEDG(4,IFCE(2,I)) IFCE(2,NFCE-2)=IEDG(4,IFCE(3,I)) IFCE(3,NFCE-2)=NEDG-1 IFCE(4,NFCE-2)=0 IFCE(1,NFCE-1)=IEDG(3,IFCE(3,I)) IFCE(2,NFCE-1)=IEDG(4,IFCE(1,I)) IFCE(3,NFCE-1)=NEDG IFCE(4,NFCE-1)=0 ELSE IFCE(1,NFCE-3)=IEDG(3,IFCE(1,I)) IFCE(2,NFCE-3)=IEDG(3,IFCE(2,I)) IFCE(3,NFCE-3)=NEDG-2 IFCE(4,NFCE-3)=0 IFCE(1,NFCE-2)=IEDG(3,IFCE(3,I)) IFCE(2,NFCE-2)=IEDG(4,IFCE(2,I)) IFCE(3,NFCE-2)=NEDG-1 IFCE(4,NFCE-2)=0 IFCE(1,NFCE-1)=IEDG(4,IFCE(1,I)) IFCE(2,NFCE-1)=IEDG(4,IFCE(3,I)) IFCE(3,NFCE-1)=NEDG IFCE(4,NFCE-1)=0 END IF ELSE IF (IEDG(1,IFCE(1,I)).EQ.IEDG(2,IFCE(2,I))) IF (IEDG(2,IFCE(1,I)).EQ.IEDG(1,IFCE(3,I))) IFCE(1,NFCE-3)=IEDG(3,IFCE(1,I)) IFCE(2,NFCE-3)=IEDG(4,IFCE(2,I)) IFCE(3,NFCE-3)=NEDG-2 IFCE(4,NFCE-3)=0 IFCE(1,NFCE-2)=IEDG(3,IFCE(2,I)) IFCE(2,NFCE-2)=IEDG(4,IFCE(3,I)) IFCE(3,NFCE-2)=NEDG-1 IFCE(4,NFCE-2)=0 IFCE(1,NFCE-1)=IEDG(3,IFCE(3,I)) IFCE(2,NFCE-1)=IEDG(4,IFCE(1,I)) IFCE(3,NFCE-1)=NEDG IFCE(4,NFCE-1)=0 ELSE IFCE(1,NFCE-3)=IEDG(3,IFCE(1,I)) IFCE(2,NFCE-3)=IEDG(4,IFCE(2,I)) IFCE(3,NFCE-3)=NEDG-2 IFCE(4,NFCE-3)=0 IFCE(1,NFCE-2)=IEDG(3,IFCE(2,I)) IFCE(2,NFCE-2)=IEDG(3,IFCE(3,I)) IFCE(3,NFCE-2)=NEDG-1 IFCE(4,NFCE-2)=0 IFCE(1,NFCE-1)=IEDG(4,IFCE(1,I)) IFCE(2,NFCE-1)=IEDG(4,IFCE(3,I)) IFCE(3,NFCE-1)=NEDG IFCE(4,NFCE-1)=0 END IF ELSE IF (IEDG(1,IFCE(1,I)).EQ.IEDG(1,IFCE(3,I))) IF (IEDG(2,IFCE(1,I)).EQ.IEDG(1,IFCE(2,I))) IFCE(1,NFCE-3)=IEDG(3,IFCE(2,I)) IFCE(2,NFCE-3)=IEDG(4,IFCE(1,I)) IFCE(3,NFCE-3)=NEDG-2 IFCE(4,NFCE-3)=0 IFCE(1,NFCE-2)=IEDG(4,IFCE(2,I)) IFCE(2,NFCE-2)=IEDG(4,IFCE(3,I)) IFCE(3,NFCE-2)=NEDG-1 IFCE(4,NFCE-2)=0 IFCE(1,NFCE-1)=IEDG(3,IFCE(1,I)) IFCE(2,NFCE-1)=IEDG(3,IFCE(3,I)) IFCE(3,NFCE-1)=NEDG IFCE(4,NFCE-1)=0 ELSE IFCE(1,NFCE-3)=IEDG(4,IFCE(1,I)) IFCE(2,NFCE-3)=IEDG(4,IFCE(2,I)) IFCE(3,NFCE-3)=NEDG-2 IFCE(4,NFCE-3)=0 IFCE(1,NFCE-2)=IEDG(3,IFCE(2,I)) IFCE(2,NFCE-2)=IEDG(4,IFCE(3,I)) IFCE(3,NFCE-2)=NEDG-1 IFCE(4,NFCE-2)=0 IFCE(1,NFCE-1)=IEDG(3,IFCE(1,I)) IFCE(2,NFCE-1)=IEDG(3,IFCE(3,I)) IFCE(3,NFCE-1)=NEDG IFCE(4,NFCE-1)=0 END IF ELSE IF (IEDG(2,IFCE(1,I)).EQ.IEDG(1,IFCE(2,I))) IFCE(1,NFCE-3)=IEDG(3,IFCE(2,I)) IFCE(2,NFCE-3)=IEDG(4,IFCE(1,I)) IFCE(3,NFCE-3)=NEDG-2 IFCE(4,NFCE-3)=0 IFCE(1,NFCE-2)=IEDG(3,IFCE(3,I)) IFCE(2,NFCE-2)=IEDG(4,IFCE(2,I)) IFCE(3,NFCE-2)=NEDG-1 IFCE(4,NFCE-2)=0 IFCE(1,NFCE-1)=IEDG(3,IFCE(1,I)) IFCE(2,NFCE-1)=IEDG(4,IFCE(3,I)) IFCE(3,NFCE-1)=NEDG IFCE(4,NFCE-1)=0 ELSE IFCE(1,NFCE-3)=IEDG(4,IFCE(1,I)) IFCE(2,NFCE-3)=IEDG(4,IFCE(2,I)) IFCE(3,NFCE-3)=NEDG-2 IFCE(4,NFCE-3)=0 IFCE(1,NFCE-2)=IEDG(3,IFCE(2,I)) IFCE(2,NFCE-2)=IEDG(3,IFCE(3,I)) IFCE(3,NFCE-2)=NEDG-1 IFCE(4,NFCE-2)=0 IFCE(1,NFCE-1)=IEDG(3,IFCE(1,I)) IFCE(2,NFCE-1)=IEDG(4,IFCE(3,I)) IFCE(3,NFCE-1)=NEDG IFCE(4,NFCE-1)=0 END IF END IF IFCE(1,NFCE )=NEDG-2 IFCE(2,NFCE )=NEDG-1 IFCE(3,NFCE )=NEDG IFCE(4,NFCE )=0 END DO C C Update pointers which will tell us what faces, edges, and vertices C were added during this step. C LFCE(2,IELA)=NFCE LEDG(2,IELA)=NEDG LVTX(2,IELA)=NVTX C C End of elaboration step. C END DO C C Copy the desired lowest and highest values to COMMON. C RLOW=DLOW RHGH=DHGH C C Compute the minimum and maximum value in the data. C RMIN=+1.E36 RMAX=-1.E36 C DO (I=1,NVTX) RMIN=MIN(RMIN,VRTX(4,I)) RMAX=MAX(RMAX,VRTX(4,I)) END DO C C Done. C RETURN C C Error exits. C 901 PRINT * , 'STOP IN GGDINI - TOO MANY VERTICES' STOP C 902 PRINT * , 'STOP IN GGDINI - TOO MANY EDGES' STOP C 903 PRINT * , 'STOP IN GGDINI - TOO MANY FACES' STOP C END FUNCTION GGDPNT (RLAT,RLON) C C Following is a copy of the PARAMETER statement from GGDINI. C PARAMETER (NELA=5,MFCE=27300,MEDG=40950,MVTX=10242,NELP=NELA+1) C C Define the common block for the package. C COMMON /GGDCMN/ NFCE,IFCE(4,MFCE),NEDG,IEDG(4,MEDG), + NVTX,VRTX(4,MVTX), + LFCE(2,NELP),LEDG(2,NELP),LVTX(2,NELP), + RLOW,RHGH,RMIN,RMAX SAVE /GGDCMN/ C C Define constants for converting from degrees to radians and back. C DATA DTOR / .017453292519943 / DATA RTOD / 57.2957795130823 / C C Compute the position of the point on the surface of the unit sphere. C X0=COS(RLAT)*COS(RLON) Y0=COS(RLAT)*SIN(RLON) Z0=SIN(RLAT) C C Find out what face that point falls above. C IFST=1 NUMB=20 ITRY=1 C REPEAT IPOS=IFST LOOP X1=VRTX(1,IEDG(1,IFCE(1,IPOS))) Y1=VRTX(2,IEDG(1,IFCE(1,IPOS))) Z1=VRTX(3,IEDG(1,IFCE(1,IPOS))) W1=VRTX(4,IEDG(1,IFCE(1,IPOS))) X2=VRTX(1,IEDG(2,IFCE(1,IPOS))) Y2=VRTX(2,IEDG(2,IFCE(1,IPOS))) Z2=VRTX(3,IEDG(2,IFCE(1,IPOS))) W2=VRTX(4,IEDG(2,IFCE(1,IPOS))) IF (IEDG(1,IFCE(2,IPOS)).NE.IEDG(1,IFCE(1,IPOS)).AND. + IEDG(1,IFCE(2,IPOS)).NE.IEDG(2,IFCE(1,IPOS))) X3=VRTX(1,IEDG(1,IFCE(2,IPOS))) Y3=VRTX(2,IEDG(1,IFCE(2,IPOS))) Z3=VRTX(3,IEDG(1,IFCE(2,IPOS))) W3=VRTX(4,IEDG(1,IFCE(2,IPOS))) ELSE X3=VRTX(1,IEDG(2,IFCE(2,IPOS))) Y3=VRTX(2,IEDG(2,IFCE(2,IPOS))) Z3=VRTX(3,IEDG(2,IFCE(2,IPOS))) W3=VRTX(4,IEDG(2,IFCE(2,IPOS))) END IF A1=Y3*Z2-Y2*Z3 B1=X2*Z3-X3*Z2 C1=X3*Y2-X2*Y3 A2=Y1*Z3-Y3*Z1 B2=X3*Z1-X1*Z3 C2=X1*Y3-X3*Y1 A3=Y2*Z1-Y1*Z2 B3=X1*Z2-X2*Z1 C3=X2*Y1-X1*Y2 IF (ITRY.EQ.1) XM=(X1+X2+X3)/3. YM=(Y1+Y2+Y3)/3. ZM=(Z1+Z2+Z3)/3. S1=SIGN(1.,A1*XM+B1*YM+C1*ZM) S2=SIGN(1.,A2*XM+B2*YM+C2*ZM) S3=SIGN(1.,A3*XM+B3*YM+C3*ZM) T1=SIGN(1.,A1*X0+B1*Y0+C1*Z0) T2=SIGN(1.,A2*X0+B2*Y0+C2*Z0) T3=SIGN(1.,A3*X0+B3*Y0+C3*Z0) EXIT IF (S1.EQ.T1.AND.S2.EQ.T2.AND.S3.EQ.T3) ELSE D1=ABS(A1*X0+B1*Y0+C1*Z0)/SQRT(A1*A1+B1*B1+C1*C1) D2=ABS(A2*X0+B2*Y0+C2*Z0)/SQRT(A2*A2+B2*B2+C2*C2) D3=ABS(A3*X0+B3*Y0+C3*Z0)/SQRT(A3*A3+B3*B3+C3*C3) DIST=MIN(D1,D2,D3) IF (IBST.EQ.0.OR.DIST.LT.DBST) IBST=IPOS DBST=DIST END IF END IF IF (IPOS.EQ.IFST+NUMB-1) IF (ITRY.EQ.1) IPOS=IFST-1 ITRY=2 IBST=0 DBST=0. ELSE IPOS=IBST X1=VRTX(1,IEDG(1,IFCE(1,IPOS))) Y1=VRTX(2,IEDG(1,IFCE(1,IPOS))) Z1=VRTX(3,IEDG(1,IFCE(1,IPOS))) W1=VRTX(4,IEDG(1,IFCE(1,IPOS))) X2=VRTX(1,IEDG(2,IFCE(1,IPOS))) Y2=VRTX(2,IEDG(2,IFCE(1,IPOS))) Z2=VRTX(3,IEDG(2,IFCE(1,IPOS))) W2=VRTX(4,IEDG(2,IFCE(1,IPOS))) IF (IEDG(1,IFCE(2,IPOS)).NE.IEDG(1,IFCE(1,IPOS)).AND. + IEDG(1,IFCE(2,IPOS)).NE.IEDG(2,IFCE(1,IPOS))) X3=VRTX(1,IEDG(1,IFCE(2,IPOS))) Y3=VRTX(2,IEDG(1,IFCE(2,IPOS))) Z3=VRTX(3,IEDG(1,IFCE(2,IPOS))) W3=VRTX(4,IEDG(1,IFCE(2,IPOS))) ELSE X3=VRTX(1,IEDG(2,IFCE(2,IPOS))) Y3=VRTX(2,IEDG(2,IFCE(2,IPOS))) Z3=VRTX(3,IEDG(2,IFCE(2,IPOS))) W3=VRTX(4,IEDG(2,IFCE(2,IPOS))) END IF EXIT END IF END IF IPOS=IPOS+1 END LOOP IFST=IFCE(4,IPOS) NUMB=4 ITRY=1 UNTIL (IFST.EQ.0) C C Compute a data value at the point. C IF ((X1-X0)**2+(Y1-Y0)**2+(Z1-Z0)**2.LT..000001) RVAL=W1 ELSE AP=Y1*(Z3-Z2)+Y2*(Z1-Z3)+Y3*(Z2-Z1) BP=X1*(Z2-Z3)+X2*(Z3-Z1)+X3*(Z1-Z2) CP=X1*(Y3-Y2)+X2*(Y1-Y3)+X3*(Y2-Y1) DP=X1*Y2*Z3-X1*Y3*Z2-X2*Y1*Z3+X3*Y1*Z2+X2*Y3*Z1-X3*Y2*Z1 TP=-DP/(AP*X0+BP*Y0+CP*Z0) X4=TP*X0 Y4=TP*Y0 Z4=TP*Z0 A1=Y3*Z2-Y2*Z3 B1=X2*Z3-X3*Z2 C1=X3*Y2-X2*Y3 TP=-(A1*X1+B1*Y1+C1*Z1)/(A1*(X4-X1)+B1*(Y4-Y1)+C1*(Z4-Z1)) X5=X1+TP*(X4-X1) Y5=Y1+TP*(Y4-Y1) Z5=Z1+TP*(Z4-Z1) W5=W2+(W3-W2)*SQRT(((X5-X2)**2+(Y5-Y2)**2+(Z5-Z2)**2)/ + ((X3-X2)**2+(Y3-Y2)**2+(Z3-Z2)**2)) RVAL=W1+(W5-W1)*SQRT(((X4-X1)**2+(Y4-Y1)**2+(Z4-Z1)**2)/ + ((X5-X1)**2+(Y5-Y1)**2+(Z5-Z1)**2)) END IF C C Adjust the value to the desired range. C GGDPNT=RLOW+(RHGH-RLOW)*(RVAL-RMIN)/(RMAX-RMIN) C C Done. C RETURN C END FUNCTION FRAN() C C Pseudo-random-number generator. C DOUBLE PRECISION X SAVE X DATA X / 2.718281828459045 / X=MOD(9821.D0*X+.211327D0,1.D0) FRAN=REAL(X) RETURN END