-------
Integer :
Integer :
Integer :
Integer :
Integer :
Integer :
Integer :
Character
Integer :
Real :
Integer :
Integer :
Integer :
Real :
Real :
Integer :
Character
Integer :
Real :
Integer :
Integer :
Integer :
Character
Integer :
Logical :
Integer :
Integer :
Integer :
Character
Real :
Logical :
EHR
EDNR
GHR_v, GHR_u
yy, mm, dd, hh ! 2-digit year, month, day, hour
' 4-digit year
! Extraterrestrial Horizontal Radiation
! Extraterrestrial Direct Normal Radiation
! Global Horizontal Radiation
! Direct Normal Radiation
! Diffuse Horizontal Radiation
_ , DNR_s , DHR_s
Total_Sky_Cover, Opaque_Sky_Cover
Dry_Bulb_Temperature, Dew_Point_Temperature
Relative Humidity
Station_Pressure
Wind_Direction
Wind Speed ! in m/s
Horizontal Visibility
Ceiling_height ! in meters
Len=09) :: Present_weather
Precipitable water
baod ! Broadband aerosol optical depth
Snow_Depth
Days_Since_Last_Snowf all
Hourly Precipitation
(Len=01) : : Hourly Precipitation Flag
Fyear, nlines
ok, is daytime
doy, j yyyy, jmm, jdd, jdoy, jhh, jv
jd_today
nbase, iv
Len=200) :: tbuf, qline
ulO, Station elevation in meters
Xfound
Integer, Dimension(:,:), Pointer : :
Type(Val_and_Flag) :: d
Character(Len(d%s)) :: data source
!Logical
!Integer
Print now
jvO, jvl
'0'
! PWT — Present Weather Table
! Used to validate Present_weather's 9 flags
! Numbers with "*" denote additions to the SAMSON
! document. See .
Character(Len=*), Parameter :: PWT = '0123456789'
Integer :: i
nlines = 0
260
-------
Empty_File = .True.
! Open the input file
If (Len_trim(InFile) == 0) Then
! No file name.
Return
End If
If (Is_vll) Then
Call IORead(jin, InFile)
data_source = T_SAMSONvll
Else
Call Decompress_and_Open_File(Jin, InFile)
data_source = T_SAMSONvlO
End If
Fyear = WFlx%yyyy
Write (ULog, *) '## File: ', Trim(InFile)
Call ToTTy('Read SAMSON vlx '//Trim(InFile))
! Assign pointers
p_every_hour => WFlx%Every_Hour_Present
! * [Isr] 25 Mar 2002 8:50 am
! + We will copy the missing years from existing years.
! Say 1990 is missing. We will copy 1961 or 1964 (if
! leap year). This will allow make_rO to complete the
! run. We will then eliminate the 1990.hvf file and the
! MET file records associated with 1990.
I
! Fid Position Element Definition
I
! 001 - 001 Indicator Indicates the header record.
I
WBAN Number Station's Weather Bureau Army Navy number.
City
261
-------
! 031 - 032 State State where the station is located (abbreviated
! to two letters).
I
Time Zone Time zone is the number of hours by which the
local standard time lags or leads Universal
Time. For example, Mountain Standard Time is
designated -7 because it lags Universal Time by
7 hours.
Latitude Latitude of the station.
N = North of equator
Degrees
Minutes
Longitude Longitude of the station.
W = West, E = East
Degrees
Minutes
Elevation
Read (gline, 9130, lostat = ios) WBAN_Number, City, State, Time_Zone, &
Latitude Let, Latitude Degrees, Latitude Minutes, &
Longitude_Let, Longitude_Degrees, Longitude_Minutes, &
Elevation meters
If (WFlx%Check_Header_Info) Then
WFlx%Head%WBAN = WBAN_Number
WFlx%Head%Text = Trim(City) //','// Trim(State)
WFlx%Head%Lat = Coords(Latitude_Let, Latitude_Degrees, Latitude_Minutes)
WFlx%Head%Lon = Coords(Longitude Let, Longitude Degrees, Longitude Minutes)
WFlx%Head%Elev = Elevation_meters
WFlx%Head%TZ = Time_Zone
End If
262
-------
If (Present(Jout)) Then
Write (Jout, 9150) WBAN_Number, City, State, Time_Zone, &
Latitude_Let, Latitude_Degrees, Latitude_Minutes, &
Longitude_Let, Longitude_Degrees, Longitude_Minutes, &
Elevation meters
! Same format as
Format (
2x, al,
2x, al,
2x, 14)
! Add the output header
Write (Jout, 9170)
Format ('~YR MO DA HR I
8 9 10
17 18 19
End If
14
15
Table 2. Data Elements in the NSRDB Hourly Data Files
(For all except the first record of each file)
IP - position in input file
OP - field number in output file
! OP IP
Element
Local Standard Time
Year
Month
Day
Hour
[Isr] Thu 19 Oct 2000 11:58:19
This field appears to be
column 96 of the input file.
The manual provides no data
to support this assertion.
[Isr] Wed Nov 21 09:33:59 2001
Observation Indicator appears
only in SAMSON v 1.0 files.
Extraterrestrial
Horizontal Radiation
Values Definition
Year of observation
Month of observation
Day of month
Hour of day in local standard time
0 = Weather observation made.
9 = Weather observation not
made or missing.
If this field = 9 OR if field
13 (wind speed) = missing
(9999. or 99.0), then
fields 6, 7, 8, 10, 11, 17,
and 18 were all modeled and
not actually observed.
Amount of solar radiation in Wh/m2
received on a horizontal surface at
the top of the atmosphere during the
60 minutes preceding the hour indicated.
263
-------
Extraterrestrial Direct
Normal Radiation
Amount of solar radiation in Wh/m2
received on a surface normal to
the sun at the top of the atmosphere
during the 60 minutes preceding the
hour indicated.
Global Horizontal Radiation
Data Value 0-1415
Flag for Data Source A-H, ?
Flag for Data Uncertainty 0-9
Total amount of direct and diffuse
solar radiation in Wh/m2 received on a
horizontal surface during the 60
minutes preceding the hour indicated.
9999 = missing data.
Direct Normal Radiation
Data Value 0-141E
Flag for Data Source A-H, 1
Flag for Data Uncertainty 0-9
Diffuse Horizontal Radiation
Data Value 0-1415
Flag for Data Source A-H, ?
Flag for Data Uncertainty 0-9
Amount of solar radiation in Wh/m2
received within a 5.7o field of view
centered on the sun, during the 60
minutes preceding the hour indicated.
9999 = missing data.
Amount of solar radiation in Wh/m2
received from the sky (excluding the
solar disk) on a horizontal surface,
during the 60 minutes preceding the
hour indicated. 9999 = missing data.
== End of SAMSON version 1.1 parameters =====
Amount of sky dome (in tenths)
covered by clouds that prevent
observing the sky or higher cl
layers. 99 = missing data.
.oud
! 08 054-058 Dry Bulb Temperature
Dry bulb temperature in degrees C.
9999. = missing data.
060-064 Dew Point Temperature
Dew point temperature in degrees C.
9999. = missing data.
Relative humidity in percent.
999 = missing data.
! 11 070-073 Station Pressure
264
-------
! 18 110-11J
Wind Direction
Wind Speed
Visibility
Ceiling Height
Hourly Precipitation
Hourly Precipitation Flag
Table
Below
Present weather
Precipitable Water
Broadband Aerosol
Optical_Depth
Snow Depth
Days Since Last Snowfall 0-88
Wind direction in degrees.
(N = 0 or 360, E = 90, S = 180,
W = 270) 999 = missing data.
missing data.
0.0-160.9 Horizontal visibility in
Precipitable water in
millimeters. 9999 = missing data.
Broadband aerosol optical depth
(broadband turbidity) on the
day indicated.
99999. = missing data.
In inches and hundredths
(See information below).
See explanation below.
Loop: Do
Read (Jin, '(a)', lostat = ios) qline
If (ios /= 0) Goto 9999
nlines = nlines + 1
265
-------
! Print the Hourly Precipitation part of the record.
! This record is available only for vl.0 files.
!If (.Not. Is_vll) Then
! Write(1001, '(Ix,a,a,Ix,a)') 'SAMSON: ', qline(l:13), Trim(qline(123
!End If
Read (qline, 9190, lostat = ios) yy, mm, dd, hh, EHR, EDNR, &
GHR_v, GHR_s, GHR_u, DNR_v, DNR_s, DNR_u, DHR_v, DHR_s, DHR_u, &
Total_Sky_Cover, Opaque_Sky_Cover, &
Dry_Bulb_Temperature, Dew_Point_Temperature, &
Relative Humidity, Station Pressure, &
Wind_Direction, Wind_Speed, &
Horizontal_Visibility, Ceilinq_heiqht, Observation_Indicator, &
Present weather, Precipitable water, &
baod, Snow Depth, Days Since Last Snowfall, &
Hourly_Precipitation, Hourly_Precipitation_Flag
9190 Format ( 4(lx,12), 2(lx,14), &
3(lx,14,lx,al,il), &
2(lx,12), 2(lx,f5.1), Ix, 13, Ix, 14, &
Ix, 13, £5.1, £6.1, 16, Ix, 11, &
a9, 14, &
£6.3, 14, 13, &
Ix, 16, al )
MI *********************** Calendar and "iv" tests commented out on Wed Dec 12 12:29:37 2001, because
MI *********************** jfj__ All tests were passed without problems.
MI *********************** #2. With the tests active, it took 3m 21s to read all files;
iii *********************** With the tests commented out, it took 1m 36s to read all files;
jd_today = jd(yyyy, mm, dd)
!!! Call Jd_to_ymd(jd_today, jyyyy, jmm, jdd)
!!! ok = (yyyy == jyyyy) .And. (mm == jmm) .And. (dd == jdd)
!!! If (.Not. ok) Then
!!! Write(ULoq,0102) '?? jd /= Jd_to_ymd ', yyyy, mm, dd
!!! Write(6,0102) '?? jd /= Jd_to_ymd ', yyyy, mm, dd
! ! ! Stop '?? jd /= Jd_to_ymd'
! ! ! End If
266
-------
doy = iDoY(yyyy, mm, dd)
Call Calend(yyyy, doy, jmm, jdd)
ok = (mm == jmm) .And. (dd == jdd)
If (.Not. ok) Then
Write(ULog,0102) '?? IDoY /= Calend
Write(6,0102) '?? IDoY /= Calend
Stop '?? IDoY /= Calend1
End If
dd
dd
24
Nhours = 2J
! Test inverse.
jdoy = (iv+Nhours-1) / Nhours
jhh = iv - (jdoy-1)*Nhours
ok = (hh == jhh) .And. (doy == jdoy)
If (.Not. ok) Then
Write(ULog,0102) '?? iv /= (doy,hh) ', yyyy, mm, dd, hh
Write(6,0102) '?? iv /= (doy,hh) ', yyyy, mm, dd, hh
Stop '?? iv /= (doy,hh)'
End If
267
-------
! Print now = ( ( j vO <= iv) . And . ( iv <= j vl) )
!If (Print_now) Then
! Write (ULog, *) yyyy, mm, dd, hh, &
! ': Days_Since_Last_Snowfall = ', Days_Since_Last_Snowfall
!End If
I I I
I I I
I I I
I I I
I I I
I I I
! ! !0202
I I I
I I I
I I I
I I I
I I I
I I I
I I I
!!!0212
I I I
I I I
! Test iv to ymdh and inverse.
Call iv_to_ymdh(iv, jyyyy, jmm
ok = (yyyy == jyyyy)
If (.Not. ok) Then
Write(ULog,0202) '
Write(6, 0202) '
Format!///, lx, a,
Stop '?? iv to ymdh1
End If
Call ymdh_to_iv(yyyy,
ok = (iv == jv) .And.
If (.Not. ok) Then
Write(ULog,0212) '?? ymdh_to_iv: ',
Write(6, 0212) '?? ymdh_to_iv: ',
Format!///, lx, a, 'iv=', 10, '/=',
Stop '?? ymdh to iv1
End If
! Tally this day and hour.
p every hour(doy,hh) = p every hour(doy,hh) + 1
!WFlx%Every_Hour_Present(doy,hh) = 1 + &
! WFlx%Every_Hour_Present(doy,hh)
! It is daylight if Rs>0 (and Rs is not missing). (Rs is an alias for GHR v.)
! "is daytime" is used to determine if other parameters are truly
! missing, e.g., if a parameter X is to be output only for daylight hours,
! then parameter X is (effectively) non-missing if Rs>0 for the same record (hour).
! The default for "is daytime" is .True., i.e., any parameter will be missing
! unless we can determine it is nighttime. This will take care of the (rare) case
! when Rs is missing.
is_daytime = .True.
lobal Horizontal Radiation
Total amount of direct and diffuse
Data Value 0-1415
Flag for Data Source A-H, ?
Flag for Data Uncertainty 0-9
solar radiation in Wh/m2 received on
a horizontal surface during the 60
minutes preceding the hour indicated.
i 9999 = missing data.
! Ref[FAO] Rs == Global Horizontal Radiation == samson GHR_v == huswo iGRAD
! 1 Watt hour == 3.6e-3 MJoule
! 1 Watt hour m^-2 == 8.59845E-02 Langley
If ((0 <= GHR_v) .And. (GHR_v <= 1415)) Then
is_daytime = (GHR_v > 0)
If (Is_vll) Then
Xparam(f_GHR)%Samson_vll(iv)%v = GHR_v
268
-------
Xparam(f GHR)%Samson vll(iv)%s = data source
Else
Xparam(f_GHR)%Samson_vlO(iv)%v = GHR_v
Xparamtf_GHR)%Samson_vlO(iv)%s = data_source
Endif
If (Verify(GHR_s, 'abcdefghABCDEFGH?') == 0) Then
If (Is_vll) Then
Xparam(f_GHR)%Samson_vll(iv)%f(1:1) = GHR_s
Else
Xparam(f_GHR)%Samson_vlO(iv)%f(1:1) = GHR_s
Endif
Else
Write (ULog, 9290) iv, yyyy, mm, dd, hh, &
'GHR_s', GHR_s, '{abcdefghABCDEFGH?}'
End If
If ((0 <= GHR_u) .And. (GHR_u <= 9)) Then
If (Is_vll) Then
Xparamtf_GHR)%Samson_vll(iv)%f(2:2) = Achar(GHR_u+iaO)
Else
Xparam(f_GHR)%Samson_vlO(iv)%f(2:2) = Achar(GHR_u+iaO)
Endif
Else
Write (ULog, 9270) iv, yyyy, mm, dd, hh, &
1GHR_u', GHR_u, '{0..9}'
End If
End If
! 01 014-017 Extraterrestrial 0-1415 Amount of solar radiation in Wh/m2
! Horizontal Radiation received on a horizontal surface at
! the top of the atmosphere during the
! 60 minutes preceding the hour indicated.
! Ref[FAO] Ra == Extraterrestrial Horizontal Radiation
! 1 Watt hour == 3.6e-3 MJoule
If ((0 <= EHR) .And. (EHR <= 1415)) Then
If (Is_vll) Then
Xparam(f_EHR)%Samson_vll(iv)%v = EHR
Xparamtf_EHR)%Samson_vll(iv)%s = data_source
Else
Xparam(f_EHR)%Samson_vlO(iv)%v = EHR
Xparam(f_EHR)%Samson_vlO(iv)%s = data_source
Endif
End If
! Extraterrestrial Direct Normal Radiation
If ((0 <= EDNR) .And. (EDNR <= 1415)) Then
If (Is_vll) Then
Xparam(f_EDNR)%Samson_vll(iv)%v = EDNR
Xparam(f_EDNR)%Samson_vll(iv)%s = data_source
Else
Xparam(f_EDNR)%Samson_vlO(iv)%v = EDNR
269
-------
Xparam(f_EDNR)%Samson_vlO(iv)%s = data_source
Endif
End If
! Direct Normal Radiation
If ((0 <= DNR_v) .And. (DNR_v <= 1415)) Then
If (Is_vll) Then
Xparam(f_DNR)%Samson_vll(iv)%v = DNR_v
Xparam(f_DNR)%Samson_vll(iv)%s = data_source
Else
Xparam(f_DNR)%Samson_vlO(iv)%v = DNR_v
Xparamtf_DNR)%Samson_vlO(iv)%s = data_source
Endif
If (Verify(DNR_s, 'abcdefghABCDEFGH?') == 0) Then
If (Is_vll) Then
Xparam(f_DNR)%Samson_vll(iv)%f(1:1) = DNR_s
Else
Xparam(f_DNR)%Samson_vlO(iv)%f(1:1) = DNR_s
Endif
Else
Write (ULog, 9290) iv, yyyy, mm, dd, hh, &
'DNR_s', DNR_s, '{abcdefghABCDEFGH?}'
End If
If ((0 <= DNR_u) .And. (DNR_u <= 9)) Then
If (Is_vll) Then
Xparam(f_DNR)%Samson_vll(iv)%f(2:2) = Achar(DNR_u+iaO)
Else
Xparam(f_DNR)%Samson_vlO(iv)%f(2:2) = Achar(DNR_u+iaO)
Endif
Else
Write (ULog, 9270) iv, yyyy, mm, dd, hh, &
1DNR_u', DNR_u, '{0..9}'
End If
End If
! Diffuse Horizontal Radiation
If ((0 <= DHR_v) .And. (DHR_v <= 1415)) Then
If (Is_vll) Then
Xparam(f_DHR)%Samson_vll(iv)%v = DHR_v
Xparam(f_DHR)%Samson_vll(iv)%s = data_source
Else
Xparam(f_DHR)%Samson_vlO(iv)%v = DHR_v
Xparam(f_DHR)%Samson_vlO(iv)%s = data_source
Endif
If (Verify(DHR_s, 'abcdefghABCDEFGH?') == 0) Then
If (Is_vll) Then
Xparam(f_DHR)%Samson_vll(iv)%f(1:1) = DHR_s
Else
Xparam(f_DHR)%Samson_vlO(iv)%f(1:1) = DHR_s
Endif
270
-------
Else
Write (ULog, 9290) iv, yyyy, mm, dd, hh, &
'DHR_s', DHR_s, '{abcdefghABCDEFGH?}'
End If
If ((0 <= DHR_u) .And. (DHR_u <= 9)) Then
If (Is_vll) Then
Xparam(f_DHR)%Samson_vll(iv)%f(2:2) = Achar(DHR_u+iaO)
Else
Xparam(f_DHR)%Samson_vlO(iv)%f(2:2) = Achar(DHR_u+iaO)
Endif
Else
Write (ULog, 9270) iv, yyyy, mm, dd, hh, &
1DHR_u', DHR_u, ' {0 . . 9} '
End If
End If
! SAMSON version 1.1 files have no more parameters
If (Is vll) Cycle Loop
! Total Sky_Cover
If ((0 <= Total_Sky_Cover) .And. (Total_Sky_Cover <= 10)) Then
Xparam(f_TSC)%Samson_vlO(iv)%v = Total_Sky_Cover
Xparam(f_TSC)%Samson_vlO(iv)%s = data_source
End If
! 07 051-052 Opaque Sky_Cover 0-10 Amount of sky dome (in tenths)
! covered by clouds that prevent
! observing the sky or higher cloud
! layers. 99 = missing data.
If ((0 <= Opaque_Sky_Cover) .And. (Opaque_Sky_Cover <= 10)) Then
Xparam(f_OSC)%Samson_vlO(iv)%v = Opaque_Sky_Cover
Xparam(f_OSC)%Samson_vlO(iv)%s = data_source
End If
! 08 054-058 Dry Bulb Temperature -70.0 to Dry bulb temperature in degrees C.
! 60.0 9999. = missing data.
If ((-70.0 <= Dry_Bulb_Temperature) .And. &
(Dry_Bulb_Temperature <= +60.0)) Then
Xparam(f_DBT)%Samson_vlO(iv)%v = Dry_Bulb_Temperature
Xparam(f_DBT)%Samson_vlO(iv)%s = data_source
End If
! Dew Point Temperature
If ((-70.0 <= Dew Point Temperature) .And. &
(Dew_Point_Temperature <= +60.0)) Then
Xparam(f_DPT)%Samson_vlO(iv)%v = Dew_Point_Temperature
Xparam(f_DPT)%Samson_vlO(iv)%s = data_source
271
-------
Relative Humidity 0-100 Relative humidity in percent.
i 999 = missing data.
If ((0 <= Relative_Humidity) .And. (Relative_Humidity <= 100)) Then
Xparam(f_RH)%Samson_vlO(iv)%v = Relative_Humidity
Xparam(f_RH)%Samson_vlO(iv)%s = data_source
End If
! 11 070-073 Station Pressure 700-1100 File: Station pressure in millibars.
! 1 millibar = 0.10000 kilopascal Internal: kilopascal
If ((700 <= Station_Pressure) .And. (Station_Pressure <= 1100)) Then
Xparam(f_SP)%Samson_vlO(iv)%v = Station_Pressure * millibar to kilopascal
Xparam(f_SP)%Samson_vlO(iv)%s = data_source
End If
Wind direction in degrees.
! (N = 0 or 360, E = 90, S = 180,
! W = 270) 999 = missing data.
If ((0 <= Wind_Direction) .And. (Wind_Direction <= 360)) Then
Xparam(f_WD)%Samson_vlO(iv)%v = Wind_Direction
Xparam(f WD)%Samson vlO(iv)%s = data source
End If
! Wind Speed will be normalized to z=10 meters. See "Data Base Project Documentation".
If ((0.0 <= Wind_Speed) .And. (Wind_Speed < +99.0)) Then
If (Wind_Speed .NotEqual. 0.0) Then
ulO = Wind Speed * 5.81 / Log(Station elevation in meters/0.03)
Xparam(f_WS)%Samson_vlO(iv)%v = ulO
Else
Xparam(f_WS)%Samson_vlO(iv)%v = Wind_Speed
End If
Xparam(f_WS)%Samson_vlO(iv)%s = data_source
End If
!! too many missing.
!! 14 083-088 Visibility 0.0-160.9 Horizontal visibility in
!! kilometers. 777.7 = unlimited
!! visibility. 99999. = missing data.
If ((0.0 <= Horizontal_Visibility) .And. &
(Horizontal_Visibility < 99999)) Then
Xparam(f_HV)%Samson_vlO(iv)%v = Horizontal_Visibility
Xparam(f HV)%Samson vlO(iv)%s = data source
If ((0.0 <= Horizontal Visibility) .And. &
(Horizontal_Visibility <= 160.9)) Then
Maximum_Horizontal_Visibility = Max(Maximum_Horizontal_Visibility, Horizontal_Visibility)
Else If (Horizontal_Visibility .Eguals. 777.7) Then
272
-------
Xparam(f_HV)%Samson_vlO(iv)%s = T_Unlimited
End If
End If
! 15 089-094 Ceiling Height 0-30450 Ceiling height in meters.
! 77777 = unlimited ceiling height.
! 88888 = cirroform.
I gggggg = missing data.
If ((0 <= Ceiling_height) .And. (Ceiling_height < 999999)) Then
Xparam(f CH)%Samson vlO(iv)%v = Ceiling height
Xparam(f CH)%Samson vlO(iv)%s = data source
If ((0 <= Ceiling_height) .And. (Ceiling_height <= 30450)) Then
Maximum_Ceiling_Height = Max(Maximum_Ceiling_Height, Ceiling_height)
Else If (Ceiling height == 77777) Then
Xparam(f_CH)%Samson_vlO(iv)%s = T_Unlimited
Else If (Ceiling_height == 88888) Then
Xparam(f_CH)%Samson_vlO(iv)%s = T_Cirroform
End If
End If
! Observation Indicator 0 or 9 0 = Weather observation made.
! 9 = Weather observation not made or missing.
! Present weather - Present weather conditions denoted by 9 indicators.
Select Case(Observation_Indicator)
Case(0, 9)
! Obs made. Save the present weather table.
Xparam(f OI)%Samson vlO(iv)%v = Observation Indicator
Xparam(f_OI)%Samson_vlO(iv)%f = Present_weather
Xparam(f_OI)%Samson_vlO(iv)%s = data_source
! Check the present weather table.
Do i = 1, 9
If (Verify(Present_weather(i:i), PWT) == 0) Cycle ! Ok.
Write (ULog, 9250) iv, yyyy, mm, dd, hh, &
1Present_weather( ',i, ' ) ' , Present_weather(i:i), PWT
Format (/, Ix, '?? Read_SAMSON_File_vlx: ', &
17, Ix, 14, 2('-',i2.2), 13, 'h, ', a, 10, a, &
1: invalid value: "', a, '", expecting: "', a, '"' /)
End Do
Case Default
Write (ULog, 9270) iv, yyyy, mm, dd, hh, &
'Observation Indicator1, Observation Indicator, '{0 9}'
Format (/, Ix, '?? Read_SAMSON_File_vlx: ', &
17, Ix, 14, 2('-',i2.2), 13, 'h, ', a, &
1: invalid value: ', 10, ', expecting: ', a, /)
End Select
! Precipitable water
If ((0 <= Precipitable_water) .And. (Precipitable_water <= 100)) Then
273
-------
Xparam(f_pH2O)%Samson_vlO(iv)%v = Precipitable_water
Xparam(f_pH2O)%Samson_vlO(iv)%s = data_source
End If
! Broadband aerosol optical_depth (broadband turbidity) on the day indicated.
If ((0.0 <= baod) .And. (baod <= 0.900)) Then
Xparam(f baod)%Samson vlO(iv)%v = baod
Xparam(f_baod)%Samson_vlO(iv)%s = data_source
End If
! Snow Depth
If ((0 <= Snow_Depth) .And. (Snow_Depth <= 100)) Then
Xparam(f_SD)%Samson_vlO(iv)%v = Snow_Depth
Xparam(f SD)%Samson vlO(iv)%s = data source
End If
! Days Since Last Snowfall
If ((0 <= Days_Since_Last_Snowfall) .And. (Days_Since_Last_Snowfall <= 88)) Then
Xparam(f_DSLS)%Samson_vlO(iv)%v = Days_Since_Last_Snowfall
Xparamtf_DSLS)%Samson_vlO(iv)%s = data_source
End If
!If (Print_now) Then
! Write (ULog, ^) yyyy, mm, dd, hh, &
! ': Xparam(f_DSLS)%Samson_vlO(iv) = ', Xparam(f_DSLS)%Samson_vlO(iv)
!End If
! 21 124-129 Hourly Precipitation 000000- In inches and hundredths
i 099999 (See information below).
I
! 130-130 Hourly Precipitation Flag See explanation below.
I
! DATA FORMAT—HOURLY PRECIPITATION
I
! It stands to reason that for most hours the non-occurrence of
! precipitation is prevalent. Therefore, in order to save space in
! the original digital file, there are entries only for:
I
! 1. The first day and hour of each month where observations were
! taken even if no precipitation occurred during that month.
I
! 2. Hours with precipitation > zero.
I
! 3. Beginning and ending hours of missing periods.
I
! 4. Beginning and ending hours of accumulating periods.
I
! 5. Beginning and ending periods of deleted data.
274
-------
! The actual precipitation data value: The data value portion is a
! six-digit integer. Units are inches and hundredths. Range =
! 000000-099999. 000000 will be used only on the first hour of each
! month unless there is precipitation during that hour, in which case
! the measured value will be provided. On other days during the
! month without precipitation, no entry will be made. 099999
! indicates that the value is unknown.
I
! Hourly Precipitation Flag:
I
! A Accumulated period and amount. An accumulated period
! indicates that the precipitation amount is correct, but
! the exact beginning and ending times are only known to
! the extent that the precipitation occurred sometime
! withinthe accumulation period. Begin accumulation data
! value will always be 099999.
I
! D Deleted Flag. Beginning and ending of a deleted period.
! A deleted value indicates that the original data were
! received, but were unreadable or clearly recognized as
! noise.
I
! M Missing Flag. (Beginning and ending of a missing
! Period.) A missing flag indicates that the data were
! not received. This flag appears on the first and last
! dayofeach month for which data were not received or
! not processed by NCDC. Prior to 1984 a missing period
! was recorded as " OOOOOM" at the beginning and ending
! hours. If precipitation occurred during the last hour
! of the missing period, the second M appears with a non-
! zero value. Beginning in 1984 the beginning and ending
! hours of the missing period are recorded as "099999M".
I
! b Blank. No Flag needed.
I
! Examples:
I
! The precipitation accumulation from 1st month day 02 to 2nd month
i Hc,ir n/i -
Accumulation begins
099999A Accumulation continues
000390A Accumulation ends
! Accumulated precipitation for 1 month only:
I
! 01
275
-------
Accumulated, deleted, and missing precipitation data through months
01 and 02:
01 0001 0100 OOOOOOb First record of the
month
Accumulation begins
Accumulation continues
Accumulation ends
Deleted data begins
Deleted data ends
Missing data
Missing data
I
! Reguired precipitation charts or forms were never
I
! 01 0001
! Again Hourly Precipitation In inches and hundredths
! Example: WBAN 14914, 1961 range: 0 - 218; i.e., 218
If ((0 <= Hourly_Precipitation) .And. (Hourly_Precipitation < 099999)) Then
Xparam(f_HP)%Samson_vlO(iv)%v = inches to cm * Hourly_Precipitation / 100.0
Xparam(f_HP)%Samson_vlO(iv)%s = data_source
End If
! Valid range: ADM E(undocumented) Blank
Xparam(f_HP)%Samson_vlO(iv)%f = Hourly_Precipitation_Flag
Select Case(Hourly_Precipitation_Flag)
Case ( 'A', 'D', 'M', 'E', ' ' )
! Ok.
Case Default
! Invalid value.
Write (ULog, 9290) iv, yyyy, mm,
'Hourly Precipitation flag1
Hourly Precipitation Flag, &
1{A D M E Blank}'
Format (/, Ix, '?? Read_SAMSON_File_vlx: ', &
17, Ix, 14, 2('-',i2.2), 13, 'h, ', a, &
1: invalid value: "', a, '"; Expecting: ', a, /)
End Select
End Do Loop
9999 Continue
Call lOClose(jin)
Empty_File = (nlines == 0)
End Subroutine Read_SAMSON_vlx
End Module SAMSON
276
-------
setup
! Last change: LSR 17 May 2002 6:11 pm
Module Setup
Use Date_Module
Use F2kCLI
Use FileStuff
Use Global Variables
Use Stats
Use Utilsl
Use loSubs
Implicit None
Subroutine InitialSetUp()
Implicit None
Character(Len=30) :: ydate = ''
Character(Len=80) :: xexe = ''
Character(Len=80) :: xversion = ''
Integer :: nn
! Get the directory delimiter
DirDelim = DirectoryDelimiter()
ydate = Unix_Date()
Call Get Command Argument(0, xexe) ! Get application name.
! Remove possible path.
nn = 1 !+ Index(xexe, DirDelim, Back=.True.)
Call FileTimeStamp(xexe, xversion)
Call STD_dir(Log_dir)
Call STD_dir(R0_root)
Call STD dir(Raw Data dir)
! Open the log file.
Call Get Log File Name(Log file, Log dir, ULog)
! Dump data for Mathematica here
Call lOWrite(Umath, 'Omath.del')
277
-------
! Header for the log file: date and files used.
Write (ULog, *) 'Executable: ', Trim(xexe(nn:))
Write (ULog, *) 'Version ..:
Write (ULog, *) 'run of ...:
Write (ULog, *) 'TimeStamp :
Write (ULog, 9130)
9130 Format (/ &
Ix, 'Notes:', /, &
Ix, '* None as of 11 Feb 2002 2:27 pm', /)
! To Do List:
Write (ULog, 9150)
9150 Format ( &
Ix, 'To Do List: ', /, &
Ix, '* Look for "***" "!!" "!!!"', /)
Call SetFieldlnfo()
Call Set_Hours_since_JdO()
ICall TestJATinteracter ( )
ICall Test_Multiples_of()
ICall Test_Azimuth()
ICall Test_NaN()
I Call Test_iv_to_ymdh()
Call FLushAll()
End Subroutine InitialSetUp
Subroutine Check Drives and Files()
Implicit None
Character(Len=l0), Dimension(2), Parameter :: xDrive = (/ 'v:\.', 'z:\.' /)
Logical :: have drive
Logical :: have_file
Integer : : nerr, j , n
I Verify that the substituted drives are present.
Do j = 1, Ubound(xDrive,1)
Inquire(File=xDrive(j), Exist=have drive)
If (.Not. have_drive) Then
n = Len_trim(xDrive(j)) - 2 I Do not display trailing "\."
Write(ULog, 9130) xDrive(j)(1:n)
278
-------
Format (Ix, '?? Check Drives and Files: Could not find drive ', a)
nerr = nerr + 1
End If
End Do
! make sure "zcat" exists.
Inquire(File=zcat, Exist=have file)
If (.Not. have_file) Then
Write (ULog, *) '?? Check_Drives_and_Files: Could not find file zcat: ', Trim(zcat)
nerr = nerr + 1
End If
! make sure metadata file exists.
Inquire(File=metadata coda, Exist=have file)
If (have_file) Then
Call lORead(Umetadata, metadata_coda)
Else
Write (ULog, *) '?? Check_Drives_and_Files: Could not find file: ', Trim(metadata_coda)
nerr = nerr + 1
End If
If (nerr > 0) Then
Stop '?? Missing drives/files'
End If
End Subroutine Check Drives and Files
I ^^ insert subroutine in
!
! Generate log file name, e.g., 'Logs\2001-10-03_090547.log';
! If Log Path is present, it must terminate with the directory delimiter, 'Logs\';
! If ULog is present, Log_Name will be opened with write access and attached to ULog,
Implicit None
Character(Len=^), Intent(Out) :: Log Name
Character(Len=*), Optional, Intent(In) :: Log_Path
Integer, Optional, Intent(Out) :: ULog
Call Date and time(Date = xdate, Time = xtime)
279
-------
! tbuf 2001-10-03_090547.log
! xdate yyyymmdd
! xtime hhmmss.sss
! 123456789-123456789-1
tbuf = 'yyyy-mm-dd_hhmmss.log'
tbuf(01:04) = xdate(l:4) ! year
tbuf(06:07) = xdate(5:6) ! Month
tbuf(09:10) = xdate(7:8) ! day of Month
tbuf(12:17) = xtime(1:6) ! hhmmss
If (Present(Log_Path)) Then
tdir = Log_Path
Call STD_dir(tdir)
Log_Name = Trim(tdir) // tbuf
Else
Log_Name = tbuf
End If
If (Present(ULog)) Then
Call lOWrite(ULog, Log_Name)
End If
End Subroutine Get Log File Name
Subroutine SetFieldlnfo()
Implicit None
! Parameter Zero: Other counts
Fieldlnfo(0)%Name = 'Year counts'
Fieldlnfo(0)%Minimum obs per day = 0 ! Unused for this field.
! FAO Ra == Extraterrestrial Horizontal Radiation
! Extraterrestrial Horizontal Radiation - Amount of solar radiation in Wh/m2
! received on a horizontal surface at the top of the atmosphere during the
! 60 minutes preceding the hour indicated.
Fieldlnfo(f_EHR)%Name = 'Extraterrestrial Horizontal Radiation [Wh/m2] (Ra)'
Fieldlnfo(f_EHR)%Minimum_obs_per_day = 3
Fieldlnfo(f EHR)%minimum value = Zero
Fieldlnfo(f_EHR)%maximum_value = 1415 * 24 ! 24 hours
Fieldlnfo(f_EDNR)%Name = 'Extraterrestrial Direct Normal Radiation [Wh/m2]'
Fieldlnfo(f_EDNR)%Minimum_obs_per_day = 3
Fieldlnfo(f_EDNR)%minimum_value = Zero
Fieldlnfo(f EDNR)%maximum value = 1415 * 24 ! 24 hours
280
-------
! MET Solar Radiation == FAO Rs == Global Horizontal Radiation
! == samson GHRv == huswo iGRAD
! Global Horizontal Radiation - Total amount of direct and diffuse solar
! radiation in Wh/m2 received on a horizontal surface during the
! 60 minutes preceding the hour indicated.
Fieldlnfo(f_GHR)%Name = 'Global Horizontal Radiation [Wh/m2] (Rs)'
Fieldlnfo(f_GHR)%Minimum_obs_per_day = 3
Fieldlnfo(f_GHR)%minimum_value = Zero
Fieldlnfo(f_GHR)%maximum_value = 1415 * 24 ! 24 hours
Fieldlnfo(f_DNR)%Name = 'Direct Normal Radiation [Wh/m2]'
Fieldlnfo(f_DNR)%Minimum_obs_per_day = 3
Fieldlnfo(f_DNR)%minimum_value = Zero
Fieldlnfo(f_DNR)%maximum_value = 1415 * 24 ! 24 hours
Fieldlnfo(f_DHR)%Name = 'Diffuse Horizontal Radiation [Wh/m2]'
Fieldlnfo(f_DHR)%Minimum_obs_per_day = 3
Fieldlnfo(f DHR)%minimum value = Zero
Fieldlnfo(f_DHR)%maximum_value = 1415 * 24 ! 24 hours
! Total Sky Cover
Fieldlnfo(f_TSC)%Name = 'Total Sky Cover [tenths]'
Fieldlnfo(f_TSC)%Minimum_obs_per_day = 3
Fieldlnfo(f_TSC)%minimum_value = Zero
Fieldlnfo(f_TSC)%maximum_value = 10
! Opaque Sky Cover - Amount of sky dome (in tenths) covered by clouds that prevent
! observing the sky or higher cloud layers.
Fieldlnfo(f_OSC)%Name = 'Opaque_Sky_Cover [tenths]'
Fieldlnfo(f OSC)%Minimum obs per day = 3
Fieldlnfo(f OSC)%minimum value = Zero
Fieldlnfo(f_OSC)%maximum_value = 10
! Dry bulb temperature in degrees C.
Fieldlnfo(f_DBT)%Name = 'Dry Bulb Temperature [°C]'
Fieldlnfo(f_DBT)%Minimum_obs_per_day = 3
Fieldlnfo(f_DBT)%minimum_value = -70.0
Fieldlnfo(f_DBT)%maximum_value = +60.0
! Dew Point Temperature in degrees C.
Fieldlnfo(f_DPT)%Name = 'Dew Point Temperature [°C]'
Fieldlnfo(f_DPT)%Minimum_obs_per_day = 3
Fieldlnfo(f_DPT)%minimum_value = -80.0
Fieldlnfo(f_DPT)%maximum_value = +60.0
! Relative humidity in percent.
Fieldlnfo(f_RH)%Name = 'Relative_Humidity [Percent]'
Fieldlnfo(f_RH)%Minimum_obs_per_day = 3
Fieldlnfo(f_RH)%minimum_value = Zero
Fieldlnfo(f RH)%maximum value = 100
281
-------
! Station pressure in kPa.
Fieldlnfo(f_SP)%Name = 'Station_Pressure [kPa]1
Fieldlnfo(f_SP)%Minimum_obs_per_day = 1
Fieldlnfo(f_SP)%minimum_value = 70.0
Fieldlnfo(f_SP)%maximum_value = 110.0
! Wind direction in degrees.
! (N = 0 or 360, E = 90, S = 180, W = 270)
Fieldlnfo(f_WD)%Name = 'Wind_Direction [Degrees]'
Fieldlnfo(f_WD)%Minimum_obs_per_day = 3
Fieldlnfo(f WD)%minimum value = Zero
Fieldlnfo(f_WD)%maximum_value = 360
! Wind Speed in meters/sec at height = 10 meters
Fieldlnfo(f_WS)%Name = 'Wind_Speed @z=10m [m/s]'
Fieldlnfo(f_WS)%Minimum_obs_per_day = Fieldlnfo(f_WD)%Minimum_obs_per_day
Fieldlnfo(f_WS)%minimum_value = Zero
Fieldlnfo(f_WS)%maximum_value = 98.9
! Horizontal Visibility
Fieldlnfo(f_HV)%Name = 'Horizontal Visibility [km]'
Fieldlnfo(f HV)%Minimum obs per day = 3
Fieldlnfo(f HV)%minimum value = Zero
Fieldlnfo(f_HV)%maximum_value = Zero ! See Process_Set. Maximum value is station-dependent
! Ceiling Height
Fieldlnfo(f_CH)%Name = 'Ceiling Height [m]'
Fieldlnfo(f_CH)%Minimum_obs_per_day = 3
Fieldlnfo(f_CH)%minimum_value = Zero
Fieldlnfo(f CH)%maximum value = Zero ! See Process Set. Maximum value is station-dependent
! Observation Indicator 0 or 9 0 = Weather observation made.
! 9 = Weather observation not made or missing.
! Present weather - Present weather conditions denoted by 9 indicators.
Fieldlnfo(f_OI)%Name = 'Observation Indicator/Present weather'
Fieldlnfo(f_OI)%Minimum_obs_per_day = 1
Fieldlnfo(f_OI)%minimum_value = Zero ! N/A for this parameter
Fieldlnfo(f OI)%maximum value = Zero ! N/A for this parameter
! Precipitable Water
Fieldlnfo(f_pH2O)%Name = 'Precipitable Water [mm]'
Fieldlnfo(f_pH2O)%Minimum_obs_per_day = 3
Fieldlnfo(f pH2O)%minimum value = Zero
Fieldlnfo(f_pH2O)%maximum_value = 100
! Broadband Aerosol Optical Depth (0.0-0.90) - Broadband aerosol optical depth
! (broadband turbidity) on the day indicated.
Fieldlnfo(f_baod)%Name = 'Broadband Aerosol Optical_Depth []'
Fieldlnfo(f_baod)%Minimum_obs_per_day = 1
Fieldlnfo(f baod)%minimum value = Zero
282
-------
Fieldlnfo(f baod)%maximum value = 0.900
Fieldlnfo(f_SD)%Name = 'Snow Depth [cm]'
Fieldlnfo(f_SD)%Minimum_obs_per_day = 3
Fieldlnfo(f_SD)%minimum_value = Zero
Fieldlnfo(f_SD)%maximum_value = 100
Fieldlnfo(f_DSLS)%Name = 'Days since last Snowfall [days]'
Fieldlnfo(f_DSLS)%Minimum_obs_per_day = 3
Fieldlnfo(f_DSLS)%minimum_value = Zero
Fieldlnfo(f_DSLS)%maximum_value = 88
! Hourly Precipitation
Fieldlnfo(f_HP)%Name = 'Hourly Precipitation [cm]'
Fieldlnfo(f_HP)%Minimum_obs_per_day = 3
Fieldlnfo(f_HP)%minimum_value = Zero
Fieldlnfo(f_HP)%maximum_value =400
Fieldlnfo(f_FAO_SG_PET)%Name = 'Eto, FAO Short Grass [mm/day]'
Fieldlnfo(f_FAO_SG_PET)%Minimum_obs_per_day = 0
Fieldlnfo(f_FAO_SG_PET)%minimum_value = -1
Fieldlnfo(f_FAO_SG_PET)%maximum_value = 35
Fieldlnfo(f_KP_FWS_Evaporation)%Name = 'K-P FWS Evaporation [mm/day]'
Fieldlnfo(f_KP_FWS_Evaporation)%Minimum_obs_per_day = 0
Fieldlnfo(f_KP_FWS_Evaporation)%minimum_value = -10
Fieldlnfo(f_KP_FWS_Evaporation)%maximum_value = +20
Fieldlnfo(f_Ep)%Name = 'Ep, Class A pan Evaporation [mm/day]'
Fieldlnfo(f Ep)%Minimum obs per day = 0
Fieldlnfo(f Ep)%minimum value = -1
Fieldlnfo(f_Ep)%maximum_value = 180
! MET field information.
MET_field(g_Precipitation)%Name = 'Precipitation [cm/day]'
MET_field(g_Pan_Evaporation)%Name = 'Pan Evaporation [cm/day]'
MET_field(g_Temperature_mean)%Name = 'Temperature mean [°C]'
MET_field(g_Wind_Speed)%Name = 'Wind Speed @10 meter [cm/s]'
MET_field(g_Solar_Radiation)%Name = 'Solar Radiation [Langleys/day]'
MET_field(g_FAO_Short_Grass)%Name = Fieldlnfo(f_FAO_SG_PET)%Name
MET_field(g_Daylight_Station_Pressure)%Name = 'Daylight Station Pressure [kiloPascal]'
MET_field(g_Daylight_Relative_Humidity)%Name = 'Daylight Relative Humidity [%]'
MET_field(g_Daylight_Opague_Sky_Cover)%Name = 'Daylight Opague Sky Cover [tenths]'
MET_field(g_Daylight_Temperature)%Name = 'Daylight Temperature [°C]'
MET_field(g_Daylight_Broadband_Aerosol)%Name = 'Daylight Broadband Aerosol [optical depth]'
MET_field(g_Daylight_Mean_Wind_Speed)%Name = 'Daylight Mean Wind Speed @ 10 meters [m/s]'
MET field(d Daylight max wind speed)%Name = 'Maximum Daylight Mean Wind Speed @10m [m/s]'
MET_field(d_Daylight_direction_of_max_wind_speed)%Name = 'Direction of Maximum Daylight Wind [degre
MET_field(d_PWS)%Name = 'Daylight Prevailing Wind_Speed @z=10m [m/s]'
MET_field(d_PWD)%Name = 'Daylight Prevailing Wind_Direction [Degrees]'
283
-------
End Subroutine SetFieldlnfo
Subroutine TestjAfinteracter ( )
Use Winteracter
Implicit None
Integer :: uu = 6
Character(Len=200) :: tbuf, wbuf, tdir
Character(Len=200), Dimension(50) :: vfiles
Integer :: nn, j, ierr
Logical :: have_dir
Write (uu, *) ' *** Winteracter *** '
tdir = 's:\ '
have_dir = lOsDirExists(tdir)
Write (uu, ^) 'Have directory ', Trim(tdir), ' ? ', have dir
! Check create directory; Get rid of the path and create
! the new directory in the current directory.
! Warning: cannot create C:\a\b\c where \a\b does not exist.
Call Temp_File_Name(tdir)
j = Index(tdir, '\', Back=.True.) + 1
tdir = '000.' // tdir(j:)
If (.Not. lOsDirExists(tdir)) Then
Write (uu, ^ ) 'Directory ', Trim(tdir), ' did not exist; creating ....'
Call IQsDirMake(tdir)
have_dir = lOsDirExists(tdir)
Write (uu, *) 'Directory ', Trim(tdir), ' exists ? ', have dir
Else
Write (uu, *) 'Bummer. Directory ', Trim(tdir), ' exists; will not create.
End If
! %temp% == "e:\TEMP"
! Note: no trailing slash.
tbuf = 'temp'
Call IQsVariable(tbuf, wbuf) ! wbuf == e:\TEMP
Write (uu, 9130) Trim(tbuf), Trim(wbuf)
9130 Format (Ix, '%', a, '% == "', a, '"')
! does not look into subdirs.
Vfiles(l) = '**None**'
tdir = 'z:\'
tbuf = '14914_*.*'
nn = 1
284
-------
Do j = 1, 2
Call lOsDirEntryType('FD')
Call lOsDirlnfo(tdir, tbuf, vfiles, nn)
ierr = InfoError(l)
Write (uu, *) 'lOsDirlnfo: nn = ', nn, '; ierr
Write (uu, *) 'Vfiles(l) == "', Vfiles(l), '"'
nn = Ubound(vfiles,1)
End Do
Stop 'Stopping at Test_Winteracter'
End Subroutine Test Winteracter
Subroutine Test_Multiples_of()
Implicit None
Integer :: iNum
Integer :: uu = 9021
!! iNum =
Call ttO( 1, 24)
Call ttO(24, 25)
Call ttO(25, 25)
Call ttO(24,
Call ttO(48,
Call ttO(48,
Call ttO( 2, J
iNum = 3
Write (uu, *)
Call ttO(l, 2) ! 0
Call ttO(2, 4) ! 1
Call ttO(2, 8) ! 2
Call ttO(27, 29) ! 1
Call ttO(26, 33) ! 3
Call ttO(26, 34) ! 3
Call ttO(33, 33) ! 1
Call ttO(33, 39) ! 3
Stop 'Stopping at Test_Multiples_of'
Contains
285
-------
Implicit None
Integer, Intent(In) :: kO, kl
Integer :: n
End Subroutine ttO
End Subroutine Test_Multiples_of
Subroutine Test Azimuth()
! Wind Direction
! Azimuth: Wind direction in degrees.
! (N = 0 or 360, E = 90, S = 180, W = 270)
I
! Theta: regular angle measured counterclockwise
! from the +x axis.
I
! Interconversion Formulae:
! Theta = Modulo(90-Azimuth, 360)
! Azimuth = Modulo(90-Theta, 360)
I
! From Mathematica, for the test case: (ArcTan(x/y) + Pi) == Azimuth
! Note that Theta = ArcTan(y/x) + Pi
286
-------
1
I
I
I
1
1
I
I
1
1
I
I
1
1
I
I
1
1
I
I
1
1
I
I
1
1
Impli
Real
Real
Real
Real
Integ
Integ
Chara
Azimuth
0.0
15.0
30.0
45. 0
60. 0
75.0
90.0
105. 0
120. 0
135.0
150.0
165. 0
180. 0
195.0
210.0
225. 0
240.0
255.0
270.0
285. 0
300. 0
315.0
330.0
345. 0
360. 0
cit None
: : xO, xl,
Tt
(
]
3'
3:
3]
3t
2E
2'
2;
2^
2]
1<
If
If
i;
i:
i;
1C
c
xinc
: : aj ! Azimuth
: : tj ! Theta
: : bj ! A
er : : npts
er : : Jout
cter (Len=2
zimuth
, j
= 200
) : : tc
Theta Azimuth from Theta
90.0 0.0
75.0 15.0
60.0 30.0
45.0 45.0
30.0 60.0
15.0 75.0
0.0 90.0
5.0 105.0
330.0 120.0
315.0 135.0
300.0 150.0
5.0 165.0
270.0 180.0
255.0 195.0
240.0 210.0
225. 0
240.0
255.0
270.0
285. 0
300. 0
315.0
330.0
345. 0
0.0 ?? ! 0 is equivalent to 360. Result is ok.
tcoda
xO = 0
xl = 360
xinc = 15
npts = Ceiling((xl-xO)/xinc)
Write (Jout, 9130) 'Azimuth', 'Theta', 'Azimuth from Theta'
Format (Ix, alO, 3x, alO, 3x, a!8)
Do j = 0, npts
aj = xO + j ^xinc
tj = Modulo(90-aj, 360)
bj = Modulo(90-tj, 360)
If (Abs(aj-bj) <= le-6) Then
tcoda = ''
287
-------
Else
tcoda = '??'
End If
Write (Jout, 9150) aj, tj, bj, tcoda
9150 Format (Ix, flO.l, 3x, flO.l, 3x, f!8.1, 4x, a)
End Do
! Zero/Zero generates a NAN
Implicit None
If (IsNaN(r2)) Then
Write (6,*) 'Zero/Zero == NAN'
Else
Write (6,*) 'Zero/Zero /= NAN'
End If
! 0 1960-12-31 25h
! 1 1961-01-01 Ih
! 273924 1990-12-31 24h
i 999999 2070-07-07 24h
! Given iv, this subroutine computes the corresponding yyyy-mm-dd hh
! and, optionally, doy (day of the year).
288
-------
Implicit None
Integer : : yyyy, mm, old, hh
Integer :: jv, jdoy, nbase, ymax, i, nperiods
Write (ULog, *)
Write (ULog, *) 'jv : ', jv
Write (ULog, *) 'nperiods...: ', nperiods
Write (ULog, ^ ) ' ymax : ' , ymax
Do yyyy =
!! Hours_since_JdO(yyyy) = (Jd(yyyy,01,01) - JdO) * Nhours
!!nbase = Hours since JdO(yyyy)
nbase = (Jd(yyyy,01,01) - JdO) * Nhours
If (jv > nbase) Exit
End Do
Write (ULog, ^ ) ' nperiods ^Nhours : ' , nperiods ^Nhours
Write (ULog, *) '(Jd(yyyy,01,01) - JdO) : ', (Jd(yyyy,01,01) - JdO)
Write (ULog, *) '(Jd(yyyy,01,01) - JdO) * Nhours...: ', (Jd(yyyy,01,01) - JdO) * Nhours
289
-------
j v = j v - nbase
j doy = (j v+Nhours-1) / Nhours
End Do
End Subroutine Test iv to ymdh
End Module Setup
290
-------
Stats
! Last change: LSR 12 Mar 2002 9:05 am
Module Stats
Use Global Variables
! Numerically stable computation of the variance.
I
! 1 n
! Sample variance = Sum (x_i-x_mean)^2
! n - 1 1=1
! Reference:
! [2] Nicholas J. Higham. 1996. Accuracy and Stability of Numerical
! Algorithms. SIAM (Society for Industrial & Applied Mathematics
! ISBN 0-89871-355-2. Page 13.
I
! Accumulate:
! Note that the updating formulae can be written:
I
! M 0 = 0
291
-------
Implicit None
Public
Real, Private, Parameter :: rZero = 0.0
Subroutine Stat_Test(Uout, xOk)
! Simple test of module "Stats"
Implicit None
Integer, Intent(In)
Integer :: k
Real
Type(Stat_Block)
! Test set #1
Real, Dimension(3)
Real
Real
Eps
Xblo<
! Initialize accumulator block
Call Stat_Initialize(Xblock)
! Add points to accumulator
Do k = 1, Ubound(Setl,l)
Call Stat_Add_Point(Xblock, Setl(k)
End Do
Test set #1
292
-------
! Compute results
Call Stat Results(Xblock)
If (xOk) Then
Write(Uout, 9130) ' Setl'
Write(Uout, 9170) Xblock%xMean
Write(Uout, 9190) Xblock%xVariance
Else
Write(Uout, 9150) 'Setl1
Write(Uout, 9170) Xblock%xMean, Msetl
Write(Uout, 9190) Xblock%xVariance, Vsetl
End If
! Initialize accumulator block
Call Stat_Initialize(Xblock)
! Add points to accumulator
Do k = 1, Ubound(Set2,l)
Call Stat_Add_Point(Xblock, Set2(k)
End Do
Test set #2
! Epsilon(Vset2) = 1.19209290E-07
! Abs(Xblock%xMean-Mset2) = 2.38418579E-06
! Abs(Xblock%xVariance-Vset2) = 5.96046448E-08
Eps = l.Oe-5 ! Kludge
xOk = (Abs(Xblock%xMean-Mset2) < Eps) .And. &
(Abs(Xblock%xVariance-Vset2) < Eps)
IWrite(Uout,*) 'Eps = ', Eps
IWrite(Uout, *) 'Epsilon(Vset2) = ', Epsilon(Vset2)
IWrite(Uout,*) 'Abs(Xblock%xMean-Mset2) = ', Abs(Xblock%xMean-Mset2)
IWrite(Uout,*) 'Abs(Xblock%xVariance-Vset2) = ', Abs(Xblock%xVariance-Vseti
293
-------
End If
End Subroutine Stat_Test
Subroutine Stat Initialize(Xblock, Header)
Implicit None
Type(Stat_Block), Intent(Out) :
Character(Len=*), Optional, Intent(In) :
Xblock%k = 0
Xblock%M k = rZero
Xblock%Q k = rZero
Xblock%xMean = -Huge(rZero)
Xblock%xVariance = -Huge(rZero)
Xblock%xmin = +Huge(rZero)
Xblock%xmax = -Huge(rZero)
If (Present(Header)) Then
Xblock%Header = Header
Else
Xblock%Header = ''
End If
End Subroutine Stat Initialize
: Xblock
: Header
Implicit None
Type(Stat_Block), Target, Intent(InOut) ::
Real,
Intent(In)
Pointer
Real,
Real
Integer, Pointer
Q k =
k = k + 1
!If (k >= 2) Then
! M_kml = M_k
! M_k = M_kml + (x_k-M_kml)/Real(k)
! Q_k = Q_k + (Real(k-1)/Real(k))*(x_k-M_
!Else
294
-------
! ! k = 1
! M_k = x_k
! Q_k = rZero
!End If
M_kml = M_k
M_k = M_kml + (x_k-M_kml)/Real(k)
Q_k = Q_k + (Real(k-1)/Real(k))*(x_k-M_kml)**2
Xblock%xmin = Min(x_k, Xblock%xmin)
Xblock%xmax = Max(x_k, Xblock%xmax)
End Subroutine Stat Add Point
Implicit None
Type(Stat_Block), Target, Intent(InOut) :: Xblock
If (Xblock%k >= 2) Then
Xblock%xVariance = Xblock%Q_k/Real(Xblock%k-l)
Else
! Variance not defined for k<=l
Xblock%xVariance = -Huge(rZero)
End If
! Output data descriptive statistics.
! Calls: Stat_Results(Xblock)
Implicit None
Integer, Intent(In) :: Uout
Type(Stat_Block) , Intent(InOut) :: Xblock
Character(Len=*), Optional, Intent(In) :: Header
Real :: std dev ! Sample standard deviation
Write(Uout, *)
If (Present(Header)) Then
295
-------
Write(Uout, 9130) Trim(Header)
Else If (Len_trim(Xblock%Header) > 0) Then
Write(Uout, 9130) Trim(Xblock%Header)
End If
9130 Format (Ix, a)
Write(Uout, 9150) Xblock%k, Xblock%xMean, Xblock%xVariance
9150 Format ( &
Ix, ' n ....: ', 10, /, &
Ix, ' Mean .: ', Ipgl4.6, /, &
Ix, ' Var ..: ', Ipgl4.6)
If (Xblock%xVariance > rZero) Then
std dev = Sqrt(Xblock%xVariance)
Write(Uout, 9170) std_dev
Format (Ix, ' StdDev: ', Ipgl4.6)
End If
End Subroutine Stat_Output
End Module Stats
296
-------
UtilsO
Use Date_Module
Use Global_Variables
Use loSubs
Use Strings
Implicit None
! Skip lines until we find a line starting with Xstring;
! All errors are fatal.
Implicit None
Character(Len=^), Intent(In) :: FileName
Character(Len=^), Intent(In) :: Xstring
Integer, Intent(In) :: Jin
Character(Len=132) :: tbuf
Integer :: ios, lien
Logical :: Found_header_line
Errors Detected = .False.
Found header line = .False.
lien = Len(Xstring)
GetHeader: Do
Read (Jin, '(a)', lostat = ios) tbuf
If (ios /= 0) Exit GetHeader
If (tbuf(l:Ilen) == Xstring) Then
! Found the header line.
Found header line = .True.
Exit
End If
End Do GetHeader
If (.Not. Found_header_line)
Errors_Detected = .True.
Write (6, 9130) Xstring, Trim(FileName)
Write (ULog, 9130) Xstring, Trim(FileName)
Format (Ix, '?? Skip_Until: Could not find "', a,
Stop '?? Skip_Until: Could not find header line'
!Return
297
-------
End If
End Subroutine Skip_Until
Subroutine Roundoff(Xval, Nval, Wformat, &
Target total, Xeps, Yval, Delta_Sum)
Statement of the problem: The problem when we distribute
hourly precipitation among hours. The total amount to be
distributed is given (Target_total, 2.3622). When the
partial amounts are printed (.e.g, ' ( f 10 • 2 ) ' ) , the total
sum of the printed values (xtotal, 2.35) is different.
This subroutine attempts to
Xval
HP(full)
Implicit None
Real, Dimension(:
Integer,
Character(Len=*),
Real,
Real,
Real, Dimension!:
Real,
tval
HP(f0.2)
Intent(In)
Intent(In)
Intent(In)
Intent(In)
Intent(In)
Pointer
Intent(Out)
Xval
Nval
Wformat
Target_total
Xeps
Yval
Delta Sum
tbuf
xtotal = Zero
n = Len_trim(Wformat)
Do i = 1, Nval
298
-------
Write(tbuf(i), Wformat(1:n)) Xval(i)
Read(tbuf(i), *) Yval(i)
xtotal = xtotal + Yval(i)
End Do
Delta_Sum = Target_total - xtotal
If (Abs(Delta_Sum) > Xeps) Then
Yval(l) = Yval(l) + Delta_Sum
End If
End Subroutine Roundoff
End Module UtilsO
299
-------
Utilsl
Use Date_Module
Use FileStuff
Use Global_Variables
Use loSubs
Use Strings
Use Winteracter
Use Utils5
Implicit None
Contains
Recursive Subroutine Display Station Info(Xstations, Xfull)
Implicit None
Type(Site_Info), Pointer :: Xstations ! Pointer to root of tree
Logical, Optional, Intent(In) :: Xfull
! Visit smaller values first.
Call Display_Station_Info(Xstations%pLeft, Xfull)
! Process and write the current node.
Call Display_This_Station(Xstations, Xfull)
! Then visit larger values.
Call Display_Station_Info(Xstations%pRight, Xfull)
End If
End Subroutine Display Station Info
Subroutine Display_This_Station(xWBAN, Xfull)
Implicit None
Type(Site_Info), Intent(In) :: xWBAN
Logical, Optional, Intent(In) :: Xfull
300
-------
Integer : : i, yyyy, mm, dd
Logical : : all info
If (Present (Xfull) ) Then
all_info = Xfull
Else
all info = .False.
End If
Write (ULog, 9150) Trim (xWBAN%WBAN )
Format (Ix, '@@@ Station Name: ', a)
Write (ULog, *) 'Latitude : ', xWBAN%Lat
Write (ULog, *) 'Longitude: ', xWBAN%Lon
Write (ULog, 9190) xWBAN%Elev
Format (Ix, 'Elev: ', Ipgl4.6)
Write (ULog, 9230) xWBAN%Nelev
Format (Ix, 'Nelev: ', 10)
Do i = 1, xWBAN%Nelev
Call Jd_to_ymd(xWBAN%Elev_Directives(i)%Julian_Day, yyyy, mm, dd)
Write (ULog, 9250) i, yyyy, mm, dd, xWBAN%Elev_Directives(i)%Elevation_meter
Format (Ix, ' [', 10, '] ', 14, '-', 12.2, '-', 12.2, ': ', Ipgl4.6)
End Do
End If
If (xWBAN%Nelev == 0) Then
Errors_Detected = .True.
Write (ULog, 9270) Trim(xWBAN%WBAN), Trim(xWBAN%Text)
Format ('?? Station without elevation data: ', 3x, a, 3x, a)
End If
End Subroutine Display_This_Station
301
-------
! Propagate modifications to Allocate_SAMSON_arrays and Deallocate_SAMSON_arrays
! *
! *
! 0.0/0.0 generates a NaN. I am using the NaN so that the program will
! generate a message if the entry is used without being set first.
! Good idea, but what about missing entries in the SAMSON file?
! Missing entries are not read ==> value is NaN.
Implicit None
! Nelements: used for testing Senegal. We need an array of an
! specific size
Integer, Optional, Intent(In) :: Nelements
If (Present(Nelements)) Then
M elems = Nelements
Else
M_elems = (jdl-jdO+1) * Nhours
End If
! Samson vlO
Allocate(Xparam(jpar)%Samson_vlO(M_elems))
Xparam(jpar)%Samson_vlO%s = T_Missing
Xparam(jpar)%Samson vlO%v = Missing Data
!Xparam(jpar)%Samson vlO%v = Zero/Zero ! NaN
Xparam(jpar)%Samson_vlO%f = ''
! Samson vll
If ((f_EHR <= jpar) .And. (jpar <= f_DHR)) Then
Allocate(Xparam(jpar)%Samson_vll(M_elems))
Xparam(jpar)%Samson_vll%s = T_Missing
Xparam(jpar)%Samson vll%v = Missing Data
!Xparam(jpar)%Samson vll%v = Zero/Zero ! NaN
Xparam(jpar)%Samson_vll%f = ''
End If
End Do
302
-------
End Subroutine Allocate SAMSON arrays
vs .
! Deallocate SAMSON arrays
! Propagate modifications to Allocate_SAMSON_arrays and Deallocate_SAMSON_arrays
Implicit None
Integer :: jpar
! Parameter zero is special:
Xparam(0)%Samson_vlO => Null() ! Unused for param zero
Xparam(0)%Samson vll => Null() ! Unused for param zero
! Deallocate arrays for each parameter.
Do jpar = 1, f_end
Deallocate(Xparam(jpar)%Samson vl0)
If ((f_EHR <= jpar) .And. (jpar <= f_DHR)) Then
Deallocate (Xparam (jpar) %S amis on vll)
End If
End Do
Implicit None
Integer, Intent(Out) :: Jin
Character(Len=*), Intent(In) :: Zfile ! e.g., d:\AHFiles\13873_91.
! Decompress data file and open it. See Assumption[5].
Call FileNameParts(Zfile, ypath, yname, ytype)
Call Temp File Name(tname, NamePrologue=yname)
303
-------
Call IORead(Jin, tname)
End Subroutine Decompress and Open File
Subroutine Compare_SAMSON_Headers(WF10, WF11)
! Compare SAMSON vl.0 and vl.1 headers.
Implicit None
Type(Filelnfo), Intent(In) :: WF10, WF11
Logical :: ok
If (WF10%Head%WBAN /= WFll%Head%WBAN) Then
Write (ULog, '(lx,4a)') '?? WBANs do not match: ', &
Trim(WF10%Head%WBAN), '; ', Trim(WFll%Head%WBAN)
End If
If (.Not. String_Eq(WF10%Head%Text, WFll%Head%Text)) Then
Write (ULog, '(lx,4a)') '?? WBANs do not match: ', &
Trim(WF10%Head%Text), '; ', Trim(WFll%Head%Text)
End If
ok = (String_Eq(WF10%Head%Lat%Letter, WFll%Head%Lat%Letter)) .And. &
(WF10%Head%Lat%degrees == WFll%Head%Lat%degrees) .And. &
(WF10%Head%Lat%minutes == WFll%Head%Lat%degrees)
If (.Not. ok) Then
Write (ULog, '(a,2(2x,a2,10,Ix,10,:,"; "))') '?? Latitudes do not match: ', &
WF10%Head%Lat%Letter, WF10%Head%Lat%degrees, WF10%Head%Lat%minutes, &
WFll%Head%Lat%Letter, WFll%Head%Lat%degrees, WFll%Head%Lat%minutes
End If
ok = (String_Eq(WF10%Head%Lon%Letter, WFll%Head%Lon%Letter)) .And. &
(WF10%Head%Lon%degrees == WFll%Head%Lon%degrees) .And. &
(WF10%Head%Lon%minutes == WFll%Head%Lon%degrees)
If (.Not. ok) Then
Write (ULog, '(a,2(2x,a2,10,Ix,10,:,"; "))') '?? Longitudes do not match: ', &
WF10%Head%Lon%Letter, WF10%Head%Lon%degrees, WF10%Head%Lon%minutes, &
WFll%Head%Lon%Letter, WFll%Head%Lon%degrees, WFll%Head%Lon%minutes
End If
If (Abs(WF10%Head%Elev - WFll%Head%Elev) > EpsO) Then
! WF10%Head%Elev /= WFll%Head%Elev
Write (ULog, '(a,Ipgl4.6,3x,Ipgl4.6)') '?? Elevations do not match: ', &
WF10%Head%Elev , WF10%Head%Elev
End If
304
-------
End If
End Subroutine Compare_SAMSON_Headers
Implicit None
Type(Site_Info), Target, Intent(In) :: xWBAN
Integer, Intent(In) :: Julian_Day
Logical, Intent(Out) :: Xfound
Real, Intent(Out) :: Elevation
If (Nelev == 0) Return
! The Julian Day limits in Elev Directives are such that
! Note
! * Jd(i) = Elev_Directives(i)%Julian_Day
! ^ Jd(Nelev+l) is defined to be +Infinity, i.e., present day
Do ielev = Nelev, 1, -1
If (Elev_Directives(ielev)%Julian_Day <= Julian_Day) Then
Xfound = .True.
Elevation = Elev_Directives(ielev)%Elevation_meter
Return
305
-------
End If
End Do
End Subroutine GetElevation
Subroutine Set Hours since JdO()
! See
Implicit None
Integer :: y4
If (HsJDO_unset) Then
HsJDO_unset = .False.
Do y4 = MinYear, MaxYear
HsJDO(y4) = (Jd(y4,01,01) - JdO) * Nhours
End Do
End If
I
! Examples: Assume JdO = Julian_Day 1961-01-01
! YYYY Hours_since_JdO
! 1961 0
! 1962 365 * Nhours
! 1963 365*2 * Nhours
! 1964 365*3 * Nhours
! 1965 (365*3 + 366) * Nhours
I ; ;
I
! See
Implicit None
Integer, Intent(In)
Integer
306
-------
End Function Hours since JdO
Implicit None
Integer, Intent(In) :: jday
Character(Len=l8) :: str Jd to ymd
! 1
! 123456789012345678
! 2451545 2000-01-01
str Jd to ymd = 'iiiiiii yyyy-mm-dd'
Call Jd to ymd (jday, j yyyy, jmm, jdd)
xstr = itoa(j mm)
str_Jd_to_ymd(14:15) = Adjustr(xstr(1:2))
If (str_Jd_to_ymd(14:14) == '') Then
str_Jd_to_ymd(14:14) = '0'
End If
xstr = itoa(jdd)
str_Jd_to_ymd(17 :18) = Adjustr(xstr(1:2))
If (str_Jd_to_ymd(17:17) == '') Then
str_Jd_to_ymd(17:17) = '0'
End If
End Function str Jd to ymd
Function str iv to ymdh(iv)
307
-------
Implicit None
Integer, Intent(In) :: iv
Character(Len=21) :: str_iv_to_ymdh
Character(Len(str_iv_to_ymdh)) :: xstr
Integer :: jyyyy, jmm, jdd, jhh
! 1 2
! 123456789012345678901
! 223676 1985-07-01 Ih
str iv to ymdh = 'iiiiii yyyy-mm-dd HHh'
Call iv_to_ymdh(iv, jyyyy, jmm, jdd, jhh)
xstr = itoa (jmm)
str_iv_to_ymdh(13:14) = Adjustr(xstr(1:2))
If (str_iv_to_ymdh(13:13) == '') Then
str_iv_to_ymdh(13:13) = '0'
End If
xstr = itoa(jdd)
str_iv_to_ymdh(16:17) = Adjustr(xstr(1:2))
If (str_iv_to_ymdh(16:16) == '') Then
str_iv_to_ymdh(16:16) = '0'
End If
End Function str iv to ymdh
308
-------
24 25 Nhours = 2J
! See
! iv Date nperiods ymax
I
! 0 1960-12-31 25h
! 1 1961-01-01 Ih
! 273924 1990-12-31 24h
I gggggg 2070-07-07 24h
Implicit None
Integer, Intent(In) :
Integer, Intent(Out) :
Integer, Optional, Intent(Out) :
Integer :: jv, jdoy, nbase
! Get the year
i v = i v
Do yyyy = MaxYear, MinYear, -1
!nbase = Hours_since_JdO(yyyy)
nbase = HsJDO(yyyy)
If (jv > nbase) Exit
End Do
If (yyyy < MinYear) Then
Stop '?? Internal error in iv_to_ymdh: yyyy < MinYear1
End If
309
-------
24 25 Nhours = 25
9124
9149
Implicit None
Integer, Intent(In)
Integer, Intent(Out)
Integer, Optional, Intent(Out)
!nbase = Hours since JdO(yyyy)
nbase = HsJDO(yyyy)
iv = (jdoy-1)*Nhours + hh + nbase
End Subroutine ymdh to iv
310
-------
! Compute the number of multiples of INum that appear
! in the interval [Kmin, Kmax], 0 < Kmin <= Kmax.
Nmultiples
0
1
2
1
! 33 33 1
! 33 39 3
Implicit None
Integer, Intent(In) :
Integer :
Integer :: i, j
Logical :: ok
Real :: rnum
! If iNum == 0 we should not be here in the first place.
If (iNum == 0) Return
311
-------
End Function Multiples_of
Function Pythag(a, b) Result(Rval)
! Computes Sgrt(a^*2 + b*^2) without destructive overflow or underflow.
I
! History:
! = [Isr] Tue Dec 18 16:37:05 2001
! . processed by "to_f90" on Tue Dec 18 16:37:05 2001
! . Adapted from SLATEC/Pythag; ported to f95
I
!***LIBRARY SLATEC
l***TYpE SINGLE PRECISION (Pythag-S)
!***AUTHOR (UNKNOWN)
!***REVISION HISTORY (YYMMDD)
! 811101 DATE WRITTEN
! 890531 Changed all specific intrinsics to generic. (WRB)
! 891214 Prologue converted to Version 4.0 format. (BAB)
! 900402 Added TYPE section. (WRB)
Implicit None
Real, Intent(In) :: a
Real, Intent(In) :: b
Real : : Rval
Real ::p,g,r,s,t
Real, Parameter :: Two = 2.0eO
Real, Parameter :: Four = 4.OeO
Real, Parameter :: eps = Tiny(Two)
! The algorithm converges in at most three iterations,
! assuming unit roundoff >= le-20.
If (g <= eps) Then ! (g == Zero) ?
Rval = p
Return
End If
312
-------
Do
r = (q/p) ** 2
t = Four + r
If (Abs(t-Four) <= eps) Then ! (t == Four) ?
Rval = p
Exit
End If
s = r / t
p = p + Two*p*s
q = q * s
End Do
End Function Pythag
Implicit None
Integer, Intent(In) :: TimeZone
Integer :: Central Meridian
Select Case(TimeZone)
Case(4) ! Letter Q, Atlantic
Central Meridian = 60 ! 60W
Case(5) ! Letter R, Eastern
Central Meridian = 75
Case(7) ! Letter T, Mountain
Central_Meridian = 105
Case(8) ! Letter U, Pacific
Central Meridian = 120
Case(10) ! Letter W, Hawaii-Aleutian
Central Meridian = 150
313
-------
Central_Meridian = -150 ! -150W is equivalent to 150E
Case Default
Central_Meridian = TimeZone * 15
End Select
End Function TimeZone to Central Meridian
Subroutine Es_and_Delta(Tc, e_s, Delta)
! Tc Temperature [°C]
! e s saturation water vapor pressure [kPa] at temperature Tc; Ref[1:36]
! Delta slope of saturation water vapor pressure [kPa/°C]
Implicit None
Real, Intent(In) :: Tc
Real, Intent (Out) :: e_s
Real, Intent (Out) :: Delta
! Propagate changes in the parameters
! to Es_and_Delta and DewPointF
! [FAO, Page 36, eq. 11]
I
! eO(T) = es_cl * Exp(es_c2 * T / (es_c3 + T))
tmpC = Exp(es_c2 * Tc / (es_c3+Tc))
e_s = es_cl * tmpC
Delta = es_cl * es_c2 * es_c3 * tmpO / (es_c3+Tc)**2
Cp * P
eps ^ lambda
314
-------
Real, Parameter :: Cp = 1.013 ! KJ Kg~-l K~-l
Real, Parameter :: Eps = 0.622 ! dimensionless
Real :: LambdaF
End Function GammaF
!Function LambdaF(Tdew) Result(LambdaV)
I
! !
I I
! ! Tdew Dew point temperature [ °C]
! ! LambdaV Latent heat (enthalpy) of vaporazation [KiloJoules/Kg]
! Implicit None
! Real, Intent(In) :: Tdew
! Real :: LambdaV
I
! LambdaV = 2501 - 2.361*Tdew
I
!End Function LambdaF
Function Wind_Speed_F(ulO, T_Height) Result(ux)
! ulO Wind Speed in [m/s] at z=10 meters
! T_Height Height indicator
! ux Wind Speed in [m/s] at z indicated by T_Height
Implicit None
Real, Intent(In) :: ulO
Integer, Intent(In) :: T_Height
Real :: ux
!
! On input, Wind Speed in meters/sec was normalized to z=10 meters.
! Height Height indicator h x
315
-------
5.81 == Ln((10-0)/O.03)
4.87
4.89
3.56
1. 05
! u_c : wind speed at height z_c (m)
! d c : zero plane displacement (m)
! z Oc: surface roughness length or roughness height (m)
I
! For Open Flat Terrain (used for Metereological Stations):
! zO = 0.03 meters
! dO = 0
! z = 10 meters (reference height)
tu = ulO / 5.81
Select Case(T_Height)
Case(T_ulO)
! Height = 10 m
ux = ulO
Case (T_u4 )
! Height = 4 m
ux = tu * 4.89
Case (T_upl )
! Height = 0.1 m
ux = tu * 1 . 05
Case Default
Write (6, *) '?? Wind_Speed_F: unknown wind indicator == ', T_Height
316
-------
Write (ULog, *) '?? Wind_Speed_F: unknown wind indicator == ', T_Height
Stop '?? Wind Speed F: unknown wind indicator1
End Select
End Function Wind_Speed_F
! Estimate dewpoint from inversion of eq. 14, FAO, Page 37.
I
! Ta Air temperature [°C]
! RH Relative humidity [%]
! Tdew Dew point temperature [°C]
Implicit None
Real, Intent(In) :: Ta
Real, Intent(In) :: RH
Real :: Tdew
! e_s saturation water vapor pressure [kPa] at temperature Tc; Ref[l:36]
! e_a Atmospheric water vapor pressure; actual water vapor pressure, [kPa]
! Delta slope of saturation water vapor pressure [kPa/°C]
Real :: e_s, Delta, e_a
! e_s = eO(Ta)
Call Es and Delta(Ta, e s, Delta) ! Delta unused.
! Propagate changes in the parameters
! to Es_and_Delta and DewPointF
! [FAO, Page 36, eq. 11]
317
-------
Subroutine Generate File names()
! Generate the names of all files associated with the current wban number.
!
!
Implicit None
Character(MaxNamLen) :: tdir, tname
Integer :: yyyy, n
Character(Len(pWBAN%WBAN)+1) :: W_Wban
Character(Len=5):: W_State
! Output directory of the form: v:\rO\\14914\
! where is the two letter abbreviation of
! state where the station is located, .e.g,
! v:\rO\AK\25501\
W_State = pWBAN%State
R0_full = Trim(RO_root) // Trim(W_State) // DirDelim // Trim(pWBAN%WBAN)
Call STD_dir(R0_full)
tdir = RO full
End Do
! name_met == v:\rO\ND\14914\wl4914.dvf
! name_txt == v:\rO\ND\14914\wl4914.txt
tname = Trim(tdir) // W_Wban(l:n)
name_met = Trim(tname) // '.dvf
name_txt = Trim(tname) // '.txt'
! Generate filenames:
! Daily Precipitation files:
! v:\precip\03940.d
! v:\precip\03940g.d ! with gaps in the yearly record.
! Hourly Precipitation files:
! v:\precip\03940.h
! v:\precip\03940g.h
318
-------
! Daily evaporation file:
! v:\evaporation\T_14914.evp
name Daily Evap = Trim(Raw Data dir) // 'evaporation1 // &
DirDelim // 'T_' // Trim(pWBAN%WBAN) // '.evp'
End Subroutine Generate File names
Subroutine Initialize Output Directory()
Character(Len(pWBAN%WBAN)) :: Wban
Character(MaxNamLen) :: tdir
Integer :: yyyy, kO, kl, tlen
Logical :: have file
! If the output directory does not exist, exit.
! Output directory of the form: v:\r0.by.State\AK\25501\
! IQsDirMake will not create directory path x:\a\b\c\
! if, for example, directories a or b do not exist.
! We have to check and create each component individually — bummer.
tdir = R0_full
Call STD_dir(tdir) ! Make sure tdir has a trailing delimiter
tlen = Len_Trim(tdir)
kl = 1
do
If (kl > tlen) Exit
If (.Not. lOsDirExists(tdir(l:kO))
Call IQsDirMake(tdir(l:kO))
If (.Not. lOsDirExists(tdir(l:kO))) Then
Write(ULog, *) '?? Could not create Directory ', Trim(tdir(1:kO)
Stop '?? Stopping -- Could not create Directory.'
!Return
319
-------
End If
End If
kl = kO + 1 ! Skip over previous delimiter
End Do
! Hourly Values Files
Do yyyy = MinYear, MaxYear
Inquire(File=name_rO(yyyy), Exist=have_flie)
If (have_file) Then
Call lOsDeleteFile(name_rO(yyyy))
End If
End Do
! Later vintage: Winteracter does not complain if
! the file does not exist.
Call lOsDeleteFile(name_met)
Call lOsDeleteFile(name txt)
! http://www.spc.noaa.gov/fag/tornado/beaufort.html
! gyre.umeoce.maine.edu/data/gomoos/php/variable description.php?variable=wind speed
! -- Warning: the m/s quantities in gyre.umeoce.maine.edu are wrong!
! http://whale.wheelock.edu/whalenet-stuff/beaufort.html
Implicit None
Real, Intent(In) :: Wind Speed ! in meters/second
Integer, Intent(Out) :: Vforce ! Beaufort Force
Character(Len=*), Intent(Out) :: Vtext ! explanation
Type :: BF
Real :: Wind_Speed ! Lower limit of range, in meters/second
Character(Len=80) :: WMO_Classification = ''
Character(Len=80) :: SeaText = ''
Character(Len=80) :: LandText = ''
End Type BF
Integer :: i
Type(BF), Dimension(0:12), Parameter :: VInfo = (/ &
BF( 0.0, 'Calm', 'Sea surface smooth and mirror-like', 'Smoke rises vertically'), &
BF( 0.3, 'Light air', 'Scaly ripples, no foam crests', &
'Smoke drift indicates wind direction, still wind vanes'), &
BF( 1.6, 'Light Breeze', 'Small wavelets, crests glassy, no breaking', &
'Wind felt on face, leaves rustle, vanes begin to move'), &
BF( 3.4, 'Gentle Breeze', 'Large wavelets, crests begin to break, scattered whitecaps', &
'Leaves and small twigs constantly moving, light flags extended '), &
320
-------
BF( 5.5, 'Moderate Breeze1, 'Small waves 1-4 ft. becoming longer, numerous whitecaps', &
'Dust, leaves, and loose paper lifted, small tree branches move'), &
BF( 8.0, 'Fresh Breeze', 'Moderate waves 4-8 ft taking longer form, many whitecaps, some spray', &
'Small trees in leaf begin to sway'), &
BF(10.8, 'Strong Breeze', 'Larger waves 8-13 ft, whitecaps common, more spray', &
'Larger tree branches moving, whistling in wires'), &
BF(13.9, 'Moderate Gale', 'Sea heaps up, waves 13-20 ft, white foam streaks off breakers', &
'Whole trees moving, resistance felt walking against wind'), &
BF(17.2, 'Fresh Gale', 'Moderately high (13-20 ft) waves of greater length, edgesS
& of crests begin to break into spindrift, foam blown in streaks', &
'Twigs broken off trees, walking against wind very difficult'), &
BF(20.8, 'Strong Gale', 'High waves (20 ft), sea begins to roll, denseS
& streaks of foam, spray may reduce visibility', &
'Slight structural damage occurs, slate blows off roofs'), &
BF(24.5, 'Storm', 'Very high waves (20-30 ft) with overhanging crests,&
& sea white with densely blown foam, heavy rolling, lowered visibility', &
'Seldom experienced on land, trees broken or uprooted, "considerable structural damage"'), &
BF(28.5, 'Violent Storm', 'Exceptionally high (30-45 ft) waves, foam patchesS
& cover sea, visibility more reduced', &
'Widespread damage, very rare occurence'), &
BF(32.7, 'Hurricane', 'Air filled with foam, waves over 45 ft, sea completelyS
& white with driving spray, visibility greatly reduced', &
'Violent destruction') /)
Do i = Ubound(VInfo,l), Lbound(VInfo,1), -1
If (VInfo(i)%Wind_Speed <= Wind_Speed) Then
Vforce = i
Vtext = VInfo(i)%LandText
Return
End If
End Do
Vforce = -1
Write (Vtext, 9130) Wind_Speed
9130 Format ('?? Beaufort_Wind_Scale: incorrect WindSpeed = ', Ipgl4.6)
End Subroutine Beaufort Wind Scale
Subroutine Meta_Data_File()
Implicit None
Integer :: j, kO, kl, uu, fO, fl
Logical :: xok
Character(Len=80) :: tbuf
321
-------
Errors_Detected = .True.
Return
End If
j = Index(pWBAN%Text, ',')
kO = j - 1 ! City
kl = j + 2 ! State
j = Index(pWBAN%TZ, '(') - 1
Write (uu, 9130) 'Station WBAN Number: ', Trim(pWBAN%WBAN)
Write (uu, 9130) 'Station Name: ', pWBAN%Text(1:kO)
Write (uu, 9130) 'Station Location (State): ', pWBAN%Text(kl:kl+1)
Write (uu, 9130) 'Station Time Zone: ', Trim(pWBAN%TZ)
9130 Format (a, a)
Write (uu, 9150) 'Station Latitude: ', &
pWBAN%Lat%Letter, pWBAN%Lat%degrees, pWBAN%Lat%minutes
Write (uu, 9150) 'Station Longitude: ', &
pWBAN%Lon%Letter, pWBAN%Lon%degrees, pWBAN%Lon%minutes
9150 Format (a, al, 15, ' degrees ', 12, ' minutes')
Write (uu, 9130) 'File generation date: ', Trim(TimeStamp)
Call Str_years(tbuf)
Write (uu, '(a,a)') 'Years present:', Trim(tbuf)
fO = 1 + Index (name_met, DirDelim, Back=.True.)
fl = Len_trim(name_met)
Write (uu, 9190) name_met(f0:f1)
9190 Format (&
1!', /, &
'! The meteorological Daily Values File "', a, &
'" has the following format:')
Call lOClose(uu)
End Subroutine Meta Data File
322
-------
! Copy contents of Uin to Uout.
Implicit None
Integer, Intent(In) :: Uin, Uout
iostatus
tbuf
Read a Line: Do
Read (Uin, '(a)', iostat=iostatus) tbuf
If (iostatus /= 0) Exit Read_a_Line
Write (Uout, '(a)') Trim(tbuf)
End Do Read a Line
323
-------
Utils2
!Use Binary Tree
!Use Date_Module
!Use FileStuff
!Use Floating Point Comparisons
!Use GetNumbers
!Use loSubs
!Use Linked_List
!Use Read_Info
!Use Reallocate Module
!Use SAMSON
!Use Strings
!Use UtilsO
!Use Utilsl
Use ETO
Use Fix_Data_Records
Use Global Variables
Use Precipitation module
Use Process_Gaps
Use Stats
Use Utils4
Implicit None
Type :: XQuadrant
Integer :: n = 0
Integer, Dimension(1:24) :: h = 0
Real : : ws_max = Zero
Integer :: ws_d!2 = 0
Real :: ws mean = Zero
End Type XQuadrant
Type :: XPerSeg
Type(XQuadrant), Dimension(:), Pointer :: Quads
Integer, Dimension(:), Pointer :: Tied List
Integer :: Tied_N
Integer :: isub, imod
Real :: rdiv
End Type XPerSeg
Do k = 1, f_SAMSON
Select Case(k)
Case(f_EHR)
Case(f_EDNR)
Case(f GHR)
324
-------
I I I
I I I
I I I
I I I
I I I
I I I
I I I
I I I
I I I
I I I
I I I
I I I
I I I
I I I
I I I
I I I
I I I
I I I
I I I
I I I
I I I
I I I
I I I
Case(f_DNR)
Case(f_DHR)
Case(f_TSC)
Case(f_OSC)
Case(f_DBT)
Case(f_DPT)
Case(f_RH)
Case(f_SP)
Case(f_WD)
Case(f_WS)
Case(f_HV)
Case(f_CH)
Case(f_pH2O)
Case(f baod)
Case(f_SD)
Case(f_DSLS)
Case(f_HP)
Case(f_OI)
Case Default
Write(6,*;
Stop
End Select
End Do
! Direct Normal Radiation
! Diffuse Horizontal Radiation
! Total Sky Cover
! Opaque Sky Cover
! Dry Bulb Temperature
! Dew Point Temperature
! Relative Humidity
! Station Pressure
! Wind Direction
! Wind Speed
! Horizontal Visibility
! Ceiling Height
! Precipitable Water
! Broadband Aerosol Optical Depth
! Snow Depth
! Days since last Snowfall
! Hourly Precipitation (value + flags)
! Observation Indicator / Present Weather
!
Integer :: yyyy, ybase, from beg, from end, to beg, to end
Logical :: okay
Integer :: k, kl, iv, jday, hh25, j, vdim
Integer :: ierr
Type(Val_and_Flag), Dimension(:), Pointer :: HV, CH, HP, DSLS, vf
! Transfer data SAMSON v 1.1 to SAMSON v 1.0
Xparam(f_GHR )%Samson_vlO = Xparam(f_GHR )%Samson_vll
Xparam(f_EHR )%Samson_vlO = Xparam(f_EHR )%Samson_vll
Xparam(f_EDNR)%Samson_vlO = Xparam(f_EDNR)%Samson_vll
Xparam(f_DNR )%Samson_vlO = Xparam(f_DNR )%Samson_vll
325
-------
Xparam(f_DHR )%Samson_vlO = Xparam(f_DHR )%Samson_vll
! Read observed ppt data
Call Read_Hourly_Ppt(okay)
If (.Not. okay) Then
Errors Detected = .True.
End If
Call Read_Daily_Ppt(okay)
If (.Not. okay) Then
Errors_Detected = .True.
End If
!
! See also
! if some years are missing, fill the missing years with data
! from an existing year. We will delete the added data when
! generating the rO files and the MET records.
Do yyyy = MinYear, MaxYear
! If the SAMSON version 1.0 file is missing, then we will
! say that everyting is missing, even if SAMSON version 1.1
! is present.
If (Year_Data(yyyy)%SAMSON_vlO == 0) Then
! For non-leap years, copy year = 1963
! For leap years, copy year = 1964
! These years are always present.
If (IsLeapYear(yyyy)) Then
ybase = 1964
Else
yb ase = 1963
End If
! My previous comment notwithstanding, verify that the
! data for the stated base year exists.
If (Year_Data(ybase)%SAMSON_vlO == 0) Then
Write(*, *) '?? Missing base year ', yyyy, '. Stopping
Write(ULog, *) '?? Missing base year ', yyyy, '. Stopping
Stop '?? Missing base year '
End If
Call ymdh to iv(ybase, mm=01, dd=01, hh=01, iv=from beg)
Call ymdh to iv(ybase, mm=12, dd=31, hh=Nhours, iv=from end)
Call ymdh_to_iv(yyyy, mm=01, dd=01, hh=01, iv=to_beg)
Call ymdh_to_iv(yyyy, mm=12, dd=31, hh=Nhours, iv=to_end)
Write (ULog, 9130) yyyy, ybase, yyyy
Format (/, &
Ix, '## Data for ', 14, &
1 is missing and was filled with data from ', 14, &
326
-------
1 for analysis purposes. ' , /, &
Ix, ' The year ', 14, &
1 will not be present in the final suite of files.'
! Copy all records.
Do k = 1, f_SAMSON
Xparam(k)%Samson vlO(to beg:to end) = &
Xparam(k)%Samson_vlO(from_beg:from_end)
End Do
Obs ppt(to beg:to end) = Obs ppt(from beg:from end)
End If
End Do
Call Standardize ppt(HP, Accum Samson, okay)
If (.Not. okay) Then
Errors_Detected = .True.
End If
If (Have_ppt_Obs_hourly_data) Then
Call Standard!ze_ppt (Obs_ppt, Accum_EI, okay)
If (.Not. okay) Then
Errors Detected = .True.
End If
End If
! Fix SAMSON performs manual correction of the data records.
! Be careful when fixing precipitation records, in particular,
! with runs of missing, deleted, or accumulated values.
! Make sure the beginning or end of such a run is not lost
! on transfer.
Call Fix_SAMSON(okay)
If (.Not. okay) Then
Errors_Detected = .True.
End If
! #1.1 Process BAOD first. For some sites, BAOD is
! missing for nighhtime (in particular, hh = 24).
! #3. Interpolate, fill gaps, and otherwise process the
! hourly data as needed.
I
! #4. THEN fill the 25th hour with the daily values.
! Compute other parameters as needed.
! Note that if we compute the daily values first
! and then fill gaps (invert steps #3 and #4), the
! interpolation scheme would have to determine if
327
-------
! the interpolation boundary or the interpolation
! range contains a "25th" hour. Do not use nor
! fill a 25th hour during interpolation.
! Missing values for BAOD reguire special treatment.
Call Process_BAOD(okay)
If (.Not. okay) Then
Errors_Detected = .True.
End If
!! Use Fill-Gaps algorithms. 12 Mar 2002 2:08 pm
!! Wind Speed: Fill missing hourly values with the monthly mean
MCall Set_Hourly_Values (f_WS, okay)
!!If (.Not. okay) Then
!! Errors_Detected = .True.
!!End If
! Set the 25th hour to NaN. We will see these on output,
! unless set first.
! 6 May 2002 3:24 pm: Warning: according to If95,
! Zero/Zero == NaN
! and Nint(Xparam(1)%Samson_vlO(2*NHours)%v) == 0,
! i.e, Nint(NaN) == 0 !?
Do k = 1, f_SAMSON
vf => Xparam(k)%Samson vlO
vdim = Ubound(vf,l)
vf(NHours:vdim:NHours)%v = Zero / Zero
End Do
Call ToTTy('Process_Set: Process_Days_since_last_Snowfall')
Call FLushAll()
Call Process_Days_since_last_Snowfall(okay)
If (.Not. okay) Then
Errors_Detected = .True.
End If
! Convert "flag" values to numbers
HREF="Onotes.txt#Note_15">
HREF="Onotes.txt#Note_16">
Write(ULog, *)
Write(ULog, 9150) '## Maximum Horizontal Visibility:
Write(ULog, 9150) '## Maximum_Ceiling_Height :
Format (Ix, a, Ipgl4.6)
Maximum Horizontal Visibility
Maximum Ceiling Height
328
-------
! Change the flag value of Unlimited visibility (777.7)
! to 110% of the maximum value attained.
! Change the flag values of Unlimited ceiling height (77777),
! Cirroform (88888) to 110% of the maximum attained.
Fieldlnfo(f_HV)%maximum_value = Maximum_Horizontal_Visibility * 1.10
Fieldlnfo(f CH)%maximum value = Maximum Ceiling Height ^ 1.10
Do iv = 1, kl
If (HV(iv)%s == TJJnlimited) Then
HV(iv)%v = Fieldlnfo(f_HV)%maximum_value
End If
Select Case(CH(iv)%s)
Case(T_Unlimited)
CH(iv)%v = Fieldlnfo(f_CH)%maximum_value
Case(T Cirroform)
CH(iv)%v = Fieldlnfo(f_CH)%maximum_value
End Select
End Do
Call ToTTy('Process_Set: Fill gaps et al.')
Call FLushAll()
! Fill gaps et al.
Do k = 1, f_SAMSON
! Gaps anywhere.
Select Case(k)
Case(f_HP)
! Hourly precipitation must be processed after
! everybody else. We will process HP manually.
!Case(f_HV) ! Horizontal Visibility
!Case(f_CH) ! Ceiling Height
! Also, be careful filling gaps when
! #1. Visibility == 777.7 = unlimited visibility.
! #2. Ceiling Height == 77777 = unlimited ceiling height, or
! == 88888 = Cirroform.
! We have this problem for 14914 (Fargo), 1965-01-01 12h-15h
Case Default
Call Find_Gaps(k)
End Select
End Do
329
-------
! Process precipitation.
Call Process Precipitation(okay)
If (.Not. okay) Then
Errors_Detected = .True.
End If
! Estimate missing Dew Point observations.
Call Fill_Dew_Point(okay)
If (.Not. okay) Then
Errors_Detected = .True.
End If
Call ToTTy('Process_Set: Compute daily values.')
Do k = 1, f_SAMSON
okay = .True.
f_EHR, & ! Extraterrestrial Horizontal Radiation
f_EDNR, & ! Extraterrestrial Direct Normal Radiation
f_GHR, & ! Global Horizontal Radiation
f_DNR, & ! Direct Normal Radiation
f_DHR) ! Diffuse Horizontal Radiation
Call Daily Values(k, T Cumulative, okay)
& ! **•
f_TSC, & ! Total Sky Cover
f_OSC, & ! Opaque Sky Cover
f_DBT, & ! Dry Bulb Temperature
f_DPT, & ! Dew Point Temperature
f_RH, & ! Relative Humidity
f_SP, & ! Station Pressure
f_HV, & ! Horizontal Visibility
f_baod, & ! Broadband Aerosol Optical Depth
f_CH, & ! Ceiling Height
f_pH2O, & ! Precipitable Water
f SD) ! Snow Depth
Case(f_HP) ! Hourly Precipitation (value + flags)
! Do nothing. The daily value was computed by
!
ICall Daily_Values(k, T_Cumulative, okay)
330
-------
Case(f_WD) ! Wind Direction
! Do nothing. Wind Direction and Wind Speed must
! be done simultaneously. See f_WS below.
! = 29 Jan 2002 4:11 pm
! ^ since we are computing mean wind speed,
! we need to fill the 25-th hour with Zero.
Call Daily_Values(k, T_Not_Applicable, okay)
Case(f_WS) ! Wind Speed
i ***
! If a particle were transported ... **** Needs work
i *** Decompose into x- and y- components ...
! The daily value is the vector sum of the hourly components.
! -- not anymore 29 Jan 2002 4:08 pm
! Call Daily_Wind_Speed(okay) ! unused, 29 Jan 2002 4:08 pm
Case(f_DSLS) ! Days since last Snowfall
! Copy the 24h value to the 25h.
hh25 = 0
Do jday = jdO, jdl ! step by day
hh25 = hh25 + NHours ! Index of the 25th hour of jday == (jday-jdO+1)*NHours
j = hh25 - 1 ! the 24th hour of the current day (jday)
Select Case(DSLS(j)%s)
Case(T_Missing, T_Not_Applicable, T_Undefined, T_Perpetual_Darkness)
! Not a valid number.
DSLS(hh25) = Val_and_Flag(T_Estimated, Zero, '')
Case Default
DSLS(hh25) = Val_and_Flag(T_Estimated, DSLS(j)%v, '')
End Select
End Do
okay = .True.
Case(f OI) ! Observation Indicator / Present Weather
Call Daily_Values(k, T_Not_Applicable, okay)
Case Default
Write(6,*) 'Process_Set: Internal error: no "CASE(f_k)" for k=', k
Stop '?? Stopping in Process_Set: Internal error: no "CASE(f_x)"'
End Select
If (.Not. okay) Then
Errors Detected = .True.
End If
Call FLushAll()
End Do
331
-------
! Rs (a + b[opaque_sky_cover]) * Ra
I =
! Rso a * Ra
I
! (a + b[opaque_sky_cover])
Call ToTTy('ETO_et_al'
Call ETO_et_al(okay)
Call FLushAll()
! Must be the last statement of the file.
Call ToTTy('Process_Set: Set_Missing_Flags')
Call Set_Missing_Flags()
End Subroutine Process Set
Implicit None
Integer, Intent(In) :: k_id
Integer, Pointer :: min obs per day
Type(Val_and_Flag), Dimension(:), Pointer :: xsd ! X-Data
! Gaps anywhere.
Select Case(k_id)
Case(f_OI)
! Observation Indicator, Present weather
! It makes no sense to fill gaps in these parameters
!Case(f_HV) ! Horizontal Visibility
!Case(f CH) ! Ceiling Height
! Also, be careful filling gaps when
! #1. Visibility == 777.7 = unlimited visibility.
! #2. Ceiling Height == 77777 = unlimited ceiling height, or
332
-------
! == 88888 = cirroform.
! We have this problem for 14914 (Fargo), 1965 1 1 12h-15h
Case Default
Call Find_Gaps(k_id)
End Select
End Subroutine Process Param
!
! Fill missing hourly values with the monthly mean
Implicit None
Integer, Intent(In) :: k_id
Logical, Intent(Out) :: Xok
Integer :: yyyy, mm, dd, hh, iv, jvO, jvl
Integer :: ierr, npoints, minN
Real :: xsum
Type(Val_and_Flag), Dimension(:), Pointer :: vf, WD
Integer, Dimension(:), Pointer :: Days_in_Month
! Wind Direction
CheckThisMonth: Do iv = JvO, jvl
If (Modulo(iv,25) == 0) Then
! Skip 25-th hour of the day.
Cycle CheckThisMonth
End If
Select Case(vf(iv)%s)
Case(T_Missing, T_Not_Applicable, T_Undefined, T_Perpetual_Darkness)
333
-------
! Not a valid number: Do nothing.
Case Default
npoints = npoints + 1
xsum = xsum + vf(iv)%v
End Select
End Do CheckThisMonth
! Done with the month. Do we have points?
If (npoints < minN) Then
ierr = ierr + 1
Write(ULog, 9130) Trim(Fieldlnfo(k_id)%Name), yyyy, mm, npoints, minN
Format(Ix, '?? Set_Hourly_Values: ', a, Ix, &
14, '-', 12.2, ': found ', 10, ' points, need ', 10)
Cycle MonthLoop
End If
! Set Values
SetThisMonth: Do iv = jvO, jvl
If (Modulo(iv,25) == 0) Then
! Skip 25-th hour of the day.
Cycle SetThisMonth
End If
Select Case(vf(iv)%s)
Case(T_Missing, T_Not_Applicable, T_Undefined, T_Perpetual_Darkness)
vf(iv)%v = xsum
vf(iv)%s = T_Estimated
vf(iv)%f = ''
! If we are dealing with Wind Speed,
! Then zero Wind Direction also.
If (k_id == f_WS) Then
WD(iv)%v = Zero
WD(iv)%s = T_Estimated
WD(iv)%f = ''
End If
End Select
End Do SetThisMonth
End Do MonthLoop
End Do
Xok = (ierr == 0)
End Subroutine Set Hourly Values
Subroutine Daily_Wind_Speed(Xok)
334
-------
! = 29 Jan 2002 4:05 pm
! * module not used. Wind Speed daily value will be
! a simple average.
! Compute Mean Wind Speed over 24-hour periods.
!*
Implicit None
Logical, Intent(Out) :: Xok
Integer :: jday, hhOl, hh24, hh25, hh
Integer :: minN, nn, ierr, ntt, jyyyy, jmm, jdd
Real :: xsum, ysum, resultant, Azimuth, Theta
Type(Val_and_Flag), Dimension!:), Pointer :: WD
Type(Val_and_Flag), Dimension(:), Pointer :: WS
Character(Len(Fieldlnfo(1)%Name)), Pointer :: txt
Type(Stat_Block) :: daily_ws
Call Stat_Initialize(daily_ws, 'Wind Speed, daily values')
ierr = 0
WD => Xparam(f_WD)%Samson_vlO ! Wind Direction in degrees
WS => Xparam(f_WS)%Samson_vlO ! Wind Speed in m/s
minN = Fieldlnfo(f WD)%Minimum obs per day
txt => Fieldlnfo(f_WD)%Name
ntt = Len_trim(txt)
! Wind Direction 0-360 Wind direction in degrees.
! (N = 0 or 360, E = 90, S = 180,
! W = 270)
I
! Wind Speed
I
! 0
! Azimuth: Wind direction in degrees.
! (N = 0 or 360, E = 90, S = 180, W = 270)
I
! Theta: regular angle measured counterclockwise
335
-------
! From Mathematica, for the test case: (ArcTan(x/y) + Pi) == Azimuth
! Note that Theta = ArcTan(y/x) + Pi
I
! Subroutine Test Azimuth tests these formulae.
Do jday = jdO, jdl ! step by day
hhOl = (jday-jdO)*Nhours + 1 ! First hour of the day
hh24 = hhOl + 23 ! Last hour of the day (24th)
hh25 = hh24 + 1 ! 25th hour of jday == (jday-jdO+1)*NHours
OneDay: Do hh = hhOl, hh24
Select Case(WD(hh)%s)
Case(T_Missing, T_Not_Applicable, T_Undefined, T_Perpetual_Darkness)
Cycle OneDay
End Select
Select Case(WS(hh)%s)
Case(T_Missing, T_Not_Applicable, T_Undefined, T_Perpetual_Darkness)
Cycle OneDay
End Select
nn = nn + 1
! Convert Azimuth (measured in degrees) to
! Theta (measured in radians)
Theta = Modulo(90-WD(hh)%v, 360) * Degrees_to_Radians
xsum = xsum + WS(hh)%v * Cos(Theta)
ysum = ysum + WS(hh)%v * Sin(Theta)
End Do OneDay
336
-------
WS(hh25)%v = resultant
WS(hh25)%s = T_Estimated
WS(hh25)%f = ''
Call Stat_Add_Point(daily_ws, resultant)
Else
ierr = ierr + 1
WD(hh25)%v = Missing_Data
WD(hh25)%s = T_Missing
WD(hh25)%f = ''
WS(hh25)%v = Missing_Data
WS(hh25)%s = T_Missing
WS(hh25)%f = ''
Call Jd to ymd(jday, j yyyy, jmm, jdd)
Write(ULog, 9130) jyyyy, jmm, jdd, txt(1:ntt), nn, minN
Format(lx, '?? Daily Values: ', 14, 2('-',12.2), Ix, a, ' found ', 10, ', need ', 10)
End If
End Subroutine Daily Wind Speed
Subroutine Process Days since last Snowfall(Xok)
!*
! Correct the field Days since last Snowfall.
! Observation Indicator (OI)
! | Present_weather (PW)
! | | Days Since Last Snowfall (DSLS)
! yy mm dd hh | |
, __ __ __ __ _ 45
61
62
63
63 ! hours 2-23
63
... 64 ! Counter increments
337
-------
! Days Since Last Snowfall counter increments on hour == 1.
SAMSON information. For a description of the SAMSON file format see
[Isr] Thu 19 Oct 2000 11:58:19
This field appears to be
column 96 of the input file.
The manual provides no data
to support this assertion.
[Isr] Wed Nov 21 09:33:59 2001
Observation Indicator appears
only in SAMSON v 1.0 files.
0 = Weather observation made.
9 = Weather observation not
made or missing.
If this field = 9 OR if field
13 (wind speed) = missing
(9999. or 99.0), then
fields 6, 7, 8, 10, 11, 17,
and 18 were all modeled and
not actually observed.
Snow-related fields of Present_weather are 4 and 5.
Field Number
Contents
I Values
I | Description
(4)
0 = Light snow
1 = Moderate snow
2 = Heavy snow
3 = Light snow pellets
4 = Moderate snow pellets
5 = Heavy snow pellets
6 = Light ice crystals
7 = Moderate ice crystals
8 = Heavy ice crystals
9 = None if Observation
Indicator element eguals
0, else unknown or
missing if Observation
338
-------
Indicator element
equals 9.
Notes:
Beginning in April 1963, any
occurrence of ice crystals
is recorded as a 7.
Implicit None
Logical, Intent(Out) :: Xok
! Relevant field numbers
Integer, Parameter :: f4 = 4
Integer, Parameter :: f5 = 5
0 = Light snow showers
1 = Moderate snow showers
2 = Heavy snow showers
3 = Light snow sguall
4 = Moderate snow sguall
5 = Heavy snow sguall
9 = None if Observation
Indicator element eguals 0,
else unknown or
missing if Observation
Indicator element eguals 9.
Integer
Integer
Integer
Logical
Integer
Real
Logical
Integer
Integer
jday, hhOl, hh24, hh25, hh
ierr, jyyyy, jmm, jdd, jhh
Observation Indicator
DSLS_unknown, observation_is_missing
nDSLS, hours_with_snow
old_DSLS
and_Flag), Dimension(:), Pointer :: OI, DSLS
339
-------
! Start with Days Since Last Snowfall (DSLS) == Missing,
! because we do not know when the last snow occurred.
I
! This is what SAMSON did for San Juan, P.R. (WBAN=11641),
! where it has not snow in recorded history.
nDSLS = -1 ! Non sensical value.
Do jday = jdO, jdl ! step by day
hhOl = (jday-jdO)*Nhours + 1 ! First hour of the day
hh24 = hhOl + 23 ! Last hour of the day (24th)
hh25 = hh24 + 1 ! 25th hour of jday == (jday-jdO+1)*NHours
!Print_now = (hhOl == jvO)
!If (print_now) Then
! Write(6,*) 'in dsls'
! Write(6,*) 'in dsls'
!End If
! If at any time during the day DSLS(hh)%v holds
! a non-missing value, we want pDSLS to point to it.
! We will use this value for the 25th hour (the daily value).
pDSLS = Tbogus
OneDay: Do hh = hhOl, hh24
! Save the old value to test cognitive dissonances.
old_DSLS = DSLS(hh)%v
If (DSLS(hh)%s /= T_Missing) Then
nDSLS = Nint(old_DSLS)
pDSLS = hh
DSLS_unknown = .False.
Else
DSLS_unknown = .True.
End If
340
-------
!If (Print_now) Then
! Write (6, *) 'in Process Days since last Snowfall1
!End If
Observation_Indicator = Nint(OI(hh)%v) ! 0 or 9
Select Case(Observation Indicator)
Case(9) ! Weather observation not made or missing.
observation_is_missing = .True.
Case(O) ! Weather observation made.
! It snow if
! a) (Present_weather, field #4 /= 9) ! snow
! OR
! b) (Present weather, field #5 /= 9) ! snow showers
observation is missing = .False.
If ((OI(hh)%f(f4:f4) /= '9') .Or. (OI(hh)%f(f5:f5) /= '9')
hours_with_snow = hours_with_snow + 1
!Call iv to ymdh(hh, jyyyy, jmm, jdd, jhh)
IWrite (ULog, 9310) hh, jyyyy, jmm, jdd, jhh, OI(hh)%f(f4:f5)
Format ('## Observation Indicator ==0: ', &
1 for ', 17, Ix, 14, '-', 12.2, '-', 12.2, 13, 'h', &
'; OI(hh)%f(f4:f5) == "', a, '"')
End If
Case Default
Call iv to ymdh(hh, jyyyy, jmm, jdd, jhh)
Write (ULog, 9150) Observation Indicator, hh, jyyyy, jmm, jdd, jhh
Format ('?? Observation Indicator not 0 nor 9: ', 10, &
1 for ', 17, Ix, 14, '-', 12.2, '-', 12.2, 13, 'h')
observation is missing = .True.
End Select
! If Days since last Snowfall is unknown, try to do something sensible.
If (DSLS_unknown) Then
! If it snow during this fraction of a day, then Days since last Snowfall == 0.
! This inference is valid regardless of the value of observation_is_missing.
If (hours_with_snow > 0) Then
nDSLS = 0
DSLS(hh)%v = nDSLS
DSLS(hh)%s = T_Estimated
DSLS(hh)%f = ''
pDSLS = hh
DSLS_unknown = .False.
Else If (observation_is_missing) Then
! Nothing we can do. It may have snow (or not).
nDSLS = -1 ! Non sensical value.
Else
! It has not snow.
If (nDSLS >= 0) Then
! We have an estimate.
341
-------
DSLS(hh)%v = nDSLS
DSLS(hh)%s = T_Estimated
DSLS(hh)%f = ''
pDSLS = hh
DSLS_unknown = .False.
Else
! We still do not know when the last snowfall occurred.
DSLS(hh)%v = Missing_Data
DSLS(hh)%s = T_Missing
DSLS(hh)%f = ''
End If
End If
Else
! Old value not unknown. Keep it. Nothing to do.
End If
!If (print_now) Then
! Write (ULog, *) str_iv_to_ymdh(hh) , DSLS(hh)
!End If
End Do OneDay
! Now assign a daily value. If we had a non-missing value
! during the day, use it.
If (pDSLS /= Tbogus) Then
! We had a good hour. Use it for the daily value
DSLS(hh25)%v = DSLS(pDSLS)%v
DSLS(hh25)%s = T_Estimated
DSLS(hh25)%f = ''
Else
! Else no good hours during this day. Missing.
DSLS(hh25)%v = Missing_Data
DSLS(hh25)%s = T_Missing
DSLS(hh25)%f = ''
End If
!If (print_now) Then
! Write (ULog, *) str_iv_to_ymdh(hh25), DSLS(hh25)
!End If
! If the value of nDSLS is known, increment it.
! The maximum value of this field is 88 days.
! See SAMSON description above.
If (nDSLS >= 0) nDSLS = Min(nDSLS + 1, 88)
342
-------
! SAMSON. BAOD is not present during nighttime. One value (measured or
! estimated) was produced and replicated for the day. See fragment
! file below (13893_61.txt):
! yy mm dd hh BAOD
! 61 1 1 1 99999.
! 61 1 1 : 99999. ! hours 2-6
! 61 1 1 7 99999.
! 61 1 1 8 .034
! 61 1 1 : .034 ! hours 9-16
! 61 1 2
! 61 1 2
I
! Algorithm:
! For every day: Starting at 1 h, find the first non-missing
! value, call it "rv". Replace all missing values during that
! day with "rv".
Implicit None
Logical, Intent(Out) :: Xok
Integer :: jday, hhOl, hh24, hh
Integer :: ierr, jyyyy, jmm, jdd
Real : : rv
Logical : : xfound
Type(Val_and_Flag), Dimension!:), Pointer :: BAOD
Do jday = jdO, jdl ! step by day
hhOl = (jday-jdO)*Nhours + 1 ! First hour of the day
hh24 = hhOl + 23 ! Last hour of the day (24th)
xfound = .False.
FindRv: Do hh = hhOl, hh24
Select Case(BAOD(hh)%s)
Case(T_Missing, T_Not_Applicable, T_Undefined, T_Perpetual_Darkness)
! Do nothing.
Case Default
rv = BAOD(hh)%v
xfound = .True.
Exit FindRv
End Select
343
-------
End Do FindRv
If (.Not. xfound) Then
ierr = ierr + 1
Call Jd to ymd(jday, j yyyy, jmm, jdd)
Write(ULog, 9130) jyyyy, jmm, jdd
Format(lx, '?? Process_BAOD: Values all missing for ', 14, 2('-',12
Cycle
End If
OneDay: Do hh = hhOl, hh24
If (BAOD(hh)%s == T_Missing) Then
BAOD(hh)%v = rv
BAOD(hh)%s = T_Estimated
BAOD(hh)%f = ''
End If
End Do OneDay
End Do
Xok = (ierr == 0)
Subroutine Daylight_Prevailing_Wind(Xok, MET_pwd, MET_pws, OnlyDaylight)
! Compute Daylight Prevailing Wind
!
! Xok — truth of "all entries in MET ^ not missing"
! MET_pwd -- prevailing wind direction
! MET_pws -- prevailing wind speed
! OnlyDaylight — Truth of "compute prevailing stuff using only daylight hours'1
! The direction of wind is measured in terms of where the air
! is coming from. A northerly wind blows air from north to south.
! A southwesterly wind blows air from the southwest to the northeast.
I
! The prevailing wind is the wind that blows most freguently across a
! particularly region. Different regions on Earth have different
! prevailing wind directions which are dependent upon the nature of
! the general circulation of the atmosphere and the latitudinal wind
! zones.
I
! Prevailing wind - The wind direction most freguently observed during a
! given period. The periods most freguently used are the observational day,
! month, season, and year. Methods of determination vary from a single count
344
-------
! of periodic observations to the computation of a wind rose.
Implicit None
Logical, Intent(Out) :: Xok
Type(Val_and_Flag), Dimension(jdO:jdl), Intent(Out) :: MET_pwd, MET_pws
Logical, Intent(In) :: OnlyDaylight
Integer
Integer
Integer
Integer
Integer
Logical
Real
i, j, n, nO, JO, jl, ierr, jb, points_in_day
ig, idim, jt, ik, nmax, ndel, min d!2
jday, hhOl, hh24, hh
JyyyYf jmm, jdd, jhh
prevailing_guadrant, Prevailing_CoordSys, ncount
is daytime, N is even
xws, xwd, median wd, mean ws, vmax
Integer, Parameter :: CoordAx dim = 6
Type(XPerSeg), Dimension(CoordAx dim), Target, Save :: CoordAx ! Coordinate Axes
Type(XQuadrant), Pointer :: qik
Logical :: have_quadrant
Integer, Dimension(:), Pointer :: p2Hour
Integer, Pointer :: p2N
Logical, Save :: first_time = .True.
Logical :: Print now = .False.
Integer :: jul dayO, jul dayl
ierr = 0
WD => Xparam(f_WD)%Samson_vlO
WS => Xparam(f_WS)%Samson_vlO
Rs => Xparam(f_Rs)%Samson_vlO
! Wind Direction in degrees
! Wind Speed in m/s
! Global Horizontal Radiation
MET_pwd = Val_and_Flag(T_Missing, Missing_Data, '')
MET_pws = Val_and_Flag(T_Missing, Missing_Data, '')
! Allocate only once.
first_time = .False.
! Allocate Quadrant/HemiPlane arrays.
! isub, rdiv, imod define the function
! i = Modulo(Floor((x-isub)/rdiv),imod) + 1
! imod = number of guadrants
! rdiv * imod = 360
345
-------
! Define 2 sets of coordinate axes (qOO, q4 5 ) , each with 4 quadrants
! qOO — standard coordinate axes.
! i = Modulo(Floor((x)/90.0),4) + 1
iq = 1
idim = 4
Allocate(CoordAx(iq)%Quads(idim) , CoordAx(iq)%Tied List(idim) )
CoordAx(iq)%isub = 0
CoordAx(iq)%imod = idim
CoordAx(iq)%rdiv = 360.0 / idim
! q45 -- axes rotated 45°
! i = Modulo(Floor((x-45)/90.0),4) + 1
iq = iq + 1
idim = 4
Allocate(CoordAx(iq)%Quads(idim), CoordAx(iq)%Tied_List(idim))
CoordAx(iq)%isub = 45
CoordAx(iq)%imod = idim
CoordAx(iq)%rdiv = 360.0 / idim
! hemiplanes: 2 "quadrants"
Do i = 1, 4
iq = iq + 1
idim = 2
Allocate(CoordAx(iq)%Quads(idim), CoordAx(iq)%Tied_List(idim))
jb = (i-1) * 45 ! jb == 0, 45, 90, 135
!i = Modulo(Floor((x-jb)/ISO.0),2) + 1
CoordAx(iq)%isub = jb
CoordAx(iq)%imod = idim
CoordAx(iq)%rdiv = 360.0 / idim
End Do
!! http://www.windpower.dk/tour/wres/rose.htm
!! Wind rose. Divide the compass into 12 sectors (European Wind
!! Atlas standard), one for each 30 degrees of the horizon.
!! A wind rose may also be drawn for 8 or 16 sectors.
!iq = iq + 1
!idim = 12
!Allocate(CoordAx(iq)%Quads(idim), CoordAx(iq)%Tied_List(idim))
ICoordAx(iq)%isub = 0
!CoordAx(iq)%imod = idim
!CoordAx(iq)%rdiv = 360.0 / idim
End If
ByJD: Do jday = jdO, jdl ! step by day
hhOl = (jday-jdO)*Nhours + 1 ! First hour of the day
hh24 = hhOl + 23 ! Last hour of the day (24th)
! #1. Find the quadrant with the most observations.
346
-------
If tie: rotate axes 4J
and try again.
q45
(4)
(1)
(3)
(2)
22
13E
The formulae below make sure quadrants are wrapped around properly,
so that the resulting quadrant is always an integer in [1, 4] .
See Mathematica notebook
If there are still ties, cut into hemiplanes:
\
\
347
-------
:= Mod[Floor[(x)/ISO.0],2]+1
:= Mod[Floor[(x-45)/ISO.0],2]+1
:= Mod[Floor[(x-90)/ISO.0],2]+1
:= Mod[Floor[(x-135)/ISO.0],2]+1
! XQuadrant(n, h, max, d!2, mean)
CoordAx(iq)%Quads = XQuadrant(0, 0, Zero, Huge(0), Zero)
CoordAx(iq)%Tied_List = 0
CoordAx(iq)%Tied_N = 0
End Do
points in day = 0
Select Case(Rs(hh)%s)
Case(T_Missing, T_Not_Applicable, T_Undefined, T_Perpetual_Darkness)
Cycle OneDay
End Select
If (OnlyDaylight) Then
! Skip non-daylight hours.
is_daytime = (Rs(hh)%v > Zero)
If (.Not. is_daytime) Then
Cycle OneDay
End If
End If
Select Case(WD(hh)%s)
Case(T_Missing, T_Not_Applicable, T_Undefined, T_Perpetual_Darkness)
Cycle OneDay
End Select
348
-------
Select Case(WS(hh)%s)
Case(T_Missing, T_Not_Applicable, T_Undefined, T_Perpetual_Darkness)
Cycle OneDay
End Select
points in day = points in day + 1
xwd = WD(hh)%v
xws = WS(hh)%v
Do iq = 1, CoordAx dim
! (1): ik = Modulo(Floor((xwd-0)/90.0),4) + 1
! (2): ik = Modulo(Floor((xwd-45)/90.0),4) + 1
! (3-6): jb = (j-1) * 45, j = 1..4, jb == 0, 45, 90, 135
! ik = Modulo(Floor( (xwd-jb)/ISO.0) ,2) + 1
! (7): ik = Modulo(Floor((xwd-0)/30.0),12) + 1
n = qik%n + 1
qik%n = n
qik%h(n) = hh
qik%ws_mean = qik%ws_mean + xws
! ndel — the absolute delta time from 12 noon
ndel = Abs(hh - hhOl +1-12)
If ((xws-qik%ws_max) > EpsO) Then
! xws > maximum so far : Found a new maximum.
qik%ws max = xws
qik%ws_d!2 = ndel
Else If (Abs(xws-qik%ws_max) < EpsO) Then
! xws == maximum wind speed
! We have a wind speed equal to the recorded maximum WS.
! If this wind speed falls in the same quadrant of the
! recorded maximum, ok. Else we cannot determine the
! most frequent quadrant uniquely.
If (ndel < qik%ws_d!2) Then
qik%ws_d!2 = ndel
End If
End If
End Do
End Do OneDay
If (points in day == 0) Then
! No data.
ierr = ierr + 1
Cycle ByJD
End If
349
-------
! Compute mean wind speeds for each quadrant.
Do iq = 1, CoordAx dim
Do ik = 1, CoordAx(iq)%imod
n = CoordAx(iq)%Quads(ik)%n
If (n == 0) Cycle
CoordAx(iq)%Quads(ik)%ws mean = CoordAx(iq)%Quads(ik)%ws mean / n
End Do
End Do
have quadrant = .False.
Do iq = 1, CoordAx dim
! prevailinq_quadrant -- the quadrant with the most observations.
! This is the prevailinq wind quadrant.
! CoordAx(iq)%Quads(prevailinq quadrant)%n —
! The number of elements the fell
! in prevailinq_quadrant.
! ncount > 1 ==> ties.
prevailinq quadrant = Maxloc(CoordAx(iq)%Quads(:)%n, Dim=l)
nmax = CoordAx(iq)%Quads(prevailinq quadrant)%n
! Tied_List is a list of the tied quadrants.
ncount = 0
Do ik = 1, CoordAx(iq)%imod
If (CoordAx(iq)%Quads(ik)%n == nmax) Then
ncount = ncount + 1
CoordAx(iq)%Tied List(ncount) = ik
End If
End Do
CoordAx(iq)%Tied_N = ncount
If (ncount == 1) Then
have_quadrant = .True.
Prevailinq_CoordSys = iq
ik = prevailinq quadrant
p2N => CoordAx(iq)%Quads(ik)%n
p2Hour => CoordAx(iq)%Quads(ik)%h
Selection_Criterion = 'quadrant with the most observations'
Exit
End If
End Do
If (.Not. have_quadrant) Then
! Amonq the tied quadrants, find the one with the larqest mean wind speed.
L41: Do iq = 1, CoordAx_dim
ncount = 0
vmax = -Huqe(Zero)
350
-------
If (CoordAx(iq)%Quads(ik)%ws mean > vmax) Then
ncount = 1
vmax = CoordAx(iq)%Quads(ik)%ws_mean
prevailing_quadrant = ik
Else If (Abs(CoordAx(iq)%Quads(ik)%ws_mean - vmax) < EpsO) Then
! CoordAx(iq)%Quads(ik)%ws mean == vmax
! A tie. Note that we cannot qive up yet because vmax
! may still change.
ncount = ncount + 1
End If
End Do L42
If (ncount == 1) Then
have quadrant = .True.
Prevailing CoordSys = iq
ik = prevailing_quadrant
p2N => CoordAx(iq)%Quads(ik)%n
p2Hour => CoordAx(iq)%Quads(ik)%h
Selection Criterion = 'tied quadrant with the largest mean wind speed1
Exit
End If
End Do L41
End If
If (.Not. have_quadrant) Then
! Among the tied quadrants, select the quadrant where the
! maximum wind speed is closest to noon.
L51: Do iq = 1, CoordAx dim
ncount = 0
min_d!2 = Huge(0)
L52: Do jt = 1, CoordAx(iq)%Tied_N
ik = CoordAx(iq)%Tied_List(jt)
If (CoordAx(iq)%Quads(ik)%ws_d!2 < min_d!2) Then
ncount = 1
min_dl2 = CoordAx(iq)%Quads(ik)%ws_d!2
prevailing_quadrant = ik
Else If (CoordAx(iq)%Quads(ik)%ws_d!2 == min_d!2) Then
! A tie. Note that we cannot give up yet because min d!2
! may still change.
ncount = ncount + 1
End If
End Do L52
If (ncount == 1) Then
have_quadrant = .True.
Prevailing CoordSys = iq
351
-------
ik = prevailing_quadrant
p2N => CoordAx(iq)%Quads(ik)%n
p2Hour => CoordAx(iq)%Quads(ik)%h
Selection_Criterion = 'tied quadrant with the maximum wind speed closest to noon'
Exit
End If
End Do L51
End If
! 25 Apr 2002 11:21 am: Debuqqinq test.
! For 26451: Anchoraqe, AK,
! hhOl == 66701, 1968-04-22 Ih
! Selection Criterion: tied quadrant with the larqest mean wind speed
!If (.Not. have_quadrant) Then
If ((.Not. have_quadrant) .Or. (Print_now)) Then
Call iv_to_ymdh(hhOl, jyyyy, jmm, jdd, jhh)
If (Print_now) Then
Write(ULoq, 9130) &
'?? Dayliqht Prevailing Wind: "Print now" at ', &
hhOl, jyyyy, jmm, jdd, jhh, 'h'
Else
Write(ULog, 9130) &
'?? Daylight Prevailing Wind: No resolution at ', &
hhOl, jyyyy, jmm, jdd, jhh, 'h'
End If
Format (/, Ix, a, 17, Ix, 14, 2('-',i2.2), 13, a)
Qerr: Do hh = hhOl, hh24
Select Case(Rs(hh)%s)
Case(T_Missing, T_Not_Applicable, T_Undefined, T_Perpetual_Darkness)
Cycle Qerr
End Select
If (OnlyDaylight) Then
! Skip non-daylight hours.
is_daytime = (Rs(hh)%v > Zero)
If (.Not. is_daytime) Then
Cycle Qerr
End If
End If
Select Case(WD(hh)%s)
Case(T_Missing, T_Not_Applicable, T_Undefined, T_Perpetual_Darkness)
Cycle Qerr
End Select
Select Case(WS(hh)%s)
352
-------
Case(T_Missing, T_Not_Applicable, T_Undefined, T_Perpetual_Darkness)
Cycle Qerr
End Select
i = hh - hhOl + 1
Write(ULog, 9150) &
', hh, jyyyy,
Format (Ix, a, 17, Ix,
End Do Qerr
Write (ULog, *) ' have guadrant: ', have guadrant
If (have_guadrant) Then
Write (ULog, *) 'Selection_Criterion
Write (ULog, *) 'Prevailing_CoordSys
Write (ULog, ^) 'prevailing guadrant
Write (ULog, *) ' p2N
nO = Min(p2N, Ubound(CoordAx(1)%Quads(1)%h,1), 24)
Write (ULog, *) ' p2Hour: ', p2Hour(l:nO)
End If
Do ig = 1, CoordAx_dim
Write(ULog, *)
Write(ULog, *) '================================'
Write(ULog, *) 'CoordAx(', ig, ')'
n = CoordAx(ig)%Tied_N
nO = Min(n, 4)
Write(ULog, *) ' %Tied_List:'
Do ik = 1, CoordAx(ig)%imod
n = CoordAx(ig)%Quads(ik)%n
nO = Min(n, Ubound(CoordAx(1)%Quads(1)%h,1
Write(ULog, *)
' ) ' %Quads( ', ik,
<) ' h: '
Trim(Selection_Criterion)
Prevailing CoordSys
prevailing guadrant
p2N
Write(ULog,
Write(ULog,
Write(ULog,
Write(ULog,
Write(ULog,
End Do
End Do
Call FLushAll()
End If
%ws_d!2:
%ws mean:
): n == ', n
CoordAx(ig)%Quads(ik)%h(l:nO)
CoordAx(ig)%Quads(ik)%ws_max
CoordAx(ig)%Quads(ik)%ws_d!2
CoordAx(ig)%Quads(ik)%ws mean
If (.Not. have_guadrant)
ierr = ierr + 1
Cycle ByJD
End If
! The prevailing wind direction is the median of the
! wind directions that fell in the most freguent guadrant.
353
-------
! (at most 24 entries).
! is problematic: The entries of WD (hh01:hh24) are (say)
! 137058:137067. On entry to Insertion_Sort_Real the addresses
! of WD are 1:24. However, the values in pindex are still in
! the range 137058:137067. So, references to WD(p2Hour(j)) will
! be out of bounds *within* the subroutine. We need to pass the
! address bounds of WD to the subroutine. This is more trouble
! than it is worth, since the insertion sort is used only once
! and the code is short. Therefore ( 2 Jul 2002 1:29 pm) code in line.
! Insertion sort is best if the number of entries is less than
! approximately 20 entries.
Do i = 1, p2N - 1
Do j = i+1, p2N
If (WD(p2Hour(j))%v < WD(p2Hour(i))%v) Then
nO = p2Hour(i)
p2Hour(i) = p2Hour(j)
p2Hour(j) = nO
End If
End Do
End Do
If (Print_now) Then
Write(ULog, *)
Write(ULog, *) '================================'
nO = Min(p2N, Ubound(CoordAx(1)%Quads(1)%h, 1) , 24)
Write(ULog, *) 'Sorted p2Hour: ', p2Hour(l:nO)
Write(ULog, *) 'Sorted WD(p2Hour(1:nO))%v: ', WD(p2Hour(1:nO))%v
Write(ULog, *)
End If
! Compute median
N_is_even = (Modulo(p2N,2) == 0)
If (N_is_even) Then
i = p2N / 2
j 0 = p2Hour(i)
jl = p2Hour(i+l)
median_wd = (WD(jO)%v + WD(jl)%v) / 2
If (Print_now) Then
Write(ULog, *) 'N is even: i, JO, jl: ', i, JO, jl
End If
Else
i = (p2N + 1) / 2
j 0 = p2Hour(i)
median_wd = WD(jO)%v
354
-------
If (Print_now) Then
Write(ULog, *) 'N is odd: i, JO: ', i, JO
End If
End If
!! The prevailing wind speed is the mean of the
!! wind speeds that fell in the prevailing guadrant.
!mean_ws = Zero
!Do i = 1, p2N
! j 0 = p2Hour(i)
! mean ws = mean ws + WS(jO)%v
!End Do
!mean_ws = mean_ws / p2N
mean ws = CoordAx(Prevailing CoordSys)%Quads(prevailing guadrant)%ws mean
If (Print_now) Then
Write(ULog, *) 'median_wd: ', median_wd
Write(ULog, *) 'mean ws..: ', mean ws
End If
MET_pwd(jday) = Val_and_Flag(T_Estimated, median_wd, '')
MET_pws(jday) = Val_and_Flag(T_Estimated, mean_ws, '')
End Subroutine Daylight Prevailing Wind
Function a eg b(A, B)
Implicit None
Type(XQuadrant), Dimension(:), Intent(In) :: A,
Logical :: a eg b
If (sa /= sb) Return
Do i = 1, sa
If (A(i)%n /= B(i)%n) Then
Return
End If
355
-------
Return
End If
End Do
a_eq_b = .True.
End Function a eq b
356
-------
UtilsS
! ! ! 1
I I I 2
! ! ! 3
MI 4
III 5
! ! ! 6
! ! ! 7
! ! ! 8
I I I 9
! ! ! 10
! ! ! 11
! ! ! 12
III 13
! ! ! 14
! ! ! 15
! ! ! 16
I I I 17
! ! ! 18
! ! ! 19
! ! ! 20
! ! ! 21
I I I 22
! ! ! 23
! ! ! 24
III 25
16- 22
24- 30
32- 40
42- 50
52- 60
62- 65
67- 70
72- 78
80- 86
88- 92
94-100
102-106
108-114
116-123
125-132
134-136
138-148
150-154
156-163
165-170
172-176
178-186
188-197
199-208
210-
Print Cols
Extraterrestrial Horizontal Radiation
Extraterrestrial Direct Normal Radiation
Global Horizontal Radiation
Direct Normal Radiation
Diffuse Horizontal Radiation
Total Sky_Cover
Opaque Sky Cover
Dry Bulb Temperature
Dew_Point Temperature
Relative Humidity
Station Pressure
Wind Direction
Wind Speed
Visibility
Ceiling Height
Observation Indicator
Present_weather
Precipitable Water
Broadband Aerosol Optical_Depth
Snow_Depth
Days Since Last Snowfall
Hourly Precipitation
FAO Short Grass PET
K-P FWS Evaporation
Pan Evaporation
!Use Binary_Tree
!Use FileStuff
!Use GetNumbers
!Use Linked_List
!Use Read_Info
!Use SAMSON
Use Strings
!Use UtilsO
!Use Utilsl
!Use Date_Module
Use Global_Variables
!Use loSubs
!Use Utils2
!Use Winteracter
Implicit None
Contains
357
-------
Subroutine Print_Cols(Col_Text, kO, kl, Last, Icount, Initialize)
! Call Print_Cols('Extraterrestrial Horizontal Radiation', kO, kl, Icount=l)
Implicit None
Character(Len=*), Intent(In) :: Col_Text
Integer, Intent(In) :: kO, kl
! The value of "Last" is unimportant.
! We check only for the presence/absence of the variable.
Logical, Optional, Intent(In) :: Last
Integer, Optional, Intent(In) :: Icount
Logical, Optional, Intent(In) :: Initialize
If (Present(Initialize)) Then
! By increasing 'uu' we will be able to store
! different seguences in different files.
uu = uu + 1
do return = .False.
Return
End If
If (do_return) Then
Return
End If
If (Present(Icount)) Then
Kcount = Icount - 1
End If
If (Present(Last)) Then
do_return = .True.
Return
End If
End Subroutine Print Cols
358
-------
Implicit None
Integer, Intent(In) :: kO
Character(Len=*), Intent(In) :: xFMT
! The value of "Last" is unimportant.
! We check only for the presence/absence of the variable.
Logical, Optional, Intent(In) :: Last
Integer, Save :: ip = 1
Logical, Save :: do return = .False.
Integer : : j 0, j1
If (do_return) Then
Return
End If
If (Present(Last)) Then
do return = .True.
Write (ULog, '(/, 2a, /)') '## Daily values file format: ', Trim(FMT_dvf)
Return
! Stop '## Stoping in "dvf_MkFMT" as reguested'
End If
JO = 2
jl = Len_trim(xFMT) - 1
If (ip > 1) Then
FMT_dvf(ip:i
ip = ip + 2
End If
9130 Format ('t', 10, ',', a)
End Subroutine dvf_MkFMT
Subroutine dvf_Cols(Col_Text, kO, kl, Last, Icount, Initialize)
! Call dvf Cols('Extraterrestrial Horizontal Radiation1, kO, kl, Icount=l)
Implicit None
Character(Len=*), Intent(In) :: Col_Text
Integer, Intent(In) :: kO, kl
359
-------
! The value of "Last" is unimportant.
! We check only for the presence/absence of the variable.
Logical, Optional, Intent(In) :: Last
Integer, Optional, Intent(In) :: Icount
Logical, Optional, Intent(In) :: Initialize
Character(132) :: qbuf
Integer, Save :: uu = 9010
Integer, Save :: Kcount = 0
Logical, Save :: do return = .False.
If (Present(Initialize)) Then
! By increasing 'uu' we will be able to store
! different sequences in different files.
uu = ULog
do_return = .False.
Return
End If
If (do_return) Then
Return
End If
If (Present(Icount)) Then
Kcount = Icount - 1
End If
If (Present(Last)) Then
do_return = .True.
Return
! Stop '## Stoping in "dvf Cols" as requested1
End If
Write (ULog, '(lx,a)') Trim(qbuf)
End Subroutine dvf Cols
Subroutine hvf_Cols(Col_Text, kO, kl, Xfmt, Last, Icount, Initialize)
! Call hvf_Cols('Extraterrestrial Horizontal Radiation', kO, kl, Icount=l)
Implicit None
Character(Len=*), Intent(In) :: Col_Text
Integer, Intent(In) :: kO, kl
Character(Len=*), Intent(In) :: Xfmt
360
-------
! The value of "Last" is unimportant.
! We check only for the presence/absence of the variable.
Logical, Optional, Intent(In) :: Last
Integer, Optional, Intent(In) :: Icount
Logical, Optional, Intent(In) :: Initialize
Character(Len=132) :: qbuf
Character(Len=50) :: wfmt, wtype
Integer : : j 0, j1
Integer, Save :: ip = 1
uu = ULog
If (Present(Initialize)) Then
do return = .False.
Return
End If
If (do_return) Then
Return
End If
If (Present(Icount)) Then
Kcount = Icount - 1
End If
If (Present(Last)) Then
do return = .True.
Write (ULog, '(/, 2a, /)') '## Hourly values file format: ', Trim(FMT_hvf)
Return
! Stop '## Stoping in "hvf_MkFMT" as requested'
End If
! We assume xFMT with leading and trailing parenthesis.
! We will not check.
wfmt = Xfmt
Call Collapse(wfmt, NewLen=j1)
JO = 2
jl = jl - 1
Select Case(wfmt(j0:j0))
Case('f')
wtype = 'Real'
Case( '!' )
wtype = 'Integer'
Case Default
361
-------
wtype = ''
End Select
If (ip > 1) Then
FMT_hvf(ip:ip+1) = ', '
ip = ip + 2
End If
End Module UtilsS
362
-------
Utils4
! Last change: LSR 16 May 2002 3:21 pm
Module Utils4
Use Global Variables
Use Date_Module
Use Utilsl
Use UtilsS
Implicit None
Implicit None
Integer, Intent(In) :: k_id
Character(Len=*), Intent(In) :: T_code ! Daily value code.
Logical, Intent(Out) :: Xok
Integer :: jday, hhOl, hh24, hh25, iv
Integer :: nn, ierr
Real :: xsum
Type(Val_and_Flag), Dimension(:), Pointer :: VF
Character(Len(Fieldlnfo(1)%Name)), Pointer :: txt
ierr = 0
VF => Xparam(k_id)%Samson_vlO
txt => Fieldlnfo(k_id)%Name
Do jday = jdO, jdl ! step by day
hhOl = (jday-jdO)*Nhours + 1 ! First hour of the day
hh24 = hhOl + 23 ! Last hour of the day (24th)
hh25 = hh24 + 1 ! 25th hour of jday == (jday-jdO+1)*NHours
If (T_code == T_Not_Applicable) Then
VF(hh25) = Val_and_Flag(T_Not_Applicable, Zero, '')
Cycle
End If
! xsum = Sum(VF(hh01:hh24)%v, Mask=(VF(hhOl:hh24)%s /= T_Missing))
! nn = Count(Mask=(VF(hhOl:hh24)%s /= T_Missing))
xsum = Zero
nn = 0
Do iv = hhOl, hh24
Select Case(VF(iv)%s)
Case(T_Missing, T_Not_Applicable, T_Undefined, T_Perpetual_Darkness)
! Do nothing.
363
-------
Case Default
xsum = xsum + VF(iv)%v
nn = nn + 1
End Select
End Do
Else If (T_code == T_Cumulative) Then
! Daily value is cumulative.
VF(hh25) = Val_and_Flag(T_Estimated, xsum, ''
End If
End Subroutine Daily Values
Subroutine Set_Missing_Flags()
!
!
! 1 n
! Sample variance = Sum (x i-x mean)^2
! n - 1 1=1
! Reference:
! [2] Nicholas J. Higham. 1996. Accuracy and Stability of Numerical
! Algorithms. SIAM (Society for Industrial & Applied Mathematics)
! ISBN 0-89871-355-2. Page 13.
I
! Accumulate:
I
! 1 k
! M_k = - •
! k
I
I k
364
-------
! Updating formulae:
! M_k = M_{k-l} +
I
! Q_l = 0
I
! Q_k = Q_{k-l} +
I
! After which:
I
! Sample Mean = M n
! Note that the updating formulae can be written:
I
! M_0 = 0
I
I
! Q 0 = 0
k
Implicit None
Integer
Integer
Integer
Integer
Integer
Integer
Logical
Real
Real
Real
iv, jpar, yyyy, mm, dd, hh
npts, k
n missing, n nan, n out of range, n undefined
n_Perpetual_Darkness
n_l_24, n_25
loop min, loop max, loop inc
erase value, Have Precip Data
minV, maxV
M_k, Q_k, M_kml, x_k
v mean, v variance, v std dev, v min, v max
Type(Val and Flag), Dimension(:), Pointer :: vf
Call Station_With_Missing_Data(pWBAN%WBAN, Have_Precip_Data)
Write (ULog, '(//)')
Loop_Jpar: Do jpar = 1, f_end
365
-------
Select Case(jpar)
Case Default
Case(f_EHR, f_EDNR, f_BAOD, f_GHR, f_DNR, f_TSC, f_OSC, &
f_DBT, f_DPT, f_RH, f_SP, f_WD, f_WS, f_HV, f_CH, &
f_pH2O, f_SD, f_DSLS)
loop_min = 1
loop_max = Ubound(vf,1)
loop inc = 1
Case( &
f_FAO_SG_PET, & ! FAO Short Grass PET, mm/day
f Ep) ! Class A pan Evaporation, mm/day
! Loop only through daily values.
loop_min = NHours
loop_max = Ubound(vf,l)
loop inc = NHours
Case( f_KP_FWS_Evaporation) ! K-P FWS Evaporation, mm/day
! 8 Feb 2002 5:41 pm: until we develop a better formulation,
! kill KP_FWS_Evaporation.
! Do nothing.
Cycle Loop_Jpar
npts = 0
n_missing = 0
n_nan = 0
npts = 0
n_l_24 = 0
n_25 = 0
n_out_of_range = 0
n undefined = 0
n Perpetual Darkness = 0
minV = Fieldlnfo(jpar)%minimum_value
maxV = Fieldlnfo(jpar)%maximum_value
! Number of points.
366
-------
npts = npts + 1
erase value = .False.
I I I I
I I I I
I I I I
!!! If (IsNaN(vf(iv)%v)) Then
! ! ! erase_value = .True.
!!! n_nan = n_nan + 1
! ! ! Call iv to ymdh(iv, yyyy, mm, old, hh)
!!! Write (ULog, 9150) Trim(Fieldlnfo(jpar)%Name), yyyy, mm, dd, hh
I I I
!!! Else
!!! Select Case(vf(iv)%s)
!!! Case(T_Missing)
! ! ! erase_value = .True.
!!! n_missing = n_missing + 1
!!! Case(T_Not_Applicable)
!!! ! Do nothing.
!!! Case (TJJndefined)
!!! n_undefined = n_undefined + 1
!!! Case(T_Perpetual_Darkness)
!!! n_Perpetual_Darkness = n_Perpetual_Darkness + 1
I I I
! ! ! Case Default
!!! in_range = ((minV <= vf(iv)%v) .And. (vf(iv)%v <= maxV))
!!! If (in_range) Then
!!! ! We have a value
! ! ! k = k + 1
!!! x_k = VF(iv)%v
!!! M_kml = M_k
!!! M_k = M_kml + (x_k-M_kml)/Real(k)
!!! Q_k = Q_k + (Real(k-1)/Real(k))*(x_k-M_kml)**2
!!! v min = Min(x k, v min)
!!! v max = Max(x k, v max)
!!! ICall Stat_Add_Point(Xp_ranges(jpar) , vf(iv)%v)
!!! Else
!!! n out of range = n out of range + 1
!!! erase_value = .True.
!!! Call iv_to_ymdh(iv, yyyy, mm, dd, hh)
!!! Write (ULog, 9170) Trim(Fieldlnfo(jpar)%Name), yyyy, mm, dd, hh, vf(i\
!!! End If
! ! ! End Select
!!! End If
! 19 Mar 2002 12:10 pm: It is faster (by about 5 seconds)
367
-------
! to perform all the checks with IFs statements, rather
! than breaking the if-sequence by using "Case Select(vf(iv)%s)"
! and then returning to IFs.
If (IsNaN(vf(iv)%v)) Then
! Produced, by example, by the operation Zero/Zero
erase_value = .True.
n_nan = n_nan + 1
Call iv to ymdh(iv, yyyy, mm, dd, hh)
Write (ULog, 9150) Trim(Fieldlnfo(jpar)%Name), &
iv, yyyy, mm, dd, hh
Format(Ix, 6x, a, ' NaN value for ', &
17, Ix, 14, '-', 12.2, '-', 12.2, 13, 'h')
Else If (vf(iv)%s == T_Missing) Then
erase_value = .True.
n_missing = n_missing + 1
If (Modulo(iv,25) /= 0) Then
n_l_24 = n_l_24 + 1
Else
n_25 = n_25 + 1
Call iv to ymdh(iv, yyyy, mm, dd, hh)
vf(iv) = Val_and_Flag(T_Estimated, Zero, '')
Write (ULog, 9160) Trim(Fieldlnfo(jpar)%Name), &
iv, yyyy, mm, dd, hh
Format(Ix, 6x, a, ' h25 missing for ', &
17, Ix, 14, '-', 12.2, '-', 12.2, 13, 'h; value zeroed.')
End If
Else If (vf(iv)%s == T_Not_Applicable) Then
! Do nothing.
Else If (vf(iv)%s == TJJndefined) Then
n undefined = n undefined + 1
Else If ((minV <= vf(iv)%v) .And. (vf(iv)%v <= maxV)) Then
! Value in range; we have a value.
k = k + 1
x_k = VF(iv)%v
M_kml = M_k
M_k = M_kml + (x_k-M_kml)/Real(k)
Q_k = Q_k + (Real(k-1)/Real(k))*(x_k-M_kml)**2
v min = Min(x k, v min)
v max = Max(x k, v max)
ICall Stat_Add_Point(Xp_ranges(jpar), vf(iv)%v)
Else
368
-------
n out of range = n out of range + 1
erase value = .True.
Call iv_to_ymdh(iv, yyyy, mm, dd, hh)
Write (ULog, 9170) Trim(Fieldlnfo(jpar)%Name), &
iv, yyyy, mm, dd, hh, vf(iv)%v
Format(Ix, 6x, a, ' Out-of-range value for ', &
17, Ix, 14, '-', 12.2, '-', 12.2, 13, 'h', 3x, Ipgl4.6)
End If
If (erase_value) Then
vf(iv)%v = 0
vf(iv)%s = T_Missing
vf(iv)%f = ''
If (jpar == f_HP) Then
Call iv to ymdh(iv, yyyy, mm, dd, hh)
Write (ULog, 9130) Trim(Fieldlnfo(jpar)%Name), iv, yyyy, mm, dd, hh
Format(Ix, 6x, a, ' missing value for ', &
17, Ix, 14, '-', 12.2, '-', 12.2, 13, 'h')
End If
End If
End Do Loop_iv
! Compute mean and variance
v_mean = M_k
If (k >= 2) Then
v_variance = Q_k / Real(k-l)
v_std_dev = Sqrt(v_variance)
Else
! Variance not defined for k<=l
v variance = Huge(Zero)
v_std_dev = Huge(Zero)
End If
Write(ULog, 9190) k, v_mean, v_variance
Format ( &
Ix, ' n ....: ', 10, /, &
Ix, ' Mean .: ', Ipgl4.6, /, &
Ix, ' Var ..: ', Ipgl4.6)
Write(ULog, 9210) v_std_dev
Format (Ix, ' StdDev: ', Ipgl4.6)
If (n missing > 0) Then
Write (ULog, 9250) Trim(Fieldlnfo(jpar)%Name), &
n_missing, npts
Format (Ix, 3x, '?? Missing values for ', a, &
' = ', 10, '/', 10)
369
-------
If (n_l_24 > 0) Write (ULog, 9270) n_l_24
If (n_25 > 0) Write (ULog, 9290) n_25
Format (ix, 3x, ' * Ih - 24h hours missing: ', 10)
Format (Ix, 3x, ' * Summary of the day (h25) missing: ', 10)
End If
If (n_nan > 0) Then
Write (ULog, 9310) Trim(Fieldlnfo(jpar)%Name), &
n_nan, npts
Format (Ix, 3x, '?? NaN values for ', a, ' = ', 10, '/', 10)
End If
If (n_out_of_range > 0) Then
Write (ULog, 9330) Trim(Fieldlnfo(jpar)%Name), &
n out of range, npts
Format (Ix, 3x, '?? Out-of-range values for ', a, '
End If
If (n_undefined > 0) Then
Write (ULog, 9350) Trim(Fieldlnfo(jpar)%Name), &
n_undefined, npts
Format (Ix, 3x, '?? Undefined values for ', a, ' =
End If
If (n_Perpetual_Darkness > 0) Then
Write (ULog, 9370) Trim(Fieldlnfo(jpar)%Name), &
n Perpetual Darkness, npts
Format (Ix, 3x, '?? Perpetual Darkness values for ', a, ' = ', 10, '/', 10)
End If
End Do Loop Jpar
End Subroutine Set_Missing_Flags
End Module Utils4
370
-------
UtilsS
! Last change: LSR 4 Jun 2002
Module Utils5
Contains
Subroutine Station_With_Missing_Data(WBAN, Have_Precip_Data)
Implicit None
Character(Len=*
Logical,
WBAN
Have_Precip_Data
Select Case(WBAN)
Case Default
! Default:
! * All stations have precipitation data
Have_Precip_Data = .True.
! STATIONS WITH LITTLE OR NO HOURLY PRECIPITATION DATA:
! See:
!
-------
& ! Cedar City UT
& ! Washington-Dulles VA
& ! Lufkin TX
'94725') ! Massena NY
! = 26 Apr 2002 2:40 pm
! ^ make all hourly precipitation missing, even when
! present. Fill daily value record (hh == 25) with
! the summary of the day.
Have_Precip_Data = .False.
Case( '22516') ! Kahului HI
! Special case. SAMSON identified this station
! "with little or no hourly ppt data" but we
! augmented the precipitation data with Earthlnfo
! hourly and summary of the day data.
Have_Precip_Data = .True.
End Select
End Subroutine Station With Missing Data
! The following stations have unfixable gaps in
! the precipitation record. We will not issue a year
! if the year's precipitation record is incomplete.
!
! This routine should be called after all the data files
! have been read and before data processing starts. This
! routine consolidates the information gotten from the
! files read variable Year_Data(yyyy)%SAMSON_vlO) with
! its internal information, e.g.,
! Issue This Year(1965:1990) = .True.
Implicit None
Logical, Intent(Out) :: Some Years Missing, All Years Missing
Integer ::
372
-------
Case( '23184') ! Prescott, AZ
! Only the years 1961-1968 are complete.
! Too few years.
Issue_This_Year = .False. ! Kill all years.
Case( '94185') ! Burns, OR
! Only the years 1983-1988 are complete.
! Too few years.
Issue This Year = .False. ! Kill all years.
!!! 4 Jun 2002 10:34 am — precipitation files
!!! corrected. We have all years.
!Case( '03927') ! Fort Worth, TX
! ! Only the years 1975-1990 are complete.
! ! Too few years.
! Issue This Year = .False. ! Kill all years.
373
-------
Do yyyy = MinYear, MaxYear
If (Year_Data(yyyy)%SAMSON_vlO == 0) Then
! SAMSON input routines determined that this year
! is missing. Do not keep this year.
Issue_This_Year(yyyy) = .False.
End If
End Do
End Subroutine Issue Years
! Output: a string suitable for:
! ## Years present: 1961-1990
! ## Years present: ** none **
! ## Years present: 1965-1988
Implicit None
Character(Len=*), Intent(Out) :: Xstr
ymax = Ubound(Issue_This_Year,1) + 1
in_gap = .False.
n_issued = .False.
Xstr = ''
tlen = 1
Do yyyy = MinYear, ymax
If (yyyy < ymax) Then
entry ok = Issue This Year(yyyy)
Else
! yyyy == ymax
! The point is not missing — it does not exist.
entry ok = .False.
End If
If (in_gap) Then
374
-------
If (entry_ok) Then
! Nothing to do.
Else
! The gap run just ended.
Write (Xstr(tlen:), '(14)') yyyy-1
tlen = Len trim(Xstr) + 1
in gap = .False.
n_issued = .True.
End If
Else
If (entry_ok) Then
! Start of a new gap
Write (Xstr(tlen:), '(Ix,14,"-")') yyyy
tlen = Len trim(Xstr) + 1
in gap = .True.
Else
! No missing value and no gap.
! Do nothing.
End If
End If
End Do
If (.Not. n_issued) Then
Write (Xstr(tlen:), '(lx,a)') '** none **'
End If
End Subroutine Str years
375
-------
Test Data for Laramie, Wyoming: average values for August 1987(BurmanandPochop 1994:221)
Element
Wind movement
Solar radiation
Relative Humidity
Air temperature
Dew-point temperature
Elevation
Station pressure
Vapor pressure deficit
Saturation vapor pressure
Actual vapor pressure
A
ap (for class A evap pan)
a
Daily pan evaporation
Daily FWS evaporation
S.I.
261 km/day
5964 Wh/m2/day
43.7%
15.2C
3.0C
2200m
78.1kPa
0.972 kPa
1.727kPa
0.755 kPa
0.1110kPa/°C
0.1225kPa/°C
0.0519 kPa/°C
7.32mm
0.06 mm
English
163 mi/day
514Ly/day
43.7%
59.4F
37.4F
7218 feet
23.04 inches Hg
0.2869 inches Hg
0.5095 inches of Hg
0.2227 inches of Hg
0.0182 inches Hg/°F
0.02007 inches Hg/°F
0.00846 inches Hg/°F
0.29 inches
0.0025 inches
376
-------
SAMSON Station Locations
Weather Bureau Army Navy (WBAN) number, station location (City and State), geographic coordinates (latitude and
longitude) and elevation (mM.S.L.) of the 234 stations available for coordinated climatological dataset for AgDisp,
PRZM, and EXAMS (see ?). Primary stations (those with measured solar radiation data for at least one year) are in bold
type; stations lacking hourly precipitation data are italicized.
WBAN
Station
State
Latitude
Longitude Elevation (m)
03103
03812
03813
03820
03822
03856
03860
03870
03927
03928
03937
03940
03945
03947
04725
04751
11641
12834
12836
12839
12842
Flagstaff
Asheville
Macon
Augusta
Savannah
Huntsville
Huntington
Greenville
Fort Worth
Wichita
Lake Charles
Jackson
Columbia
Kansas City
Binghamton
Bradford
San Juan
Daytona Beach
Key West
Miami
Tampa
AZ
NC
GA
GA
GA
AL
WV
SC
TX
KS
LA
MS
MO
MO
NY
PA
PR
FL
FL
FL
FL
35.1
35.4
32.7
33.4
32.1
34.7
38.4
34.9
32.8
37.7
30.1
32.3
38.8
39.3
42.2
41.8
18.4
29.2
24.6
25.8
28.0
-111.7
-82.5
-83.7
-82.0
-81.2
-86.8
-82.6
-82.2
-97.1
-97.4
-93.2
-90.1
-92.2
-94.7
-76.0
-78.6
-66.0
-81.1
-81.8
-80.3
-82.5
2135
661
110
45
16
190
255
296
164
408
o
3
101
270
315
499
600
19
12
1
2
3
377
-------
WBAN
12844
12912
12916
12917
12919
12921
12924
12960
13722
13723
13729
13733
13737
13739
13740
13741
13748
13781
13865
13866
13873
13874
13876
13877
13880
13881
13882
13883
13889
13891
Station
West Palm
Beach
Victoria
New Orleans
Port Arthur
Brownsville
San Antonio
Corpus Christi
Houston
Raleigh
Greensboro
Elkins
Lynchburg
Norfolk
Philadelphia
Richmond
Roanoke
Wilmington
Wilmington
Meridian
Charleston
Athens
Atlanta
Birmingham
Bristol
Charleston
Charlotte
Chattanooga
Columbia
Jacksonville
Knoxville
State
FL
TX
LA
TX
TX
TX
TX
TX
NC
NC
wv
VA
VA
PA
VA
VA
NC
DE
MS
WV
GA
GA
AL
TN
SC
NC
TN
SC
FL
TN
Latitude
26.7
28.9
30.0
30.0
25.9
29.5
27.8
30.0
35.9
36.1
38.9
37.3
36.9
39.9
37.5
37.3
34.3
39.7
32.3
38.4
34.0
33.7
33.6
36.5
32.9
35.2
35.0
34.0
30.5
35.8
Longitude
-80.1
-96.9
-90.3
-94.0
-97.4
-98.5
-97.5
-95.4
-78.8
-80.0
-79.9
-79.2
-76.2
-75.3
-77.3
-80.0
-77.9
-75.6
-88.8
-81.6
-83.3
-84.4
-86.8
-82.4
-80.0
-80.9
-85.2
-81.1
-81.7
-84.0
Elevation (m)
6
32
3
7
6
242
13
33
134
270
594
279
9
9
50
358
9
24
94
290
244
315
192
459
12
234
210
69
9
299
378
-------
WBAN
13893
13894
13895
13897
13957
13958
13959
13962
13963
13964
13966
13967
13968
13970
13985
13994
13995
13996
14607
14733
14734
14735
14737
14739
14740
14742
14745
14751
14764
14765
Station
Memphis
Mobile
Montgomery
Nashville
Shreveport
Austin
Waco
Abilene
Little Rock
Fort Smith
Wichita Falls
Oklahoma
City
Tulsa
Baton Rouge
Dodge City
St. Louis
Springfield
Topeka
Caribou
Buffalo
Newark
Albany
Allentown
Boston
Hartford
Burlington
Concord
Harrisburg
Portland
Providence
State
TN
AL
AL
TN
LA
TX
TX
TX
AR
AR
TX
OK
OK
LA
KS
MO
MO
KS
ME
NY
NJ
NY
PA
MA
CT
VT
NH
PA
ME
RI
Latitude
35.1
30.7
32.3
36.1
32.5
30.3
31.6
32.4
34.7
35.3
34.0
35.4
36.2
30.5
37.8
38.8
37.2
39.1
46.9
42.9
40.7
42.8
40.7
42.4
41.9
44.5
43.2
40.2
43.7
41.7
Longitude
-90.0
-88.3
-86.4
-86.7
-93.8
-97.7
-97.2
-99.7
-92.2
-94.4
-98.5
-97.6
-95.9
-91.2
-100.0
-90.4
-93.4
-95.6
-68.0
-78.7
-74.2
-73.8
-75.4
-71.0
-72.7
-73.2
-71.5
-76.9
-70.3
-71.4
Elevation (m)
87
67
62
180
79
189
155
534
81
141
314
397
206
23
787
172
387
270
190
215
9
89
117
5
55
104
105
106
19
19
379
-------
WBAN
Station
State
Latitude Longitude Elevation (m)
14768
14771
14777
14778
14820
14821
14826
14827
14836
14837
14839
14840
14842
14847
14848
14850
14852
14860
14891
14895
14898
14913
14914
14918
14920
14922
14923
14925
14926
14933
Rochester
Syracuse
Wilkes-Barre
Williamsport
Cleveland
Columbus
Flint
Fort Wayne
Lansing
Madison
Milwaukee
Muskegon
Peoria
Sault Ste. Marie
South Bend
Traverse City
Youngstown
Erie
Mansfield
Akron
Green Bay
Duluth
Fargo
International
Falls
La Crosse
Minneapolis
Moline
Rochester
Saint Cloud
Des Moines
NY
NY
PA
PA
OH
OH
MI
IN
MI
WI
WI
MI
IL
MI
IN
MI
OH
PA
OH
OH
WI
MN
ND
MN
WI
MN
IL
MN
MN
IA
43.1
43.1
41.3
41.3
41.4
40.0
43.0
41.0
42.8
43.1
43.0
43.2
40.7
46.5
41.7
44.7
41.3
42.1
40.8
40.9
44.5
46.8
46.9
48.6
43.9
44.9
41.5
43.9
45.6
41.5
-77.7
-76.1
-75.7
-77.1
-81.9
-82.9
-83.7
-85.2
-84.6
-89.3
-87.9
-86.3
-89.7
-84.4
-86.3
-85.6
-80.7
-80.2
-82.5
-81.4
-88.1
-92.2
-96.8
-93.4
-91.3
-93.2
-90.5
-92.5
-94.1
-93.7
169
124
289
243
245
254
233
252
256
262
211
191
199
221
236
192
361
225
395
377
214
432
274
361
205
255
181
402
313
294
380
-------
WBAN
Station
State
Latitude Longitude Elevation (m)
14935
14936
14940
14941
14943
14944
14991
21504
22516
22521
22536
23023
23034
23042
23044
23047
23050
23061
23065
23066
23129
23153
23154
23155
23160
23161
23169
23174
23183
23184
Grand Island
Huron
Mason City
Norfolk
Sioux City
Sioux Falls
Eau Claire
Hilo
Kahului
Honolulu
Lihue
Midland/Odessa
San Angelo
Lubbock
El Paso
Amarillo
Albuquerque
Alamosa
Goodland
Grand Junction
Long Beach
Tonopah
Ely
Bakersfield
Tucson
Daggett
Las Vegas
Los Angeles
Phoenix
Prescott
NE
SD
IA
NE
IA
SD
WI
HI
HI
HI
HI
TX
TX
TX
TX
TX
NM
CO
KS
CO
CA
NV
NV
CA
AZ
CA
NV
CA
AZ
AZ
41.0
44.4
43.2
42.0
42.4
43.6
44.9
19.7
20.9
21.3
22.0
31.9
31.4
33.7
31.8
35.2
35.1
37.5
39.4
39.1
33.8
38.1
39.3
35.4
32.1
34.9
36.1
33.9
33.4
34.7
-98.3
-98.2
-93.3
-97.4
-96.4
-96.7
-91.5
-155.1
-156.4
-157.9
-159.4
-102.2
-100.5
-101.8
-106.4
-101.7
-106.6
-105.9
-101.7
-108.5
-118.2
-117.1
-114.9
-119.1
-110.9
-116.8
-115.2
-118.4
-112.0
-112.4
566
393
373
471
336
435
273
11
15
5
45
871
582
988
1194
1098
1619
2297
1124
1475
17
1653
1906
150
779
588
664
32
339
1531
381
-------
WBAN
23185
23188
23232
23234
23273
24011
24013
24018
24021
24023
24025
24027
24028
24029
24033
24089
24090
24121
24127
24128
24131
24143
24144
24146
24153
24155
24156
24157
24221
24225
Station
Reno
San Diego
Sacramento
San
Francisco
Santa Maria
Bismarck
Minot
Cheyenne
Lander
North Platte
Pierre
Rock Springs
Scottsbluff
Sheridan
Billings
Casper
Rapid City
Elko
Salt Lake City
Winnemucca
Boise
Great Falls
Helena
Kalispell
Missoula
Pendleton
Pocatello
Spokane
Eugene
Medford
State
NV
CA
CA
CA
CA
ND
ND
WY
WY
NE
SD
WY
NE
WY
MT
WY
SD
NV
UT
NV
ID
MT
MT
MT
MT
OR
ID
WA
OR
OR
Latitude
39.5
32.7
38.5
37.6
34.9
46.8
48.3
41.2
42.8
41.1
44.4
41.6
41.9
44.8
45.8
42.9
44.1
40.8
40.8
40.9
43.6
47.5
46.6
48.3
46.9
45.7
42.9
47.6
44.1
42.4
Longitude Elevation (m)
-119.8
-117.2
-121.5
-122.4
-120.5
-100.8
-101.3
-104.8
-108.7
-100.7
-100.3
-109.1
-103.6
-107.0
-108.5
-106.5
-103.1
-115.8
-112.0
-117.8
-116.2
-111.4
-112.0
-114.3
-114.1
-118.9
-112.6
-117.5
-123.2
-122.9
1341
9
8
5
72
502
522
1872
1696
849
526
2056
1206
1209
1088
1612
966
1547
1288
1323
874
1116
1188
904
972
456
1365
721
109
396
382
-------
WBAN
24227
24229
24230
24232
24233
24243
24283
24284
25308
25339
25501
25503
25624
25713
26411
26415
26425
26451
26510
26528
26533
26615
26616
26617
27502
41415
93037
93058
93129
93193
Station
Olympia
Portland
Redmond/Bend
Salem
Seattle/Tacoma
Yakima
Arcata
North Bend
Annette
Yakutat
Kodiak
King Salmon
Cold Bay
St Paul Is.
Fairbanks
Big Delta
Gulkana
Anchorage
Mcgrath
Talkeetna
Settles
Bethel
Kotzebue
Nome
Barrow
Guam
Colorado Springs
Pueblo
Cedar City
Fresno
State
WA
OR
OR
OR
WA
WA
CA
OR
AK
AK
AK
AK
AK
AK
AK
AK
AK
AK
AK
AK
AK
AK
AK
AK
AK
PI
CO
CO
UT
CA
Latitude
47.0
45.6
44.3
44.9
47.5
46.6
41.0
43.4
55.0
59.5
57.8
58.7
55.2
57.2
64.8
64.0
62.2
61.2
63.0
62.3
66.9
60.8
66.9
64.5
71.3
13.6
38.8
38.3
37.7
36.8
Longitude
-122.9
-122.6
-121.2
-123.0
-122.3
-120.5
-124.1
-124.3
-131.6
-139.7
-152.3
-156.7
-162.7
-170.2
-147.9
-145.7
-145.5
-150.0
-155.6
-150.1
-151.5
-161.8
-162.6
-165.4
-156.8
-144.8
-104.7
-104.5
-113.1
-119.7
Elevation (m)
61
12
940
61
122
325
69
5
34
9
34
15
29
7
138
388
481
35
103
105
205
46
5
7
4
110
1881
1439
1712
100
383
-------
WBAN
Station
State
Latitude Longitude Elevation (m)
93721
93729
93730
93738
93805
93814
93815
93817
93819
93820
93821
93822
93842
93987
94008
94018/23062
94185
94224
94240
94702
94725
94728/14732
94746
94814
94822
94823
94830
94846
94847
94849
Baltimore
Cape Hatteras
Atlantic City
Sterling (Washington-Dulles Airpt.)
Tallahassee/Apalachicola
Covington
Dayton
Evansville
Indianapolis
Lexington
Louisville
Springfield
Columbus
Lufkin
Glasgow
Boulder/Denver
Burns
Astoria
Quillayute
Bridgeport
Massena
New York (LaGuardia Airpt.)
Worcester
Houghton
Rockford
Pittsburgh
Toledo
Chicago
Detroit
Alpena
MD
NC
NJ
VA
FL
KY
OH
IN
IN
KY
KY
IL
GA
TX
MT
CO
OR
OR
WA
CT
NY
NY
MA
MI
IL
PA
OH
IL
MI
MI
39.2
35.3
39.5
39.0
30.4
39.1
39.9
38.1
39.7
38.0
38.2
39.8
32.5
31.2
48.2
39.8
43.6
46.2
48.0
41.2
44.9
40.8
42.3
47.2
42.2
40.5
41.6
41.8
42.4
45.1
-76.7
-75.6
-74.6
-77.5
-84.4
-84.7
-84.2
-87.5
-86.3
-84.6
-85.7
-89.7
-85.0
-94.8
-106.6
-104.9
-119.1
-123.9
-124.6
-73.1
-74.9
-73.9
-71.9
-88.5
-89.1
-80.2
-83.8
-87.8
-83.0
-83.6
47
2
20
82
21
271
306
118
246
301
149
187
136
96
700
1610
1271
7
55
2
63
11
301
329
221
373
211
190
191
210
384
-------
WBAN Station State Latitude Longitude Elevation (m)
94860 Grand Rapids MI 42.9 -85.5 245
94910 Waterloo IA 42.6 -92.4 265
94918/14942 Omaha NE 41.3 -95.9 298
385
-------
References
Allen, R. G. 1996. Assessing integrity of weather data for reference evapotranspiration estimation.
Journal of Irrigation and Drainage Engineering 122:97-106.
Allen, R. G. 2000. REF-ET: Reference Evapotranspiration Calculation Software for FAO and ASCE
Standardized Equations, Version 2.0. University of Idaho, Kimberly.
Allen, R. G., L. S. Pereira, D. Raes, and M. Smith. 1998. Crop evapotranspiration: Guidelines for
computing crop water requirements. FAO Irrigation and Drainage Paper 56, Food and
Agriculture Organization of the United Nations, Rome.
ASCE. 1996. Hydrology Handbook (Manual No. 28), Second edition. American Society of Civil
Engineers, New York.
Brutsaert, W. 1982. Evaporation into the Atmosphere: Theory, History, and Applications. D. Reidel
Publishing Co., Dordrecht, The Netherlands.
Burman, R., andL. O. Pochop. 1994. Evaporation, Evapotranspiration and Climatic Data. Elsevier,
Amsterdam.
Burns, L. A. 2000. Exposure Analysis Modeling System (Exams): User Manual and System
Documentation. EPA/600/R-00/081, U.S. Environmental Protection Agency, Office of
Research and Development, National Exposure Research Laboratory, Research Triangle
Park, North Carolina, USA.
Carousel, R. F., J. C. Imhoff, P. R. Hummel, J. M. Cheplick, and A. S. Donigian, Jr. 2005. PRZM-3,
A Model for Predicting Pesticide and Nitrogen Fate in the Crop Root and Unsaturated Soil
Zones: Users Manual for Release 3.12.2. Athens, Georgia.
EPA. 2000. Meteorological Monitoring Guidance for Regulatory Modeling Applications. EPA-
454/R-99-005, U.S. Environmental Protection Agency, Office of Air Quality Planning and
Standards, Research Triangle Park.
Farnsworth, R. K., E. S. Thompson, and E. L. Peck. 1982. Evaporation Atlas for the Contiguous 48
United States. NOAA Technical Report NWS 33, Office of Hydrology, National Weather
Service, Washington.
Hanson, C. L., K. A. Cumming, D. A. Woolhiser, and C. W. Richardson. 1993. Program for Daily
Weather Simulation. Water Resources Investigations Rep. 93-4018,U.S. Geological Survey,
Denver.
Hanson, C. L., K. A. Cumming, D. A. Woolhiser, and C. W. Richardson. 1994. Microcomputer
Program for Daily Weather Simulation in the Contiguous United States. Agricultural
Research Service ARS-114, U.S. Department of Agriculture, Boise.
Harbeck, G. E., andF. W. Kennon. 1954. Lake Hefner Studies Technical Report. Professional Paper
269, U.S. Geological Survey, Washington.
Harrison, L. P. 1963. Fundamental concepts and definitions relating to humidity. Pages 3-80 in A.
Wexler, editor. Humidity and Moisture. Reinhold Publishing Company, New York.
Johnson, G. L., C. L. Hanson, S. P. Hardegree, and E. B. Ballard. 1996. Stochastic weather
simulation: Overview and analysis of two commonly used models. Journal of Applied
Meteorology 35:1878-1896.
Kohler, M. A., T. J. Nordenson, and W. E. Fox. 1955. Evaporation from Pans and Lakes. Research
Paper 38, U.S. Department of Commerce, Weather Bureau, Washington.
Kohler, M. A., and L. H. Parmele. 1967. Generalized estimates of free-water evaporation. Water
386
-------
Resources Research 3:997-1005.
Lamoreux, W. W. 1962. Modern evaporation formulae adapted to computer use. Monthly Weather
Review 90:26-28.
Merkel, W. H. 1988. Appendix A, Revision of Reservoir Operations Study Computer Program &
User Manual. Technical Release 210-VI-TR-19A, U.S. Department of Agriculture, Soil
Conservation Service, Washington.
Nicks, A. D., and G. A. Gander. 1994. CLIGEN: A weather generator for climate inputs to water
resource and other models, in Proceedings of the Fifth International Conference on
Computers in Agriculture. American Society of Agricultural Engineers, Orlando.
Penman, H. L. 1948. Natural evaporation from open water, bare soil and grass. Proceedings of the
Royal Society of London, Series A: Mathematical and Physical Sciences 193:120-145.
Penman, H. L. 1963. Vegetation and Hydrology. Technical Communication 53, Commonwealth
Agricultural Bureaux, Farnham Royal, Bucks, England.
Richardson, C. W., and D. A. Wright. 1984. WGEN: A Model for Generating Daily Weather
Variables. Agricultural Research Service ARS-8, U.S. Department of Agriculture,
Washington.
Semenov, M. A., andE. M. Barrow. 1997. Use of a stochastic weather generator in the development
of climate change scenarios. Climate Change 35:397-414.
Semenov, M. A., R. J. Brooks, E. M. Barrow, and C. W. Richardson. 1998. Comparison of the
WGEN and LARS-WG stochastic weather generators for diverse climates. Climate Research
10:95-107.
387
-------
vvEPA
United States
Environmental Protection
Agency
Office of Research
and Development (8101R)
Washington, DC 20460
Official Business
Penalty for Private Use
$300
EPA 600/R-07/053
May 2007
www.epa.gov
Please make all necessary changes on the below label,
detach or copy, and return to the address in the upper
left-hand corner.
If you do not wish to receive these reports CHECK HERE
D; detach, or copy this cover, and return to the address in
the upper left-hand corner.
PRESORTED STANDARD
POSTAGE & FEES PAID
EPA
PERMIT No. G-35
Recycled/Recyclable
Printed with vegetable-based ink on
paper that contains a minimum of
50% post-consumer fiber content
processed chlorine free
-------