PROGRAM DBCOMP                                                    00001020
C
C     dbcomp_mars_1500.f in /tmp6/purucker/mgs_science_mag
C      Calculates the Br field on a regular 1 degree
C      grid at 200 km from a equivalent source magnetization
C      distribution. The Br field is calculated only at 
C      or near locations where there are observations.
C                                                                               
      parameter (ntot=11562)
      REAL P(ntot),CC(ntot),CS(ntot),SS(ntot),bfield(ntot)                     
      REAL FLAT(ntot),FLON(ntot),CTH(ntot),STH(ntot)   
      Real br(361,181)                        
C                                                                               
      Open(40,
     &  file='br_spoabgt3_br_it87_mars_rdipole_200km_1500cf.xyz',
     &     status='new',access='sequential',form='formatted')
C
      PI = 3.1415927                                                            
      arc = 0.0174533
      r = 3393.5
      cutoff= 1500.
      t3 = 3.0
      t2 = 2.0
      scaling = 100.0
      pad = 99999.0
C
      Do 888 m=1,361
        Do 889 n=1,181
          br(m,n)=pad
 889    continue
 888  continue
C
      xyz3=200.
C
C     READ IN THE DIPOLE SET AND ASSOCIATED DATA SETS                           
C        COMPUTE VOL. FROM DIP-SPACE.                                   00000340
C                                                                       00000350
      itot = 0
      NP= 11550                                                                 
      DSP = 1.89                                                                
      VOL= (59.2**2)*(DSP**2)*40.0  
C                                                                               
      Open(15,
     & file='mars.spoabgt3.br.it87.cc.978.radial.dipole.xyz',
     &        status='old',access='sequential',form='formatted')
      Do 18 i=1,np
      READ(15,92) FLON(I),FLAT(I),P(I),ss(i),cc(i),cs(i)
 92   Format(2F10.2,F15.4,3F10.3)
 18   CONTINUE                                                                  
C                                                                               
      T=R+XYZ3                                     
C                                                   
C     CALCULATE FIELD at observation locations and nearby                      
C                                                                               
      Open(30,file='combined_spo_ab.dat',form='formatted',status='old')
C                                                                               
C     ITERATE OVER LATITUDES                                                    
C                                                                               
 33   Read(30,*,end=99) xyz2re,xyz1re,i10,r10,r11,r12,i11,r13
      If (r13.eq.0.0) go to 33
      itot = itot+1
      If (mod(itot,500).eq.0) print *,itot
C
      Do 55 m1=-2,2,1
        do 56 m2=-2,2,1   
          ix1=int(xyz1re-m1)+91
          xyz1 = xyz1re-real(m1)
          If (xyz2re.ge.180) xyz2re=xyz2re-360.
          ix2=int(xyz2re-m2)+181
          If (ix2.lt.1) ix2=ix2+360
          If (ix2.gt.361) ix2=ix2-360
          xyz2 = xyz2re-real(m2)
          If (xyz2.lt.-180.0) xyz2=xyz2+360.
          If (xyz2.ge.180.0) xyz2=xyz2-360.
          If (br(ix2,ix1).ne.pad) go to 56
        
      TCR=(90.-XYZ1)*ARC                             
      CTC=COS(TCR)                                      
      STC=SIN(TCR)                                       
C                                                                               
      FR = 0.0                                                          
      FT = 0.0                                                          
      FP = 0.0                                                          
C                                                                               
C     SUM THE CONTRIBUTIONS FROM THE INDIVIDUAL DIPOLES                         
C                                                                               
      DO 57 I=1,NP                                         
        CTD=sin(flat(I)*arc)
        STD=cos(flat(i)*arc)
        DPH=(XYZ2-FLON(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 = SQRT(E)
        If (c.gt.1500.0) go to 57
        C = C * E
        C= (scaling*vol)/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=DD*D11+BD*D12+CD*D13
        DFDP2=DD*D21+BD*D22+CD*D23
        DFDP3=DD*D31+BD*D32+CD*D33
 
        FR=P(I)*DFDP1 + FR
        FT=P(I)*DFDP2 + FT
        FP=P(I)*DFDP3 + FP
   57   CONTINUE 
        br(ix2,ix1)=fr
   56  Continue
   55 Continue
C
      Go to 33
  99  continue 
      Do 81 m=1,361
        Do 82 n=1,181
          xyz1=real(n-91)
          xyz2=-180.+(m-1)
          Write(40,883) xyz2,xyz1,br(m,n)  
 883      Format(3F10.2)
 82     continue
 81   Continue
C                                                                               
      End