Changeset 25295

Timestamp:

07/20/20 23:24:15 (5 years ago)

Author:

adhikari

Message:

CHG: we extract Love kernels at planetary depth.

Location:

issm/trunk-jpl/src/c/modules/FourierLoveCorex

Files:

• love_numbers.f90 (modified) (15 diffs)
• lovenb_sub.f90 (modified) (11 diffs)
• model.f90 (modified) (1 diff)

issm/trunk-jpl/src/c/modules/FourierLoveCorex/love_numbers.f90

@@ -90 +90 @@
  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 +96 @@
  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 +108 @@
  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 +130 @@
  ! 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 +156 @@
 
  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 +191 @@
  !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 +216 @@
 
  bc(:,:)=0.d0
+
  
  
@@ -227 +232 @@
         endif
  deg=degmin-1
-
 
 
@@ -242 +246 @@
   !write(*,*)
   !!write(101,*), fr, T
-
-
-
 
  do fr=1,nfreq
@@ -289 +290 @@
    !print*,(deg)
    bc(:,:)=0.d0
+   lovekernels(:,1)=0.d0
    call boundary_conditions_matrix(bc,indx,nbc)
         !print*,'bc ok', layerrap
@@ -295 +297 @@
    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 +329 @@
      end do
      radbc(n) =.false.
+     radbc(n+1) =.true.
      !write(101,*) ' New start interface: ', radius(n+1)/1.d3,' km'
      if (display) then
@@ -323 +335 @@
      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 +349 @@
 
     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 +364 @@
     !!! 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))

issm/trunk-jpl/src/c/modules/FourierLoveCorex/lovenb_sub.f90

@@ -42 +42 @@
         
    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 +60 @@
   !    and typing the matrix boundary conditions 
   ibc = 0   ! Counter of interfaces with boundary conditions 
-
+ 
   do i = ifirst,nlayer 
   
@@ -81 +81 @@
       !!! >>> 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 +99 @@
                         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 +149 @@
   end do
 
-
   !-- Internal sphere
 
@@ -165 +164 @@
                   bc(nbc,nbc)=1.d0
         end if
-        
 
         !print*, bc, imag_eval
@@ -317 +315 @@
   logical :: ok
 
-
   bcsav(:,:) = bc(:,:)
   rads=0.d0
@@ -323 +320 @@
   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 +366 @@
         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 +377 @@
         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<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
@@ -426 +397 @@
 
 
-        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,1)*ra*go_surf*valini
-                                lovel1 = sc(nbc -ibc*6 -3 +3,1)*ra*go_surf*valini
-                                lovek1 = sc(nbc -ibc*6 -3 +5,1)*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
-                                        sumh = sumh+bcsav(ii+1,jj+ic)*sc(jj+ic,1)
-                                        suml = suml+bcsav(ii+3,jj+ic)*sc(jj+ic,1)
-                                        sumk = sumk+bcsav(ii+5,jj+ic)*sc(jj+ic,1)
-                                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 = zabs(loveh1s/loveh)
-                                if (layerrap > 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)

issm/trunk-jpl/src/c/modules/FourierLoveCorex/model.f90

@@ -65 +65 @@
  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