[cig-commits] r6702 - geodyn/3D/MAG/trunk/src
wei at geodynamics.org
wei at geodynamics.org
Thu Apr 26 17:01:21 PDT 2007
Author: wei
Date: 2007-04-26 17:01:20 -0700 (Thu, 26 Apr 2007)
New Revision: 6702
Modified:
geodyn/3D/MAG/trunk/src/cmbcoeff.f
geodyn/3D/MAG/trunk/src/nl.f
geodyn/3D/MAG/trunk/src/par.bnch1mv
geodyn/3D/MAG/trunk/src/unscramble.f
Log:
More changes in cmbcoeff.f, it work.
Modified: geodyn/3D/MAG/trunk/src/cmbcoeff.f
===================================================================
--- geodyn/3D/MAG/trunk/src/cmbcoeff.f 2007-04-26 23:47:56 UTC (rev 6701)
+++ geodyn/3D/MAG/trunk/src/cmbcoeff.f 2007-04-27 00:01:20 UTC (rev 6702)
@@ -20,6 +20,7 @@
dimension la(nlma),ma(nlma)
dimension alm(lmax,lmax),blm(lmax,lmax)
dimension glm(lmax,lmax),hlm(lmax,lmax)
+ dimension cc(2.*nlma)
c
c constants
c
@@ -40,28 +41,10 @@
cdepth=rc-ri ! Depth of outer core
rcre=rc/re ! core/surface radius ratio
rcsqinv=(cdepth/rc)**2 ! dim-less inverse squared core radius
-c
-c harmonic factors
-c
- fl=float(l)
- flp2=fl + 2.
- fl2p1=(2.*fl) + 1.
- fm=float(m)
- fact=sqrt(fl2p1*pi4inv)
- fact1=((-1.)**fm)*fl*fact
- if(m.gt.0) fact1=fact1*sqrt(2.)
- fact2=rcsqinv*(rcre**flp2)
c
-c define MAG harmonic-Gauss harmonic conversion factors
-c conalm,conblm
+c write header to cc file
c
- conalm=fact1
- conblm=-fact1
-
-c
-c write header
-c
write(21,2100) nlma,lmax,minc,r(1),r(kc),time/tscale
2100 format(/, 2x,"nlma=",i3,2x,"lmax=",i3,2x,"minc=",i3,2x,
$ "r(1)=",f7.4,2x,"r(kc)=",f7.4,2x,"time/tscale=",f9.6)
@@ -73,30 +56,22 @@
c write(21,2102) (dw(lm,kc),lm=1,nlma)
c write(21,2102) (z(lm,kc),lm=1,nlma)
2101 format(256(1X,f9.5))
- 2102 format(256(1X,f9.3))
+c 2102 format(256(1X,f9.3))
c processing unscramble start here
+c write header to cg file
write(22,2200) nlma,lmax,minc,r(1),r(kc),time/tscale
2200 format(/, 2x,"nlma=",i3,2x,"lmax=",i3,2x,"minc=",i3,2x,
$ "r(1)=",f7.4,2x,"r(kc)=",f7.4,2x,"time/tscale=",f9.6)
-c define intermediate indices
-c defined in param.f
-c mmax=(lmax/minc)*minc
-c nmaf=mmax+1
-c nlaf=lmax+1
-c print out mmax nmaf and nlaf
-c write(22,2203) mmax,nmaf,nlaf
-c 2203 format(/,2x,3i3)
-
c define the unscramble array mclm(lm)
lm=0
do 35 mc=1,nmaf,minc
do 31 lc=mc,nlaf
lm=lm+1
mclm(lm)=mc
- write(22,2201) lm, mclm(lm)
- 2201 format(/, 2x,i3,2x,i4)
+c write(22,2201) lm, mclm(lm)
+c 2201 format(/, 2x,i3,2x,i4)
31 continue
35 continue
@@ -109,8 +84,8 @@
tl3=-ma(lm)*(lmax+1)/minc
la(lm)=tl1+tl2+tl3
c PRINT lm, la(lm), ma(lm) HERE
- write(22,2202) lm, la(lm), ma(lm)
- 2202 format(/, 2x,i3,2x,i4,2x,i4)
+c write(22,2202) lm, la(lm), ma(lm)
+c 2202 format(/,2x,i3,2x,i4,2x,i4)
36 continue
c
@@ -118,12 +93,13 @@
c and assign the new indices
do 37 lm=1,nlma
-cc READ (cc-file) c1,c2
+ c1=real(b(lm,1))
+ c2=aimag(b(lm,1))
c assign new indices
- l=la(lm)
- m=ma(lm)
- alm(l,m)=c1
- blm(l,m)=c2
+ l=la(lm)
+ m=ma(lm)
+ alm(l,m)=c1
+ blm(l,m)=c2
c Convertion starts here
c the following code converts between Gauss coefficients (glm, hlm)
@@ -155,21 +131,43 @@
c conductivity, and uses the depth of the outer core as the
c length scale.
+c
+c harmonic factors
+c
+ fl=float(l)
+ flp2=fl + 2.
+ fl2p1=(2.*fl) + 1.
+ fm=float(m)
+ fact=sqrt(fl2p1*pi4inv)
+ fact1=((-1.)**fm)*fl*fact
+ if(m.gt.0) fact1=fact1*sqrt(2.)
+ fact2=rcsqinv*(rcre**flp2)
+
+c
+c define MAG harmonic-Gauss harmonic conversion factors
+c conalm,conblm
+c
+ conalm=fact1
+ conblm=-fact1
+
if (l .le. 0 .or. m .lt. 0 .or. m .gt. l) then
write(6,'(''bad l or m in getgauss'')')
return
endif
- if (id .gt. 0) then ! form Gauss coeffs in nT
+c if (id .gt. 0) then ! form Gauss coeffs in nT
glm(l,m)=anano*escale*fact2*conalm*alm(l,m)
hlm(l,m)=anano*escale*fact2*conblm*blm(l,m)
- return
- else ! form dimensionless fully normalized potential coeffs
- alm(l,m)=glm(l,m)/(anano*escale*fact2*conalm)
- blm(l,m)=hlm(l,m)/(anano*escale*fact2*conblm)
- return
- endif
+ write(22,2203) l,m,alm(l,m),blm(l,m),glm(l,m),hlm(l,m)
+ 2203 format(/,2x,2i3,2x,2(f9.5),2x,2(f15.5))
+c return
+c else ! form dimensionless fully normalized potential coeffs
+c alm(l,m)=glm(l,m)/(anano*escale*fact2*conalm)
+c blm(l,m)=hlm(l,m)/(anano*escale*fact2*conblm)
+c return
+c endif
+
37 continue
return
Modified: geodyn/3D/MAG/trunk/src/nl.f
===================================================================
--- geodyn/3D/MAG/trunk/src/nl.f 2007-04-26 23:47:56 UTC (rev 6701)
+++ geodyn/3D/MAG/trunk/src/nl.f 2007-04-27 00:01:20 UTC (rev 6702)
@@ -115,7 +115,7 @@
c
character chd*3,chg*3
character*72 movefile,movmfile,movafile,logfile,grafile,
- $ lsfile,lpfile,ccfile
+ $ lsfile,lpfile,ccfile, cgfile
c---------------------------------------------------------------
c
ifirst=1
@@ -138,6 +138,7 @@
movafile='ma.'//outfile
movmfile='mm.'//outfile
ccfile='cc.'//outfile
+ cgfile='cg.'//outfile
lpfile='lp.'//outfile
open(15,file=logfile,status='new',form='formatted')
open(16,file=lsfile,status='new',form='formatted')
@@ -151,6 +152,7 @@
$ open(20,file=movmfile,status='new',form='formatted')
if(imovopt.ge.1000)
$ open(21,file=ccfile,status='new',form='formatted')
+ open(22,file=cgfile,status='new',form='formatted')
c
c write header for movie file
c
Modified: geodyn/3D/MAG/trunk/src/par.bnch1mv
===================================================================
--- geodyn/3D/MAG/trunk/src/par.bnch1mv 2007-04-26 23:47:56 UTC (rev 6701)
+++ geodyn/3D/MAG/trunk/src/par.bnch1mv 2007-04-27 00:01:20 UTC (rev 6702)
@@ -30,8 +30,8 @@
ktopv=2,
kbotv=2,
difamp=0,
- imovopt=1111,
- tmovstart=2.0,
+ imovopt=1101,
+ tmovstart=0.02,
tmovstep=2.E-3,
iframes=100,
icour=4 /
Modified: geodyn/3D/MAG/trunk/src/unscramble.f
===================================================================
--- geodyn/3D/MAG/trunk/src/unscramble.f 2007-04-26 23:47:56 UTC (rev 6701)
+++ geodyn/3D/MAG/trunk/src/unscramble.f 2007-04-27 00:01:20 UTC (rev 6702)
@@ -1,3 +1,4 @@
+ program unscramble
c formulas in unscramble.f were taken from the prep.f and other mag files
c it is used to unscramble the cc-files, converting the cc-file
c output to arrays in harmonic degree l and order m.
@@ -4,11 +5,18 @@
c It defines arrays al(lm) and am(lm) for lm=1,nlma,
c which contain the l and m values we need to convert the cc-file contents
+ integer lm, mc, lc
+ real mclm(lm),la(lm),ma(lm)
+ paramter (lmax, minc, nlma)
+
+c open cc file for reading
+
+c open (21,
c input lmax, minc, and nlma from cc file header HERE:
- lmax=??
- minc=??
- nlma=??
+ lmax=153
+ minc=4
+ nlma=32
c define intermediate indices
mmax=(lmax/minc)*minc
@@ -38,14 +46,15 @@
c if this works, then read the contents of the cc-file block in pairs,
c and assign the new indices
- do 37 lm=1,nlma
- READ (cc-file) c1,c2
+cc do 37 lm=1,nlma
+cc READ (cc-file) c1,c2
c assign new indices
- l=la(lm)
- m=ma(lm)
- anew(l,m)=c1
- bnew(l,m)=c2
- 37 continue
+cc l=la(lm)
+cc m=ma(lm)
+cc anew(l,m)=c1
+cc bnew(l,m)=c2
+cc 37 continue
-
\ No newline at end of file
+
+ end
More information about the cig-commits
mailing list