* * marsabs_edcgnnems.f in /home/purucker/mgs_science_mag ********************************************************************* * * 2. NAME DERIVATION - Equivalent Source MAPping - Whole world * 1.89 degree dipole spacing - Conjugate Gradient * technique. Icosahedral tesselation of sphere. * * 3. FUNCTION - Solves the design equation using the Conjugate * gradient technique. Sparse storage mode obs (m) > dipoles (n) * 11550 dipoles Row Preconditioning and Weights * * 4. LANGUAGE - Standard Fortran 77. * * 5. IMPORTANT VARIABLES - * NX: THE MAXIMUM NUMBER OF DIPOLES ALLOWED. * MX: The maximum number of observations allowed * These need to be modified at the front of routines * esmap2 and linbcgd in parameter statements. * * 7. FILE REFERENCES - * Unit 11: Input. Main field model. * Unit 18: Output. Decomposed matrices. * Unit 20: Input. Dipole locations (icosahedron subdivision) * Unit 32: Input. Observations (B,X,Y,Z) of field with locs * Unit 33: Input. Observations (B,X,Y,Z) of field with locs * Unit 40: Output. Dipole parameters (Magnetization or suscep) * Unit 45: Output. Fit to observations and error. * * 8. HARDWARE/SOFTWARE ENVIRONMENT - Works on PC running Linux * No special libraries required, Single precision version * * 9. HISTORY - * * VERSION DATE AUTHOR DESCRIPTION * _______ ____ ______ ___________ * * 1.0 2000 M. Purucker PC version * ******************************************************************* * C Open(11,file='gsfc1283.coef',status='old', C & access='sequential',form='formatted') C Open(18,file='mars2.edcgnnema',status='old', & form='unformatted') C Open(20,file='icos2radial.mars',status='old', & access='sequential',form='formatted') C Open(40, & file='test.xyz', & status='new', & form='formatted') C Open(45, & file='test1.xyz', & status='new', & form='formatted') C C For formatted data C Open(32,file='mgsab.dat', & status='old',access='sequential',form='formatted') Open(33,file='98spo12.bin', & status='old',access='sequential',form='formatted') C C Program switches C CC Print *,'Are you entering Scalar(1) or Vector(2) data?' CC Read(5,669) isvdata 669 Format(I1) isvdata = 2 C CC Uses preconditioning if you say (2) CC No preconditioning if you say (1) CC Read(5,669) incerr incerr=2 C CC Print *,'Is this 7F10.2 data (1) or observations (2)?' CC Read(5,669) idtype idtype = 1 C CC Print *,'Using magnetization (1) or susceptibility (2)?' CC Read(5,669) imagsus imagsus = 1 C it17 = 2 If (isvdata.eq.1) then Call SDATA(it17,incerr,idtype,ci) Else if (isvdata.eq.2) then Call VDATA(it17,idtype) End if C Call esmap2(ci,imagsus,incerr) C End Subroutine esmap2(ci,imagsus,incerr) C PARAMETER (mx=40546,NX=11550,nmax=25000000) C C MX is the number of observations C NX is the number of dipoles C DIMENSION W(NX),E(NX),var(nx) DIMENSION SS(NX),CC(NX),CS(NX),DW(MX),D(NX) DIMENSION CTH(NX),STH(NX),P(NX),RTHICK(NX) DIMENSION BD(NX),XL(3),cdist(NX),sr(nx) DIMENSION DFDP1(NX),DFDP2(NX),DFDP3(NX),DFDP4(NX) DIMENSION THICK(NX),ALAT(NX),ALON(NX) Dimension srow(mx) Real sc(nmax) Integer ijc(nmax) C C SPACG = DIPOLE SPACING IN DEGREES C LONGRNG = LONGITUDE RANGE specified in Unit 10 C 1 = -180 TO 180 C 2 = 0 TO 360 C Thresh = Retaining dipoles with thresh km. of observation C ithcv = 0 longrng = 2 spacg = 1.89 C C PRINT OUT THE Critical Information C WRITE(6,1) SPACG,LONGRNG 1 FORMAT(/,1X,'Dipole spacing = ',F10.2,/,1x, * 'LONGITUDE RANGE = ',I4) C WRITE(6,814) 814 FORMAT(1X,'CONSTANT THICKNESS CRUST OF 40 KM ASSUMED.') C C R = Mars RADIUS (KM) C R = 3393.5 C C SELECT THE LOCATION OF THE NP DIPOLES C Rewind 20 CALL SELECT(nx,ALAT,ALON,NP, * THICK,SS,CC,CS,CTH,STH,BD,LONGRNG) C C INITIALIZE THE GEOMAGNETIC FIELD MODEL C CC CALL FIDD(0.,0.,R,1980.,A,B,C,W1,0.0) C 4 FORMAT(1X,'LATITUDE LONGITUDE ALTITUDE OBSERVATION',3X, * 'STD.ERR') C 62 FORMAT(/,15X,'TOTAL FIELD ANOMALY VALUES') 63 FORMAT(/,15X,'DELTA X ANOMALY VALUES') 64 FORMAT(/,15X,'DELTA Y ANOMALY VALUES') 65 FORMAT(/,15X,'DELTA Z ANOMALY VALUES') 35 FORMAT(1X,6F10.3) C C Read the accumulated matrices from disk C Read(18) mp Read(18) np Read(18) sc Read(18) ijc Read(18) dw C Write(6,889) ijc(ijc(1)-1)-1 889 Format(/,1x,'Size of ijc and sc vectors is ',I20) Print *,'mp = ',mp,' np = ',np C scmin=0.0 scmax=0.0 Do 304 m=1,ijc(ijc(1)-1)-1 If (sc(m).ne.sc(m)) print *,sc(m),m,' error' If (abs(sc(m)).gt.scmax) scmax=sc(m) If (abs(sc(m)).lt.scmin) scmin=sc(m) 304 Continue Print *,'sc min and max = ',scmin,scmax C dwmin=0.0 dwmax=0.0 Do 305 m=1,mx If (abs(dw(m)).gt.dwmax) dwmax=dw(m) If (abs(dw(m)).lt.dwmin) dwmin=dw(m) 305 Continue Print *,'dw min and max = ',dwmin,dwmax itol = 1 tol = 0.000001 itmax = 200 C C Precondition the matrix C unless you are using weights C If (incerr.eq.2) then Do 377 m=1,mp srow(m) = 0.0 Do 378 k=ijc(m),ijc(m+1)-1 C C Sum over the rows C srow(m) = srow(m) + abs(sc(k)) 378 Continue 377 Continue C Do 379 m=1,mp dw(m) = dw(m)/srow(m) Do 380 k=ijc(m),ijc(m+1)-1 sc(k)=sc(k)/srow(m) 380 Continue 379 Continue End if C C Solve using a cg technique designed especially for a design matrix C Call linbcgd(mp,np,dw,p,itol,tol,itmax,iter,err,sc,ijc) C Print *,'Conjugate gradient complete after ',iter,' iterations' Print *,' with an error of ',err C C WRITE OUT THE DIPOLES AND CALCULATE THEIR STANDARD DEVIATION C If (imagsus.eq.2) then WRITE(6,39) 39 FORMAT(//,1X,'LATITUDE LONGITUDE SUSC*THICKNESS') Else if (imagsus.eq.1) then Write(6,882) 882 Format(//,1x,'Latitude Longitude Magnetization') End if C SD1 = 0.0 SD2 = 0.0 c DO 40 I=1,NP SD1 = SD1 + P(I)**2 SD2 = SD2 + P(I) C IF (MOD(I,600).EQ.0) THEN WRITE(6,45) ALAT(I),ALON(I),p(i) END IF WRITE(40,45) ALON(I),ALAT(I),P(i),ss(i),cc(i),cs(i) 45 FORMAT(2F10.2,F15.4,3F10.3) 40 CONTINUE C SD3 = SQRT((NP*SD1-SD2**2)/(NP*(NP-1))) If (imagsus.eq.2) then WRITE(6,178) SD3 178 FORMAT(/,1X,'STD.DEVIATION OF SUSC*thick = ',F10.4,/) Else if (imagsus.eq.1) then Write(6,172) SD3 172 Format(/,1x,'STD.Deviation of magnetization = ',F10.4,/) End if C C NOW COMPARE THE FIELDS CALCULATED FROM THE DIPOLES C WITH THE OBSERVED FIELDS C WRITE(6,234) 234 FORMAT(1X,'LATITUDE LONGITUDE ALTITUDE OBSERVATION',3X, * 'CALCULATED FIELD Std.Dev.on Bin') C NBA=0 NXA=0 NYA=0 NZA=0 ICOUNT = 0 SXYI1 = 0.0 SXI1 = 0.0 SYI1 = 0.0 S2XI1 = 0.0 S2YI1 = 0.0 SXYI2 = 0.0 SXI2 = 0.0 SYI2 = 0.0 S2XI2 = 0.0 S2YI2 = 0.0 SXYI3 = 0.0 SXI3 = 0.0 SYI3 = 0.0 S2XI3 = 0.0 S2YI3 = 0.0 SXI4 = 0.0 SYI4 = 0.0 S2XI4 = 0.0 SXYI4 = 0.0 S2YI4 = 0.0 SF1 = 0.0 S2F1 = 0.0 SF2 = 0.0 S2F2 = 0.0 SF3 = 0.0 S2F3 = 0.0 SF4 = 0.0 S2F4 = 0.0 C 195 CALL RDATA(NX,ALON,NP,SS,CC,CS,CTH,STH,THICK,FR,FT,FP,YC,P, * DFDP1,DFDP2,DFDP3,DFDP4,Y,SPACG,ICOMP,IEND,XL,BD, * STDERR,LONGRNG,cdist,imagsus) C IF (IEND.EQ.1) GO TO 199 C IF (LONGRNG.EQ.2.AND.XL(2).LT.0.0) XL(2)=360.+XL(2) C IF (ICOMP.EQ.1.AND.NBA.EQ.0) WRITE(6,62) IF (ICOMP.EQ.2.AND.NXA.EQ.0) WRITE(6,63) IF (ICOMP.EQ.3.AND.NYA.EQ.0) WRITE(6,64) IF (ICOMP.EQ.4.AND.NZA.EQ.0) WRITE(6,65) ICOUNT = ICOUNT + 1 C C ACCUMULATE THE INFORMATION NEEDED TO CALCULATE THE C CORRELATION COEFFICIENT AND THE SIGMA C IF (ICOMP.EQ.1) THEN CALCF=YC DELTA1 = Y-CALCF NBA=NBA+1 SF1 = SF1 + DELTA1 S2F1 = S2F1 + DELTA1**2 SXYI1=SXYI1+(CALCF*Y) SXI1=SXI1+CALCF SYI1=SYI1+Y S2XI1=S2XI1+(CALCF**2) S2YI1=S2YI1+(Y**2) ELSE IF (ICOMP.EQ.2) THEN CALCF=-FT DELTA2 = Y-CALCF NXA=NXA+1 SF2 = SF2 + DELTA2 S2F2 = S2F2 + DELTA2**2 SXYI2=SXYI2+(CALCF*Y) SXI2=SXI2+CALCF SYI2=SYI2+Y S2XI2=S2XI2+(CALCF**2) S2YI2=S2YI2+(Y**2) ELSE IF (ICOMP.EQ.3) THEN CALCF=FP DELTA3 = Y-CALCF NYA=NYA+1 SF3 = SF3 + DELTA3 S2F3 = S2F3 + DELTA3**2 SXYI3=SXYI3+(CALCF*Y) SXI3=SXI3+CALCF SYI3=SYI3+Y S2XI3=S2XI3+(CALCF**2) S2YI3=S2YI3+(Y**2) ELSE IF (ICOMP.EQ.4) THEN CALCF=-FR DELTA4 = Y-CALCF NZA=NZA+1 SF4 = SF4+DELTA4 S2F4=S2F4+DELTA4**2 SXYI4=SXYI4+(CALCF*Y) SXI4=SXI4+CALCF SYI4=SYI4+Y S2XI4=S2XI4+(CALCF**2) S2YI4=S2YI4+(Y**2) END IF C IF (MOD(ICOUNT,1000).EQ.0.or.icount.eq.1) * WRITE(6,35) XL(1),XL(2),XL(3),Y,CALCF,stderr Write(45,49) XL(1),XL(2),XL(3),Y,Calcf,stderr 49 Format(6F10.3) C GO TO 195 199 CONTINUE C Rewind (55) C IF (NBA.GT.0) THEN CCB=((REAL(NBA)*SXYI1)-(SXI1*SYI1))/ * SQRT((REAL(NBA)*S2XI1-SXI1**2)*(REAL(NBA)*S2YI1-SYI1**2)) SIG1 = SQRT((NBA*S2F1-SF1**2)/(NBA*(NBA-1))) WRITE(6,726) CCB,SIG1 726 FORMAT(1X,'COR. COEF. FOR SCALAR = ',F8.3, * ' SIGMA ON DEL = ',F8.4) END IF IF (NXA.GT.0) THEN CCXin = (REAL(NXA)*S2XI2-SXI2**2)*(REAL(NXA)*S2YI2-SYI2**2) SIG2ins = (NXA*S2F2-SF2**2)/(NXA*(NXA-1)) If (sig2ins.gt.0.0.and.ccxin.gt.0.0) then SIG2 = SQRT(sig2ins) CCX=((REAL(NXA)*SXYI2)-(SXI2*SYI2))/sqrt(ccxin) WRITE(6,727) CCX,SIG2 727 FORMAT(1X,'COR. COEF. FOR X COMP = ',F8.3, * ' SIGMA ON DEL = ',F8.4) Else Write(6,797) ccxin,sig2ins,nxa 797 Format(E15.4,E15.4,I6) End if END IF IF (NYA.GT.0) THEN ccyin = (REAL(NYA)*S2XI3-SXI3**2)*(REAL(NYA)*S2YI3-SYI3**2) SIG3ins = (NYA*S2F3-SF3**2)/(NYA*(NYA-1)) If (sig3ins.gt.0.0.and.ccyin.gt.0.0) then CCY=((REAL(NYA)*SXYI3)-(SXI3*SYI3))/sqrt(ccyin) SIG3 = SQRT(sig3ins) WRITE(6,728) CCY,SIG3 728 FORMAT(1X,'CORRELATION COEFFICIENT FOR DELTA Y DATA = ',F8.3, * ' SIGMA = ',F8.4) Else Write(6,797) ccyin,sig3ins,nya End if END IF IF (NZA.GT.0) THEN cczin = (REAL(NZA)*S2XI4-SXI4**2)*(REAL(NZA)*S2YI4-SYI4**2) SIG4ins = (NZA*S2F4-SF4**2)/(NZA*(NZA-1)) If (sig4ins.gt.0.0.and.cczin.gt.0.0) then CCZ=((REAL(NZA)*SXYI4)-(SXI4*SYI4))/sqrt(cczin) SIG4 = SQRT(sig4ins) WRITE(6,729) CCZ,SIG4 729 FORMAT(1X,'COR. COEF. FOR Z COMP = ',F8.3, * ' SIGMA ON DEL = ',F8.4) Else Write(6,797) cczin,sig4ins,nza End if END IF C 99 Continue Return END SUBROUTINE FUN(NX,X,NP,FR,FT,FP,YC,P,SS,CC,CS,STH,CTH, * ALON,THICK,DFDP1,DFDP2,DFDP3,DFDP4,SPACG,BF, * cdist,imagsus) C C THE PURPOSE OF THIS SUBROUTINE IS TO DETERMINE THE MAGNETIC FIELD C AT A POINT P1 DUE TO A MAGNETIC DIPOLE AT A POINT P2. C C THE SUBROUTINE ALSO CALCULATES THE PARTIAL DERIVATIVES, DFDP(L), C WHICH FORM THE ELEMENTS OF THE MATRIX A C SAVE DIMENSION DFDP1(*),DFDP2(*),DFDP3(*),DFDP4(*) DIMENSION X(3),SS(*),CC(*),CS(*) DIMENSION STH(*),CTH(*),THICK(*),BF(*) DIMENSION ALON(*),P(*),cdist(*) REAL ONE,TWO,THREE C C Mars Radius = R DATA R/3393.5/ DATA PI/3.1415927/ DATA ARC /1.745329E-2/ Data scaling /100./ C ONE=X(1) TWO=X(2) THREE=R+X(3) C C CALCULATE THE MAIN FIELD AT THE OBSERVATION FOR USE C IN CALCULATING THE TOTAL FIELD PARTIAL C C CALL FIDD(ONE,TWO,THREE,1980.,A,B,C,W,0.0) C DC=-C/W C BC=-A/W C EC=+B/W C T=R+X(3) TCR=(90.-X(1))*ARC CTC=COS(TCR) STC=SIN(TCR) FR=0.0 FT=0.0 FP=0.0 YC=0.0 DO 57 I=1,NP DFDP1(I)=0.0 DFDP2(I)=0.0 DFDP3(I)=0.0 DFDP4(I)=0.0 CTD=STH(I) STD=CTH(I) DPH=(X(2)-ALON(I))*ARC CDP=COS(DPH) SDP=SIN(DPH) DD=SS(I) BD=CC(I) CD=CS(I) A1=CTD*CTC+STD*STC*CDP B1=CTC*STD-STC*CTD*CDP C1=STC*SDP A2=-STC*CTD+CTC*STD*CDP B2=-STC*STD-CTC*CTD*CDP C2=CTC*SDP A3=-STD*SDP B3=CTD*SDP C3=CDP E = R*R + T*T - 2.0*R*T*A1 C C THE SQUARE ROOT OF E IS THE DISTANCE BETWEEN THE DIPOLE AT C THE EARTH'S SURFACE AND THE OBSERVATION AT SATELLITE ALTITUDE C C = SQRT(E) cdist(i) = c C = C * E C C 1 degree of latitude on Mars equal approx. 59.2 km C VOL=(59.2**2)*(SPACG**2)*THICK(I) C C Switch between magnetization and susceptibility C If (imagsus.eq.1) then C C Magnetization: The scaling factor takes into account the mu C nought and the conversion from nT to T C C= (scaling*vol)/c C Else if (imagsus.eq.2) then C C Susceptibility C susc = bf(i)/(4.0*pi) c=(susc*vol)/c C End if C F1= T*A1-R F2= -T*B1 F3= T*C1 D1=T-R*A1 D2=-R*A2 D3=-R*A3 D1=3.*D1 D2=3.*D2 D3=3.*D3 F1=F1/E F2=F2/E F3=F3/E D11=C*(D1*F1-A1) D12=C*(D1*F2+B1) D13=C*(D1*F3-C1) D21=C*(D2*F1-A2) D22=C*(D2*F2+B2) D23=C*(D2*F3-C2) D31=C*(D3*F1-A3) D32=C*(D3*F2+B3) D33=C*(D3*F3-C3) C C FR = RADIAL COMPONENT OF FIELD. C FT = COLATITUDE COMPONENT OF FIELD. C FP = LONGITUDE COMPONENT OF FIELD. C DFDP1(I)=DD*D11+BD*D12+CD*D13 DFDP2(I)=DD*D21+BD*D22+CD*D23 DFDP3(I)=DD*D31+BD*D32+CD*D33 C DFDP4(I)=DC*DFDP1(I)+BC*DFDP2(I)+EC*DFDP3(I) FR=P(I)*DFDP1(I) + FR FT=P(I)*DFDP2(I) + FT FP=P(I)*DFDP3(I) + FP C YC=P(I)*DFDP4(I) + YC 57 CONTINUE RETURN END SUBROUTINE RDATA(NX,ALON,NP,SS,CC,CS,CTH,STH,THICK,FR,FT,FP,YC, * P,DFDP1,DFDP2,DFDP3,DFDP4,Y,SPACG,ICOMP,IEND,X,BD, * STDERR,LONGRNG,cdist,imagsus) SAVE C C READS IN THE DATA C DIMENSION X(3),STH(*),CTH(*) DIMENSION SS(*),CC(*),CS(*) DIMENSION DFDP1(*),DFDP2(*),DFDP3(*),DFDP4(*) DIMENSION THICK(*),BD(*) DIMENSION ALON(*),P(*),cdist(*) REAL FLAT,FLON,FALT,FFLD DATA PAD,PAD1,PAD2,PAD3,PAD4 /99999.,99991.,99992.,99993.,99994./ C C READ AN OBSERVATION C 48 READ(55,END=52) FLAT,FLON,FALT,FFLD,STDERR IEND = 0 10 FORMAT(5F10.3) IF (FLAT.EQ.PAD1) THEN ICOMP=1 GO TO 48 ELSE IF (FLAT.EQ.PAD2) THEN ICOMP=2 GO TO 48 ELSE IF (FLAT.EQ.PAD3) THEN ICOMP=3 GO TO 48 ELSE IF (FLAT.EQ.PAD4) THEN ICOMP=4 GO TO 48 ELSE IF (FLAT.EQ.PAD) THEN GO TO 48 END IF IF (STDERR.EQ.0.0.OR.STDERR.EQ.PAD) STDERR=1.0 C X(1)=FLAT X(2)=FLON X(3)=FALT Y = FFLD C IF (LONGRNG.EQ.2) THEN IF (X(2).LT.0.0) X(2)=X(2)+360.0 ELSE IF (LONGRNG.EQ.1) THEN IF (X(2).GT.180.0) X(2) = X(2)-360.0 END IF C C THE THICKNESS IS SET TO A CONSTANT 40 KM HERE C DO 59 K4=1,NP THICK(K4)=40.0 59 CONTINUE C C CALCULATE THE PARTIAL DERIVATIVE RELATING THE OBSERVATION C TO EACH DIPOLE C CALL FUN(nx,X,NP,FR,FT,FP,YC,P,SS,CC,CS,STH,CTH,ALON, * THICK,DFDP1,DFDP2,DFDP3,DFDP4,SPACG,BD, * cdist,imagsus) C RETURN 52 IEND = 1 RETURN END SUBROUTINE SELECT(NX,ALAT,ALON,NP, * THICK,SS,CC,CS,CTH,STH,BD,LONGRNG) SAVE C C SELECTS THE LOCATION OF THE DIPOLES AND CALCULATES C OR READS THEIR PARAMETERS. C C O U T P U T S C ALAT: Latitude of dipole location I C ALON: Longitude of dipole locations I C THICK: Thickness of crust at location I C NP: Number of dipoles C SS: Z DIRECTION COSINE OF DIPOLE LOCATION I C CC: X DIRECTION COSINE OF DIPOLE LOCATION I C CS: Y DIRECTION COSINE OF DIPOLE LOCATION I C CTH: COSINE OF DIPOLE LATITUDE LOCATION I C STH: SIN OF DIPOLE LATITUDE LOCATION I C BD: TOTAL FIELD AT DIPOLE LATITUDE LOCATION I C DIMENSION THICK(*),ALAT(*),ALON(*) DIMENSION SS(*),CC(*),CS(*),CTH(*),STH(*),BD(*) DATA ARC /1.745329E-2/ C WRITE(6,59) 59 FORMAT(1X,'LATITUDE LONGITUDE THICKNESS ') C J = 1 DO 390 I=1,11550 READ(20,82) A1,B1,C1,D1,E1,F1 82 FORMAT(6F10.3) IF (LONGRNG.EQ.2.AND.B1.LT.0.0) B1=360.+B1 ALAT(J) = A1 ALON(J) = B1 THICK(J) = C1 SS(J) = D1 CC(J) = E1 CS(J) = F1 CTH(J) = COS(A1*ARC) STH(J) = SIN(A1*ARC) IF (MOD(J,1000).EQ.0) * WRITE(6,45) A1,B1,THICK(J) 45 FORMAT(1X,3F10.3) J=J+1 390 CONTINUE C NP = J-1 WRITE(6,57) NP 57 FORMAT(/,1X,I6,' DIPOLES HAVE BEEN SELECTED.',//) RETURN END SUBROUTINE FIDD(DLAT,DLONG,ALT1,TM,X,Y,Z,F,DST) SAVE DATA JJ /0/ ,MM/0/ ,NMX /14/ ,NEXT /0/ DATA IO/11/ DATA IDST/1/ DATA LL/1/ ALT = ALT1 CALL FID(IO,JJ,MM,NEXT,IDST,DLAT,DLONG,ALT,TM,DST,NMX,LL,X,Y,Z,F) JJ = 1 RETURN END SUBROUTINE FID (IU,J,MM,NEXT,IDST,DLAT,DLONG,Q1,TM,DST,NMX,L,X,Y, * Z,F) SAVE C*********************************************************************** C J.EQ.0 INPUTS LATITUDE & Q1=ALTITUDE (KM) RELATIVE TO ELLIPSOID C (GEODETIC COORDINATES) C J.EQ.0 OUTPUT FIELD COMPONENTS NORTH,EAST,VERTICAL C IN GEODETIC COORDINATES C C J.NE.0 LAT.&LONG IN SPHERICAL COORDINATES,Q1=GEOCENTRIC RADIUS(KM) C J.NE.0 OUTPUT FLD COMPONENTS NORTH,EAST,VERTICAL IN SPHERICAL COOR C C MM.EQ.0 USE DEFAULT VALUES AE=6378.16,FLAT=298.25 C MM.NE.0 INPUT VALUES FOR AE,FLAT ON FIRST CALL TO FDG C C NEXT.EQ.0 DO NOT READ INPUT VALUES FOR EXTERNAL FIELD PARAMETERS C WHEN L IS GREATER THAN 0 C NEXT.EQ.0 DO NOT EVALUATE EXTERNAL FIELD FROM MODEL C NEXT.NE.0 READ INPUT VALUES FOR EXTERNAL FIELD PARAMETERS WHEN C L GREATER 0 C NEXT.NE.0 EVALUATE EXTERNAL FIELD MODEL C C IDST.EQ.0 DO NOT EVALUATE INDUCED COEFFICIENTS C IDST.EQ.1 EVALUATE INDUCED COEFFICIENTS C C DLAT GEODETIC LATITUDE IN DEGREES WHEN J=0 C GEOCENTRIC LATITUDE IN DEGREES WHEN J=1 C C DLONG LONGITUDE IN DEGREES C C Q1 GEODETIC ALTITUDE (KM) WHEN J=0 C GEOCENTRIC RADIUS (KM) WHEN J=1 C C NMX MAX DEG OF MODEL EVALUATION C C DST DST VALUE C C NMAX MAXIMUM DEGREE AND ORDER OF CONSTANT TERMS OF FIELD MODEL C NMAXT " " " FIRST ORDER TIME " " " C NMAXTT " " " SECOND " " " " " C NMXTTT " " " THIRD " " " " " C C K.EQ.0 FIELD MODEL COEFFICIENTS SCHMIDT NORMALIZED C K.NE.0 FIELD MODEL COEFFICIENTS GAUSS NORMALIZED C C TZERO EPOCH TIME FOR FIELD MODEL COEFFICIENTS C C ABAR MEAN RADIUS USED IN FIELD MODEL POTENTIAL EXPANSION C (DEFAULT = 6371.2) C C MODEXT.EQ.0 NO EXTERNAL FIELD SOLVED WITH MODEL C MODEXT.NE.0 EXTERNAL FIELD SOLVED WITH MODEL C C MODIND.EQ.0 NO INDUCED COEFFS SOLVED WITH MODEL C MODIND.NE.0 INDUCED COEFFS SOLVED WITH MODEL C C L.EQ.0 EVALUATE FIELD C L.GT.0 READ IN FIELD MODEL AND EVALUATE FIELD C L.LE.0 EVALUATE FIELD AT OLD TIME C C*********************************************************************** DIMENSION Q(31,31),TG(31,31) DIMENSION G(31,31),GT(31,31),SHMIT(31,31),AID(33) DIMENSION GTTT(8,8),GTT(31,31) DATA IFRST/0/ DATA AE,FLAT/6378.16,298.25/ DATA TLAST/0./ DATA TABAR/6371.2/ IF (IFRST) 110,100,110 C C C EQUATORIAL EARTH RADIUS AND FLATTENING FACTOR C C USED IN GEODETIC-GEOCENTRIC COORDINATES. C C C C THE MODEL ITSELF IS INDEPENDENT OF THOSE C C PARAMETERS C C C 100 IF (MM.NE.0) READ(IU,101) AE,FLAT 101 FORMAT(1X,2F6.1) WRITE(6,112) 112 FORMAT(//5X,'SUBROUTINE FID HAS BEEN ADJUSTED TO READ AND EVALUATE * LARGE EXTERNAL MODELS') WRITE(6,109) AE,FLAT 109 FORMAT(//5X,'CONSTANTS USED : '/,22X,'EQUATORIAL EARTH RADIUS ', *F8.3/,22X,'EARTH RECIPROCAL FLATTENING ',F6.1//) IFRST=1 FLAT=1. -1./FLAT E1=0. E2=0. E3=0. ALFA1=0. ALFA2=0. ALFA3=0. ALFA4=0. A2=AE**2 A4=AE**4 B2=(AE*FLAT)**2 A2B2=A2*(1.-FLAT**2) A4B4=A4*(1.-FLAT**4) 110 IF (L) 19,1,2 1 IF (TM-TLAST) 17,19,17 2 READ (IU,3) NMAX,NMAXT,NMAXTT,NMXTTT,MODEXT,K,TZERO,ABAR,MODIND, * (AID(I),I=1,12) 3 FORMAT(4I2,2I2,2F6.1,I2,12A4,A2) IF(ABAR.EQ.0.) ABAR=TABAR READ(IU,103) (AID(I),I=14,33) 103 FORMAT(20A4) L=0 WRITE (6,104) (AID(I),I=1,33) 104 FORMAT (25X,12A4,A2/5X,20A4//) WRITE(6,105) NMAX,NMAXT,NMAXTT,NMXTTT,MODEXT,K,TZERO,ABAR,NEXT 105 FORMAT(5X,'FIELD MODEL ORDER (',I2,',',I2,',',I2,',',I2,')'/, *5X,'EXTERNAL FIELD SOLVED WITH MODEL ( 0-NO;.GT.0-DEGREE)',I2/, *5X,'NORMALIZATION (K=0-SCHMIDT ; K.NE.0-GAUSS)',I2/, *5X,'FIELD MODEL EPOCH ',F6.1/, *5X,'FIELD MODEL MEAN RADIUS ',F6.1/, *5X,'EVALUATE EXTERNAL FIELD TO DEGREE',I2//) MAXN=0 TEMP=0. 5 READ (IU,6) N,M,GNM,HNM,GTNM,HTNM,GTTNM,HTTNM 6 FORMAT (2I3,6F11.4) IF (N.LE.0) GO TO 7 MAXN=(MAX0(N,MAXN)) G(N,M)=GNM GT(N,M)=GTNM GTT(N,M)=GTTNM TEMP=AMAX1(TEMP,ABS(GTNM)) IF (M.EQ.1) GO TO 5 G(M-1,N)=HNM GT(M-1,N)=HTNM GTT(M-1,N)=HTTNM GO TO 5 7 IF (NMXTTT.EQ.0) GO TO 107 106 READ(IU,6)N,M,GTTTNM,HTTTNM IF(N.EQ.0) GO TO 107 IF(N.GT.8) STOP 106 GTTT(N,M)=GTTTNM IF (M.EQ.1) GO TO 106 GTTT(M-1,N)=HTTTNM GO TO 106 107 CONTINUE C READ EXTERNAL FIELD IF(MODEXT.NE.0) THEN 30 READ(IU,616) N,M,QNM,SNM 616 FORMAT(2I3,2F12.5) IF(N .LE. 0) GO TO 31 Q(N,M) = QNM IF(M .EQ. 1) GO TO 30 Q(M-1,N) = SNM GO TO 30 END IF 31 CONTINUE IF(MODIND.NE.0.AND.IDST.NE.0) READ(IU,102)ALFA1,ALFA2,ALFA3, * ALFA4 102 FORMAT(6X,4F11.4) CC WRITE(6,8) 8 FORMAT(6H0 N M,6X,1HG,10X,1HH,9X,2HGT,9X,2HHT,8X,3HGTT, * 8X,3HHTT,7X,4HGTTT,7X,4HHTTT//) DO 12 N=2,MAXN DO 12 M=1,N MI=M-1 IF (M.EQ.1) GO TO 10 CC IF (N.GT.NMXTTT) WRITE(6,9) N,M,G(N,M),G(MI,N), CC * GT(N,M),GT(MI,N),GTT(N,M),GTT(MI,N) CC IF (N.LE.NMXTTT) WRITE(6,9) N,M,G(N,M),G(MI,N),GT(N,M), CC * GT(MI,N),GTT(N,M),GTT(MI,N),GTTT(N,M),GTTT(MI,N) 9 FORMAT(2I3,8F11.4) GO TO 12 10 CONTINUE CC IF (N.GT.NMXTTT) WRITE(6,11)N,M,G(N,M),GT(N,M), CC * GTT(N,M) CC IF (N.LE.NMXTTT) WRITE(6,11)N,M,G(N,M),GT(N,M), CC * GTT(N,M),GTTT(N,M) 11 FORMAT(2I3,F11.4,11X,F11.4,11X,F11.4,11X,F11.4) 12 CONTINUE IF (MODEXT.NE.0) THEN C WRITE(6,108) DO 32 N = 2,MODEXT DO 32 M = 1,N IF (M.EQ.1) THEN SNM = 0.0 ELSE SNM = Q(M-1,N) END IF WRITE(6,6) N,M,Q(N,M),SNM 32 CONTINUE END IF CC IF(IDST.NE.0) WRITE(6,111) ALFA1,ALFA2,ALFA3,ALFA4 111 FORMAT(//5X,'INDUCED COEFFS, ',4F10.4,//) C108 FORMAT(//5X,8HEXTFLD, /) 13 FORMAT (1H1) IF (TEMP.EQ.0.) L=-1 14 IF (K.NE.0) GOTO17 SHMIT(1,1)=-1. DO 15 N=2,MAXN SHMIT(N,1)=SHMIT(N-1,1)*FLOAT(2*N-3)/FLOAT(N-1) SHMIT(1,N)=0. JJ=2 DO 15 M=2,N SHMIT(N,M)=SHMIT(N,M-1)* * SQRT(FLOAT((N-M+1)*JJ)/FLOAT(N+M-2)) SHMIT(M-1,N)=SHMIT(N,M) 15 JJ=1 C WRITE(6,300) C300 FORMAT(' FID SHMIT') DO 16 N=2,MAXN DO 16 M=1,N G(N,M)=G(N,M)*SHMIT(N,M) GT(N,M)=GT(N,M)*SHMIT(N,M) GTT(N,M)=GTT(N,M)*SHMIT(N,M) IF (NMXTTT.GT.0.AND.N.LE.8) * GTTT(N,M)=GTTT(N,M)*SHMIT(N,M) IF (M.EQ.1) GO TO 16 G(M-1,N)=G(M-1,N)*SHMIT(M-1,N) GT(M-1,N)=GT(M-1,N)*SHMIT(M-1,N) GTT(M-1,N)=GTT(M-1,N)*SHMIT(M-1,N) IF (NMXTTT.GT.0.AND.N.LE.8) * GTTT(M-1,N)=GTTT(M-1,N)*SHMIT(M-1,N) 16 CONTINUE IF (MODEXT .NE. 0) THEN DO 33 N = 2,MODEXT DO 33 M = 1,N Q(N,M) = Q(N,M)*SHMIT(N,M) IF(M .EQ. 1) GO TO 33 Q(M-1,N) = Q(M-1,N)*SHMIT(M-1,N) 33 CONTINUE END IF C C IF THE TIME AT WHICH THE FIELD IS TO BE CALCULATED DIFFERS C FROM THE EPOCH TIME FOR THE FIELD MODEL COEFFICIENTS, C THE FOLLOWING CODE WILL CALCULATE THE CORRECTION C 17 T=TM-TZERO DO 18 N=1,MAXN DO 18 M=1,N TGX=0. THX=0. IF(M.EQ.1) GO TO 270 IF(N.GT.NMXTTT) GO TO 210 TGX=GTTT(N,M)*T THX=GTTT(M-1,N)*T 210 IF(N.GT.NMAXTT) GO TO 220 TGX=(TGX + GTT(N,M))*T THX= (THX + GTT(M-1,N))*T 220 IF(N.GT.NMAXT) GO TO 230 TGX=(TGX + GT(N,M))*T THX=(THX+GT(M-1,N))*T 230 TGX=TGX+G(N,M) THX=THX+G(M-1,N) TG(N,M)=TGX TG(M-1,N)=THX GO TO 18 270 CONTINUE IF(N.GT.NMXTTT) GO TO 240 TGX=GTTT(N,M)*T 240 IF(N.GT.NMAXTT) GO TO 250 TGX=(TGX+GTT(N,M))*T 250 IF(N.GT.NMAXT) GO TO 260 TGX=(TGX+GT(N,M))*T 260 TGX= TGX+G(N,M) TG(N,M)=TGX 18 CONTINUE TLAST=TM C C END OF SV SECTION C 19 DLATR=DLAT/57.2957795D0 SINLA=SIN(DLATR) RLONG=DLONG/57.2957795D0 CPH=COS(RLONG) SPH=SIN(RLONG) IF (J.EQ.0) GO TO 20 C C Q1 IS GEOCENTRIC RADIUS WHEN J=1 C R=Q1 CT=SINLA GO TO 21 20 SINLA2=SINLA**2 C C Q1 IS GEODETIC ALTITUDE WHEN J=0 C ALT=Q1 C COSLA2=1.-SINLA2 DEN2=A2-A2B2*SINLA2 DEN=SQRT(DEN2) FAC=(((Q1*DEN)+A2)/((Q1*DEN)+B2))**2 CT=SINLA/SQRT(FAC*COSLA2+SINLA2) R=SQRT(Q1*(Q1+2.*DEN)+(A4-A4B4*SINLA2)/DEN2) 21 ST=SQRT(1.-CT**2) NMAX=MIN0(NMX,MAXN) NEXTF=NEXT DSTT=DST IIDST=IDST CALL MAGF(IIDST,ALFA1,ALFA2,ALFA3,ALFA4,DSTT,ST,CT,SPH,CPH, * R,NMAX,BT,BP,BR,B,ABAR,NEXTF,Q,TG) Y=BP F=B IF (J) 22,23,22 22 X=-BT Z=-BR RETURN C C TRANSFORMS FIELD TO GEODETIC DIRECTIONS C 23 SIND=SINLA*ST-SQRT(COSLA2)*CT COSD=SQRT(1.0-SIND**2) X=-BT*COSD-BR*SIND Z=BT*SIND-BR*COSD RETURN END C SUBROUTINE MAGF(IDST,ALFA1,ALFA2,ALFA3,ALFA4,DST,ST,CT,SPH,CPH, * R,NMAX,BT,BP,BR,B,ABAR,NEXT,Q,G) SAVE DIMENSION P(0:31,31),DP(0:31,31),CONST(31,31),SP(31),CP(31), * FM(31),DXDQ(31,31),DXDS(31,31),DYDQ(31,31),DYDS(31,31), * DZDQ(31,31),DZDS(31,31),FN(31) DIMENSION Q(31,31),G(31,31) DATA NCORE/14/ IF (P(1,1).EQ.1.0) GO TO 3 1 P(1,1)=1. DP(1,1)=0. SP(1)=0. CP(1)=1. DO 2 N=2,NMAX FN(N)=N DO 2 M=1,N FM(M)=M-1 2 CONST(N,M)=FLOAT((N-2)**2-(M-1)**2)/FLOAT((2*N-3)*(2*N-5)) 3 SP(2)=SPH CP(2)=CPH DO 4 M=3,NMAX SP(M)=SP(2)*CP(M-1)+CP(2)*SP(M-1) 4 CP(M)=CP(2)*CP(M-1)-SP(2)*SP(M-1) AOR=ABAR/R AR=AOR**2 BTC=0.0 BPC=0.0 BRC=0.0 BC=0.0 BT=0. BP=0. BR=0. IF (IDST.EQ.0) GO TO 12 GBAR=G(2,1) G(2,1)=GBAR + ALFA1*DST E1BAR = Q(2,1) E2BAR = Q(2,2) E3BAR = Q(1,2) Q(2,1) = Q(2,1) + ALFA2*DST Q(2,2) = Q(2,2) + ALFA3*DST Q(1,2) = Q(1,2) + ALFA4*DST 12 CONTINUE DO 8 N=2,NMAX AR=AOR*AR DO 8 M=1,N IF (N-M) 6,5,6 5 P(N,N)=ST*P(N-1,N-1) DP(N,N)=ST*DP(N-1,N-1)+CT*P(N-1,N-1) GO TO 7 6 P(N,M)=CT*P(N-1,M)-CONST(N,M)*P(N-2,M) C C NOTE : CONST(2,1)=0 C DP(N,M)=CT*DP(N-1,M)-ST* * P(N-1,M)-CONST(N,M)*DP(N-2,M) 7 PAR=P(N,M)*AR IF (M.EQ.1) GO TO 9 TEMP=G(N,M)*CP(M)+G(M-1,N)*SP(M) BP=BP-(G(N,M)*SP(M)-G(M-1,N)*CP(M))*FM(M)*PAR GO TO 10 9 TEMP=G(N,M)*CP(M) 10 BT=BT+TEMP*DP(N,M)*AR BR=BR-TEMP*FN(N)*PAR IF (N.GT.NCORE) GO TO 8 BTC=BT BRC=BR BPC=BP 8 CONTINUE BP=BP/ST BPC=BPC/ST BNEXT=SQRT(BT*BT + BP*BP + BR*BR) IF (NEXT.GT.0) THEN MONO = 2 SIND = 0.0 COSD = 1.0 CX = -BT CY = BP CZ = -BR ROA= 1.0/AOR RB= (ROA)**(2*(MONO-2)+1) ROA2= ROA*ROA DO 11 N= MONO,NEXT FNC= N-1 RB= RB*ROA2 DO 11 M= 1,N FMC= M-1 P(N,M)= P(N,M)*RB DP(N,M)= DP(N,M)*RB TEMP= -FNC*P(N,M)*SIND - DP(N,M)*COSD DXDQ(N,M)= TEMP*CP(M) DXDS(N,M)= TEMP*SP(M) TEMP= FMC*P(N,M)/ST DYDQ(N,M)= -TEMP*SP(M) DYDS(N,M)= TEMP*CP(M) TEMP= -FNC*P(N,M)*COSD + DP(N,M)*SIND DZDQ(N,M)= TEMP*CP(M) DZDS(N,M)= TEMP*SP(M) IF (M .EQ. 1) THEN CX= CX + Q(N,M)*DXDQ(N,M) CY= CY + Q(N,M)*DYDQ(N,M) CZ= CZ + Q(N,M)*DZDQ(N,M) ELSE CX= CX + Q(N,M)*DXDQ(N,M) + Q(M-1,N)*DXDS(N,M) CY= CY + Q(N,M)*DYDQ(N,M) + Q(M-1,N)*DYDS(N,M) CZ= CZ + Q(N,M)*DZDQ(N,M) + Q(M-1,N)*DZDS(N,M) END IF 11 CONTINUE BRC = BRC + (-CZ - BR) BPC = BPC + ( CY - BP) BTC = BTC + (-CX - BT) BT = -CX BP = CY BR = -CZ END IF B=SQRT(BT*BT+BP*BP+BR*BR) C C **** B(14 - 30) *** C BC=SQRT(BTC*BTC + BRC*BRC + BPC*BPC) BTC=BT - BTC BPC=BP - BPC BRC=BR - BRC BC= B - BC IF(IDST.EQ.0) RETURN Q(2,1) = E1BAR Q(2,2) = E2BAR Q(1,2) = E3BAR G(2,1)=GBAR RETURN End SUBROUTINE linbcgd(m,n,b,x,itol,tol,itmax,iter,err,sc,ijc) INTEGER ijc(*),iter,itmax,itol,n Real err,tol,b(*),x(*),sc(*),EPS C C MX is the number of observations C NX is the number of dipoles C PARAMETER (EPS=1.d-14,mx=40546,nx=11550) Real xold(nx) INTEGER j Real ak,akden,bk,bkden,bknum,bnrm,xnrm, * znrm,snrm Real p(NX),rp(NX),r(NX),s(mx),w(mx) C iter=0 call dsprstx(sc,ijc,b,r,m,n) do 11 j=1,m s(j) = b(j) 11 Continue err = sderr(m,s) xtemp = 0.0 C Write(70) iter,xtemp,err Print *,'Initial SD on fit = ',err C Write(6,882) 882 Format(/,1x,'Iteration',4x,'SD on fit',4x,'Solution RMS') C Do 932 j=1,n p(j) = r(j) 932 continue 100 if (iter.le.itmax) then iter=iter+1 call dsprsax(sc,ijc,p,w,m,n) aknum = 0.d0 akden = 0.d0 do 12 j=1,n aknum=aknum+r(j)*r(j) 12 continue Do 433 j=1,m akden = akden+w(j)*w(j) 433 continue ak = aknum/akden do 16 j=1,n x(j)=x(j)+ak*p(j) 16 continue Do 434 j=1,m s(j)=s(j)-ak*w(j) 434 continue err = sderr(m,s) xtemp=sderr(n,x) Do 24 j=1,n rp(j) = r(j) 24 Continue Call dsprstx(sc,ijc,s,r,m,n) If (iter.gt.1) then If (err.gt.errold) then Do 432 j=1,n x(j) = xold(j) 432 Continue iter = iter-1 err = errold Print *,'Iteration ending because error has increased' Return End if End if Write(6,714) iter,err,xtemp 714 Format(1x,I5,4x,2E12.3) C Write(70) iter,xtemp,err errold = err Do 431 j=1,n xold(j) = x(j) 431 Continue bknum = 0.d0 bkden = 0.d0 Do 20 j=1,n bknum = bknum + r(j)*r(j) bkden = bkden + rp(j)*rp(j) 20 Continue bk = bknum/bkden Do 22 j=1,n p(j) = r(j)+bk*p(j) 22 Continue if(err.gt.tol) goto 100 endif return END FUNCTION sderr(n,x) Integer n Real x(n),sx sx2 = 0.0 Do 11 i=1,n sx2 = x(i)**2 + sx2 11 Continue sderr = sqrt(sx2/real(n)) return end FUNCTION snrm(n,sx,itol) INTEGER n,itol,i,isamax Real sx(n),snrm if (itol.le.3)then snrm=0. do 11 i=1,n snrm=snrm+sx(i)**2 11 continue snrm=sqrt(snrm) else isamax=1 do 12 i=1,n if(abs(sx(i)).gt.abs(sx(isamax))) isamax=i 12 continue snrm=abs(sx(isamax)) endif return END Subroutine SDATA(it17,incerr,idtype,ci) C C Selects scalar obs. and writes them to an 'esmap-style' C file on units 55 (corresponding to Region 1) C Dimension z(720,357),zerr(720,357) C PAD1=99991.0 PAD=99999.0 icount = 0 ict = 0 ci = 0.0 C WRITE(55) PAD1,PAD1,PAD1,PAD1,PAD1 100 Continue C 558 Continue C If (idtype.eq.2) then Go to 200 Else if (idtype.eq.1) then C C 7F10.2 data C Read(32,901) ((z(i,j),j=1,357),i=1,720) 901 Format(7F10.2) If (incerr.eq.1) then Read(36,901) ((zerr(i,j),j=1,357),i=1,720) End if C alt = 400.0 alat = -89. alon = 0.0 dkn = z(361,1) ster=1.0 If (incerr.eq.1) then ster = zerr(361,1) ci = ci + (dkn/ster)**2 Else ci = ci + dkn**2 End if Write(55) alat,alon,alt,dkn,ster icount = icount + 1 C alat = 89. alon = 0.0 dkn = z(361,357) If (incerr.eq.1) then ster=zerr(361,357) ci = ci + (dkn/ster)**2 Else ci = ci + dkn**2 End if Write(55) alat,alon,alt,dkn,ster icount = icount + 1 C Read(20,82) a1,b1,c1,d1,e1,f1,g1 82 Format(7F10.3) Do 501 n=1,11560 Read(20,82) dlat,dlon,c1,d1,e1,f1,g1 latix = nint((dlat+89.)*2) + 1 lonix = nint((dlon+180.)*2) + 1 If (dlon.eq.180.) lonix = 1 dkn = z(lonix,latix) If (incerr.eq.2) then ci = ci + dkn**2 Else if (incerr.eq.1) then ster = zerr(lonix,latix) ci = ci + (dkn/ster)**2 End if alat = ((latix-1)/2.)-89. alon = ((lonix-1)/2.)-180. icount = icount+1 Write(55) alat,alon,alt,dkn,ster 501 Continue Rewind 20 Go to 90 Else Go to 558 End if C C Binned data C 200 Read(32,end=90) alat,alon,rdk,dkn,ster C alt = rdk - 6371.2 ci = ci + dkn**2 Write(55) alat,alon,alt,dkn,ster icount = icount + 1 C C Return for more data C Go to 200 90 Continue C Print *,'Input data set complete' Print *,icount,' observations read' Write(6,882) ci 882 Format(1x,'Initial cost function = ',E15.4) Write(59) ci,icount C C Indicate that the files are now closed C WRITE(55) PAD,PAD,PAD,PAD,PAD Close (55) Rewind (55) C Return End Subroutine VDATA(it17,idtype) C C Selects vector obs. and writes them to an 'esmap-style' C file on units 55 (corresponding to Region 1) C Pad4=99994.0 Pad3=99993.0 PAD1=99991.0 PAD=99999.0 icount = 0 C WRITE(55) PAD4,PAD4,PAD4,PAD4,PAD4 Do 501 n=1,49300 Read(32,*) alon,alat,alt,br,bt,bp,number,sigma_r, & sigma_t,sigma_p z = -br xyz1=alat+0.5 xyz2=alon+0.5 xyz3=alt+5. If (sigma_r.eq.0.0) go to 501 icount = icount+1 Write(55) xyz1,xyz2,xyz3,z,sigma_r 501 Continue Do 504 n=1,16927 Read(33,*) alon,alat,alt,br,bt,bp,number,sigma_r, & sigma_t z = -br xyz1=alat+0.5 xyz2=alon+0.5 xyz3=alt+5. If (sigma_r.eq.0.0) go to 504 icount = icount+1 Write(55) xyz1,xyz2,xyz3,z,sigma_r 504 Continue Write(55) pad,pad,pad,pad,pad Rewind(32) Rewind(33) Rewind(55) C Go to 59 C Rewind(32) WRITE(55) PAD2,PAD2,PAD2,PAD2,PAD2 Do 502 n=1,49300 Read(32,*) alon,alat,alt,br,bt,bp,number,sigma_r, & sigma_t,sigma_p x = -bt xyz1=alat+0.5 xyz2=alon+0.5 xyz3=alt+5. If (sigma_t.eq.0.0) sigma_t=5.0 icount = icount+1 Write(55) xyz1,xyz2,xyz3,x,sigma_t 502 Continue Write(55) pad,pad,pad,pad,pad C Rewind(32) WRITE(55) PAD3,PAD3,PAD3,PAD3,PAD3 Do 503 n=1,49300 Read(32,*) alon,alat,alt,br,bt,bp,number,sigma_r, & sigma_t,sigma_p x = bp xyz1=alat+0.5 xyz2=alon+0.5 xyz3=alt+5. If (sigma_p.eq.0.0) sigma_p=5.0 icount = icount+1 Write(55) xyz1,xyz2,xyz3,x,sigma_p 503 Continue Write(55) pad,pad,pad,pad,pad C 59 Print *,'Input data set complete' Print *,icount,' observations read' C Return End SUBROUTINE dsprsax(sa,ija,x,b,m,n) INTEGER n,ija(*) Real b(*),sa(*),x(*) INTEGER i,k if (ija(1).ne.m+2) then Print *,'mismatched vector and matrix in sprsax' Stop End if Do 10 i=1,m b(i) = 0.0 10 continue do 12 i=1,m do 11 k=ija(i),ija(i+1)-1 b(i)=b(i)+sa(k)*x(ija(k)) 11 continue 12 continue return END SUBROUTINE dsprstx(sa,ija,x,b,m,n) INTEGER n,ija(*) Real b(*),sa(*),x(*) INTEGER i,j,k if (ija(1).ne.m+2) pause 'mismatched vector and matrix in sprstx' do 10 i=1,n b(i) = 0.0 10 continue do 13 i=1,m do 12 k=ija(i),ija(i+1)-1 j=ija(k) b(j)=b(j)+sa(k)*x(i) 12 continue 13 continue return END