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