Subject: Re: Year 68 Virus
From: Jim Rosinski (rosinski@bearmtn-e0.cgd.ucar.edu)
Date: Mon Sep 14 1998 - 15:44:13 MDT
CCM Users;
We have come up with a fix for the "68 year" problem originally brought up by
Kevin Werner. The problem he encountered (on a SUN) was that 32-bit integer
computations would overflow at approximately year 68 of the simulation in the
routine which computes calendar day information. The fix involves a few
lines of changes to the routine calendr.F in CCM3.6 (and CCM3.6.6). A few of
the computations are now forced to be done using 64-bit integer arithmetic.
We have tested the fix on all platforms on which the CCM is supported and
gotten bit-for-bit answers as compared to the old code (for a one-day run).
An offline driver test was run out to 7000 years of simulation and no
problems were found.
The fix we are providing applies to CCM3.6 and CCM3.6.6. The modified code
is appended to the end of this mail message, and is also available on the
CCM3 news page: http://www.cgd.ucar.edu/cms/ccm3/news/
For those who wish to retrofit the changes into
CCM3.2, you will need to change caldyi.F in addition to calendr.F.
Please let us know if you encounter further problems.
Regards,
Jim Rosinski
CCM Core Group
-----------------------------Cut Here-----------------------------------------
#include <misc.h>
subroutine calendr(nstep ,dtime ,mdbase ,msbase ,mbdate ,
& mbsec ,mdcur ,mscur ,mcdate ,mcsec ,
& calday )
!-----------------------------------------------------------------------
!
! Compute current date and day information for history tape header.
! Compute current julian day (including fraction) for sun angle
! computations and time interpolation of boundary dataset information.
! One year is defined as 365 days.
!
! Computational notes:
!
! 86400 is the number of seconds in 1 day.
!
! Dividing an integer by 10**n has the effect of right-shifting the
! decimal digits n positions (ex: 861231/100 = 008612).
!
! mod(integer,10**n) has the effect of extracting the rightmost
! n decimal digits of the integer (ex: mod(861231,10000) = 1231).
!
!---------------------------Code history--------------------------------
!
! Original version: J. Rosinski
! Standardized: J. Rosinski, June 1992
! Reviewed: J. Kiehl, B. Boville, August 1992
! Reviewed: J. Kiehl, April 1996
! Reviewed: B. Boville, April 1996
!
!-----------------------------------------------------------------------
!
! $Id: calendr.F,v 1.5 1998/09/14 20:46:44 rosinski Exp $
! $Author: rosinski $
!
implicit none
!
! Input arguments
!
integer nstep ! current time step (0, 1, ...)
real dtime ! length of time step (seconds)
integer mdbase ! base day of run (e.g., 0)
integer msbase ! base seconds of base day (e.g., 0)
integer mbdate ! base date (yyyymmdd format) of run (e.g., 000901)
integer mbsec ! base seconds of base date (e.g., 0)
!
! Output arguments
!
integer mdcur ! current day (0, 1, ...)
integer mscur ! current seconds of current day (0, ..., 86400)
integer mcdate ! current date (yyyymmdd format) (e.g., 021105)
integer mcsec ! current seconds of current date (0, ..., 86400)
real calday ! current julian day including fraction
!
!---------------------------Local workspace-----------------------------
!
integer*8 nstep8 ! current time step (0, 1, ...)
integer*8 nsecs ! run time in seconds
integer ndays ! ndays + nyears = number of day changes since start
integer nyears ! ndays + nyears = number of day changes since start
integer mcyear ! current year of current date (0-99)
integer mcmnth ! current month of current date (1-12)
integer mcday ! current day of current date (1-31)
integer jday ! Julian day (1-365)
integer mbmnth ! base month (1-12) validity check
integer mbday ! base day of base date (1-31) validity check
integer ndm(12) ! number of days per month
integer jdcon(12) ! convert month index to julian day
integer*8 nspdy ! number of seconds per day
parameter (nspdy=86400)
!
save ndm,jdcon
data ndm/31,28,31,30,31,30,31,31,30,31,30,31/
data jdcon/0,31,59,90,120,151,181,212,243,273,304,334/
!
!-----------------------------------------------------------------------
!
nstep8 = nstep
!
! Check validity of input data
!
mbmnth = mod(mbdate,10000)/100
mbday = mod(mbdate,100)
if (mbmnth.lt.1 .or. mbmnth.gt.12) then
write(6,*)' CALENDR: Invalid base month input:',mbmnth
call endrun
end if
if (mbday.lt.1 .or. mbday.gt.ndm(mbmnth)) then
write(6,*)' CALENDR: Invalid base day of base date input:',mbday
call endrun
end if
if (msbase.lt.0 .or. msbase.ge.nspdy) then
write(6,*)' CALENDR: Invalid base seconds(msbase):',msbase
call endrun
end if
if (mbsec.lt.0 .or. mbsec.ge.nspdy) then
write(6,*)' CALENDR: Invalid base seconds(mbsec):',mbsec
call endrun
end if
!
nsecs = nstep8*nint(dtime)
!
! First: current day, seconds
!
mdcur = mdbase + (nsecs+msbase)/nspdy
mscur = mod((nsecs+msbase),nspdy)
!
! Next: current date, seconds.
! Uncommenting the next line will have the effect of modifying
! nsecs to include base day and base seconds. This is the way
! the current date was computed in CCM1.
!
! nsecs = nsecs + mdbase*nspdy + msbase
ndays = (nsecs+mbsec)/nspdy
nyears = ndays/365
ndays = mod(ndays,365)
!
! Current seconds of current date
!
mcsec = mod(nsecs+mbsec,nspdy)
!
! Initialize current year, month, day.
!
mcyear = mbdate/10000 + nyears
mcmnth = mod(mbdate,10000)/100
mcday = mod(mbdate,100) + ndays
!
! Now loop through months, converting yyyy, mm, and ddd to yyyymmdd.
! ex: 791235 becomes 800104. 190001370 becomes 19010105.
!
10 if (mcday.gt.ndm(mcmnth)) then
mcday = mcday - ndm(mcmnth)
mcmnth = mcmnth + 1
if (mcmnth.eq.13) then ! add a year
mcyear = mcyear + 1
mcmnth = 1
end if
go to 10
end if
mcdate = mcyear*10000 + mcmnth*100 + mcday
!
! Convert current month, day, seconds to Julian day + fraction. Multiply
! nspdy by floating point 1 in order to ensure real*8 arithmetic. As of
! 9/14/98 this only makes a difference on Linux.
!
jday = jdcon(mcmnth) + mcday
calday = float(jday) + float(mcsec)/(1.*nspdy)
return
end
This archive was generated by hypermail 2b27 : Thu Jun 01 2000 - 09:06:55 MDT