Index: sm/trunk-jpl/src/c/modules/FourierLoveCorex/love_numbers.f90.bak =================================================================== --- /issm/trunk-jpl/src/c/modules/FourierLoveCorex/love_numbers.f90.bak (revision 22007) +++ (revision ) @@ -1,362 +1,0 @@ - program lnb_setup - - use lnb_param - use util - - integer :: source_type, fr, i - double precision :: T,cst, fluid_tscale,g - double complex, allocatable :: hs(:,:), ls(:,:), ks(:,:) - double precision, allocatable :: frequencies(:) - - degmax=180 - degmin=1 - nfreq=500 - nlayer=6 - allocate( radbc(nlayer+1), soliddim(nlayer) ) - allocate( roc(nlayer), lac(nlayer), muc(nlayer), muc2(nlayer), vic(nlayer), vic2(nlayer)) - allocate(radius(nlayer+1), burgers(nlayer), frequencies(nfreq)) - allocate( hs(nfreq,degmax+1), ls(nfreq,degmax+1), ks(nfreq,degmax+1) ) - - pi=dacos(-1.d0) - - allow_layer_del=.true. - - - radius(1) = 100.d3; radius(2) = 1221.5d3; radius(3) = 3480.d3; - radius(4) = 5701.d3; radius(5) = 5951.d3; radius(6) = 6271.d3; radius(7) = 6371.d3 - vic(1) = 0.d0; vic(2) = 0.d0; vic(3) = 2.d22; vic(4) = 5.d20; vic(5) = 5.d20; vic(6) = 1.d40 - incompressible=.false. - !ra=6378.137d0 - muc(1) = 1.d6; muc(2) = 0.d0; muc(3) = 2.20425d11 - muc(4) = 0.870257d11; muc(5) = 0.870257d11; muc(6) = 0.510268d11 - roc(1) = 10987.9d0; roc(2) = 10986.9d0; roc(3) = 4903.58d0 - roc(4) = 3628.29d0; roc(5) = 3628.29d0; roc(6) = 3113.94d0 - roc(6)=3300.d0 - - - !roc(6)=3300.d0 - muc2=muc;vic2=vic;lac=1e11 - !burgers(:)=.false. - soliddim(:) = .true. - soliddim(2) = .false. - ra=radius(nlayer+1) - !ra = 6378.137d3 - rb = 3480.d3 - rc = 1221.5d3 - - go0=10 - mu0=1e11 - r0=ra - - 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 - - source_type=11 - cst = 365.25d0*24.d0*3600.d0*1000.d0 - T=8*cst - burgers(1:6)=.false. - print*,'ok',T, burgers(1:6), ' ', soliddim(1:nlayer), radbc(1:nlayer),ra, rb, rc - - - do i=1,nlayer - print*, radius(i+1), roc(i), muc(i), lac(i), vic(i) - print*, '' - enddo - benchmark_spada=.false. - !frequencies=1E-5/cst*2*pi - - do fr=1,nfreq - 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 - frequencies(fr)=2.d0*pi*(10.d0**(fluid_tscale)/cst) -! it is meant to approximate the theoretical freq=0.d0 that should be assessed for fft - elseif (fr <= nfreq/2+1) then - frequencies(fr)=(dble(fr-1)/T)*2.d0*pi - else - frequencies(fr)=(-dble(nfreq+1-fr)/T)*2.d0*pi - end if - enddo - - display=.true. - call love_numbers(frequencies,source_type,hs,ls,ks) - !print*, hs - do i=1,degmax - do fr=1,nfreq - write(17,*), real(hs(fr,i+1)),aimag(hs(fr,i+1)), real(ls(fr,i+1)),aimag(ls(fr,i+1)), real(ks(fr,i+1)),aimag(ks(fr,i+1)) - enddo - enddo - - -end program - -!==================== - subroutine love_numbers(frequencies,source_type,hs,ls,ks) -!==================== - - use lnb_param - 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 - 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 :: i,j,k,n, fr,m,l, source_type - double complex :: hs(nfreq,degmax+1), ls(nfreq,degmax+1), ks(nfreq,degmax+1) - - - double precision :: cpu_time1(100) - integer :: cpu_count - - - pi=acos(-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 - - - - ! Test number of interfaces and of boundary conditions - nbc = 0 - do i=1,nlayer+1 - if (radbc(i)) nbc=nbc+6 - enddo - - - ifmin=source_type - ifmax=source_type - - - print*,'model init' - call model_init2(nbc,nlayer) - print*, 'done' - - nfext = 13 ! Number of potential excitation sources - allocate( bc(nbc,nbc), indx(nbc) ) - allocate( f(nfext,nbc), 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' - - - 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 - - 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='(A23, 10X,A2,I5,A1,I5)',ADVANCE='NO'), repeat(char(8),23),'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=complex(0.d0,frequencies(fr)) - !print*,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 - 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) - - ! Automatic reduction of the number of layers when the attenuation with depth becomes too strong - if (allow_layer_del.eqv..true.) then - - do while ((layerrap<=epsdb).and.(nbc>12).or.(isnan(layerrap)) ) - - !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. - !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 ) - allocate( bc(nbc,nbc), indx(nbc) ) - allocate( f(nfext,nbc), sourcs(nfext) ) - - - 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) - - end do - end if - - !-- Saving Love numbers - - - hs(fr,deg+1) = loveh - ks(fr,deg+1) = lovek - 1.d0 - ls(fr,deg+1) = lovel - - !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,real(hs(fr,deg+1)), aimag(hs(fr,deg+1)),real(ks(fr,deg+1)), aimag(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 - -end subroutine