module model
 
 implicit none

 contains

  subroutine model_init2(nbc)
!--------------------- 
 use lnb_param

  double precision :: ro,g,drad,la,mu,ro1,ro2,x
  logical :: solid,test
  integer :: nbc,i,n,j



  !- Calculus parameters    
  pi=dacos(-1.d0)
  epsdb=epsilon(1.d0)
  
  !write(101,*), 'Epsilon machine = ',epsdb
  !write(101,*), 'Constants:'
  !write(101,*), ' pi = ',pi
  !write(101,*), ' G  = ',GG 
  !write(101,*), 'Dimensioned scales:'
  !write(101,*), ' rigidity mu0 = ',mu0
  !write(101,*), ' gravity  go0 = ',go0
  !write(101,*), ' distance r0  = ',r0
  
 

 
 !call earth_homogeneous_layers_prembased(radius,nlayer,roc,lac,muc,g,solid) !yields the Earth layers radii, density, lamé parameters, gravity and solid/liquid state boolean identifier
 
 !if (benchmark_spada) then
 ! muc(1) = 1.d6; muc(2) = 0.d0; muc(3) = 2.28340d11 
 ! muc(4) = 1.05490d11; muc(5) = 0.70363d11; muc(6) = 0.50605d11
 ! roc(1) = 10751.d0; roc(2) = 10750.d0; roc(3) = 4978.d0 
 ! roc(6)=3300.d0
 !else
  !muc(1) = 1.d6; muc(2) = 0.d0; muc(3) = 2.20425d11 
  !muc(4) = 0.870257d11; muc(5) = 0.870257d11; muc(6) = 0.510268d11
  !roc(1) = 10987.9d0; roc(2) = 10986.9d0; roc(3) = 4903.58d0 
  !roc(4) = 3628.29d0; roc(5) = 3628.29d0; roc(6) = 3113.94d0
 !endif

 !vic2=vic/vicratio
 !muc2=muc/mucratio

 if (display) then
	print*, '  Layer top boundary        rho                       la                        mu                        eta' 
 do i = 1,nlayer
  ! print*, radius(i+1),roc(i),lac(i),muc(i),vic(i)
  if (burgers(i)) then
  print*, radius(i+1),roc(i),lac(i),muc(i),vic(i),muc2(i),vic2(i), 'B'
  else
  print*, radius(i+1),roc(i),lac(i),muc(i),vic(i)
  endif
 end do

 endif
 !write(101,*),


 call earth_nlayers(1.d0,ro,la,mu,g,solid)
   go_surf = g
   go0 = g !!! LC+SA: 07/20/2020. 
   ro_mean = 3.d0/4.d0*go_surf/(pi*GG*ra)
	if (display) then
		print*, 'Surface Gravity: ', go_surf, ' m.s^-2, mean density: ', ro_mean, ' kg.m^-3'
		print*, 'Center of the Earth approximated at r=', radius(1)/1e3, ' km'
	endif 
  
  ! Test number of interfaces and of boundary conditions
  n = 0
  !!! >>> SA: 01/27/2018 => following changes made to retrieve love numbers at depth 
  !do i=1,nlayer+1 
 !if (radbc(i)) n=n+1
  !enddo
  n = nlayer+1 
  !!! >>> SA 
  if (nbc/=(6*n)) then
 !write(101,*),
 print*, 'Error nbc should be ',6*n,' instead of ',nbc
 stop
  endif
  
  ! Test solid-liquid boundaries
 test = .false.
  do i = 1,nlayer
 if (radbc(i)) test = .true.
 if (test) then
 if ((.not.soliddim(i)).and.((.not.radbc(i)).or.(.not.radbc(i+1)))) then
  print*, 'Error in model definition: one fluid layer is not set with correct boundary conditions'
  stop
 endif
 endif
  enddo
  
  !-- Display a few quantities
  !write(101,*), ' surface radius = ',ra
  !write(101,*), ' CMB radius = ',rb
  !write(101,*), ' ICB radius = ',rc
  !write(101,*), ' go surface = ',go_surf
  !write(101,*), ' mean density = ',ro_mean 
  
  !-- Diplsay ICB conditions
  !x=rc/ra
  !write(101,*), 
  !call earth_nlayers_viscoelas(x-epsdb,ro1,complex(la,0.d0),complex(mu,0d0),g,solid)
  !write(101,*), 'ICB, g = ',g
  !write(101,*), 'ICB, ro(Inner Core) = ', ro1
  !call earth_nlayers_viscoelas(x,ro2,complex(la,0.d0),complex(mu,0d0),g,solid)
  !write(101,*), 'ICB, g = ',g
  !write(101,*), 'ICB, ro(Outer Core) = ', ro2
  !go_icb = g
  
  !-- Display CMB conditions
  !x=rb/ra
  !write(101,*),
  !call earth_nlayers_viscoelas(x-epsdb,ro1,complex(la,0.d0),complex(mu,0d0),g,solid)
  !write(101,*), 'CMB, g = ',g
  !write(101,*), 'CMB, ro(Outer Core) = ', ro1
  !call earth_nlayers_viscoelas(x,ro2,complex(la,0.d0),complex(mu,0d0),g,solid)
  !write(101,*), 'CMB, g = ',g
  !write(101,*), 'CMB, ro(Mantle) = ', ro2
  !go_cmb = g
 
 end subroutine model_init2


 subroutine earth_homogeneous_layers_prembased(radin,nl,rot,lat,mut,GR,solid)
 use lnb_param
   
   double precision :: x0,GR
   logical :: solid
   integer                 :: i,j,nl
   doubleprecision      :: la,mu,ro,vp,vs,RT
   doubleprecision      :: t1,t2,t3,t4,r1,r2
   doubleprecision, dimension(14) :: r
   doubleprecision, dimension(13,4) :: d,p,s
   doubleprecision, dimension(:), pointer :: radin
   doubleprecision, dimension(:), pointer :: rot,lat,mut
    doubleprecision, dimension(nl+1) :: rad
   
   
   RT = ra
   r(1) = 0.d0;     r(2) = 1221.5d0; r(3) = 3480.d0; r(4) = 3630.d0
   r(5) = 5600.d0;  r(6) = 5701.d0;  r(7) = 5771.d0; r(8) = 5971.d0
   r(9) = 6151.d0; r(10) = 6291.d0; r(11) = 6346.6d0
   r(12) = 6356.d0; r(13) = 6368.d0; r(14) = ra
  
!   if ((r(14)*1.d3) /= ra) stop'Problem in prem: Earth radius /= ra !' 
!
  d(:,:) = 0.d0
  d(1,1) = 13.0885d0;               d(1,3) = -8.8381d0 
  d(2,1) = 12.5815d0; d(2,2) = -1.2638d0; d(2,3) = -3.6426d0; d(2,4) = -5.5281d0
  d(3,1) = 7.9565d0 ; d(3,2) = -6.4761;   d(3,3) = 5.5283d0;  d(3,4) = -3.0807d0
  d(4,1) = 7.9565d0 ; d(4,2) = -6.4761;   d(4,3) = 5.5283d0;  d(4,4) = -3.0807d0
  d(5,1) = 7.9565d0 ; d(5,2) = -6.4761;   d(5,3) = 5.5283d0;  d(5,4) = -3.0807d0
  d(6,1) = 5.3197d0 ; d(6,2) = -1.4836d0
  d(7,1) = 11.2494d0; d(7,2) = -8.0298d0
  d(8,1) = 7.1089d0 ; d(8,2) = -3.8045d0
  d(9,1) = 2.6910d0 ; d(9,2) = 0.6924d0
  d(10,1) = 2.6910d0; d(10,2) = 0.6924d0
  d(11,1) = 2.9d0  
  d(12,1) = 2.6d0  

! ocean
  if (soliddim(13).eqv..false.) then
   d(13,1) = 1.02d0 

! continental
  else
  d(13,1) = d(12,1)
  end if

  p(:,:) = 0.d0
  p(1,1) = 11.2622d0 ; p(1,3) = -6.3640d0
  p(2,1) = 11.0487d0 ; p(2,2) = -4.0362d0; p(2,3)  = 4.8023d0; p(2,4) = -13.5732d0
  p(3,1) = 15.3891d0 ; p(3,2) = -5.3181d0; p(3,3)  = 5.5242d0; p(3,4) = -2.5514d0
  p(4,1) = 24.952d0 ; p(4,2)  = -40.4673d0; p(4,3) = 51.4832d0; p(4,4) = -26.6419d0
  p(5,1) = 29.2766d0 ; p(5,2) = -23.6027d0; p(5,3) = 5.5242d0; p(5,4) = -2.5514d0
  p(6,1) = 19.0957d0 ; p(6,2)  = -9.8672d0
  p(7,1) = 39.7027d0 ; p(7,2)  = -32.6166d0
  p(8,1) = 20.3926d0 ; p(8,2)  = -12.2569d0
  p(9,1) = 4.1875d0 ; p(9,2)  = 3.9382d0
  p(10,1) = 4.1875d0 ; p(10,2) = 3.9382d0
  p(11,1) = 6.8d0 
  p(12,1) = 5.8d0
!
! ocean
  if (soliddim(13).eqv..false.) then
  p(13,1) = 1.45d0 
!
! continental
  else
  p(13,1) = p(12,1)
  end if
!----
!
  s(:,:) = 0.d0
!
  s(1,1) = 3.6678d0; s(1,3) = -4.4475d0

  s(3,1) = 6.9254d0; s(3,2) = 1.4672d0; s(3,3) = -2.0834d0; s(3,4) = 0.9783d0
  s(4,1) = 11.1671d0; s(4,2) = -13.7818d0; s(4,3) = 17.4575d0; s(4,4) = -9.2777d0
  s(5,1) = 22.3459d0; s(5,2) = -17.2473d0; s(5,3) = -2.0834d0; s(5,4) = 0.9783d0
  s(6,1) = 9.9839d0; s(6,2) = -4.9324
  s(7,1) = 22.3512d0; s(7,2) = -18.5856d0 
  s(8,1) = 8.9496d0; s(8,2) = -4.4597
  s(9,1) = 2.1519d0; s(9,2) = 2.3481d0
  s(10,1) = 2.1519d0; s(10,2) = 2.3481d0
  s(11,1) = 3.9d0 
  s(12,1) = 3.2d0 
!
! ocean (please don't modify)
  if (soliddim(13).eqv..false.) then
!
! continental
  else
  s(13,1) = s(12,1)
  end if
!
!
  r(:) = r(:)*1.d3
  
  !- handling the first layer : central sphere
  rad = radin
  rad(1) = 0.d0
  
  do j = 1,nl
  
 ro = 0.d0
 vp = 0.d0
 vs = 0.d0

 do i = 1,13

  r1 = 0.d0
  r2 = 0.d0
  if ((rad(j) > r(i)).and.(rad(j) <= r(i+1))) then
   if (rad(j+1) <= r(i+1)) then
    r2 = rad(j+1)
    r1 = rad(j)
   else
    r2 = r(i+1)
    r1 = rad(j)
   end if
  else if (rad(j) <= r(i)) then
   if ((rad(j+1) > r(i)).and.(rad(j+1) <= r(i+1))) then
    r2 = rad(j+1)
    r1 = r(i)
   else if (rad(j+1) > r(i+1)) then
    r2 = r(i+1)
    r1 = r(i)
   end if
  end if

  t1 = d(i,1)/3.d0
  t2 = d(i,2)/(Rt*4.d0)
  t3 = d(i,3)/((Rt**2)*5.d0)
  t4 = d(i,4)/((Rt**3)*6.d0)
  ro =  ro + t1*(r2**3) + t2*(r2**4) + t3*(r2**5) + t4*(r2**6) - &
   ( t1*(r1**3) + t2*(r1**4) + t3*(r1**5) + t4*(r1**6) )
     
  t1 = p(i,1)/3.d0
  t2 = p(i,2)/(Rt*4.d0)
  t3 = p(i,3)/((Rt**2)*5.d0)
  t4 = p(i,4)/((Rt**3)*6.d0)
  vp =  vp + t1*(r2**3) + t2*(r2**4) + t3*(r2**5) + t4*(r2**6) - &
   ( t1*(r1**3) + t2*(r1**4) + t3*(r1**5) + t4*(r1**6) )
     
  t1 = s(i,1)/3.d0
  t2 = s(i,2)/(Rt*4.d0)
  t3 = s(i,3)/((Rt**2)*5.d0)
  t4 = s(i,4)/((Rt**3)*6.d0)
  vs =  vs + t1*(r2**3) + t2*(r2**4) + t3*(r2**5) + t4*(r2**6) - &
   ( t1*(r1**3) + t2*(r1**4) + t3*(r1**5) + t4*(r1**6) )

 end do
 ro = ro*3 / (rad(j+1)**3-rad(j)**3)
 vp = vp*3 /(rad(j+1)**3-rad(j)**3)
 vs = vs*3 / (rad(j+1)**3-rad(j)**3)
 mu = ro*vs**2
 la = ro*vp**2 - 2.d0*mu
 ro = ro*1.d3
 la = la*1.d9
 mu = mu*1.d9

 !print*, rad(j+1), ro, la , mu
 rot(j) = ro
 lat(j) = la
 mut(j) = mu
  end do
   
   
 end subroutine earth_homogeneous_layers_prembased


 subroutine earth_nlayers(x0,ro,la,mu,GR,solid)
 use lnb_param
!---------------------------------------
! Gets planet properties at non-dimentionalized radius x0:
!ro: density
!la: Lame lambda constant
!mu: Lame mu constant (shear modulus)
!GR: gravity
!solid: boolean, true if solid, otherwise liquid
    integer                 :: i,j
    doubleprecision      :: la,mu,CST,GR
    doubleprecision      :: x,x0,xx0,ro
    doubleprecision      ::r1,r2
    logical       :: solid
    double precision, dimension(nlayer+1) :: r
    


    solid = .true.

    r(:) = radius(:)
    r(1) = 0.d0
    x    = x0   * ra
    CST = 4.d0*pi*GG/3.d0
    if ( x > r(nlayer+1) ) x = r(nlayer+1)
    
    GR = 0.d0
    
    do i = 1,nlayer

 if (x > r(i)) then 

  r2 = r(i+1)
  r1 = r(i)
  ro = roc(i)

  if (x <= r(i+1)) then

   GR = GR + ro*(x**3-r1**3)
   if (x > epsdb) GR = GR*CST/(x**2)

  else
   GR = GR + ro*(r2**3-r1**3)
  end if
 end if

 if ((x >= r(i)).and.(x < r(i+1))) then    
  la = lac(i)
  mu = muc(i)
  ro = roc(i)
  if ( soliddim(i).eqv..false. ) solid = .false.

 end if

    end do

    if (x>=r(nlayer+1)) then
 la = lac(nlayer)
 mu = muc(nlayer)
 ro = roc(nlayer)
  if ( soliddim(nlayer).eqv..false. ) solid = .false.

     end if

 end subroutine earth_nlayers
 
 subroutine earth_nlayers_viscoelas(x0,ro,la,mu,GR,solid)
 use lnb_param
!---------------------------------------
! Defines the earth rheology at non-dimensionalized radius x0 and angular frequency freq (global/public variable)

    integer                 :: i
    double complex   :: la, mu, cst
    double precision      :: GR,la0,mu00,ka,mu1,mu2,vi2
    double precision      :: x,x0,xx0,ro,vi
    logical       :: notfound,solid,burg



  call earth_nlayers(x0,ro,la0,mu00,GR,solid)
  
  x    = x0 * ra
 do i = 1,nlayer
 if ((x >= radius(i)).and.(x < radius(i+1))) then
 vi = vic(i) 
 vi2 = vic2(i)
 mu2=muc2(i)
 burg=burgers(i)
 end if
 end do

  if (burg.eqv..true.) then
 mu1=mu00

 if ((vi*freq) /= 0.d0) then
  ka=la0 + 2.d0/3.d0*mu00
  !print*,mu1,mu2,vi,vi2
   mu=mu1*freq*(freq+mu2/vi2)/((freq+mu2/vi2)*(freq+mu1/vi)+mu1/vi2*freq)
  !la=4*la0*mu**2/(2*mu*(3*la0+2*mu)-3*la0*mu)
  la=ka-2.d0/3.d0*mu
  !print*,freq,mu,la
 else
  mu = mu00
  la = la0
 endif
  else
   ka = la0 + 2.d0/3.d0*mu00
   cst = vi*freq
   if (cst /= 0.d0) then
  mu = mu00/(1.d0+mu00/cst)
  la = (la0 + mu00*ka/cst)/(1.d0 + mu00/cst)
   else
  mu = mu00
  la = la0
   endif
  endif
   
  !
  !!write(101,*), freq,cst,la0,la,mu00,mu,vi

 end subroutine earth_nlayers_viscoelas
 
 


end module model
