Mapping code

Program mapmaker

use sep
  use random_mod
  use globecalc

implicit none
  real                              :: backaz,foraz,distance,ilat_tmp,ilon_tmp
  real                              :: ptest,rtest,ntest,stest, pi
  real,dimension (:)  , pointer     :: lat,lon,shaz,temp,regime,plate,subregion,datalat,datalon,datashaz
  integer,dimension (:)  , pointer  :: subr
  real,dimension (:,:), pointer     :: azdiff
  integer                           :: ilat,ilon,idat,pos1,pos2,k,maxk
  integer                           :: n1,m1,m2,o1,o2,d1,d2,stat,t,r,s

call sep_init ()

call from_history("n1",n1)

allocate ( lat(n1),lon(n1),shaz(n1),regime(n1),plate(n1),subregion(n1) )
  allocate ( subr(n1) )

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

call from_param( "o1", o1)
  call from_param( "o2", o2)
  call from_param( "d1", d1)
  call from_param( "d2", d2)
  m1= abs((2*o1)/d1)
  m2= abs((2*o2)/d2)

call from_param( "platetest",       ptest,0.)
  call from_param( "regimetest",      rtest,0.)
  call from_param( "subregiontest",   stest,0.)
  call from_param( "fixnorms",        ntest,1.)
  !--all data are used when tests /=0

if( exist_file("randaz") ) then
     call to_history( "n1", n1, "randaz")
     call to_history( "o1",  1, "randaz")
     call to_history( "d1",  1, "randaz")
  end if
  call to_history("n1",m1+1)
  call to_history("n2",m2+1)
  call to_history("o1",o1)
  call to_history("o2",o2)
  call to_history("d1",d1)
  call to_history("d2",d2)
  call to_history("label1","latitude (deg)")
  call to_history("label2","longitude(deg)")
  call to_history("label3","lack of intersections")

call sep_close()
  pi=3.141593
  !---Convert degrees lat/long to radians---------------
  lat= lat * (pi/180) 
  lon= lon * (pi/180)
  shaz=shaz* (pi/180)
  !-NF regime should consider Sh=SH-90, not SH--------------   
  if (ntest/=0) where ( int(regime)==1 ) shaz=abs(asin(sin( shaz(idat)-(pi/2.) )))             
  !--- change ocean subcode to -1 from 0
  where ( (subregion)==0. ) subregion=-2.
  subr=int(subregion)
  write(0,*) 'region test?',rtest,'plate test?', ptest,'subregion test?', stest,'fix normals?', ntest
  !---Test with random azimuths at each data location---------
  if( exist_file("randaz")) then
     allocate(temp(n1))
     call random_number(temp)
     datashaz=temp*pi
  end if

if ((rtest==1.3).or.(rtest==2.3)) then
     regime=int(regime)
     rtest=int(rtest)
  end if
  !------count how many data points going to use
  k=0
  do idat=1,n1
    if (shaz(idat) <=2*pi) then !---null values in data file =999
    if ( (rtest==0).or.(regime(idat)==rtest) ) then!-only consider one stress regime
    if ( (ptest==0).or.(ptest==plate(idat)) ) then !-only consider one plate
    if ( (stest==0).or.( (stest==4).and.(subregion(idat)>0) ).or.(stest==mod(subr(idat),3)+1) ) then
       k=k+1
    end if   
    end if    
    end if  
    end if
  end do
  maxk=k 
  allocate(datalat(maxk))
  allocate(datalon(maxk))
  allocate(datashaz(maxk))
  write(0,*) 'You are about to consider ', maxk,' points of 10921.'
  k=0
  do idat=1,n1
    if (shaz(idat) <=2*pi) then !---null values in data file =999
    if ( (rtest==0).or.(regime(idat)==rtest) ) then!-only consider one stress regime
    if ( (ptest==0).or.(ptest==plate(idat)) ) then !-only consider one plate
    if ( (stest==0).or.( (stest==4).and.(subregion(idat)>0) ).or.(stest==mod(subr(idat),3)+1) ) then
       k=k+1
       datalat(k)=lat(idat)
       datalon(k)=lon(idat)
       datashaz(k)=shaz(idat)
    end if   
    end if    
    end if  
    end if
  end do
  maxk=k
  !---Test with random azimuths at each data location---------
  if( exist_file("randaz")) then
     allocate(temp(n1))
     call random_number(temp)
     datashaz=temp*pi
     call sep_write(datashaz*(180/pi),"randaz")
  end if
  !---prepare a file w/ data locations used when discarding with tests
  if (exist_file("datafile") ) then  
     call sep_init()
     call to_history("n1",maxk,"datafile")
     call to_history("n2",   3,"datafile")
     call sep_write (datalon*(180/pi),  "datafile")
     call sep_write (datalat*(180/pi),  "datafile")
     call sep_write (datashaz*(180/pi), "datafile")
     call sep_close()
  end if
! goto 100   ! for only making actually used data file

write(0,*) 'Now mapping. Please take your seat.'
  allocate ( azdiff(m1+1,m2+1) )
  azdiff=0  !--safety zeroing
  do ilat=o1,o1+m1,d1
  write(0,*) 'done ',(ilat-o1)/d1,' of ',m1
  do ilon=o2,o2+m2,d2
     ilat_tmp = ilat * (pi/180)!----temporary radians of model grid points
     ilon_tmp = ilon * (pi/180)
     do k=1,maxk
                    !---finally to the meat of the program----------------
                    call calc( datalat(k),datalon(k),ilat_tmp,ilon_tmp,distance,backaz,foraz)
                    pos1=ilat-o1    ;       pos2=ilon-o2
                    azdiff(pos1,pos2) = azdiff(pos1,pos2) + asin(sin( (datashaz(k)-backaz)))**2
       end do  ! idata
  end do  ! ilon
  end do  ! idatalat

call sep_write(azdiff)

100 k=2

end program mapmaker