C Program does statistical analysis on a table of Mars rock data.  Produces   
C  number density and areal coverage functions versus rock diameter for       
C  rock population at the Viking Lander 1 or Lander 2 sites.                  
                                                                              
                                                                              
C Usage                                                                       
C =====                                                                       
C                    rockana                                                  
C                                                                             
C   There are no arguments.                                                   
C                                                                             
C   To compile and run under Solaris, rename the source file to rockana.f     
C   and type                                                                  
C                      f77 -u -c rockana.f                                    
C                      f77 -o rockana rockana.o                               
C                                                                             
C   Move the data file into the same directory.  Its name should be the       
C   first (and only) line in file  rockana.in  (also in the same directory)   
C                                                                             
C   Then type                                                                 
C                      rockana                                                
                                                                              
                                                                              
C Input File                                                                  
C ==========                                                                  
C   The input file is  rockana.in  and it contains only the name of the file  
C   of input rock data.  If the file name is not VL0ROCKS.TAB, VL1ROCKS.TAB,  
C   or VL2ROCKS.TAB, the program will fail at about line 290.  You must add   
C   a new option to the IF statement giving the areas of the drift, far       
C   field, and near field surfaces.                                           
                                                                              
                                                                              
C Output Files                                                                
C ============                                                                
C   There are four output files.  The file  rockana.out  is the most readable 
C   and probably the most useful.  For each surface type and burial state it  
C   shows:                                                                    
C         (a) the surface area covered by rocks in each size range, where     
C             rocks have been binned according to their average diameter      
C         (b) the average axial ratio of rocks in each bin (computed for      
C             each pair of axes)                                              
C   For each burial state it repeats the above by summing over surface        
C   types (drift material, far field, and near field).  For each surface      
C   type it sums over burial states.  And, finally, it gives both the         
C   raw number distribution and area covered as a function of rock diameter   
C   as well as normalized values (per square meter) after adjusting far       
C   field statistics to have the same distribution of small rocks as were     
C   found in the near field.                                                  
C                                                                             
C   rockana.tab  is a variation on the format above which allows simple       
C   conversion to PDS TABLE objects.                                          
C                                                                             
C   rockana.arrays  is (a) and (b) above for each surface type and burial     
C   state; there have been no accumulations over surface types or burial      
C   states.                                                                   
C                                                                             
C   rockana.db  is an expansion of the original data file (VLnROCKS.TAB)      
C   showing the individual axial ratios and areas for each rock.              
                                                                              
                                                                              
C Advisory                                                                    
C ========                                                                    
C   This program is not guaranteed to do anything.  After limited testing,    
C   it executed under both SunOS and Solaris, and it appeared to reproduce    
C   results reported by Moore and Keller.  It was designed to do a very       
C   simple task initially and grew, without much regard for software          
C   aesthetics, into a much larger package.                                   
                                                                              
                                                                              
                                                                              
      integer*4       ndim             !Max number of records in input file   
      parameter       ( ndim = 500 )                                          
      integer*4       ndbins           !Number of diameter sorting bins       
      parameter       ( ndbins = 18 )                                         
      integer*4       mdbins                                                  
      parameter       ( mdbins = ndbins - 2 )                                 
      real*4          dbins(mdbins)    !Boundaries between bins (m);          
                                       !first bin is junk                     
                                       ! See initialization below.            
                                                                              
      character*4     id(ndim)         !Rock ID                               
      character*1     ty(ndim)         !Surface type (N,F,D, or X)            
      real*4          x(ndim)          !X position (m)                        
      real*4          y(ndim)          !Y position (m)                        
      real*4          z(ndim)          !Z position (m)                        
      real*4          l(ndim)          !Length (m)                            
      real*4          w(ndim)          !Width (m)                             
      real*4          h(ndim)          !Height (m)                            
      real*4          h2(ndim)         !Adjusted height (m)                   
      real*4          ang(ndim)        !Azimuth (deg)                         
      character*1     bur(ndim)        !Burial state (P,H,B, or U)            
      character*24    notes(ndim)      !Name, other notes                     
      character*80    filename         !Input data file name                  
      integer*4       irec             !Counter for reading records           
      integer*4       nrec             !Number of records in input file       
                                                                              
      real*4          dlw(ndim)        !Effective diameter from length        
                                       ! and width (m)                        
      real*4          dwh(ndim)        !Effective diameter from width         
                                       ! and height (m)                       
      real*4          dlh(ndim)        !Effective diameter from length        
                                       ! and height (m)                       
                                                                              
      real*4          arlw(ndim)       !L-W axial ratio                       
      real*4          arwh(ndim)       !W-H axial ratio                       
      real*4          arlh(ndim)       !L-H axial ratio                       
                                                                              
      real*4          area(ndim)       !L*W area                              
                                                                              
C In the following arrays the indices are:                                    
C          1.  sorting by average diameter                                    
C          2.  surface type                                                   
C                   1 - drift                                                 
C                   2 - far field (but not in drift material)                 
C                   3 - near field                                            
C                   4 - excluded                                              
C          3.  burial state:                                                  
C                   1 - perched                                               
C                   2 - partially buried                                      
C                   3 - buried                                                
C                   4 - unknown                                               
C          4.  quantity                                                       
C                   1 - number (count or quantity^0)                          
C                   2 - quantity                                              
C                   3 - quantity^2                                            
                                                                              
      real*4          srlw(ndbins,4,4,3)    !Cumulative L-W axial ratio       
                                            ! based on dlw                    
      real*4          srwh(ndbins,4,4,3)    !Cumulative W-H axial ratio       
                                            ! based on dwh                    
      real*4          srlh(ndbins,4,4,3)    !Cumulative L-H axial ratio       
                                            !based on dlh                     
                                                                              
      real*4          av_arlw(ndbins,4,4,3) !Cumulative L-W axial ratio       
                                            ! based on dlw                    
      real*4          av_arwh(ndbins,4,4,3) !Cumulative W-H axial ratio       
                                            ! based on dwh                    
      real*4          av_arlh(ndbins,4,4,3) !Cumulative L-H axial ratio       
                                            ! based on dlh                    
                                                                              
      real*4          sarea(ndbins,4,4,2)   !Cumulative area                  
                                                                              
                                                                              
C In the following arrays the indices are same except the first is omitted    
                                                                              
      real*4          trlw(4,4,3)        !Avg/SD/N L-W ratio over all sizes   
      real*4          trwh(4,4,3)        !Avg/SD/N W-H ratio over all sizes   
      real*4          trlh(4,4,3)        !Avg/SD/N L-H ratio over all sizes   
      real*4          tarea(4,4,2)       !Summed area/number over all sizes   
                                                                              
                                                                              
C Resume miscellaneous declarations                                           
                                                                              
      integer*4       nsrlw              !Integer version of srlw count       
      integer*4       nsrwh              !Integer version of srwh count       
      integer*4       nsrlh              !Integer version of srlh count       
      integer*4       nsarea             !Integer version of sarea count      
                                                                              
      integer*4       ntrlw                                                   
      integer*4       ntrwh                                                   
      integer*4       ntrlh                                                   
      integer*4       ntarea                                                  
                                                                              
      integer*4       i0                                                      
      integer*4       i1, i1a, i1b, i1c                                       
      integer*4       i2, i3, i4                                              
      real*4          pi                 !Constant                            
      real*4          rt2                !Constant                            
                                                                              
      integer*4       i                  !Miscellaneous index                 
                                                                              
      character*1     lparen, rparen, comma, dquote                           
      character*3     plusminus                                               
                                                                              
      character*14    avg_diam(ndbins)                                        
      character*17    surf_type(4)                                            
      character*23    bur_state(4)                                            
                                                                              
      real*4          surface_area(3)    !Area within each surface type       
      real*4          SA                 !Surface area used for normalization 
                                                                              
      real*4          n_per_bin(ndbins,3)                                     
      real*4          n_cum(3)                                                
      real*4          area_per_bin(ndbins,3)                                  
      real*4          area_cum(3)                                             
                                                                              
      real*4          d_thresh           !Cut-off diameter for filling far    
                                         ! field with near field statistics   
      integer*4       n_d_thresh         !Index in dbins for threshold        
                                                                              
      real*4          ff_interrock_area  !Area in far field not covered       
                                         !by rocks larger than d_thresh       
      real*4          nf_interrock_area  !Area in near field not covered      
                                         !by rocks larger than d_thresh       
      real*4          factor             !ratio of FF to NF interrock areas   
                                                                              
      real*4          comp(ndbins,4)     !composite number/area dist function 
      real*4          comp_cum(4)        !cumulative function                 
                                                                              
      real*4          ar_surf_av(3,ndbins,2,3)                                
      real*4          ar_sarea(ndbins,4,2)                                    
                                                                              
      real*4          t_ar_surf_av(3,3)                                       
      real*4          t_ar_sarea(2)                                           
                                                                              
      integer*4       lw_n                                                    
      integer*4       wh_n                                                    
      integer*4       lh_n                                                    
                                                                              
      real*4          lw_av                                                   
      real*4          wh_av                                                   
      real*4          lh_av                                                   
                                                                              
      real*4          lw_sd                                                   
      real*4          wh_sd                                                   
      real*4          lh_sd                                                   
                                                                              
                                                                              
C Initializations                                                             
                                                                              
      pi  = 4.*atan(1.0)                                                      
      rt2 = sqrt(2.)                                                          
                                                                              
      dquote = char(34)                                                       
      lparen = char(40)                                                       
      rparen = char(41)                                                       
      comma  = char(44)                                                       
      plusminus = "+/-"                                                       
                                                                              
      d_thresh = 0.25                                                         
                                                                              
c      do i = 1,ndbins-2             !Equal bins of 0.05 from 0.05 (RAS)      
c         dbins(i) = 0.05*i                                                   
c      end do                                                                 
                                                                              
      dbins(1) = 0.25/32             !RT2 bins starting from 1/128 m (M&K)    
      do i = 2,16                                                             
         dbins(i) = rt2*dbins(i-1)                                            
      end do                                                                  
                                                                              
      do i=1,ndbins-3                                                         
         write(avg_diam(i+2),'(f5.3,"<=D<",f5.3)') dbins(i),dbins(i+1)        
      end do                                                                  
      avg_diam(1) = "   undefined  "                                          
      write(avg_diam(000002),'("       D<",f5.3)') dbins(1)                   
      write(avg_diam(ndbins),'(f5.3,"<=D      ")') dbins(ndbins-2)            
                                                                              
      surf_type(1) = "DRIFT AREA (I2=1)"                                      
      surf_type(2) = "FAR  FIELD (I2=2)"                                      
      surf_type(3) = "NEAR FIELD (I2=3)"                                      
      surf_type(4) = " EXCLUDED  (I2=4)"                                      
                                                                              
      bur_state(1) = "     PERCHED     (I3=1)"                                
      bur_state(2) = "PARTIALLY BURIED (I3=2)"                                
      bur_state(3) = "     BURIED      (I3=3)"                                
      bur_state(4) = "     UNKNOWN     (I3=4)"                                
                                                                              
      open (unit   = 16,                                                      
     *      file   = 'rockana.out',                                           
     *      status = 'unknown')                                               
                                                                              
      open (unit   = 17,                                                      
     *      file   = 'rockana.db',                                            
     *      status = 'unknown')                                               
                                                                              
      open (unit   = 18,                                                      
     *      file   = 'rockana.arrays',                                        
     *      status = 'unknown')                                               
                                                                              
      open (unit   = 19,                                                      
     *      file   = 'rockana.tab',                                           
     *      status = 'unknown')                                               
                                                                              
      open (unit   = 30,                                                      
     *      file   = 'rockana.in',                                            
     *      status = 'old')                                                   
                                                                              
      read (30,'(a)') filename                                                
      write(16,'("Input from data file ",(a))') filename(01:40)               
                                                                              
      close (30)                                                              
                                                                              
      open (unit   = 31,                                                      
     *      file   = filename,                                                
     *      status = 'old')                                                   
                                                                              
                                                                              
C Define surface areas for each surface type per Henry Moore message of       
C 1996-04-05.  If you want to use this program with other input files,        
C you will have to add to the IF statement below.  Note that                  
C VL0ROCKS.TAB is VL1ROCKS.TAB with two outcrops omitted.                     
                                                                              
      if      (filename(01:12) .eq. "VL0ROCKS.TAB") then                      
         surface_area(1) = 15.27                         !Drift               
         surface_area(2) = 57.85                         !Far Field with      
                                                         ! outcrop area       
C         surface_area(2) = 57.85 - 3.7643                !Far Field less     
                                                          ! outcrops          
         surface_area(3) = 10.48                         !Near Field          
      else if (filename(01:12) .eq. "VL1ROCKS.TAB") then                      
         surface_area(1) = 15.27                         !Drift               
         surface_area(2) = 57.85                         !Far Field           
         surface_area(3) = 10.48                         !Near Field          
      else if (filename(01:12) .eq. "VL2ROCKS.TAB") then                      
         surface_area(1) = 28.09                         !Drift               
         surface_area(2) = 42.83                         !Far Field           
         surface_area(3) = 12.75                         !Near Field          
      else                                                                    
         write(6,'("Unknown file ... cannot get surface areas")')             
         stop                                                                 
      end if                                                                  
                                                                              
                                                                              
C Find threshold bin                                                          
                                                                              
      n_d_thresh = ndbins                                                     
      do i=ndbins-2,1,-1                                                      
         if (int(100.*d_thresh+0.001) .lt.                                    
     *       int(100.*dbins(i)+0.001)) n_d_thresh = i+1                       
      end do                                                                  
                                                                              
      write(16,'("Requested diameter threshold: ",f10.4,/,                    
     *           "Assigned diameter threshold:  ",f10.4,/,                    
     *           "Threshold bin:                ",i10)')                      
     *   d_thresh,dbins(n_d_thresh-2),n_d_thresh                              
                                                                              
C Enter loop to read data records; calculate                                  
C                adjusted height                                              
C                average diameter (3 ways)                                    
C                axial ratios                                                 
C                area                                                         
                                                                              
      irec = 0                                                                
   10 continue                                                                
      irec = irec + 1                                                         
      read (31,'(1x,a4,3x,a1,2x,6(f5.2,1x),f5.0,2x,a1,3x,a21)',end=20)        
     * id(irec),ty(irec),x(irec),y(irec),z(irec),                             
     * l(irec),w(irec),h(irec),ang(irec),bur(irec),notes(irec)                
                                                                              
      if (bur(irec) .eq. 'P') then                                            
        h2(irec) = h(irec)                                                    
      else if (bur(irec) .eq. 'H') then                                       
        h2(irec) = 2*h(irec)                                                  
      else                                                                    
        h2(irec) = -999.                                                      
      end if                                                                  
                                                                              
      dlw(irec) = sqrt(l(irec)*w(irec))                                       
      dwh(irec) = sqrt(w(irec)*h2(irec))                                      
      dlh(irec) = sqrt(l(irec)*h2(irec))                                      
                                                                              
      arlw(irec) = min(l(irec),w(irec))/max(l(irec),w(irec))                  
      arwh(irec) = min(w(irec),h2(irec))/max(w(irec),h2(irec))                
      arlh(irec) = min(l(irec),h2(irec))/max(l(irec),h2(irec))                
                                                                              
      area(irec) = l(irec)*w(irec)*pi/4.                                      
                                                                              
      if (irec .eq. 1) write(17,'("IREC  ID TY    X     Y     Z    ",         
     *     " L     W     H   ANG BUR   H2    DLW   DWH   DLH ",               
     *     " ARLW  ARWH  ARLH   AREA ")')                                     
      write(17,'(i4,1x,a4,1x,a1,6f6.2,f6.0,1x,a1,1x,7f6.2,f8.4)')             
     *  irec,id(irec),ty(irec),                                               
     *  x(irec),y(irec),z(irec),l(irec),w(irec),h(irec),                      
     *  ang(irec),bur(irec),h2(irec),dlw(irec),dwh(irec),dlh(irec),           
     *  arlw(irec),arwh(irec),arlh(irec),area(irec)                           
                                                                              
      go to 10                                                                
                                                                              
   20 continue                                                                
      nrec = irec - 1                                                         
                                                                              
      write(16,'("Read ",i5," input records ...")') nrec                      
                                                                              
C Calculate axial ratios                                                      
                                                                              
   30 continue                                                                
                                                                              
      do i1=1,ndbins                        !Rock size index                  
        do i2=1,4                      !Type index                            
          do i3=1,4                    !Burial index                          
            do i4=1,3                  !count,sum vals, sum val^2             
              srlw(i1,i2,i3,i4)    = 0.                                       
              srwh(i1,i2,i3,i4)    = 0.                                       
              srlh(i1,i2,i3,i4)    = 0.                                       
              av_arlw(i1,i2,i3,i4) = 0.                                       
              av_arwh(i1,i2,i3,i4) = 0.                                       
              av_arlh(i1,i2,i3,i4) = 0.                                       
              if (i4 .le. 2) sarea(i1,i2,i3,i4) = 0.                          
            end do                                                            
          end do                                                              
        end do                                                                
      end do                                                                  
                                                                              
      do 40 irec = 1,nrec                                                     
                                                                              
        i1a = 1                                                               
        if ((dlw(irec) .ge. 0.) .and. (dlw(irec) .lt. dbins(1))) i1a = 2      
        do i = 3,ndbins-1                                                     
           if ((dlw(irec) .ge. dbins(i-2)) .and.                              
     *         (dlw(irec) .lt. dbins(i-1))) i1a = i                           
        end do                                                                
        if (dlw(irec) .ge. dbins(ndbins-2)) i1a = ndbins                      
                                                                              
        i1b = 1                                                               
        if ((dwh(irec) .ge. 0.) .and. (dwh(irec) .lt. dbins(1))) i1b = 2      
        do i = 3,ndbins-1                                                     
           if ((dwh(irec) .ge. dbins(i-2)) .and.                              
     *         (dwh(irec) .lt. dbins(i-1))) i1b = i                           
        end do                                                                
        if (dwh(irec) .ge. dbins(ndbins-2)) i1b = ndbins                      
                                                                              
        i1c = 1                                                               
        if ((dlh(irec) .ge. 0.) .and. (dlh(irec) .lt. dbins(1))) i1c = 2      
        do i = 3,ndbins-1                                                     
           if ((dlh(irec) .ge. dbins(i-2)) .and.                              
     *         (dlh(irec) .lt. dbins(i-1))) i1c = i                           
        end do                                                                
        if (dlh(irec) .ge. dbins(ndbins-2)) i1c = ndbins                      
                                                                              
                                                                              
        if      (ty(irec) .eq. "D") then                                      
           i2 = 1                                                             
        else if ((ty(irec) .eq. "F") .or. (ty(irec) .eq. "O")) then           
           i2 = 2                                                             
        else if (ty(irec) .eq. "N") then                                      
           i2 = 3                                                             
        else if (ty(irec) .eq. "X") then                                      
           i2 = 4                                                             
        else                                                                  
           write(6,'("Unexpected rock type ty(",i3.3,")=",a1)')               
     *       irec,ty(irec)                                                    
           stop                                                               
        end if                                                                
                                                                              
        if      (bur(irec) .eq. "P") then                                     
           i3 = 1                                                             
        else if (bur(irec) .eq. "H") then                                     
           i3 = 2                                                             
        else if (bur(irec) .eq. "B") then                                     
           i3 = 3                                                             
        else if (bur(irec) .eq. "U") then                                     
           i3 = 4                                                             
        else                                                                  
           write(6,'("Unexpected burial state bur(",i3.3,")=",a1)')           
     *       irec,bur(irec)                                                   
           stop                                                               
        end if                                                                
                                                                              
        srlw(i1a,i2,i3,1) = srlw(i1a,i2,i3,1) + 1.                            
        srlw(i1a,i2,i3,2) = srlw(i1a,i2,i3,2) + arlw(irec)                    
        srlw(i1a,i2,i3,3) = srlw(i1a,i2,i3,3) + arlw(irec)*arlw(irec)         
                                                                              
        srwh(i1b,i2,i3,1) = srwh(i1b,i2,i3,1) + 1.                            
        srwh(i1b,i2,i3,2) = srwh(i1b,i2,i3,2) + arwh(irec)                    
        srwh(i1b,i2,i3,3) = srwh(i1b,i2,i3,3) + arwh(irec)*arwh(irec)         
                                                                              
        srlh(i1c,i2,i3,1) = srlh(i1c,i2,i3,1) + 1.                            
        srlh(i1c,i2,i3,2) = srlh(i1c,i2,i3,2) + arlh(irec)                    
        srlh(i1c,i2,i3,3) = srlh(i1c,i2,i3,3) + arlh(irec)*arlh(irec)         
                                                                              
        sarea(i1a,i2,i3,1) = sarea(i1a,i2,i3,1) + 1.                          
        sarea(i1a,i2,i3,2) = sarea(i1a,i2,i3,2) + area(irec)                  
                                                                              
                                                                              
 40   continue                                                                
                                                                              
      do i1=1,ndbins                                                          
         do i2=1,4                                                            
            do i3=1,4                                                         
                                                                              
               if (srlw(i1,i2,i3,1) .gt. 0.5) then                            
                 av_arlw(i1,i2,i3,2) = srlw(i1,i2,i3,2)/srlw(i1,i2,i3,1)      
               else                                                           
                 av_arlw(i1,i2,i3,2) = 0.                                     
               end if                                                         
                                                                              
               if (srlw(i1,i2,i3,1) .gt. 0.5) then                            
                 av_arlw(i1,i2,i3,3) =                                        
     *                       sqrt(srlw(i1,i2,i3,3)/srlw(i1,i2,i3,1) -         
     *                         av_arlw(i1,i2,i3,2)*av_arlw(i1,i2,i3,2))       
               else                                                           
                 av_arlw(i1,i2,i3,3) = 0.                                     
               end if                                                         
                                                                              
               if (srwh(i1,i2,i3,1) .gt. 0.5) then                            
                 av_arwh(i1,i2,i3,2) = srwh(i1,i2,i3,2)/srwh(i1,i2,i3,1)      
               else                                                           
                 av_arwh(i1,i2,i3,2) = 0.                                     
               end if                                                         
                                                                              
               if (srwh(i1,i2,i3,1) .gt. 0.5) then                            
                  av_arwh(i1,i2,i3,3) =                                       
     *                      sqrt(srwh(i1,i2,i3,3)/srwh(i1,i2,i3,1) -          
     *                        av_arwh(i1,i2,i3,2)**2)                         
               else                                                           
                 av_arwh(i1,i2,i3,3) = 0.                                     
               end if                                                         
                                                                              
               if (srlh(i1,i2,i3,1) .gt. 0.5) then                            
                 av_arlh(i1,i2,i3,2) = srlh(i1,i2,i3,2)/srlh(i1,i2,i3,1)      
               else                                                           
                 av_arlh(i1,i2,i3,2) = 0.                                     
               end if                                                         
                                                                              
               if (srlh(i1,i2,i3,1) .gt. 0.5) then                            
                 av_arlh(i1,i2,i3,3) =                                        
     *                      sqrt(srlh(i1,i2,i3,3)/srlh(i1,i2,i3,1) -          
     *                        av_arlh(i1,i2,i3,2)**2)                         
               else                                                           
                 av_arlh(i1,i2,i3,3) = 0.                                     
               end if                                                         
                                                                              
            end do                                                            
         end do                                                               
      end do                                                                  
                                                                              
                                                                              
C Output average axial rations to unit 16                                     
                                                                              
      do i2=1,4                                                               
         write(16,'("*********************************************",          
     *                "*******************************************")')        
         do i3=1,4                                                            
            write(16,'(/,8x,"SURFACE IS: ",a17,4x,                            
     *                   "BURIAL STATE IS:",a23)')                            
     *         surf_type(i2),bur_state(i3)                                    
         write(16,'(/,"                      AREA       ",                    
     *                "MIN(L,W)/MAX(L,W)  MIN(W,H)/MAX(W,H)  "                
     *                "MIN(L,H)/MAX(L,H)")')                                  
            write(16,'(  "I1    AVG_DIAM                 ",                   
     *                   "   AVG+/-StD    N",                                 
     *                   "     AVG+/-StD    N",                               
     *                      "     AVG+/-StD    N")')                          
                                                                              
            do i=1,3                                                          
               trlw(i2,i3,i) = 0.                                             
               trwh(i2,i3,i) = 0.                                             
               trlh(i2,i3,i) = 0.                                             
               if (i .ne. 3) tarea(i2,i3,i) = 0.                              
            end do                                                            
                                                                              
            do i1=1,ndbins                                                    
               nsrlw  = srlw(i1,i2,i3,1)                                      
               nsrwh  = srwh(i1,i2,i3,1)                                      
               nsrlh  = srlh(i1,i2,i3,1)                                      
               nsarea = sarea(i1,i2,i3,1)                                     
                                                                              
               write(16,'(i2,1x,a14,f8.4,1x,a1,i3.3,a1,                       
     *              3(2x,f4.2,a3,f4.2,1x,a1,i3.3,a1))')                       
     *           i1,avg_diam(i1),                                             
     *           sarea(i1,i2,i3,2),lparen,nsarea,rparen,                      
     *           av_arlw(i1,i2,i3,2),plusminus,av_arlw(i1,i2,i3,3),           
     *           lparen,nsrlw,rparen,                                         
     *           av_arwh(i1,i2,i3,2),plusminus,av_arwh(i1,i2,i3,3),           
     *           lparen,nsrwh,rparen,                                         
     *           av_arlh(i1,i2,i3,2),plusminus,av_arlh(i1,i2,i3,3),           
     *           lparen,nsrlh,rparen                                          
                                                                              
               if (i1 .eq. 1 ) then                                           
                 write(19,'(i1,a1,i1,a1,i2,a1,"     ",a1,"     ",a1,          
     *              f7.4,a1,i3,a1,                                            
     *              3(f4.2,a1,f4.2,a1,i3,a1))')                               
     *             i2,comma,i3,comma,i1,comma,comma,comma,                    
     *             sarea(i1,i2,i3,2),comma,nsarea,comma,                      
     *             av_arlw(i1,i2,i3,2),comma,av_arlw(i1,i2,i3,3),             
     *             comma,nsrlw,comma,                                         
     *             av_arwh(i1,i2,i3,2),comma,av_arwh(i1,i2,i3,3),             
     *             comma,nsrwh,comma,                                         
     *             av_arlh(i1,i2,i3,2),comma,av_arlh(i1,i2,i3,3),             
     *             comma,nsrlh                                                
               else                                                           
                 write(19,'(i1,a1,i1,a1,i2,a1,a5,a1,a5,                       
     *              a1,f7.4,a1,i3,a1,                                         
     *              3(f4.2,a1,f4.2,a1,i3,a1))')                               
     *             i2,comma,i3,comma,i1,comma,avg_diam(i1)(01:05),            
     *             comma,avg_diam(i1)(10:14),                                 
     *             comma,sarea(i1,i2,i3,2),comma,nsarea,comma,                
     *             av_arlw(i1,i2,i3,2),comma,av_arlw(i1,i2,i3,3),             
     *             comma,nsrlw,comma,                                         
     *             av_arwh(i1,i2,i3,2),comma,av_arwh(i1,i2,i3,3),             
     *             comma,nsrwh,comma,                                         
     *             av_arlh(i1,i2,i3,2),comma,av_arlh(i1,i2,i3,3),             
     *             comma,nsrlh                                                
               end if                                                         
                                                                              
               if (i1 .ne. 1) then                                            
                  do i=1,3                                                    
                                                                              
                     if (i.ne.3) tarea(i2,i3,i) = tarea(i2,i3,i) +            
     *                                              sarea(i1,i2,i3,i)         
                     trlw(i2,i3,i) = trlw(i2,i3,i) + srlw(i1,i2,i3,i)         
                     trwh(i2,i3,i) = trwh(i2,i3,i) + srwh(i1,i2,i3,i)         
                     trlh(i2,i3,i) = trlh(i2,i3,i) + srlh(i1,i2,i3,i)         
                  end do                                                      
               end if                                                         
                                                                              
            end do                                                            
                                                                              
            if (trlw(i2,i3,1) .ge. 0.5) then                                  
               trlw(i2,i3,2) = trlw(i2,i3,2)/trlw(i2,i3,1)                    
               trlw(i2,i3,3) = sqrt(trlw(i2,i3,3)/trlw(i2,i3,1) -             
     *                          trlw(i2,i3,2)*trlw(i2,i3,2))                  
            end if                                                            
                                                                              
            if (trwh(i2,i3,1) .ge. 0.5) then                                  
               trwh(i2,i3,2) = trwh(i2,i3,2)/trwh(i2,i3,1)                    
               trwh(i2,i3,3) = sqrt(trwh(i2,i3,3)/trwh(i2,i3,1) -             
     *                          trwh(i2,i3,2)*trwh(i2,i3,2))                  
            end if                                                            
                                                                              
            if (trlh(i2,i3,1) .ge. 0.5) then                                  
               trlh(i2,i3,2) = trlh(i2,i3,2)/trlh(i2,i3,1)                    
               trlh(i2,i3,3) = sqrt(trlh(i2,i3,3)/trlh(i2,i3,1) -             
     *                          trlh(i2,i3,2)*trlh(i2,i3,2))                  
            end if                                                            
                                                                              
            ntrlw = trlw(i2,i3,1)                                             
            ntrwh = trwh(i2,i3,1)                                             
            ntrlh = trlh(i2,i3,1)                                             
            ntarea = tarea(i2,i3,1)                                           
                                                                              
            write(16,'(3x,"    Totals    ",f8.4,1x,a1,i3.3,a1,                
     *              3(2x,f4.2,a3,f4.2,1x,a1,i3.3,a1),/)')                     
     *           tarea(i2,i3,2),lparen,ntarea,rparen,                         
     *           trlw(i2,i3,2),plusminus,trlw(i2,i3,3),                       
     *            lparen,ntrlw,rparen,                                        
     *           trwh(i2,i3,2),plusminus,trwh(i2,i3,3),                       
     *            lparen,ntrwh,rparen,                                        
     *           trlh(i2,i3,2),plusminus,trlh(i2,i3,3),                       
     *            lparen,ntrlh,rparen                                         
                                                                              
            write(19,'(i1,a1,i1,a1,"  ",a1,"0.001",a1,"9.999",a1,             
     *              f7.4,a1,i3,a1,                                            
     *              3(f4.2,a1,f4.2,a1,i3,a1),/)')                             
     *           i2,comma,i3,comma,comma,comma,comma,                         
     *           tarea(i2,i3,2),comma,ntarea,comma,                           
     *           trlw(i2,i3,2),comma,trlw(i2,i3,3),                           
     *            comma,ntrlw,comma,                                          
     *           trwh(i2,i3,2),comma,trwh(i2,i3,3),                           
     *            comma,ntrwh,comma,                                          
     *           trlh(i2,i3,2),comma,trlh(i2,i3,3),                           
     *            comma,ntrlh                                                 
                                                                              
         end do                                                               
      end do                                                                  
                                                                              
C Output counts and sums to unit 18                                           
                                                                              
      do i2=1,4                                                               
         write(18,'("*********************************************",          
     *                "*******************************************")')        
         do i3=1,4                                                            
            write(18,'(/,7x,"SURFACE IS: ",a17,4x,                            
     *                   "BURIAL STATE IS:",a23)')                            
     *         surf_type(i2),bur_state(i3)                                    
            write(18,'(/,"                      AREA              ARLW",      
     *                   "               ARWH  ",                             
     *                      "             ARLH     ")')                       
            write(18,'(  "I1    AVG_DIAM                 ",                   
     *                   "    SUM SUM^2   N",                                 
     *                   "      SUM SUM^2   N",                               
     *                      "      SUM SUM^2   N")')                          
            do i1=1,ndbins                                                    
               nsrlw  = srlw(i1,i2,i3,1)                                      
               nsrwh  = srwh(i1,i2,i3,1)                                      
               nsrlh  = srlh(i1,i2,i3,1)                                      
               nsarea = sarea(i1,i2,i3,1)                                     
                                                                              
               write(18,'(i2,1x,a14,f8.4,1x,a1,i3.3,a1,                       
     *              3(2x,f5.2,1x,f5.2,1x,a1,i3.3,a1))')                       
     *           i1,avg_diam(i1),                                             
     *           sarea(i1,i2,i3,2),lparen,nsarea,rparen,                      
     *           srlw(i1,i2,i3,2),srlw(i1,i2,i3,3),                           
     *           lparen,nsrlw,rparen,                                         
     *           srwh(i1,i2,i3,2),srwh(i1,i2,i3,3),                           
     *           lparen,nsrwh,rparen,                                         
     *           srlh(i1,i2,i3,2),srlh(i1,i2,i3,3),                           
     *           lparen,nsrlh,rparen                                          
                                                                              
            end do                                                            
            write(18,'()')                                                    
         end do                                                               
      end do                                                                  
                                                                              
                                                                              
                                                                              
C Compute the axial ratio statistics summed over the three main surface       
C  types (drift, far field, near field) for perched and partially buried rocks
C        i0 = 1      LW axial ratio                                           
C        i0 = 2      WH axial ratio                                           
C        i0 = 3      LH axial ratio                                           
                                                                              
      do i3=1,2                                                               
                                                                              
         write(16,'(//,2x,"AXIAL RATIO -- AVERAGE OVER DRIFT, ",              
     *                    "FAR, AND NEAR FIELDS")')                           
         write(16,'(49x,"BURIAL STATE IS:",a23)') bur_state(i3)               
         write(16,'(/,"                      AREA       ",                    
     *                "MIN(L,W)/MAX(L,W)  MIN(W,H)/MAX(W,H)  "                
     *                "MIN(L,H)/MAX(L,H)")')                                  
         write(16,'(  "I1    AVG_DIAM                 ",                      
     *                   "   AVG+/-StD    N",                                 
     *                   "     AVG+/-StD    N",                               
     *                      "     AVG+/-StD    N")')                          
                                                                              
         t_ar_sarea(1) = 0.                                                   
         t_ar_sarea(2) = 0.                                                   
                                                                              
         do i0 = 1,3                                                          
            do i4 = 1,3                                                       
               t_ar_surf_av(i0,i4) = 0.                                       
            end do                                                            
         end do                                                               
                                                                              
         do i1=1,ndbins                                                       
            ar_sarea(i1,i3,1) = sarea(i1,1,i3,1) + sarea(i1,2,i3,1) +         
     *                               sarea(i1,3,i3,1)                         
            ar_sarea(i1,i3,2) = sarea(i1,1,i3,2) + sarea(i1,2,i3,2) +         
     *                               sarea(i1,3,i3,2)                         
            t_ar_sarea(1) = t_ar_sarea(1) + ar_sarea(i1,i3,1)                 
            t_ar_sarea(2) = t_ar_sarea(2) + ar_sarea(i1,i3,2)                 
                                                                              
            do i4=1,3                                                         
               do i0 = 1,3                                                    
                  ar_surf_av(i0,i1,i3,i4) = 0.                                
               end do                                                         
               do i2=1,3                                                      
                  ar_surf_av(1,i1,i3,i4) = ar_surf_av(1,i1,i3,i4) +           
     *                                        srlw(i1,i2,i3,i4)               
                  ar_surf_av(2,i1,i3,i4) = ar_surf_av(2,i1,i3,i4) +           
     *                                        srwh(i1,i2,i3,i4)               
                  ar_surf_av(3,i1,i3,i4) = ar_surf_av(3,i1,i3,i4) +           
     *                                        srlh(i1,i2,i3,i4)               
               end do                                                         
               do i0 = 1,3                                                    
                  t_ar_surf_av(i0,i4) = t_ar_surf_av(i0,i4) +                 
     *                                    ar_surf_av(i0,i1,i3,i4)             
               end do                                                         
            end do                                                            
                                                                              
            lw_n  = int(ar_surf_av(1,i1,i3,1) + 0.0001)                       
            lw_av = ar_surf_av(1,i1,i3,2)/ar_surf_av(1,i1,i3,1)               
            lw_sd = sqrt(ar_surf_av(1,i1,i3,3)/ar_surf_av(1,i1,i3,1) -        
     *                lw_av*lw_av)                                            
                                                                              
            wh_n  = int(ar_surf_av(2,i1,i3,1) + 0.0001)                       
            wh_av = ar_surf_av(2,i1,i3,2)/ar_surf_av(2,i1,i3,1)               
            wh_sd = sqrt(ar_surf_av(2,i1,i3,3)/ar_surf_av(2,i1,i3,1) -        
     *                wh_av*wh_av)                                            
                                                                              
            lh_n  = int(ar_surf_av(3,i1,i3,1) + 0.0001)                       
            lh_av = ar_surf_av(3,i1,i3,2)/ar_surf_av(3,i1,i3,1)               
            lh_sd = sqrt(ar_surf_av(3,i1,i3,3)/ar_surf_av(3,i1,i3,1) -        
     *                lh_av*lh_av)                                            
                                                                              
            nsarea = int(ar_sarea(i1,i3,1) + 0.00001)                         
                                                                              
            write(16,'(i2,1x,a14,f8.4,1x,a1,i3.3,a1,                          
     *              3(2x,f4.2,a3,f4.2,1x,a1,i3.3,a1))')                       
     *           i1,avg_diam(i1),                                             
     *           ar_sarea(i1,i3,2),lparen,nsarea,rparen,                      
     *           lw_av,plusminus,lw_sd,                                       
     *           lparen,lw_n,rparen,                                          
     *           wh_av,plusminus,wh_sd,                                       
     *           lparen,wh_n,rparen,                                          
     *           lh_av,plusminus,lh_sd,                                       
     *           lparen,lh_n,rparen                                           
         end do                                                               
                                                                              
         lw_n  = int(t_ar_surf_av(1,1) + 0.0001)                              
         lw_av = t_ar_surf_av(1,2)/t_ar_surf_av(1,1)                          
         lw_sd = sqrt(t_ar_surf_av(1,3)/t_ar_surf_av(1,1) -                   
     *                lw_av*lw_av)                                            
                                                                              
         wh_n  = int(t_ar_surf_av(2,1) + 0.0001)                              
         wh_av = t_ar_surf_av(2,2)/t_ar_surf_av(2,1)                          
         wh_sd = sqrt(t_ar_surf_av(2,3)/t_ar_surf_av(2,1) -                   
     *                wh_av*wh_av)                                            
                                                                              
         lh_n  = int(t_ar_surf_av(3,1) + 0.0001)                              
         lh_av = t_ar_surf_av(3,2)/t_ar_surf_av(3,1)                          
         lh_sd = sqrt(t_ar_surf_av(3,3)/t_ar_surf_av(3,1) -                   
     *                lh_av*lh_av)                                            
                                                                              
         nsarea = int(t_ar_sarea(1) + 0.00001)                                
                                                                              
         write(16,'(3x,"    Totals    ",f8.4,1x,a1,i3.3,a1,                   
     *              3(2x,f4.2,a3,f4.2,1x,a1,i3.3,a1))')                       
     *           t_ar_sarea(2),lparen,nsarea,rparen,                          
     *           lw_av,plusminus,lw_sd,                                       
     *           lparen,lw_n,rparen,                                          
     *           wh_av,plusminus,wh_sd,                                       
     *           lparen,wh_n,rparen,                                          
     *           lh_av,plusminus,lh_sd,                                       
     *           lparen,lh_n,rparen                                           
                                                                              
      end do                                                                  
                                                                              
                                                                              
C Find area as a function of bin over all rock types and burials              
                                                                              
                                                                              
                                                                              
      write(16,'(//,"*********************************************",          
     *             "*******************************************")')           
      write(16,'(  "        --------- DRIFT -------    ------ FA",            
     *             "R FIELD ------    ----- NEAR FIELD ------ ")')            
      write(16,'(  "       ",3(" Surface Area:",f6.2," m^2",3x))')            
     *   (surface_area(i),i=1,3)                                              
      write(16,'(/,"  D>=    NUM   NUM   AREA   AREA    NUM   NU",            
     *             "M   AREA   AREA    NUM   NUM   AREA   AREA ")')           
      write(16,'(  "        (BIN) (CUM) (BIN)  (CUM)   (BIN) (CU",            
     *             "M) (BIN)  (CUM)   (BIN) (CUM) (BIN)  (CUM) ")')           
      write(19,'(" ")')                                                       
                                                                              
      do i=1,3                                                                
         n_cum(i) = 0.                                                        
         area_cum(i) = 0.                                                     
      end do                                                                  
                                                                              
      do i1=ndbins,2,-1                                                       
         do i=1,3                                                             
            n_per_bin(i1,i) = 0.                                              
            area_per_bin(i1,i) = 0.                                           
         end do                                                               
         do i2=1,3                                                            
            do i3 = 1,4                                                       
               n_per_bin(i1,i2) = n_per_bin(i1,i2) +                          
     *                           sarea(i1,i2,i3,1)                            
               area_per_bin(i1,i2) = area_per_bin(i1,i2) +                    
     *                           sarea(i1,i2,i3,2)                            
            end do                                                            
            n_cum(i2) = n_cum(i2) +                                           
     *                            n_per_bin(i1,i2)                            
            area_cum(i2) = area_cum(i2) +                                     
     *                                area_per_bin(i1,i2)                     
         end do                                                               
         write(16,'(1x,a5,3(1x,2f6.0,2f7.4)))') avg_diam(i1)(01:05),          
     *    (n_per_bin(i1,i),n_cum(i),                                          
     *    area_per_bin(i1,i),area_cum(i),i=1,3)                               
         write(19,'(1x,a5,a1,3(f4.0,a1,f4.0,a1,f6.4,a1,f6.4,a1))')            
     *    avg_diam(i1)(01:05),comma,                                          
     *    (n_per_bin(i1,i),comma,n_cum(i),comma,                              
     *    area_per_bin(i1,i),comma,area_cum(i),comma,i=1,2),                  
     *    (n_per_bin(i1,i),comma,n_cum(i),comma,                              
     *    area_per_bin(i1,i),comma,area_cum(i),i=3,3)                         
      end do                                                                  
                                                                              
                                                                              
      ff_interrock_area = surface_area(2)                                     
      nf_interrock_area = surface_area(3)                                     
      do i1 = ndbins,n_d_thresh,-1                                            
         ff_interrock_area = ff_interrock_area - area_per_bin(i1,2)           
         nf_interrock_area = nf_interrock_area - area_per_bin(i1,3)           
      end do                                                                  
                                                                              
      factor = ff_interrock_area/nf_interrock_area + 1.                       
                                                                              
      write(16,'(/,"Far Field Inter-Rock Area:  ",f10.4,/,                    
     *             "Near Field Inter-Rock Area: ",f10.4,/,                    
     *             "FF/NF Scale Factor:         ",f10.4)')                    
     *             ff_interrock_area,nf_interrock_area,factor                 
                                                                              
      write(16,'(//,"*********************************************",          
     *             "*******************************************")')           
      write(16,'(  "        --- NUMBER DIST ---    --- AREA DIST ---",        
     *             "  -- NORM NUM DIST --  --NORM AREA DIST-")')              
      write(16,'(  "  D>=      (BIN)    (CUM)        (BIN)   (CUM)  ",        
     *             "   (BIN)    (CUM)        (BIN)   (CUM)  ")')              
      write(19,'(" ")')                                                       
                                                                              
      comp_cum(1) = 0.                                                        
      comp_cum(2) = 0.                                                        
                                                                              
      SA = surface_area(1) + surface_area(2) + surface_area(3)                
                                                                              
      do i1=ndbins,2,-1                                                       
                                                                              
                                                                              
         if (i1 .ge. n_d_thresh) then                                         
            comp(i1,1) = n_per_bin(i1,1) + n_per_bin(i1,2) +                  
     *                      n_per_bin(i1,3)                                   
            comp(i1,2) = area_per_bin(i1,1) + area_per_bin(i1,2) +            
     *                      area_per_bin(i1,3)                                
         else                                                                 
            comp(i1,1) = n_per_bin(i1,1) + n_per_bin(i1,3)*factor             
            comp(i1,2) = area_per_bin(i1,1) + area_per_bin(i1,3)*factor       
         end if                                                               
                                                                              
         comp(i1,3) = comp(i1,1)/SA                                           
         comp(i1,4) = comp(i1,2)/SA                                           
                                                                              
         comp_cum(1) = comp_cum(1) + comp(i1,1)                               
         comp_cum(2) = comp_cum(2) + comp(i1,2)                               
         comp_cum(3) = comp_cum(3) + comp(i1,3)                               
         comp_cum(4) = comp_cum(4) + comp(i1,4)                               
                                                                              
         write(16,'(1x,a5,2(1x,2f9.2,5x,2f8.4)))') avg_diam(i1)(01:05),       
     *    comp(i1,1),comp_cum(1),comp(i1,2),comp_cum(2),                      
     *    comp(i1,3),comp_cum(3),comp(i1,4),comp_cum(4)                       
                                                                              
         write(19,'(1x,a5,a1,2(f7.2,a1,f7.2,a1,f7.4,a1,f7.4,a1)))')           
     *    avg_diam(i1)(01:05),comma,                                          
     *    comp(i1,1),comma, comp_cum(1),comma,                                
     *    comp(i1,2),comma, comp_cum(2), comma,                               
     *    comp(i1,3),comma,comp_cum(3),comma,                                 
     *    comp(i1,4),comma,comp_cum(4)                                        
                                                                              
      end do                                                                  
                                                                              
                                                                              
      stop                                                                    
      end                                                                     
                                                                              
