 
 !subroutine lnb_setup {{{
 subroutine lnb_setup(  ksr, ksi, hsr, hsi, lsr, lsi, lksr, lksi, dnfreq,  frequencies,  ddegmax, ddegmin, &
          dgo0,  dr0,  dmu0, dallow_layer_del, source_type, ddisplay, dnlayer,  dradius,    &
          dvic,  dlac,  dmuc, dvic2,  dmuc2,  droc, dburgers,  dsoliddim)

 use lnb_param
 use util
 implicit none

 integer :: dnfreq, ddegmax, ddegmin
 integer :: dnlayer
 integer :: source_type, fr, i, IERR
 double precision :: T,cst, fluid_tscale,g
 double precision, dimension(dnfreq,ddegmax+1) :: hsr,hsi,lsr,lsi,ksr,ksi
 double complex, dimension(dnfreq,ddegmax+1) :: hs,ls,ks
 double precision, dimension(dnfreq) :: frequencies

 double complex :: lks(dnfreq,(ddegmax+1)*(dnlayer+1)*6) 
 double precision, dimension(dnfreq,(ddegmax+1)*(dnlayer+1)*6) :: lksr, lksi 
 
 double precision   :: dmu0,dgo0,dr0
 logical :: dallow_layer_del, ddisplay

 double precision, dimension(dnlayer) :: droc,dlac,dmuc,dvic,dvic2,dmuc2 
 double precision, dimension(dnlayer+1) :: dradius 
 double precision, dimension(dnlayer) :: dsoliddim,dburgers 

 !assign local variables to public variables
 mu0=dmu0; go0=dgo0; r0=dr0; allow_layer_del=dallow_layer_del; 
 nfreq=dnfreq;degmax=ddegmax;degmin=ddegmin;nlayer=dnlayer; 
 allocate(roc(nlayer),lac(nlayer),muc(nlayer),vic(nlayer),vic2(nlayer),muc2(nlayer))
 roc=droc; lac=dlac; muc=dmuc; vic=dvic; vic2=dvic2; muc2=dmuc2; 
 allocate(radius(nlayer+1),soliddim(nlayer),burgers(nlayer),radbc(nlayer+1))
 radius=dradius; 

	display=ddisplay
	!print*,display
	do i=1,nlayer
	if (NINT(dsoliddim(i))==1) then
		soliddim(i)=.true.
	else
		soliddim(i)=.false.
	endif 
	if (NINT(dburgers(i))==1) then
		burgers(i)=.true.
	else
		burgers(i)=.false.
	endif 
	enddo

 pi=dacos(-1.d0)

 ra=radius(nlayer+1)
 rb = 3480.d3
 rc = 1221.5d3
 r0=ra
 cst = 365.25d0*24.d0*3600.d0*1000.d0 

 benchmark_spada=.false.
 !frequencies=1E-5/cst*2*pi

	hs(:,:)=(0.d0,0.d0)
	ks=hs
	ls=hs
 
 call love_numbers(frequencies,source_type,hs,ls,ks,lks)

 hsr=dble(hs)
 hsi=dimag(hs)
 ksr=dble(ks)
 ksi=dimag(ks)
 lsr=dble(ls)
 lsi=dimag(ls)
 !!! love kernels 
 lksr=dble(lks)
 lksi=dimag(lks)

end subroutine
! }}}

!subroutine love_numbers {{{
 subroutine love_numbers(frequencies,source_type,hs,ls,ks,lks)

 use model
 use lovenb_sub

 double complex, dimension(:,:), pointer :: bc
 double complex, dimension(:,:), pointer :: f
 character*40,     dimension(:),   pointer :: sourcs
 integer,          dimension(:),   pointer :: indx,indx2
 integer :: nbc,nfext,ifreq,ntheta2,nphi2,pp 
 double precision :: cst,prec, delta!,h1,h2,k1,k2,l1,l2,dh1,dh2,dk1,dk2,dl1,dl2
 double precision :: he,ke,le
 double precision :: frequencies(nfreq), fluid_tscale
 double complex :: loveh,lovel,lovek
 logical :: logi
 integer :: nbc_init
 integer :: i,j,k,n, fr,m,l, source_type
 double complex :: hs(nfreq,degmax+1), ls(nfreq,degmax+1), ks(nfreq,degmax+1)
 !double complex, dimension(nbc_init,1) :: lovekernels
 double complex, dimension(:,:), allocatable :: lovekernels
 double complex :: lks(nfreq,(degmax+1)*(nlayer+1)*6) 

 double precision :: cpu_time1(100)
 integer :: cpu_count

 pi=dacos(-1.d0)

  open(unit=101, file='lastrun_log')

 cpu_count=0
 cpu_count=cpu_count+1;call cpu_time(cpu_time1(cpu_count))
 !write(101,*), 'Done !', cpu_time1(cpu_count), 's'
 !write(101,*),''


 !double precision, dimension(:), pointer :: vech,vecl,veck,vecd,vec

 radbc(:)=.false.
 radbc(1)=.true.
 radbc(nlayer+1)=.true.
 do i=2,nlayer
  if ((.not.soliddim(i)).or.(.not.soliddim(i-1))) then
   radbc(i)=.true.
  end if
 end do
 

 ! Test number of interfaces and of boundary conditions
  nbc = 0
  lks=0.d0
  !!! >>> SA: 01/27/2018 => following changes are made to retrieve love numbers at depth 
  !do i=1,nlayer+1 
 !if (radbc(i)) nbc=nbc+6
  !enddo
  nbc_init = (nlayer+1)*6
  nbc=nbc_init
  !!! >>> SA 

 ifmin=source_type
 ifmax=source_type
 

  if (display) then
	 print*,'model init'
 endif


 call model_init2(nbc)  


  if (display) then
	 print*, 'done'
 endif

 nfext = 13 ! Number of potential excitation sources
 allocate( lovekernels(nbc_init,1) )
 allocate( bc(nbc,nbc_init), indx(nbc_init) )
 allocate( f(nfext,nbc_init), sourcs(nfext) )

 sourcs(1)='ICB --Volumetric Potential' 
 sourcs(2)='ICB --Pressure'
 sourcs(3)='ICB --Loading'
 sourcs(4)='ICB --Tangential Traction'
 sourcs(5)='CMB --Volumetric Potential'
 sourcs(6)='CMB --Pressure'
 sourcs(7)='CMB --Loading'
 sourcs(8)='CMB --Tangential Traction'
 sourcs(9)='SURF--Volumetric Potential'
 sourcs(10)='SURF--Pressure'
 sourcs(11)='SURF--Loading'
 sourcs(12)='SURF--Tangential Traction'

 
 if (display)  print*, 'source_type =', source_type, ', ', sourcs(source_type)
 if (source_type<9) then
  print*,'Error: Internal loading not supported yet, please input source_type between 9 and 12'
  print*,'Reference:'
  do i=1,12
   print*,i,sourcs(i)
  enddo
 return
 endif


! Print/write control
 !display=.false. ! Printing on the terminal
 
 ! Calculation optimization (avoids underflows at high degree)
 !allow_layer_del = .true.     ! Do we allow deletion of the central layers in the calculation if necessary?
 layerrap = 1.d0  ! Max absolute ratio between love numbers at the top and bottom of the last layer

 ! Mode search parameters
 freq0 = -1.d10
 logi = .true.
 cst = 365.25d0*24.d0*3600.d0*1000.d0 
 prec = 1.d-8
 firstmode = 0.d0

 
 cpu_count=cpu_count+1;call cpu_time(cpu_time1(cpu_count))
 !write(101,*), 'Earth model initialization :', cpu_time1(cpu_count) -cpu_time1(cpu_count-1), 's'
 !write(101,*),''
 if (display) then
 write(*,*), 'Earth model initialization :', cpu_time1(cpu_count) -cpu_time1(cpu_count-1), 's'
 write(*,*),''
 endif

! Spherical Harmonics initialization



! -- SH degree loop


 bc(:,:)=0.d0

 
 
  !write(101,*), 'Calculating impulse response'
 if (display) then
 write(*,*), 'Calculating impulse response'
 endif
 	if (maxval(frequencies).ne.0 .and. display) then
		print*, 'Frequency ranging from ', minval(dabs(frequencies)), &
		' to ', maxval(frequencies),' s^-1'
	endif 

	if (frequencies(1)==0 .and. display) then
		print*,'Elastic calculation'
	endif
 deg=degmin-1



 do l = degmin,degmax 
 deg=deg+1
 !write(*,*)
  !write(*,*)
  !write(101,fmt='(A19, 10X,A2,I3,A1,I3)',ADVANCE='NO'), repeat(char(8),19),'l=',deg,'/',degmax
  if (display) then
  write(*,fmt='(A10,A2,I5,A1,I5)',ADVANCE='NO'),repeat(char(32),10),'l=',deg,'/',degmax
  endif

  !write(*,*)
  !!write(101,*), fr, T

 do fr=1,nfreq

  !!write(101,*),fr


  !the way Fourier Transform algorithms works : they calculate for fmin=0 to fmax=1/dt
  !by aliasing for this algo, f=1/(2*dt) to f=1/dt is the same as f=-1/(2*dt) to f=0
  !but physics-wise, the Earth response at high frequency is not what we are looking for
  !the [0 ; 1/(2*dt)] interval and its negative counterpart are what we are looking for
  !therefore in our love number calculation the frenquency set must be [0 ; 1/(2*dt) ] U [-1/(2*dt) 0[

  !if (fr==1) then ! the variable freq is, in fact, the pulsation = 2*pi*frequency
  ! fluid_tscale=-575.502 / (deg+91.1765) -0.176471 ! this formula sets fluid time scale = 1e6 kyr at degree 2, 1e5 kyr at degree 20 and 1e3 kyr at degree 90, which allows no error growth and still the fluid number as the benchmark at all degree (tested up to degmax=147)
   !print*,fluid_tscale
  ! freq=complex(0.d0,1.d0)*2.d0*pi*(10.d0**(fluid_tscale)/cst) ! empiric time period that changes according to harmonic degree to ensure getting fluid response without crazy error growing for high degree
! it is meant to approximate the theoretical freq=0.d0 that should be assessed for fft
  !elseif (fr <= nfreq/2+1) then
  ! freq=complex(0.d0,1.d0)*(dble(fr-1)/T)*2.d0*pi
  !else 
  ! freq=complex(0.d0,1.d0)*(-dble(nfreq+1-fr)/T)*2.d0*pi
  !end if
  
  freq=dcmplx(0.d0,frequencies(fr))
  !print*,fr,freq

  !if (deg==degmin) write(31,*),fr,real(freq)*cst,aimag(freq)*cst,time(fr)
  !freq=dble(fr)/T
  !--  Elastic Love number calculation
  !freq=2.d0*pi/(1e4*cst)*complex(0.d0,1.d0)

   !if (deg==2) then ! for rotationnal feedback
   !ifmin=9;ifmax=9 ! Sets tidal calculation
   !bc(:,:)=0.d0
   !call boundary_conditions_matrix(bc,indx,nbc)
   !call external_forcing(deg,f,sourcs,nfext,nbc) 
   !call solution(deg,bc,indx,f,sourcs,nfext,nbc,loveh,lovel,lovek,delta)
   !k2tidal(fr)=lovek-1.d0
   !h2tidal(fr)=loveh
   !l2tidal(fr)=lovel
   !ifmin=11;ifmax=11 ! Sets back loading calculation
   !endif
   !print*,(deg)
   bc(:,:)=0.d0
   lovekernels(:,1)=0.d0
   call boundary_conditions_matrix(bc,indx,nbc)
	!print*,'bc ok', layerrap
   call external_forcing(deg,f,sourcs,nfext,nbc) 
	!print*,'forcing ok', layerrap
   call solution(deg,bc,indx,f,sourcs,nfext,nbc,loveh,lovel,lovek,delta,lovekernels)

   !do i=1,nbc
   !   lks(fr,deg*nbc+i) = lovekernels(i)
   !end do

	!print*,'first_sol', layerrap, epsdb

   ! Automatic reduction of the number of layers when the attenuation with depth becomes too strong
   allow_layer_del=.true.
   !allow_layer_del=.false.
   if (allow_layer_del.eqv..true.) then

      print*, layerrap, epsdb, nbc 
      !do while ((layerrap<=epsdb).and.(nbc>12*6).or.(isnan(layerrap)) )
      do while ((layerrap<=epsdb).and.(nbc>12).or.(isnan(layerrap)) )
     !do while ((layerrap<=epsdb).and.(nbc>12*(nlayer+1)).or.(isnan(layerrap)) ) ! SA: for love Kernels
 !		print*,'trying to delete layer', layerrap, epsdb, nbc
     !write(101,*)
     !write(101,*) 'Rapport Nombre de Love surface/profondeur faible : ', layerrap 
     !write(101,*) ' Changement d''interface de debut d''integration' 

     if (display) then
     write(*,*)
     write(*,*) 'Surface/Depth Love number ratio small: ', layerrap 
     write(*,*) ' Changing the interface where the integration starts'
     endif

     nbc = nbc-6
     n = 1
     do while (.not.radbc(n))
      n = n+1
     end do
     radbc(n) =.false.
     radbc(n+1) =.true.
     !write(101,*) ' New start interface: ', radius(n+1)/1.d3,' km'
     if (display) then
     write(*,*) ' New start interface: ', radius(n+1)/1.d3,' km'
     endif
        
     deallocate( bc, indx, f, sourcs, lovekernels )
     allocate( bc(nbc,nbc), indx(nbc) )
     allocate( f(nfext,nbc), sourcs(nfext) )
    
     allocate( lovekernels(nbc,1) )

     bc(:,:)=0.d0
     lovekernels(:,1)=0.d0

     call boundary_conditions_matrix(bc,indx,nbc)
     call external_forcing(deg,f,sourcs,nfext,nbc) 
     call solution(deg,bc,indx,f,sourcs,nfext,nbc,loveh,lovel,lovek,delta,lovekernels)

    end do

      !do i=1,nbc
      !   pp = ((nlayer+1)*6 - nbc)/6 
      !   lks(fr, deg*nbc + pp*6 + i) = lovekernels(i)
      !end do

   end if

   !-- Saving Love numbers
    hs(fr,deg+1) = loveh
    ks(fr,deg+1) = lovek - 1.d0
    ls(fr,deg+1) = lovel
    
    !!! and love kernels 
    do i=1,nbc
      !lks(fr,deg*nbc+i) = lovekernels(i,1)
      lks(fr,(deg+1)*nbc_init+i-nbc) = lovekernels(i,1)
    end do

	!print*,dble(hs(fr,deg+1)), dimag(hs(fr,deg+1)),dble(ks(fr,deg+1)), dimag(ks(fr,deg+1))

    !if (fr==1) then ! if this is supposed to be fluid response
    !  ! then cut off the imaginary part, which is inherited from approximating t=infinity to a few million years (see above the setting of freq)
    ! hs(fr,deg+1)=real(hs(fr,deg+1))
    ! ks(fr,deg+1)=real(ks(fr,deg+1))
    ! ls(fr,deg+1)=real(ls(fr,deg+1))
    !endif
!~~~~~~~~~~~~~~TEST ZONE~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 if (benchmark_spada) then
 write(17,*), deg,fr,dble(hs(fr,deg+1)), dimag(hs(fr,deg+1)),dble(ks(fr,deg+1)), dimag(ks(fr,deg+1))
 endif
!~~~~~~~~~~~~~~END TEST ZONE~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


  !write(*,*) 'Elastic Love Numbers:  ',ke,he,le 
  !write(ifile1,*) deg,ke,he,le
 end do
 
 end do


 cpu_count=cpu_count+1;call cpu_time(cpu_time1(cpu_count))
 !write(101,*), 'Earth model initialization :', cpu_time1(cpu_count) -cpu_time1(cpu_count-1), 's'
 !write(101,*),''
 if (display) then
 write(*,*), 'Love number calc done :', cpu_time1(cpu_count) -cpu_time1(cpu_count-1), 's'
 write(*,*),''
 endif

 deallocate( bc, indx, f, sourcs )
end subroutine
! }}} 

