Index: /issm/trunk-jpl/src/c/modules/FourierLoveCorex/love_numbers.f90 =================================================================== --- /issm/trunk-jpl/src/c/modules/FourierLoveCorex/love_numbers.f90 (revision 25294) +++ /issm/trunk-jpl/src/c/modules/FourierLoveCorex/love_numbers.f90 (revision 25295) @@ -90,5 +90,5 @@ character*40, dimension(:), pointer :: sourcs integer, dimension(:), pointer :: indx,indx2 - integer :: nbc,nfext,ifreq,ntheta2,nphi2 + 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 @@ -96,8 +96,9 @@ 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((nlayer+1)*6,1) :: lovekernels + !double complex, dimension(nbc_init,1) :: lovekernels + double complex, dimension(:,:), allocatable :: lovekernels double complex :: lks(nfreq,(degmax+1)*(nlayer+1)*6) @@ -107,6 +108,6 @@ 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)) @@ -129,9 +130,11 @@ ! 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 = (nlayer+1)*6 + nbc_init = (nlayer+1)*6 + nbc=nbc_init !!! >>> SA @@ -153,6 +156,7 @@ nfext = 13 ! Number of potential excitation sources - allocate( bc(nbc,nbc), indx(nbc) ) - allocate( f(nfext,nbc), sourcs(nfext) ) + 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' @@ -187,5 +191,5 @@ !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 @@ -212,4 +216,5 @@ bc(:,:)=0.d0 + @@ -227,5 +232,4 @@ endif deg=degmin-1 - @@ -242,7 +246,4 @@ !write(*,*) !!write(101,*), fr, T - - - do fr=1,nfreq @@ -289,4 +290,5 @@ !print*,(deg) bc(:,:)=0.d0 + lovekernels(:,1)=0.d0 call boundary_conditions_matrix(bc,indx,nbc) !print*,'bc ok', layerrap @@ -295,10 +297,19 @@ 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 - do while ((layerrap<=epsdb).and.(nbc>12).or.(isnan(layerrap)) ) + 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,*) @@ -318,4 +329,5 @@ end do radbc(n) =.false. + radbc(n+1) =.true. !write(101,*) ' New start interface: ', radius(n+1)/1.d3,' km' if (display) then @@ -323,10 +335,13 @@ endif - deallocate( bc, indx, f, sourcs ) + 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) @@ -334,4 +349,10 @@ 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 @@ -343,10 +364,7 @@ !!! and love kernels do i=1,nbc - lks(fr,deg*nbc+i) = lovekernels(i,1) + !lks(fr,deg*nbc+i) = lovekernels(i,1) + lks(fr,(deg+1)*nbc_init+i-nbc) = lovekernels(i,1) end do - !!! for ks... 6th variable is y5 => love_k - do i = 1,nbc/6 - lks(fr,deg*nbc+6*i) = lks(fr,deg*nbc+6*i) - (radius(i)/ra)**deg - end do !print*,dble(hs(fr,deg+1)), dimag(hs(fr,deg+1)),dble(ks(fr,deg+1)), dimag(ks(fr,deg+1)) Index: /issm/trunk-jpl/src/c/modules/FourierLoveCorex/lovenb_sub.f90 =================================================================== --- /issm/trunk-jpl/src/c/modules/FourierLoveCorex/lovenb_sub.f90 (revision 25294) +++ /issm/trunk-jpl/src/c/modules/FourierLoveCorex/lovenb_sub.f90 (revision 25295) @@ -42,10 +42,10 @@ if ((dlog10(1.d0/ra**dint(dble(deg)/10.d0)))<-250.d0) then - valini=1.d-300!*complex(1.d0,1.d0) + !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) + valini= 1.d0!*complex(1.d0,1.d0) endif - ! Which among the model interfaces is the first to display boundary conditions ifirst = 0 @@ -60,5 +60,5 @@ ! and typing the matrix boundary conditions ibc = 0 ! Counter of interfaces with boundary conditions - + do i = ifirst,nlayer @@ -81,5 +81,5 @@ !!! >>> SA: 01/27/2018 => changes made to retrieve love numbers at depth !!!if (radbc(i)) then - ystart(:)=0.d0 + ystart(:)= 0.d0 ystart(j)= valini k=ibc+1 @@ -99,5 +99,5 @@ jj = 6*k+j+is-3 do kk=1,ny - bc(ii+kk,jj) = ystart(kk)*one + bc(ii+kk,jj) = ystart(kk)*one !/valini end do !!!else @@ -149,5 +149,4 @@ end do - !-- Internal sphere @@ -165,5 +164,4 @@ bc(nbc,nbc)=1.d0 end if - !print*, bc, imag_eval @@ -317,5 +315,4 @@ logical :: ok - bcsav(:,:) = bc(:,:) rads=0.d0 @@ -323,4 +320,11 @@ 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 @@ -362,8 +366,10 @@ call ZGESV(nbc,1,bc,lda,ipiv,sc,ldb,info) - if (info.ne.0) then + 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 @@ -371,47 +377,12 @@ delta = (1.d0-dble(n+1)/dble(n)*(lovek-1.d0)+2.d0/dble(n)*loveh) - !!! >>> at depth [by SA on 1/31/2018] - !!! sc outputs as follows (so scale it accordingly) - !!! y4 y2 y6 y1 y3 y5 - !!! at center of mass of the earth (top row) => at surface (bottom row) - !print *, nbc - do j = 1,nbc/6 - !!! y4, y2, y6 are not scaled yet =>>> ask lambert ==>>> - !love_kernels((j-1)*6+1,1) = sc((j-1)*6+1,1) ! y4 NOT SCALED YET - !love_kernels((j-1)*6+2,1) = sc((j-1)*6+2,1) ! y2 NOT SCALED YET - !love_kernels((j-1)*6+3,1) = sc((j-1)*6+3,1) ! y6 NOT SCALED YET - love_kernels((j-1)*6+4,1) = sc((j-1)*6+4,1)*ra*go_surf ! y1 => love_h - love_kernels((j-1)*6+5,1) = sc((j-1)*6+5,1)*ra*go_surf ! y3 => love_l - love_kernels((j-1)*6+6,1) = sc((j-1)*6+6,1)*ra*go0 ! y5 => love_k - !if (deg==2) then - ! print *, real(love_kernels((j-1)*6+1,1)), real(love_kernels((j-1)*6+4,1)) - !end if - end do - ! some reshuffling is needed for y4, y2, y6... - do j = 2,nbc/6 - !!! y4, y2, y6 are not scaled yet =>>> ask lambert ==>>> - love_kernels((j-2)*6+1,1) = sc((j-1)*6+1,1) ! y4 NOT SCALED YET - love_kernels((j-2)*6+2,1) = sc((j-1)*6+2,1) ! y2 NOT SCALED YET - love_kernels((j-2)*6+3,1) = sc((j-1)*6+3,1) ! y6 NOT SCALED YET - end do - !!! these are surface values... - love_kernels(nbc-5,1) = sc(1,1) ! y4 NOT SCALED YET - love_kernels(nbc-4,1) = sc(2,1) ! y4 NOT SCALED YET - love_kernels(nbc-3,1) = sc(3,1) ! y4 NOT SCALED YET - - !print *, ' ************* ' - !if (deg==2) then - ! do ibc=1,nbc/6 - !print *, real(sc( (ibc-1)*6+1:(ibc-1)*6+6 , 1)) - ! print *, real(love_kernels( (ibc-1)*6+1:(ibc-1)*6+6 , 1)) - ! end do - !print *, ' ************* ' - !do ibc=1,nbc/6 - ! print *, real(love_kernels( (ibc-1)*6+1:(ibc-1)*6+6 , 1)) - !end do - !print *, ' ************* ' - !print*, ra*go_surf, ra*go0, mu0 - !end if - !!! >>> SA + !print *, ra, go0, go_surf + + do j=1,nbc + if (j zabs(lovel1s/lovel)) layerrap = zabs(lovel1s/lovel) - if (layerrap > zabs((lovek1s-(rads/ra)**deg)/(lovek-1.d0))) & - layerrap = zabs((lovek1s-(rads/ra)**deg)/(lovek-1.d0)) - endif - - loveh1s = loveh1 - lovel1s = lovel1 - lovek1s = lovek1 - rads = radius(j) - !if (radius(j)==rb) then - ! write(*,*) 'CMB ',n,loveh1,lovel1,lovek1-(rb/ra)**deg - !else if (radius(j)==rc) then - ! write(*,*) 'ICB ',n,loveh1,lovel1,lovek1-(rc/ra)**deg - !else - ! if (j<10) then - ! write(*,*) 'ITF-',j,n,loveh1,lovel1,lovek1-(radius(j)/ra)**deg - ! else - ! write(*,*) 'ITF-',j,n,loveh1,lovel1,lovek1-(radius(j)/ra)**deg - ! end if - !endif - end if - end do - - end do - end if - + !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) Index: /issm/trunk-jpl/src/c/modules/FourierLoveCorex/model.f90 =================================================================== --- /issm/trunk-jpl/src/c/modules/FourierLoveCorex/model.f90 (revision 25294) +++ /issm/trunk-jpl/src/c/modules/FourierLoveCorex/model.f90 (revision 25295) @@ -65,4 +65,5 @@ 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