[cig-commits] [commit] devel: improved the source code for attenuation models: fixed a double to single precision conversion statement, replaced Qmu = 0 with Qmu = 4999 in model_1dref.f90 because Qmu = 0 makes no sense (41f2c30)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Mon Apr 21 14:00:48 PDT 2014
Repository : ssh://geoshell/specfem3d_globe
On branch : devel
Link : https://github.com/geodynamics/specfem3d_globe/compare/46b00e6b38b82d267194090fb713f38485f162cc...9419f59565e2229712f8b5429d09a2dc53fcec62
>---------------------------------------------------------------
commit 41f2c3022e0dd2fbb42d39b9be44879edbb8ba23
Author: Dimitri Komatitsch <komatitsch at lma.cnrs-mrs.fr>
Date: Mon Apr 21 22:57:16 2014 +0200
improved the source code for attenuation models: fixed a double to single precision conversion statement, replaced Qmu = 0 with Qmu = 4999 in model_1dref.f90 because Qmu = 0 makes no sense
>---------------------------------------------------------------
41f2c3022e0dd2fbb42d39b9be44879edbb8ba23
flags.guess | 17 --
src/meshfem3D/get_model.f90 | 41 +++-
src/meshfem3D/meshfem3D_models.f90 | 2 +-
src/meshfem3D/model_1dref.f90 | 377 ++++++++++++++++++-------------------
src/meshfem3D/model_s362ani.f90 | 68 ++-----
5 files changed, 225 insertions(+), 280 deletions(-)
diff --git a/flags.guess b/flags.guess
index 0a61a76..3c0c454 100644
--- a/flags.guess
+++ b/flags.guess
@@ -17,16 +17,6 @@
# or -mcmodel=medium -shared-intel (if you use the Intel ifort / icc compiler)
# to the configure options of CFLAGS, FCFLAGS and LDFLAGS otherwise the compiler will display an error
# message (for instance 'relocation truncated to fit: R\_X86\_64\_PC32 against .bss' or something similar);
-# BEWARE that using -mcmodel=medium -shared-intel is known for currently leading to incorrect seismograms,
-# at least when flag ATTENUATION is on (and maybe even without), at least in the case of the Intel ifort / icc compiler.
-# This likely comes from intrinsic functions such as size() that return and integer8 instead of an integer4, thus leading
-# to incorrect results when used in function calls is the -i8 flag is not added to the compiler options;
-# however, when adding -i8 the code does not compile because the MPI calls then refuse to compile (they need integer4 as arguments).
-# It seems that this problem is triggered when 3D models are used and mostly (only?) when attenuation is on, but not triggered for 1D Earth model benchmarks
-# even with attenuation (at least not triggered for EXAMPLES/benchmarks/attenuation_benchmark_GJI_2002_versus_normal_modes
-# nor for Anne Sieminski's three 1D benchmarks), which makes it uneasy to locate.
-# Since most current users run the code with less than 2 GB of static memory per core, we have not investigated that problem carefully for now,
-# and thus recommend that you do NOT use these compiler flags until this problem is located and fixed.
#
###########################################################################################################################
###########################################################################################################################
@@ -90,13 +80,6 @@ case $my_FC in
# I/O throughput lingers at 2.5 MB/s, with it it can increase to ~44 MB/s
# However it does not make much of a difference on NFS mounted volumes or with SFS 3.1.1 / Lustre 1.6.7.1
DEF_FFLAGS="-O3 -check nobounds -xHost -fpe0 -ftz -assume buffered_io -assume byterecl -align sequence -vec-report0 -std03 -diag-disable 6477 -implicitnone -gen-interfaces -warn all" # -mcmodel=medium -shared-intel
-##############################################################################################################################
-##############################################################################################################################
-##############################################################################################################################
-# adding -mcmodel=medium -shared-intel can/will lead to problems, see the detailed comment at the beginning of this file
-##############################################################################################################################
-##############################################################################################################################
-##############################################################################################################################
# useful for debugging...
# for debugging: change -O3 -check nobounds to -check all -debug -g -O0 -fp-stack-check -traceback -ftrapuv
#
diff --git a/src/meshfem3D/get_model.f90 b/src/meshfem3D/get_model.f90
index 57ce6d4..14f8762 100644
--- a/src/meshfem3D/get_model.f90
+++ b/src/meshfem3D/get_model.f90
@@ -305,18 +305,41 @@
if(ATTENUATION) then
if(ATTENUATION_3D .or. ATTENUATION_1D_WITH_3D_STORAGE) then
- do i_sls = 1,N_SLS
- tau_e_store(i,j,k,i_sls,ispec) = tau_e(i_sls)
- enddo
- Qmu_store(i,j,k,ispec) = Qmu
- else
- ! store values from mid-point for whole element
- if( i == NGLLX/2 .and. j == NGLLY/2 .and. k == NGLLZ/2 ) then
+
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ do i_sls = 1,N_SLS
+ tau_e_store(i,j,k,i_sls,ispec) = sngl(tau_e(i_sls))
+ enddo
+ Qmu_store(i,j,k,ispec) = sngl(Qmu)
+ else
do i_sls = 1,N_SLS
- tau_e_store(1,1,1,i_sls,ispec) = tau_e(i_sls)
+ tau_e_store(i,j,k,i_sls,ispec) = tau_e(i_sls)
enddo
- Qmu_store(1,1,1,ispec) = Qmu
+ Qmu_store(i,j,k,ispec) = Qmu
endif
+
+ else
+
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ ! store values from mid-point for whole element
+ if( i == NGLLX/2 .and. j == NGLLY/2 .and. k == NGLLZ/2 ) then
+ do i_sls = 1,N_SLS
+ tau_e_store(1,1,1,i_sls,ispec) = sngl(tau_e(i_sls))
+ enddo
+ Qmu_store(1,1,1,ispec) = sngl(Qmu)
+ endif
+ else
+ ! store values from mid-point for whole element
+ if( i == NGLLX/2 .and. j == NGLLY/2 .and. k == NGLLZ/2 ) then
+ do i_sls = 1,N_SLS
+ tau_e_store(1,1,1,i_sls,ispec) = tau_e(i_sls)
+ enddo
+ Qmu_store(1,1,1,ispec) = Qmu
+ endif
+ endif
+
endif
endif
diff --git a/src/meshfem3D/meshfem3D_models.f90 b/src/meshfem3D/meshfem3D_models.f90
index 818a44f..d4e7444 100644
--- a/src/meshfem3D/meshfem3D_models.f90
+++ b/src/meshfem3D/meshfem3D_models.f90
@@ -81,7 +81,7 @@
call model_ak135_broadcast(myrank,CRUSTAL)
case(REFERENCE_MODEL_1DREF)
- call model_1dref_broadcast(myrank,CRUSTAL)
+ call model_1dref_broadcast(CRUSTAL)
case(REFERENCE_MODEL_SEA1D)
call model_sea1d_broadcast(myrank,CRUSTAL)
diff --git a/src/meshfem3D/model_1dref.f90 b/src/meshfem3D/model_1dref.f90
index 3fa4689..775a561 100644
--- a/src/meshfem3D/model_1dref.f90
+++ b/src/meshfem3D/model_1dref.f90
@@ -54,7 +54,7 @@
integer, parameter :: NR_REF = 750
! model_1dref_variables
- double precision, dimension(:), allocatable :: &
+ double precision, dimension(NR_REF) :: &
Mref_V_radius_ref,Mref_V_density_ref, &
Mref_V_vpv_ref,Mref_V_vph_ref, &
Mref_V_vsv_ref,Mref_V_vsh_ref, &
@@ -66,7 +66,7 @@
!-------------------------------------------------------------------------------------------------
!
- subroutine model_1dref_broadcast(myrank,CRUSTAL)
+ subroutine model_1dref_broadcast(CRUSTAL)
! standard routine to setup model
@@ -74,25 +74,8 @@
implicit none
- integer :: myrank
logical :: CRUSTAL
- ! local parameters
- integer :: ier
-
- ! allocates model arrays
- allocate(Mref_V_radius_ref(NR_REF), &
- Mref_V_density_ref(NR_REF), &
- Mref_V_vpv_ref(NR_REF), &
- Mref_V_vph_ref(NR_REF), &
- Mref_V_vsv_ref(NR_REF), &
- Mref_V_vsh_ref(NR_REF), &
- Mref_V_eta_ref(NR_REF), &
- Mref_V_Qkappa_ref(NR_REF), &
- Mref_V_Qmu_ref(NR_REF), &
- stat=ier)
- if( ier /= 0 ) call exit_MPI(myrank,'error allocating Mref_V arrays')
-
! all processes will define same parameters
call define_model_1dref(CRUSTAL)
@@ -4404,194 +4387,194 @@
104.d0 /)
Mref_V_Qmu_ref( 181 : 210 ) = (/ &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 /)
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 /)
Mref_V_Qmu_ref( 211 : 240 ) = (/ &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 /)
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 /)
Mref_V_Qmu_ref( 241 : 270 ) = (/ &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 /)
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 /)
Mref_V_Qmu_ref( 271 : 300 ) = (/ &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 /)
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 /)
Mref_V_Qmu_ref( 301 : 330 ) = (/ &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 /)
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 /)
Mref_V_Qmu_ref( 331 : 360 ) = (/ &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
- 0.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
+ 4999.d0 , &
355.d0 , &
355.d0 /)
diff --git a/src/meshfem3D/model_s362ani.f90 b/src/meshfem3D/model_s362ani.f90
index 20a0527..673d954 100644
--- a/src/meshfem3D/model_s362ani.f90
+++ b/src/meshfem3D/model_s362ani.f90
@@ -224,11 +224,11 @@
logical upper,upper_650
logical lower,lower_650
- real(kind=4), parameter :: r0=6371.
- real(kind=4), parameter :: rmoho=6371.0 - 24.4 ! subtracting the thickness here
- real(kind=4), parameter :: r670=6371. - 670. ! subtracting the thickness here
- real(kind=4), parameter :: r650=6371. - 650. ! subtracting the thickness here
- real(kind=4), parameter :: rcmb=3480.0
+ real(kind=4), parameter :: r0 = 6371.
+ real(kind=4), parameter :: rmoho = 6371.0 - 24.4 ! subtracting the thickness here
+ real(kind=4), parameter :: r670 = 6371. - 670. ! subtracting the thickness here
+ real(kind=4), parameter :: r650 = 6371. - 650. ! subtracting the thickness here
+ real(kind=4), parameter :: rcmb = 3480.0
integer :: i,nspl,nskip,nlower,nupper,iker,lstr
@@ -782,7 +782,6 @@
!-------------------------------------------------------------------------------------------------
!
-
subroutine gt3dmodl(lu,targetfile, &
maxhpa,maxker,maxcoe, &
numhpa,numker,numcoe,lmxhpa, &
@@ -911,7 +910,6 @@
!-------------------------------------------------------------------------------------------------
!
-
subroutine rd3dmodl(lu,filename,ierror, &
nmodkern,nhorpar,ityphpar, &
ihorpar,lmaxhor,ncoefhor, &
@@ -1073,16 +1071,6 @@
integer, intent(in) :: nver
integer, intent(out) :: ncon
-! Daniel Peter: original define
-!
-! real(kind=4) verlat(1)
-! real(kind=4) verlon(1)
-! real(kind=4) verrad(1)
-!
-! integer icon(1)
-! real(kind=4) con(1)
-
-! Daniel Peter: avoiding out-of-bounds errors
real(kind=4), intent(in) :: verlat(nver)
real(kind=4), intent(in) :: verlon(nver)
real(kind=4), intent(in) :: verrad(nver)
@@ -1137,7 +1125,7 @@
!-------------------------------------------------------------------------------------------------
!
-! --- evaluate perturbations in percent
+! --- evaluate perturbations
subroutine model_s362ani_subshsv(xcolat,xlon,xrad,dvsh,dvsv,dvph,dvpv)
@@ -1190,22 +1178,9 @@
call ylm(y,x,lmax,ylmcof(1,ihpa),wk1,wk2,wk3)
else if(itypehpa(ihpa) == 2) then
numcof=numcoe(ihpa)
-
-! originally called
- call splcon(y,x,numcof,xlaspl(1,ihpa), &
+ call splcon(y,x,numcof,xlaspl(1,ihpa), &
xlospl(1,ihpa),radspl(1,ihpa), &
nconpt(ihpa),iconpt(1,ihpa),conpt(1,ihpa))
-
-! making sure of array bounds
-!! note from DK DK, July 2013: the code below is correct but it is exactly the same as the code
-!! note from DK DK, July 2013: above because xlaspl(1:numcof,ihpa) is contiguous in memory, and thus
-!! note from DK DK, July 2013: using a pointer to it such as xlaspl(1,ihpa) in the call seems fine as well;
-!! note from DK DK, July 2013: and some compilers could create expensive and useless memory copies for each call
-!! note from DK DK, July 2013: with the new syntax below
-! call splcon(y,x,numcof,xlaspl(1:numcof,ihpa), &
-! xlospl(1:numcof,ihpa),radspl(1:numcof,ihpa), &
-! nconpt(ihpa),iconpt(1,ihpa),conpt(1,ihpa))
-
else
write(6,"('problem 1')")
endif
@@ -1259,7 +1234,6 @@
enddo ! --- ieval
! evaluate perturbations in vsh and vsv
-
if(ish == 1) then
vsh3drel=valu(1)+0.5*valu(2)
else if(ish == 0) then
@@ -1270,8 +1244,7 @@
enddo ! --- by ish
-! evaluate perturbations in percent
-
+! evaluate perturbations
dvsh=vsh3drel
dvsv=vsv3drel
dvph=0.55*dvsh ! scaling used in the inversion
@@ -1319,22 +1292,9 @@
call ylm(y,x,lmax,ylmcof(1,ihpa),wk1,wk2,wk3)
else if(itypehpa(ihpa) == 2) then
numcof=numcoe(ihpa)
-
-! originally called
- call splcon(y,x,numcof,xlaspl(1,ihpa), &
+ call splcon(y,x,numcof,xlaspl(1,ihpa), &
xlospl(1,ihpa),radspl(1,ihpa), &
nconpt(ihpa),iconpt(1,ihpa),conpt(1,ihpa))
-
-! making sure of array bounds
-!! note from DK DK, July 2013: the code below is correct but it is exactly the same as the code
-!! note from DK DK, July 2013: above because xlaspl(1:numcof,ihpa) is contiguous in memory, and thus
-!! note from DK DK, July 2013: using a pointer to it such as xlaspl(1,ihpa) in the call seems fine as well;
-!! note from DK DK, July 2013: and some compilers could create expensive and useless memory copies for each call
-!! note from DK DK, July 2013: with the new syntax below
-! call splcon(y,x,numcof,xlaspl(1:numcof,ihpa), &
-! xlospl(1:numcof,ihpa),radspl(1:numcof,ihpa), &
-! nconpt(ihpa),iconpt(1:maxver,ihpa),conpt(1:maxver,ihpa))
-
else
write(6,"('problem 1')")
endif
@@ -1744,10 +1704,9 @@
!
! 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)
+ real(kind=4), dimension((LMAX+1)**2) :: Y
real(kind=4), parameter :: RADIAN = 57.2957795
@@ -1766,9 +1725,7 @@
! 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)
- CALL legndr(THETA,L,L,WK1(1:L+1),WK2(1:L+1),WK3(1:L+1))
+ CALL legndr(THETA,L,L,WK1,WK2,WK3)
FAC=(1.,0.)
DFAC=CEXP(CMPLX(0.,PHI))
@@ -1797,8 +1754,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
+ real(kind=4) :: X(M+1),XP(M+1),XCOSEC(M+1)
double precision :: SMALL,sumval,COMPAR,CT,ST,FCT,COT,X1,X2,X3,F1,F2,XM,TH
More information about the CIG-COMMITS
mailing list