Subject: Array-Length phenomena (bugs?)
From: Peter Paul Smolka (smolka@uni-muenster.de)
Date: Fri May 29 1998 - 15:19:56 MDT
Dear Dr. Kluszek,
Dear ccm-users,
while porting ccm3 to another platform I got into contact with
two potentially significant phenomena (bugs). Before there are
any consequences note that runtime-systems may behave differently
(so your runs may not be affected).
Summary:
I. Insufficient length of real-header array in Subr. RDHDR
and its potential workaround.
II. Small phenomenon in subr. settau
III. Length mismatch of arry u3 in call to inidat from subr.
inital.
End of summary.
(1) sigapb has the length 2*plev+1 = 37 (with plev=18).
(2) from subr. inital subr. rdhdr is called. Here the 4th but last
array is sigapb.
(3) IN subr. rdhdr this is hedr(*), which means hedr gets
a length of 37 that is practically: hedr(37) (!)
This can be checked for example with a do loop running from 1 to
40 before reading to some value. This do-loop will cause a
controlled crash (if compiler options set in accordance) after
element 37 (for example if inserted in subr. rdhdr)
(4) As hedi(32) (integer header) is 239, not only according to
the file but also according to manual p. 125 the real-header
implicite-read attempts to read beyond 37 which causes a
subscript-range error, see code: (hedr(i),i=1,hedi(32))
Or in other words: If other runtime-systems do not crash
what is with the coefficients 33ff and/or where does the input
go (consequences for results)?
It is not excluded that some run-time systems may flood the elements in
case of subscript-range error automatically into the other arrays,
so not only sigapb is filled but also siga etc. The old IBM G0 / G1
compilers (around 1981) had this bug, which sometimes was used as
a feature as well.
End of description phenomenon I.
***********************************
My suggested solution for this (subr. RDHDR can then be left as it is,
which is practical as it is called also from elsewhere).
Additional array HDRPP in Subr. Inital which is allocated with a length
of LPP=3 * (2*PLEV+1) + 2*PLAT (no. of args deduced from length of
array in common-block comhdr
The CALT means C OLD, the CPP: Comments starting changes by me (Initials)
if (masterproc) then
C
C Attach initial dataset
C
call atchbnd(ninit ,ncdata ,'U' ,lcroot)
C
C Read initial data
C
write (6,*) 'ninit=',ninit,' ncdata=',ncdata
write (6,*) 'plenhi=',plenhi
write (6,*) 'Before call rdhdr plev=',plev, 'plat=',plat,
* ' length header=',LPP
LPP=3 * (2*PLEV+1) + 2*PLAT
ALLOCATE (HDRPP(LPP),STAT=ISTAT)
IF(ISTAT.GT.0) GOTO 2001
CALT call rdhdr(ninit ,plenhi ,plenhis ,plenhc ,plenhcs ,
CALT $ plenhr ,comhdi ,comhdc ,sigapb ,mflds ,
CALT $ mcflds ,ier )
call rdhdr(ninit ,plenhi ,plenhis ,plenhc ,plenhcs ,
$ plenhr ,comhdi ,comhdc ,HDRPP ,mflds ,
$ mcflds ,ier )
CPP now move coefficients from HDRPP into sigapb etc according
CPP to description of comhdr and manual p. 128
CPP code it as primitive as possible (please apologize the style)
CPP common /comhdr/
CPP $ sigapb(2*plev+1) ,siga(2*plev+1) ,sigb(2*plev+1) ,
CPP $ hdlat(plat) ,hdwt(plat)
nfrom=0
nend=2*plev+1
do 2002 npp=1,nend
nfrom=nfrom+1
sigapb(npp)=hdrpp(nfrom)
write (6,20021) nfrom,npp,sigapb(npp)
20021 format ('n-header=',I5,' n-sig=',I5,' sigapb=',E16.7)
2002 continue
do 2003 npp=1,nend
nfrom=nfrom+1
siga(npp)=hdrpp(nfrom)
write (6,20031) nfrom,npp,siga(npp)
20031 format ('n-header=',I5,' n-sig=',I5,' siga=',E16.7)
2003 continue
do 2004 npp=1,nend
nfrom=nfrom+1
sigb(npp)=hdrpp(nfrom)
write (6,20041) nfrom,npp,sigb(npp)
20041 format ('n-header=',I5,' n-sig=',I5,' sigb=',E16.7)
2004 continue
do 2005 npp=1,plat
nfrom=nfrom+1
hdlat(npp)=hdrpp(nfrom)
write (6,20051) nfrom,npp,hdlat(npp)
20051 format ('n-header=',I5,' n-hdl=',I5,' hdlat=',E16.7)
2005 continue
do 2006 npp=1,plat
nfrom=nfrom+1
hdwt(npp)=hdrpp(nfrom)
write (6,20061) nfrom,npp,hdwt(npp)
20061 format ('n-header=',I5,' n-hdl=',I5,' hdwt=',E16.7)
2006 continue
deallocate (hdrpp,stat=istat)
if (ier.ne.0) then
write(6,*)'INITAL:Error returned from RDHDR for record ',ier
call endrun
end if
end if
C#line 18490
By the way: Indirect deduction of sequence correctly as from overall
logic that seems to be the most likely one.
II. Subr. Settau, negative gravity-wave speeds (some)
Having done so, the results of the provided test-data-set look exactly
as they should until (including) manual p. 191 reference-temperatures
for semi-implicit scheme (e.g. headers are checked correctly, a similar
machine epsilon is computed (really a good idea what they did), reference-pressures
are exactly(!) as in handbook BUT:
Gravity-Wave Phase speeds (M/S)... are partially negative,
same with Gravity wave equiv. depths, so ccm3 aborts controlled.
Following this leads to the computation of zci/zcr and the comments
inserted by NCAR into the respective routines.
Temporary workaround: zci and zcr set to the values prescribed in the
manual om p. 192 (as the order of magnitude was correct)
with error-msg for safety (hessenburg and eigenvalue optimization
later).
III. Inidat called from inital
The ccm3 proceed to inidat called from inital.
Here another potentially significant array mismatch occurres:
Array u3, length mismatch
C#line 18703
write (6,*) 'Before call inidat'
call inidat(ps(1,1,1),u3(i1,1,j1,1),v3(i1,1,j1,1),
.......
The variables i1 and j1 have the values:
i1= 2 j1= 3 based on:
nxpt= 1 jindmx= 1
(see
parameter(nxpt = 1)
parameter(jintmx = 1)
parameter(plond = plon + 1 + 2*nxpt)
parameter(platd = plat + 2*nxpt + 2*jintmx)
parameter(plevd = plev*(3 + pcnst))
parameter(i1 = 1 + nxpt)
parameter(j1 = 1 + nxpt + jintmx)
beglat= 1 endlatex= 68 platd= 68
u3(plond,plevd,beglatex:endlatex,2)
plond= 131 plevd= 72 beg:endlatex= 1: 682
u3(i1,1,j1,1) i1= 2 j1= 3
IN Inidat however the declaration is:
C
C#line 17642 "ccm.f"
C------------------------------Arguments--------------------------------
C
C Input/Output arguments
C
real*8 ps(plond,plat,2) ! surface pressure
real*8 u3(plond,plevd,platd,2) ! u-wind
That is: plond and platd have here the normal values, which is:
higher than that of i1 and j1.
All this can be checked with own write stmts, so do NOT stop any
running simulations, especially as due to the static nature of
most fortran-arrays some runtime systems may (or may not) floodfill
the numbers somehow (correctly) as long as the size is declared once
correctly.
Sorry for any inconvenice caused, but while porting to a non-Cray
environment I put of course for safety all debugging options
on and check critical passages manually.
Anyway: The NCAR people did a absolutely impressive job especially
as the code reflects the data-processing-history of then the past
10 years (and who else beyond NCAR puts the code for public discussion?)
Most sincerely, Peter Smolka
P.S.: And maybe some minor suggestions for future versions:
Declaring all reals explicitely to REAL*8 would be nice
(not all compilers have the Autodbl option). The other changes
are then only calling function DBLE( ) in all MAX( ) and MIN( )
Function calls that have real constants, changing the sign call to
DSIGN and as far as I remember changing one AMAX1 to DMAX1.
(a world without Cray-Pointers would even be nicer....).
And for distribution:
The Cray etc. as it is and maybe in addition for the rest of the world:
And for the rest of the world: A SST-driven T42
(and maybe T21) version code, so the UNIX step can be avoided.
(that is only 2 files more, maybe in case of SOM 2 other files more).
This makes ccm easier to handle as the special code for
parallel-computing simply does not exist any more and some
architectures do the scheduling (if applicable) even automatically.
Dr. Peter P. Smolka
University Muenster
Geological Institute
Corrensstr. 24
D-48149 Muenster
Tel.: +49/251/833-3989 +49/2533/4401
Fax: +49/251/833-3968 +49/2433/4401
Try to leave the earth as you would like to meet her when you come.
This archive was generated by hypermail 2b27 : Thu Jun 01 2000 - 09:06:17 MDT