[cig-commits] [commit] devel, master: cosmetic changes in src/meshfem3D/model_s362ani.f90 (9daa3eb)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Nov 6 08:10:37 PST 2014
Repository : https://github.com/geodynamics/specfem3d_globe
On branches: devel,master
Link : https://github.com/geodynamics/specfem3d_globe/compare/bc58e579b3b0838a0968725a076f5904845437ca...be63f20cbb6f462104e949894dbe205d2398cd7f
>---------------------------------------------------------------
commit 9daa3eb6da30840366978605801996068421af99
Author: Dimitri Komatitsch <komatitsch at lma.cnrs-mrs.fr>
Date: Sun Apr 13 15:32:08 2014 +0200
cosmetic changes in src/meshfem3D/model_s362ani.f90
>---------------------------------------------------------------
9daa3eb6da30840366978605801996068421af99
src/meshfem3D/model_s362ani.f90 | 151 +++++++++++++++-------------------------
1 file changed, 55 insertions(+), 96 deletions(-)
diff --git a/src/meshfem3D/model_s362ani.f90 b/src/meshfem3D/model_s362ani.f90
index 8ce76a2..c65d050 100644
--- a/src/meshfem3D/model_s362ani.f90
+++ b/src/meshfem3D/model_s362ani.f90
@@ -83,8 +83,8 @@
! standard routine to setup model
- use constants
use model_s362ani_par
+ use constants
implicit none
@@ -185,7 +185,8 @@
numvar,ivarkern,varstr, &
refmdl,kerstr,hsplfl,dskker,ierror)
else
- write(6,"('the model ',a,' does not exits')") modeldef(1:len_trim(modeldef))
+ write(6,"('model ',a,' does not exist')") modeldef(1:len_trim(modeldef))
+ stop 'model does not exist in s362_ani'
endif
! check arrays
@@ -284,8 +285,7 @@
enddo
call chebyfun(u2,13,chebyshev2)
do i=1+nskip+nupper,nskip+nupper+nlower
- dvercof(i)=(chebyshev2(i-nskip-nupper)- &
- chebyshev(i-nskip-nupper))/ddep
+ dvercof(i)=(chebyshev2(i-nskip-nupper) - chebyshev(i-nskip-nupper))/ddep
enddo
endif
else if(string(1:13) == 'WDC+SHSVWM20A') then
@@ -372,11 +372,8 @@
! vercof(23)=1.
! vercof(24)=1.
! vercof(25)=1.
- else if( &
- (string(1:lstr) == 'WDC+ANI_362_U6L8'.and.lstr == 16) &
- .or. &
- (string(1:lstr) == 'WDC+ANI_362_U6L8_TOPO'.and.lstr == 21) &
- ) then
+ else if( (string(1:lstr) == 'WDC+ANI_362_U6L8'.and.lstr == 16) &
+ .or. (string(1:lstr) == 'WDC+ANI_362_U6L8_TOPO'.and.lstr == 21) ) then
if(upper) then
nspl=6
splpts(1)=24.4
@@ -439,11 +436,8 @@
vercof(30)=1.
vercof(31)=1.
vercof(32)=1.
- else if( &
- (string(1:lstr) == 'WDC+ANI_362_U6L8_650'.and.lstr == 20) &
- .or. &
- (string(1:lstr) == 'WDC+ANI_362_U6L8_TOPO_650'.and.lstr == 25) &
- ) then
+ else if( (string(1:lstr) == 'WDC+ANI_362_U6L8_650'.and.lstr == 20) &
+ .or. (string(1:lstr) == 'WDC+ANI_362_U6L8_TOPO_650'.and.lstr == 25) ) then
if(upper_650) then
nspl=6
splpts(1)=24.4
@@ -472,8 +466,7 @@
vercof(1)=1.
vercof(22)=1.
vercof(23)=1.
- else if(string(1:lstr) == 'WDC+WM_362_U6L8_650' &
- .and.lstr == 19) then
+ else if(string(1:lstr) == 'WDC+WM_362_U6L8_650' .and.lstr == 19) then
if(upper_650) then
nspl=6
splpts(1)=24.4
@@ -579,11 +572,8 @@
vercof(34)=1.
vercof(35)=1.
vercof(36)=1.
- else if( &
- (string(1:lstr) == 'WDC+U8L8_I1D_650'.and.lstr == 16) &
- .or. &
- (string(1:lstr) == 'WDC+U8L8_I3D_650'.and.lstr == 16) &
- ) then
+ else if( (string(1:lstr) == 'WDC+U8L8_I1D_650'.and.lstr == 16) &
+ .or. (string(1:lstr) == 'WDC+U8L8_I3D_650'.and.lstr == 16) ) then
if(upper_650) then
nspl=8
splpts(1)=24.4
@@ -772,10 +762,7 @@
1.00196657023780,1.0015515913133,1.0012554932754,1.0010368069141, &
1.00087070107920,1.0007415648034 /
- if(kmax > 13)then
- write(*,"(' kmax exceeds the limit in chebyfun')")
- stop
- endif
+ if(kmax > 13) stop ' kmax exceeds the limit in chebyfun'
f(0)=1.0
f(1)=u
@@ -875,9 +862,7 @@
dskker(i)=desckern(i)
do j=1,ncoefhor(ihpakern(i))
coe(j,i)=coef(j,i)
-! if(j == 1) then
-! write(6,"(e12.4)") coe(j,i)
-! endif
+! if(j == 1) write(6,"(e12.4)") coe(j,i)
enddo
enddo
@@ -970,8 +955,7 @@
if ( ios /= 0 ) then
write(IMAIN,*) 'error opening "', trim(filename), '": ', ios
call flush_IMAIN()
- ! stop
- call exit_MPI(0, 'error model s362ani')
+ call exit_MPI(0, 'error in model s362ani')
endif
do while (ios == 0)
@@ -982,7 +966,7 @@
substr=string(17:lstr)
ifst=1
ilst=len_trim(substr)
- do while (substr(ifst:ifst) == ' '.and.ifst < ilst)
+ do while (substr(ifst:ifst) == ' ' .and. ifst < ilst)
ifst=ifst+1
enddo
if(ilst-ifst <= 0) then
@@ -994,7 +978,7 @@
substr=string(12:len_trim(string))
ifst=1
ilst=len_trim(substr)
- do while (substr(ifst:ifst) == ' '.and.ifst < ilst)
+ do while (substr(ifst:ifst) == ' ' .and. ifst < ilst)
ifst=ifst+1
enddo
if(ilst-ifst <= 0) then
@@ -1059,8 +1043,7 @@
lmaxhor(idummy)=0
ncoefhor(idummy)=ncoef
do i=1,ncoef
- read(lu,*) ixlspl(i,idummy),xlaspl(i,idummy), &
- xlospl(i,idummy),xraspl(i,idummy)
+ read(lu,*) ixlspl(i,idummy),xlaspl(i,idummy),xlospl(i,idummy),xraspl(i,idummy)
enddo
endif
else if(string(1:4) == 'STRU'.and.string(9:9) == ':') then
@@ -1154,8 +1137,7 @@
!-------------------------------------------------------------------------------------------------
!
-
-! --- evaluate perturbations in per cent
+! --- evaluate perturbations in percent
subroutine model_s362ani_subshsv(xcolat,xlon,xrad,dvsh,dvsv,dvph,dvpv)
@@ -1201,6 +1183,7 @@
y=90.0-xcolat
x=xlon
+
do ihpa=1,numhpa
if(itypehpa(ihpa) == 1) then
lmax=lmxhpa(ihpa)
@@ -1228,7 +1211,7 @@
endif
enddo
-! --- evaluate 3-D perturbations in velocity and anisotropy
+! evaluate 3-D perturbations in velocity and anisotropy
valu(1)=0. ! --- velocity
valu(2)=0. ! --- anisotropy
@@ -1240,14 +1223,12 @@
lstr=len_trim(varstr(ivarkern(iker)))
vstr=(varstr(ivarkern(iker)))
if(ieval == 1) then
- if(vstr(1:lstr) == 'UM (SH+SV)*0.5,'.or. &
- vstr(1:lstr) == 'LM (SH+SV)*0.5,'.or. &
+ if(vstr(1:lstr) == 'UM (SH+SV)*0.5,'.or. vstr(1:lstr) == 'LM (SH+SV)*0.5,'.or. &
vstr(1:lstr) == 'EA (SH+SV)*0.5,') then
isel=1
endif
else if(ieval == 2) then
- if(vstr(1:lstr) == 'UM SH-SV,'.or. &
- vstr(1:lstr) == 'LM SH-SV,'.or. &
+ if(vstr(1:lstr) == 'UM SH-SV,'.or. vstr(1:lstr) == 'LM SH-SV,'.or. &
vstr(1:lstr) == 'EA SH-SV,') then
isel=1
endif
@@ -1259,19 +1240,16 @@
ihpa=ihpakern(iker)
nylm=(lmxhpa(ihpakern(iker))+1)**2
do i=1,nylm
- value=value+vercof(iker)*ylmcof(i,ihpa) &
- *coe(i,iker)
+ value=value+vercof(iker)*ylmcof(i,ihpa)*coe(i,iker)
enddo
else if(itypehpa(ihpakern(iker)) == 2) then
ihpa=ihpakern(iker)
do i=1,nconpt(ihpa)
iver=iconpt(i,ihpa)
- value=value+vercof(iker)*conpt(i,ihpa) &
- *coe(iver,iker)
+ value=value+vercof(iker)*conpt(i,ihpa)*coe(iver,iker)
enddo
else
- write(6,"('problem 2')")
- stop
+ stop 'problem 2'
endif ! --- itypehpa
endif ! --- vercof(iker) /= 0.
endif ! --- isel == 1
@@ -1280,24 +1258,24 @@
valu(ieval)=value
enddo ! --- ieval
-! --- evaluate perturbations in vsh and vsv
+! evaluate perturbations in vsh and vsv
if(ish == 1) then
vsh3drel=valu(1)+0.5*valu(2)
else if(ish == 0) then
vsv3drel=valu(1)-0.5*valu(2)
else
- stop 'something wrong'
+ stop 'something is wrong in model_s362ani_subshsv'
endif
enddo ! --- by ish
-! --- evaluate perturbations in percent
+! evaluate perturbations in percent
dvsh=vsh3drel
dvsv=vsv3drel
- dvph=0.55*dvsh ! --- scaling used in the inversion
- dvpv=0.55*dvsv ! --- scaling used in the inversion
+ dvph=0.55*dvsh ! scaling used in the inversion
+ dvpv=0.55*dvsv ! scaling used in the inversion
end subroutine model_s362ani_subshsv
@@ -1331,7 +1309,7 @@
! -------------------------------------
-! --- contributing horizontal basis functions at xlat,xlon
+! contributing horizontal basis functions at xlat,xlon
y=90.0-xcolat
x=xlon
@@ -1357,13 +1335,12 @@
! xlospl(1:numcof,ihpa),radspl(1:numcof,ihpa), &
! nconpt(ihpa),iconpt(1:maxver,ihpa),conpt(1:maxver,ihpa))
-
else
write(6,"('problem 1')")
endif
enddo
-! --- evaluate topography (depression) in km
+! evaluate topography (depression) in km
valu(1)=0. ! --- 410
valu(2)=0. ! --- 650
@@ -1398,8 +1375,7 @@
value=value+conpt(i,ihpa)*coe(iver,iker)
enddo
else
- write(6,"('problem 2')")
- stop
+ stop 'problem 2'
endif ! --- itypehpa
endif ! --- isel == 1
enddo ! --- end of do iker=1,numker
@@ -1454,17 +1430,11 @@
interval=np
endif
- if(interval == 0) then
-! write(6,"('low value:',2f10.3)") x,xarr(1)
- else if(interval > 0.and.interval < np) then
-! write(6,"('bracket:',i5,3f10.3)") interval,xarr(interval),x,xarr(interval+1)
- else
-! write(6,"('high value:',2f10.3)") xarr(np),x
- endif
-
do ib=1,np
+
val=0.
vald=0.
+
if(ib == 1) then
r1=(x-xarr(1))/(xarr(2)-xarr(1))
@@ -1488,6 +1458,7 @@
r13d=-1./(xarr(2)-xarr(1))
if(interval == ib.or.interval == 0) then
+
if(iflag == 0) then
val=r1*r4*r10 + r2*r5*r10 + r2*r6*r11 +r13**3
vald=r1d*r4*r10+r1*r4d*r10+r1*r4*r10d
@@ -1495,8 +1466,7 @@
vald=vald+r2d*r6*r11+r2*r6d*r11+r2*r6*r11d
vald=vald+3.*r13d*r13**2
else if(iflag == 1) then
- val=0.6667*(r1*r4*r10 + r2*r5*r10 + r2*r6*r11 &
- + 1.5*r13**3)
+ val=0.6667*(r1*r4*r10 + r2*r5*r10 + r2*r6*r11 + 1.5*r13**3)
vald=r1d*r4*r10+r1*r4d*r10+r1*r4*r10d
vald=vald+r2d*r5*r10+r2*r5d*r10+r2*r5*r10d
vald=vald+r2d*r6*r11+r2*r6d*r11+r2*r6*r11d
@@ -1565,14 +1535,10 @@
vald=vald+r1d*r4*r9+r1*r4d*r9+r1*r4*r9d
vald=vald+r2d*r5*r9+r2*r5d*r9+r2*r5*r9d
if(iflag == 1) then
- val=val+0.3333*(rr1*rr4*rr10 + rr2*rr5*rr10 + &
- rr2*rr6*rr11)
- vald=vald+0.3333*(rr1d*rr4*rr10+rr1*rr4d*rr10+ &
- rr1*rr4*rr10d)
- vald=vald+0.3333*(rr2d*rr5*rr10+rr2*rr5d*rr10+ &
- rr2*rr5*rr10d)
- vald=vald+0.3333*(rr2d*rr6*rr11+rr2*rr6d*rr11+ &
- rr2*rr6*rr11d)
+ val=val+0.3333*(rr1*rr4*rr10 + rr2*rr5*rr10 + rr2*rr6*rr11)
+ vald=vald+0.3333*(rr1d*rr4*rr10+rr1*rr4d*rr10+ rr1*rr4*rr10d)
+ vald=vald+0.3333*(rr2d*rr5*rr10+rr2*rr5d*rr10+ rr2*rr5*rr10d)
+ vald=vald+0.3333*(rr2d*rr6*rr11+rr2*rr6d*rr11+ rr2*rr6*rr11d)
endif
else if(interval == ib) then
val=r1*r4*r10 + r2*r5*r10 + r2*r6*r11
@@ -1581,8 +1547,7 @@
vald=vald+r2d*r6*r11+r2*r6d*r11+r2*r6*r11d
if(iflag == 1) then
val=val+0.3333*rr2*rr6*rr12
- vald=vald+0.3333*(rr2d*rr6*rr12+rr2*rr6d*rr12+ &
- rr2*rr6*rr12d)
+ vald=vald+0.3333*(rr2d*rr6*rr12+rr2*rr6d*rr12+ rr2*rr6*rr12d)
endif
else if(interval == ib+1) then
val=r2*r6*r12
@@ -1590,6 +1555,7 @@
else
val=0.
endif
+
else if(ib == np-1) then
rr1=(x-xarr(np-2))/(xarr(np)-xarr(np-2))
@@ -1644,8 +1610,7 @@
vald=vald+r2d*r5*r9+r2*r5d*r9+r2*r5*r9d
if(iflag == 1) then
val=val+0.3333*rr1*rr3*rr7
- vald=vald+0.3333*(rr1d*rr3*rr7+rr1*rr3d*rr7+ &
- rr1*rr3*rr7d)
+ vald=vald+0.3333*(rr1d*rr3*rr7+rr1*rr3d*rr7+ rr1*rr3*rr7d)
endif
else if(interval == ib.or.interval == np) then
val=r1*r4*r10 + r2*r5*r10 + r2*r6*r11
@@ -1653,18 +1618,15 @@
vald=vald+r2d*r5*r10+r2*r5d*r10+r2*r5*r10d
vald=vald+r2d*r6*r11+r2*r6d*r11+r2*r6*r11d
if(iflag == 1) then
- val=val+0.3333*(rr1*rr3*rr8 + rr1*rr4*rr9 + &
- rr2*rr5*rr9)
- vald=vald+0.3333*(rr1d*rr3*rr8+rr1*rr3d*rr8+ &
- rr1*rr3*rr8d)
- vald=vald+0.3333*(rr1d*rr4*rr9+rr1*rr4d*rr9+ &
- rr1*rr4*rr9d)
- vald=vald+0.3333*(rr2d*rr5*rr9+rr2*rr5d*rr9+ &
- rr2*rr5*rr9d)
+ val=val+0.3333*(rr1*rr3*rr8 + rr1*rr4*rr9 + rr2*rr5*rr9)
+ vald=vald+0.3333*(rr1d*rr3*rr8+rr1*rr3d*rr8+ rr1*rr3*rr8d)
+ vald=vald+0.3333*(rr1d*rr4*rr9+rr1*rr4d*rr9+ rr1*rr4*rr9d)
+ vald=vald+0.3333*(rr2d*rr5*rr9+rr2*rr5d*rr9+ rr2*rr5*rr9d)
endif
else
val=0.
endif
+
else if(ib == np) then
r1=(x-xarr(np-2))/(xarr(np)-xarr(np-2))
@@ -1703,8 +1665,7 @@
vald=vald+r2d*r5*r9+r2*r5d*r9+r2*r5*r9d
vald=vald+3.*r13d*r13**2
else if(iflag == 1) then
- val=0.6667*(r1*r3*r8 + r1*r4*r9 + r2*r5*r9 + &
- 1.5*r13**3)
+ val=0.6667*(r1*r3*r8 + r1*r4*r9 + r2*r5*r9 + 1.5*r13**3)
vald=r1d*r3*r8+r1*r3d*r8+r1*r3*r8d
vald=vald+r1d*r4*r9+r1*r4d*r9+r1*r4*r9d
vald=vald+r2d*r5*r9+r2*r5d*r9+r2*r5*r9d
@@ -1778,13 +1739,12 @@
complex TEMP,FAC,DFAC
- !real(kind=4) WK1(1),WK2(1),WK3(1),Y(1),XLAT,XLON
-
integer :: LMAX
!
! WK1,WK2,WK3 SHOULD BE DIMENSIONED AT LEAST (LMAX+1)*4
!
+! real(kind=4) WK1(1),WK2(1),WK3(1),Y(1),XLAT,XLON
real(kind=4),dimension(LMAX+1) :: WK1,WK2,WK3
real(kind=4) :: XLAT,XLON
real(kind=4),dimension((LMAX+1)**2) :: Y !! Y should go at least from 1 to fac(LMAX)
@@ -1806,8 +1766,8 @@
! index L goes from 0 to LMAX
L=IL1-1
- !! see legndr(): WK1,WK2,WK3 should go from 1 to L+1
- !CALL legndr(THETA,L,L,WK1,WK2,WK3)
+ ! see legndr(): WK1,WK2,WK3 should go from 1 to L+1
+ ! CALL legndr(THETA,L,L,WK1,WK2,WK3)
CALL legndr(THETA,L,L,WK1(1:L+1),WK2(1:L+1),WK3(1:L+1))
FAC=(1.,0.)
@@ -1838,6 +1798,7 @@
implicit none
!real(kind=4) :: X(2),XP(2),XCOSEC(2) !! X, XP, XCOSEC should go from 1 to M+1
+ real(kind=4) :: X(M+1),XP(M+1),XCOSEC(M+1) !! X, XP, XCOSEC should go from 1 to M+1
double precision :: SMALL,sumval,COMPAR,CT,ST,FCT,COT,X1,X2,X3,F1,F2,XM,TH
@@ -1847,8 +1808,6 @@
real(kind=4) :: THETA,DSFL3,COSEC,SFL3
- real(kind=4) :: X(M+1),XP(M+1),XCOSEC(M+1) !! X, XP, XCOSEC should go from 1 to M+1
-
!!!!!! illegal statement, removed by Dimitri Komatitsch DFLOAT(I)=FLOAT(I)
sumval=0.D0
@@ -1869,7 +1828,7 @@
XP(I)=0.
enddo
- IF(L > 1.AND.ABS(THETA) > 1.E-5) GO TO 3
+ IF(L > 1.AND.ABS(THETA) > 1.E-5) goto 3
X(1)=FCT
IF(L == 0) RETURN
X(1)=CT*FCT
More information about the CIG-COMMITS
mailing list