c isd_process.f c Program to process ISD data, calculate apparent temperatures, update c database files and put on ftp site c --------------------------------------------------------------------- c Type declaration character cY*4,cM*2,cMeta(500)*5,cStation(500)*18,cStat*18, : cSt(500)*2,cStID*5,cS*2,cOtype*5,cSyr*4,cSmon*2, : cMxfl(31)*1,cMnfl(31)*1,cMxS(31)*1,cMnS(31)*1, : cPath*28,cJunk*23,cStnid*5,cMon(12)*3 integer iDtarr(3),iStation(500),iTz(500),iMlength(12),iL(31), : iEmax(31),iEmin(31),iMaxT(31),iMinT(31),iEtemp(31,24), : iN(500,2),iFyear(6,500),iDays(6,500,1940:2020,12), : iMcount(3),iMhwave(3),iCount(6),iHwave(3),iLyear(6,12), : iVal(12),i4count(3) real rAlt(500),rPress(31,24),rTemp(31,24),rDp(31,24), : rEtemp(31,24),rAthresh(500,6),rFtemp(31,24),rEval(31,3), : rTval(31,3) data iMlength / 31,28,31,30,31,30,31,31,30,31,30,31 / data cMon / 'Jan','Feb','Mar','Apr','May','Jun', : 'Jul','Aug','Sep','Oct','Nov','Dec' / cPath = '/cmb1/monitoring/heat_index/' c --------------------------------------------------------------------- c Get current date and build input file name. ISD data is listed in c UTC dates and times so must grep from two monthly files to aquire c a full month of local time data call iDate(iDtarr) iYear = iDtarr(3) * iMonth = iDtarr(2) - 3 iMonth = iDtarr(2) - 2 if (iMonth .lt. 1) then iYear = iYear - 1 iMonth = 12 + iMonth endif iSyr = iDtarr(3) * iSmon = iDtarr(2) - 2 iSmon = iDtarr(2) - 1 if (iSmon .lt. 1) then iSyr = iSyr - 1 iSmon = 12 + iSmon endif ****************************************** * iMonth = 8 * iSmon = 9 ****************************************** write(cY,'(i4.4)') iYear write(cM,'(i2.2)') iMonth write(cSyr,'(i4.4)') iSyr write(cSmon,'(i2.2)') iSmon if ((mod(iYear,4) .eq. 0) .and. (iMonth .eq. 2)) then iMlen = 29 else iMlen = iMlength(iMonth) endif c --------------------------------------------------------------------- c Make sure data files are available * inquire(file=cPath//'data/ISD'//cY//cM//'-V030.dat',exist=lExist) inquire(file=cPath//'data/ISH'//cMon(iMonth)//cY//'.dat', : exist=lExist) if (lExist .eq. .false.) then write(6,*) 'First months data file does not exist' go to 1001 endif * inquire(file=cPath//'data/ISD'//cSyr//cSmon//'-V030.dat', * : exist=lExist) inquire(file=cPath//'data/ISH'//cMon(iSmon)//cSyr//'.dat', : exist=lExist) if (lExist .eq. .false.) then write(6,*) 'Second months data file does not exist' go to 1001 endif c --------------------------------------------------------------------- c Make backups of output files call system('cp -f '//cPath//'output/*.ts '//cPath// : 'output/backup') call system('cp -f '//cPath//'app_temps/*.apptemp '//cPath// : 'app_temps/backup') c --------------------------------------------------------------------- c Read metadata file and put into array open(2,file=cPath//'programs/metadata.dat',status='old') n = 0 do 100,i=1,500 read(2,1,end=1000) iMet,cStat,cS,rLat,rLong,iAlt,iTime 1 format(i5,3x,a18,a2,1x,2f7.2,42x,i6,i4) n = n + 1 write(cMeta(n),'(i5.5)') iMet do 101,ii=1,20 if (cStat(ii:ii) .eq. ',') then cStation(n) = cStat(1:ii-1) iStation(n) = ii-1 go to 3000 else if (cStat(ii:ii) .eq. ' ') then cStat(ii:ii) = '_' endif 101 continue 3000 cSt(n) = cS iTz(n) = iTime 100 continue 1000 iMeta = n close (2) c --------------------------------------------------------------------- c Open and read threshold values (plus number of obs in 1961 - 90 c period) open(2,file=cPath//'programs/app_temp_thr.dat',status='unknown') do 105,i=1,iMeta read(2,10,end=1010) iN(i,1),(rAthresh(i,j),j=1,3),iN(i,2), : (rAthresh(i,j),j=4,6) 10 format(27x,2(i5,3f7.2)) 105 continue 1010 close (2) c --------------------------------------------------------------------- c Get existing overall tables into arrays open(30,file=cPath//'output/ave_app_temp_1day.ts',status='old') open(31,file=cPath//'output/max_app_temp_1day.ts',status='old') open(32,file=cPath//'output/min_app_temp_1day.ts',status='old') open(33,file=cPath//'output/ave_app_temp_4day.ts',status='old') open(34,file=cPath//'output/max_app_temp_4day.ts',status='old') open(35,file=cPath//'output/min_app_temp_4day.ts',status='old') do 500,i=1,6 c --------------------------------------------------------------------- c Initialize arrays do 501,j=1,500 iFyear(i,j) = -9 do 502,k=1940,2020 do 503,l=1,12 iDays(i,j,k,l) = -9 503 continue 502 continue 501 continue n = 1 lFirst = .true. c --------------------------------------------------------------------- c Read overall tables into arrays do 510,j=1,50000 read(29+i,31,end=5000) cStnid,iY,(iVal(k),k=1,12) 31 format(a5,1x,i4,1x,12i5) if (lFirst .eq. .true.) then iFyear(i,n) = iY lFirst = .false. endif if (cStnid .eq. cMeta(n)) then do 511,k=1,12 iDays(i,n,iY,k) = iVal(k) 511 continue else n = n + 1 if (cStnid .eq. cMeta(n)) then iFyear(i,n) = iY do 512,k=1,12 iDays(i,n,iY,k) = iVal(k) 512 continue else write(6,*) 'Error reading existing overall files ',i,' ',n stop endif endif 510 continue 5000 rewind 29+i 500 continue c --------------------------------------------------------------------- c Step through metadata array and grep data from ISD file do 110,i=1,iMeta * do 110,i=1,1 write(6,*) 'Processing '//cMeta(i),' ', : cStation(i)(1:iStation(i)),' ',cSt(i) call system('rm '//cPath//'programs/grep.dat') * call system('grep '//cMeta(i)//cY//cM//' '//cPath//'data/ISD'// * : cY//cM//'-V030.dat > '//cPath//'programs/grep.dat') call system('grep '//cMeta(i)//cY//cM//' '//cPath//'data/ISH'// : cMon(iMonth)//cY//'.dat > '//cPath//'programs/grep.dat') * call system('grep '//cMeta(i)//cSyr//cSmon//' '//cPath// * : 'data/ISD'//cSyr//cSmon//'-V030.dat >> '//cPath// * : 'programs/grep.dat') call system('grep '//cMeta(i)//cSyr//cSmon//' '//cPath// : 'data/ISH'//cMon(iSmon)//cSyr//'.dat >> '//cPath// : 'programs/grep.dat') c --------------------------------------------------------------------- c Check that data exists call system('rm '//cPath//'programs/size.dat') call system('wc -l '//cPath//'programs/grep.dat > '// : cPath//'programs/size.dat') open(2,file=cPath//'programs/size.dat',status='old') read(2,11) iSize 11 format(i8) close (2) * if (iSize .lt. 1) go to 110 c --------------------------------------------------------------------- c Open station heat index and apparent temperature files and read c through to end so that new data can be appended open(20,file=cPath//'app_temps/'//cStation(i)(1:iStation(i))// : '_'//cSt(i)//'.apptemp',status='unknown') do 310,j=1,10000 read(20,4,end=3100) cJunk 4 format(a23) read(cJunk(18:21),'(i4)') iPyear read(cJunk(22:23),'(i2)') iPmon if ((iPyear .eq. iYear) .and. (iPmon .eq. iMonth)) then backspace (20) go to 3100 endif 310 continue c --------------------------------------------------------------------- c Initialize count variables 3100 iMnum = 0 iYnum = 0 do 520,m=1,3 iMcount(m) = -9 iMhwave(m) = -9 iCount(m) = -9 iHwave(m) = -9 520 continue c --------------------------------------------------------------------- c Open and read existing station data files open(21,file=cPath//'output/'//cStation(i)(1:iStation(i))// : '_'//cSt(i)//'_1day_ave.ts',status='old') open(22,file=cPath//'output/'//cStation(i)(1:iStation(i))// : '_'//cSt(i)//'_1day_max.ts',status='old') open(23,file=cPath//'output/'//cStation(i)(1:iStation(i))// : '_'//cSt(i)//'_1day_min.ts',status='old') open(24,file=cPath//'output/'//cStation(i)(1:iStation(i))// : '_'//cSt(i)//'_4day_ave.ts',status='old') open(25,file=cPath//'output/'//cStation(i)(1:iStation(i))// : '_'//cSt(i)//'_4day_max.ts',status='old') open(26,file=cPath//'output/'//cStation(i)(1:iStation(i))// : '_'//cSt(i)//'_4day_min.ts',status='old') do 530,j=1,6 do 531,k=1,4 read(20+j,4,end=530) cJunk 531 continue do 532,k=iFyear(j,i),iYear-1 read(20+j,4,end=530) cJunk 532 continue if (iMonth .gt. 1) then read(20+j,51,end=530) (iLyear(j,k),k=1,12) 51 format(5x,12i5) else do 533,k=1,12 iLyear(j,k) = -9 533 continue endif rewind (20+j) 530 continue c --------------------------------------------------------------------- c Open grepped data file and initialize variables open(2,file=cPath//'programs/grep.dat',status='old') rAlt(i) = -9999 do 210,iD=1,31 do 211,iT=1,24 rPress(iD,iT) = -999.9 rDp(iD,iT) = -999.9 rTemp(iD,iT) = -999.9 rEtemp(iD,iT) = -999.9 rFtemp(iD,iT) = -999.9 iEtemp(iD,iT) = 99999 211 continue 210 continue iYnum = 0 do 139,k=1,6 iCount(k) = 0 139 continue do 138,k=1,3 i4count(k) = 0 iHwave(k) = 0 138 continue iMnum = 0 do 240,kk=1,3 iMcount(kk) = 0 iMhwave(kk) = 0 240 continue c --------------------------------------------------------------------- c Read grepped data into arrays of slp, RH and temperature do 111,j=1,50000 read(2,2,end=1002) iCh,cStID,iY,iM,iDay,iTime,cOtype, : iAlt,iTemp,iDp,iMslp 2 format(i4,6x,a5,i4,2i2,i4,14x,a5,i5,36x,3(i5,1x)) if (cOtype .ne. 'FM-15') go to 111 iD = iDay c --------------------------------------------------------------------- c Put observation into array of local time values if ((iAlt .lt. 9999) .and. (rAlt(i) .eq. -9999)) : rAlt(i) = float(iAlt) if (mod(iTime,100) .gt. 30) then iT = iTime/100 + 1 else iT = iTime/100 endif if ((iM .eq. iSmon) .and. : ((iDay .gt. 1) .or. : ((iDay .eq. 1) .and. (iT .gt. iTz(i))))) go to 1002 if ((iT .lt. 0) .or. (iT .gt. 24)) go to 111 if (iT .eq. 24) then iD = iD + 1 iT = 0 endif c --------------------------------------------------------------------- c Convert to local time iT = iT - iTz(i) if (iT .lt. 0) then iT = 24 + iT if (iM .eq. iSmon) then iD = iMlen else iD = iD - 1 endif endif iT = iT + 1 if ((iD .lt. 1) .or. (iD .gt. iMlen)) go to 111 c --------------------------------------------------------------------- c Write to monthly arrays if (rAlt(i) .eq. -9999) go to 111 if (iMslp .lt. 90000) then rPress(iD,iT) = iMslp/10. rPress(iD,iT) = rPress(iD,iT)/exp(rAlt(i)/8000.) rPress(iD,iT) = rPress(iD,iT)/10. endif if (iTemp .lt. 9000) then rTemp(iD,iT) = iTemp/10. endif if (iDp .lt. 9000) then rDp(iD,iT) = iDp/10. endif 111 continue 1002 close (2) c --------------------------------------------------------------------- c Step through monthly values and use subroutine to calculate hourly c apparent temperatures lWrite = .false. do 120,iD=1,iMlen iL(iD) = iD cMxfl(iD) = '0' cMnfl(iD) = '0' cMxS(iD) = ' ' cMnS(iD) = ' ' do 121,iT=1,24 if ((rPress(iD,iT) .gt. 0) .and. : (rTemp(iD,iT) .gt. -900) .and. : (rDp(iD,iT) .gt. -900)) then call etmp(rTemp(iD,iT),rDp(iD,iT),rPress(iD,iT), : rEtemp(iD,iT)) c --------------------------------------------------------------------- c Convert to Fahrenheit and write to a new array, then compile an array c of integer values rFtemp(iD,iT) = rEtemp(iD,iT)*9/5 + 32 endif if (rFtemp(iD,iT) .gt. -900) then if (rFtemp(iD,iT) .ge. 0) then if (rFtemp(iD,iT)-int(rFtemp(iD,iT)) .gt. 0.5) then iEtemp(iD,iT) = int(rFtemp(iD,iT)) + 1 else iEtemp(iD,iT) = int(rFtemp(iD,iT)) endif else if (rFtemp(iD,iT)-int(rFtemp(iD,iT)) .lt. : (-1)*0.5) then iEtemp(iD,iT) = int(rFtemp(iD,iT)) - 1 else iEtemp(iD,iT) = int(rFtemp(iD,iT)) endif endif endif 121 continue c --------------------------------------------------------------------- c Ascertain daily max and min values, and times of occurrence iEmax(iD) = -999 iMaxT(iD) = 99 iEmin(iD) = 999 iMinT(iD) = 99 rEmax = -999.9 rEmin = 999.9 rMax = -999.9 rMin = 999.9 do 122,iT=1,24 if (iEtemp(iD,iT) .ge. 900) go to 122 lWrite = .true. if (iEtemp(iD,iT) .gt. iEmax(iD)) then iEmax(iD) = iEtemp(iD,iT) iMaxT(iD) = iT - 1 endif if (iEtemp(iD,iT) .lt. iEmin(iD)) then iEmin(iD) = iEtemp(iD,iT) iMinT(iD) = iT - 1 endif if (rEtemp(iD,iT) .gt. rEmax) rEmax = rEtemp(iD,iT) if (rEtemp(iD,iT) .lt. rEmin) rEmin = rEtemp(iD,iT) if (rTemp(iD,iT) .gt. rMax) rMax = rTemp(iD,iT) if (rTemp(iD,iT) .lt. rMin) rMin = rTemp(iD,iT) 122 continue c --------------------------------------------------------------------- c Evaluate signs and QC indicators if (iEmax(iD) .eq. -999) then iEmax(iD) = 99999 cMxfl(iD) = '9' endif if (iEmin(iD) .eq. 999) then iEmin(iD) = 99999 cMnfl(iD) = '9' endif if (iEmax(iD) .lt. 0) then cMxS(iD) = '-' endif if (iEmin(iD) .lt. 0) then cMnS(iD) = '-' endif c --------------------------------------------------------------------- c Perform exceedence counts if ((rEmax .gt. -990) .and. (rEmin .lt. 990)) then iYnum = iYnum + 1 iMnum = iMnum + 1 rEave = (rEmax + rEmin)/2. else rEave = -999.9 rEmin = -999.9 endif rEval(iD,1) = rEave rEval(iD,2) = rEmax rEval(iD,3) = rEmin if ((rMax .gt. -990) .and. (rMin .lt. 990)) then rAve = (rMax + rMin)/2. else rAve = -999.9 endif rTval(iD,1) = rAve rTval(iD,2) = rMax rTval(iD,3) = rMin do 144,m=1,3 if (rEval(iD,m) .gt. rAthresh(i,m)) then iCount(m) = iCount(m) + 1 i4count(m) = i4count(m) + 1 iMcount(m) = iMcount(m) + 1 else if (i4count(m) .ge. 4) then iHwave(m) = iHwave(m) + 1 iMhwave(m) = iMhwave(m) + 1 endif i4count(m) = 0 endif if (rTval(iD,m) .gt. rAthresh(i,m+3)) : iCount(m+3) = iCount(m+3) + 1 144 continue 120 continue if (iMnum .eq. 0) then do 223,m=1,3 iMcount(m) = -9 iMhwave(m) = -9 223 continue endif if (iYnum .eq. 0) then do 150,m=1,3 iCount(m) = -9 iHwave(m) = -9 150 continue endif do 155,m=1,3 iLyear(m,iMonth) = iMcount(m) iLyear(m+3,iMonth) = iMhwave(m) 155 continue c --------------------------------------------------------------------- c Write past and new data to exceedence files do 320,iA=1,6 do 331,k=1,4 read(20+iA,4) cJunk 331 continue do 311,iB=iFyear(iA,i),iYear-1 read(20+iA,4) cJunk write(29+iA,31) cMeta(i),iB,(iDays(iA,i,iB,k),k=1,12) 311 continue write(29+iA,31) cMeta(i),iYear,(iLyear(iA,k),k=1,12) write(20+iA,35) iYear,(iLyear(iA,k),k=1,12) 35 format(i4,1x,12i5) 320 continue c --------------------------------------------------------------------- c Write to apparent temperature values to output if (lWrite .eq. .true.) then write(20,3) 'DLY000',cMeta(i),'ETMXF ',iYear,iMonth,'XXXX', : iMlen,(iL(l),iMaxT(l),cMxS(l),iEmax(l),cMxfl(l),l=1,iMlen) 3 format(a6,a5,a6,i4.4,i2.2,a4,i3.3,31(2i2.2,a1,i5.5,1x,a1)) write(20,3) 'DLY000',cMeta(i),'ETMNF ',iYear,iMonth,'XXXX', : iMlen,(iL(l),iMinT(l),cMnS(l),iEmin(l),cMnfl(l),l=1,iMlen) endif do 650,j=1,7 close (19+j) 650 continue c --------------------------------------------------------------------- c FTP to ftp0.ncdc.noaa.gov call system('rm '//cPath//'programs/ftpout') call system('echo "open ftp0.ncdc.noaa.gov" > '// : cPath//'programs/ftpout') call system('echo "user anonymous Trevor.Wallis@noaa.gov" >> '// : cPath//'programs/ftpout') call system('echo "cd /pub/data/heatstress/stationdata" >> '// : cPath//'programs/ftpout') call system('echo "lcd '//cPath//'output" >> '// : cPath//'programs/ftpout') call system('echo "delete '//cStation(i)(1:iStation(i))// : '_'//cSt(i)//'_1day_ave.ts" >> '//cPath//'programs/ftpout') call system('echo "put '//cStation(i)(1:iStation(i))// : '_'//cSt(i)//'_1day_ave.ts" >> '//cPath//'programs/ftpout') call system('echo "delete '//cStation(i)(1:iStation(i))// : '_'//cSt(i)//'_1day_max.ts" >> '//cPath//'programs/ftpout') call system('echo "put '//cStation(i)(1:iStation(i))// : '_'//cSt(i)//'_1day_max.ts" >> '//cPath//'programs/ftpout') call system('echo "delete '//cStation(i)(1:iStation(i))// : '_'//cSt(i)//'_1day_min.ts" >> '//cPath//'programs/ftpout') call system('echo "put '//cStation(i)(1:iStation(i))// : '_'//cSt(i)//'_1day_min.ts" >> '//cPath//'programs/ftpout') call system('echo "delete '//cStation(i)(1:iStation(i))// : '_'//cSt(i)//'_4day_ave.ts" >> '//cPath//'programs/ftpout') call system('echo "put '//cStation(i)(1:iStation(i))// : '_'//cSt(i)//'_4day_ave.ts" >> '//cPath//'programs/ftpout') call system('echo "delete '//cStation(i)(1:iStation(i))// : '_'//cSt(i)//'_4day_max.ts" >> '//cPath//'programs/ftpout') call system('echo "put '//cStation(i)(1:iStation(i))// : '_'//cSt(i)//'_4day_max.ts" >> '//cPath//'programs/ftpout') call system('echo "delete '//cStation(i)(1:iStation(i))// : '_'//cSt(i)//'_4day_min.ts" >> '//cPath//'programs/ftpout') call system('echo "put '//cStation(i)(1:iStation(i))// : '_'//cSt(i)//'_4day_min.ts" >> '//cPath//'programs/ftpout') call system('echo "lcd '//cPath//'app_temps" >> '// : cPath//'programs/ftpout') call system('echo "cd /pub/data/heatstress/app_temp" >> '// : cPath//'programs/ftpout') call system('echo "delete '//cStation(i)(1:iStation(i))// : '_'//cSt(i)//'.apptemp" >> '//cPath//'programs/ftpout') call system('echo "put '//cStation(i)(1:iStation(i))// : '_'//cSt(i)//'.apptemp" >> '//cPath//'programs/ftpout') call system('echo "quit" >> '//cPath//'programs/ftpout') call system('ftp -in < '//cPath//'programs/ftpout') 110 continue c --------------------------------------------------------------------- c Close overall threshold exceedance files and FTP to ftp0.ncdc.noaa.gov do 175,i=1,6 close (29+i) 175 continue call system('rm '//cPath//'programs/ftpout') call system('echo "open ftp0.ncdc.noaa.gov" > '// : cPath//'programs/ftpout') call system('echo "user anonymous Trevor.Wallis@noaa.gov" >> '// : cPath//'programs/ftpout') call system('echo "cd /pub/data/heatstress/stationdata/'// : 'monthly" >> '//cPath//'programs/ftpout') call system('echo "lcd '//cPath//'output" >> '// : cPath//'programs/ftpout') call system('echo "delete ave_app_temp_1day.ts" >> '// : cPath//'programs/ftpout') call system('echo "put ave_app_temp_1day.ts" >> '// : cPath//'programs/ftpout') call system('echo "delete ave_app_temp_4day.ts" >> '// : cPath//'programs/ftpout') call system('echo "put ave_app_temp_4day.ts" >> '// : cPath//'programs/ftpout') call system('echo "delete max_app_temp_1day.ts" >> '// : cPath//'programs/ftpout') call system('echo "put max_app_temp_1day.ts" >> '// : cPath//'programs/ftpout') call system('echo "delete max_app_temp_4day.ts" >> '// : cPath//'programs/ftpout') call system('echo "put max_app_temp_4day.ts" >> '// : cPath//'programs/ftpout') call system('echo "delete min_app_temp_1day.ts" >> '// : cPath//'programs/ftpout') call system('echo "put min_app_temp_1day.ts" >> '// : cPath//'programs/ftpout') call system('echo "delete min_app_temp_4day.ts" >> '// : cPath//'programs/ftpout') call system('echo "put min_app_temp_4day.ts" >> '// : cPath//'programs/ftpout') call system('echo "quit" >> '//cPath//'programs/ftpout') call system('ftp -in < '//cPath//'programs/ftpout') c --------------------------------------------------------------------- c Compress input data to avoid multiple processing * call system('compress '//cPath//'data/ISD'//cY//cM//'-V030.dat') call system('compress '//cPath//'data/ISH'//cMon(iMonth)//cY// : '.dat') c --------------------------------------------------------------------- c Quit program 1001 stop end *********************************************************************** c Dian Gaffen's subroutine for calculating apparent temperature c Modified to accept dew points rather than RH subroutine etmp(t,dp,p,ta) F= 1.0007 + P * 3.46E-6 Tk=t+273.15 ESAT = F * 6.1121 * EXP ((17.502*(Tk-273.15))/(Tk-32.18)) RH = exp((17.27*237.7*(dp-t))/((237.7+t)*(237.7+dp))) e = ESAT*RH EKPA = e/10. ta = -1.3+0.92*T+2.2*EKPA return end