module lovenb_sub

 use util
 use lnb_param

 implicit none
 
 contains
 
 ! subroutine boundary_conditions_matrix {{{ 
 subroutine boundary_conditions_matrix(bc,indx,nbc)

  implicit none
  double complex :: la, mu
  double precision :: xmin,xmax,eps,hini,hmin
  double precision :: x,ro1,ro2,g,d,one
  double complex, dimension (:,:), pointer :: bc
  integer, dimension(:), pointer   :: indx
  double complex, dimension(6)   :: ystart,dydx
  double complex, dimension(6,6) :: ystep
  integer :: ny,i,j,k,ii,jj,kk,nok,nbad,is,nbc,ibc,ifirst,nstep
  logical :: solid,vali
  
  !-- Initialisation of boundary condition matrix

  indx(:) = 0
  
  !-- Integration 

  eps= 1.d-10
  hini=1.d-10
  hmin=0.d0

  ! Underflows risk tests in the numerical integration
  ! The starting ystart are chosen to be different than 1 to limit under and overflows in the Runge-Kutta
  if (-300.d0<dlog10(tiny(ra))) then 
	write(*,*) 'Warning problem in the dimensioning of the starting values in Runge-Kutta'
	write(*,*) 'log tiny(ra) = ', dlog(tiny(ra))
	stop
  endif

	
   if ((dlog10(1.d0/ra**dint(dble(deg)/10.d0)))<-250.d0) then
	!valini=1.d-300!*complex(1.d0,1.d0)
	   valini=10.d0**(-dble(150)/dble(nbc/6-1))!*complex(1.d0,1.d0)
   else
	   valini= 1.d0!*complex(1.d0,1.d0)
   endif

  ! Which among the model interfaces is the first to display boundary conditions
  ifirst = 0
  i = 0
  do while (ifirst==0)
	i = i+1
	if (radbc(i)) ifirst = i
  end do
  
  
  !-- Layer iteration for the solutions propagation 
  !    and typing the matrix boundary conditions 
  ibc = 0   ! Counter of interfaces with boundary conditions 
 
  do i = ifirst,nlayer 
  
   xmin= radius(i)/ra
   xmax= radius(i+1)/ra - epsdb
   nstep=100
  
   if (soliddim(i)) then
   		ny = 6
   		is = 0
   		one= 1.d0
   else	
   		ny = 2
   		is = 4
   		one= -1.d0
   end if 	
   
   do j = 1,ny
		! Typing ystart 
      !!! >>> SA: 01/27/2018 => changes made to retrieve love numbers at depth 
		!!!if (radbc(i)) then
   		ystart(:)= 0.d0
			ystart(j)= valini
 			k=ibc+1
 		!!!else
 		!!!	ystart(:)=ystep(j,:)
 		!!!end if
 		
		! Numerical Integration 
		!print*,j
		call propagFDS(ystart,xmin,xmax,nstep)
 		!call odeint(ystart,ny,xmin,xmax,eps,hini,hmin,nok,nbad)
		!print*,j
 		
 		! Boundary Condition matrix - propagation part 
		!!!if (radbc(i+1)) then
 			ii = 6*k+is
 			jj = 6*k+j+is-3
 			do kk=1,ny
            bc(ii+kk,jj) = ystart(kk)*one !/valini 
			end do
		!!!else
		!!!	ystep(j,:)=ystart(:) 
 		!!!end if
   end do
	
	
   !!!if (radbc(i)) then
   
	! Boundary Condition matrix - solid regions
	if (soliddim(i)) then
		one = -1.d0
		if (i>1)then 
			if (.not.soliddim(i-1)) one = 1.d0
		end if
		do  j=1,6
			bc(j+6*ibc,j+6*ibc+3) = one
		end do
		
	! Boundary Condition matrix - liquid regions
	else
		x = radius(i)/ra
		call earth_nlayers_viscoelas(x,ro2,la,mu,g,solid)
		ii = 6*ibc+1
		bc(ii,ii+3) = -1.d0
		bc(ii,ii+4+3) = -go0/g
		bc( ii+1,ii+3)=-ro2*g*ra/mu0
		bc( ii+2,ii+1+3)=-1.d0
		bc(ii+5,ii+3)= 4.d0*pi*GG*ro2*ra/go0
		bc(ii+4,ii+4+3)=-1.d0
		bc(ii+5,ii+5+3)=-1.d0
		x = radius(i+1)/ra
		call earth_nlayers_viscoelas(x-epsdb,ro1,la,mu,g,solid)
		call earth_nlayers_viscoelas(x,ro2,la,mu,g,solid)
		ii = 6*(ibc+1)+1
		bc(ii,ii-1)=-1.d0
		bc(ii,ii+1)=bc(ii+4,ii+1)*go0/g ! b(17,14) solution integration 1 of z5 CMB
		bc(ii,ii+2)=bc(ii+4,ii+2)*go0/g ! b(17,15) solution integration 2 of z5 CMB
		 					  ! b(13,..) y1 CMB
		bc(ii+1,ii-1)=-ro1*g*ra/mu0
		bc(ii+2,ii)=-1.d0
		bc(ii+5,ii-1)= 4.d0*pi*GG*ro1*ra/go0
	end if	
	ibc = ibc+1
	
   !! end if
	
  end do

  !-- Internal sphere

  call sphere(ifirst,deg,bc)

  !-- Surface conditions
	bc(nbc-5,nbc-2)=-1.d0
	bc(nbc-3,nbc-1)=-1.d0
	bc(nbc-1,nbc)=-1.d0
	bc(nbc,nbc)=dble(deg+1)

  !-- Degree 1 special case
	if(deg==1) then
		  bc(nbc,:)=0.d0
		  bc(nbc,nbc)=1.d0
	end if

	!print*, bc, imag_eval
 end subroutine boundary_conditions_matrix
  ! }}} 
 
 ! subroutine sphere {{{
 subroutine sphere(i,n,bc)
!------------------------

	implicit none
	double complex :: la,mu
	integer :: n,i
	double precision :: x,e,r,ro,g,cst,xsen,eps
	double complex, dimension(:,:),pointer :: bc

	logical :: solid
	

	eps= 1.d-10
	x = radius(i)/ra
	e = x
	call earth_nlayers_viscoelas(x-eps,ro,la,mu,g,solid)
	
	xsen=(x/e)**n
	r=x*ra
	cst = 4.d0*pi*GG*ro

	!write(89,*),xsen,r,ra,cst,mu0,mu,n	


	bc(1,1)=xsen*r/ra
	bc(1,2)=xsen/(r*ra)
	bc(1,3)=0.d0

	bc(2,1)=(2.d0*mu*dble(n*n-n-3)/dble(n) + cst/3.d0*ro*r**2)*xsen/mu0
	bc(2,2)=(2.d0*mu*dble(n-1)/r**2 + cst/3.d0*ro)*xsen/mu0
	bc(2,3)=-ro*xsen/mu0

	bc(3,1)=dble(n+3)/dble(n*(n+1))*xsen*r/ra
	bc(3,2)=1.d0/dble(n)*xsen/(r*ra)
	bc(3,3)=0.d0

	bc(4,1)=2.d0*mu*dble(n+2)/dble(n+1)*xsen/mu0
	bc(4,2)=2.d0*mu*dble(n-1)/dble(n)*xsen/(r**2*mu0)
	bc(4,3)=0.d0

	bc(5,1)=0.d0
	bc(5,2)=0.d0
	bc(5,3)=xsen/(go0*ra)

	bc(6,1)=-cst*xsen*r/go0
	bc(6,2)=-cst*xsen/(r*go0)
	bc(6,3)=dble(n)*xsen/(r*go0)


 end subroutine sphere
 ! }}} 
 
 ! subroutine external_forcing {{{
 subroutine external_forcing(n,source,sourcs,nfext,nbc)
!-----------------------------------

   implicit none
   integer :: n
   integer :: nfext,nbc
   character*40,     dimension(:),   pointer :: sourcs
   double complex, dimension(:,:), pointer :: source
   
   source(:,:)=0.d0
   sourcs(:)=' '
	  
   ! In Case of a Inner Core - Outer Core - Mantle planet and Boundary conditions on these 3 interfaces
   if (nbc == 24) then 
   
	!-- Sources at the ICB
   
	sourcs(1)='ICB --Volumetric Potential'
	source(1,12)=dble(n)/(rc*go0)
	source(1,11)=1.d0/(ra*go0)
	
	sourcs(2)='ICB --Pressure'
	source(2, 8)=-ro_mean/mu0

	sourcs(3)='ICB --Loading'
	source(3, 8)=-ro_mean*dble(2*n+1)/(3.d0*mu0)*ra/rc
	source(3,12)= dble(2*n+1)/(rc*go0)

	sourcs(4)='ICB --Tangential Traction'
	source(4,10)= ro_mean/mu0

	!--Sources at the CMB

	sourcs(5)='CMB --Volumetric Potential'
	source(5,14)=-ro_mean/mu0*ra/rb
	source(5,18)= dble(2*n+1)/(rb*go0)

	sourcs(6)='CMB --Pressure'
	source(6,14)=-ro_mean/mu0

	sourcs(7)='CMB --Loading'
	source(7,14)=-ro_mean*dble(2*n+1)/(3.d0*mu0)*ra/rb
	source(7,18)= dble(2*n+1)/(rb*go0)

	sourcs(8)='CMB --Tangential Traction'
	source(8,16)=-ro_mean/mu0
	
    endif

    !--Sources at the surface

	sourcs(9)='SURF--Volumetric Potential'
	source(9,nbc)=dble(2*n+1)/(ra*go0)

	sourcs(10)='SURF--Pressure'
	source(10,nbc-4)=-ro_mean/mu0

	sourcs(11)='SURF--Loading'
	source(11,nbc-4)=-ro_mean*dble(2*n+1)/(3.d0*mu0)
	source(11,nbc)= dble(2*n+1)/(ra*go0)

	sourcs(12)='SURF--Tangential Traction'
	source(12,nbc-2)= ro_mean/mu0

	!--Special case n=1

	if (n.eq.1) then
	  source(11,nbc)=0.d0
	  source(9,nbc)=0.d0
	end if
	
 
 end subroutine external_forcing
 ! }}} 
 
 ! subroutine solution {{{ 
 subroutine solution(n,bc,indx,f,sourcs,nfext,nbc,loveh,lovel,lovek,delta,love_kernels)
!-------------------------------------------------

  integer :: n,nfext,nbc,i,j,ic,ibc,ii,jj,icmin,ifirst,info,lda,ldb
  double complex, dimension(:,:), pointer :: bc,f
  integer,      dimension(:), pointer :: indx
  character*40, dimension(:), pointer :: sourcs	
  double complex, dimension(nbc,1) :: sc, love_kernels
  double complex :: loveh,lovel,lovek,loveh1,lovel1,lovek1,loveh1s,lovel1s,lovek1s
  double complex :: sumh,suml,sumk,d,val,sumy2,sumy4,sumy6,rads
  double precision :: delta
  double complex, dimension(6) :: icbh,icbl,icbk,y2srf,y4srf,y6srf
  double complex, dimension(nbc,nbc) :: bcsav
  integer :: ipiv(nbc)
  logical :: ok

  bcsav(:,:) = bc(:,:)
  rads=0.d0
  lda=nbc
  ldb=nbc

 ! open(unit=370,file='../temp') 
 !  do i=1,nbc
 !     do j=1,nbc
 !        write(370,*), real(bc(i,j))
 !     end do 
 !  end do
 ! close(370) 

  ifirst = 0
  i = 0
  do while (ifirst==0)
	i = i+1
	if (radbc(i)) ifirst = i
  end do
! do i = 1,nbc 
!	do j = 1,nbc
!	if (bc(i,j)/=0.d0) print*, i,j,bc(i,j)
!	enddo
! enddo
! write(*,*)
  
  !-- LU decomposition and determinant
    ok = .true.
   
  !call ludcmp(bc,nbc,nbc,indx,d)
  !do i=1,nbc
  !   d=d*bc(i,i)
  !end do
  !if (display) write(*,*)
  !write(*,*),dlog(freq),'   Matrix det = ',d/dabs(d)*dlog(dabs(d))
  !if (display) write(*,*),dlog10(freq*365.25d0*24.d0*3600.d0*1000.d0),'   Matrix det = ',d
  !write(ifile5,*) freq,d

!-- Resolution

  if (ok) then
  do i = ifmin,ifmax
  	
	do j=1,nbc
	  sc(j,1)=f(i,j)
	end do

	!if (display) write(*,*) 'Source: ',sourcs(i)
	!call lubksb(bc,nbc,nbc,indx,sc)
	call ZGESV(nbc,1,bc,lda,ipiv,sc,ldb,info)

  	if (info.ne.0) then 
		print*, 'Error in ZGESV : LAPACK linear equation solver couldn''t resolve the system'
	end if
   !!! >>> y1, y3, y5 at surface (by default) 
	!loveh = sc(nbc-2,1)*ra*go0 !go_surf
	!lovel = sc(nbc-1,1)*ra*go0 !go_surf
	loveh = sc(nbc-2,1)*ra*go_surf
	lovel = sc(nbc-1,1)*ra*go_surf
	lovek = sc(nbc,1)*ra*go0
	delta = (1.d0-dble(n+1)/dble(n)*(lovek-1.d0)+2.d0/dble(n)*loveh)

   !print *, ra, go0, go_surf

	do j=1,nbc
      if (j<nbc-2) then ! no need for valini restoration for surface. 
         sc(j,1) = sc(j,1)*valini**((nbc-j-3)/6+1)
      endif
	   love_kernels(j,1)=sc(j,1) !*valini
	end do

	sumy2=0.d0
	sumy4=0.d0
	sumy6=0.d0

	do ic=1,6
	  sumy2 = sumy2+bcsav(nbc-4,ic+nbc-9)*sc(ic+nbc-9,1)
	  sumy4 = sumy4+bcsav(nbc-2,ic+nbc-9)*sc(ic+nbc-9,1)
	  sumy6 = sumy6+bcsav(nbc,ic+nbc-9)*sc(ic+nbc-9,1)
	end do


   !print*, '***********' 
   !print*, radbc
   !print*, '***********' 

   loveh1 = sc(4,1)
   lovel1 = sc(6,1)
   lovek1 = sc(8,1) - (radius(nlayer-nbc/6+3)/ra)**deg/(go0*ra)

   loveh1s = sc(nbc-2,1)
   lovel1s = sc(nbc-1,1)
   lovek1s = sc(nbc,1) - 1.d0/(go0*ra)

   ! ratio of love_number_at_depth vs love_number_at_surface was compared 
   ! the comparison is suited for two solid-layers. 
   ! Not sure whether it works as well when we have a liquid layer at depth. 
   layerrap = zabs(loveh1/loveh1s) 
   if (zabs(lovel1/lovel1s) < layerrap) then 
      layerrap = zabs(lovel1/lovel1s) 
   endif 
   if (zabs(lovek1/lovek1s) < layerrap) then 
      layerrap = zabs(lovek1/lovek1s) 
   endif 

   end do
   end if
  
   !open(unit=371,file='../temp2') 
   !do i=1,nbc
   !   write(371,*), real(sc(i,1))
   !end do
   !close(371) 

 !301 format(a4,5x,'n=',i3,5x,'h=',f14.10,5x,'l=',f14.10,5x,'k=',f14.10)
 301 format(a4,5x,'n=',i5,5x,'h=',d14.7,5x,'l=',d14.7,5x,'k=',d14.7)
 305 format(a4,i1,4x,'n=',i5,5x,'h=',d14.7,5x,'l=',d14.7,5x,'k=',d14.7)
 306 format(a4,i2,3x,'n=',i5,5x,'h=',d14.7,5x,'l=',d14.7,5x,'k=',d14.7)
 302 format(i5,5x,f14.10,5x,f14.10,5x,f14.10)
 303 format(i5,5x,f14.10,5x,f14.10,5x,f14.10,5x,f14.10)
 304 format(i5,5x,d20.10,5x,d20.10,5x,d20.10)

 end subroutine solution 
 ! }}}
 
  ! subroutine solution_nbloutput {{{ 
  subroutine solution_nbloutput(n,bc,indx,f,sourcs,nfext,nbc,loveh,lovel,lovek,delta)
!-------------------------------------------------

  integer :: n,nfext,nbc,i,j,ic,ibc,ii,jj,icmin,ifirst
  double precision, dimension(:,:), pointer :: bc,f
  integer,      dimension(:), pointer :: indx
  character*40, dimension(:), pointer :: sourcs	
  double precision, dimension(nbc) :: sc
  double precision :: loveh,lovel,lovek,loveh1,lovel1,lovek1,loveh1s,lovel1s,lovek1s
  double precision :: sumh,suml,sumk,d,val,delta,sumy2,sumy4,sumy6,rads
  double precision, dimension(6) :: icbh,icbl,icbk,y2srf,y4srf,y6srf
  double precision, dimension(nbc,nbc) :: bcsav
  logical :: ok

  
  !do ic=1,6
!	if (nbc==18) then 
!		icbh(ic) = bc(1,ic)
!		icbl(ic) = bc(3,ic)
!		icbk(ic) = bc(5,ic)
!!	elseif (nbc==24) then 
!		icbh(ic) = bc(7,ic+3)
!		icbl(ic) = bc(9,ic+3)
!		icbk(ic) = bc(11,ic+3)
!	end if
!	y2srf(ic) = bc(nbc-4,ic+nbc-9)
!	y4srf(ic) = bc(nbc-2,ic+nbc-9)
!	y6srf(ic) = bc(nbc,ic+nbc-9)
 ! end do
  
  bcsav(:,:) = bc(:,:)
  rads=0
  ifirst = 0
  i = 0
  do while (ifirst==0)
	i = i+1
	if (radbc(i)) ifirst = i
  end do
 !do i = 1,nbc 
!	do j = 1,nbc
!	if (bc(i,j)/=0.d0) print*, i,j,bc(i,j)
!	enddo
 !enddo
 !write(*,*)
  
  !-- LU decomposition and determinant
    ok = .true.
   
  call ludcmp(bc,nbc,nbc,indx,d)
  do i=1,nbc
     d=d*bc(i,i)
  end do
  write(*,*)
  !write(*,*),dlog(freq),'   Matrix det = ',d/dabs(d)*dlog(dabs(d))
 
!-- Resolution
 
  if (ok) then
  do i = ifmin,ifmax
  	
	do j=1,nbc
	  sc(j)=f(i,j)
	end do

	write(*,*) 'Source: ',sourcs(i)
	call lubksb(bc,nbc,nbc,indx,sc)

	loveh = sc(nbc-2)*ra*go_surf
	lovel = sc(nbc-1)*ra*go_surf
	lovek = sc(nbc)*ra*go0
	delta = 1.d0-dble(n+1)/dble(n)*(lovek-1.d0)+2.d0/dble(n)*loveh
	
	sumy2=0.d0
	sumy4=0.d0
	sumy6=0.d0

	do ic=1,6
	  sumy2 = sumy2+bcsav(nbc-4,ic+nbc-9)*sc(ic+nbc-9)
	  sumy4 = sumy4+bcsav(nbc-2,ic+nbc-9)*sc(ic+nbc-9)
	  sumy6 = sumy6+bcsav(nbc,ic+nbc-9)*sc(ic+nbc-9)
	end do
	
	!write(*,*) 'SURF',n,loveh,lovel,lovek-1.d0
	
	ibc = 0
	do  j = nlayer,1,-1
		if (radbc(j)) then 
			ibc=ibc+1
			if (soliddim(j)) then
				loveh1 = sc(nbc - ibc*6 -3 +1)*ra*go_surf*valini
				lovel1 = sc(nbc -ibc*6 -3 +3)*ra*go_surf*valini
				lovek1 = sc(nbc -ibc*6 -3 +5)*ra*go0*valini
			else 
				sumh=0.d0
				suml=0.d0
				sumk=0.d0
				ii = nbc - (ibc+1)*6
				jj = nbc - (ibc+1)*6 -3
				if (j==ifirst) then 
					icmin = 4
				else
					icmin = 1
				end if
				do ic=icmin,6
					!print*, ii+1,jj+ic
					sumh = sumh+bcsav(ii+1,jj+ic)*sc(jj+ic)
					suml = suml+bcsav(ii+3,jj+ic)*sc(jj+ic)
					sumk = sumk+bcsav(ii+5,jj+ic)*sc(jj+ic)
				end do
				loveh1 = sumh*ra*go_surf*valini
				lovel1 = suml*ra*go_surf*valini
				lovek1 = sumk*ra*go0*valini
			end if			
			if (j==ifirst) then
				layerrap = dabs(loveh1s/loveh)
				if (layerrap > dabs(lovel1s/lovel)) layerrap = dabs(lovel1s/lovel)
				if (layerrap > dabs((lovek1s-(rads/ra)**deg)/(lovek-1.d0))) &
				layerrap = dabs((lovek1s-(rads/ra)**deg)/(lovek-1.d0)) 
			endif

			loveh1s = loveh1
			lovel1s = lovel1
			lovek1s = lovek1
			rads = radius(j)
			!if (radius(j)==rb) then 
			!	write(*,301) 'CMB ',n,loveh1,lovel1,lovek1-(rb/ra)**deg
			!else if (radius(j)==rc) then
			!	write(*,301) 'ICB ',n,loveh1,lovel1,lovek1-(rc/ra)**deg
			!else
			!	if (j<10) then
			!		write(*,305) 'ITF-',j,n,loveh1,lovel1,lovek1-(radius(j)/ra)**deg
			!	else
			!		write(*,306) 'ITF-',j,n,loveh1,lovel1,lovek1-(radius(j)/ra)**deg
			!	end if
			!endif
		end if
	end do
			
  end do
  end if
  
   !write(*,*),dlog10(aimag(freq)*365.d0*24.d0*3600.d0*1000.d0),'   Matrix det = ',d, lovek
  
 end subroutine solution_nbloutput 
 ! }}}

end module lovenb_sub
