* * edcgnnema.f in /u/purucker/fortran/esmap ********************************************************************* * * 2. NAME DERIVATION - Equivalent Source MAPping - Whole world * 3.4 degree dipole spacing - Conjugate Gradient * technique. Icosahedral tesselation of sphere. * Using DESIGN matrix instead of normal matrix * * 3. FUNCTION - Accumulates the design matrix and the rhs * in sparse storage mode, data (m) > dipoles (n) * 3612 dipoles. X, Y, Z or scalar components * * 4. LANGUAGE - Standard Fortran 77. * * 5. IMPORTANT VARIABLES - * NX: THE MAXIMUM NUMBER OF DIPOLES ALLOWED. * MX: The maximum number of observations allowed * Change as needed in the parameter statement at the * front of esmap2 * * 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 36: Input. Error associated with observation * Unit 40: Output. Dipole parameters (Magnetization or suscep) * * 8. HARDWARE/SOFTWARE ENVIRONMENT - Works on Cray * No special libraries required. * * 9. HISTORY - * * VERSION DATE AUTHOR DESCRIPTION * _______ ____ ______ ___________ * * 1.0 1995 M. Purucker Cray version * ******************************************************************* * Open(11,file='gsfc1283.coef',status='old',access='sequential', & form='formatted') C Open(18,file='orsted_3612.edcgnnema',status='new', & form='unformatted') C Open(20,file='icosfid_3612',status='old',access='sequential', & form='formatted') C C For unformatted data C C Open(32,file='dxz400.ricard1',status='old', C & access='sequential',form='unformatted') C C For formatted data C Open(32,file='obs_model_oersted_03c_indrem.xyzz',status='old', & access='sequential',form='formatted') CC Open(36,file='b1005.pgmg.lesserr.mkmp',status='old', CC & 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 = 1 C CC Print *,'Include (1) or omit (2) an error file?' 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 = 2 C CC Print *,'Enter threshold value in km:' CC Read(5,671) thresh thresh = 1500. C it17 = 2 If (isvdata.eq.1) then Call SDATA(it17,incerr,idtype,ci,icount) Else if (isvdata.eq.2) then Call VDATA(it17,idtype,icount) End if C Call esmap2(ci,imagsus,thresh,icount) C End Subroutine esmap2(ci,imagsus,thresh,mp) C PARAMETER (mx=3612,NX=3612,nmax=1400000) 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) DIMENSION DFDP1(NX),DFDP2(NX),DFDP3(NX),DFDP4(NX) DIMENSION THICK(NX),ALAT(NX),ALON(NX) 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 = 3.4 C C PRINT OUT THE Critical Information C WRITE(6,1) SPACG,LONGRNG 1 FORMAT(/,1X,'Dipole spacing = ',F10.2,/,1x, * 'LONGITUDE RANGE = ',I4) write(6,779) thresh 779 Format(1x,'Dipoles retained out to ',F10.2,' km') C WRITE(6,814) 814 FORMAT(1X,'CONSTANT THICKNESS CRUST OF 40 KM ASSUMED.') C C R = EARTH RADIUS (KM) C R = 6371.2 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 ylat = 0.0 ylon = 0.0 ytime = 1980.0 ydst = 0.0 CALL FIDD(ylat,ylon,R,ytime,A,B,C,W1,ydst) C C READ IN AN OBSERVATION AND CALCULATE THE PARTIALS C WRITE(6,4) 4 FORMAT(1X,'LATITUDE LONGITUDE ALTITUDE OBSERVATION',3X, * 'STD.ERR') C NBA=0 NXA=0 NYA=0 NZA=0 ICOUNT = 0 ijc(1) = mp + 2 k = mp+1 C 5 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 9 C IF (ICOMP.EQ.1.AND.NBA.EQ.0) WRITE(6,62) 62 FORMAT(/,15X,'TOTAL FIELD ANOMALY VALUES') IF (ICOMP.EQ.2.AND.NXA.EQ.0) WRITE(6,63) 63 FORMAT(/,15X,'DELTA X ANOMALY VALUES') IF (ICOMP.EQ.3.AND.NYA.EQ.0) WRITE(6,64) 64 FORMAT(/,15X,'DELTA Y ANOMALY VALUES') IF (ICOMP.EQ.4.AND.NZA.EQ.0) WRITE(6,65) 65 FORMAT(/,15X,'DELTA Z ANOMALY VALUES') C ICOUNT = ICOUNT + 1 If (mod(icount,2000).eq.0.or.icount.eq.1) * Write(6,35) xl(1),xl(2),xl(3),y,stderr 35 FORMAT(1X,5F10.3) C C ACCUMULATE DW, SC, and IJC C IF (ICOMP.EQ.1) THEN C C SCALAR C NBA = NBA+1 DY = Y If (stderr.ne.1.) then DW(icount)=DY*(1./STDERR) DO 14 J=1,np If (cdist(j).le.thresh) then k=k+1 If (k.gt.nmax) Print *,'nmax too small.' sc(k) = dfdp4(j)*(1./stderr) ijc(k) = j end if 14 Continue ijc(icount+1) = k+1 Else DW(icount)=DY DO 217 J=1,np If (cdist(j).le.thresh) then k=k+1 If (k.gt.nmax) Print *,'nmax too small.' sc(k) = dfdp4(j) ijc(k) = j end if 217 Continue ijc(icount+1)=k+1 End if C ELSE IF (ICOMP.EQ.2) THEN C C DELTA X C NXA = NXA+1 DY = -Y If (stderr.ne.1.) then DW(icount)=DY*(1./STDERR) DO 80 J=1,np If (cdist(j).le.thresh) then k=k+1 If (k.gt.nmax) Print *,'nmax too small.' sc(k) = dfdp2(j)*(1./stderr) ijc(k) = j end if 80 Continue ijc(icount+1) = k+1 Else DW(icount)=DY DO 81 J=1,np If (cdist(j).le.thresh) then k=k+1 If (k.gt.nmax) Print *,'nmax too small.' sc(k) = dfdp2(j) ijc(k) = j end if 81 Continue ijc(icount+1)=k+1 End if C ELSE IF (ICOMP.EQ.3) THEN C C DELTA Y C NYA = NYA+1 DY = Y If (stderr.ne.1.) then DW(icount)=DY*(1./STDERR) DO 82 J=1,np If (cdist(j).le.thresh) then k=k+1 If (k.gt.nmax) Print *,'nmax too small.' sc(k) = dfdp3(j)*(1./stderr) ijc(k) = j end if 82 Continue ijc(icount+1) = k+1 Else DW(icount)=DY DO 83 J=1,np If (cdist(j).le.thresh) then k=k+1 If (k.gt.nmax) Print *,'nmax too small.' sc(k) = dfdp3(j) ijc(k) = j end if 83 Continue ijc(icount+1)=k+1 End if C ELSE IF (ICOMP.EQ.4) THEN C C DELTA Z C NZA = NZA+1 DY = -Y If (stderr.ne.1.) then DW(icount)=DY*(1./STDERR) DO 84 J=1,np If (cdist(j).le.thresh) then k=k+1 If (k.gt.nmax) Print *,'nmax too small.' sc(k) = dfdp1(j)*(1./stderr) ijc(k) = j end if 84 Continue ijc(icount+1) = k+1 Else DW(icount)=DY DO 85 J=1,np If (cdist(j).le.thresh) then k=k+1 If (k.gt.nmax) Print *,'nmax too small.' sc(k) = dfdp1(j) ijc(k) = j end if 85 Continue ijc(icount+1)=k+1 End if C END IF C C RETURN FOR MORE DATA C GO TO 5 9 REWIND (55) C ialen = ijc(ijc(1)-1)-1 C WRITE(6,182) NXA,NYA,NZA,NBA,ialen 182 FORMAT(/,1X,I7,' X OBSERVATIONS',/,1X,I7,' Y OBSERVATIONS', * /,1X,I7,' Z OBSERVATIONS',/,1X,I7,' SCALAR OBSERVATIONS', * /,1x,'Length of ijc and sc vectors = ',I13) C C Save the results to disk C Write(18) mp Write(18) np Write(18) sc Write(18) ijc Write(18) dw 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 DATA R/6371.2/ DATA PI/3.1415927/ DATA ARC /1.745329E-2/ Data scaling /100./ Data ytime /1980.0/ Data ydst /0.0/ 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 CALL FIDD(ONE,TWO,THREE,ytime,A,B,C,W,ydst) DC=-C/W BC=-A/W 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 VOL=(111.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 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 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 Z X ', * ' Y TOTAL FIELD',/) C J = 1 DO 390 I=1,3612 READ(20,82) A1,B1,C1,D1,E1,F1,G1 82 FORMAT(7F10.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 BD(J) = G1 CTH(J) = COS(A1*ARC) STH(J) = SIN(A1*ARC) IF (MOD(J,1000).EQ.0) * WRITE(6,45) A1,B1,THICK(J),D1,E1,F1,G1 45 FORMAT(1X,7F10.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 SDATA(it17,incerr,idtype,ci,icount) C C Selects scalar obs. and writes them to an 'esmap-style' C file on units 55 (corresponding to Region 1) C PAD1=99991.0 PAD=99999.0 icount = 0 ict = 0 ci = 0.0 alt = 400.0 ster = 1.0 C WRITE(55) PAD1,PAD1,PAD1,PAD1,PAD1 C 200 Read(32,*,end=90) alat,alon,dobs,dmod C dkn = dobs-dmod 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,icount) C C Selects vector obs. and writes them to an 'esmap-style' C file on units 55 (corresponding to Region 1) C Dimension z(721,361),x(721,361) C Pad4=99994.0 Pad3=99993.0 Pad2=99992.0 PAD1=99991.0 PAD=99999.0 icount = 0 ict = 0 alt = 400. ster = 1.0 C WRITE(55) PAD2,PAD2,PAD2,PAD2,PAD2 C 558 Continue C If (idtype.eq.2) then Go to 200 Else if (idtype.eq.1) then C C Gridded binary data C Read(32) x,z C Do 501 n=1,11562 Read(20,82) dlat,dlon,c1,d1,e1,f1,g1 82 Format(7F10.3) latix = nint((dlat+90.)*2) + 1 lonix = nint((dlon+180.)*2) + 1 If (dlon.eq.180.) lonix = 1 dkn = x(lonix,latix) alat = ((latix-1)/2.)-90. alon = ((lonix-1)/2.)-180. icount = icount+1 If (alat.eq.90.0) alat = 89.9 If (alat.eq.-90.0) alat = -89.9 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 If (alat.eq.pad2) go to 200 If (alat.eq.pad) go to 90 icount = icount + 1 C alt = rdk - 6371.2 C C Determine which regions the data falls in C C To be determined C Return for more data C Go to 200 90 Continue Write(55) pad,pad,pad,pad,pad C WRITE(55) PAD4,PAD4,PAD4,PAD4,PAD4 C 758 Continue If (idtype.eq.2) then Go to 209 Else if (idtype.eq.1) then C C gridded binary data C Do 701 n=1,11562 Read(20,82) dlat,dlon,c1,d1,e1,f1,g1 latix = nint((dlat+90.)*2) + 1 lonix = nint((dlon+180.)*2) + 1 If (dlon.eq.180.) lonix = 1 dkn = z(lonix,latix) alat = ((latix-1)/2.)-90. alon = ((lonix-1)/2.)-180. icount = icount+1 If (alat.eq.90.0) alat = 89.9 If (alat.eq.-90.0) alat = -89.9 Write(55) alat,alon,alt,dkn,ster 701 Continue Rewind 20 Go to 823 End if C C Binned data C 209 Read(32,end=823) alat,alon,rdk,dkn,ster If (alat.eq.pad4) go to 209 If (alat.eq.pad) Go to 823 icount = icount + 1 C alt = rdk - 6371.2 C C Determine which regions the data falls in C C To be determined C C Return for more data C Go to 209 823 Continue Write(55) pad,pad,pad,pad,pad Close (55) Rewind (55) C Print *,'Input data set complete' Print *,icount,' observations read' C Return End