!===============================================================================
! One of Andy Nowacki's Fortran utility modules for dealing with seismic
! anisotropy and other problems.
!
! Andy Nowacki <andy.nowacki@bristol.ac.uk>
!
! See the file LICENCE for licence details.
!===============================================================================
   module density_1d
!------------------------------------------------------------------------------
!
! Returns density for some global Earth models, given a depth.
!
!
!  ** size constants
      integer, parameter, private :: i4 = selected_int_kind(9) ; ! long int
      integer, parameter, private :: r4 = selected_real_kind(6,37) ; ! SP
      integer, parameter, private :: r8 = selected_real_kind(15,307) ; ! DP
      
!  ** precision selector
      integer, parameter, private :: rs = r8
      
!  ** maths constants and other useful things
      real(rs), parameter, private :: pi = 3.141592653589793238462643 ;
      real(rs), parameter, private :: to_rad = 1.74532925199433e-002 ;  
      real(rs), parameter, private :: to_deg = 57.2957795130823e0 ;  
      real(rs), parameter, private :: to_km = 111.194926644559 ;      

      real(rs), parameter, private :: big_number = 10.e36 ;      
      
   contains
   
!===============================================================================
   function ak135den(depth)
!===============================================================================
!  Uses AK135 model

   implicit none
   
   real(rs),intent(in)  :: depth
   real(rs)  :: dep(136),den(136),ak135den,den1,den2,dep1,dep2
   integer   :: i
   logical   :: twodepths
   data dep / &
		 0.0,   20.0,   20.0,   35.0,   35.0, &
		77.5,  120.0,  165.0,  210.0,  210.0, &
	   260.0,  310.0,  360.0,  410.0,  410.0, &
	   460.0,  510.0,  560.0,  610.0,  660.0, &
	   660.0,  710.0,  760.0,  809.5,  859.0, &
	   908.5,  958.0, 1007.5, 1057.0, 1106.5, &
	  1156.0, 1205.5, 1255.0, 1304.5, 1354.0, &
	  1403.5, 1453.0, 1502.5, 1552.0, 1601.5, &
	  1651.0, 1700.5, 1750.0, 1799.5, 1849.0, &
	  1898.5, 1948.0, 1997.5, 2047.0, 2096.5, &
	  2146.0, 2195.5, 2245.0, 2294.5, 2344.0, &
	  2393.5, 2443.0, 2492.5, 2542.0, 2591.5, &
	  2640.0, 2690.0, 2740.0, 2740.0, 2789.7, &
	  2839.3, 2891.5, 2891.5, 2939.3, 2989.7, &
	  3040.0, 3090.3, 3140.7, 3191.0, 3241.3, &
	  3291.7, 3342.0, 3392.3, 3442.6, 3493.0, &
	  3543.3, 3593.6, 3644.0, 3694.3, 3744.6, &
	  3795.0, 3845.3, 3895.6, 3945.9, 3996.3, &
	  4046.6, 4096.9, 4147.3, 4197.6, 4247.9, &
	  4298.3, 4348.6, 4398.9, 4449.3, 4499.6, &
	  4549.9, 4600.3, 4650.6, 4700.9, 4801.6, &
	  4851.9, 4902.2, 4952.6, 5002.9, 5053.2, &
	  5103.6, 5153.5, 5153.5, 5204.6, 5255.3, &
	  5306.0, 5356.8, 5407.5, 5458.2, 5508.9, &
	  5559.6, 5610.3, 5661.0, 5711.7, 5813.2, &
	  5863.9, 5914.6, 5965.3, 6016.0, 6066.7, &
	  6117.4, 6168.1, 6218.9, 6269.6, 6320.3, &
	  6371.0 /
	  
   data den / &
	   2.44900,  2.44900,  2.71420,  2.71420,  3.29760, &
	   3.29940,  3.30130,  3.34870,  3.39600,  3.39600, &
	   3.46520,  3.53430,  3.60340,  3.67260,  3.79760, &
	   3.86120,  3.92480,  3.98850,  4.05210,  4.11580, &
	   4.33930,  4.38960,  4.43990,  4.47010,  4.50290, &
	   4.53500,  4.56640,  4.59700,  4.62700,  4.65630, &
	   4.68490,  4.71300,  4.74040,  4.76730,  4.79360, &
	   4.81950,  4.84480,  4.87000,  4.89420,  4.91820, &
	   4.94230,  4.96530,  4.98820,  5.01090,  5.03330, &
	   5.05530,  5.07720,  5.09900,  5.12060,  5.14220, &
	   5.16380,  5.18480,  5.20610,  5.22700,  5.24810, &
	   5.26980,  5.29070,  5.31220,  5.33380,  5.35600, &
	   5.37760,  5.39990,  5.42240,  5.42240,  5.42380, &
	   5.42510,  5.42650,  9.91450,  9.99420, 10.07220, &
	  10.14850, 10.22330, 10.29640, 10.36790, 10.43780, &
	  10.50620, 10.57310, 10.63850, 10.70230, 10.76470, &
	  10.82570, 10.88520, 10.94340, 11.00010, 11.05550, &
	  11.10950, 11.16230, 11.21370, 11.26390, 11.31270, &
	  11.36040, 11.40690, 11.45210, 11.49620, 11.53910, &
	  11.58090, 11.62160, 11.66120, 11.69980, 11.73730, &
	  11.77370, 11.80920, 11.84370, 11.87720, 11.94140, &
	  11.97220, 12.00010, 12.03110, 12.05930, 12.08670, &
	  12.11330, 12.13910, 12.70370, 12.72890, 12.75300, &
	  12.77600, 12.79800, 12.81880, 12.83870, 12.85740, &
	  12.87510, 12.89170, 12.90720, 12.92170, 12.94740, &
	  12.95860, 12.96880, 12.97790, 12.98590, 12.99290, &
	  12.99880, 13.00360, 13.00740, 13.01000, 13.01170, &
	  13.01220 /

!  Ensure depth is in km and not too large
   if (depth > 6371.0) then
      write(*,'(a,f9.1,a)') '  ak135den: depth ',depth,' is too large.'
      stop
   endif

!  Look for the specific point first
   do i=1,136
      if (depth == dep(i)) then 
         ak135den = den(i)
         twodepths = .false.
         exit
      else
!  Otherwise, find the points above and below this depth
         if (depth > dep(i) .and. depth < dep(i+1)) then
            dep1 = dep(i)
            den1 = den(i)
            dep2 = dep(i+1)
            den2 = den(i+1)
            twodepths = .true.
            exit
         endif
      endif
   enddo
   
   if (twodepths) then
!  Extrapolate linearly to find density
      ak135den = den1 + (depth-dep1)*(den2-den1)/(dep2-dep1)
   endif
   
   return
   
   end function ak135den
!------------------------------------------------------------------------------

!===============================================================================
   function premden(depth)
!===============================================================================
!  PREM
   
   implicit none
   
   real(rs),intent(in) :: depth
   real(rs)  :: dep(88),den(88),premden,den1,den2,dep1,dep2
   integer   :: i
   logical   :: twodepths
   data dep /        0.00,   15.00,   15.00,   24.40,   24.40, &
		   40.00,   60.00,   80.00,  115.00,  150.00,  185.00, &
		  220.00,  220.00,  265.00,  310.00,  355.00,  400.00, &
		  400.00,  450.00,  500.00,  550.00,  600.00,  635.00, &
		  670.00,  670.00,  721.00,  771.00,  871.00,  971.00, &
		 1071.00, 1171.00, 1271.00, 1371.00, 1471.00, 1571.00, &
		 1671.00, 1771.00, 1871.00, 1971.00, 2071.00, 2171.00, &
		 2271.00, 2371.00, 2471.00, 2571.00, 2671.00, 2741.00, &
		 2771.00, 2871.00, 2891.00, 2891.00, 2971.00, 3071.00, &
		 3171.00, 3271.00, 3371.00, 3471.00, 3571.00, 3671.00, &
		 3771.00, 3871.00, 3971.00, 4017.00, 4171.00, 4271.00, &
		 4371.00, 4471.00, 4571.00, 4671.00, 4771.00, 4871.00, &
		 4971.00, 5071.00, 5149.50, 5149.50, 5171.00, 5271.00, &
		 5371.00, 5471.00, 5571.00, 5671.00, 5771.00, 5871.00, &
		 5971.00, 6071.00, 6171.00, 6271.00, 6371.00/
   
   data den /       2.60000,  2.60000,  2.90000,  2.90000,  3.38076, &
		  3.37906,  3.37688,  3.37471,  3.37091,  3.36710,  3.36330, &
		  3.35950,  3.43578,  3.46264,  3.48951,  3.51639,  3.54325, &
		  3.72378,  3.78678,  3.84980,  3.91282,  3.97584,  3.98399, &
		  3.99214,  4.38071,  4.41241,  4.44317,  4.50372,  4.56307, &
		  4.62129,  4.67844,  4.73460,  4.78983,  4.84422,  4.89783, &
		  4.95073,  5.00299,  5.05469,  5.10590,  5.15669,  5.20713, &
		  5.25729,  5.30724,  5.35706,  5.40681,  5.45657,  5.49145, &
		  5.50642,  5.55641,  5.56645,  9.90349, 10.02940, 10.18134, &
		 10.32726, 10.46727, 10.60152, 10.73012, 10.85321, 10.97091, &
		 11.08335, 11.19067, 11.29298, 11.39042, 11.48311, 11.57119, &
		 11.65478, 11.73401, 11.80900, 11.87990, 11.94682, 12.00989, &
		 12.06924, 12.12500, 12.16634, 12.76360, 12.77493, 12.82501, &
		 12.87073, 12.91211, 12.94912, 12.98178, 13.01009, 13.03404, &
		 13.05364, 13.06888, 13.07977, 13.08630, 13.08848 /

!  Ensure depth is in km and not too large
   if (depth > 6371.0) then
      write(*,'(a,f9.1,a)') '  premden: depth ',depth,' is too large.'
      stop
   endif

!  Look for the specific point first
   do i=1,88
      if (depth == dep(i)) then 
         premden = den(i)
         twodepths = .false.
         exit
      else
!  Otherwise, find the points above and below this depth
         if (depth > dep(i) .and. depth < dep(i+1)) then
            dep1 = dep(i)
            den1 = den(i)
            dep2 = dep(i+1)
            den2 = den(i+1)
            twodepths = .true.
            exit
         endif
      endif
   enddo
   
   if (twodepths) then
!  Extrapolate linearly to find density
      premden = den1 + (depth-dep1)*(den2-den1)/(dep2-dep1)
   endif
   
   return
   
   end function premden
!------------------------------------------------------------------------------
   
   end module density_1d
