Program dbind C C Beginning with the magnetic crustal thickness model icosfid2.f2c C or a modified crustal thickness model like icosfid2.f2c_namod C C The resulting thickness is then used with susceptibilities for C the continents of 0.035 SI and for the oceans of 0.04 SI. C Produces total field output at the locations of the icosahedal basis C Parameter (ntotl=11562) REAL P(ntotl),CC(ntotl),CS(ntotl),SS(ntotl) REAL FLAT(ntotl),FLON(ntotl),CTH(ntotl),STH(ntotl) REAL cr(ntotl) Real sd(ntotl),mc(ntotl) Integer kflag(ntotl),kqz(ntotl) Character*40 gfield,dfile1,ofile,ozfile,ofiles C DATA ARC,R,T3,T2,YC/.0174533,6371.2,3.0,2.0,0.0/ PI = 3.1415927 xyz3 = 400.0 C C Induced field direction file C Open(11,file='gsfc1283.coef',status='old', & access='sequential',form='formatted') zero=0.0 tm=1980. Call Fidd(zero,zero,R,tm,a,b,c,w,zero) C C Output files C Open(40,file='tfsemm0+mf.allmap1am.xyz',status='new', & form='formatted') Open(42,file='tfsemm0+mf.allmap1am_360.xyz',status='new', & form='formatted') C C READ IN THE DIPOLE SET AND ASSOCIATED DATA SETS C COMPUTE VOL. FROM DIP-SPACE. C NP=11562 C C DSP is dipole spacing in degrees C DSP = 1.89 C C Vol is in cubic km C VOL= 494572.5*DSP**2 C CC Print *,'Is this a SEMM-0 (0) or dipole (1) set:' CC Read(5,881) isd isd = 0 881 Format(I1) C C Input file of magnetic crustal thickness C Open(15,file='icosfid2.f2c_namod',form='formatted') i=1 59 Read(15,558,End=69) Flat(i),Flon(i),ss(i),cc(i),cs(i), & f,cr(i),sd(i),kflag(i),kqz(i),mc(i) 558 Format(2F9.2,3F7.3,F7.0,2F6.1,2I3,F7.1) th = flat(i) sth(i) = sin(th*arc) cth(i) = cos(th*arc) If (kflag(i).ge.8) then susc = 0.04 Else if (kflag(i).le.3) then susc = 0.035 End if P(i) = (F*mc(i)*susc)/(4.*pi*40.) i = i+1 Go to 59 69 Mp = i-1 Print *,mp,' dipoles' C T=R+XYZ3 C C CALCULATE FIELD C DO 100 IC1=1,MP XYZ1 = flat(ic1) If (xyz1.eq.-90.0) xyz1=-89.9 If (xyz1.eq.90.0) xyz1=89.9 TCR=(90.-XYZ1)*ARC CTC=COS(TCR) STC=SIN(TCR) C XYZ2 = flon(ic1) YCR = 0.0 YCT = 0.0 YCP = 0.0 C C SUM THE CONTRIBUTIONS FROM THE INDIVIDUAL DIPOLES C DO 57 I=1,NP DPH=(XYZ2-FLON(I))*ARC CDP = COS(DPH) SDP = SIN(DPH) A1=STH(I)*CTC+CTH(I)*STC*CDP B1=CTC*CTH(I)-STC*STH(I)*CDP C1=STC*SDP A2=-STC*STH(I)+CTC*CTH(I)*CDP B2=-STC*CTH(I)-CTC*STH(I)*CDP A3=-CTH(I)*SDP B3=STH(I)*SDP E=R*R+T*T-T2*T*R*A1 C=SQRT(E) G=(T*(SS(I)*A1-CC(I)*B1+CS(I)*C1)-SS(I)*R)*T3/E GRS=G*R+SS(I) FR=A1*GRS-G*T-CC(I)*B1+CS(I)*C1 FT=A2*GRS-CC(I)*B2+CS(I)*CTC*SDP FP=A3*GRS-CC(I)*B3+CS(I)*CDP DFDP1 = -P(I)/(C*E) YCR=YCR+FR*DFDP1 YCT=YCT+FT*DFDP1 YCP=YCP+FP*DFDP1 57 CONTINUE YCR=YCR*VOL YCT=YCT*VOL YCP=YCP*VOL Call Fidd(xyz1,xyz2,t,tm,a,b,c,w,zero) dc=-c/w bc=-a/w ec= b/w tf=ycr*dc+yct*bc+ycp*ec If (mod(ic1,1000).eq.1) Print *,ic1,' complete' Write(40,85) flat(ic1),flon(ic1),tf qlon=mod(360.+flon(ic1),360.) Write(42,85) flat(ic1),qlon,tf 85 Format(3F9.2) 100 CONTINUE C End