c Program GenerateStationParameters c c Purpose: Using data at ftp://hydrolab.arsusda.gov/pub/nicks/countries c generate the monthly temperature and precipitation statistics c for a Cligen ".par" file. (Incidentally .par files posted at c that website appear to be grossly incorrect.) c c Note(s): For purpose of calculating P(W|W) & P(W|D) we are assuming c the data is order by date. Missing days are simply ignored. c A day with any precip, regardless of amount, is considered c a wet day. To eliminate conditions that prohibit transition c between dry and wet states, P(W|D)=0 and P(W|W)=1 are not c permitted. Skewness is without units: the value calculated c using tenths of a millimeter is _identical_ to that for c hundredths of an inch -- really. c c Author: Charles R. Meyer c c Date: 9/03/2003 - 9/09/2003; 9/17/2003 - 9/18/2003; 9/22/2003; 07/24/06. c c Version: 1.005 c Changes: c - Modified Data Read line to show which line is being read, c when a read error occurs. c CRM -- 07/24/06 & 8/02/06. c c Version: 1.004 c Changes: c - Corrected a "divide-by-zero" which occurs when a month has no c precipitation for any day in the period of record. In that c case, P(W|W) is set to zero; P(W|D) is set to 0.01; and the c amount of precip per event is set to 0.01 . c CRM -- 10/29/03. c c Version: 1.003 c Changes: c - Added the string "TYPE= 0" to line 2 of the .TOP file. c From the SCS TR-55, this parameter is used in single storm runs c to indicate which of four climates occurs at this site: c 1 & 2 (SCS 1 & 1A) - Pacific maritime with wet winters and dry c summers. 4 (SCS 3) - Areas where large 24-hr rainfall amounts c result from tropical storms in the Gulf of Mexico and the c Atlantic coast. 3 (SCS 2) - the rest of the U. S. Later in c FindMatch, this value gets replaced with that of the surrogate. c CRM -- 9/24/03. c c Inputs: c id_sta -- Station (numeric) ID c id_cty -- Country (numeric) ID c stanam -- Station Name c lat_d -- Station Latitude (degrees) c lat_m -- Station Latitude (minutes) c lon_d -- Station Longitude (degrees) c lon_m -- Station Longitude (minutes) c elev -- Station Elevation (m) c yr -- year of observation c mo -- month of observation c da -- day of observation c tmax -- daily maximum temperature (C) c tmin -- daily minimum temperature (C) c pcp -- daily precipitation (mm) c c Derived Values: c latt -- Station Latitude (in decimal format) c lonn -- Station Longitude (in decimal format) c y_pcp -- Precip Yesterday (0.0=dry; 1.0=wet) c dd_knt -- count of dry day followed by a dry day c dw_knt -- count of dry day followed by a wet day c wd_knt -- count of wet day followed by a dry day c ww_knt -- count of wet day followed by a wet day c tx_knt -- count of days having max temperature value c tn_knt -- count of days having min temperature value c p_ww -- P(W|W); ie, Monthly probability of wet after wet c p_dw -- P(W|D); ie, Monthly probability of wet after dry c pw -- Monthly probability of precip c precip -- total precip for the month c monpcp -- Monthly average precip per event c t_max -- Monthly average maximum temperature c t_min -- Monthly average minimum temperature c sd_pcp -- Monthly standard deviation for precip c sd_tmx -- Monthly standard deviation for maximum temperature c sd_tmn -- Monthly standard deviation for minimum temperature c c Used for population-based ('n', not 'n-1') calculations of SD c tempx -- Sum of max temp values for each month c tempn -- Sum of min temp values for each month c precip -- Sum of precip values for each month c tempx2 -- Sum of max-precip-values squared for each month c tempn2 -- Sum of min-precip-values squared for each month c precp2 -- Sum of precip-values squared for each month c integer id_sta, id_cty, lat_d, lat_m, lon_d, lon_m integer yr, mo, da integer elev integer dd_knt(12),dw_knt(12),wd_knt(12),ww_knt(12) integer years integer lineno real rdays(12) real tmax, tmin, pcp, y_pcp real latt, lonn real p_ww(12), p_dw(12), precip(12), monpcp(12), pw(12) real t_max(12), t_min(12), tempx(12), tempn(12), skewx(12) real precp2(12), tempx2(12), tempn2(12) real sd_pcp(12), sd_tmx(12), sd_tmn(12) real tx_knt(12), tn_knt(12) character*3 st_code character*20 infile, outfile character*21 country character*25 city character*46 stanam, string character*80 buffer c data y_pcp /-1./ data dd_knt /12*0/ data dw_knt /12*0/ data wd_knt /12*0/ data ww_knt /12*0/ data tx_knt /12*0./ data tn_knt /12*0./ data precip /12*0./ data tempx /12*0./ data tempn /12*0./ data precp2 /12*0./ data tempx2 /12*0./ data tempn2 /12*0./ data sd_tmx /12*0./ data sd_tmn /12*0./ data sd_pcp /12*0./ data skewx /12*0./ c 1100 format(i2,i3,a46,i3,i2,2x,i3,i2,i6) 1101 format(3i2,2(f5.1,2x),f5.1) 1102 format(2x,i2,16x,f5.1) 1103 format(a80) c lineno = 0 write(*,*) write(*,*) "Welcome to the program GenStPar version 1.005." write(*,*) " updated August 2006" write(*,*) write(*,*) "This program builds the top of a Cligen station" write(*,*) "parameter file from the international daily" write(*,*) "temperature and precipitation data available on" write(*,*) "the webpage at" write(*,*) "ftp://hydrolab.arsusda.gov/pub/nicks/countries" write(*,*) write(*,*) "Please enter the name of your raw data file." write(*,*) "Example: 76151.GDS" write(*,*) read(*,'(a20)') infile open(unit=10, file=infile, status='old', err=9992) read(10,1100) id_sta,id_cty,stanam,lat_d,lat_m,lon_d,lon_m,elev lineno = lineno + 1 c write(string,'(a46)') stanam c-----eliminate leading blanks from stanam i_rej = 0 11 continue i_rej = i_rej+1 if(i_rej.lt.46 .and. string(i_rej:i_rej).eq.' ' ) goto 11 write(stanam,'(a)') string(i_rej:46) c c-----use first 21 characters of station name as the country: last 23 as city read(stanam,'(a21,a25)') country,city c-----try to find a "match" so 3-letter code can be substituted for country open(unit=11, file='WEPP_CountryCodes.txt',status='old',err=9993) 12 continue read(11,'(a3,1x,a21)',end=13) st_code, string if(string(1:21-i_rej) .eq. country(1:21-i_rej)) then c-------alter stanam to format used by WEPP GUI interface write(stanam,'(a25,1x,a3)') city, st_code else goto 12 endif c 13 continue write(outfile,'(a20)') infile i = 21 14 continue i = i-1 if(i.gt.1 .and. outfile(i:i).ne.'.' ) goto 14 outfile(i+1:i+3)='top' c open(unit=12, file=outfile, status='unknown', err=9994) c c ---- Convert Longitudes to E/W pos/neg longitude: c-----eastern longitudes are positive up to 180 if (lon_d .gt. 180) then lon_d = 360 - lon_d lon_m = 60 - lon_m c-----western longitudes are negative else if (lon_d .gt. 0) then lon_d = -1 * lon_d lon_m = -1 * lon_m endif latt = lat_d + lat_m/60.0 lonn = lon_d + lon_m/60.0 c 10 continue CC read(10,1101,end=9996,err=9991) yr,mo,da,tmax,tmin,pcp read(10,1103,end=9996) buffer lineno = lineno + 1 C if(lineno.le.10) write(*,*) buffer read(buffer,1101,err=9991) yr,mo,da,tmax,tmin,pcp C if(lineno.le.10) write (*,*) yr,mo,da,tmax,tmin,pcp if(tmax .gt. -99.) then c----------temperatures are in tenths of a degree C. tempx(mo) = tempx(mo) + tmax tempx2(mo) = tempx2(mo) + tmax**2 tx_knt(mo) = tx_knt(mo) + 1. endif c----------temperatures (C) if(tmin .gt. -99.) then tempn(mo) = tempn(mo) + tmin tempn2(mo) = tempn2(mo) + tmin**2 tn_knt(mo) = tn_knt(mo) + 1. endif c----------precip (mm) if(pcp .gt. 0.0) then precip(mo) = precip(mo) + pcp precp2(mo) = precp2(mo) + pcp**2 if(y_pcp .gt. 0.0) then ww_knt(mo) = ww_knt(mo) + 1 else dw_knt(mo) = dw_knt(mo) + 1 endif c--------non-missing days else if(pcp .eq. 0.0) then if(y_pcp .gt. 0.0) then wd_knt(mo) = wd_knt(mo) + 1 else dd_knt(mo) = dd_knt(mo) + 1 endif endif c--------update yesterday's precip state for nonmissing data if(pcp .ge. 0.0) y_pcp = pcp goto 10 c goto 9996 c------ ERRORS: 9991 continue write(*,*) ' Error reading input data on line:', lineno C write(*,*) lineno write(*,1103) buffer goto 9999 9992 continue write(*,*) ' File:',infile,' not found.' write(*,*) ' It needs to be in the same directory (folder) as prog &ram GenStPar.' goto 9999 9993 continue write(*,*) ' File: WEPP_CountryCodes.txt not found.' write(*,*) ' It needs to be in the same directory (folder) as prog &ram GenStPar.' goto 9999 9994 continue write(*,*) ' File:',outfile,' cannot be opened.' write(*,*) ' It should be in the same directory (folder) as progra &m GenStPar.' goto 9999 9995 continue write(*,*) ' Cannot Write to File:',outfile write(*,*) ' It should be in the same directory (folder) as progra &m GenStPar.' goto 9999 c 9996 continue c------Calculate Means & SD's for Temperatures & Precipitation do 700, i=1, 12 temp_d1 = float(ww_knt(i) + wd_knt(i)) temp_d2 = float(dw_knt(i) + dd_knt(i)) c--------trap divide by zero if(temp_d1 .lt. 0.00001) temp_d1 = 0.00001 if(temp_d2 .lt. 0.00001) temp_d2 = 0.00001 p_ww(i) = ww_knt(i) / temp_d1 p_dw(i) = dw_knt(i) / temp_d2 c--------trap conditions that prevent transition between states if(p_ww(i) .gt. 0.99) p_dw(i) = 0.99 if(p_dw(i) .lt. 0.01) p_dw(i) = 0.01 c--------monthly average precip per event rdays(i) = float(ww_knt(i) + dw_knt(i)) if(rdays(i).gt.0.0 .and. precip(i).gt.0.0) then monpcp(i) = precip(i) / (rdays(i) * 25.4) else monpcp(i) = 0.01 endif c--------monthly average temperatures t_max(i) = tempx(i)*9/(5*tx_knt(i))+32.0 t_min(i) = tempn(i)*9/(5*tn_knt(i))+32.0 c--------monthly SD precip per event (using 'n', not 'n-1'; ie, "population") C sd_pcp(i) = sqrt(precp2(i)/rdays(i)-(precip(i)/rdays(i))**2)/ C & 25.4 c--------monthly SD temperatures (using 'n', not 'n-1') sd_tmx(i) = sqrt(tempx2(i)/tx_knt(i) - (tempx(i)/tx_knt(i)) & **2)*9/5.0 sd_tmn(i) = sqrt(tempn2(i)/tn_knt(i) - (tempn(i)/tn_knt(i)) & **2)*9/5.0 700 continue years = int(temp_d1 + temp_d2) / 31 c c------Calculate SD for Monthly Precipitation (using 'n-1') rewind(10) read(10,*) 20 continue read(10,1102,end=9997) mo,pcp if(pcp .gt. 0.0) then if(rdays(mo) .gt. 1.0) then sd_pcp(mo)=sd_pcp(mo)+((pcp/25.4-monpcp(mo))**2)/ & (rdays(mo)-1) else sd_pcp(mo) = 0.0 endif endif goto 20 c c-----finish SD calculation 9997 continue do 21, i=1, 12 sd_pcp(i) = sqrt(sd_pcp(i)) 21 continue c c------Calculate Skewness for Monthly Precipitation rewind(10) read(10,*) 30 continue read(10,1102,end=9998) mo,pcp if(pcp .gt. 0.0) then skewx(mo) = skewx(mo)+((pcp/25.4-monpcp(mo))/sd_pcp(mo))**3 endif goto 30 c c-----finish skewness calculation 9998 continue do 31, i=1, 12 if(rdays(i).gt.2) then skewx(i) = rdays(i)*skewx(i)/((rdays(i)-1)*(rdays(i)-2)) else skewx(i) = 0.0 endif 31 continue c write(12,2000,err=9995) stanam write(12,2001,err=9995) latt,lonn,float(years),0 c-----elevation is converted from feet to meters write(12,2002,err=9995) float(elev)*3.28,0.0,0.0 write(12,2003,err=9995) monpcp write(12,2004,err=9995) sd_pcp write(12,2005,err=9995) skewx write(12,2006,err=9995) p_ww write(12,2007,err=9995) p_dw write(12,2008,err=9995) t_max write(12,2009,err=9995) t_min write(12,2010,err=9995) sd_tmx write(12,2011,err=9995) sd_tmn c 2000 format(1x,a46) c2001 format(" LATT=",f7.2," LONG=",f7.2," YEARS=",f4.0) 2001 format(" LATT=",f7.2," LONG=",f7.2," YEARS=",f4.0," TYPE=",i2) c2002 format(" ELEVATION =",f6.0) 2002 format(" ELEVATION =",f6.0," TP5 =",f5.2," TP6=",f5.2) 2003 format(" MEAN P ",12f6.2) 2004 format(" S DEV P",12f6.2) 2005 format(" SKEW P ",12f6.2) 2006 format(" P(W|W) ",12f6.2) 2007 format(" P(W|D) ",12f6.2) 2008 format(" TMAX AV",12f6.2) 2009 format(" TMIN AV",12f6.2) 2010 format(" SD TMAX",12f6.2) 2011 format(" SD TMIN",12f6.2) write(*,*) write(*,*) "Congratulations. You have completed a run of the" write(*,*) "GenStPar program. Next you will need to find a" write(*,*) "surrogate station to supply the climate parameters" write(*,*) "missing from the GDS datasets. To find a matching" write(*,*) "US station you can now run the program FindMatch." write(*,*) "It will suggest 10 US stations that should be good" write(*,*) "surrogates. US data is available from the Cligen" write(*,*) "website http://horizon.nserl.purdue.edu/Cligen/" write(*,*) "When you find an acceptable station, simply copy its" write(*,*) "file, and use a text editor to replace the top 12" write(*,'(" lines with the file that we just built: ",a)')outfile c 9999 continue close(10) close(11) close(12) end