Variogram code

Program variogram

use sep
  use bin1
  use globecalc

implicit none
  real                              :: az1,az2,ov,dv,a,b,c,d,pi
  real,dimension (:)  , pointer     :: lat,lon,shaz,temp,anglesum,ang,delta,ones,fold
  integer                           :: idat,idat2,nv,n1,k,stat

call sep_init ()

call from_history("n1",n1)

allocate (lat(n1),lon(n1),shaz(n1))

call sep_read (lat)
  call sep_read (lon)
  call sep_read (shaz)
  !write(0,*) 'lat is ', lat
  !write(0,*) 'lon is ', lon
  !write(0,*) 'az is  ', shaz

call from_param( "nv", nv)
  call from_param( "ov", ov)
  call from_param( "dv", dv)

call to_history("n1",nv)
  call to_history("n2", 1)
  call to_history("o1",ov)
  call to_history("d1",dv)

call to_history("label2","orthogonality")
  call to_history("label1","kilometers")

call sep_close()
  pi=3.141593
!------------quadrant tests for the subtraction ------
!  k=1
!  if (k/=0) then
!     a=asin(sin( 0.17453295-pi )) !10 in radians
!     b=asin(sin( 3.3161259 -pi )) !190 in radians
!     a=asin(sin( 10*(pi/180) ))                 ! = 0.17453295 
!     b=asin(sin(-10*(pi/180) ))                 ! =-0.17453295
!     c=asin(sin(170*(pi/180) ))                 ! = 0.174532562
!     d=asin(sin(190*(pi/180) ))                 ! =-0.174533218
!     write(0,*) '10 -90=',a*(180/pi),' 190 -90=',b*(180/pi),'   170=',c,'   190=',d
!
!  else   !---to use, must also take out end if at end of program--
!-----------------------------------

lat= lat * (pi/180)             !---Convert degrees lat/long to radians --
     lon= lon * (pi/180)
     shaz=shaz* (pi/180)

write(0,*) 'Doing the variance bit'
     allocate(delta ( n1*(n1-1)/2) )
     allocate(ang   ( n1*(n1-1)/2) )
     delta=0.  ; ang=0.
     !     write(0,*) 'allocated'
     !     write(0,*) 'size of ang is ',size(ang),'delta is ',size(delta)
     k=0
     do idat =1,n1-1 
        !        write(0,*) 'Finished ',k,' of',size(ang),' loops.'
        if( shaz(idat) <=180 ) then                 !null values in data file =999
           do idat2=idat+1,n1
              k=k+1
              if( shaz(idat2) <=180 ) then
                 call calc( lat(idat),lon(idat),lat(idat2),lon(idat2),delta( (idat-1)*n1+(idat2-1) ),az1,az2)
                 ang(k) = sqrt( abs(asin(sin( shaz(idat)-shaz(idat2) ))))
              end if
           end do
        end if
     end do
     !---------------------------------------------
     !---------------putting data into distance bins ---------------
     delta= delta * 6371                             !  convert to km's: rad. of earth=6371 km's
     allocate( anglesum(nv) )
     anglesum=0                                      !  initialize to zero
     !     write(0,*)'size of anglesum is ', size(anglesum)
     call bin1_init( nv, ov, dv, delta )             !  bins should go to 0.5*circ. of earth
     stat = bin1_lop( .true., .true.,anglesum, ang)  !  anglesum is output
     !--------------to normalize anlgesum,  fold=Btrans (ones) -----------
     allocate( ones(n1) ) 
     allocate( fold(nv) )
     ones=1.; fold=0.
     call bin1_init( nv, ov, dv, delta )
     stat = bin1_lop( .true., .false., fold, ones)
     !  write(0,*) 'anglesum at 1 is ', anglesum
     where(fold/=0)                                 !  mind not dividing by zero
        anglesum=anglesum/fold                      !  divide by number of data point in the bin
     elsewhere
        anglesum=0
     end where
     !  write(0,*) 'fold at 1 is ', fold
     !  write(0,*) 'anglesum is now', anglesum

call sep_write( anglesum)
!  end if   !---for when you don't really want to run the program
end program variogram