| Data and Documentation | Chronology and News | Program Status | Publications | Related Data and Resources | Contact Points |
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@noaa.gov
Updated: Jun 5, 2014 20:10:45 UTC
http://www.icoads.noaa.gov/lmrlib.html