C C icosmap_marsab.f C C Tesselates a sphere and calculates field parameters at the C tesselation points. For input to equivalent source routines. C Parameter (mtot=11562) Dimension qlat(mtot),qlon(mtot) Dimension denlon(mtot),cc(mtot),cs(mtot) Dimension ss(mtot),bfield(mtot) C C Initialize the main field model C C Open(11,file='gsfc1283.coef', C & status='old',access='sequential',form='formatted') C Call Fidd(0.0,0.0,6371.2,1980.,a,b,c,w1,0.0) C Open(31,file='icos2radial.mars', & status='new',access='sequential',form='formatted') C nedge = 35 C C An nedge of C 5 produces 162 dipoles. C 6 produces 252 dipoles C 7 produces 362 dipoles with a separation of 12 degrees C 20 produces a separation of 3.4 degrees and 3612 dipoles C An nedge of C 35 produces a separation of 1.89 degrees and 11562 dipoles C Call icosen(nedge,qlat,qlon,ntot) C C WRITE OUT THE LOCATION OF THE DIPOLES C sthick=40.0 DO 125 I=1,ntot qlat(i) = 90. - qlat(i) If (qlon(i).gt.180.0) then denlon(i) = qlon(i) - 360.0 Else denlon(i) = qlon(i) End if tlat = qlat(i) If (tlat.eq.90.0) tlat=89.9 If (tlat.eq.-90.0) tlat=-89.9 tlon = denlon(i) C Call Fidd(tlat,tlon,6371.2,1980.,x,y,z,f,0.0) C C Induced in main field case C CC ss(i) = -z/f CC cc(i) = -x/f CC cs(i) = y/f C C Radial dipole case C ss(i) = -1.0 cc(i) = 0.0 cs(i) = 0.0 C bfield(i) = f IF (MOD(I,750).EQ.0) WRITE(6,127) I,qLAT(I),DENLON(I), & ss(i),cc(i),cs(i) 127 FORMAT(1X,I10,5F10.2) If (qlat(i).gt.-88.and.qlat(i).lt.88.) then Write(31,55) qlat(i),denlon(i),sthick,ss(i), & cc(i),cs(i) 55 Format(6F10.3) End if 125 CONTINUE C Stop End Subroutine icosen(nedge,qlat,qlon,ntot) C C************************************************************* C THIS PROGRAM GENERATES THE LATITUDE AND LONGITUDE LOCATIONS C ON A SPHERICAL ICOSAHEDRON OF DATA POINTS DISTRIBUTED IN AN C EQUAL SPACING, EQUAL AREA ARRANGEMENT AS IN VESTINE ET AL C GEOMAGNETISM AND GEOELECTRICITY, 15(1963), PP. 73-89 C************************************************************* Parameter (mtot=11562) REAL LLAT,LLONG,NPOLLAT,NPOLLON,SPOLLON,SPOLLAT, * LONSPAC,LATSPAC Real qlat(mtot),qlon(mtot) Real atheta(181),aphi(181,361),qtemp(mtot) Integer nphi(361),ind(mtot) C******************************************************* C THE NUMBER OF POINTS PER EDGE INCLUDES BOTH VERTICES C OF THE TRIANGLE SIDE. ONE IS SUBTRACTED FROM NEDGE TO C GIVE # OF POINTS PER SIDE NOT INCLUDING BOTH VERTICES. C******************************************************* LATSPAC = 180.0/((NEDGE-1)*3) ntot = 1 C******* THESE ARE THE POLE LOCATIONS * NPOLLAT =0.0 NPOLLON = 0.0 SPOLLAT = 180.0 SPOLLON = 0.0 C************************************** qlat(ntot) = npollat qlon(ntot) = npollon ntot = ntot+1 101 FORMAT(F10.3,F10.3) C*************************************************************** C AN ASSUMPTION IS MADE THAT THERE ARE THREE LATITUDE C REGIONS COMPOSING THE ICOSAHEDRON: TWO POLE CAPS CONTAINING C 5 SPHERICAL TRIANGLES, AND THE MIDDLE BAND COMPOSED OF 10 C SPHERICAL TRIANGLES. THE FIRST REGION IS ASSUMED TO GO FROM C 0 TO 60 DEGREES LATITUDE, THE SECOND 60 TO 120, THE THIRD C 120 TO 180 DEGREES LATITUDE. C************************************************************** C*** THE NORTH POLE REGION *********************** DO 10 N = 1, ((NEDGE-1)*3) IF(N.LE.(NEDGE-1)) THEN LONSPAC = 360.0/(REAL(N)*5.0) LLAT = REAL(N)*LATSPAC DO 20 J = 1, (N*5) LLONG = LONSPAC * REAL(J) qlat(ntot) = llat qlon(ntot) = llong ntot = ntot + 1 20 CONTINUE END IF C**** SOUTH POLAR REGION ********************** IF(N.GE.(2*(NEDGE-1))) THEN IF(N.EQ.(3*(NEDGE-1))) GO TO 1000 NN=((NEDGE-1)*3)-N LONSPAC = 360.0/(REAL(NN)*5) LLAT = REAL(N)*LATSPAC DO 30 J = 1, (NN*5) LLONG = LONSPAC * REAL(J) + 36.0 qlat(ntot) = llat qlon(ntot) = llong ntot = ntot + 1 30 CONTINUE END IF C**** THE MIDDLE LATITUDE REGION *********************** IF(N.GT.(NEDGE-1).AND.N.LT.(2*(NEDGE-1))) THEN LONSPAC = 360.0/((NEDGE-1)*5) LLAT = REAL(N)*LATSPAC DO 40 J = 1, (NEDGE-1)*5 LLONG = (LONSPAC*J) + (N-(NEDGE-1))*(0.5*LONSPAC) qlat(ntot) = llat qlon(ntot) = llong ntot = ntot + 1 40 CONTINUE C**************************************************************** C THE NESTED J COUNTER LOOPS IN EACH REGION GENERATE THE C LONGITUDE LOCATIONS FOR A PARTICULAR LATITUDE. THE LATITUDE C IS INCREMENTED BY THE N COUNTER LOOP WHICH INCREMENTS THE C THE LATITUDE SPACING WHICH IS DEPENDENT ON THE CHOICE OF NEDGE. C THE VALUES OF LATITUDE ARE GIVEN IN COLATITUDE AND THE VALUES C OF LONGITUDE ARE GIVEN IN EAST LONGITUDE WITH WRAP AROUND PAST C 360 DEGREES. C**************************************************************** END IF 10 CONTINUE C 1000 qlat(ntot) = spollat qlon(ntot) = spollon C Print *,ntot,' total points on the surface of the sphere' C C Reorder the data points by increasing colatitude C Call reordr(ntot,2,qlat,qlon,dummy,dummy,ind) C i = 1 Do 300 m=1,ntot If (qlon(m).gt.180.0) qlon(m) = qlon(m) - 360.0 If (m.eq.1) then atheta(1) = qlat(m) aphi(1,1) = qlon(m) j = 1 Else if (qlat(m-1).ne.qlat(m)) then nphi(i) = j C C Now reorder within a latitude belt C Call reordr(j,1,qtemp,dummy,dummy,dummy,ind) Do 350 m2=1,j aphi(i,m2) = qtemp(m2) 350 Continue i = i+1 atheta(i) = qlat(m) j = 1 qtemp(j) = qlon(m) If (m.eq.ntot) then nphi(i) = 1 aphi(i,1) = qlon(m) End if Else j = j + 1 qtemp(j) = qlon(m) End if CC Write(11,101) qlat(m),qlon(m) 300 Continue nbands = i C Print *,nbands,' total latitude bands' C Return END SUBROUTINE REORDR (N,IFLAG, A,B,C,D, IND) INTEGER N, IFLAG, IND(N) REAL A(N), B(N), C(N), D(N) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS SUBROUTINE USES A QUICK SORT TO REORDER THE REAL C ARRAY A INTO INCREASING ORDER. A RECORD OF THE PERMUTA- C TIONS APPLIED TO A IS STORED IN IND, AND THESE PERMUTA- C TIONS MAY BE APPLIED TO ONE, TWO, OR THREE ADDITIONAL C VECTORS BY THIS ROUTINE. ANY OTHER VECTOR V MAY BE PER- C MUTED IN THE SAME FASHION BY CALLING SUBROUTINE PERMUT C WITH N, IND, AND V AS PARAMETERS. C A SET OF NODES (X(I),Y(I),Z(I)) AND DATA VALUES W(I) C MAY BE PREPROCESSED BY REORDR FOR INCREASED EFFICIENCY IN C THE TRIANGULATION ROUTINE TRMESH. THIS PREPROCESSING IS C ALSO USEFUL FOR DETECTING DUPLICATE NODES. NOTE THAT THE C FIRST THREE NODES MUST NOT BE COLLINEAR ON INPUT TO C TRMESH. EITHER X, Y, OR Z MAY BE USED AS THE SORT KEY C (ASSOCIATED WITH A). IF THE NODAL COORDINATES ARE IN C TERMS OF LATITUDE AND LONGITUDE, IT IS USUALLY ADVAN- C TAGEOUS TO SORT ON LONGITUDE WITH THE SAME INTERCHANGES C APPLIED TO LATITUDE AND THE DATA VALUES. D WOULD BE A C DUMMY PARAMETER IN THIS CASE. C C INPUT PARAMETERS - N - NUMBER OF NODES. C C IFLAG - NUMBER OF VECTORS TO BE PER- C MUTED. C IFLAG .LE. 0 IF A, B, C, AND C D ARE TO REMAIN UNALTERED. C IFLAG = 1 IF ONLY A IS TO BE C PERMUTED. C IFLAG = 2 IF A AND B ARE TO BE C PERMUTED. C IFLAG = 3 IF A, B, AND C ARE C TO BE PERMUTED. C IFLAG .GE. 4 IF A, B, C, AND D C ARE TO BE PERMUTED. C C A,B,C,D - VECTORS OF LENGTH N TO BE C SORTED (ON THE COMPONENTS OF A) C OR DUMMY PARAMETERS, DEPENDING C ON IFLAG. C C IND - VECTOR OF LENGTH .GE. N. C C N, IFLAG, AND ANY DUMMY PARAMETERS ARE NOT ALTERED BY THIS C ROUTINE. C C OUTPUT PARAMETERS - A,B,C,D - SORTED OR UNALTERED VECTORS, C DEPENDING ON IFLAG. C C IND - SEQUENCE OF INDICES 1,...,N C PERMUTED IN THE SAME FASHION C AS THE REAL VECTORS. THUS C THE ORDERING MAY BE APPLIED C TO A VECTOR V AND STORED IN C W BY SETTING W(I) = C V(IND(I)), OR V MAY BE OVER- C WRITTEN WITH THE ORDERING BY C A CALL TO PERMUT. C C MODULES REFERENCED BY REORDR - QSORT, PERMUT C C*********************************************************** C INTEGER NN, NV C C LOCAL PARAMETERS - C C NN = LOCAL COPY OF N C NV = LOCAL COPY OF IFLAG C NN = N NV = IFLAG CALL QSORT (NN,A, IND) IF (NV .LE. 0) RETURN CALL PERMUT (NN,IND, A ) IF (NV .EQ. 1) RETURN CALL PERMUT (NN,IND, B ) IF (NV .EQ. 2) RETURN CALL PERMUT (NN,IND, C ) IF (NV .EQ. 3) RETURN CALL PERMUT (NN,IND, D ) RETURN END SUBROUTINE PERMUT (N,IP, A ) INTEGER N, IP(N) REAL A(N) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS ROUTINE APPLIES A SET OF PERMUTATIONS TO A VECTOR. C C INPUT PARAMETERS - N - LENGTH OF A AND IP. C C IP - VECTOR CONTAINING THE SEQUENCE OF C INTEGERS 1,...,N PERMUTED IN THE C SAME FASHION THAT A IS TO BE PER- C MUTED. C C A - VECTOR TO BE PERMUTED. C C N AND IP ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETER - A - REORDERED VECTOR REFLECTING THE C PERMUTATIONS DEFINED BY IP. C C MODULES REFERENCED BY PERMUT - NONE C C*********************************************************** C INTEGER NN, K, J, IPJ REAL TEMP C C LOCAL PARAMETERS - C C NN = LOCAL COPY OF N C K = INDEX FOR IP AND FOR THE FIRST ELEMENT OF A IN A C PERMUTATION C J = INDEX FOR IP AND A, J .GE. K C IPJ = IP(J) C TEMP = TEMPORARY STORAGE FOR A(K) C NN = N IF (NN .LT. 2) RETURN K = 1 C C LOOP ON PERMUTATIONS C 1 J = K TEMP = A(K) C C APPLY PERMUTATION TO A. IP(J) IS MARKED (MADE NEGATIVE) C AS BEING INCLUDED IN THE PERMUTATION. C 2 IPJ = IP(J) IP(J) = -IPJ IF (IPJ .EQ. K) GO TO 3 A(J) = A(IPJ) J = IPJ GO TO 2 3 A(J) = TEMP C C SEARCH FOR AN UNMARKED ELEMENT OF IP C 4 K = K + 1 IF (K .GT. NN) GO TO 5 IF (IP(K) .GT. 0) GO TO 1 GO TO 4 C C ALL PERMUTATIONS HAVE BEEN APPLIED. UNMARK IP. C 5 DO 6 K = 1,NN 6 IP(K) = -IP(K) RETURN END SUBROUTINE QSORT (N,X, IND) INTEGER N, IND(N) REAL X(N) C C*********************************************************** C C ROBERT RENKA C OAK RIDGE NATL. LAB. C (615) 576-5139 C C THIS SUBROUTINE USES AN ORDER N*LOG(N) QUICK SORT TO C SORT THE REAL ARRAY X INTO INCREASING ORDER. THE ALGOR- C ITHM IS AS FOLLOWS. IND IS INITIALIZED TO THE ORDERED C SEQUENCE OF INDICES 1,...,N, AND ALL INTERCHANGES ARE C APPLIED TO IND. X IS DIVIDED INTO TWO PORTIONS BY PICKING C A CENTRAL ELEMENT T. THE FIRST AND LAST ELEMENTS ARE COM- C PARED WITH T, AND INTERCHANGES ARE APPLIED AS NECESSARY SO C THAT THE THREE VALUES ARE IN ASCENDING ORDER. INTER- C CHANGES ARE THEN APPLIED SO THAT ALL ELEMENTS GREATER THAN C T ARE IN THE UPPER PORTION OF THE ARRAY AND ALL ELEMENTS C LESS THAN T ARE IN THE LOWER PORTION. THE UPPER AND LOWER C INDICES OF ONE OF THE PORTIONS ARE SAVED IN LOCAL ARRAYS, C AND THE PROCESS IS REPEATED ITERATIVELY ON THE OTHER C PORTION. WHEN A PORTION IS COMPLETELY SORTED, THE PROCESS C BEGINS AGAIN BY RETRIEVING THE INDICES BOUNDING ANOTHER C UNSORTED PORTION. C C INPUT PARAMETERS - N - LENGTH OF THE ARRAY X. C C X - VECTOR OF LENGTH N TO BE SORTED. C C IND - VECTOR OF LENGTH .GE. N. C C N AND X ARE NOT ALTERED BY THIS ROUTINE. C C OUTPUT PARAMETER - IND - SEQUENCE OF INDICES 1,...,N C PERMUTED IN THE SAME FASHION AS X C WOULD BE. THUS, THE ORDERING ON C X IS DEFINED BY Y(I) = X(IND(I)). C C MODULES REFERENCED BY QSORT - NONE C C INTRINSIC FUNCTIONS CALLED BY QSORT - IFIX, FLOAT C C*********************************************************** C C NOTE -- IU AND IL MUST BE DIMENSIONED .GE. LOG(N) WHERE C LOG HAS BASE 2. C C*********************************************************** C INTEGER IU(21), IL(21) INTEGER M, I, J, K, L, IJ, IT, ITT, INDX REAL R, T C C LOCAL PARAMETERS - C C IU,IL = TEMPORARY STORAGE FOR THE UPPER AND LOWER C INDICES OF PORTIONS OF THE ARRAY X C M = INDEX FOR IU AND IL C I,J = LOWER AND UPPER INDICES OF A PORTION OF X C K,L = INDICES IN THE RANGE I,...,J C IJ = RANDOMLY CHOSEN INDEX BETWEEN I AND J C IT,ITT = TEMPORARY STORAGE FOR INTERCHANGES IN IND C INDX = TEMPORARY INDEX FOR X C R = PSEUDO RANDOM NUMBER FOR GENERATING IJ C T = CENTRAL ELEMENT OF X C IF (N .LE. 0) RETURN C C INITIALIZE IND, M, I, J, AND R C DO 1 I = 1,N 1 IND(I) = I M = 1 I = 1 J = N R = .375 C C TOP OF LOOP C 2 IF (I .GE. J) GO TO 10 IF (R .GT. .5898437) GO TO 3 R = R + .0390625 GO TO 4 3 R = R - .21875 C C INITIALIZE K C 4 K = I C C SELECT A CENTRAL ELEMENT OF X AND SAVE IT IN T C IJ = I + IFIX(R*FLOAT(J-I)) IT = IND(IJ) T = X(IT) C C IF THE FIRST ELEMENT OF THE ARRAY IS GREATER THAN T, C INTERCHANGE IT WITH T C INDX = IND(I) IF (X(INDX) .LE. T) GO TO 5 IND(IJ) = INDX IND(I) = IT IT = INDX T = X(IT) C C INITIALIZE L C 5 L = J C C IF THE LAST ELEMENT OF THE ARRAY IS LESS THAN T, C INTERCHANGE IT WITH T C INDX = IND(J) IF (X(INDX) .GE. T) GO TO 7 IND(IJ) = INDX IND(J) = IT IT = INDX T = X(IT) C C IF THE FIRST ELEMENT OF THE ARRAY IS GREATER THAN T, C INTERCHANGE IT WITH T C INDX = IND(I) IF (X(INDX) .LE. T) GO TO 7 IND(IJ) = INDX IND(I) = IT IT = INDX T = X(IT) GO TO 7 C C INTERCHANGE ELEMENTS K AND L C 6 ITT = IND(L) IND(L) = IND(K) IND(K) = ITT C C FIND AN ELEMENT IN THE UPPER PART OF THE ARRAY WHICH IS C NOT LARGER THAN T C 7 L = L - 1 INDX = IND(L) IF (X(INDX) .GT. T) GO TO 7 C C FIND AN ELEMENT IN THE LOWER PART OF THE ARRAY WHCIH IS C NOT SMALLER THAN T C 8 K = K + 1 INDX = IND(K) IF (X(INDX) .LT. T) GO TO 8 C C IF K .LE. L, INTERCHANGE ELEMENTS K AND L C IF (K .LE. L) GO TO 6 C C SAVE THE UPPER AND LOWER SUBSCRIPTS OF THE PORTION OF THE C ARRAY YET TO BE SORTED C IF (L-I .LE. J-K) GO TO 9 IL(M) = I IU(M) = L I = K M = M + 1 GO TO 11 C 9 IL(M) = K IU(M) = J J = L M = M + 1 GO TO 11 C C BEGIN AGAIN ON ANOTHER UNSORTED PORTION OF THE ARRAY C 10 M = M - 1 IF (M .EQ. 0) RETURN I = IL(M) J = IU(M) C 11 IF (J-I .GE. 11) GO TO 4 IF (I .EQ. 1) GO TO 2 I = I - 1 C C SORT ELEMENTS I+1,...,J. NOTE THAT 1 .LE. I .LT. J AND C J-I .LT. 11. C 12 I = I + 1 IF (I .EQ. J) GO TO 10 INDX = IND(I+1) T = X(INDX) IT = INDX INDX = IND(I) IF (X(INDX) .LE. T) GO TO 12 K = I C 13 IND(K+1) = IND(K) K = K - 1 INDX = IND(K) IF (T .LT. X(INDX)) GO TO 13 IND(K+1) = IT GO TO 12 END