Commit 27ec9c4a authored by Yann CAPDEVILLE's avatar Yann CAPDEVILLE
Browse files

little modifcation to fix issue with prem 0T2

 On branch Yann
 Your branch is up to date with 'origin/Yann'.
 Changes to be committed:
	modified:   Changlog
	modified:   nms/get_fctp.f90
	modified:   nms/global_main.f90
	modified:   nms/module_modes.f90
	modified:   nms/util_fctp.f90
	modified:   yannos/util_minos.f90
parent d969aa5c
......@@ -91,3 +91,5 @@ tagging rel-1-4
and one local approximation (corresponding to what is measuring
a sensor?)
documentation update
10/07/20 fixing little bug for 0T2 in prem. Not sure if it has concequences for other model.
adding the option to dump gravity potential in get_fctp
......@@ -33,7 +33,7 @@ program get_fctp
write(*,*) "Usage : "//trim(adjustl(prog_name))//" eigenfile_prefix function_code n l"
write(*,*) "With : "
print*, "- eigenfile_prefix: eigen function file prefix"
print*, "- function_code: 1=U,11,Up,2=V,22=Vp,3=W,33=Wp"
print*, "- function_code: 1=U,11,Up,2=V,22=Vp,3=W,33=Wp,4=P,44=Pp"
print*, "- n: radial order"
print*, "- l: angular order"
print*, "Remark: if n<0 then all the available n will be output"
......@@ -65,6 +65,10 @@ program get_fctp
write(name,'(i4.4,"U",i4.4)') n,l
case(11)
write(name,'(i4.4,"Up",i4.4)') n,l
case(4)
write(name,'(i4.4,"P",i4.4)') n,l
case(44)
write(name,'(i4.4,"Pp",i4.4)') n,l
case(2)
print*,'Warning, output will be multiplied by sqrt(l(l+1))'
write(name,'(i4.4,"V",i4.4)') n,l
......
......@@ -37,6 +37,7 @@ module global_main_param
logical :: sources_available=.false.,stacking_sources,source_delay,source_force
doubleprecision, parameter :: SEXP=1.d18
character(len=4), dimension(: ), allocatable :: stations_name
!character(len=6), dimension(: ), allocatable :: stations_name
!temps:
integer :: NBT
real(DP):: dt,t0,amp_source
......
......@@ -27,6 +27,7 @@ module modes
!
integer :: length_qnl,length_qnl_all,nbchunk
integer, parameter :: memorymax=512
! integer, parameter :: memorymax=128000000
integer, dimension(:), allocatable :: qnldeb_chunk,qnlfin_chunk
integer :: qnldeb,qnlfin,iqnld
!
......@@ -332,6 +333,7 @@ contains
write(name,'("U",a1,"_",a4)') compr(ic),stations_name(ir)
else
write(name,'("U",a1,"_",a4)') comp(ic) ,stations_name(ir)
!write(name,'("U",a1,"_",a6)') comp(ic) ,stations_name(ir)
endif
open(17,file=name)
write(17,'(f8.1,1x,e14.7)') ((i-1)*dt,trace_time(i,ic),i=1,NBT)
......
......@@ -777,10 +777,15 @@ contains
use global_main_param, only: RA
implicit none
integer, intent(in) :: unit,n,l,s
real(DP), dimension(4,nkS) :: xdum
real(DP), dimension(:,:), allocatable :: xdum
integer :: i,iflag
if (.not.files_openS) stop 'write_modeS_ascii: use open_fctp1 with S type first'
if (s>22) then
allocate(xdum(6,nkS))
else
allocate(xdum(4,nkS))
endif
call readS(n,l,xdum,iflag)
!print*,maxval(abs(xdum)),ubound(xdum)
if (s.eq.1) then
......@@ -799,6 +804,14 @@ contains
do i=1,nkS
write(unit,*)(RA-rayS(i)),sqrt(float(l*(l+1)))*xdum(4,i)
enddo
else if (s.eq.4) then
do i=1,nkS
write(unit,*)(RA-rayS(i)),xdum(5,i)
enddo
else if (s.eq.44) then
do i=1,nkS
write(unit,*)(RA-rayS(i)),xdum(6,i)
enddo
endif
!-----------------------------------------------------------------
end subroutine write_modeS_ascii
......@@ -979,7 +992,7 @@ contains
implicit none
integer, intent(in) :: l,n
integer, intent(out):: iflag
real(DP), dimension(4,nkS), intent(out) :: xdum
real(DP), dimension(:,:), intent(out) :: xdum
!
real(DP) :: WW,QQ,GC
real(DP), dimension(nkS*6) :: tt
......@@ -989,8 +1002,12 @@ contains
nb=nkS
!nve!est la longueur necessaire pour atteindre la profondeur
!desire pour v, v' u, u' ( on ne lit pas p et p', si on en a besoin il
!ecrire une autre routine de lecture.)
nvec=int(nb*4)
!ecrire une autre routine de lecture.)
if (ubound(xdum,dim=1)==6) then
nvec=int(nb*6)
else
nvec=int(nb*4)
endif
if (l<=LLmaxS.and.n<=NNmaxS) then
if (irecS(l,n)/=0) then
!reading and killing discontinuities
......@@ -1001,6 +1018,12 @@ contains
xdum(3,i)=tt(4*(nkS-i)+3)
xdum(4,i)=tt(4*(nkS-i)+4)
enddo
if (ubound(xdum,dim=1)==6) then
do i=1,nkS
xdum(5,i)=tt((nkS-i)+4*nkS+1)
xdum(6,i)=tt((nkS-i)+5*nkS+1)
enddo
endif
iflag=0
else
iflag=2
......
......@@ -1052,7 +1052,6 @@ end subroutine get_fctpL
endif
else
if(knsw==1) then
!test
knt=kount-1
if(jcom/=4.and.l/=1)then
irem=mod(knt,2)
......@@ -1060,6 +1059,7 @@ end subroutine get_fctpL
(irem.ne.0.and.det.gt.0.d0))) knt=knt+1
endif
endif
if (l==2) knt=knt-1 !unkonw problem in prem for l=2 Toroidal modes ..
endif
else
!l=0
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment