diff --git a/Changlog b/Changlog
index fdafadd93c2136a5ed2519df1ab84a5af5880d50..32cb69396862856987890e643074bb8de2ee25b6 100644
--- a/Changlog
+++ b/Changlog
@@ -98,3 +98,5 @@ tagging rel-1-4
31/10/21 hsrc is working for point sources. srcinv seems to work.
01/11/21 commit and tag srcinv01 before eavy changes (Nyquist time step etc)
03/11/21 commit and tag srcinv02 before eavy changes (multi layer mesh and good source description)
+05/05/22 commit changes done in April 2022 about psi source type (instead of f). Not sure it is finished
+ tagging rel-hsrc-1
diff --git a/nms/cutsrc.f90 b/nms/cutsrc.f90
index de4ce2f018e65ae9f9fe4ea1175adae3dfb41338..e5d36b6ebe5971aa47ca37dd40e76ef5447fd2a6 100644
--- a/nms/cutsrc.f90
+++ b/nms/cutsrc.f90
@@ -8,7 +8,7 @@ program cutsrc
character(len=100) :: name
type(meshinv) :: geom
real(DP), dimension(:), pointer :: coef
- integer :: ianswer,it,nn,idir,i
+ integer :: ianswer,it,nn,idir,i,irg
real(DP) :: r,theta,phi,dt,tt,ca,cb,rr,dr
real(DP),dimension(3) :: fsrc
@@ -33,7 +33,7 @@ program cutsrc
open(13,file='f_p.txt')
do it=1,nn
tt=(it-1)*dt
- call get_src_val(geom,coef,tt,r,theta,phi,fsrc,2)
+ call get_src_val(geom,coef,tt,irg,r,theta,phi,fsrc,2)
write(11,*) tt,fsrc(1)
write(12,*) tt,fsrc(2)
write(13,*) tt,fsrc(3)
@@ -48,13 +48,13 @@ program cutsrc
read*,tt,ca,cb
ca=ca*deg2rad
cb=cb*deg2rad
- dr=(geom%r2-geom%r1)/(nn-1)
+ dr=(geom%ring(geom%nring)%r2-geom%ring(1)%r1)/(nn-1)
open(11,file='f_r.txt')
open(12,file='f_t.txt')
open(13,file='f_p.txt')
do i=1,nn
- rr=geom%r1+(i-1)*dr
- call get_src_val(geom,coef,tt,rr,ca,cb,fsrc,2)
+ rr=geom%ring(1)%r1+(i-1)*dr
+ call get_src_val(geom,coef,tt,irg,rr,ca,cb,fsrc,2)
write(11,*) rr,fsrc(1)
write(12,*) rr,fsrc(2)
write(13,*) rr,fsrc(3)
@@ -70,7 +70,7 @@ program cutsrc
open(13,file='f_p.txt')
do i=1,nn
theta=pi/2+(i-1)*dt-geom%dtheta/2.d0
- call get_src_val(geom,coef,tt,ca,theta,cb,fsrc,2)
+ call get_src_val(geom,coef,tt,irg,ca,theta,cb,fsrc,2)
write(11,*) theta*rad2deg,fsrc(1)
write(12,*) theta*rad2deg,fsrc(2)
write(13,*) theta*rad2deg,fsrc(3)
@@ -86,7 +86,7 @@ program cutsrc
open(13,file='f_p.txt')
do i=1,nn
phi=(i-1)*dt-geom%dphi/2.d0
- call get_src_val(geom,coef,tt,ca,cb,phi,fsrc,2)
+ call get_src_val(geom,coef,tt,irg,ca,cb,phi,fsrc,2)
write(11,*) phi*rad2deg,fsrc(1)
write(12,*) phi*rad2deg,fsrc(2)
write(13,*) phi*rad2deg,fsrc(3)
diff --git a/nms/global_main.f90 b/nms/global_main.f90
index f7143da4150d202a47f22365002c82323fe222b6..47532c733194da4b8c0b718a1de3b6e76d24ef6a 100644
--- a/nms/global_main.f90
+++ b/nms/global_main.f90
@@ -92,13 +92,15 @@ contains
integer, dimension(:), optional, intent(in) :: is2rad_
real(DP), dimension(:), optional, intent(in) :: srcrad_
!
- integer :: i,ic,jj,NBPT
+ integer :: i,ic,jj,NBPT,NBC
+ real(DP) :: t,p
logunit=logunit_
!
if (present(srccoord_)) then
srcinv=.true.
NBPT=ubound(srccoord_,dim=2)
- NBE=3*NBPT
+ NBC=6 !6 moment tensor component possible.
+ NBE=NBC*NBPT
NBRS=NBRS_
allocate(srcrad(NBRS),is2rad(NBE))
srcrad(:)=srcrad_(:)
@@ -110,10 +112,14 @@ contains
half_duration(:)=0.d0
t_delay_sources(:)=0.d0
do i=1,NBPT
- do ic=1,3
- jj=ic+(i-1)*3
+ do ic=1,NBC
+ jj=ic+(i-1)*NBC
is2rad(jj)=is2rad_(i)
+!rotate xyz to rtp here?
Mtmp(ic,jj)=1.d0
+ t=srccoord_(2,i)*deg2rad
+ p=srccoord_(3,i)*deg2rad
+ call convert_M_to_spherical(Mtmp(:,jj),t,p)
coord_sources(:,jj)=srccoord_(:,i)
write(sources_name(jj),'("G",i5.5,"_",i1)') i,ic
enddo
@@ -196,7 +202,7 @@ contains
write(logunit,'(a)') eigenfileS
write(logunit,'(a)') eigenfileT
if (srcinv) then
- source_force=.true.
+ source_force=.false. !using moment tensor now
channel=VEL
amp_source=1.d0
else
@@ -673,6 +679,84 @@ contains
close(11)
enddo
end subroutine load_station_response
+!-------------------------------------
+ subroutine convert_M_to_spherical(M,theta,phi)
+ implicit none
+ doubleprecision,dimension(:), intent(inout) :: M
+ doubleprecision, intent(in) :: theta,phi
+ doubleprecision, dimension(3,3) :: Mt
+ !
+ Mt(1,1)=M(1)
+ Mt(2,2)=M(2)
+ Mt(3,3)=M(3)
+ Mt(1,2)=M(4)
+ Mt(2,1)=M(4)
+ Mt(1,3)=M(5)
+ Mt(3,1)=M(5)
+ Mt(2,3)=M(6)
+ Mt(3,2)=M(6)
+ call c_tensor2(Mt,theta,phi)
+ M(1)=Mt(1,1)
+ M(2)=Mt(2,2)
+ M(3)=Mt(3,3)
+ M(4)=Mt(1,2)
+ M(5)=Mt(1,3)
+ M(6)=Mt(2,3)
+ end subroutine convert_M_to_spherical
+!-------------------------------------
+ subroutine c_tensor2(M,theta,phi)
+! cartesian to spherical conversion
+!-------------------------------------
+ implicit none
+ doubleprecision,dimension(:,:), intent(inout) :: M
+ doubleprecision, intent(in) :: theta,phi
+!
+ doubleprecision,dimension(3,3) :: CCS,CSC,T
+ doubleprecision :: ct,st,cp,sp
+ integer :: i,j,k
+!
+ ct=dcos(theta)
+ st=dsin(theta)
+ cp=dcos(phi )
+ sp=dsin(phi )
+!
+! passage de catersien a spherique
+!
+ Ccs(1,1)= st*cp; Ccs(1,2)=st*sp; Ccs(1,3)= ct
+ Ccs(2,1)= ct*cp; Ccs(2,2)=ct*sp; Ccs(2,3)= -st
+ Ccs(3,1)=-sp ; Ccs(3,2)=cp ; Ccs(3,3)= 0.0d0
+!
+! passage de spherique a cartesien
+!
+ Csc(1,1)= st*cp; Csc(1,2)= ct*cp; Csc(1,3)= -sp
+ Csc(2,1)= st*sp; Csc(2,2)= ct*sp; Csc(2,3)= cp
+ Csc(3,1)= ct ; Csc(3,2)=-st ; Csc(3,3)= 0.0d0
+!
+! M*Ccs
+!
+ do i=1,3
+ do j=1,3
+ T(i,j)=0.0d0
+ do k=1,3
+ T(i,j)=T(i,j)+M(i,k)*Csc(k,j)
+ enddo
+ enddo
+ enddo
+!
+! Csc*M*Csc
+!
+ do i=1,3
+ do j=1,3
+ M(i,j)=0.0d0
+ do k=1,3
+ M(i,j)=M(i,j)+Ccs(i,k)*T(k,j) !tR.M.R
+ enddo
+ enddo
+ enddo
+!-------------------------------------
+ end subroutine c_tensor2
+!-------------------------------------
+
!------------------------------------------------------
end module global_main_param
!------------------------------------------------------
diff --git a/nms/module_modele.f90 b/nms/module_modele.f90
index c25500463989ef39e8d3f6998e0905f8ff9cb7b1..be6d01ffe67d37d861a6dbd8f486e2ef92d383cb 100644
--- a/nms/module_modele.f90
+++ b/nms/module_modele.f90
@@ -2,6 +2,7 @@
module earth_modele
!--------------------------------------------------------------------------
implicit none
+ private:: locate
!
!----------------------------------------------------------
type mod_der
@@ -33,16 +34,17 @@ module earth_modele
! eta : F/(A-2L)
!----------------------------
character(len=80) :: titre
- integer :: nbcou
+ integer :: nbcou,nbi
integer :: ifanis,ifdeck,nic,noc,nbceau
real :: tref
type(mod_der) :: der
- doubleprecision, dimension(:), pointer :: r,rho,vpv,vsv,qkappa,qshear,vph,vsh,eta
+ doubleprecision, dimension(:), pointer :: r,rho,vpv,vsv,qkappa,qshear,vph,vsh,eta,radi
+ integer, dimension(:), pointer :: inti
logical :: allocated
!----------------------------
end type modele
!----------------------------
-
+
!--------------------------------------------------------------------------
contains
!--------------------------------------------------------------------------
@@ -121,7 +123,7 @@ contains
character(*), intent(in) :: name
type(modele), intent(out) :: mod
!
- integer :: i,unit=127,icode
+ integer :: i,unit=127,icode,nbi
!
call init_modele(mod)
open(unit,file=name,iostat=icode,status='old',action='read')
@@ -148,6 +150,28 @@ contains
stop 'read_modele: prob de de lecture'
endif
close(unit)
+ nbi=0
+ do i=1,mod%nbcou-1
+ if ( (mod%r(i+1)-mod%r(i))/mod%r(i+1)<1.d-7 ) then
+ nbi=nbi+1
+ endif
+ enddo
+ mod%nbi=nbi
+ allocate(mod%radi(0:nbi+1),mod%inti(0:nbi+1))
+ mod%inti(0) =1
+ mod%radi(0) =0.
+ mod%inti(nbi+1)=mod%nbcou
+ mod%radi(nbi+1)=mod%r(mod%nbcou)
+!
+ nbi=0
+ do i=1,mod%nbcou-1
+ if ( (mod%r(i+1)-mod%r(i))/mod%r(i+1)<1.d-7 ) then
+ nbi=nbi+1
+ mod%radi(nbi)=mod%r(i)
+ mod%inti(nbi)=i
+ endif
+ enddo
+
!
100 format(a80)
105 format(f8.0, 3f9.2, 2f9.1, 2f9.2, f9.5)
@@ -269,6 +293,107 @@ contains
!--------------------------------------------------------------------------
end function modele_different
!--------------------------------------------------------------------------
+!--------------------------------------------------------------------------
+ subroutine reinterpa_interface(A,NBA,RA,nbi,inti,ri,is,B,NBB,RB,flag)
+!flag: zero deriv:
+!--------------------------------------------------------------------------
+ use def_gparam
+ use module_spline
+ implicit none
+ integer, intent(in) :: NBA,NBB,nbi,is
+ logical, intent(in) :: flag
+ integer, dimension(0:), intent(in) :: inti
+ real(DP), dimension(0:), intent(in) :: ri
+ real(DP), dimension(:), intent(in) :: A,RA,RB
+ real(DP), dimension(:), intent(out) :: B
+!
+ integer :: i,nban,nd,ka,kb
+ real(DP), dimension(NBA) :: y2
+ real(DP) :: yp1,ypn,rad
+
+ do nd=1,nbi+1
+ ka=inti(nd-1)+1
+ kb=inti(nd)
+ if (nd==nbi+1.and.flag) then
+ ypn=0.d0
+ else
+ ypn=(A(kb )-A(kb-1))/(RA(kb )-RA(kb-1))
+ endif
+ if (nd==1.and.flag) then
+ yp1=0.d0
+ else
+ yp1=(A(ka+1)-A(ka ))/(RA(ka+1)-RA(ka ))
+ endif
+ call spline(RA(ka:kb),A(ka:kb),yp1,ypn,y2(ka:kb))
+ enddo
+ do i=1,NBB
+ if (RB(i)<RA(is)) then
+ B(i)=A(is)
+ else
+ rad=RB(i)
+ nd=locate(ri,nbi,rad)
+ ka=inti(nd-1)+1
+ kb=inti(nd)
+ B(i)=splint(RA(ka:kb),A(ka:kb),y2(ka:kb),RB(i))
+ endif
+ enddo
+!--------------------------------------------------------------------------
+ end subroutine reinterpa_interface
+!--------------------------------------------------------------------------
+!-----------------------------------------------------------------
+ subroutine get_1dmodel_value(mod,rho,Vsh,Vsv,Vph,Vpv,eta,Qs,NBV,rad)
+!-----------------------------------------------------------------
+ use def_gparam
+ use module_spline
+ implicit none
+ integer, intent(in) :: NBV
+ type(modele), intent(in) :: mod
+ real(DP), dimension(:), intent(in ) :: rad
+ real(DP), dimension(:), intent(out) :: rho,Vsh,Vph,Vsv,Vpv,Qs,eta
+!
+ call reinterpa_interface(mod%rho ,mod%nbcou,mod%r,mod%nbi,mod%inti,mod%radi,1,rho,NBV,rad,.true.)
+ call reinterpa_interface(mod%vsh ,mod%nbcou,mod%r,mod%nbi,mod%inti,mod%radi,1,vsh,NBV,rad,.true.)
+ call reinterpa_interface(mod%vsv ,mod%nbcou,mod%r,mod%nbi,mod%inti,mod%radi,1,vsv,NBV,rad,.true.)
+ call reinterpa_interface(mod%vph ,mod%nbcou,mod%r,mod%nbi,mod%inti,mod%radi,1,vph,NBV,rad,.true.)
+ call reinterpa_interface(mod%vpv ,mod%nbcou,mod%r,mod%nbi,mod%inti,mod%radi,1,vpv,NBV,rad,.true.)
+ call reinterpa_interface(mod%eta ,mod%nbcou,mod%r,mod%nbi,mod%inti,mod%radi,1,eta,NBV,rad,.true.)
+ call reinterpa_interface(mod%qshear,mod%nbcou,mod%r,mod%nbi,mod%inti,mod%radi,1,qs ,NBV,rad,.true.)
+
+
+!--------------------------------------------------------------------------
+ end subroutine get_1dmodel_value
+!--------------------------------------------------------------------------
+
+!-----------------------------------------------------------------
+ integer function locate(inter,n,x)
+!-----------------------------------------------------------------
+ use def_gparam
+ implicit none
+ integer, intent(in) :: n
+ real(DP), dimension(0:n+1), intent(in) :: inter
+ real(DP), intent(in) :: x
+!
+ integer :: i
+ i=0
+ if (abs (x-inter(0))/max(x,inter(0)) <1.E-10_DP) then
+ i=1
+ else if ( abs (x-inter(n+1))/max(x,inter(n+1)) <1.E-10_DP) then
+ i=n
+ else
+ do while ( x < inter(i) .or. x > inter(i+1) )
+ i=i+1
+ if (i > n) then
+ print*,'i=',i,'n=',n
+ print*,'x=',x
+ print*,'inter=',inter
+ stop 'locate: failed to locate x in inter!'
+ endif
+ enddo
+ endif
+ locate=i+1
+!-----------------------------------------------------------------
+ end function locate
+!-----------------------------------------------------------------
!--------------------------------------------------------------------------
end module earth_modele
diff --git a/nms/module_modes.f90 b/nms/module_modes.f90
index f67e147f551e5d458324e0c9efbb8ea4284972b5..8d5c21e346cf60eacda0dc855059836a519c481f 100644
--- a/nms/module_modes.f90
+++ b/nms/module_modes.f90
@@ -169,10 +169,10 @@ contains
real(NBRS)*((Nmax_all+1)*(LmaxS+1)*4+ &
(maxval(NmaxT(:))+1)*(LmaxT+1)*2)*8/1024./1024. &
,' Mb/proc)'
- allocate(UUs(NBE,0:Nmax_allS,0:LmaxS),UUps(NBE,0:Nmax_allS,0:LmaxS) &
- ,VVs(NBE,0:Nmax_allS,0:LmaxS),VVps(NBE,0:Nmax_allS,0:LmaxS) )
+ allocate(UUs(NBRS,0:Nmax_allS,0:LmaxS),UUps(NBRS,0:Nmax_allS,0:LmaxS) &
+ ,VVs(NBRS,0:Nmax_allS,0:LmaxS),VVps(NBRS,0:Nmax_allS,0:LmaxS) )
allocate(WW(1,0:Nmax_allT,0:LmaxT),WWp(1,0:Nmax_allT,0:LmaxT))
- allocate(WWs(NBE,0:Nmax_allT,0:LmaxT),WWps(NBE,0:Nmax_allT,0:LmaxT))
+ allocate(WWs(NBRS,0:Nmax_allT,0:LmaxT),WWps(NBRS,0:Nmax_allT,0:LmaxT))
UU=0.0_DP;VV=0.0_DP;UUp=0.0_DP;VVp=0.0_DP;WW=0.0_DP;WWp=0.0_DP
!
print*,'reading and interpolating eigenfunction ...'
@@ -1408,6 +1408,7 @@ contains
!!$ endif
irr=ir2rad(ir)
irs=is2rad(is)
+! print*,'coucou',ir,irr,is,irs
tr=rec_coord(2,ir)
pr=rec_coord(3,ir)
ts=src_coord(2,is)
diff --git a/nms/module_srcinv.f90 b/nms/module_srcinv.f90
index 7f195ebb8a446304905baf495b35cdcbdf7aa771..67d63d137f58758abba1b765c39c2fc1e56da6ee 100644
--- a/nms/module_srcinv.f90
+++ b/nms/module_srcinv.f90
@@ -5,28 +5,48 @@ module module_srcinv
!-----------------------------------------------------------------------------
use def_gparam
use time_function
+ use earth_modele
implicit none
private
public:: meshinv, generate_invmesh,write_mesh2geo,xyz2rtp,init_srcinv,build_green_fcts,write_src,read_src &
,sum_green,read_data,build_ATCA,invert_ATCA,get_pt_time_function,comp_misfit,get_src_val,rotate2center,sum_green_fcts
!-------------------------
- type meshinv
- real(DP) :: phi_center,theta_center !theta and phi
- real(DP) :: dtheta,dphi,r1,r2,t0,time_support,fmax,fmaxb,tmin,tmax,dtinv,dttrace
- real(DP) :: dtdata,t0data,damp_coef,timetapper,dtout
- integer :: ngllh_shape,ngllh_int,ngllv_shape,ngllv_int,NBTinv,NBP,NBPs,NBTtrace,NBTdata,NBR,nbcomp,NBTout
+ type mesh_ring
+ real(DP) :: r1,r2
+ integer :: ngllh_shape,ngllh_int,ngllv_shape,ngllv_int,NBPs,NBPsc,NBPi,NBPic,NBradic
real(DP), dimension(:,:,:,:), pointer :: mesh_shape,mesh_int
+ real(DP), dimension(:,:,:,:,:), pointer :: deriv
real(DP), dimension(:,:,:,:,:,:), pointer :: shape_fcts
- real(DP), dimension(:,:) , pointer :: mesh_int_rtp
+ real(DP), dimension(:,:,:,:,:,:,:), pointer :: dshape_fcts
real(DP), dimension(:,:,:), pointer :: w_int
+!NPPs: number of shape paramter in one ring
+!NBPsc: NumBer of shape Paramter cumulated up to ring i-1
+!NPPi: number of integration points in one ring
+!NBPic: NumBer of integration points cumulated up to ring i-1
+!NBradic: cumulative number of integration points in the previous rings
+ end type mesh_ring
+ type meshinv
+ real(DP) :: phi_center,theta_center !theta and phi
+ real(DP) :: dtheta,dphi,t0,time_support,fmax,fmaxb,tmin,tmax,dtinv,dttrace
+ real(DP) :: t0data,damp_coef,timetapper,dtout
+ integer :: NBCM,NBTinv,NBP,NBPs,NBPsc,NBTtrace,NBR,nbcomp,NBTout,nring,nbrad_int,NBPint
+!NBCM: source number of components (in normal mode. For psi inv, =6)
+!NBPint: total number if integration points
+!NBPsc : total number of shape functions
+!nbrad_int: total number of different radius for integration points
+ real(DP), dimension(:,:) , pointer :: mesh_int_rtp
+ real(DP), dimension(:,:,:) , pointer :: CIJ
real(DP), dimension(:,:,:,:), pointer :: green
real(DP), dimension(:,:), pointer :: ATCA
real(DP), dimension(:), pointer :: mask,damping_factor,mesh_rad
integer, dimension(:), pointer :: iglob2irad
- character(len=100) :: datadir
+ type(mesh_ring), dimension(:), pointer :: ring
+ character(len=100) :: datadir,modelename
type(timefunction) :: base
logical :: applymask,save_green,time_side_damping
+ type(modele) :: modele1d
end type meshinv
+
!------------------------
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
contains
@@ -37,18 +57,43 @@ contains
use def_gparam
implicit none
type(meshinv), intent(inout) :: geom
+ integer :: i
!
open(11,file='srcinv.dat')
read(11,*) ! center of the inversion boxe (phi,theta) in degree
read(11,*) geom%phi_center,geom%theta_center
read(11,*) ! box extents (dphi,dtheta) in degree
read(11,*) geom%dphi,geom%dtheta
- read(11,*) ! box inversion radius bounds r1,r2 in m
- read(11,*) geom%r1,geom%r2
- read(11,*) ! inversion degrees in horizontal and vertical directions
- read(11,*) geom%ngllh_shape,geom%ngllv_shape
- read(11,*) ! integration de degrees in horizontal and vertical directions
- read(11,*) geom%ngllh_int,geom%ngllv_int
+ read(11,*) !enter the number of mesh rings
+ read(11,*) geom%nring
+ allocate(geom%ring(geom%nring))
+ read(11,*) !enter the nring+1 radius
+ do i=1,geom%nring
+ read(11,*) geom%ring(i)%r1
+ if (i>1) geom%ring(i-1)%r2=geom%ring(i)%r1
+ enddo
+ read(11,*) geom%ring(geom%nring)%r2
+ read(11,*) ! for each ring: inversion degrees in horizontal and vertical directions
+ read(11,*) ! then: integration de degrees in horizontal and vertical directions
+ do i=1,geom%nring
+ if (i==1) then
+ geom%ring(i)%NBPsc=0
+ geom%ring(i)%NBPic=0
+ geom%ring(i)%NBradic=0
+ else
+ geom%ring(i)%NBPsc=geom%ring(i-1)%NBPsc+geom%ring(i-1)%NBPs
+ geom%ring(i)%NBPic=geom%ring(i-1)%NBPic+geom%ring(i-1)%NBPi
+ geom%ring(i)%NBradic=geom%ring(i-1)%NBradic+geom%ring(i-1)%ngllv_int !cumulative number of integration points in the previous rings
+ endif
+ read(11,*) geom%ring(i)%ngllh_shape,geom%ring(i)%ngllv_shape,geom%ring(i)%ngllh_int,geom%ring(i)%ngllv_int
+ geom%ring(i)%NBPs=geom%ring(i)%ngllh_shape**2*geom%ring(i)%ngllv_shape
+ geom%ring(i)%NBPi=geom%ring(i)%ngllh_int**2 *geom%ring(i)%ngllv_int
+ geom%nbrad_int=geom%nbrad_int+geom%ring(i)%ngllv_int
+ enddo
+ geom%NBPint=geom%ring(geom%nring)%NBPic+geom%ring(geom%nring)%NBPi !total number of integration points
+ geom%NBPsc =geom%ring(geom%nring)%NBPsc+geom%ring(geom%nring)%NBPs !total number of shape points
+ read(11,*) !1D modele file name
+ read(11,'(a)') geom%modelename
read(11,*) ! t0, time duration and maximum frequency
read(11,*) geom%t0,geom%time_support,geom%fmax
read(11,*) !data directory
@@ -83,7 +128,7 @@ contains
geom%dphi =geom%dphi *deg2rad
geom%dtheta =geom%dtheta *deg2rad
! number of free parameters:
- geom%NBPs=geom%ngllh_shape**2*geom%ngllv_shape*3 !number of parameter per time slice
+ geom%NBPs=geom%NBPsc*3 !number of parameter per time slice
geom%NBP=geom%NBTinv*geom%NBPs
print*,'Number of space parameters x components:',geom%NBPs
print*,'Number of time parameters per space parameter:',geom%NBTinv
@@ -93,73 +138,84 @@ contains
close(11)
allocate(geom%damping_factor( geom%NBP))
geom%damping_factor(:)=1._DP
+ geom%NBCM=6
call get_shape_fcts(geom)
+ call read_modele(geom%modelename,geom%modele1d)
!----------------------------------------------------------------------------
end subroutine init_srcinv
+!----------------------------------------------------------------------------
!----------------------------------------------------------------------------
subroutine get_shape_fcts(geom)
+!----------------------------------------------------------------------------
use def_gparam
use funaro
implicit none
type(meshinv), intent(inout) :: geom
!
- integer :: ii,jj,kk,i,j,k,N,Ni
+ integer :: ii,jj,kk,i,j,k,N,Ni,irg
real(DP), dimension(:), allocatable :: ET,VN,Q,ETi,VNi
- real(DP), dimension(geom%ngllh_int,geom%ngllh_shape) :: hx
- real(DP), dimension(geom%ngllv_int,geom%ngllv_shape) :: hv
- if (geom%ngllv_shape==1.and.geom%ngllh_shape==1.and.geom%ngllv_int==1.and.geom%ngllh_int==1) then
- allocate(geom%shape_fcts(geom%ngllh_int,geom%ngllh_int,geom%ngllv_int &
- ,geom%ngllh_shape,geom%ngllh_shape,geom%ngllv_shape))
- geom%shape_fcts=1._DP
+ real(DP), dimension(:,:), allocatable :: hx,hv,dhx,dhv
+ do irg=1,geom%nring
+ if (geom%ring(irg)%ngllv_shape==1.and.geom%ring(irg)%ngllh_shape==1.and.geom%ring(irg)%ngllv_int==1.and.geom%ring(irg)%ngllh_int==1) then
+ allocate(geom%ring(irg)%shape_fcts(geom%ring(irg)%ngllh_int,geom%ring(irg)%ngllh_int,geom%ring(irg)%ngllv_int &
+ ,geom%ring(irg)%ngllh_shape,geom%ring(irg)%ngllh_shape,geom%ring(irg)%ngllv_shape))
+ geom%ring(irg)%shape_fcts=1._DP
else
+ allocate (hx(geom%ring(irg)%ngllh_int,geom%ring(irg)%ngllh_shape),hv(geom%ring(irg)%ngllv_int,geom%ring(irg)%ngllv_shape) )
+ allocate (dhx(geom%ring(irg)%ngllh_int,geom%ring(irg)%ngllh_shape),dhv(geom%ring(irg)%ngllv_int,geom%ring(irg)%ngllv_shape) )
! x and y:
- N=geom%ngllh_shape-1
+ N=geom%ring(irg)%ngllh_shape-1
allocate(ET(0:N),VN(0:N),Q(0:N))
call ZELEGL(N,ET,VN)
- if (geom%ngllh_int>0) then
- Ni=geom%ngllh_int-1
+ if (geom%ring(irg)%ngllh_int>0) then
+ Ni=geom%ring(irg)%ngllh_int-1
allocate(ETi(0:Ni),VNi(0:Ni))
call ZELEGL(Ni,ETi,VNi)
- do ii=1,geom%ngllh_shape
+ do ii=1,geom%ring(irg)%ngllh_shape
Q(:)=0.d0
Q(ii-1)=1.d0
- do i=1,geom%ngllh_int
- call INLEGL(N,ET,VN,Q,ETi(i-1),hx(i,ii) )
+ do i=1,geom%ring(irg)%ngllh_int
+ call INLEGL (N,ET,VN,Q,ETi(i-1), hx(i,ii) )
+ call INLEGL2(N,ET,VN,Q,ETi(i-1),dhx(i,ii) )
enddo
enddo
deallocate(ETi,VNi)
endif
deallocate(ET,VN,Q)
! z:
- N=geom%ngllv_shape-1
+ N=geom%ring(irg)%ngllv_shape-1
allocate(ET(0:N),VN(0:N),Q(0:N))
call ZELEGL(N,ET,VN)
- Ni=geom%ngllv_int-1
- if (geom%ngllv_int>0) then
+ Ni=geom%ring(irg)%ngllv_int-1
+ if (geom%ring(irg)%ngllv_int>0) then
allocate(ETi(0:Ni),VNi(0:Ni))
call ZELEGL(Ni,ETi,VNi)
- do ii=1,geom%ngllv_shape
+ do ii=1,geom%ring(irg)%ngllv_shape
Q(:)=0.d0
Q(ii-1)=1.d0
- do i=1,geom%ngllv_int
- call INLEGL(N,ET,VN,Q,ETi(i-1),hv(i,ii) )
+ do i=1,geom%ring(irg)%ngllv_int
+ call INLEGL (N,ET,VN,Q,ETi(i-1),hv (i,ii) )
+ call INLEGL2(N,ET,VN,Q,ETi(i-1),dhv(i,ii) )
enddo
enddo
deallocate(ETi,VNi)
endif
deallocate(ET,VN,Q)
!
- if (geom%ngllh_int>1) then
- allocate(geom%shape_fcts(geom%ngllh_int,geom%ngllh_int,geom%ngllv_int &
- ,geom%ngllh_shape,geom%ngllh_shape,geom%ngllv_shape))
+ if (geom%ring(irg)%ngllh_int>1) then
+ allocate(geom%ring(irg)%shape_fcts(geom%ring(irg)%ngllh_int,geom%ring(irg)%ngllh_int,geom%ring(irg)%ngllv_int &
+ ,geom%ring(irg)%ngllh_shape,geom%ring(irg)%ngllh_shape,geom%ring(irg)%ngllv_shape))
!
- do kk=1,geom%ngllv_shape
- do jj=1,geom%ngllh_shape
- do ii=1,geom%ngllh_shape
- do k=1,geom%ngllv_int
- do j=1,geom%ngllh_int
- do i=1,geom%ngllh_int
- geom%shape_fcts(i,j,k,ii,jj,kk)=hx(i,ii)*hx(j,jj)*hv(k,kk)
+ do kk=1,geom%ring(irg)%ngllv_shape
+ do jj=1,geom%ring(irg)%ngllh_shape
+ do ii=1,geom%ring(irg)%ngllh_shape
+ do k=1,geom%ring(irg)%ngllv_int
+ do j=1,geom%ring(irg)%ngllh_int
+ do i=1,geom%ring(irg)%ngllh_int
+ geom%ring(irg)%shape_fcts(i,j,k,ii,jj,kk)=hx(i,ii)*hx(j,jj)*hv(k,kk)
+ geom%ring(irg)%dshape_fcts(1,i,j,k,ii,jj,kk)=dhx(i,ii)* hx(j,jj)* hv(k,kk)
+ geom%ring(irg)%dshape_fcts(2,i,j,k,ii,jj,kk)= hx(i,ii)*dhx(j,jj)* hv(k,kk)
+ geom%ring(irg)%dshape_fcts(3,i,j,k,ii,jj,kk)= hx(i,ii)* hx(j,jj)*dhv(k,kk)
enddo
enddo
enddo
@@ -168,23 +224,29 @@ contains
enddo
endif
endif
+ deallocate(hv,hx)
+ enddo
!----------------------------------------------------------------------------
end subroutine get_shape_fcts
!----------------------------------------------------------------------------
- integer function ipnumber(geom,ic,i,j,k,type)
+ integer function ipnumber(geom,ic,irg,i,j,k,type)
!type: 1: int, 2 : shape
!----------------------------------------------------------------------------
implicit none
type(meshinv), intent(in) :: geom
- integer, intent(in) :: ic,i,j,k,type
+ integer, intent(in) :: ic,irg,i,j,k,type
+ !
+ integer :: ngllh,ngllv,istart,nbc
!
select case(type)
case (1)
+!integration points
! ipnumber=ic+(k-1)*3+(j-1)*3*geom%ngllh_int+(i-1)*3*geom%ngllh_int**2
- ipnumber=ic+(i-1)*3+(j-1)*3*geom%ngllh_int+(k-1)*3*geom%ngllh_int**2
+ ipnumber=ic+((i-1)+(j-1)*geom%ring(irg)%ngllh_int+(k-1)*geom%ring(irg)%ngllh_int**2+geom%ring(irg)%NBPic)*geom%NBCM
case (2)
+!shape functions
! ipnumber=ic+(k-1)*3+(j-1)*3*geom%ngllh_shape+(i-1)*3*geom%ngllh_shape**2
- ipnumber=ic+(i-1)*3+(j-1)*3*geom%ngllh_shape+(k-1)*3*geom%ngllh_shape**2
+ ipnumber=ic+((i-1)+(j-1)*geom%ring(irg)%ngllh_shape+(k-1)*geom%ring(irg)%ngllh_shape**2+geom%ring(irg)%NBPsc)*3
case default
STOP 'upnumber: type can only be 1 or 2'
end select
@@ -192,12 +254,12 @@ contains
end function ipnumber
!----------------------------------------------------------------------------
!-------------------------------------------------
- subroutine get_pt_time_function(geom,coef,ic,ix,iy,iz,tf,cc_)
+ subroutine get_pt_time_function(geom,coef,ic,irg,ix,iy,iz,tf,cc_)
!-------------------------------------------------
implicit none
type(meshinv), intent(inout) :: geom
real(DP), dimension(:), intent(in) :: coef
- integer, intent(in) :: ic,ix,iy,iz
+ integer, intent(in) :: ic,ix,iy,iz,irg
real(DP), dimension(:), pointer, intent(out) :: tf
real(DP), dimension(:), optional, intent(out) :: cc_
!
@@ -207,7 +269,7 @@ contains
!
allocate(tf(geom%NBTtrace))
do i=1,geom%NBTinv
- ii=ipnumber(geom,ic,ix,iy,iz,2)+(i-1)*geom%NBPs
+ ii=ipnumber(geom,ic,irg,ix,iy,iz,2)+(i-1)*geom%NBPs
cc(i)=coef(ii)
enddo
if (present(cc_)) cc_=cc
@@ -234,14 +296,19 @@ contains
type(meshinv), intent(inout) :: geom
!
real (DP), dimension(:,:,:,:), pointer :: green
- integer :: logunit,i,j,k,ii,jj,kk,ic,ip,ipint,it,ir,ips,icp,idelta
+ integer :: logunit,i,j,k,ii,jj,kk,ic,ip,ipint,it,ir,ips,icp,idelta,irg,a
real(DP) :: tshift,t,ta,tb,coef,time_damp
+ real(DP), dimension(:), allocatable :: Mc
character(len=30) :: name
logunit=1327
!
print*,'Preparing green functions ...'
open(logunit,file='srcinv.log')
- call init_global_main(logunit,geom%mesh_int_rtp,geom%ngllv_int,geom%mesh_rad,geom%iglob2irad)
+! do i=1,ubound(geom%iglob2irad,dim=1)
+! print*,i,geom%iglob2irad(i)
+! enddo
+! stop 'test'
+ call init_global_main(logunit,geom%mesh_int_rtp,geom%nbrad_int,geom%mesh_rad,geom%iglob2irad)
geom%NBTout=NBT
geom%dtout=dt
!over reading NBT and dt from nms.dat with Nyquist dt
@@ -250,7 +317,7 @@ contains
NBT=int(2.d0**(int(log(dble(NBT))/log(2.d0))+1))
dt=(geom%NBTout-1)*geom%dtout/(NBT-1)
! dt=dt*1.5
-! NBT=int((geom%NBTout-1)*geom%dtout/dt)+1
+ ! NBT=int((geom%NBTout-1)*geom%dtout/dt)+1
print*,'Actual dt and NBT used:',dt,NBT
geom%NBTtrace=NBT
geom%dttrace=dt
@@ -282,20 +349,27 @@ contains
allocate(geom%green(geom%NBTtrace,3,geom%NBR,geom%NBP))
geom%green(:,:,:,:)=0.d0
ip=0
-!$OMP PARALLEL DO COLLAPSE(4) DEFAULT(SHARED) PRIVATE(kk,jj,ii,ic,ip,ir,icp,k,j,i,ipint)
- do kk=1,geom%ngllv_shape
- do jj=1,geom%ngllh_shape
- do ii=1,geom%ngllh_shape
- do ic=1,3
- ip=ipnumber(geom,ic,ii,jj,kk,2)
- do ir=1,geom%NBR
-! do icp=1,3
- do k=1,geom%ngllv_int
- do j=1,geom%ngllh_int
- do i=1,geom%ngllh_int
- ipint=ipnumber(geom,ic,i,j,k,1)
- geom%green(:,:,ir,ip)=geom%green(:,:,ir,ip) &
- +green(:,:,ir,ipint)*geom%w_int(i,j,k)*geom%shape_fcts(i,j,k,ii,jj,kk)
+ do irg=1,geom%nring
+!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(kk,jj,ii,ic,ip,ir,icp,k,j,i,ipint,Mc)
+ allocate(Mc(6*geom%ring(irg)%ngllh_int*geom%ring(irg)%ngllh_int*geom%ring(irg)%ngllv_int))
+!$OMP DO COLLAPSE(4)
+ do kk=1,geom%ring(irg)%ngllv_shape
+ do jj=1,geom%ring(irg)%ngllh_shape
+ do ii=1,geom%ring(irg)%ngllh_shape
+ do ic=1,3
+ ip=ipnumber(geom,ic,irg,ii,jj,kk,2)
+ call get_Mc(geom,ic,irg,ii,jj,kk,Mc)
+! do ir=1,geom%NBR
+ do k=1,geom%ring(irg)%ngllv_int
+ do j=1,geom%ring(irg)%ngllh_int
+ do i=1,geom%ring(irg)%ngllh_int
+ do a=1,geom%NBCM
+ ipint=ipnumber(geom,a,irg,i,j,k,1)
+ geom%green(:,:,:,ip)=geom%green(:,:,:,ip) &
+ +green(:,:,:,ipint)*Mc(ipint)
+ enddo
+! geom%green(:,:,ir,ip)=geom%green(:,:,ir,ip) &
+! +green(:,:,ir,ipint)*geom%ring(irg)%w_int(i,j,k)*geom%ring(irg)%shape_fcts(i,j,k,ii,jj,kk)
enddo
enddo
enddo
@@ -304,8 +378,10 @@ contains
enddo
enddo
enddo
+!$OMP END DO
+ deallocate(Mc)
+!$OMP END PARALLEL
enddo
-!$OMP END PARALLEL DO
!$OMP PARALLEL DO COLLAPSE(2) DEFAULT(SHARED) PRIVATE(it,ips,ip,tshift,ic,ir)
do it=2,geom%NBTinv
@@ -390,14 +466,14 @@ contains
real (DP), dimension(:,:,:,:), pointer :: green,green_shift
real(DP), dimension(:,:,:), allocatable :: traces
real(DP), dimension(:), allocatable :: tmp
- integer :: logunit,i,j,k,ii,jj,kk,ic,ip,ipint,it,ir,ips,icp,idelta
+ integer :: logunit,i,j,k,ii,jj,kk,ic,ip,ipint,it,ir,ips,icp,idelta,irg
real(DP) :: tshift,t,ta,tb,coef,time_damp
character(len=30) :: name
logunit=1327
!
print*,'Preparing green functions ...'
open(logunit,file='srcinv.log')
- call init_global_main(logunit,geom%mesh_int_rtp,geom%ngllv_int,geom%mesh_rad,geom%iglob2irad)
+ call init_global_main(logunit,geom%mesh_int_rtp,geom%nbrad_int,geom%mesh_rad,geom%iglob2irad)
geom%NBTout=NBT
geom%dtout=dt
!over reading NBT and dt from nms.dat with Nyquist dt
@@ -439,49 +515,50 @@ contains
print*,'summing (1) ...'
allocate(geom%green(geom%NBTtrace,3,geom%NBR,geom%NBPs) &
,green_shift(geom%NBTtrace,3,geom%NBR,geom%NBPs) )
- if (geom%ngllv_shape==geom%ngllv_int.and.geom%ngllh_shape==geom%ngllh_int) then
+ geom%green(:,:,:,:)=0.d0
+ do irg=1,geom%nring
+ if (geom%ring(irg)%ngllv_shape==geom%ring(irg)%ngllv_int.and.geom%ring(irg)%ngllh_shape==geom%ring(irg)%ngllh_int) then
!$OMP PARALLEL DO COLLAPSE(5) DEFAULT(SHARED) PRIVATE(ir,ic,k,j,i,ip)
- do ir=1,geom%NBR
- do ic=1,3
- do k=1,geom%ngllv_int
- do j=1,geom%ngllh_int
- do i=1,geom%ngllh_int
- ip=ipnumber(geom,ic,i,j,k,1)
- geom%green(:,:,ir,ip)=green(:,:,ir,ip)*geom%w_int(i,j,k)
+ do ir=1,geom%NBR
+ do ic=1,3
+ do k=1,geom%ring(irg)%ngllv_int
+ do j=1,geom%ring(irg)%ngllh_int
+ do i=1,geom%ring(irg)%ngllh_int
+ ip=ipnumber(geom,ic,irg,i,j,k,1)
+ geom%green(:,:,ir,ip)=green(:,:,ir,ip)*geom%ring(irg)%w_int(i,j,k)
+ enddo
enddo
enddo
enddo
enddo
- enddo
!$OMP END PARALLEL DO
- else
- geom%green(:,:,:,:)=0.d0
- ip=0
+ else
+ ip=0
!$OMP PARALLEL DO COLLAPSE(4) DEFAULT(SHARED) PRIVATE(kk,jj,ii,ic,ip,ir,icp,k,j,i,ipint)
- do kk=1,geom%ngllv_shape
- do jj=1,geom%ngllh_shape
- do ii=1,geom%ngllh_shape
- do ic=1,3
- ip=ipnumber(geom,ic,ii,jj,kk,2)
- do ir=1,geom%NBR
- ! do icp=1,3
- do k=1,geom%ngllv_int
- do j=1,geom%ngllh_int
- do i=1,geom%ngllh_int
- ipint=ipnumber(geom,ic,i,j,k,1)
- geom%green(:,:,ir,ip)=geom%green(:,:,ir,ip) &
- +green(:,:,ir,ipint)*geom%w_int(i,j,k)*geom%shape_fcts(i,j,k,ii,jj,kk)
+ do kk=1,geom%ring(irg)%ngllv_shape
+ do jj=1,geom%ring(irg)%ngllh_shape
+ do ii=1,geom%ring(irg)%ngllh_shape
+ do ic=1,3
+ ip=ipnumber(geom,ic,irg,ii,jj,kk,2)
+ do ir=1,geom%NBR
+ do k=1,geom%ring(irg)%ngllv_int
+ do j=1,geom%ring(irg)%ngllh_int
+ do i=1,geom%ring(irg)%ngllh_int
+ ipint=ipnumber(geom,ic,irg,i,j,k,1)
+ geom%green(:,:,ir,ip)=geom%green(:,:,ir,ip) &
+ +green(:,:,ir,ipint)*geom%ring(irg)%w_int(i,j,k)*geom%ring(irg)%shape_fcts(i,j,k,ii,jj,kk)
+ enddo
enddo
enddo
- enddo
! enddo
+ enddo
enddo
enddo
enddo
enddo
- enddo
!$OMP END PARALLEL DO
- endif
+ endif
+ enddo
print*,'summing (2) ...'
traces(:,:,:)=0._DP
do it=1,geom%NBTinv
@@ -543,6 +620,62 @@ contains
!----------------------------------------------------------------------------
end subroutine sum_green_fcts
!----------------------------------------------------------------------------
+!----------------------------------------------------------------------------
+ subroutine get_Mc(geom,ic,irg,ii,jj,kk,Mc)
+!for epsi
+! 1 = xx
+! 2 = yy
+! 3 = zz
+! 4 = xy
+! 5 = xz
+! 6 = yz
+!----------------------------------------------------------------------------
+ use def_gparam
+ implicit none
+ type(meshinv), intent(in) :: geom
+ integer, intent(in) :: ic,irg,ii,jj,kk
+ real(DP), dimension(:), intent(out) :: Mc
+!
+ integer :: a,i,j,k,ipint,b,irad,ijk
+ real(DP), parameter:: c2=dsqrt(2.d0)/2.d0
+ real(DP), dimension(6) :: epsi
+!
+ ijk=geom%ring(irg)%NBPic
+ do k=1,geom%ring(irg)%ngllv_int
+ do j=1,geom%ring(irg)%ngllh_int
+ do i=1,geom%ring(irg)%ngllh_int
+ ijk=ijk+1
+ do a=1,geom%NBCM
+ ipint=ipnumber(geom,a,irg,i,j,k,1)
+ epsi(:)=0.d0
+ select case(ic) !psi component
+ case (1) !x
+ epsi(1)=geom%ring(irg)%deriv(1,1,i,j,k)*geom%ring(irg)%dshape_fcts(1,i,j,k,ii,jj,kk)
+ epsi(4)=geom%ring(irg)%deriv(2,1,i,j,k)*geom%ring(irg)%dshape_fcts(1,i,j,k,ii,jj,kk)*c2
+ epsi(5)=geom%ring(irg)%deriv(3,1,i,j,k)*geom%ring(irg)%dshape_fcts(1,i,j,k,ii,jj,kk)*c2
+ case (2) !y
+ epsi(2)=geom%ring(irg)%deriv(2,2,i,j,k)*geom%ring(irg)%dshape_fcts(2,i,j,k,ii,jj,kk)
+ epsi(4)=geom%ring(irg)%deriv(1,2,i,j,k)*geom%ring(irg)%dshape_fcts(2,i,j,k,ii,jj,kk)*c2
+ epsi(6)=geom%ring(irg)%deriv(3,2,i,j,k)*geom%ring(irg)%dshape_fcts(2,i,j,k,ii,jj,kk)*c2
+ case (3) !z
+ epsi(3)=geom%ring(irg)%deriv(3,3,i,j,k)*geom%ring(irg)%dshape_fcts(3,i,j,k,ii,jj,kk)
+ epsi(5)=geom%ring(irg)%deriv(2,3,i,j,k)*geom%ring(irg)%dshape_fcts(1,i,j,k,ii,jj,kk)*c2
+ epsi(6)=geom%ring(irg)%deriv(3,3,i,j,k)*geom%ring(irg)%dshape_fcts(1,i,j,k,ii,jj,kk)*c2
+ end select
+ irad=geom%iglob2irad(ijk)
+ Mc(ipint)=0.d0
+ do b=1,6
+ Mc(ipint)=Mc(ipint)+geom%CIJ(b,a,irad)*epsi(b)
+ enddo
+ Mc(ipint)=Mc(ipint)*geom%ring(irg)%w_int(i,j,k) !w_int already contains the jacobian
+ enddo
+ enddo
+ enddo
+ enddo
+
+!----------------------------------------------------------------------------
+ end subroutine get_Mc
+!----------------------------------------------------------------------------
!----------------------------------------------------------------------------
subroutine build_ATCA(geom)
!----------------------------------------------------------------------------
@@ -597,7 +730,7 @@ contains
end subroutine build_ATCA
!----------------------------------------------------------------------------
!----------------------------------------------------------------------------
- subroutine get_src_val(geom,coef,time,r,theta,phi,fsrc,icode)
+ subroutine get_src_val(geom,coef,time,irg,r,theta,phi,fsrc,icode)
! icode: 2 equator centered coordinates
! 1 true center coordinates
!----------------------------------------------------------------------------
@@ -607,25 +740,25 @@ contains
real(DP), dimension(:), intent(in) :: coef
real(DP), intent(in) :: time,r,theta,phi
real(DP), dimension(:), intent(out) :: fsrc
- integer, intent(in) :: icode
+ integer, intent(in) :: icode,irg
!
integer :: it,N,ic,i,j,k,ip
real(DP) :: rtil,xi,eta,t,cc
real(DP), dimension(:), allocatable :: ET,VN,Q
- real(DP), dimension(geom%ngllh_shape):: hx,hy
- real(DP), dimension(geom%ngllv_shape):: hz
+ real(DP), dimension(geom%ring(irg)%ngllh_shape):: hx,hy
+ real(DP), dimension(geom%ring(irg)%ngllv_shape):: hz
real(DP), dimension(3,geom%NBTinv) :: tt
!get reference cub coordinates:
- call get_ref_cub_coord(geom,phi,theta,r,xi,eta,rtil,icode)
+ call get_ref_cub_coord(geom,irg,phi,theta,r,xi,eta,rtil,icode)
! hNi values:
- N=geom%ngllh_shape-1
+ N=geom%ring(irg)%ngllh_shape-1
if (N==0) then
hx(1)=1._DP
hy(1)=1._DP
else
allocate(ET(0:N),VN(0:N),Q(0:N))
call ZELEGL(N,ET,VN)
- do i=1,geom%ngllh_shape
+ do i=1,geom%ring(irg)%ngllh_shape
Q(:)=0.d0
Q(i-1)=1.d0
call INLEGL(N,ET,VN,Q,xi ,hx(i) )
@@ -633,13 +766,13 @@ contains
enddo
deallocate(ET,VN,Q)
endif
- N=geom%ngllv_shape-1
+ N=geom%ring(irg)%ngllv_shape-1
if (N==0) then
hz(1)=1._DP
else
allocate(ET(0:N),VN(0:N),Q(0:N))
call ZELEGL(N,ET,VN)
- do i=1,geom%ngllv_shape
+ do i=1,geom%ring(irg)%ngllv_shape
Q(:)=0.d0
Q(i-1)=1.d0
call INLEGL(N,ET,VN,Q,rtil,hz(i) )
@@ -653,11 +786,11 @@ contains
if (abs(geom%tmin+(it-1)*geom%dtinv-time)/geom%tmax < 1.d-7) then
!no time interpolation
fsrc(:)=0.d0
- do k=1,geom%ngllv_shape
- do i=1,geom%ngllh_shape
- do j=1,geom%ngllh_shape
+ do k=1,geom%ring(irg)%ngllv_shape
+ do i=1,geom%ring(irg)%ngllh_shape
+ do j=1,geom%ring(irg)%ngllh_shape
do ic=1,3
- ip=ipnumber(geom,ic,i,j,k,2)+(it-1)*geom%NBPs
+ ip=ipnumber(geom,ic,irg,i,j,k,2)+(it-1)*geom%NBPs
fsrc(ic)=fsrc(ic)+coef(ip)*hx(i)*hy(j)*hz(k)
!if (i==10.and.k==27.and.ic==1) write(701,*) j,coef(ip)
enddo
@@ -668,11 +801,11 @@ contains
!time interpolation needed
tt(:,:)=0._DP
do it=1,geom%NBTinv
- do k=1,geom%ngllv_shape
- do i=1,geom%ngllh_shape
- do j=1,geom%ngllh_shape
+ do k=1,geom%ring(irg)%ngllv_shape
+ do i=1,geom%ring(irg)%ngllh_shape
+ do j=1,geom%ring(irg)%ngllh_shape
do ic=1,3
- ip=ipnumber(geom,ic,i,j,k,2)+(it-1)*geom%NBPs
+ ip=ipnumber(geom,ic,irg,i,j,k,2)+(it-1)*geom%NBPs
tt(ic,it)=tt(ic,it)+coef(ip)*hx(i)*hy(j)*hz(k)
enddo
enddo
@@ -699,75 +832,97 @@ contains
implicit none
type(meshinv), intent(inout) :: geom
!
- real(DP), dimension(geom%ngllh_shape) :: xgllh_shape
- real(DP), dimension(geom%ngllv_shape) :: xgllv_shape
- real(DP), dimension(geom%ngllh_int) :: xgllh_int,wgllh_int
- real(DP), dimension(geom%ngllv_int) :: xgllv_int,wgllv_int
- integer :: i,j,k,a,ijk
- real(DP) :: r,nu,eta,X,Y,isdelta,Jac
- !
- call def_xgll(xgllh_shape,geom%ngllh_shape)
- call def_xgll(xgllv_shape,geom%ngllv_shape)
- if (geom%ngllv_int>0) then
- call def_xgll(xgllh_int, geom%ngllh_int )
- call def_xgll(xgllv_int, geom%ngllv_int )
+ real(DP), dimension(:), allocatable :: xgllh_shape
+ real(DP), dimension(:), allocatable :: xgllv_shape
+ real(DP), dimension(:), allocatable :: xgllh_int,wgllh_int
+ real(DP), dimension(:), allocatable :: xgllv_int,wgllv_int
+ integer :: i,j,k,a,ijk,irg
+ real(DP) :: r,nu,eta,X,Y,isdelta,Jac,eps
+ real(DP), dimension(3) :: xx
!
- call def_wgll(wgllh_int, geom%ngllh_int )
- call def_wgll(wgllv_int, geom%ngllv_int )
- allocate(geom%mesh_int (3,geom%ngllh_int,geom%ngllh_int,geom%ngllv_int))
- allocate(geom%mesh_int_rtp(3,geom%ngllh_int*geom%ngllh_int*geom%ngllv_int))
- allocate(geom%w_int(geom%ngllh_int,geom%ngllh_int,geom%ngllv_int))
- allocate(geom%mesh_rad(geom%ngllv_int),geom%iglob2irad(geom%ngllh_int*geom%ngllh_int*geom%ngllv_int))
+ eps=1.e-8_DP*geom%ring(geom%nring)%r2/2.d0 !small value to catch interfaces
+ allocate(geom%mesh_rad(geom%nbrad_int),geom%iglob2irad(geom%NBPint))
+ allocate(geom%mesh_int_rtp(3,geom%NBPint))
+ do irg=1,geom%nring
+ allocate(xgllh_shape(geom%ring(irg)%ngllh_shape),xgllv_shape(geom%ring(irg)%ngllv_shape))
+ allocate(xgllh_int(geom%ring(irg)%ngllh_int),wgllh_int(geom%ring(irg)%ngllh_int))
+ allocate(xgllv_int(geom%ring(irg)%ngllv_int),wgllv_int(geom%ring(irg)%ngllv_int))
+ call def_xgll(xgllh_shape,geom%ring(irg)%ngllh_shape)
+ call def_xgll(xgllv_shape,geom%ring(irg)%ngllv_shape)
+ if (geom%ring(irg)%ngllv_int>0) then
+ call def_xgll(xgllh_int, geom%ring(irg)%ngllh_int )
+ call def_xgll(xgllv_int, geom%ring(irg)%ngllv_int )
+!
+ call def_wgll(wgllh_int, geom%ring(irg)%ngllh_int )
+ call def_wgll(wgllv_int, geom%ring(irg)%ngllv_int )
+ allocate(geom%ring(irg)%mesh_int (3,geom%ring(irg)%ngllh_int,geom%ring(irg)%ngllh_int,geom%ring(irg)%ngllv_int))
+ allocate(geom%ring(irg)%deriv (3,3,geom%ring(irg)%ngllh_int,geom%ring(irg)%ngllh_int,geom%ring(irg)%ngllv_int))
+ allocate(geom%ring(irg)%w_int(geom%ring(irg)%ngllh_int,geom%ring(irg)%ngllh_int,geom%ring(irg)%ngllv_int))
endif
- allocate(geom%mesh_shape(3,geom%ngllh_shape,geom%ngllh_shape,geom%ngllv_shape))
- if (geom%ngllv_shape==1.and.geom%ngllh_shape==1.and.geom%ngllv_int==1.and.geom%ngllh_int==1) then
- geom%mesh_shape(1,1,1,1)=(geom%r2+geom%r1)/2._DP
- geom%mesh_shape(2,1,1,1)=geom%theta_center
- geom%mesh_shape(3,1,1,1)=geom%phi_center
- geom%mesh_int_rtp(1,1)=geom%mesh_shape(1,1,1,1)
+ allocate(geom%ring(irg)%mesh_shape(3,geom%ring(irg)%ngllh_shape,geom%ring(irg)%ngllh_shape,geom%ring(irg)%ngllv_shape))
+ if (geom%ring(irg)%ngllv_shape==1.and.geom%ring(irg)%ngllh_shape==1.and.geom%ring(irg)%ngllv_int==1.and.geom%ring(irg)%ngllh_int==1) then
+ geom%ring(irg)%mesh_shape(1,1,1,1)=(geom%ring(irg)%r2+geom%ring(irg)%r1)/2._DP
+ geom%ring(irg)%mesh_shape(2,1,1,1)=geom%theta_center
+ geom%ring(irg)%mesh_shape(3,1,1,1)=geom%phi_center
+ geom%mesh_int_rtp(1,1)=geom%ring(irg)%mesh_shape(1,1,1,1)
!_rtp must be in degree:
- geom%mesh_int_rtp(2:3,1)=geom%mesh_shape(2:3,1,1,1)*rad2deg
- geom%w_int(1,1,1)=1.d0
-! call rotate2center(geom%theta_center,geom%phi_center,geom%mesh_shape(:,1,1,1))
+ geom%mesh_int_rtp(2:3,1)=geom%ring(irg)%mesh_shape(2:3,1,1,1)*rad2deg
+ geom%ring(irg)%w_int(1,1,1)=1.d0
+! call rotate2center(geom%ring(irg)%theta_center,geom%ring(irg)%phi_center,geom%ring(irg)%mesh_shape(:,1,1,1))
else
- do k=1,geom%ngllv_shape
- r=(xgllv_shape(k)*(geom%r2-geom%r1)+geom%r2+geom%r1)/2.d0
- do j=1,geom%ngllh_shape
+ do k=1,geom%ring(irg)%ngllv_shape
+ r=(xgllv_shape(k)*(geom%ring(irg)%r2-geom%ring(irg)%r1)+geom%ring(irg)%r2+geom%ring(irg)%r1)/2.d0
+ do j=1,geom%ring(irg)%ngllh_shape
nu=xgllh_shape(j)*geom%dtheta/2._DP
Y=tan(nu)
- do i=1,geom%ngllh_shape
+ do i=1,geom%ring(irg)%ngllh_shape
eta=xgllh_shape(i)*geom%dphi/2._DP
X=tan(eta)
isdelta=1._DP/sqrt(1+X**2+Y**2)
- geom%mesh_shape(1,i,j,k)=r *isdelta
- geom%mesh_shape(2,i,j,k)=r*X*isdelta
- geom%mesh_shape(3,i,j,k)=r*Y*isdelta
- call rotate2center(geom%theta_center,geom%phi_center,geom%mesh_shape(:,i,j,k))
+ geom%ring(irg)%mesh_shape(1,i,j,k)=r *isdelta
+ geom%ring(irg)%mesh_shape(2,i,j,k)=r*X*isdelta
+ geom%ring(irg)%mesh_shape(3,i,j,k)=r*Y*isdelta
+ call rotate2center(geom%theta_center,geom%phi_center,geom%ring(irg)%mesh_shape(:,i,j,k))
enddo
enddo
enddo
- if (geom%ngllv_int>0) then
- ijk=0
- do k=1,geom%ngllv_int
- r=(xgllv_int(k)*(geom%r2-geom%r1)+geom%r2+geom%r1)/2.d0
- do j=1,geom%ngllh_int
+ if (geom%ring(irg)%ngllv_int>0) then
+ ijk=geom%ring(irg)%NBPic
+ do k=1,geom%ring(irg)%ngllv_int
+ r=(xgllv_int(k)*(geom%ring(irg)%r2-geom%ring(irg)%r1)+geom%ring(irg)%r2+geom%ring(irg)%r1)/2.d0
+ if (k==1) r=r+eps !catching interfaces
+ if (k==geom%ring(irg)%ngllv_int) r=r-eps !catching interfaces
+ do j=1,geom%ring(irg)%ngllh_int
nu=xgllh_int(j)*geom%dtheta/2._DP
Y=tan(nu)
- do i=1,geom%ngllh_int
+ do i=1,geom%ring(irg)%ngllh_int
eta=xgllh_int(i)*geom%dphi/2._DP
X=tan(eta)
isdelta=1._DP/sqrt(1+X**2+Y**2)
- geom%mesh_int(1,i,j,k)=r *isdelta
- geom%mesh_int(2,i,j,k)=r*X*isdelta
- geom%mesh_int(3,i,j,k)=r*Y*isdelta
+ xx(1)=r *isdelta
+ xx(2)=r*X*isdelta
+ xx(3)=r*Y*isdelta
+ geom%ring(irg)%mesh_int(1,i,j,k)=xx(1)
+ geom%ring(irg)%mesh_int(2,i,j,k)=xx(2)
+ geom%ring(irg)%mesh_int(3,i,j,k)=xx(3)
+ geom%ring(irg)%deriv(1,1,i,j,k) = -xx(2)/( xx(1)**2 + xx(2)**2 ) ! xi_x
+ geom%ring(irg)%deriv(2,1,i,j,k) = xx(1)/( xx(1)**2 + xx(2)**2 ) ! xi_y
+ geom%ring(irg)%deriv(3,1,i,j,k) = 0._DP !xi_z
+ geom%ring(irg)%deriv(1,2,i,j,k) = -xx(3)/( xx(1)**2 + xx(3)**2 ) !et_x
+ geom%ring(irg)%deriv(2,2,i,j,k) = 0._DP ! et_y
+ geom%ring(irg)%deriv(3,2,i,j,k) = xx(1)/( xx(1)**2 + xx(3)**2 ) ! et_z
+ geom%ring(irg)%deriv(1,3,i,j,k) = xx(1)/r ! ga_x
+ geom%ring(irg)%deriv(2,3,i,j,k) = xx(2)/r ! ga_y
+ geom%ring(irg)%deriv(3,3,i,j,k) = xx(3)/r ! ga_z
Jac=r**2*(1+X**2)*(1+Y**2)/(1+X**2+Y**2)**1.5
- geom%w_int(i,j,k)=wgllh_int(i)*wgllh_int(j)*wgllh_int(k)*Jac
- call rotate2center(geom%theta_center,geom%phi_center,geom%mesh_int(:,i,j,k))
+ geom%ring(irg)%w_int(i,j,k)=wgllh_int(i)*wgllh_int(j)*wgllv_int(k)*Jac
+ call rotate2center(geom%theta_center,geom%phi_center,geom%ring(irg)%mesh_int(:,i,j,k))
ijk=ijk+1
- call xyz2rtp(geom%mesh_int(:,i,j,k),geom%mesh_int_rtp(1,ijk),geom%mesh_int_rtp(2,ijk) &
+ call xyz2rtp(geom%ring(irg)%mesh_int(:,i,j,k),geom%mesh_int_rtp(1,ijk),geom%mesh_int_rtp(2,ijk) &
,geom%mesh_int_rtp(3,ijk),deg=.true.)
- if (j==1.and.i==1) geom%mesh_rad(k)=geom%mesh_int_rtp(1,ijk)
- geom%iglob2irad(ijk)=k
+ if (j==1.and.i==1) geom%mesh_rad(k+geom%ring(irg)%NBradic)=geom%mesh_int_rtp(1,ijk)
+ geom%iglob2irad(ijk)=k+geom%ring(irg)%NBradic
+! print*,'hooo',ijk,geom%iglob2irad(ijk)
! geom%mesh_int_rtp(1,ijk)=(RA-geom%mesh_int_rtp(1,ijk))/1000.d0 !depth in meter
! print*,geom%mesh_int(:,i,j,k),geom%mesh_int_rtp(:,ijk)
! stop
@@ -776,21 +931,31 @@ contains
enddo
endif
endif
-!
+ deallocate(xgllh_shape,xgllv_shape)
+ deallocate(xgllh_int,wgllh_int)
+ deallocate(xgllv_int,wgllv_int)
+ enddo
+ call get_CIJvalues(geom)
+ !
+! do ijk=1,geom%NBPint
+! print*,'raaa',ijk,geom%iglob2irad(ijk)
+! enddo
end subroutine generate_invmesh
!------------------------------------------------------------
subroutine write_mesh2geo(geom)
implicit none
type(meshinv), intent(inout) :: geom
- integer :: i,j,k,l
+ integer :: irg,i,j,k,l
open(11,file="mesh_shape.geo",status="replace")
write(11,421) 1,0.,0.,0.
l=1
- do k=1,geom%ngllv_shape
- do j=1,geom%ngllh_shape
- do i=1,geom%ngllh_shape
- l=l+1
- write(11,421) l,geom%mesh_shape(1,i,j,k),geom%mesh_shape(2,i,j,k),geom%mesh_shape(3,i,j,k)
+ do irg=1,geom%nring
+ do k=1,geom%ring(irg)%ngllv_shape
+ do j=1,geom%ring(irg)%ngllh_shape
+ do i=1,geom%ring(irg)%ngllh_shape
+ l=l+1
+ write(11,421) l,geom%ring(irg)%mesh_shape(1,i,j,k),geom%ring(irg)%mesh_shape(2,i,j,k),geom%ring(irg)%mesh_shape(3,i,j,k)
+ enddo
enddo
enddo
enddo
@@ -951,7 +1116,7 @@ contains
end subroutine rotate2centerm1
!-------------------------------------------------
!-------------------------------------------------
- subroutine get_ref_cub_coord(geom,phi,theta,r,xi,eta,rtilde,icode)
+ subroutine get_ref_cub_coord(geom,irg,phi,theta,r,xi,eta,rtilde,icode)
!get reference cube point coordinates centered on (phi=0.,theta=90)
!icode : 1 in the rotated to center frame
! 2 in the reference frame
@@ -961,7 +1126,7 @@ contains
type(meshinv), intent(inout) :: geom
real(DP), intent(in) :: r,theta,phi
real(DP), intent(out) :: xi,eta,rtilde
- integer, intent(in) :: icode
+ integer, intent(in) :: icode,irg
!
real(DP), dimension(3) :: pt
real(DP) :: sdelta,X,Y
@@ -978,7 +1143,7 @@ contains
sdelta=r/pt(1)
X=pt(2)/r*sdelta
Y=pt(3)/r*sdelta
- rtilde=(2.*r-geom%r1-geom%r2)/(geom%r2-geom%r1)
+ rtilde=(2.*r-geom%ring(irg)%r1-geom%ring(irg)%r2)/(geom%ring(irg)%r2-geom%ring(irg)%r1)
xi =atan(X)*2._DP/geom%dphi
eta =atan(Y)*2._DP/geom%dtheta
!-------------------------------------------------
@@ -1095,12 +1260,16 @@ contains
type(meshinv), intent(in) :: geom
real(DP), dimension(:), intent(in) :: coef
!
- integer :: i,j,k,ii,jj,kk
+ integer :: i,j,k,ii,jj,kk,irg
!
! open(11,file=name,form='unformatted',status='replace',action='write')
open(11,file=name,status='replace',action='write')
- write(11,*) geom%NBP,geom%NBPs,geom%NBTinv,geom%dtinv,geom%tmin,geom%tmax,geom%fmax,geom%fmaxb,geom%ngllv_shape,geom%ngllh_shape
- write(11,*) geom%ngllv_int, geom%ngllh_int,geom%phi_center,geom%theta_center,geom%dphi ,geom%dtheta ,geom%r1,geom%r2
+ write(11,*) geom%NBP,geom%NBPs,geom%NBTinv,geom%dtinv,geom%tmin,geom%tmax,geom%fmax,geom%fmaxb,geom%nring
+ write(11,*) geom%phi_center,geom%theta_center,geom%dphi ,geom%dtheta
+ do irg=1,geom%nring
+ write(11,*) geom%ring(irg)%ngllv_shape,geom%ring(irg)%ngllh_shape
+ write(11,*) geom%ring(irg)%ngllv_int, geom%ring(irg)%ngllh_int ,geom%ring(irg)%r1,geom%ring(irg)%r2
+ enddo
write(11,*) geom%t0,geom%time_support,geom%fmax
write(11,*) (coef(i),i=1,geom%NBP)
close(11)
@@ -1115,22 +1284,39 @@ contains
character(len=*), intent(in) :: name
type(meshinv), intent(out) :: geom
real(DP), dimension(:), pointer, intent(out) :: coef
- integer :: i
+ integer :: i,irg
!
! open(11,file=name,form='unformatted',status='old',action='read')
open(11,file=name,status='old',action='read')
- read(11,*) geom%NBP,geom%NBPs,geom%NBTinv,geom%dtinv,geom%tmin,geom%tmax,geom%fmax,geom%fmaxb,geom%ngllv_shape,geom%ngllh_shape
- read(11,*) geom%ngllv_int, geom%ngllh_int,geom%phi_center,geom%theta_center,geom%dphi ,geom%dtheta ,geom%r1,geom%r2
- read(11,*) geom%t0,geom%time_support,geom%fmax
+ read(11,*) geom%NBP,geom%NBPs,geom%NBTinv,geom%dtinv,geom%tmin,geom%tmax,geom%fmax,geom%fmaxb,geom%nring
+ read(11,*) geom%phi_center,geom%theta_center,geom%dphi ,geom%dtheta
+ allocate(geom%ring(geom%nring))
+ do irg=1,geom%nring
+ read(11,*) geom%ring(irg)%ngllv_shape,geom%ring(irg)%ngllh_shape
+ read(11,*) geom%ring(irg)%ngllv_int, geom%ring(irg)%ngllh_int ,geom%ring(irg)%r1,geom%ring(irg)%r2
+ enddo
if (geom%t0<geom%time_support/2) print*, 'Warning: t0 should be greater than half duration'
geom%tmin=max(0.d0,geom%t0-geom%time_support/2.d0)
geom%tmax=geom%t0+geom%time_support/2.d0
geom%dtinv =1.d0/geom%fmax/2.2d0
geom%NBTinv=int((geom%tmax-geom%tmin)/geom%dtinv)+1
geom%tmax=geom%tmin+(geom%NBTinv-1)*geom%dtinv !rounding tmax to an exact mutiple of dtinv
- geom%fmaxb=geom%fmax*1.3
- allocate(coef(geom%NBP))
+ geom%fmaxb=geom%fmax*1.3
+ read(11,*) geom%t0,geom%time_support,geom%fmax
read(11,*) (coef(i),i=1,geom%NBP)
+
+!!$ read(11,*) geom%NBP,geom%NBPs,geom%NBTinv,geom%dtinv,geom%tmin,geom%tmax,geom%fmax,geom%fmaxb,geom%ngllv_shape,geom%ngllh_shape
+!!$ read(11,*) geom%ngllv_int, geom%ngllh_int,geom%phi_center,geom%theta_center,geom%dphi ,geom%dtheta ,geom%r1,geom%r2
+!!$ read(11,*) geom%t0,geom%time_support,geom%fmax
+!!$ if (geom%t0<geom%time_support/2) print*, 'Warning: t0 should be greater than half duration'
+!!$ geom%tmin=max(0.d0,geom%t0-geom%time_support/2.d0)
+!!$ geom%tmax=geom%t0+geom%time_support/2.d0
+!!$ geom%dtinv =1.d0/geom%fmax/2.2d0
+!!$ geom%NBTinv=int((geom%tmax-geom%tmin)/geom%dtinv)+1
+!!$ geom%tmax=geom%tmin+(geom%NBTinv-1)*geom%dtinv !rounding tmax to an exact mutiple of dtinv
+!!$ geom%fmaxb=geom%fmax*1.3
+!!$ allocate(coef(geom%NBP))
+!!$ read(11,*) (coef(i),i=1,geom%NBP)
close(11)
call get_shape_fcts(geom)
!-----------------------------------------------------------
@@ -1218,20 +1404,30 @@ contains
subroutine read_data(geom,data)
!-------------------------------------------------------------
use def_gparam
- use global_main_param, only:rotation,stations_name,comp,compr,NBT,dt
+ use global_main_param, only:rotation,stations_name,comp,compr
implicit none
type(meshinv), intent(inout) :: geom
real(DP), dimension(:,:,:), pointer, intent(out) :: data
!
real(DP), dimension(:), allocatable :: tmp
- real(DP) :: bid
- integer :: i,ir,ic
+ real(DP) :: bid,dtdata,t1
+ integer :: i,ir,ic,iend,NBTdata
character(len=100) :: name
-!to update:
- geom%NBTdata=NBT
- geom%dtdata=dt
-!
- allocate(data(geom%NBTdata,3,geom%NBR))
+!get DATA dt and NBT:
+ write(name,'("U",a1,"_",a4)') compr(1),stations_name(1)
+ open(17,file=add_dir(name,geom%datadir),status='old',action='read')
+ iend=0
+ NBTdata=0
+ do while (iend==0)
+ read(17,'(f8.1,1x,e14.7)',iostat=iend) t1,bid
+ if (iend>0) STOP 'IO error in read_data'
+ if (iend==0) NBTdata=NBTdata+1
+ if (NBTdata==1) dtdata=t1
+ if (NBTdata==2) dtdata=t1-dtdata
+ enddo
+ close(17)
+ print*,'NBT and dt for data files:',NBTdata,dtdata
+ allocate(data(geom%NBTtrace,3,geom%NBR))
do ir=1,geom%NBR
do ic=1,3
if (rotation) then
@@ -1239,14 +1435,12 @@ contains
else
write(name,'("U",a1,"_",a4)') comp(ic) ,stations_name(ir)
endif
- allocate(tmp(NBT))
+ allocate(tmp(NBTdata))
open(17,file=add_dir(name,geom%datadir),status='old',action='read')
- read(17,'(f8.1,1x,e14.7)') (bid,tmp(i),i=1,geom%NBTdata)
-!ajouter appodisation:
-
+ read(17,'(f8.1,1x,e14.7)') (bid,tmp(i),i=1,NBTdata)
close(17)
- if (abs((geom%dtdata-dt)/(geom%dtdata+dt))>1.d-6.or.geom%NBTdata/=NBT) then
- call resample(tmp,geom%t0data,geom%dtdata,geom%NBTdata,data(:,ic,ir),dt,NBT)
+ if (abs((dtdata-geom%dttrace)/(dtdata+geom%dttrace))>1.d-6.or.geom%NBTtrace/=NBTdata) then
+ call resample(tmp,0.d0,dtdata,NBTdata,data(:,ic,ir),geom%dttrace,geom%NBTtrace)
else
data(:,ic,ir)=tmp(:)
endif
@@ -1314,7 +1508,51 @@ function add_dir(name,dir)
!-------------------------------------------------------------
end function add_dir
!-------------------------------------------------------------
+!-------------------------------------------------------------
+subroutine get_CIJvalues(geom)
+!-------------------------------------------------------------
+ use earth_modele
+ implicit none
+ type(meshinv), intent(inout) :: geom
+!
+ real(DP), dimension(geom%nbrad_int) :: rho,Vsh,Vsv,Vph,Vpv,eta,Qs
+ real(DP) :: A,C,F,L,N
+ integer :: i,k,j
+!
+ allocate(geom%CIJ(6,6,geom%nbrad_int))
+ geom%CIJ(:,:,:)=0._DP
+!mesh_rad has been designed to catch interfaces (by moving rads by eps)
+ call get_1dmodel_value(geom%modele1d,rho,Vsh,Vsv,Vph,Vpv,eta,Qs,geom%nbrad_int,geom%mesh_rad)
+ do i=1,geom%nbrad_int
+! write(111,*) geom%mesh_rad(i),vph(i)
+ A=rho(i)*vph(i)**2
+ C=rho(i)*vpv(i)**2
+ L=rho(i)*vsv(i)**2
+ N=rho(i)*vsh(i)**2
+ F=eta(i)*(A-2.d0*L)
+!123=rtp:
+ geom%CIJ(1,1,i)=C
+ geom%CIJ(2,2,i)=A
+ geom%CIJ(3,3,i)=A
+ geom%CIJ(1,2,i)=F
+ geom%CIJ(1,3,i)=F
+ geom%CIJ(2,3,i)=A-2.d0*N
+ geom%CIJ(4,4,i)=2.d0*N !convention de kelvin et non de volt
+ geom%CIJ(5,5,i)=2.d0*L
+ geom%CIJ(6,6,i)=2.d0*L
+!symetrie
+ do k=2,6
+ do j=1,k-1
+ geom%CIJ(k,j,i)=geom%CIJ(j,k,i)
+ enddo
+ enddo
+ enddo
+! close(111); stop 'test'
+!
+!-------------------------------------------------------------
+end subroutine get_CIJvalues
+!-------------------------------------------------------------
- !-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
end module module_srcinv
!-----------------------------------------------------------------------------
diff --git a/nms/speclib.f90 b/nms/speclib.f90
index ed083ae5a8cb036afc68eb802614ef46a121f508..b251d4812b1deaf2b803c609d06a114ee5a4d03f 100644
--- a/nms/speclib.f90
+++ b/nms/speclib.f90
@@ -2,7 +2,7 @@
MODULE FUNARO
!====================
!
- public :: def_xgll, def_wgll, def_heta, def_deriv,def_derivx ,INLEGL,ZELEGL
+ public :: def_xgll, def_wgll, def_heta, def_deriv,def_derivx ,INLEGL,INLEGL2,ZELEGL
private
!
CONTAINS
diff --git a/nms/srcdirect.f90 b/nms/srcdirect.f90
index e0771afea5e7fbf6cff349cb0f6c466a7f1175ca..9172b6f1cbab5e4e2ac7b5423096b76a561250c4 100644
--- a/nms/srcdirect.f90
+++ b/nms/srcdirect.f90
@@ -21,7 +21,7 @@ program srcdirect
! print*,'Summing'
! call sum_green(geom,coef,'direct')
!
- call get_pt_time_function(geom,coef,1,2,2,2,test)
+ call get_pt_time_function(geom,coef,1,1,2,2,2,test)
open(11,file='fctpb.txt')
do i=1,geom%NBTtrace
write(11,*) (i-1)*geom%dttrace,test(i)
diff --git a/nms/srcinv.f90 b/nms/srcinv.f90
index 1443187b7435d06b17c6e7047dde105c390c80f1..77b20effb1e853c8b6de08d604e901af20ab8023 100644
--- a/nms/srcinv.f90
+++ b/nms/srcinv.f90
@@ -28,7 +28,7 @@
call invert_ATCA(geom,data,coef)
call comp_misfit(geom,data,coef)
allocate(cc(geom%NBTinv))
- call get_pt_time_function(geom,coef,1,1,1,1,test,cc_=cc)
+ call get_pt_time_function(geom,coef,1,1,1,1,1,test,cc_=cc)
do i=1,geom%NBTinv
write(421,*) (i-1)*geom%dtinv+geom%tmin,cc(i)
enddo
@@ -38,14 +38,14 @@
enddo
close(11)
deallocate(test)
- call get_pt_time_function(geom,coef,2,1,1,1,test)
+ call get_pt_time_function(geom,coef,2,1,1,1,1,test)
open(11,file='fctp2.txt')
do i=1,geom%NBTtrace
write(11,*) (i-1)*geom%dttrace,test(i)
enddo
close(11)
deallocate(test)
- call get_pt_time_function(geom,coef,3,1,1,1,test)
+ call get_pt_time_function(geom,coef,3,1,1,1,1,test)
open(11,file='fctp3.txt')
do i=1,geom%NBTtrace
write(11,*) (i-1)*geom%dttrace,test(i)
diff --git a/nms/util_fctp.f90 b/nms/util_fctp.f90
index 762c10ed7338b4ccf684f4cbf3f0b427aa2aa4a5..6745dd5b6ac3a937e649923a1613917611ef1aaf 100644
--- a/nms/util_fctp.f90
+++ b/nms/util_fctp.f90
@@ -457,12 +457,12 @@ contains
integer, intent(in) :: Lmax
real(DP), intent(out):: fmax
integer :: l
-!on ne cherche que sur le mode fondamentale spheroidale
+!on ne cherche que sur le mode fondamentale spheroidale
+ l=Lmax+1
if (l>LLmaxS) then
print*,'get_fmaxc: Attention, l> LLmaxS:',l,LLmaxS
fmax=omtabS(LLmaxS,0)/2._DP/PI
else
- l=Lmax+1
do while(irecS(l,0)==0)
l=l-1
if (l<0) stop 'get_fmaxc: Ca va pas ....!'