Index: sm/trunk-jpl/src/c/modules/FourierLoveCorex/FourierLoveCorex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/FourierLoveCorex/FourierLoveCorex.cpp (revision 22057) +++ (revision ) @@ -1,33 +1,0 @@ -/*!\file FourierLoveCorex - * \brief: core solution module for love analysis. Based on Caron's fortran core. - */ - -#include "./FourierLoveCorex.h" - -#include "../../classes/classes.h" -#include "../../shared/shared.h" -#include "../../toolkits/toolkits.h" -#include "../InputUpdateFromConstantx/InputUpdateFromConstantx.h" - -extern "C" { - int lnb_setup_( IssmDouble* LoveKr, IssmDouble* LoveKi, IssmDouble* LoveHr, IssmDouble* LoveHi, IssmDouble* LoveLr, - IssmDouble* LoveLi, int* pnfreq, IssmDouble* frequencies, int* psh_nmax, int* psh_nmin, - IssmDouble* pg0, IssmDouble* pr0, IssmDouble* pmu0,bool* pallow_layer_deletion, int* pforcing_type, bool* pverbosemod, - int* pnumlayers, IssmDouble* radius, IssmDouble* viscosity, IssmDouble* lame_lambda, IssmDouble* lame_mu, - IssmDouble* burgers_viscosity, IssmDouble* burgers_mu, IssmDouble* density, IssmDouble* isburgers, IssmDouble* issolid); - -} - -void FourierLoveCorex( IssmDouble* LoveKr, IssmDouble* LoveKi, IssmDouble* LoveHr, IssmDouble* LoveHi, IssmDouble* LoveLr, - IssmDouble* LoveLi, int nfreq, IssmDouble* frequencies, int sh_nmax, int sh_nmin, - IssmDouble g0, IssmDouble r0, IssmDouble mu0,bool allow_layer_deletion, int forcing_type, bool verbosemod, - int numlayers, IssmDouble* radius, IssmDouble* viscosity, IssmDouble* lame_lambda, IssmDouble* lame_mu, - IssmDouble* burgers_viscosity, IssmDouble* burgers_mu, IssmDouble* density, IssmDouble* isburgers, IssmDouble* issolid){ - - /*Call fortran driver: */ - lnb_setup_( LoveKr, LoveKi, LoveHr, LoveHi, LoveLr, - LoveLi, &nfreq, frequencies, &sh_nmax, &sh_nmin, - &g0, &r0, &mu0, &allow_layer_deletion, &forcing_type, &verbosemod, - &numlayers, radius, viscosity, lame_lambda, lame_mu, - burgers_viscosity, burgers_mu, density, isburgers, issolid); -} Index: sm/trunk-jpl/src/c/modules/FourierLoveCorex/FourierLoveCorex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/FourierLoveCorex/FourierLoveCorex.h (revision 22057) +++ (revision ) @@ -1,17 +1,0 @@ -/*!\file: FourierLoveCorex.h - * \brief header file for ... - */ - -#ifndef _FOURIERLOVECOREX_H -#define _FOURIERLOVECOREX_H - -#include "../../classes/classes.h" - -/* local prototypes: */ -void FourierLoveCorex( IssmDouble* LoveKr, IssmDouble* LoveKi, IssmDouble* LoveHr, IssmDouble* LoveHi, IssmDouble* LoveLr, - IssmDouble* LoveLi,int nfreq, IssmDouble* frequencies, int sh_nmax, int sh_nmin, IssmDouble g0, IssmDouble r0, - IssmDouble mu0,bool allow_layer_deletion, int forcing_type, bool verbosemod, int numlayers, IssmDouble* radius, - IssmDouble* viscosity,IssmDouble* lame_lambda, IssmDouble* lame_mu, IssmDouble* burgers_viscosity, IssmDouble* burgers_mu, - IssmDouble* density,IssmDouble* isburgers, IssmDouble* issolid); - -#endif /* _FOURIERLOVECOREX_H */ Index: sm/trunk-jpl/src/c/modules/FourierLoveCorex/lnb_param.f90 =================================================================== --- /issm/trunk-jpl/src/c/modules/FourierLoveCorex/lnb_param.f90 (revision 22057) +++ (revision ) @@ -1,35 +1,0 @@ -module lnb_param - - implicit none - - public - - !Physical Constants - double precision, parameter :: GG = 6.67259e-11! gravitational constant - double precision :: ra,rb,rc ! radii of the surface, CMB and ICB - double precision :: pi,ro_mean,go_surf,go_icb,go_cmb !some constants - ! GG --> gravitational constant - ! ra, rb, rc --> Surface, CMB and ICB radius - ! ro_mean --> mean density of the planet - ! go_surf --> mean surface gravity - ! go_icb/cmb --> gravity at the Inner Core/Core-Mantle Boundary - double precision :: mu0,go0,r0 !3 normalization constants - -!Numerical Constants - double precision :: epsdb, layerrap !used in detecting underflows in the lovenb integration - double complex :: valini,freq ! integration constant and angular frequency being processed - logical :: allow_layer_del,display, benchmark_spada !do we allow deletion of layers in case of underflow, enable some prints in terminal - -!GIA problem setup - integer :: deg,ifmin,ifmax ! deg=SH degree, ifmin/ifmax=which kind of forcing is applied (usually we want loading) - logical :: incompressible !choose which mode of calculation is activated - -!Earth Model - integer :: nfreq, degmax, degmin - integer :: nlayer! number of layers - double precision, dimension(:), pointer :: roc,lac,muc,vic,vic2,muc2 ! Earth layers parameters - double precision, dimension(:), pointer :: radius ! Earth layers parameters - double precision :: vicratio,mucratio !2 Burgers long-term to short-term ratios - logical, dimension(:), pointer :: soliddim,radbc,burgers ! is this layer solid (if not it's liquid), does it have specific boundary conditions and does it have a Burgers rheology (if not, it uses Maxwell) - -end module lnb_param Index: sm/trunk-jpl/src/c/modules/FourierLoveCorex/love_numbers.f90 =================================================================== --- /issm/trunk-jpl/src/c/modules/FourierLoveCorex/love_numbers.f90 (revision 22057) +++ (revision ) @@ -1,364 +1,0 @@ - subroutine lnb_setup( ksr, ksi, hsr, hsi, lsr, lsi, dnfreq, frequencies, ddegmax, ddegmin, & - dgo0, dr0, dmu0, dallow_layer_del, source_type, ddisplay, dnlayer, dradius, & - dvic, dlac, dmuc, dvic2, dmuc2, droc, dburgers, dsoliddim) - use lnb_param - use util - implicit none - - integer :: source_type, fr, i, IERR - integer :: dnfreq, ddegmax, ddegmin - integer :: dnlayer - - double precision :: T,cst, fluid_tscale,g - double precision, dimension(dnfreq,ddegmax+1) :: hsr,hsi,lsr,lsi,ksr,ksi - double complex, dimension(dnfreq,ddegmax+1) :: hs,ls,ks - double precision, dimension(dnfreq) :: frequencies - - double precision :: dmu0,dgo0,dr0 - logical :: dallow_layer_del, ddisplay - double precision, dimension(dnlayer) :: droc,dlac,dmuc,dvic,dvic2,dmuc2 - double precision, dimension(dnlayer+1) :: dradius - double precision, dimension(dnlayer) :: dsoliddim,dburgers - - !assign local variables to public variables - mu0=dmu0; go0=dgo0; r0=dr0; allow_layer_del=dallow_layer_del; - nfreq=dnfreq;degmax=ddegmax;degmin=ddegmin;nlayer=dnlayer; - allocate(roc(nlayer),lac(nlayer),muc(nlayer),vic(nlayer),vic2(nlayer),muc2(nlayer)) - roc=droc; lac=dlac; muc=dmuc; vic=dvic; vic2=dvic2; muc2=dmuc2; - allocate(radius(nlayer+1),soliddim(nlayer),burgers(nlayer),radbc(nlayer+1)) - radius=dradius; - - display=ddisplay - - do i=1,nlayer - if (NINT(dsoliddim(i))==1) then - soliddim(i)=.true. - else - soliddim(i)=.false. - endif - if (NINT(dburgers(i))==1) then - burgers(i)=.true. - else - burgers(i)=.false. - endif - enddo - - pi=dacos(-1.d0) - - ra=radius(nlayer+1) - rb = 3480.d3 - rc = 1221.5d3 - r0=ra - cst = 365.25d0*24.d0*3600.d0*1000.d0 - - benchmark_spada=.false. - !frequencies=1E-5/cst*2*pi - -! hs(:,:)=complex(0.d0,0.d0) - hs(:,:)=cmplx(0.d0,0.d0) - ks=hs - ls=hs - - call love_numbers(frequencies,source_type,hs,ls,ks) - - hsr=dble(hs) - hsi=dimag(hs) - ksr=dble(ks) - ksi=dimag(ks) - lsr=dble(ls) - lsi=dimag(ls) - -end subroutine - -!==================== - subroutine love_numbers(frequencies,source_type,hs,ls,ks) -!==================== - - 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 - - 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 - - - ! 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 - - - if (display) then - print*,'model init' - endif - - - call model_init2(nbc) - - - if (display) then - print*, 'done' - endif - - 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' - - - if (display) 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 - if (maxval(frequencies).ne.0 .and. display) then - print*, 'Frequency ranging from ', minval(dabs(frequencies)), & - ' to ', maxval(frequencies),' s^-1' - endif - - if (frequencies(1)==0 .and. display) then - print*,'Elastic calculation' - 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=cmplx(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=cmplx(0.d0,1.d0)*(dble(fr-1)/T)*2.d0*pi - !else - ! freq=cmplx(0.d0,1.d0)*(-dble(nfreq+1-fr)/T)*2.d0*pi - !end if - - freq=cmplx(0.d0,frequencies(fr)) - !print*,fr,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)*cmplx(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) - !print*,'bc ok', layerrap - call external_forcing(deg,f,sourcs,nfext,nbc) - !print*,'forcing ok', layerrap - call solution(deg,bc,indx,f,sourcs,nfext,nbc,loveh,lovel,lovek,delta) - - !print*,'first_sol', layerrap, epsdb - - ! 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)) ) - ! print*,'trying to delete layer', layerrap, epsdb, nbc - !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 - - !print*,dble(hs(fr,deg+1)), dimag(hs(fr,deg+1)),dble(ks(fr,deg+1)), dimag(ks(fr,deg+1)) - - !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 - - deallocate( bc, indx, f, sourcs ) -end subroutine Index: sm/trunk-jpl/src/c/modules/FourierLoveCorex/lovenb_sub.f90 =================================================================== --- /issm/trunk-jpl/src/c/modules/FourierLoveCorex/lovenb_sub.f90 (revision 22057) +++ (revision ) @@ -1,609 +1,0 @@ -module lovenb_sub - - use util - use lnb_param - - implicit none - - contains - - subroutine boundary_conditions_matrix(bc,indx,nbc) -!------------------------------------- - - implicit none - double complex :: la, mu - double precision :: xmin,xmax,eps,hini,hmin - double precision :: x,ro1,ro2,g,d,one - double complex, dimension (:,:), pointer :: bc - integer, dimension(:), pointer :: indx - double complex, dimension(6) :: ystart,dydx - double complex, dimension(6,6) :: ystep - integer :: ny,i,j,k,ii,jj,kk,nok,nbad,is,nbc,ibc,ifirst,nstep - logical :: solid,vali - - !-- Initialisation of boundary condition matrix - - indx(:) = 0 - - !-- Integration - - eps= 1.d-10 - hini=1.d-10 - hmin=0.d0 - - ! Underflows risk tests in the numerical integration - ! The starting ystart are chosen to be different than 1 to limit under and overflows in the Runge-Kutta - if (-300.d0 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 - - - !301 format(a4,5x,'n=',i3,5x,'h=',f14.10,5x,'l=',f14.10,5x,'k=',f14.10) - 301 format(a4,5x,'n=',i5,5x,'h=',d14.7,5x,'l=',d14.7,5x,'k=',d14.7) - 305 format(a4,i1,4x,'n=',i5,5x,'h=',d14.7,5x,'l=',d14.7,5x,'k=',d14.7) - 306 format(a4,i2,3x,'n=',i5,5x,'h=',d14.7,5x,'l=',d14.7,5x,'k=',d14.7) - 302 format(i5,5x,f14.10,5x,f14.10,5x,f14.10) - 303 format(i5,5x,f14.10,5x,f14.10,5x,f14.10,5x,f14.10) - 304 format(i5,5x,d20.10,5x,d20.10,5x,d20.10) - - end subroutine solution - - - subroutine solution_nbloutput(n,bc,indx,f,sourcs,nfext,nbc,loveh,lovel,lovek,delta) -!------------------------------------------------- - - integer :: n,nfext,nbc,i,j,ic,ibc,ii,jj,icmin,ifirst - double precision, dimension(:,:), pointer :: bc,f - integer, dimension(:), pointer :: indx - character*40, dimension(:), pointer :: sourcs - double precision, dimension(nbc) :: sc - double precision :: loveh,lovel,lovek,loveh1,lovel1,lovek1,loveh1s,lovel1s,lovek1s - double precision :: sumh,suml,sumk,d,val,delta,sumy2,sumy4,sumy6,rads - double precision, dimension(6) :: icbh,icbl,icbk,y2srf,y4srf,y6srf - double precision, dimension(nbc,nbc) :: bcsav - logical :: ok - - - !do ic=1,6 -! if (nbc==18) then -! icbh(ic) = bc(1,ic) -! icbl(ic) = bc(3,ic) -! icbk(ic) = bc(5,ic) -!! elseif (nbc==24) then -! icbh(ic) = bc(7,ic+3) -! icbl(ic) = bc(9,ic+3) -! icbk(ic) = bc(11,ic+3) -! end if -! y2srf(ic) = bc(nbc-4,ic+nbc-9) -! y4srf(ic) = bc(nbc-2,ic+nbc-9) -! y6srf(ic) = bc(nbc,ic+nbc-9) - ! end do - - bcsav(:,:) = bc(:,:) - rads=0 - ifirst = 0 - i = 0 - do while (ifirst==0) - i = i+1 - if (radbc(i)) ifirst = i - end do - !do i = 1,nbc -! do j = 1,nbc -! if (bc(i,j)/=0.d0) print*, i,j,bc(i,j) -! enddo - !enddo - !write(*,*) - - !-- LU decomposition and determinant - ok = .true. - - call ludcmp(bc,nbc,nbc,indx,d) - do i=1,nbc - d=d*bc(i,i) - end do - write(*,*) - !write(*,*),dlog(freq),' Matrix det = ',d/dabs(d)*dlog(dabs(d)) - -!-- Resolution - - if (ok) then - do i = ifmin,ifmax - - do j=1,nbc - sc(j)=f(i,j) - end do - - write(*,*) 'Source: ',sourcs(i) - call lubksb(bc,nbc,nbc,indx,sc) - - loveh = sc(nbc-2)*ra*go_surf - lovel = sc(nbc-1)*ra*go_surf - lovek = sc(nbc)*ra*go0 - delta = 1.d0-dble(n+1)/dble(n)*(lovek-1.d0)+2.d0/dble(n)*loveh - - sumy2=0.d0 - sumy4=0.d0 - sumy6=0.d0 - - do ic=1,6 - sumy2 = sumy2+bcsav(nbc-4,ic+nbc-9)*sc(ic+nbc-9) - sumy4 = sumy4+bcsav(nbc-2,ic+nbc-9)*sc(ic+nbc-9) - sumy6 = sumy6+bcsav(nbc,ic+nbc-9)*sc(ic+nbc-9) - end do - - !write(*,*) 'SURF',n,loveh,lovel,lovek-1.d0 - - 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)*ra*go_surf*valini - lovel1 = sc(nbc -ibc*6 -3 +3)*ra*go_surf*valini - lovek1 = sc(nbc -ibc*6 -3 +5)*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 - !print*, ii+1,jj+ic - sumh = sumh+bcsav(ii+1,jj+ic)*sc(jj+ic) - suml = suml+bcsav(ii+3,jj+ic)*sc(jj+ic) - sumk = sumk+bcsav(ii+5,jj+ic)*sc(jj+ic) - 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 = dabs(loveh1s/loveh) - if (layerrap > dabs(lovel1s/lovel)) layerrap = dabs(lovel1s/lovel) - if (layerrap > dabs((lovek1s-(rads/ra)**deg)/(lovek-1.d0))) & - layerrap = dabs((lovek1s-(rads/ra)**deg)/(lovek-1.d0)) - endif - - loveh1s = loveh1 - lovel1s = lovel1 - lovek1s = lovek1 - rads = radius(j) - !if (radius(j)==rb) then - ! write(*,301) 'CMB ',n,loveh1,lovel1,lovek1-(rb/ra)**deg - !else if (radius(j)==rc) then - ! write(*,301) 'ICB ',n,loveh1,lovel1,lovek1-(rc/ra)**deg - !else - ! if (j<10) then - ! write(*,305) 'ITF-',j,n,loveh1,lovel1,lovek1-(radius(j)/ra)**deg - ! else - ! write(*,306) 'ITF-',j,n,loveh1,lovel1,lovek1-(radius(j)/ra)**deg - ! end if - !endif - end if - end do - - end do - end if - - !write(*,*),dlog10(aimag(freq)*365.d0*24.d0*3600.d0*1000.d0),' Matrix det = ',d, lovek - - - - !301 format(a4,5x,'n=',i3,5x,'h=',f14.10,5x,'l=',f14.10,5x,'k=',f14.10) - 301 format(a4,5x,'n=',i5,5x,'h=',d14.7,5x,'l=',d14.7,5x,'k=',d14.7) - 305 format(a4,i1,4x,'n=',i5,5x,'h=',d14.7,5x,'l=',d14.7,5x,'k=',d14.7) - 306 format(a4,i2,3x,'n=',i5,5x,'h=',d14.7,5x,'l=',d14.7,5x,'k=',d14.7) - 302 format(i5,5x,f14.10,5x,f14.10,5x,f14.10) - 303 format(i5,5x,f14.10,5x,f14.10,5x,f14.10,5x,f14.10) - 304 format(i5,5x,d20.10,5x,d20.10,5x,d20.10) - - end subroutine solution_nbloutput - - - -end module lovenb_sub Index: sm/trunk-jpl/src/c/modules/FourierLoveCorex/model.f90 =================================================================== --- /issm/trunk-jpl/src/c/modules/FourierLoveCorex/model.f90 (revision 22057) +++ (revision ) @@ -1,422 +1,0 @@ -module model - - implicit none - - contains - - subroutine model_init2(nbc) -!--------------------- - use lnb_param - - double precision :: ro,g,drad,la,mu,ro1,ro2,x - logical :: solid,test - integer :: nbc,i,n,j - - - - !- Calculus parameters - pi=dacos(-1.d0) - epsdb=epsilon(1.d0) - - !write(101,*), 'Epsilon machine = ',epsdb - !write(101,*), 'Constants:' - !write(101,*), ' pi = ',pi - !write(101,*), ' G = ',GG - !write(101,*), 'Dimensioned scales:' - !write(101,*), ' rigidity mu0 = ',mu0 - !write(101,*), ' gravity go0 = ',go0 - !write(101,*), ' distance r0 = ',r0 - - - - - !call earth_homogeneous_layers_prembased(radius,nlayer,roc,lac,muc,g,solid) !yields the Earth layers radii, density, lamé parameters, gravity and solid/liquid state boolean identifier - - !if (benchmark_spada) then - ! muc(1) = 1.d6; muc(2) = 0.d0; muc(3) = 2.28340d11 - ! muc(4) = 1.05490d11; muc(5) = 0.70363d11; muc(6) = 0.50605d11 - ! roc(1) = 10751.d0; roc(2) = 10750.d0; roc(3) = 4978.d0 - ! roc(6)=3300.d0 - !else - !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 - !endif - - !vic2=vic/vicratio - !muc2=muc/mucratio - - if (display) then - print*, ' Layer top boundary rho la mu eta' - do i = 1,nlayer - ! print*, radius(i+1),roc(i),lac(i),muc(i),vic(i) - if (burgers(i)) then - print*, radius(i+1),roc(i),lac(i),muc(i),vic(i),muc2(i),vic2(i), 'B' - else - print*, radius(i+1),roc(i),lac(i),muc(i),vic(i) - endif - end do - - endif - !write(101,*), - - - call earth_nlayers(1.d0,ro,la,mu,g,solid) - go_surf = g - ro_mean = 3.d0/4.d0*go_surf/(pi*GG*ra) - if (display) then - print*, 'Surface Gravity: ', go_surf, ' m.s^-2, mean density: ', ro_mean, ' kg.m^-3' - print*, 'Center of the Earth approximated at r=', radius(1)/1e3, ' km' - endif - - ! Test number of interfaces and of boundary conditions - n = 0 - do i=1,nlayer+1 - if (radbc(i)) n=n+1 - enddo - if (nbc/=(6*n)) then - !write(101,*), - print*, 'Error nbc should be ',6*n,' instead of ',nbc - stop - endif - - ! Test solid-liquid boundaries - test = .false. - do i = 1,nlayer - if (radbc(i)) test = .true. - if (test) then - if ((.not.soliddim(i)).and.((.not.radbc(i)).or.(.not.radbc(i+1)))) then - print*, 'Error in model definition: one fluid layer is not set with correct boundary conditions' - stop - endif - endif - enddo - - !-- Display a few quantities - !write(101,*), ' surface radius = ',ra - !write(101,*), ' CMB radius = ',rb - !write(101,*), ' ICB radius = ',rc - !write(101,*), ' go surface = ',go_surf - !write(101,*), ' mean density = ',ro_mean - - !-- Diplsay ICB conditions - !x=rc/ra - !write(101,*), - !call earth_nlayers_viscoelas(x-epsdb,ro1,cmplx(la,0.d0),cmplx(mu,0d0),g,solid) - !write(101,*), 'ICB, g = ',g - !write(101,*), 'ICB, ro(Inner Core) = ', ro1 - !call earth_nlayers_viscoelas(x,ro2,cmplx(la,0.d0),cmplx(mu,0d0),g,solid) - !write(101,*), 'ICB, g = ',g - !write(101,*), 'ICB, ro(Outer Core) = ', ro2 - !go_icb = g - - !-- Display CMB conditions - !x=rb/ra - !write(101,*), - !call earth_nlayers_viscoelas(x-epsdb,ro1,cmplx(la,0.d0),cmplx(mu,0d0),g,solid) - !write(101,*), 'CMB, g = ',g - !write(101,*), 'CMB, ro(Outer Core) = ', ro1 - !call earth_nlayers_viscoelas(x,ro2,cmplx(la,0.d0),cmplx(mu,0d0),g,solid) - !write(101,*), 'CMB, g = ',g - !write(101,*), 'CMB, ro(Mantle) = ', ro2 - !go_cmb = g - - end subroutine model_init2 - - - subroutine earth_homogeneous_layers_prembased(radin,nl,rot,lat,mut,GR,solid) - use lnb_param - - double precision :: x0,GR - logical :: solid - integer :: i,j,nl - doubleprecision :: la,mu,ro,vp,vs,RT - doubleprecision :: t1,t2,t3,t4,r1,r2 - doubleprecision, dimension(14) :: r - doubleprecision, dimension(13,4) :: d,p,s - doubleprecision, dimension(:), pointer :: radin - doubleprecision, dimension(:), pointer :: rot,lat,mut - doubleprecision, dimension(nl+1) :: rad - - - RT = ra - r(1) = 0.d0; r(2) = 1221.5d0; r(3) = 3480.d0; r(4) = 3630.d0 - r(5) = 5600.d0; r(6) = 5701.d0; r(7) = 5771.d0; r(8) = 5971.d0 - r(9) = 6151.d0; r(10) = 6291.d0; r(11) = 6346.6d0 - r(12) = 6356.d0; r(13) = 6368.d0; r(14) = ra - -! if ((r(14)*1.d3) /= ra) stop'Problem in prem: Earth radius /= ra !' -! - d(:,:) = 0.d0 - d(1,1) = 13.0885d0; d(1,3) = -8.8381d0 - d(2,1) = 12.5815d0; d(2,2) = -1.2638d0; d(2,3) = -3.6426d0; d(2,4) = -5.5281d0 - d(3,1) = 7.9565d0 ; d(3,2) = -6.4761; d(3,3) = 5.5283d0; d(3,4) = -3.0807d0 - d(4,1) = 7.9565d0 ; d(4,2) = -6.4761; d(4,3) = 5.5283d0; d(4,4) = -3.0807d0 - d(5,1) = 7.9565d0 ; d(5,2) = -6.4761; d(5,3) = 5.5283d0; d(5,4) = -3.0807d0 - d(6,1) = 5.3197d0 ; d(6,2) = -1.4836d0 - d(7,1) = 11.2494d0; d(7,2) = -8.0298d0 - d(8,1) = 7.1089d0 ; d(8,2) = -3.8045d0 - d(9,1) = 2.6910d0 ; d(9,2) = 0.6924d0 - d(10,1) = 2.6910d0; d(10,2) = 0.6924d0 - d(11,1) = 2.9d0 - d(12,1) = 2.6d0 - -! ocean - if (soliddim(13).eqv..false.) then - d(13,1) = 1.02d0 - -! continental - else - d(13,1) = d(12,1) - end if - - p(:,:) = 0.d0 - p(1,1) = 11.2622d0 ; p(1,3) = -6.3640d0 - p(2,1) = 11.0487d0 ; p(2,2) = -4.0362d0; p(2,3) = 4.8023d0; p(2,4) = -13.5732d0 - p(3,1) = 15.3891d0 ; p(3,2) = -5.3181d0; p(3,3) = 5.5242d0; p(3,4) = -2.5514d0 - p(4,1) = 24.952d0 ; p(4,2) = -40.4673d0; p(4,3) = 51.4832d0; p(4,4) = -26.6419d0 - p(5,1) = 29.2766d0 ; p(5,2) = -23.6027d0; p(5,3) = 5.5242d0; p(5,4) = -2.5514d0 - p(6,1) = 19.0957d0 ; p(6,2) = -9.8672d0 - p(7,1) = 39.7027d0 ; p(7,2) = -32.6166d0 - p(8,1) = 20.3926d0 ; p(8,2) = -12.2569d0 - p(9,1) = 4.1875d0 ; p(9,2) = 3.9382d0 - p(10,1) = 4.1875d0 ; p(10,2) = 3.9382d0 - p(11,1) = 6.8d0 - p(12,1) = 5.8d0 -! -! ocean - if (soliddim(13).eqv..false.) then - p(13,1) = 1.45d0 -! -! continental - else - p(13,1) = p(12,1) - end if -!---- -! - s(:,:) = 0.d0 -! - s(1,1) = 3.6678d0; s(1,3) = -4.4475d0 - - s(3,1) = 6.9254d0; s(3,2) = 1.4672d0; s(3,3) = -2.0834d0; s(3,4) = 0.9783d0 - s(4,1) = 11.1671d0; s(4,2) = -13.7818d0; s(4,3) = 17.4575d0; s(4,4) = -9.2777d0 - s(5,1) = 22.3459d0; s(5,2) = -17.2473d0; s(5,3) = -2.0834d0; s(5,4) = 0.9783d0 - s(6,1) = 9.9839d0; s(6,2) = -4.9324 - s(7,1) = 22.3512d0; s(7,2) = -18.5856d0 - s(8,1) = 8.9496d0; s(8,2) = -4.4597 - s(9,1) = 2.1519d0; s(9,2) = 2.3481d0 - s(10,1) = 2.1519d0; s(10,2) = 2.3481d0 - s(11,1) = 3.9d0 - s(12,1) = 3.2d0 -! -! ocean (please don't modify) - if (soliddim(13).eqv..false.) then -! -! continental - else - s(13,1) = s(12,1) - end if -! -! - r(:) = r(:)*1.d3 - - !- handling the first layer : central sphere - rad = radin - rad(1) = 0.d0 - - do j = 1,nl - - ro = 0.d0 - vp = 0.d0 - vs = 0.d0 - - do i = 1,13 - - r1 = 0.d0 - r2 = 0.d0 - if ((rad(j) > r(i)).and.(rad(j) <= r(i+1))) then - if (rad(j+1) <= r(i+1)) then - r2 = rad(j+1) - r1 = rad(j) - else - r2 = r(i+1) - r1 = rad(j) - end if - else if (rad(j) <= r(i)) then - if ((rad(j+1) > r(i)).and.(rad(j+1) <= r(i+1))) then - r2 = rad(j+1) - r1 = r(i) - else if (rad(j+1) > r(i+1)) then - r2 = r(i+1) - r1 = r(i) - end if - end if - - t1 = d(i,1)/3.d0 - t2 = d(i,2)/(Rt*4.d0) - t3 = d(i,3)/((Rt**2)*5.d0) - t4 = d(i,4)/((Rt**3)*6.d0) - ro = ro + t1*(r2**3) + t2*(r2**4) + t3*(r2**5) + t4*(r2**6) - & - ( t1*(r1**3) + t2*(r1**4) + t3*(r1**5) + t4*(r1**6) ) - - t1 = p(i,1)/3.d0 - t2 = p(i,2)/(Rt*4.d0) - t3 = p(i,3)/((Rt**2)*5.d0) - t4 = p(i,4)/((Rt**3)*6.d0) - vp = vp + t1*(r2**3) + t2*(r2**4) + t3*(r2**5) + t4*(r2**6) - & - ( t1*(r1**3) + t2*(r1**4) + t3*(r1**5) + t4*(r1**6) ) - - t1 = s(i,1)/3.d0 - t2 = s(i,2)/(Rt*4.d0) - t3 = s(i,3)/((Rt**2)*5.d0) - t4 = s(i,4)/((Rt**3)*6.d0) - vs = vs + t1*(r2**3) + t2*(r2**4) + t3*(r2**5) + t4*(r2**6) - & - ( t1*(r1**3) + t2*(r1**4) + t3*(r1**5) + t4*(r1**6) ) - - end do - ro = ro*3 / (rad(j+1)**3-rad(j)**3) - vp = vp*3 /(rad(j+1)**3-rad(j)**3) - vs = vs*3 / (rad(j+1)**3-rad(j)**3) - mu = ro*vs**2 - la = ro*vp**2 - 2.d0*mu - ro = ro*1.d3 - la = la*1.d9 - mu = mu*1.d9 - - !print*, rad(j+1), ro, la , mu - rot(j) = ro - lat(j) = la - mut(j) = mu - end do - - - end subroutine earth_homogeneous_layers_prembased - - - subroutine earth_nlayers(x0,ro,la,mu,GR,solid) - use lnb_param -!--------------------------------------- -! Gets planet properties at non-dimentionalized radius x0: -!ro: density -!la: Lame lambda constant -!mu: Lame mu constant (shear modulus) -!GR: gravity -!solid: boolean, true if solid, otherwise liquid - integer :: i,j - doubleprecision :: la,mu,CST,GR - doubleprecision :: x,x0,xx0,ro - doubleprecision ::r1,r2 - logical :: solid - double precision, dimension(nlayer+1) :: r - - - - solid = .true. - - r(:) = radius(:) - r(1) = 0.d0 - x = x0 * ra - CST = 4.d0*pi*GG/3.d0 - if ( x > r(nlayer+1) ) x = r(nlayer+1) - - GR = 0.d0 - - do i = 1,nlayer - - if (x > r(i)) then - - r2 = r(i+1) - r1 = r(i) - ro = roc(i) - - if (x <= r(i+1)) then - - GR = GR + ro*(x**3-r1**3) - if (x > epsdb) GR = GR*CST/(x**2) - - else - GR = GR + ro*(r2**3-r1**3) - end if - end if - - if ((x >= r(i)).and.(x < r(i+1))) then - la = lac(i) - mu = muc(i) - ro = roc(i) - if ( soliddim(i).eqv..false. ) solid = .false. - - end if - - end do - - if (x>=r(nlayer+1)) then - la = lac(nlayer) - mu = muc(nlayer) - ro = roc(nlayer) - if ( soliddim(nlayer).eqv..false. ) solid = .false. - - end if - - end subroutine earth_nlayers - - subroutine earth_nlayers_viscoelas(x0,ro,la,mu,GR,solid) - use lnb_param -!--------------------------------------- -! Defines the earth rheology at non-dimensionalized radius x0 and angular frequency freq (global/public variable) - - integer :: i - double complex :: la, mu, cst - double precision :: GR,la0,mu00,ka,mu1,mu2,vi2 - double precision :: x,x0,xx0,ro,vi - logical :: notfound,solid,burg - - - - call earth_nlayers(x0,ro,la0,mu00,GR,solid) - - x = x0 * ra - do i = 1,nlayer - if ((x >= radius(i)).and.(x < radius(i+1))) then - vi = vic(i) - vi2 = vic2(i) - mu2=muc2(i) - burg=burgers(i) - end if - end do - - if (burg.eqv..true.) then - mu1=mu00 - - if ((vi*freq) /= 0.d0) then - ka=la0 + 2.d0/3.d0*mu00 - !print*,mu1,mu2,vi,vi2 - mu=mu1*freq*(freq+mu2/vi2)/((freq+mu2/vi2)*(freq+mu1/vi)+mu1/vi2*freq) - !la=4*la0*mu**2/(2*mu*(3*la0+2*mu)-3*la0*mu) - la=ka-2.d0/3.d0*mu - !print*,freq,mu,la - else - mu = mu00 - la = la0 - endif - else - ka = la0 + 2.d0/3.d0*mu00 - cst = vi*freq - if (cst /= 0.d0) then - mu = mu00/(1.d0+mu00/cst) - la = (la0 + mu00*ka/cst)/(1.d0 + mu00/cst) - else - mu = mu00 - la = la0 - endif - endif - - ! - !!write(101,*), freq,cst,la0,la,mu00,mu,vi - - end subroutine earth_nlayers_viscoelas - - - - -end module model Index: sm/trunk-jpl/src/c/modules/FourierLoveCorex/util.f90 =================================================================== --- /issm/trunk-jpl/src/c/modules/FourierLoveCorex/util.f90 (revision 22057) +++ (revision ) @@ -1,226 +1,0 @@ -module util - - use model - use lnb_param - - contains - -! -!_______DERIVS_________________________________________________________________ -! - subroutine derivs(x,y,dydx) - - implicit none - double complex :: la, mu - double precision :: x,ro,g - double complex :: flm,rlm,rm0,rg0,frh,fgr,fn - logical :: solid - integer :: iy,id,n2,ny - double complex, dimension(6) :: y,dydx,f - double complex, dimension(6,6) :: dy - - call earth_nlayers_viscoelas(x,ro,la,mu,g,solid) - !print*,la,mu -! -!_______Expressions -! - flm=(la+2.d0*mu) - rlm=(3.d0*la+2.d0*mu)/(la+2.d0*mu) - rm0=mu/mu0 - rg0=g/go0 - frh=ro*g*ra/mu0 - fgr=4.d0*Pi*GG*ro*ra/go0 - fn=dble(deg*(deg+1)) - -! -!_______Systeme differentiel yi dans un solide -! - if(solid) then - - ny = 6 - - dy(1,1)= (-2.d0*la/flm)/x - dy(1,2)= mu0/flm - dy(1,3)= (la*fn/flm)/x - dy(1,4)= 0.d0 - dy(1,5)= 0.d0 - dy(1,6)= 0.d0 - - dy(2,1)= 4.d0*(-frh+rm0*rlm/x)/x + ro*freq**2*ra/mu0 - dy(2,2)=(-4.d0*mu/flm)/x - dy(2,3)= fn*(frh-2.d0*rm0*rlm/x)/x - dy(2,4)= fn/x - dy(2,5)= 0.d0 - dy(2,6)= -frh/rg0 - - dy(3,1)= -1.d0/x - dy(3,2)= 0.d0 - dy(3,3)= 1.d0/x - dy(3,4)= 1.d0/rm0 - dy(3,5)= 0.d0 - dy(3,6)= 0.d0 - - dy(4,1)= (frh-2.d0*rm0*rlm/x)/x - dy(4,2)= ( -la/flm)/x - dy(4,3)= (2.d0*rm0*(la*(2.d0*fn-1.d0)+2.d0*mu*(fn-1.d0))/flm)/(x*x) + ro*freq**2*ra/mu0 - dy(4,4)= -3.d0/x - dy(4,5)= -(frh/rg0)/x - dy(4,6)= 0.d0 - - dy(5,1)= fgr - dy(5,2)= 0.d0 - dy(5,3)= 0.d0 - dy(5,4)= 0.d0 - dy(5,5)= 0.d0 - dy(5,6)= 1.d0 - - dy(6,1)= 0.d0 - dy(6,2)= 0.d0 - dy(6,3)= -(fgr*fn)/x - dy(6,4)= 0.d0 - dy(6,5)= fn/(x*x) - dy(6,6)= -2.d0/x - - else - - ny = 2 - - dy(:,:)= 0.d0 - - dy(1,1)= fgr/rg0 - dy(1,2)= 1.d0 - dy(2,1)= (-4.d0*(fgr/rg0)+fn/x)/x - dy(2,2)= -2.d0/x-fgr/rg0 - - end if -! -!_______Derivees des yi -! - dydx(:)=0.d0 - do id=1,ny - do iy=1,ny - dydx(id)=dydx(id)+dy(id,iy)*y(iy) - end do - end do - - end subroutine derivs - - - subroutine propagFDS(ystart,xmin,xmax,nstep) - implicit none - double precision :: xmin,xmax,dr,x - integer :: nstep,i - double complex, dimension(6) :: ystart,dydx,y,y1 - - dr = (xmax -xmin)/nstep - y = ystart - x = xmin - - do i = 1,nstep - - call derivs(x,y,dydx) - y1 = y + dr*dydx - x = x + dr - y = y1 - - end do - - ystart = y - - end subroutine propagFDS - -!_____LUDCMP___________________________________________________________________ -! - SUBROUTINE ludcmp(a,n,np,indx,d) - IMPLICIT DOUBLE PRECISION (a-h,o-z) - INTEGER n,np,indx(n),NMAX - REAL*8 d,a(np,np),TINY - PARAMETER (NMAX=500,TINY=1.0D-30) - INTEGER i,imax,j,k - REAL*8 aamax,dum,sum,vv(NMAX) - d=1d0 - do 12 i=1,n - aamax=0d0 - do 11 j=1,n - if (dabs(a(i,j)).gt.aamax) aamax=dabs(a(i,j)) -11 continue - if (aamax.eq.0d0) stop 'singular matrix in ludcmp' - vv(i)=1./aamax -12 continue - do 19 j=1,n - do 14 i=1,j-1 - sum=a(i,j) - do 13 k=1,i-1 - sum=sum-a(i,k)*a(k,j) -13 continue - a(i,j)=sum -14 continue - aamax=0d0 - do 16 i=j,n - sum=a(i,j) - do 15 k=1,j-1 - sum=sum-a(i,k)*a(k,j) -15 continue - a(i,j)=sum - dum=vv(i)*dabs(sum) - if (dum.ge.aamax) then - imax=i - aamax=dum - endif -16 continue - if (j.ne.imax)then - do 17 k=1,n - dum=a(imax,k) - a(imax,k)=a(j,k) - a(j,k)=dum -17 continue - d=-d - vv(imax)=vv(j) - endif - indx(j)=imax - if(a(j,j).eq.0d0)a(j,j)=TINY - if(j.ne.n)then - dum=1d0/a(j,j) - do 18 i=j+1,n - a(i,j)=a(i,j)*dum -18 continue - endif -19 continue - return - END SUBROUTINE ludcmp -! -!_____LUBKSB___________________________________________________________________ -! - SUBROUTINE lubksb(a,n,np,indx,b) - IMPLICIT DOUBLE PRECISION (a-h,o-z) - INTEGER n,np,indx(n) - REAL*8 a(np,np),b(np) - INTEGER i,ii,j,ll - REAL*8 sum - ii=0 - do 12 i=1,n - ll=indx(i) - sum=b(ll) - b(ll)=b(i) - if (ii.ne.0)then - do 11 j=ii,i-1 - sum=sum-a(i,j)*b(j) -11 continue - else if (sum.ne.0d0) then - ii=i - endif - b(i)=sum -12 continue - do 14 i=n,1,-1 - sum=b(i) - do 13 j=i+1,n - sum=sum-a(i,j)*b(j) -13 continue - b(i)=sum/a(i,i) -14 continue - return - END SUBROUTINE lubksb - - - -end module util