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. 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 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 and (electronic documents). C-----------------------------------------------------------------------3456789 c program libtst subroutine libtst 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) 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) 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