[cig-commits] r5266 - in geodyn/3D/MoSST/trunk: . mosst-data
mosst-data/fftw_v10 mosst-dynamo
wei at geodynamics.org
wei at geodynamics.org
Tue Nov 14 15:34:51 PST 2006
Author: wei
Date: 2006-11-14 15:34:50 -0800 (Tue, 14 Nov 2006)
New Revision: 5266
Added:
geodyn/3D/MoSST/trunk/mosst-data/
geodyn/3D/MoSST/trunk/mosst-data/fftw_v10/
geodyn/3D/MoSST/trunk/mosst-data/fftw_v10/sphere.fftw_v10.1000
geodyn/3D/MoSST/trunk/mosst-dynamo/
geodyn/3D/MoSST/trunk/mosst-dynamo/make_modules
geodyn/3D/MoSST/trunk/mosst-dynamo/make_v10
geodyn/3D/MoSST/trunk/mosst-dynamo/mod_geoms.f
geodyn/3D/MoSST/trunk/mosst-lib/
Removed:
geodyn/3D/MoSST/trunk/lib/
geodyn/3D/MoSST/trunk/src/
Log:
Added directory mosst-data; added new module: mod_geoms.f; added make scripts: make_modules and make_v10; rearranged file structrue for the trunk. MoSST compiles.
Added: geodyn/3D/MoSST/trunk/mosst-data/fftw_v10/sphere.fftw_v10.1000
===================================================================
(Binary files differ)
Property changes on: geodyn/3D/MoSST/trunk/mosst-data/fftw_v10/sphere.fftw_v10.1000
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Copied: geodyn/3D/MoSST/trunk/mosst-dynamo (from rev 5193, geodyn/3D/MoSST/trunk/src)
Added: geodyn/3D/MoSST/trunk/mosst-dynamo/make_modules
===================================================================
--- geodyn/3D/MoSST/trunk/src/make_modules 2006-11-07 19:33:31 UTC (rev 5193)
+++ geodyn/3D/MoSST/trunk/mosst-dynamo/make_modules 2006-11-14 23:34:50 UTC (rev 5266)
@@ -0,0 +1,14 @@
+ex mod_params.f << EOF
+1,\$s/Lmax_v=15,mmax_v=10/Lmax_v=33,mmax_v=21/
+1,\$s/Lmax_m=15,mmax_m=10/Lmax_m=33,mmax_m=21/
+1,\$s/Lmax_t=15,mmax_t=10/Lmax_t=33,mmax_t=21/
+1,\$s/miner=1/miner=2/
+1,\$s/nmaxo=31,nmaxi=20,nmaxm=20/nmaxo=39,nmaxi=35,nmaxm=19/
+1,\$s/^cdia/ /
+1,\$s/^cdiv/ /
+wq
+EOF
+
+ifort -r8 -static -c mod_params.f mod_geoms.f mod_fields.f mod_artdis.f mod_anomaly.f mod_dataio.f mod_matrices.f -lintelsub90d -lblas -L$HOME/lib -L/usr/lib -nofor_main -O2
+
+
Property changes on: geodyn/3D/MoSST/trunk/mosst-dynamo/make_modules
___________________________________________________________________
Name: svn:executable
+ *
Added: geodyn/3D/MoSST/trunk/mosst-dynamo/make_v10
===================================================================
--- geodyn/3D/MoSST/trunk/src/make_v10 2006-11-07 19:33:31 UTC (rev 5193)
+++ geodyn/3D/MoSST/trunk/mosst-dynamo/make_v10 2006-11-14 23:34:50 UTC (rev 5266)
@@ -0,0 +1,8 @@
+ex evolutions.f << EOF
+1,\$s/^cts1/ /
+wq
+EOF
+
+
+ifort -r8 -static -o mode10 mod_artdis.o mod_fields.o mod_matrices.o mod_anomaly.o mod_dataio.o mod_params.o mod_geoms.o mosst_cig.f evolutions.f params_io.f matrices.f miscs.f bcs.f solvers.f time_integ.f forces.f -lintelsub90d -lfftw3 -lblas -lm -L$HOME/lib -L/opt/fftw3.0.1/lib -L/usr/lib
+
Property changes on: geodyn/3D/MoSST/trunk/mosst-dynamo/make_v10
___________________________________________________________________
Name: svn:executable
+ *
Added: geodyn/3D/MoSST/trunk/mosst-dynamo/mod_geoms.f
===================================================================
--- geodyn/3D/MoSST/trunk/src/mod_geoms.f 2006-11-07 19:33:31 UTC (rev 5193)
+++ geodyn/3D/MoSST/trunk/mosst-dynamo/mod_geoms.f 2006-11-14 23:34:50 UTC (rev 5266)
@@ -0,0 +1,788 @@
+!
+! This group of the moduels defines the geomatric properties of the
+! model, including spatial geometry and parity in spectral representation
+!
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+c
+c This module defines the Chebyshev polynomials; the radial collocation
+c points in the inner core, outer core and D"-layer; the radially
+c finite difference formular coefficients.
+c
+c The subroutines contained in this module assign the values to the
+c defined vectors. They can only be executed after the parameters
+c in MOD_SYSPARAM are assigned values.
+c
+c CHEB: The chebyschev polynomials T_n;
+c DCH1: The first order derivative of T_n;
+c DCH2: The second order derivative of T_n;
+c DCH3: The third order derivative of T_n;
+c DCH4: The fourth order derivative of T_n;
+c
+c RR: The radial collocation points in the outer core;
+c GG: The radial derivative coefficients;
+c DR1: The coefficients for first order differece in radius;
+c (in the outer core);
+c d/dr f(r_i) = dr1(i,1)*f_(i-2) + dr1(i,2)*f_(i-1)+
+c dr1(i,3)*f_i+dr1(i,4)*f_(i+1)+
+c dr1(i,5)*f_(i+2);
+c
+c RI: The radial collocation points in the inner core;
+c GI: The radial derivative coefficients in the inner core;
+c DRI1: Similar to DR1, but for the inner core;
+c
+c RD: The radial collocation points in the D"-layer;
+c GD: The radial derivative coefficients in the D"-layer;
+c DRD1: Similar to DR1, but for the D"-layer;
+c
+c DHR: finite difference in radius;
+c
+c CFM: The coefficients for the finite difference approximations
+c of the radial derivatives near the ICB and the CMB; they
+c are used for the linear matrices of the induction and
+c thermal equations;
+c
+c Author: Weijia Kuang
+c Date: Jan., 2001
+c COPY RIGHT: THIS CODE CAN NOT BE DUPLICATED OR USED WITHOUT
+c THE PERMISSION OF THE AUTHOR.
+c
+
+ MODULE mod_radgeom
+
+ use mod_dimparam
+ use mod_sysparam
+
+ implicit none
+
+ real (kind=8) aa,bb,hno,hni,hnm
+ real (kind=8) cfm(5,8)
+
+ real (kind=8), allocatable, dimension(:,:) :: cheb,dch1,dch2,
+ & dch3,dch4
+ real (kind=8), allocatable, dimension(:) :: rr,gg,ri,gi,rd,gd,dhr
+ real (kind=8), allocatable, dimension(:,:) :: dr1,dri1,drd1
+
+ CONTAINS
+
+c The subroutines that define the radial points and Chebyshev
+c polynomials
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ SUBROUTINE radgeom_init
+
+ if (.not. allocated(cheb)) then
+ allocate(cheb(0:nmaxo,0:nmaxo))
+ endif
+ if (.not. allocated(dch1)) then
+ allocate(dch1(0:nmaxo,0:nmaxo))
+ endif
+ if (.not. allocated(dch2)) then
+ allocate(dch2(0:nmaxo,0:nmaxo))
+ endif
+ if (.not. allocated(dch3)) then
+ allocate(dch3(0:nmaxo,0:nmaxo))
+ endif
+ if (.not. allocated(dch4)) then
+ allocate(dch4(0:nmaxo,0:nmaxo))
+ endif
+
+ if (.not. allocated(rr)) then
+ allocate(rr(0:nmaxo))
+ endif
+ if (.not. allocated(gg)) then
+ allocate(gg(0:nmaxo))
+ endif
+ if (.not. allocated(dr1)) then
+ allocate(dr1(0:nmaxo,5))
+ endif
+ if (.not. allocated(ri)) then
+ allocate(ri(0:nmaxi))
+ endif
+ if (.not. allocated(gi)) then
+ allocate(gi(0:nmaxi))
+ endif
+ if (.not. allocated(dri1)) then
+ allocate(dri1(0:nmaxi,5))
+ endif
+ if (.not. allocated(rd)) then
+ allocate(rd(0:nmaxm))
+ endif
+ if (.not. allocated(gd)) then
+ allocate(gd(0:nmaxm))
+ endif
+ if (.not. allocated(drd1)) then
+ allocate(drd1(0:nmaxm,5))
+ endif
+ if (.not. allocated(dhr)) then
+ allocate(dhr(nmxo1))
+ endif
+
+ END SUBROUTINE radgeom_init
+
+***********************************************************************
+***********************************************************************
+
+ SUBROUTINE radgeom
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+c T_n(x_i)= cos(n th_i);
+c x_i = cos(th_i);
+c r_i = a x_i + b;
+c th_i = (N-i) PI/N, for i = 0,1,...,N;
+c
+c cheb(n,i) = T_n(x_i);
+c dch1(n,i) = d/dr T_n(x_i);
+c dch2(n,i) = d^2/dr^2 T_n(x_i);
+c dch3(n,i) = d^3/dr^3 T_n(x_i);
+c dch4(n,i) = d^4/dr^4 T_n(x_i);
+c
+c d/dr f_i = dr1(i,1) f_{i-2} + dr1(i,2) f_{i-1} +
+c dr1(i,3) f_i + dr1(i,4) f_{i+1) + dr1(i,5) f_{i+2}
+c
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+ implicit none
+ integer i,k,n
+ real (kind=8) x,dps,c1,c2,c3,c4
+
+ call radgeom_init
+
+ n = nmaxo
+ aa = 0.5*(1.0-rio)
+ bb = 0.5*(1.0+rio)
+
+ rr = 0.0
+ gg = 0.0
+ dr1 = 0.0
+ ri = 0.0
+ gi = 0.0
+ dri1 = 0.0
+ rd = 0.0
+ gd = 0.0
+ drd1 = 0.0
+
+c
+c Chebyshev polynomials and radial collocation points in the outer
+c core
+c
+
+ do i = 0,n
+
+ cheb(:,i) = 0.0
+ dch1(:,i) = 0.0
+ dch2(:,i) = 0.0
+ dch3(:,i) = 0.0
+ dch4(:,i) = 0.0
+
+ x = cos((n-i)*pi/(1.0*n))
+ rr(i)= aa*x+bb
+ gg(i)= -aa*sin((n-i)*pi/(1.0*n))
+
+ cheb(0,i) = 1.0
+ cheb(1,i) = x
+ do k = 2,n
+ cheb(k,i) = 2.0*x*cheb(k-1,i)-cheb(k-2,i)
+ enddo
+
+ dch1(0,i) = 0.0
+ dch1(1,i) = 1.0
+ do k = 2,n
+ dch1(k,i) = 2.0*x*dch1(k-1,i)-dch1(k-2,i)
+ & +2.0*cheb(k-1,i)
+ enddo
+
+ dch2(0,i) = 0.0
+ dch2(1,i) = 0.0
+ do k = 2,n
+ dch2(k,i) = 2.0*x*dch2(k-1,i)-dch2(k-2,i)
+ & +4.0*dch1(k-1,i)
+ enddo
+
+ dch3(0,i) = 0.0
+ dch3(1,i) = 0.0
+ do k = 2,n
+ dch3(k,i) = 2.0*x*dch3(k-1,i)-dch3(k-2,i)
+ & +6.0*dch2(k-1,i)
+ enddo
+
+ dch4(0,i) = 0.0
+ dch4(1,i) = 0.0
+ do k = 2,n
+ dch4(k,i) = 2.0*x*dch4(k-1,i)-dch4(k-2,i)
+ & +8.0*dch3(k-1,i)
+ enddo
+
+ enddo
+
+ c1= 1.0/aa
+ c2= (1.0/aa)**2
+ c3= (1.0/aa)**3
+ c4= (1.0/aa)**4
+
+ dch1 = c1*dch1
+ dch2 = c2*dch2
+ dch3 = c3*dch3
+ dch4 = c4*dch4
+
+c
+c Radial derivative coefficients in the outer core
+c
+
+ dps = -pi/n
+ hno = 1.0/dps
+ do i = 2,n-2
+ x = 1.0/(12.0*gg(i)*dps)
+ dr1(i,1) = x
+ dr1(i,2) = -8.0*x
+ dr1(i,4) = 8.0*x
+ dr1(i,5) = -x
+ enddo
+
+ i = 1
+ x = 1.0/(12.0*gg(i)*dps)
+ dr1(i,1)= -3.0*x
+ dr1(i,2)= -10.0*x
+ dr1(i,3)= 18.0*x
+ dr1(i,4)= -6.0*x
+ dr1(i,5)= x
+
+ i = n-1
+ x = 1.0/(12.0*gg(i)*dps)
+ dr1(i,1)= -x
+ dr1(i,2)= 6.0*x
+ dr1(i,3)= -18.0*x
+ dr1(i,4)= 10.0*x
+ dr1(i,5)= 3.0*x
+
+ c1= rr(1)-rr(0)
+ c2= rr(2)-rr(0)
+ c3= rr(3)-rr(0)
+ c4= rr(4)-rr(0)
+ dr1(0,5)= -c1*c2*c3/(c4*(c4-c1)*(c4-c2)*(c4-c3))
+ dr1(0,4)= -dr1(0,5)*c4*c4*(c4-c1)*(c4-c2)/(c3*c3
+ & *(c3-c1)*(c3-c2))
+ dr1(0,3)= -(dr1(0,4)*c3**3*(c3-c1)+dr1(0,5)*c4**3
+ & *(c4-c1))/(c2**3*(c2-c1))
+ dr1(0,2)= (1.0-dr1(0,3)*c2-dr1(0,4)*c3-dr1(0,5)*
+ & c4)/c1
+ dr1(0,1)= -(dr1(0,2)+dr1(0,3)+dr1(0,4)+dr1(0,5))
+
+ c1= rr(n-1)-rr(n)
+ c2= rr(n-2)-rr(n)
+ c3= rr(n-3)-rr(n)
+ c4= rr(n-4)-rr(n)
+ dr1(n,1)= -c1*c2*c3/(c4*(c4-c1)*(c4-c2)*(c4-c3))
+ dr1(n,2)= -dr1(n,1)*c4*c4*(c4-c1)*(c4-c2)/(c3*c3
+ & *(c3-c1)*(c3-c2))
+ dr1(n,3)= -(dr1(n,2)*c3**3*(c3-c1)+dr1(n,1)*c4**3
+ & *(c4-c1))/(c2**3*(c2-c1))
+ dr1(n,4)= (1.0-dr1(n,3)*c2-dr1(n,2)*c3-dr1(n,1)*
+ & c4)/c1
+ dr1(n,5)= -(dr1(n,4)+dr1(n,3)+dr1(n,2)+dr1(n,1))
+
+c
+c Radial collocation points and the radial derivative coefficients
+c in the inner core
+c
+
+ n = nmaxi
+ dps = 1.0/n
+ hni = 1.0/dps
+ c1 = (rio-rco)
+
+ do i = 0,n
+ x = i*dps
+ ri(i)= rio-c1*(1.0-x)**2
+ gi(i)= 2.0*c1*(1.0-x)
+ enddo
+
+ do i = 2,n-2
+ x = 1.0/(12.0*gi(i)*dps)
+ dri1(i,1) = x
+ dri1(i,2) = -8.0*x
+ dri1(i,4) = 8.0*x
+ dri1(i,5) = -x
+ enddo
+
+ i = 0
+ x = 1.0/(12.0*gi(i)*dps)
+ dri1(i,1)= -25.0*x
+ dri1(i,2)= 48.0*x
+ dri1(i,3)= -36.0*x
+ dri1(i,4)= 16.0*x
+ dri1(i,5)= -3.0*x
+
+ i = 1
+ x = 1.0/(12.0*gi(i)*dps)
+ dri1(i,1)= -3.0*x
+ dri1(i,2)= -10.0*x
+ dri1(i,3)= 18.0*x
+ dri1(i,4)= -6.0*x
+ dri1(i,5)= x
+
+ i = n-1
+ x = 1.0/(12.0*gi(i)*dps)
+ dri1(i,1)= -x
+ dri1(i,2)= 6.0*x
+ dri1(i,3)= -18.0*x
+ dri1(i,4)= 10.0*x
+ dri1(i,5)= 3.0*x
+
+ c1= ri(n-1)-ri(n)
+ c2= ri(n-2)-ri(n)
+ c3= ri(n-3)-ri(n)
+ c4= ri(n-4)-ri(n)
+ dri1(n,1)= -c1*c2*c3/(c4*(c4-c1)*(c4-c2)*(c4-c3))
+ dri1(n,2)= -dri1(n,1)*c4*c4*(c4-c1)*(c4-c2)/(c3*c3
+ & *(c3-c1)*(c3-c2))
+ dri1(n,3)= -(dri1(n,2)*c3**3*(c3-c1)+dri1(n,1)*c4**3
+ & *(c4-c1))/(c2**3*(c2-c1))
+ dri1(n,4)= (1.0-dri1(n,3)*c2-dri1(n,2)*c3-dri1(n,1)*
+ & c4)/c1
+ dri1(n,5)= -(dri1(n,4)+dri1(n,3)+dri1(n,2)+dri1(n,1))
+
+c
+c Radial collocation points and the radial derivative coefficients
+c in the D"-layer
+c
+
+ n = nmaxm
+ dps = 1.0/n
+ hnm = 1.0/dps
+ c1 = (rdo-1.0)
+ do i = 0,n
+ x = i*dps
+ rd(i)= 1.0+c1*x**2
+ gd(i)= 2.0*c1*x
+ enddo
+
+ do i = 2,n-2
+ x = 1.0/(12.0*gd(i)*dps)
+ drd1(i,1) = x
+ drd1(i,2) = -8.0*x
+ drd1(i,4) = 8.0*x
+ drd1(i,5) = -x
+ enddo
+
+ c1= rd(1)-rd(0)
+ c2= rd(2)-rd(0)
+ c3= rd(3)-rd(0)
+ c4= rd(4)-rd(0)
+ drd1(0,5)= -c1*c2*c3/(c4*(c4-c1)*(c4-c2)*(c4-c3))
+ drd1(0,4)= -drd1(0,5)*c4*c4*(c4-c1)*(c4-c2)/(c3*c3
+ & *(c3-c1)*(c3-c2))
+ drd1(0,3)= -(drd1(0,4)*c3**3*(c3-c1)+drd1(0,5)*c4**3
+ & *(c4-c1))/(c2**3*(c2-c1))
+ drd1(0,2)= (1.0-drd1(0,3)*c2-drd1(0,4)*c3-drd1(0,5)*
+ & c4)/c1
+ drd1(0,1)= -(drd1(0,2)+drd1(0,3)+drd1(0,4)+drd1(0,5))
+
+ i = 1
+ x = 1.0/(12.0*gd(i)*dps)
+ drd1(i,1)= -3.0*x
+ drd1(i,2)= -10.0*x
+ drd1(i,3)= 18.0*x
+ drd1(i,4)= -6.0*x
+ drd1(i,5)= x
+
+ i = n-1
+ x = 1.0/(12.0*gd(i)*dps)
+ drd1(i,1)= -x
+ drd1(i,2)= 6.0*x
+ drd1(i,3)= -18.0*x
+ drd1(i,4)= 10.0*x
+ drd1(i,5)= 3.0*x
+
+ i = n
+ x = 1.0/(12.0*gi(i)*dps)
+ dri1(i,1)= 3.0*x
+ dri1(i,2)= -16.0*x
+ dri1(i,3)= 36.0*x
+ dri1(i,4)= -48.0*x
+ dri1(i,5)= 25.0*x
+
+c
+c Radial mesh size in the outer core
+c
+
+ do k = 1,nmaxo-1
+ dhr(k+1) = min(abs(rr(k+1)-rr(k)),abs(rr(k)-rr(k-1)))
+ enddo
+ dhr(1) = abs(rr(1)-rr(0))
+ dhr(nmxo1)= abs(rr(nmaxo)-rr(nmaxo-1))
+
+c
+c The radial derivative coefficients near the CMB and the ICB
+c
+
+ i = nmaxi-2
+ cfm(5,1)= (ri(i+2)-ri(i))/(2*ri(i+2)-ri(i)-
+ & ri(i+1))
+ cfm(4,1)= 1.0-cfm(5,1)
+ cfm(3,1)= cfm(5,1)/(ri(i+2)-ri(i+1))*((ri(i+1)-
+ & ri(i))/(ri(i+2)-ri(i)))**2
+ cfm(2,1)= (1-cfm(3,1)*(ri(i+2)-ri(i)))/(ri(i+1)
+ & -ri(i))
+ cfm(1,1)= -(cfm(2,1)+cfm(3,1))
+
+ i = nmaxi-1
+ cfm(5,2)= 1.0/(2.0+(ri(i+1)-ri(i))/(ri(i)-ri(i-1)))
+ cfm(4,2)= 1.0-cfm(5,2)
+ cfm(3,2)= cfm(5,2)*(1.0/(ri(i+1)-ri(i-1))+
+ & 2.0/(ri(i+1)-ri(i)))
+ cfm(1,2)= (1-cfm(3,2)*(ri(i+1)-ri(i)))/
+ & (ri(i-1)-ri(i))
+ cfm(2,2)= -(cfm(1,2)+cfm(3,2))
+
+ do i = 0,1
+ cfm(5,3+i)= 1.0/(2.0-(rr(i+1)-rr(i))/(rr(i+2)-rr(i)))
+ cfm(4,3+i)= 1.0-cfm(5,3+i)
+ cfm(3,3+i)= cfm(5,3+i)/(rr(i+2)-rr(i+1))*((rr(i+1)-
+ & rr(i))/(rr(i+2)-rr(i)))**2
+ cfm(2,3+i)= (1-cfm(3,3+i)*(rr(i+2)-rr(i)))/(rr(i+1)-
+ & rr(i))
+ cfm(1,3+i)= -(cfm(2,3+i)+cfm(3,3+i))
+ cfm(5,7+i)= 1.0/(2.0-(rd(i+1)-rd(i))/(rd(i+2)-rd(i)))
+ cfm(4,7+i)= 1.0-cfm(5,7+i)
+ cfm(3,7+i)= cfm(5,7+i)/(rd(i+2)-rd(i+1))*((rd(i+1)-
+ & rd(i))/(rd(i+2)-rd(i)))**2
+ cfm(2,7+i)= (1-cfm(3,7+i)*(rd(i+2)-rd(i)))/(rd(i+1)-
+ & rd(i))
+ cfm(1,7+i)= -(cfm(2,7+i)+cfm(3,7+i))
+ enddo
+
+ i = nmaxo-2
+ cfm(5,5)= 1.0/(2.0-(rr(i+1)-rr(i))/(rr(i+2)-rr(i)))
+ cfm(4,5)= 1.0-cfm(5,5)
+ cfm(3,5)= cfm(5,5)/(rr(i+2)-rr(i+1))*((rr(i+1)-rr(i))/
+ & (rr(i+2)-rr(i)))**2
+ cfm(2,5)= (1.0-cfm(3,5)*(rr(i+2)-rr(i)))/(rr(i+1)-rr(i))
+ cfm(1,5)= -(cfm(2,5)+cfm(3,5))
+
+ i = nmaxo-1
+ cfm(5,6)= 1.0/(2.0+(rr(i+1)-rr(i))/(rr(i)-rr(i-1)))
+ cfm(4,6)= 1.0-cfm(5,6)
+ cfm(3,6)= cfm(5,6)*(1.0/(rr(i+1)-rr(i-1))+
+ & 2.0/(rr(i+1)-rr(i)))
+ cfm(1,6)= (1-cfm(3,6)*(rr(i+1)-rr(i)))/
+ & (rr(i-1)-rr(i))
+ cfm(2,6)= -(cfm(1,6)+cfm(3,6))
+
+ END SUBROUTINE radgeom
+
+********************************************************************
+********************************************************************
+
+ END MODULE mod_radgeom
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+c
+c This module defines the spherical collocation points; the
+c recurrence and horizontal derivative coefficients; the local
+c meridional mesh size (used for CFL condition);
+c
+c ASLG: The associate Legendre polynomials;
+c TH: The cosine at collocation points in colatitude;
+c SINS: The cosine at collocation points in colatitude;
+c GAUWT: Gaussian weight for spherical transform;
+! TABLE: array for FFT transforms:
+! for IMSL FFT: the array for initializxation
+! for FFTW:
+! TABLE(1): forward transform
+! TABLE(2): backward transform
+!
+c WTTFR: initializing arrays for IMSL FFT routines;
+c
+c CLM: the coefficients for the recurrence relations and
+c the derivatives of Y_lm:
+c sin(th) d/dth f = \sum f1^{Lm} Y_L^m;
+c [sin(th) d/dth + 2 cos(th)] f = \sum f2^{Lm} Y_L^m;
+c [sin(th) d/dth - cos(th)] f = \sum f3^{Lm} Y_L^m;
+c cos(th) f = \sum f4^{Lm} Y_L^m;
+c f1^{Lm} = clm(L,m,1) f^{(L-1)m} - clm(L,m,2) f^{[L+1]m};
+c f2^{Lm} = clm(L,m,3) f^{(L-1)m} - clm(L,m,4) f^{[L+1]m};
+c f3^{Lm} = clm(L,m,5) f^{(L-1)m} - clm(L,m,6) f^{[L+1]m}
+c f4^{Lm} = clm(L,m,7) f^{(L-1)m} + clm(L,m,8) f^{[L+1]m}
+c LL: the angular momentum operator coefficients:
+c LL(k) = k(k+1);
+c
+c DHT: finite difference in co-lattitude;
+c
+c Author: Weijia Kuang
+c Date: Jan., 2001
+c COPY RIGHT: THIS CODE CAN NOT BE DUPLICATED OR USED WITHOUT
+c THE PERMISSION OF THE AUTHOR.
+c
+
+ MODULE mod_sphgeom
+
+ use mod_dimparam
+ use mod_sysparam
+
+ implicit none
+
+ integer, allocatable :: LL(:)
+
+ real (kind=8), allocatable :: aslg(:,:,:),clm(:,:,:)
+ real (kind=8), allocatable, dimension(:) :: th,gauwt,sins,dht
+
+! for IMSL FFT, uncomment the following decalaration
+
+! real (kind=8), allocatable :: table(:)
+
+! for FFTW, uncomment the following declaration
+
+ integer (kind=8) table(2)
+
+ CONTAINS
+
+c The subroutines that define the arrays
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ SUBROUTINE sphgeom_init
+
+ if (.not. allocated(LL)) then
+ allocate(LL(0:Lmax1))
+ endif
+ if (.not. allocated(aslg)) then
+ allocate(aslg(0:Lmaxa,0:mmaxa,ntmax))
+ endif
+ if (.not. allocated(th)) then
+ allocate(th(ntmax))
+ endif
+ if (.not. allocated(gauwt)) then
+ allocate(gauwt(ntmax))
+ endif
+ if (.not. allocated(sins)) then
+ allocate(sins(ntmax))
+ endif
+ if (.not. allocated(clm)) then
+ allocate(clm(0:Lmax1,0:mmax,8))
+ endif
+ if (.not. allocated(dht)) then
+ allocate(dht(ntmax))
+ endif
+
+! for IMSL FFT, uncoomment the following commands
+
+! if (.not. allocated(table))
+! allocate(table(4*npmax+30))
+! endif
+
+ END SUBROUTINE sphgeom_init
+
+********************************************************************
+********************************************************************
+
+ SUBROUTINE sphgeom
+
+ implicit none
+ integer i,np2,L,m
+ real (kind=8) one,c1
+ real (kind=8) x(npmax,ntmax)
+ complex (kind=8) y(0:npmax/2,ntmax)
+ real (kind=8) tempfftin(npmax)
+ complex( kind=8) tempfftout(npmax/2+1)
+
+ np2 = npmax/2
+ one = 1.0
+
+ call sphgeom_init
+
+ aslg = 0.0
+ th = 0.0
+ gauwt = 0.0
+ sins = 0.0
+ table = 0.0
+ dht = 0.0
+ clm = 0.0
+ LL = 0
+
+ call gauleg(-one,one,th,gauwt,ntmax)
+ do i = 1,ntmax
+ call aslegend(aslg(0,0,i),th(i),Lmaxa,mmaxa,2)
+ sins(i) = sqrt(one-th(i)**2)
+ enddo
+
+c call dfftri(npmax,table)
+
+c Initialize fftw plan, the last parameter 0 is for ESTIMATE
+
+ call dfftw_plan_dft_r2c_1d(table(1),npmax,tempfftin,
+ & tempfftout,0)
+ call dfftw_plan_dft_c2r_1d(table(2),npmax,tempfftout,
+ & tempfftin,0)
+
+ do i = 1,ntmax-1
+ dht(i) = abs(th(i+1)-th(i))
+ enddo
+ dht(ntmax) = dht(ntmax-1)
+
+ do m = 0,mmax
+ do L = m+1,Lmax1
+ c1 = sqrt(one*(L-m)*(L+m)/(one*(2*L-1)*(2*L+1)))
+ clm(L,m,1) = (L-1.0)*c1
+ clm(L,m,3) = (L+1.0)*c1
+ clm(L,m,5) = (L-2.0)*c1
+ clm(L,m,7) = c1
+ enddo
+ do L = m,Lmax
+ c1 = sqrt(one*(L-m+1)*(L+m+1)/
+ & (one*(2*L+1)*(2*L+3)))
+ clm(L,m,2) = (L+2.0)*c1
+ clm(L,m,4) = L*c1
+ clm(L,m,6) = (L+3.0)*c1
+ clm(L,m,8) = c1
+ enddo
+ enddo
+
+ do L = 0,Lmax1
+ LL(L) = L*(L+1)
+ enddo
+
+ END SUBROUTINE sphgeom
+
+********************************************************************
+********************************************************************
+
+ END MODULE mod_sphgeom
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+c
+c This module defines the parity conditions for velocity field
+c presentations; the radial dimensional variations of the magnetic
+c field according to electrical conductivity of the boundaries.
+c
+c LSYM: The symmetry parity indicies. For given m >= 1,
+c L = m, m+1,..., m+2*LSYM(m)+1;
+c For m = 0,
+c L = 1, 2,..., 2*LSYM(0);
+c
+c LSYM1: The upper limit of m for given L;
+c
+c KDM: The dimensional indices for velocity field
+c representation in (L,m);
+c
+c NMB: The total dimension of the magnetic field in radius;
+c NMBIC: The radial dimension of the magnetic field in the
+c inner core;
+c NMBDP: The radial dimension of the magnetic field in the
+c D"-layer;
+c
+c Author: Weijia Kuang
+c Date: Mar., 2002
+c COPY RIGHT: THIS CODE CAN NOT BE DUPLICATED OR USED WITHOUT
+c THE PERMISSION OF THE AUTHOR.
+c
+
+ MODULE mod_parity
+
+ use mod_dimparam
+ use mod_optparam
+
+ implicit none
+
+ integer nmb,nmbic,nmbdp
+ integer, allocatable, dimension(:) :: lsym_v,lsym_m,lsym_t,
+ & kdm
+
+ CONTAINS
+
+c The subroutines that define the above parameters
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ SUBROUTINE parity_init
+
+ if (.not. allocated(lsym_v)) then
+ allocate(lsym_v(0:mmax_v))
+ endif
+ if (.not. allocated(lsym_m)) then
+ allocate(lsym_m(0:Lmax_m))
+ endif
+ if (.not. allocated(lsym_t)) then
+ allocate(lsym_t(0:Lmax_t))
+ endif
+ if (.not. allocated(kdm)) then
+ allocate(kdm(0:mmax_v))
+ endif
+
+ END SUBROUTINE parity_init
+
+***********************************************************************
+***********************************************************************
+
+ SUBROUTINE parity
+
+ implicit none
+ integer i,m
+
+ call parity_init
+
+ i = mod(Lmax_v,2)
+ if (i. eq. 0) then
+ lsym_v(0) = Lmax_v/2
+ else
+ lsym_v(0) = (Lmax_v-1)/2
+ endif
+ do m = 1,mmax_v
+ i = mod(Lmax_v-m+1,2)
+ if (i .eq. 0) then
+ lsym_v(m) = (Lmax_v-m-1)/2
+ else
+ lsym_v(m) = (Lmax_v-m-2)/2
+ endif
+ enddo
+
+ do i = 0,Lmax_m
+ lsym_m(i) = min0(i,mmax_m)
+ enddo
+ do i = 0,Lmax_t
+ lsym_t(i) = min0(i,mmax_t)
+ enddo
+
+ kdm(0) = 2*lsym_v(0)*nmxo1
+ do m = 1,mmax_v
+ kdm(m)= kdm(m-1)+2*(lsym_v(m)+1)*nmxo1
+ enddo
+
+ nmbic = nmx4
+ nmbdp = nmx5
+ if (kicbb .le. 1) nmbic = 0
+ if (kcmbb .le. 1) nmbdp = 0
+
+ nmb = nmx3+nmbic+nmbdp
+
+ END SUBROUTINE parity
+
+********************************************************************
+********************************************************************
+
+ END MODULE mod_parity
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
Copied: geodyn/3D/MoSST/trunk/mosst-lib (from rev 5193, geodyn/3D/MoSST/trunk/lib)
More information about the cig-commits
mailing list