ICOADS Web information page (Thursday, 05-Jun-2014 20:10:45 UTC):

Software demo: {lmrlib}


This Fortran library, {lmrlib}, was mostly developed by 2000, as the ICOADS project increasingly recognized the need for a central shared resource of reusable code, containing well documented and thoroughly validated subroutines and functions for historical units conversions and similar requirements encountered during translation of data into the common ICOADS formats.

The name "lmr" refers (in retrospect not the best choice) to the now-obsolete Long Marine Report (binary) format (documentation here), which still however retains some structural characteristics of our current International Maritime Meteorological Archive (IMMA) format (documentation here).

In addition, Philip Brohan (UK Met Office) has begun developing a shared online library for the code he uses for IMMA analysis, which is available at https://github.com/oldweather/IMMA ("pull requests" welcome, ref. http://oss-watch.ac.uk/resources/pullrequest). This library includes a Perl implementation of {lmrlib}, which together with its test suite is believed to be functionally equivalent to the original Fortran {lmrlib}.

Following is the full Fortran source code of {lmrlib}, plus its enclosing Unix script commands, with a few portions of the code highlighted, or added explanatory comments interspersed, in red:

cat > lmrlib.f <<\EOR
C=============================================================================C
C Comprehensive Ocean-Atmosphere Data Set (COADS):   Fortran 77 Program+Shell C
C Filename:level: lmrlib:01G                                 17 November 2004 C
C Function: Tools to assist conversions into LMR6   Author: S.Woodruff et al. C
C=============================================================================C
C Software Revision Information (previous version: 17 Aug 2000, level 01F):
C Remove print statements and remaining stop from {ixdtnd,rxnddt}.
C-----------------------------------------------------------------------3456789
C Software documentation for the (modifiable) library test routine {libtst},
C and for the remaining (invariant) user-interface routines in the library.
C
C Functionality:  This is a library of tools to assist conversions from other
C formats into LMR6, whose functions are individually described in the comments
C at the beginning of each subroutine or function.  The following naming and
C other conventions have been applied:
C     1) Subprogram (function or subroutine) names begin with a letter
C     indicating their type or function:
C          a) "f" for type real (floating point)
C          b) "i" for type integer
C          c) "r" for subroutines
C          d) "t" for test and other software development subroutines.
C     2) The second letter of subprogram names indicates:
C          a) "x" for subprograms for units and related conversions.  Units
C          "to" and "from" generally compose the remainder of the name as
C          2-letter abbreviations in that order (e.g., ixbfkt is an integer
C          function that maps Beaufort numbers into speed midpoints in kts).
C          b) "w" for subprograms for corrections.  Presently these are
C          limited to barometric corrections, with "bp" in letters 3-4.
C          c) "p" for subprograms producing only printed results.
C          d) "m" for miscellaneous subprograms.
C     3) Each subprogram includes specific declarations of the type of each
C     input and output argument (implicit and default typings are not used
C     for external communication; they may be used internally to subprograms).
C     4) Input/output arguments (or function returns) that are most naturally
C     represented as integers, are typed as such.  Users will need to make
C     transformations from/to real variables, as needed.
C     5) Errors detected within the library generally result in a printed
C     statement including the name of the library routine involved, and then
C     a stop.  Exceptions are {ixdtnd,rxnddt}, which return -1 instead of stopping,
C     and {ix32dd,ixdcdd}, which return a user-supplied missing value.

! We found error handling (above) to be a tricky issue. Having the subroutine
! issue a STOP in some production situations could be problematic. However
! developing other means of exception handling could introduce a lot more
! complexity. 

C     6) Portions of the library can be exhaustively tested by running the
C     library test routine {libtst}, which should yield no output (unless
C     errors are printed).  
C Contents:  Following are the routines included, and their broader groupings:
C      barometric conversions:
C           {fxmmmb}     millimeters Hg to millibars
C           {fxmbmm}     millibars to millimeters Hg
C           {fxeimb}     inches (English) Hg to millibars
C           {fxmbei}     millibars to inches (English) Hg
C           {fxfimb}     inches (French) Hg to millibars
C           {fxmbfi}     millibars to inches (French) Hg, plus
C                        entry points {fxfim0,fxfim1} used by {tpbpfi}
C           {fwbptf}     correction value for temperature (Fahrenheit)
C           {fwbptc}     correction value for temperature (Celsius)
C           {fwbptg}     correction value for temperature (generalized)
C           {fwbpgv}     correction value for gravity
C           {tpbpfi}     test: print {fxfimb,fxfim0,fxfim1,fxmbfi} results
C           {tpbpt1}     test: print {fwbptg} results (UKMO, 1969; List, 1966)
C           {tpbpg1}     test: print {fwbpgv} results (List, 1966, Table 47B)
C           {tpbpg2}     test: print {fwbpgv} gmode intercomparisons
C           {tpbpmb}     test: print results for e.g. mercurial corrections
C      cloud conversions:
C           {ixt0ok}     tenths (of sky clear  ) to oktas (of sky covered)
C           {ixt1ok}     tenths (of sky covered) to oktas (of sky covered)
C           {tpt0t1}     test: print {ixt0ok,ixt1ok} results
C      temperature conversions:
C           {fxtftc}     Fahrenheit to Celsius
C           {fxtctf}     Celsius to Fahrenheit
C           {fxtktc}     Kelvins to Celsius
C           {fxtctk}     Celsius to Kelvins
C           {fxtrtc}     Reaumur to Celsius
C           {fxtctr}     Celsius to Reaumur
C           {tptcrf}     test: print {fxtftc,fxtctf,fxtrtc,fxtctr} results

! We highlight here that this library does not (yet anyway) include any
! humidity/related calculations, e.g. of dew point temperature. Those important
! requirements are however discussed in the separate Software demo {hum_calc}.

C      wind conversions:
C           {fxuvdd}     vector (u,v) components to direction (degrees)
C           {fxuvvv}     vector (u,v) components to velocity
C           {rxdvuv}     direction and speed to vector (u,v) components
C           {fxktms}     knots to m/s (ref. international nautical mile)
C           {fxmskt}     m/s to knots (ref. international nautical mile)
C           {fxk0ms}     knots to m/s (ref. US nautical mile)
C           {fxmsk0}     m/s to knots (ref. US nautical mile)
C           {fxk1ms}     knots to m/s (ref. Admiralty nautical mile)
C           {fxmsk1}     m/s to knots (ref. Admiralty nautical mile)
C           {ixbfkt}     0-12 Beaufort number to knots (WMO code 1100)
C           {fxbfms}     0-12 Beaufort number to m/s (WMO code 1100)
C           {ix32dd}     32-point direction abbreviation into code and degrees
C           {ixdcdd}     32-point direction code into degrees
C           {txdvuv}     test: {fxuvdd,fxuvvv,rxdvuv}
C           {tpktms}     test: print conversion factors
C      time conversions:
C           {rxltut}     local standard hour and "Julian" day into UTC
C           {ixdtnd}     date to days ("Julian") since 1 Jan 1770
C           {rxnddt}     days ("Julian") since 1 Jan 1770 to date
C           {txltut}     test: {rxltut}
C           {txdtnd}     test: {ixdtnd,rxnddt}
C      miscellaneous:
C           {rpepsi}     print machine epsilon
C           {rpepsd}     double precision version of {rpepsi}
C           {imrnde}     round to nearest even integer
C Machine dependencies:  None known.
C For more information: See <soft_info> and <soft_lmr> (electronic documents).
C-----------------------------------------------------------------------3456789
c      program libtst
      subroutine libtst

! This user-adjustable subroutine invokes all internal testing "t" routines,
! which is recommended e.g. if the library is ported to a different computer
! system. Also, the subroutine statement can be commented out and program
! statement activated, to turn this into a main program+subprograms, for stand-
! alone testing convenience. Note that "testing" in many cases refers to
! replication of published values (e.g. from the Smithsonian Met. Tables).

c-----temporary code for testing purposes
      implicit integer(a-e,g-z)
c-----barometric conversions:  print
      call tpbpfi
      call tpbpt1
      call tpbpg1(1)
      call tpbpg2(.false.)
      call tpbpmb
c-----cloud conversions:       print
      call tpt0t1
c-----temperature conversions: print
      call tptcrf
c-----wind conversions:        test
c      call txdvuv
c-----                         print
      call tpktms
c-----time conversions:        test
c      call txltut
c      call txdtnd
c-----miscellaneous:           print
      call rpepsi
      call rpepsd
      end
C=============================================================================C
C WARNING:  Code beyond this point should not require any modification.       C
C=============================================================================C
C=======================================================================3456789
c-----barometric conversions--------------------------------------------3456789
C=======================================================================3456789
      real function fxmmmb(mm)
c-----Convert barometric pressure in (standard) millimeters of mercury (mm)
c     to millibars (hPa), e.g., fxmmmb(760.) = 1013.25 (one atmosphere)
c     (List, 1966, p. 13).
c-----sdw, 23 Nov 1999.
c-----sdw, 25 Feb 2000: changes to comments.
c     References:
c     List, R.J., 1966: Smithsonian Meteorological Tables.
c           Smithsonian Institution, Washington, DC, 527 pp.
c     WMO (World Meteorological Organization), 1966: International
c           Meteorological Tables, WMO-No.188.TP.94.
      real mm
c-----factor from List (1966), p. 13 and Table 11; also in WMO (1966).
      fxmmmb = mm * 1.333224
      return
      end
C-----------------------------------------------------------------------3456789
      real function fxmbmm(mb)
c-----Convert barometric pressure in millibars (hPa; mb) to (standard)
c     millimeters of mercury.  Numerical inverse of {fxmmmb} (see for
c     background).  Note: This method yields better numerical agreement
c     in cross-testing against that routine than the factor 0.750062.
c-----sdw, 26 Jan 2000.
c-----sdw, 25 Feb 2000: changes to comments.
      real mb
      fxmbmm = mb / 1.333224
      return
      end
C-----------------------------------------------------------------------3456789
      real function fxeimb(ei)
c-----Convert barometric pressure in (standard) inches (English) of
c     mercury (in) to millibars (hPa), e.g., fxeimb(29.9213) = 1013.25
c     (one atmosphere) (List, 1966, p. 13).
c-----sdw, 26 Jan 2000.
c-----sdw, 25 Feb 2000: changes to comments.
c     References:
c     List, R.J., 1966: Smithsonian Meteorological Tables.
c           Smithsonian Institution, Washington, DC, 527 pp.
c     WMO (World Meteorological Organization), 1966: International
c           Meteorological Tables, WMO-No.188.TP.94.
      real ei
c-----factor from List (1966), Table 9.  Note: a slightly different factor 
c     33.8639 appears also on p. 13 of List (1966), and in WMO (1966).  Tests 
c     (32-bit Sun f77) over a wide range of pressure values (25.69"-31.73",
c     approximately equivalent to ~870-1074.6 mb) indicated that the choice
c     of constant made no numeric difference when data were converted to mb
c     and then rounded to tenths, except for two cases of 0.1 mb difference
c     (25.79" = 873.3/.4 mb; and 26.23" = 888.2/.3 mb).  If 1 mm = 1.333224
c     mb and 1" = 25.4 mm, then 25.4mm = 33.86389 to 7 significant digits.
      fxeimb = ei * 33.86389
      return
      end
C-----------------------------------------------------------------------3456789
      real function fxmbei(mb)
c-----Convert barometric pressure in millibars (hPa; mb) to (standard)
c     inches (English) of mercury.  Numerical inverse of {fxeimb} (see for
c     background).  Note: This method yields better numerical agreement
c     in cross-testing against that routine than the factor 0.0295300.
c-----sdw, 26 Jan 2000.
      real mb
      fxmbei = mb / 33.86389
      return
      end
C-----------------------------------------------------------------------3456789
      real function fxfimb(fi)
c-----Convert barometric pressure in inches (French) of mercury (fi) to
c     millibars (hPa).  Paris, instead of French, inches are referred
c     to in Lamb (1986), but these appear to be equivalent units.  Note:
c     data in lines (twelve lines per inch) or in inches and lines need
c     to be converted to inches (plus any decimal fraction).  Entry points
c     {fxfim0,fxfim1}, which are called by {tpbpfi} (see for background)
c     are not recommended for use in place of {fxfimb}.
c-----sdw, 26 Jan 2000.
c     References:
c     IMC (International Meteorological Committee), 1890: International
c           Meteorological Tables, published in Conformity with a Resolution
c           of the Congress of Rome, 1879.  Gauthier-Villars et Fils, Paris.
c     Lamb, H.H., 1986: Ancient units used by the pioneers of meteorological
c           instruments.  Weather, 41, 230-233.
      real fi
c-----factor for conversion of French inches to mm (IMC, 1890, p. B.2);
c     mm are then converted to mb via {fxmmmb}
      fxfimb = fxmmmb(fi * 27.069953)
      return
c-----factor for conversion of French inches to English inches (ibid.);
c     inches are then converted to mb via {fxeimb}
      entry fxfim0(fi)
      fxfim0 = fxeimb(fi * 1.0657653)
      return
c-----direct conversion factor to mb from Lamb (1986), Table 1 footnote.
      entry fxfim1(fi)
      fxfim1 = fi * 36.1
      return
      end
C-----------------------------------------------------------------------3456789
      real function fxmbfi(mb)
c-----Convert barometric pressure in millibars (hPa; mb) to inches (French)
c     of mercury.  Numerical inverse of {fxfimb} (see for background).
c-----sdw, 26 Jan 2000.
      real mb
      fxmbfi = fxmbmm(mb) / 27.069953
      return
      end
C-----------------------------------------------------------------------3456789
      real function fwbptc(bp,tc)
c-----Correction value of barometric pressure (in mm or mb; standard
c     temperature of scale 0C) (bp) for temperature in Celsius (tc)
c     (see {fwbpgt} for additional background).
c-----sdw, 23 Nov 1999.
c-----sdw, 25 Feb 2000: changes and additions to comments.
c     Reference:
c     List, R.J., 1966: Smithsonian Meteorological Tables.
c           Smithsonian Institution, Washington, DC, 527 pp.
      real bp,tc,m,l
c-----constants m and l from List (1966), p. 136.
      data m/0.0001818/,l/0.0000184/
      fwbptc = -bp * ( ((m-l)*tc) / (1.+(m*tc)) )
      return
      end
C-----------------------------------------------------------------------3456789
      real function fwbptf(bp,tf)
c-----Correction value of barometric pressure (in inches; standard
c     temperature of scale 62F) (bp) for temperature in Fahrenheit (tf)
c     (see {fwbpgt} for additional background).
c-----sdw, 23 Nov 1999.
c-----sdw, 25 Feb 2000: changes and additions to comments.
c     Reference:
c     List, R.J., 1966: Smithsonian Meteorological Tables.
c           Smithsonian Institution, Washington, DC, 527 pp.
      real bp,tf,m,l
c-----constants m and l from List (1966), p. 137.
      data m/0.000101/,l/0.0000102/
      fwbptf = -bp * ( ((m*(tf-32.))-(l*(tf-62.))) / (1.+m*(tf-32.)) )
      return
      end
C-----------------------------------------------------------------------3456789
      real function fwbptg(bp,t,u)
c-----Correction value (generalized) of barometric pressure (bp) for
c     temperature (t), depending on units (u):
c                                          standard temperature:
c          u  bp          t           of scale (ts)  of mercury (th)
c          -  ----------  ----------  -------------  -------------------
c          1  mm or mb    Celsius      0C             0C
c          2  Eng. in.    Fahrenheit  62F (16.667C)  32F (0C) (pre-1955)    
c          3  Eng. in.    Fahrenheit  32F (0C)       32F (0C) (1955-)
c          4  French in.  Reaumur     13R (16.25C)    0R (0C)

! Development of these functions for early barometric corrections led to some
! questions, not well resolved yet, about early data units and exactly how they
! should be handled, e.g. red highlighted comments here. 

c     The returned {fwbptg} value is in the same units as, and is to be
c     added to, bp.  Establishment of 0C/32F as the standard temperature
c     for both scale and mercury as implemented under u=1 and 3 became
c     effective 1 Jan 1955 under WMO International Barometer Conventions
c     (WBAN, 12 App.1.4.1--2; see also WMO, 1966 and UKMO, 1969).  List
c     (1966), p. 139 states that "the freezing point of water is universally
c     adopted as the standard temperature of the mercury, to which all
c     readings are to be reduced," but for English units uses only 62F for
c     the standard temperature of the scale.  Note: Results under u=4, and
c     the utilized settings of constants l and m, have not been verified
c     against published values, if available.  IMC (1890, p. B.24) states
c     that in "old Russian barometer readings expressed in English half lines
c     (0.05 in) the mercury and the scale were set to the same temperature
c     62F."  UKMO (1969, p. 5) states that for Met. Office barometers prior
c     to 1955 reading in millibars the standard temperature was 285K (12C).
c     This routine does not handle these, or likely other historical cases.
c-----sdw, 23 Nov 1999.
c-----sdw, 25 Feb 2000: Add u=3-4, plus changes and additions to comments.
c     References:
c     IMC (International Meteorological Committee), 1890: International
c           Meteorological Tables, published in Conformity with a Resolution
c           of the Congress of Rome, 1879.  Gauthier-Villars et Fils, Paris.
c     List, R.J., 1966: Smithsonian Meteorological Tables.
c           Smithsonian Institution, Washington, DC, 527 pp.
c     UKMO (UK Met. Office), 1969: Marine Observer's Handbook (9th ed.).
c           HMSO, London, 152 pp.
c     US Weather Bureau, Air Weather Service, and Naval Weather Service,
c           1963: Federal Meteorological Handbook No. 8--Manual of Barometry
c           (WBAN), Volume 1 (1st ed.).  US GPO, Washington, DC.
c     WMO (World Meteorological Organization), 1966: International
c           Meteorological Tables, WMO-No.188.TP.94.
      real bp,t,m(4),l(4),ts(4),th(4)
      integer u
c-----constants ts and th are from List (1966), pp. 136-137 (u=1-2); WBAN
c     12 App.1.4.1--3 (u=3); and IMC (1890), p. B.24 (u=4).
      data ts/0.0,62.,32.,13./,th/0.0,32.,32.,0.0/
c-----constants m and l are from List (1966), pp. 136-137 (u=1-3) and WBAN,
c     pp. 5-4 and 5-5 (for metric and English units).  For u=4, the u=1
c     constants were multiplied by 5/4 (after List, 1966, p. 137).
      data m/0.0001818, 0.000101, 0.000101, 0.000227/
      data l/0.0000184,0.0000102,0.0000102,0.0000230/
c-----test u for valid range
      if(u.lt.1.or.u.gt.4) then
         print *,'fwbptg error. invalid u=',u
         stop
      endif
      fwbptg = -bp * ( ((m(u)*(t-th(u)))-(l(u)*(t-ts(u))))
     +               / (1.+(m(u)*(t-th(u)))) )
      return
      end
C-----------------------------------------------------------------------3456789
      real function fwbpgv(bp,rlat,gmode)
c-----Correction value of barometric pressure (bp) for gravity depending on
c     latitude (rlat), with constants set depending on gmode (for COADS, we
c     adopt gmode=1 for 1955-forward, and gmode=2 for data prior to 1955):
c           g1 (equation 1)   g2 (equation 2)   Comment
c           ---------------   ---------------   -----------------------------
c     1 =          g45               g0         yields List (1966), Table 47B
c     2 =          g0                g0         follows GRAVCOR (pre-1955)
c     3 =          g45               g45        (of unknown utility)
c     The returned {fwbpgv} value is in the same units as, and is to be added
c     to, bp (units for bp are unspecified; Table 47B has columns for inches,
c     millimeters, and millibars).  Usage of g0 and g45 as implemented under
c     gmode=1 became effective 1 Jan 1955 under WMO International Barometer
c     Conventions: g45 is a "best" estimate of acceleration of gravity at 45
c     deg latitude and sea level, and g0 is the value of standard (normal)
c     gravity "to which reported barometric data in mm or inches of mercury
c     shall refer, but it does not represent the value of gravity at latitude
c     45 deg, at sea level" (WBAN, 12 App.1.4.1--2; see also List, 1966, pp.
c     3-4, and WMO, 1966).  For example, UK Met. Office MK I (MK II) barometers
c     issued before (starting) 1 January 1955 were graduated to read correctly
c     when the value of gravity was g45 (g0) (UKMO, 1969).  As shown by test
c     routines {tpbpg1,tpbpg2}, gmode=2 and 3 yield virtually the same results.
c-----sdw, 4 Dec 1999 (after GRAVCOR, 13 Jan 1999 e-mail fr T. Basnett, UKMO).
c-----sdw, 2 Feb 2000: changes and additions to comments.
c-----sdw, 25 Feb 2000: changes and additions to comments.
c     References:
c     List, R.J., 1966: Smithsonian Meteorological Tables.
c           Smithsonian Institution, Washington, DC, 527 pp.
c     UKMO (UK Met. Office), 1969: Marine Observer's Handbook (9th ed.).
c           HMSO, London, 152 pp.
c     US Weather Bureau, Air Weather Service, and Naval Weather Service,
c           1963: Federal Meteorological Handbook No. 8--Manual of Barometry
c           (WBAN), Volume 1 (1st ed.).  US GPO, Washington, DC.
c     WMO (World Meteorological Organization), 1966: International
c           Meteorological Tables, WMO-No.188.TP.94.
      integer gmode
      real bp,rlat,rlatr,pi,g45,g0,g1,g2,a,b,c
      data pi/3.14159265358979323846264338327950288/
c-----g45 from List (1966), p. 488 ("best" sea-level gravity at latitude 45)
      data g45/980.616/
c-----g0  from List (1966), p. 200 ("standard" acceleration of gravity)
      data  g0/980.665/
c-----check latitude 
      if(rlat.lt.-90..or.rlat.gt.90.) then
         print *,'fwbpgv error. invalid rlat=',rlat
         stop
      endif
c-----check gmode, and set g1 and g2
      if(gmode.eq.1) then
         g1 = g45
         g2 = g0
      else if(gmode.eq.2) then
         g1 = g0
         g2 = g0
      else if(gmode.eq.3) then
         g1 = g45
         g2 = g45
      else
         print *,'fwbpgv error. invalid gmode=',gmode
         stop
      endif
c-----convert degrees to radians
      rlatr  = rlat * (pi/180.)
c-----List (1966), p. 488, equation 1 (c is the local acceleration of gravity)
      a      =      0.0000059 * (cos(2.0*rlatr)**2)
      b      = 1. - 0.0026373 *  cos(2.0*rlatr)
      c      = g1 * (a + b)
c-----List (1966), p. 202, equation 2
      fwbpgv = ((c - g2)/g2) * bp
      return
      end
C-----------------------------------------------------------------------3456789
      subroutine tpbpfi
c-----Test: print {fxfimb,fxfim0,fxfim1,fxmbfi} results.  Using as input
c     the Paris lines, calculate the whole mb, column of Table 1 of Lamb
c     (1986) using the three solutions.  Calling {fxmbfi} inverts {fxfimb}.
c     The output results (32-bit Sun f77) from {fxfimb} provided slightly
c     better agreement with Table 1 than {fxfim0} (in both cases differences
c     never exceeded +0.1 mb), whereas {fxfim1} agreed less well (differences
c     up to +0.4).  Table 1 was always less than or equal to other results.
c-----sdw, 26 Jan 2000.
c-----sdw, 25 Feb 2000: change to comments.
c     Reference:
c     Lamb, H.H., 1986: Ancient units used by the pioneers of meteorological
c           instruments.  Weather, 41, 230-233.
      implicit integer(a-e,g-z)
      write(*,'(/,"tpbpfi output:")')
      write(*,1)
 1    format('   Paris:      millibars:               French inches:',/
     +      ,'------------   ----------------------   --------------',/
     +      ,'lines inches   fxfimb  fxfim0  fxfim1   fxfimb')
c-----lines is Paris lines from Table 1, transformed to fi0 real inches; fmb,
c     fm0, and fm1 contain the three solutions of converting French inches
c     to mb; and fi1 holds the inverse calculation of mb to French inches.
      do 100 lines=310,349
         fi0 = float(lines)/12.0
         fmb = fxfimb(fi0)
         fm0 = fxfim0(fi0)
         fm1 = fxfim1(fi0)
         fi1 = fxmbfi(fmb)
         write(*,2) lines,fi0,fmb,fm0,fm1,fi1
 2       format(i5,1x,f6.3,1x,3(2x,f6.1),3x,f6.3)
 100  continue
      return
      end
C-----------------------------------------------------------------------3456789
      subroutine tpbpt1
c-----Test routine for {fwbptg}.  Re-tabulate Tables 1-2 of UKMO (1969),
c     and example pages from Tables 45 and 46 of List (1966).  Agreement
c     with UKMO (1969) is approximate (different formulae may have been
c     used to construct the UKMO tables).  Exact (or near exact) numeric
c     agreement was found with List (1966): Table 45 example pages; Table
c     46 p. 166, and the first column of pp. 177 and 179.  Note: Following
c     p. 166, Table 46, except for its first column, appears to be incorrect
c     and/or incorrectly labelled.  In the implementation of this routine,
c     e.g., the p. 177 column label "0.8" under 760 mm was assumed to mean
c     768 mm instead of 760.8 mm.  But all the Table 45 numbers except the
c     first column of each page were still low compared to the re-tabulated
c     results.
c-----sdw, 25 Feb 2000.
c     References:
c     List, R.J., 1966: Smithsonian Meteorological Tables.
c           Smithsonian Institution, Washington, DC, 527 pp.
c     UKMO (UK Met. Office), 1969: Marine Observer's Handbook (9th ed.).
c           HMSO, London, 152 pp.
      implicit integer(a-e,g-z)
      real bi(11),bj1(10),bj2(10),bm1(10),bm2(10),bm3(10)
c-----pressure values in in for UKMO (1969) Tables 1-2.
      data bi/26.,26.5,27.,27.5,28.,28.5,29.,29.5,30.,30.5,31./
c-----pressure values in in for List (1966) Table 45, p. 151.
      data bj1/26.,26.2,26.4,26.6,26.8,27.,27.2,27.4,27.6,27.8/
c-----pressure values in in for List (1966) Table 45, p. 157.
      data bj2/28.,28.2,28.4,28.6,28.8,29.,29.2,29.4,29.6,29.8/
c-----pressure values in mm for List (1966) Table 46, p. 166 (1st 10 cols)
      data bm1/440.,450.,460.,470.,480.,490.,500.,510.,520.,530./
c-----pressure values in mm for List (1966) Table 46, p. 177.
      data bm2/760.,762.,764.,766.,768.,770.,772.,774.,776.,778./
c-----pressure values in mm for List (1966) Table 46, p. 179.
      data bm3/800.,802.,804.,806.,808.,810.,812.,814.,816.,818./
      write(*,'(/,"tpbpt1 output:")')
      write(*,1)
 1    format('UKMO (1969) Tables 1 (u=2) and 2 (u=3):')
c-----for {fwbptg} u=2 and 3
      do 300 u=2,3
         write(*,2) u,bi
 2       format('therm (F)  u=',i1,20x,' Barometer Reading (inches)',/
     +         ,'---------     ',11(1x,f6.1))
c-----for each whole deg F temperature value (Table 2 goes only to 90)
         do 200 jt=40,100
            write(*,3) jt,(fwbptg(bi(j),float(jt),u),j=1,11)
 3          format(i6,8x,11(1x,f6.3))
 200     continue
 300  continue
c-----for {fwbptg} u=2
      u=2
      write(*,4)
 4    format('List (1966) Table 45, pp. 151 and 157 (u=2):')
      write(*,5) u,bj1
 5    format('therm (F)  u=',i1,20x,' Barometer Reading (inches)',/
     +      ,'---------     ',10(1x,f6.1))
c-----for each half deg F temperature value (p. 151)
      jj = 0
      do 400 jt=505,750,5
         jj = jj + 1
         write(*,6) float(jt)/10.
     +             ,(fwbptg(bj1(j),float(jt)/10.,u),j=1,10)
 6       format(f6.1,8x,10(1x,f6.3))
         if(mod(jj,5).eq.0.and.jt.ne.750) print *,' '
 400  continue
c-----for each half deg F temperature value (p. 157)
      write(*,5) u,bj2
      jj = 0
      do 500 jt=755,1000,5
         jj = jj + 1
         write(*,6) float(jt)/10.
     +             ,(fwbptg(bj2(j),float(jt)/10.,u),j=1,10)
         if(mod(jj,5).eq.0.and.jt.ne.1000) print *,' '
 500  continue
c-----for {fwbptg} u=1
      u=1
      write(*,7)
 7    format('List (1966) Table 46 pp. 166, 177, and 179 (u=1):')
      write(*,8) u,bm1
 8    format('therm (C)  u=',i1,20x,' Barometer Reading (mm)',/
     +      ,'---------     ',10(1x,f6.1))
c-----for each whole deg C temperature value (p. 166)
      jj = 0
      do 600 jt=0,35
         jj = jj + 1
         write(*,9) jt,(fwbptg(bm1(j),float(jt),u),j=1,10)
 9       format(i6,8x,10(1x,f6.2))
         if(mod(jj,5).eq.0) print *,' '
 600  continue
c-----for each whole deg C temperature value (p. 177)
      write(*,8) u,bm2
      jj = 0
      do 700 jt=0,45
         jj = jj + 1
         write(*,9) jt,(fwbptg(bm2(j),float(jt),u),j=1,10)
         if(mod(jj,5).eq.0) print *,' '
 700  continue
c-----for each whole deg C temperature value (p. 179)
      write(*,8) u,bm3
      jj = 0
      do 800 jt=0,45
         jj = jj + 1
         write(*,9) jt,(fwbptg(bm3(j),float(jt),u),j=1,10)
         if(mod(jj,5).eq.0) print *,' '
 800  continue
      end
C-----------------------------------------------------------------------3456789
      subroutine tpbpg1(gmode)
c-----Test routine for {fwbpgv}.  Re-tabulate Table 47B from List (1966);
c     For reasons discussed in {fwbpgv}, gmode=1 yields good agreement
c     with Table 47B (using 32-bit Sun f77, there were nine differences
c     of no more than 0.1 pressure unit).  Table 47B calculated using
c     gmode=2 agreed almost perfectly (one 0.1 difference) with gmode=3,
c     but there were extensive differences versus gmode=1.
c-----sdw, 2 Feb 2000
c     Reference:
c     List, R.J., 1966: Smithsonian Meteorological Tables.
c           Smithsonian Institution, Washington, DC, 527 pp.
      implicit integer(a-e,g-z)
      integer gmode
      integer lat1,lat2,nin,nmm,nmb
      real uin(5),umm(6),umb(4),temp(6)
      data lat1,lat2/90,0/
c-----number of tables values for each units selection
      data nin,nmm,nmb/5,6,4/
c-----table values: inches
      data uin/  26.,  27.,  28.,  29.,  30./
c-----table values: millimeters
      data umm/ 680., 700., 720., 740., 760., 780./
c-----table values: millibars
      data umb/ 900., 950.,1000.,1050./
      write(*,'(/,"tpbpg1 output:")')
c-----inches: write header
      write(*,1) 'inches:     ',gmode,(nint(uin(j)),j=1,nin)
 1    format('lat   ',a,' gmode=',i1,/
     +      ,'---  ',6(i8))
c-----        for each latitude
      do 200 jlat=lat1,lat2,-2
         klat = jlat
 100     continue
         do 101 j=1,nin
            temp(j) = fwbpgv(uin(j),float(klat),gmode)
 101     continue
         write(*,2) klat,(temp(j),j=1,nin)
 2       format(i3,2x,6(f8.3))
         if(mod(klat,10).eq.2) print *,' '
         if(klat.eq.46) then
            klat = 45
            goto 100
         endif
 200  continue
c-----mm    : write header
      write(*,1) 'millimeters:',gmode,(nint(umm(j)),j=1,nmm)
c-----        for each latitude
      do 400 jlat=lat1,lat2,-2
         klat = jlat
 300     continue
         do 301 j=1,nmm
            temp(j) = fwbpgv(umm(j),float(klat),gmode)      
 301     continue
         write(*,3) klat,(temp(j),j=1,nmm)
 3       format(i3,2x,6(f8.2))
         if(mod(klat,10).eq.2) print *,' '
         if(klat.eq.46) then
            klat = 45
            goto 300
         endif
 400  continue
c-----mb    : write header
      write(*,1) 'millibars:  ',gmode,(nint(umb(j)),j=1,nmb)
c-----        for each latitude
      do 600 jlat=lat1,lat2,-2
         klat = jlat
 500     continue
         do 501 j=1,nmb
            temp(j) = fwbpgv(umb(j),float(klat),gmode)      
 501     continue
         write(*,3) klat,(temp(j),j=1,nmb)
         if(mod(klat,10).eq.2) print *,' '
         if(klat.eq.46) then
            klat = 45
            goto 500
         endif
 600  continue
      end
C-----------------------------------------------------------------------3456789
      subroutine tpbpg2(verbose)
c-----Test routine for {fwbpgv}.  Intercompare the corrections (c) produced by
c     {fwbpgv} in gmodes 1-3 over the range of latitudes (rlat) and barometric
c     pressure (bp) (870.0:1074.6).  If verbose is true, smaller ranges are
c     used and detailed results printed.  Note: Do not call more than once lest
c     proper initialization not occur.
c-----sdw, 4 Dec 1999.
c-----sdw, 2 Feb 2000: add initial write statement etc., rename (fr {tpbpgv}).
      logical verbose
      integer j,jlat,jbp,jlat1,jlat2,jbp1,jbp2,ir12,ir13,ir23
      integer ncomb,nx12,nx13,nx23,mx12,mx13,mx23
      real c(3),cn(3),cx(3),bpc(3)
      real cx12,cx13,cx23,px12,px13,px23
      real bx12,bx13,bx23,qx12,qx13,qx23
      data jlat1,jlat2/-900,900/,jbp1,jbp2/8700,10746/
      data ncomb/0/,nx12,nx13,nx23/3*0/,mx12,mx13,mx23/3*0/
      data ir12,ir13,ir23/3*0/
      data cn/3*999999.9/
      data cx,cx12,cx13,cx23/6*-999999.9/
      write(*,'(/,"tpbpg2 output:")')
      if(verbose) then
         jlat1=-900
         jlat2=-899
         jbp1 =8890
         jbp2 =8896
      endif
c-----for each latitude and pressure value
      do 500 jlat=jlat1,jlat2,1
      do 500 jbp = jbp1, jbp2,1
         rlat = float(jlat)/10.
         bp   = float(jbp )/10.
c-----count the number of rlat/bp combinations
         ncomb = ncomb + 1
         do 200 j=1,3
            c(j) = fwbpgv(bp,rlat,j)
            if(c(j).lt.cn(j)) cn(j) = c(j)
            if(c(j).gt.cx(j)) cx(j) = c(j)
 200     continue
c-----absolute c differences
         ax12 = abs(c(1)-c(2))
         ax13 = abs(c(1)-c(3))
         ax23 = abs(c(2)-c(3))
c-----carry max difference
         if(ax12.gt.cx12) cx12 = ax12
         if(ax13.gt.cx13) cx13 = ax13
         if(ax23.gt.cx23) cx23 = ax23
c-----count number of reversed differences
         if(c(1)-c(2).lt.0.) ir12 = ir12 + 1
         if(c(1)-c(3).lt.0.) ir13 = ir13 + 1
         if(c(2)-c(3).lt.0.) ir23 = ir23 + 1
c-----count differences >= ~0.05 (allowing for finite precision)
         if((ax12+0.001).ge.0.05) nx12 = nx12 + 1
         if((ax13+0.001).ge.0.05) nx13 = nx13 + 1
         if((ax23+0.001).ge.0.05) nx23 = nx23 + 1
c-----add c values to bp, round to tenths
         do 300 j=1,3
            bpc(j) = nint((bp + c(j))*10.)/10.
 300     continue
c-----absolute corrected bp differences
         bx12 = abs(bpc(1)-bpc(2))
         bx13 = abs(bpc(1)-bpc(3))
         bx23 = abs(bpc(2)-bpc(3))
c-----count differences >= ~0.1 (allowing for finite precision)
         if((bx12+0.01).ge.0.1) mx12 = mx12 + 1
         if((bx13+0.01).ge.0.1) mx13 = mx13 + 1
         if((bx23+0.01).ge.0.1) mx23 = mx23 + 1
c-----print out detailed results
         if(verbose) write(*,400)
     +   rlat,bp,c,ax12,ax13,ax23,cx12,cx13,cx23,ir12,ir13,ir23
     +  ,nx12,nx13,nx23,bpc,bx12,bx13,bx23,mx12,mx13,mx23
 400  format('verbose=true:          fwbpgv gmode ='/
     +      ,' rlat      bp   1  (g45,g0)   2   (g0,g0)   3 (g45,g45)',/    
     +      ,'----- -------   -----------   -----------   -----------',/
     +      ,f5.1,1x,f7.1,3x,f11.5,3x,f11.5,3x,f11.5,/
     +      ,'                correction differences =',/
     +      ,'                    1-2           1-3           2-3',/
     +      ,'                -----------   -----------   -----------',/
     +      ,'abs corr diff:  ',f11.5,3x,f11.5,3x,f11.5,/
     +      ,'max corr diff:  ',f11.5,3x,f11.5,3x,f11.5,/
     +      ,'no. neg. diff:  ',i11,3x,i11,3x,i11,/
     +      ,'no. >= ~0.05:   ',i11,3x,i11,3x,i11,/
     +      ,'                bp+correction =',/
     +      ,'                1  (g45,g0)   2   (g0,g0)   3 (g45,g45)',/    
     +      ,'                -----------   -----------   -----------',/
     +      ,'round->10ths:   ',f11.1,3x,f11.1,3x,f11.1,/
     +      ,'                bp+corr differences =',/
     +      ,'                    1-2           1-3           2-3',/
     +      ,'                -----------   -----------   -----------',/
     +      ,'abs diff:       ',f11.1,3x,f11.1,3x,f11.1,/
     +      ,'no. >= ~0.1:    ',i11,3x,i11,3x,i11,/,75('='))
 500  continue
c-----calculate percentages
      px12 = float(nx12)/float(ncomb)*100.
      px13 = float(nx13)/float(ncomb)*100.
      px23 = float(nx23)/float(ncomb)*100.
      qx12 = float(mx12)/float(ncomb)*100.
      qx13 = float(mx13)/float(ncomb)*100.
      qx23 = float(mx23)/float(ncomb)*100.
c-----write results
      write(*,601)
      write(*,602) cn,cx,cx12,cx13,cx23,ir12,ir13,ir23,nx12,nx13,nx23
     +            ,px12,px13,px23,mx12,mx13,mx23,qx12,qx13,qx23,ncomb
 601  format('c verbose=false:        fwbpgv gmode =')
 602  format('c             1  (g45,g0)   2   (g0,g0)   3 (g45,g45)',/    
     +      ,'c             -----------   -----------   -----------',/
     +      ,'c minimum:    ',f11.5,3x,f11.5,3x,f11.5,/   
     +      ,'c maximum:    ',f11.5,3x,f11.5,3x,f11.5,/,'c',/
     +      ,'c             absolute value of corr diff =',/
     +      ,'c                 1-2           1-3           2-3',/
     +      ,'c             -----------   -----------   -----------',/
     +      ,'c maximum:    ',f11.5,3x,f11.5,3x,f11.5,/
     +      ,'c no. neg dif:',i11,3x,i11,3x,i11,/,'c',/
     +      ,'c             combos w/ abs corr diff >= ~0.05 =',/
     +      ,'c                 1-2           1-3           2-3',/
     +      ,'c             -----------   -----------   -----------',/
     +      ,'c number:     ',i11,3x,i11,3x,i11,/
     +      ,'c percentage: ',f11.2,3x,f11.2,3x,f11.2,/,'c',/
     +      ,'c             combos w/ abs bp+corr diff >= ~0.1 =',/
     +      ,'c                 1-2           1-3           2-3',/
     +      ,'c             -----------   -----------   -----------',/
     +      ,'c number:     ',i11,3x,i11,3x,i11,/
     +      ,'c percentage: ',f11.2,3x,f11.2,3x,f11.2,/,'c',/
     +      ,'c total number of combinations:',i15)
      end
C-----------------------------------------------------------------------3456789
      subroutine tpbpmb
c-----Test: recalculate and print results for examples of mercurial barometer
c     corrections given in WBAN, pp. 3-5 and 3-6; and in UKMO (1969) p. 7.
c     The corrections are performed by two methods: (a) per aforementioned
c     WBAN examples; (b) and according to apparent UKMO method.  Corrections
c     are repeated using gmode values 1-2 of {fwbpgv}.  See WBAN chapter 5
c     for theoretical and approximate correction formulae, and its App.1.4.1
c     for International Barometer Conventions taking effect 1 January 1955.
c-----sdw, 25 Feb 2000.
c     References:
c     UKMO (UK Met. Office), 1969: Marine Observer's Handbook (9th ed.).
c           HMSO, London, 152 pp.
c     US Weather Bureau, Air Weather Service, and Naval Weather Service,
c           1963: Federal Meteorological Handbook No. 8--Manual of Barometry
c           (WBAN), Volume 1 (1st ed.).  US GPO, Washington, DC.
      implicit integer(a-e,g-z)
      real flat(4),bp(4),bt(4),ci(4),cg,ct,cs(4),b0i,b0t,b0g,bp1,bp2
     +    ,bp3(4)
      integer u(4),gmode
      character*11 eg(4)
c-----text describing example
      data eg/'WBAN p.3-5','WBAN p.3-5','UKMO pp.6-7','UKMO p.7'/
c-----latitude (flat)
      data flat/55.37, 4.03, 51.00, 27.00/
c-----barometric pressure (bp)
      data  bp/29.805,1006.8,30.240,1017.3/
c-----barometric pressure: corrected published values
      data bp3/1009.0,1001.2,30.215,1013.5/
c-----attached themometer (bt)
      data  bt/  55.0,  30.0,  60.0, 298.0/
c-----instrumental (index) correction (ci)
      data  ci/-0.015,+0.300,+0.005,+0.300/
c-----reduction of pressure to sea level (cs)
c     Note: calculation not available in {lmrlib}; enter from references
      data  cs/+0.048,+1.900,+0.039,+1.800/ 
c-----units (u): 1=mb/deg C, 2=inches/deg F
      data u/2,1,2,1/
      write(*,'(/,"tpbpmb output:")')
c-----2nd UKMO example: transform absolute (K) bt to degrees Celsius
c     Note: {fwbptg} doesn't handle Kelvin; try converting to Celsius
      bt(4)  = fxtktc(bt(4))
c-----1st UKMO example: transform inches bp3 (corrected value) to mb
      bp3(3) = fxeimb(bp3(3))
      write(*,1)
 1    format('note: original units (u): 1=mb/deg C, 2=in/deg F (<1955)')
c-----for each gmode
      do 900 gmode = 1,2
      write(*,2) gmode
 2    format(/,'============= fwbpgv gmode = ',i1,' =============')
c-----method (a): this is the method of WBAN pp. 3-5 and 3-6, where all
c     correction calculations use the observed pressure value
      write(*,3) 'method (a): calculations use observed press (WBAN)'
 3    format(a,/
     +,'   lat  barometric  attach   index    temp    grav'
     +,'  SLevel  corrected  corrected  corrected difference',/
     +,'u         pressure  therm.    corr    corr    corr'
     +,'    corr   (orig u)        hPa  published  hPa - pub',/
     +,'- -----  ---------  ------  ------  ------  ------'
     +,'  ------  ---------  ---------  ---------  ---------')
      do 400 j=1,4
c-----calculate gravity and temperature corrections
         cg = fwbpgv(bp(j),flat(j),gmode) 
         ct = fwbptg(bp(j),bt(j),u(j))
c-----correct to bp1 in original units
         bp1 = bp(j) + ci(j) + cg + ct + cs(j)
c-----if applicable, convert inches to mb
         if(u(j).eq.2) then
            bp2 = fxeimb(bp1)
         else
            bp2 = bp1
         endif
         write(*,4) u(j),flat(j),bp(j),bt(j),ci(j),ct,cg,cs(j),bp1,bp2
     +             ,bp3(j),(bp2-bp3(j)),eg(j)
 4       format(i1,1x,f5.2,2x,f9.3,2x,f6.1,2x,f6.3,2x,f6.3,2x,f6.3,2x
     +         ,f6.3,4(2x,f9.3),1x,a)
 400  continue
c-----method (b): this appears to be the method of UKMO (1969)
      write(*,3) 'method (b): incremental calculations (UKMO)'
      do 800 j=1,4
c     add index   correction yielding b0i, which is used as input to the
c     temperature correction yielding b0t, which is used as input to the
c     the gravity correction yielding b0g, which is used as input to the
c     height      correction yielding bp1 (else as method a)
         b0i = bp(j) + ci(j)
         ct  = fwbptg(b0i,bt(j),u(j))
         b0t = b0i + ct
         cg  = fwbpgv(b0t,flat(j),gmode) 
         b0g = b0t + cg
         bp1 = b0g + cs(j)
c-----if applicable, convert inches to mb
         if(u(j).eq.2) then
            bp2 = fxeimb(bp1)
         else
            bp2 = bp1
         endif
         write(*,4) u(j),flat(j),bp(j),bt(j),ci(j),ct,cg,cs(j),bp1,bp2
     +             ,bp3(j),(bp2-bp3(j)),eg(j)
 800  continue
 900  continue
      end
C=======================================================================3456789
c-----cloud conversions-------------------------------------------------3456789
C=======================================================================3456789
      integer function ixt0ok(t0)
c-----Convert "proportion of sky clear" in tenths (t0), to oktas (eighths
c     of sky covered; WMO code 2700).  The t0 code, specified in Maury
c     (1854), was documented for use, e.g., for US Marine Meteorological
c     Journals (1878-1893).  The dates of transition to instead reporting
c     "proportion of sky covered" (t1, as handled by {ixt1ok}) may have
c     varied nationally.  Following shows the mappings of t0/t1 to oktas
c     as provided by these routines ({tpt0t1} output):
c     10ths clear (t0)   oktas   10ths cover (t1)   oktas
c     ----------------   -----   ----------------   -----
c                   10       0                  0       0
c                    9       1                  1       1
c                    8       2                  2       2
c                    7       2                  3       2
c                    6       3                  4       3
c                    5       4                  5       4
c                    4       5                  6       5
c                    3       6                  7       6
c                    2       6                  8       6
c                    1       7                  9       7
c                    0       8                 10       8
c-----sdw, 26 Jan 2000.
c     Reference:
c     Maury, M.F., 1854: Maritime Conference held at Brussels for devising
c           a uniform system of meteorological observations at sea, August
c           and September, 1853.  Explanations and Sailing Directions to
c           Accompany the Wind and Current Charts, 6th Ed., Washington, DC,
c           pp. 54-88.
      integer t0,t1
c-----check validity of t0
      if(t0.lt.0.or.t0.gt.10) then
         print *,'ixt0ok error. illegal t0=',t0
         stop
      endif  
c-----convert tenths of "sky clear" (t0) to tenths of "sky covered" (t1)
c     (Note: assumption: no known basis in documentation)
      t1 = 10 - t0 
c-----convert tenths of "sky covered" to oktas
      ixt0ok = ixt1ok(t1)
      return
      end
C=======================================================================3456789
      integer function ixt1ok(t1)
c-----Convert tenths (of sky covered) (t1), to oktas (eighths of sky
c     covered; WMO code 2700).  This implements the mapping of tenths
c     to oktas shown below (left-hand columns) from NCDC (1968), section
c     4.5, scale 7.  In contrast, the right-hand columns show a reverse
c     mapping of "code no." (referring to oktas in the synoptic code)
c     back to tenths from Riehl (1947) (the justifications for the two
c     approaches are not known):
c           oktas  <-  tenths     |    code no.  ->  tenths
c           -----      -------    |    --------      -------
c             0         0         |        0           0
c             1         1         |        1           0
c             2         2 or 3    |        2           1
c             3         4         |        3           2.5
c             4         5         |        4           5
c             5         6         |        5           7.5
c             6         7 or 8    |        6           9
c             7         9         |        7          10 
c             8        10         |        8          10
c             9        obscured
c     Input t1 values must be limited to 0-10; "obscured" is not handled.
c-----sdw, 26 Jan 2000.
c-----sdw, 24 Jul 2000: correction of comment.
c     References:
c     NCDC (National Climatic Data Center), 1968: TDF-11 Reference Manual.
c           NCDC, Asheville, NC.
c     Riehl, 1947: Diurnal variation of cloudiness over the subtropical
c           Atlantic Ocean.  Bull. Amer. Meteor. Soc., 28, 37-40.
      integer t1,ok(0:10)
      data ok/0,1,2,2,3,4,5,6,6,7,8/
c-----check validity of t1
      if(t1.lt.0.or.t1.gt.10) then
         print *,'ixt1ok error. illegal t1=',t1
         stop
      endif  
c-----convert from tenths to oktas
      ixt1ok = ok(t1)
      return
      end
C-----------------------------------------------------------------------3456789
      subroutine tpt0t1
c-----test: print {ixt0ok,ixt1ok} results.  
c-----sdw, 26 Jan 2000.
      implicit integer(a-e,g-z)
      write(*,'(/,"tpt0t1 output:")')
      write(*,1)
 1    format('10ths clear (t0)   oktas   10ths cover (t1)   oktas',/
     +      ,'----------------   -----   ----------------   -----')
      do 100 t1=0,10
         t0  = 10 - t1
         ok0 = ixt0ok(t0)
         ok1 = ixt1ok(t1)
         write(*,2) t0,ok0,t1,ok1
 2       format(14x,i2,7x,i1,17x,i2,7x,i1)
 100  continue
      return
      end
C=======================================================================3456789
c-----temperature conversions-------------------------------------------3456789
C=======================================================================3456789
      real function fxtftc(tf)
c-----Convert temperature in degrees Fahrenheit (tc) to degrees Celsius.
c-----sdw, 26 Jan 2000 (replaces 1 Jul 1998 version from colib5s.01J).
c     Reference:
c     List, R.J., 1966: Smithsonian Meteorological Tables.
c           Smithsonian Institution, Washington, DC, 527 pp.
      real tf
c-----equation from List (1966), Table 2 (p. 17).
      fxtftc = (5.0/9.0) * (tf - 32.0)
      return
      end
C-----------------------------------------------------------------------3456789
      real function fxtctf(tc)
c-----Convert temperature in degrees Celsius (tc) to degrees Fahrenheit.
c-----sdw, 26 Jan 2000 (replaces 1 Jul 1998 version from colib5s.01J).
c     Reference:
c     List, R.J., 1966: Smithsonian Meteorological Tables.
c           Smithsonian Institution, Washington, DC, 527 pp.
      real tc
c-----equation from List (1966), Table 2 (p. 17).
      fxtctf = ((9.0/5.0) * tc) + 32.0
      return
      end
C-----------------------------------------------------------------------3456789
      real function fxtktc(tk)
c-----Convert temperature in Kelvins (tk) to degrees Celsius.
c-----Adapted from colib5s.01J function {cvtkc} (1984); sdw, 1 Jul 1998.
      real tk
      if(tk.lt.0.0) then
         print *,'fxtktc error. negative input tk=',tk
         stop
      endif
      fxtktc = tk - 273.15
      return
      end
C-----------------------------------------------------------------------3456789
      real function fxtctk(tc)
c-----Convert temperature in degrees Celsius (tc) to Kelvins.
c-----Adapted from colib5s.01J function {cvtck} (1984); sdw, 1 Jul 1998.
      real tc
      fxtctk = tc + 273.15
      if(fxtctk.lt.0.0) then
         print *,'fxtctk error. negative output=',fxtctk
         stop
      endif
      return
      end
C-----------------------------------------------------------------------3456789
      real function fxtrtc(tr)
c-----Convert temperature in degrees Reaumur (tc) to degrees Celsius.
c-----sdw, 26 Jan 2000 (replaces 1 Jul 1998 version from colib5s.01J).
c     Reference:
c     List, R.J., 1966: Smithsonian Meteorological Tables.
c           Smithsonian Institution, Washington, DC, 527 pp.
      real tr
c-----equation from List (1966), Table 2 (p. 17).
      fxtrtc = (5.0/4.0) * tr
      return
      end
C-----------------------------------------------------------------------3456789
      real function fxtctr(tc)
c-----Convert temperature in degrees Celsius (tc) to degrees Reaumur.
c-----sdw, 26 Jan 2000 (replaces 1 Jul 1998 version from colib5s.01J).
c     Reference:
c     List, R.J., 1966: Smithsonian Meteorological Tables.
c           Smithsonian Institution, Washington, DC, 527 pp.
      real tc
c-----equation from List (1966), Table 2 (p. 17).
      fxtctr = (4.0/5.0) * tc
      return
      end
C-----------------------------------------------------------------------3456789
      subroutine tptcrf
c-----test: print {fxtftc,fxtctf,fxtrtc,fxtctr} results.  Using as input
c     the Celsius column, calculate the Reaumur and Fahrenheit columns, of
c     Table 2 of Lamb (1986).  Finite precision (32-bit Sun f77) apparently
c     introduces some 0.1 degree differences shown below (units and tenths
c     positions of original Table 2 values listed to right in parentheses).
c     Also tabulated are conversions from Reaumur and Fahrenheit back to
c     Celsius:  
c             C   C -> R   C -> F    R -> C   F -> C
c        ------   ------   ------    ------   ------
c         100.0     80.0    212.0     100.0    100.0
c          62.5     50.0    144.5      62.5     62.5
c          53.5     42.8    128.3      53.5     53.5
c          51.9     41.5    125.4      51.9     51.9
c          38.9     31.1    102.0      38.9     38.9
c          35.0     28.0     95.0      35.0     35.0
c          32.2     25.8     90.0      32.2     32.2
c          31.1     24.9     88.0      31.1     31.1
c          23.3     18.6     73.9(4.0) 23.3     23.3
c          15.6     12.5     60.1(0.0) 15.6     15.6
c          14.4     11.5     57.9(8.0) 14.4     14.4
c          11.7      9.4(9.3)53.1(3.0) 11.7     11.7
c          10.0      8.0     50.0      10.0     10.0
c           3.9      3.1     39.0       3.9      3.9
c           1.1      0.9     34.0       1.1      1.1
c           0.0      0.0     32.0       0.0      0.0
c          -3.3     -2.6     26.1(6.0) -3.3     -3.3
c          -3.9     -3.1     25.0      -3.9     -3.9
c          -7.8     -6.2(6.3)18.0      -7.8     -7.8
c         -11.7     -9.4     10.9(1.0)-11.7    -11.7
c         -12.2     -9.8     10.0     -12.2    -12.2
c         -18.7    -15.0     -1.7(1.6)-18.7    -18.7
c         -19.7    -15.8     -3.5(3.4)-19.7    -19.7
c         -21.1    -16.9     -6.0     -21.1    -21.1
c         -23.8    -19.0    -10.8     -23.8    -23.8
c-----sdw, 26 Jan 2000.
c     Reference:
c     Lamb, H.H., 1986: Ancient units used by the pioneers of meteorological
c           instruments.  Weather, 41, 230-233.
      implicit integer(a-e,g-z)
      real tc(25),tf(25),tr(25)
c-----input Celsius values
      data tc/100.0, 62.5, 53.5, 51.9, 38.9, 35.0, 32.2, 31.1, 23.3
     +       , 15.6, 14.4, 11.7, 10.0,  3.9,  1.1,  0.0, -3.3, -3.9
     +       , -7.8,-11.7,-12.2,-18.7,-19.7,-21.1,-23.8/
      write(*,'(/,"tptcrf output:")')
      write(*,1)
 1    format('Celsius, Reaumur, and Fahrenheit conversions:',/
     +      ,'        C   C -> R   C -> F    R -> C   F -> C',/
     +      ,'   ------   ------   ------    ------   ------')
c-----print Reaumur and Fahrenheit, plus reverse calculations of Celsius
      do 100 j=1,25
         tr(j) = fxtctr(tc(j))
         tf(j) = fxtctf(tc(j))
         write(*,2) tc(j),tr(j),tf(j),fxtrtc(tr(j)),fxtftc(tf(j))
 2       format(3(3x,f6.1),1x,2(3x,f6.1))
 100  continue
      return
      end
C=======================================================================3456789
c-----wind conversions--------------------------------------------------3456789
C=======================================================================3456789
      real function fxuvdd(u,v)
c-----Convert wind vector eastward and northward components (u,v) to
c     direction (from) in degrees (clockwise from 0 degrees North).
c-----Adapted from colib5s.01J function {dduv} (1984); sdw, 1 Jul 1998.
      real u,v,a
      if (u.eq.0.0.and.v.eq.0.0) then
         a = 0.0
      else
         a = atan2(v,u)*(180.0/3.14159265358979323846264338327950288)
      endif
      fxuvdd = 270.0 - a
      if(fxuvdd.ge.360.0) fxuvdd = fxuvdd - 360.0
      return 
      end
C-----------------------------------------------------------------------3456789
      real function fxuvvv(u,v)
c-----Convert wind vector eastward and northward components (u,v) to
c     velocity.
c-----Adapted from colib5s.01J function {vvuv} (1984); sdw, 1 Jul 1998.
      real u,v
      fxuvvv = sqrt(u**2 + v**2)
      return
      end
C-----------------------------------------------------------------------3456789
      subroutine rxdvuv(dd,vv,u,v)
c-----Convert wind direction (dd; clockwise from 0 degrees North) and
c     velocity (vv) to vector eastward and northward components (u,v).
c-----Adapted from colib5s.01J function {uvddvv} (1984); sdw, 1 Jul 1998.
      real dd,vv,u,v,ang
      ang = dd * (3.14159265358979323846264338327950288/180.0)
      u = -vv * sin(ang)
      v = -vv * cos(ang)
      return
      end
C-----------------------------------------------------------------------3456789
      real function fxktms(kt)
c-----Convert from knots (kt; with respect to the international nautical
c     mile) to meters per second (see {tpktms} for details).
c-----Adapted from colib5s.01J function {cvskm} (1984); sdw, 26 Jun 1998.
      real kt
      fxktms = kt * 0.51444444444444444444
      return
      end
C-----------------------------------------------------------------------3456789
      real function fxmskt(ms)
c-----Convert from meters per second (ms) to knots (with respect to the
c     international nautical mile) (see {tpktms} for details).
c-----Adapted from colib5s.01J function {cvsmk} (1984); sdw, 26 Jun 1998.
      real ms
      fxmskt = ms *  1.9438444924406047516
      return
      end
C-----------------------------------------------------------------------3456789
      real function fxk0ms(k0)
c-----Convert from knots (k0; with respect to the U.S. nautical mile) to
c     meters per second (see {tpktms} for details).
c-----sdw, 26 Jun 1998.
      real k0
      fxk0ms = k0 * 0.51479111111111111111
      return
      end
C-----------------------------------------------------------------------3456789
      real function fxmsk0(ms)
c-----Convert from meters per second (ms) to knots (with respect to the
c     U.S. nautical mile) (see {tpktms} for details).
c-----sdw, 26 Jun 1998.
      real ms
      fxmsk0 = ms *  1.9425354836481679732
      return
      end
C-----------------------------------------------------------------------3456789
      real function fxk1ms(k1)
c-----Convert from knots (k1; with respect to the Admiralty nautical mile)
c     to meters per second (see {tpktms} for details).
c-----sdw, 26 Jun 1998.
      real k1
      fxk1ms = k1 * 0.51477333333333333333
      return
      end
C-----------------------------------------------------------------------3456789
      real function fxmsk1(ms)
c-----Convert from meters per second (ms) to knots (with respect to the
c     Admiralty nautical mile) (see {tpktms} for details).
c-----sdw, 26 Jun 1998.
      real ms
      fxmsk1 = ms *  1.9426025694156651471
      return
      end
C-----------------------------------------------------------------------3456789
      integer function ixbfkt(bf)

! While a simple transformation, this is an important function in {lmrlib}.
! Note however the interesting question about midpoint value 18 below.

c-----Convert from Beaufort force 0-12 (bf) to "old" (WMO code 1100)
c     midpoint in knots.  From NCDC (1968), conversion scale 5 (sec.
c     4.4).  Note: Midpoint value 18 looks questionable, but appeared
c     originally in UKMO (1948).
c     References:
c     NCDC (National Climatic Data Center), 1968: TDF-11 Reference Manual.
c           NCDC, Asheville, NC.
c     UKMO (UK Met. Office), 1948: International Meteorological Code
c           Adopted by the International Meteorological Organisation,
c           Washington, 1947 (Decode for the Use of Shipping, effective
c           from 1st January, 1949).  Air Ministry, Meteorological Office,
c           HM Stationary Office, London, 39 pp.
c-----sdw, 16 Jun 1998.
c-----sdw, 26 Jan 2000: complete UKMO reference.
c-----sdw, 17 Aug 2000: routine name in correct error message.
      integer bf,kt(0:12)
      data kt/0,2,5,9,13,18,24,30,37,44,52,60,68/
      if(bf.lt.0.or.bf.gt.12) then
         print *,'ixbfkt error.  bf=',bf
         stop
      endif
      ixbfkt = kt(bf)
      return
      end
C-----------------------------------------------------------------------3456789
      real function fxbfms(bf)
c-----Convert from Beaufort force 0-12 (bf) to "old" (WMO code 1100)
c     midpoint in meters per second.  From Slutz et al. (1985) supp.
c     K, Table K5-5 (p. K29).  See {ixbfkt} for additional background.
c     Reference:
c     Slutz, R.J., S.J. Lubker, J.D. Hiscox, S.D. Woodruff, R.L. Jenne,
c           D.H. Joseph, P.M. Steurer, and J.D. Elms, 1985: Comprehensive
c           Ocean-Atmosphere Data Data Set; Release 1.  NOAA
c           Environmental Research Laboratories, Climate Research
c           Program, Boulder, Colo., 268 pp. (NTIS PB86-105723).
c-----sdw, 16 Jun 1998
      real ms(0:12)
      integer bf
      data ms/0.,1.,2.6,4.6,6.7,9.3,12.3,15.4,19.,22.6,26.8,30.9,35./
      if(bf.lt.0.or.bf.gt.12) then
         print *,'fxbfms error.  bf=',bf
         stop
      endif
      fxbfms = ms(bf)
      return
      end
C-----------------------------------------------------------------------3456789
      integer function ix32dd(c32,dc,imiss)
c-----Convert 4-character 32-point wind direction abbreviation c32 into
c     degrees, or return imiss if unrecognized; also return numeric code
c     1-32 (or imiss) in dc (see {ixdcdd} for background).  Recognized
c     abbreviations are in cwd, with these characteristics: left-justified,
c     upper-case, with trailing blank fill, and where "X" stands for "by".
c     NOTE: No constraint is placed on imiss (it could overlap with data).
c-----sdw, 17 Aug 2000.
      implicit integer(a-e,g-z)
      character*4 c32,cwd(32)
      data cwd/'NXE ','NNE ','NEXN','NE  ','NEXE','ENE ','EXN ','E   ',
     +         'EXS ','ESE ','SEXE','SE  ','SEXS','SSE ','SXE ','S   ',
     +         'SXW ','SSW ','SWXS','SW  ','SWXW','WSW ','WXS ','W   ',
     +         'WXN ','WNW ','NWXW','NW  ','NWXN','NNW ','NXW ','N   '/
      ix32dd = imiss
      do 500 j=1,32
         if(c32.eq.cwd(j)) then
            ix32dd = ixdcdd(j,imiss)
            dc     = j
            return
         endif
 500  continue
      return
      end
C-----------------------------------------------------------------------3456789
      integer function ixdcdd(dc,imiss)
c-----Convert 32-point wind direction numeric code dc into degrees, or
c     return imiss if dc is out of range 1-32.  Release 1, Table F2-1
c     defines the mapping of code dc to degrees in dwd.
c     NOTE: No constraint is placed on imiss (it could overlap with data).
c-----sdw, 17 Aug 2000.
      implicit integer(a-e,g-z)
      dimension dwd(32)
      data dwd/    11,    23,    34,    45,    56,    68,    79,    90,
     +            101,   113,   124,   135,   146,   158,   169,   180,
     +            191,   203,   214,   225,   236,   248,   259,   270,
     +            281,   293,   304,   315,   326,   338,   349,   360/
      if(dc.ge.1.and.dc.le.32) then
         ixdcdd = dwd(dc)
      else
         ixdcdd = imiss
      endif   
      return
      end
C-----------------------------------------------------------------------3456789
      subroutine txdvuv
c-----Test {fxuvdd,fxuvvv,rxdvuv}.
c-----sdw, 1 Jul 1998.
      real u,v,dd,vv
      do 200 iu=-50,50
      do 200 iv=-50,50
         u = float(iu)
         v = float(iv)
         dd = fxuvdd(u,v)
         vv = fxuvvv(u,v)
         call rxdvuv(dd,vv,u,v)
         if(nint(u).ne.iu.or.nint(v).ne.iv) then
            print *,'txdvuv error. iu=',iu,', iv=',iv
     +      ,', dd = fxuvdd(u,v) =',dd
     +      ,', vv = fxuvvv(u,v) =',vv
     +      ,', u=',u,', v=',v
            stop
         endif
 200  continue
      end
C-----------------------------------------------------------------------3456789
      subroutine tpktms
c-----Calculate (in double precision) and print factors for conversions
c     from knots, (with reference to the international, U.S., and Admiralty
c     nautical miles) to meters per second, and for the reverse conversions.
c     One knot is defined as 1 nautical mile (nm) per hour, and there are
c     3600 seconds per hour (sph), and thus: nm/sph meters per second.  We
c     calculate the reverse factor as the inverse: 1/(nm/sph).  (Note: 20
c     significant digits were retained in transferring factors to routines.)
c-----sdw, 1 Jul 1998.
c-----sdw, 26 Jan 2000: add initial write statement, plus fix comment typos.
      double precision sph,fkt,fms,fft
      sph = 3600.d0
c-----Calculate and print factors for {fxktms/fxmskt}, where the international
c     nautical mile (nm) is defined as 1852 meters (WMO, 1966, Table 1-1).
c     References:
c     List, R.J., 1966: Smithsonian Meteorological Tables.
c           Smithsonian Institution, Washington, DC, 527 pp.
c     WMO (World Meteorological Organization), 1966: International
c           Meteorological Tables, WMO-No.188.TP.94.
      write(*,'(/,"tpktms output:")')
      print *,'tpktms.  results for international nautical mile:'
      fkt = 1852.d0/sph
      print *,'fkt=',fkt
      fms = 1.d0/fkt
      print *,'fms=',fms
      print *,'fxktms(1.)=',fxktms(1.)
      print *,'fxmskt(1.)=',fxmskt(1.)
      print *,'fxktms(45.)=',fxktms(45.)
c-----Calculate and print factors for {fxk0ms/fxmsk0}, where the U.S.
c     nautical mile is defined as 1853.248 meters = 6080.21 feet (List,
c     1966, p. 3).
      print *,'tpktms.  results for U.S. nautical mile:'
      fkt = 1853.248d0/sph
      print *,'fkt=',fkt
      fms = 1.d0/fkt
      print *,'fms=',fms
      print *,'fxk0ms(1.)=',fxk0ms(1.)
      print *,'fxmsk0(1.)=',fxmsk0(1.)
      print *,'fxk0ms(45.)=',fxk0ms(45.)
c-----Calculate and print factors for {fxk1ms/fxmsk1}, where the Admiralty
c     nautical mile is defined as 6080 feet (List, 1966, p. 3), and where
c     the factor fft for conversion from ft to m is from List, 1966, p. 2.
      print *,'tpktms.  results for Admiralty nautical mile:'
      fft = 0.3048d0
      fkt = (6080.d0*fft)/sph
      print *,'fkt=',fkt
      fms = 1.d0/fkt
      print *,'fms=',fms
      print *,'fxk1ms(1.)=',fxk1ms(1.)
      print *,'fxmsk1(1.)=',fxmsk1(1.)
      print *,'fxk1ms(45.)=',fxk1ms(45.)
      return
      end
C=======================================================================3456789
c-----time conversions--------------------------------------------------3456789
C=======================================================================3456789
      subroutine rxltut(ihr,idy,elon,uhr,udy)
c-----Convert local standard hour (ihr; in hundredths 0-2399) and "Julian"
c     day (i.e., any incrementable integer date) (idy) into coordinated
c     universal time (UTC) hour (uhr) and day (udy; decremented if the
c     dateline is crossed), using longitude (elon; in hundredths of degrees
c     0-35999, measured east of Greenwich).  Notes: a) Strict time zones,
c     including the International Date Line, are not employed.  b) In all
c     cases the western (eastern) boundary of each time zone is inclusive
c     (exclusive), i.e., 7.50W-7.49E, 7.50E-22.49E, ..., 172.50E-172.51W.
c-----sjw and sdw, 18 Jun 1998.
      integer ihr,idy,elon,uhr,udy
     +,wlon,dhr
      if(ihr.lt.0.or.ihr.gt.2399) then
         print *,'error rxltut. ihr=',ihr
         stop
      else if(elon.lt.0.or.elon.gt.35999) then
         print *,'error rxltut. elon=',elon
         stop
      endif
      wlon = 36000 - elon
      udy = idy
      dhr = (wlon + 749)/1500
      uhr = ihr + dhr*100
      if(uhr.ge.2400)then
         udy = udy + 1
         uhr = uhr - 2400
      endif
      if(wlon.ge.18000) udy = udy - 1
      return
      end
C-----------------------------------------------------------------------3456789
      integer function ixdtnd(iday,imonth,iyear)
c-----Convert from date (iday,imonth,iyear) to number of days since
c     1 Jan 1770.
c-----sjl, 17 Jun 1998.
c-----sjw and sdw, 27 Apr 1999: return ixdtnd=-1 if date is invalid.
c-----sdw, 26 Jan 2000: remove outdated variable ierr/comment.
c-----sjl, 17 Nov 2004: remove print statements.
      integer iday,imonth,iyear
     +,year,days(12)
      data days/31,28,31,30,31,30,31,31,30,31,30,31/
      logical leap
      leap(year) = mod(year,4).eq.0 .and. mod(year,100).ne.0
     + .or. mod(year,400).eq.0
c
      ixdtnd = -1
      if (iyear.lt.1770 .or. imonth.lt.1 .or. imonth.gt.12
     + .or. iday.lt.1 .or. iday.gt.days(imonth)
     + .and. (imonth.ne.2 .or. .not.leap(iyear) .or. iday.ne.29)) return
      ndays = 0
      do 190 year = 1770,iyear-1
        ndays = ndays + 365
        if (leap(year)) ndays = ndays + 1
  190 continue
      do 290 month = 1,imonth-1
        ndays = ndays + days(month)
        if (month.eq.2 .and. leap(iyear)) ndays = ndays + 1
  290 continue
      ndays = ndays + iday-1
      ixdtnd = ndays
      end
C-----------------------------------------------------------------------3456789
      subroutine rxnddt(ndays,iday,imonth,iyear)
c-----Convert from number of days (ndays) since 1 Jan 1770 to
c     date (iday,imonth,iyear).
c-----sjl, 17 Jun 1998.
c-----sjl, 17 Nov 2004: remove print statement and return
c     iday=-1, imonth=-1, and iyear=-1 if ndays is invalid.
      integer ndays,iday,imonth,iyear
     +,year,days(12)
      data days/31,28,31,30,31,30,31,31,30,31,30,31/
      logical leap
      leap(year) = mod(year,4).eq.0 .and. mod(year,100).ne.0
     + .or. mod(year,400).eq.0
c
      iday = -1
      imonth = -1
      iyear = -1
      if (ndays.lt.0) return
      mdays = ndays
      iyear = 1770
  100 continue
      n = 365
      if (leap(iyear)) n = n + 1
      if (mdays - n.ge.0) then
        mdays = mdays - n
        iyear = iyear + 1
        goto 100
      endif
      imonth = 1
  200 continue
      n = days(imonth)
      if (imonth.eq.2 .and. leap(iyear)) n = n + 1
      if (mdays - n.ge.0) then
        mdays = mdays - n
        imonth = imonth + 1
        goto 200
      endif
      iday = mdays + 1
      end
C-----------------------------------------------------------------------3456789
      subroutine txltut
c-----Test {rxltut}.
c-----sjl, 22 Jun 1998.
      integer ihr,idy,elon,uhr,udy
     +,hr,dy
      parameter(idy=366)
      do 190 ihr=0,2399
      do 190 elon=0,35999
        hr=ihr-(elon+750)/1500*100
        dy=idy
        if (elon.le.18000) then
          if (hr.lt.0) then
            hr=hr+2400
            dy=dy-1
          endif
        else
          hr=hr+2400
          if (hr.ge.2400) then
            hr=hr-2400
            dy=dy+1
          endif
        endif
        call rxltut(ihr,idy,elon,uhr,udy)
        if (uhr.ne.hr .or. udy.ne.dy) then
          print *,'error txltut.  ihr=',ihr,', idy=',idy,', elon=',elon
     +    ,', uhr=',uhr,', udy=',udy
          stop
        endif
  190 continue
      end
C-----------------------------------------------------------------------3456789
      subroutine txdtnd
c-----Test {ixdtnd,rxnddt}.
c-----sjl, 17 Jun 1998.
        i = 0
        do 190 iyear = 1770,2024
          do 180 imonth = 1,12
            do 170 iday = 1,31
              j = ixdtnd(iday,imonth,iyear)
              if (j.ne.i) then
                print *,'txdtnd error.  i=',i,', j=',j
                stop
              endif
              call rxnddt(j,jday,jmonth,jyear)
              if (jday.ne.iday .or. jmonth.ne.imonth
     +         .or. jyear.ne.iyear) then
                print *,'txdtnd error.  iday=',iday,', imonth=',imonth
     +          ,', iyear=',iyear,', jday=',jday,', jmonth=',jmonth
     +          ,', jyear=',jyear
                stop
              endif
              i = i+1
              if (iday.eq.28 .and. imonth.eq.2) then
                if (mod(iyear,400).eq.0) goto 170
                if (mod(iyear,100).eq.0) goto 180
                if (mod(iyear,4).eq.0) goto 170
                goto 180
              endif
              if (iday.eq.29 .and. imonth.eq.2) goto 180
              if (iday.eq.30 .and. (imonth.eq.4 .or. imonth.eq.6
     +         .or. imonth.eq.9 .or. imonth.eq.11)) goto 180
170         continue
180       continue
190     continue
        end
C=======================================================================3456789
c-----miscellaneous-----------------------------------------------------3456789
C=======================================================================3456789
      subroutine rpepsi
c-----Print calculated machine epsilon, i.e., the smallest power of 2,
c     2**no = ep, such that 1+2**no>1
c-----sdw, 15 Jul 1983; sdw minor revisions, 26 Jun 1998.
c-----sdw, 26 Jan 2000: add initial write statement. 
      write(*,'(/,"rpepsi output:")')
      EP=1./2.
      NO=-1 
 100  IF(1.+EP.LE.1.) GOTO 200
      EP=EP/2. 
      NO=NO-1 
      GOTO 100
 200  CONTINUE
      NO=NO+1 
      EP=EP*2.
      PRINT 300,NO,2.**NO,2.**NO,EP,EP
 300  FORMAT('rpepsi. 2.**',I4,7X,  '=',G30.14,'=',Z35,' (hex)',/ 
     +      ,'rpepsi. ep by division =',G30.14,'=',Z35,' (hex)') 
      return
      END 
C-----------------------------------------------------------------------3456789
      subroutine rpepsd
c-----Double precision version of {rpepsi} (see for background).
c-----sdw, 26 Jun 1998.
c-----sdw, 26 Jan 2000: add initial write statement. 
      double precision EP
      write(*,'(/,"rpepsd output:")')
      EP=1.d0/2.d0
      NO=-1 
 100  IF(1.d0+EP.LE.1.d0) GOTO 200
      EP=EP/2.d0 
      NO=NO-1 
      GOTO 100
 200  CONTINUE
      NO=NO+1 
      EP=EP*2.d0
      PRINT 300,NO,2.d0**NO,2.d0**NO,EP,EP
 300  FORMAT('rpepsd. 2.**',I4,7X,  '=',G30.14,'=',Z35,' (hex)',/ 
     +      ,'rpepsd. ep by division =',G30.14,'=',Z35,' (hex)') 
      return
      END 
C-----------------------------------------------------------------------3456789
      integer function imrnde(x)
c-----Round real x into integer imrnde such that a fractional part of x
c     of exactly 0.5 results in rounding to the nearest even integer.
c-----Adapted from colib5s.01J function {iround} (1984); sdw, 25 Jun 1998.
c-----sdw, 26 Jan 2000: fixed typo in comments.
      imrnde = x
      dx = x - imrnde
      jx = mod(imrnde,2)
      if ((dx.lt.-(.5)).or.(dx.eq.-(.5).and.jx.eq.-1))
     +imrnde = imrnde - 1
      if ((dx.gt.  .5) .or.(dx.eq.  .5 .and.jx.eq. 1))
     +imrnde = imrnde + 1
      return
      end
EOR
rm lmrlib.o
gfortran -c lmrlib.f

[Documentation and Software][Links to additional]


U.S. National Oceanic and Atmospheric Administration hosts the icoads website privacy disclaimer
Document maintained by icoads&#64;noaa.gov
Updated: Jun 5, 2014 20:10:45 UTC
http://icoads.noaa.gov/lmrlib.html