program temp_ret c c program reads gathers postion information from the COORDINATES file c reads raw brightness temperature data aquired via buffer data on the c NCEP IBM, creates a TIMES file, and temperature retreivals for each of c the NOAA satellites (15,16,18). c c c INPUT: COORDINATES - This file contains a sequential list of c active storms their most current position c and their position 12-h earlier along with c heading, speed, intensity, and adeck storm c number and basin information. c BUFFER DATA - Three files will be generated for each storm c in the COORDINATES file with the following c naming convention: c noaa15.out c noaa16.out c noaa18.out c DUMP FILE - This file contains specific information about c the storm including the date time group, c the geographic location, and the ATCF storm c number. c COEFFICIENT FILES - Assorted files containing temperature c retrieval coefficents. **************************************************************************** c c OUTPUT: c Temperature Retreivals - Three files will be generated for each storm c in the COORDINATES file with the following c naming convention: c noaa15.ret c noaa16.ret c noaa18.ret c Log files for each c noaa15.log c noaa16.log c noaa18.log c c **************************************************************************** * * * Written By J. Knaff * CIRA * March 2003 * * Last Modified: December 13, 2005 * * History: * April 10, 2003, changed the documentation at the top * made modifications to M. Goldenbergs code * land is now coded as 0 (damn wierd!!!) * December 13, 2005, J.Knaff replaced NOAA17 with NOAA18 * **************************************************************************** parameter(ip=10000,nup=30) integer idate15(ip),itime15(ip),isat15(ip),isen15(ip), . iscan15(ip),ielev15(ip),isatz15(ip),ls15(ip) integer idate16(ip),itime16(ip),isat16(ip),isen16(ip), . iscan16(ip),ielev16(ip),isatz16(ip),ls16(ip) integer idate18(ip),itime18(ip),isat18(ip),isen18(ip), . iscan18(ip),ielev18(ip),isatz18(ip),ls18(ip) integer times(2),dates(2) integer inum,iy,im,id,i,dir0,spd0,vtmax,numatcf,vtmax12,irmw integer iy2,im2,id2,hh2,mm2,ss2,jd2,iy1,im1,id1,hh1,mm1,ss1 integer sattime,satdate,satjday,idt,imt,iyt,ijult c real rlat15(ip),rlon15(ip),satzen15(ip),solzen15(ip),tb15(ip,15), . temp15(ip,40),tpw15(ip),clw15(ip) real rlat16(ip),rlon16(ip),satzen16(ip),solzen16(ip),tb16(ip,15), . temp16(ip,40),tpw16(ip),clw16(ip) real rlat18(ip),rlon18(ip),satzen18(ip),solzen18(ip),tb18(ip,15), . temp18(ip,40),tpw18(ip),clw18(ip) real lats(2),lons(2) real closest,strmlatc,strmlonc,strmlat,strmlon, . satclat1,satclon1,satclat2,satclon2,satclat,satclon, . dist2cen,rad,tdiff,tpw,clw real bl(15),dummy(15),tret(40),bt(15) c character*6 stnum character*8 cdir character*10 stname character*12 dumpfile character*11 coorfile character*2 cbasin character*10 cdata15, cdata16, cdata18 character*10 coutp15, coutp16, coutp18 character*10 clogf15, clogf16, clogf18 c logical open15,open16,open18,opendump c c open dump files, one for each entry in the COORDINATES file c do i=1,10 write(dumpfile,50)i 50 format('dumpjb.in_',i2.2) opendump=.false. inquire(file=dumpfile,exist=opendump) if (opendump)then inumd=i open (1,file=dumpfile,status='old') read(1,51)idtg 51 format(i10) read(1,*)istuff read(1,52)cdir 52 format(a8) close (1) endif enddo c c c open COORDINATE file c coorfile='COORDINATES' open (10,file=coorfile,status='old') c c read information from the coordinate files and loop through storms c 1 read(10,30,err=888) inum,iy,im,id,it1,it2,times(2),lats(2), . lons(2),lats(1),lons(1),dir0,spd0,irmw,vtmax,vtmax12, . stnum,stname 30 format(i2.2,1x,i4.4,2i2.2,1x,i4.4,i3.3,i2.2,4f7.1,1x, . 5(i3.3,1x),a6,1x,a10) if(inum.eq.inumd)goto 31 goto 1 31 continue close(10) c c internal read to strip atcf information out of the coordinate files c which will be used to access the ascii buffer data. c read(stnum,45)cbasin,numatcf 45 format(a2,i2) c c c convert date to julian day c ijul=md2jul(id,im,iy) c c create dates(2) c dates(2)=iy*10000+im*100+id c c determine the t-12h position time and date c times(2)=times(2) ! nearest synoptic hour in this application c times(1)=times(2)-12 if (times(1).lt.0)then times(1)=24+times(1) ijult=ijul-1 call jdayi(ijult,iyt,imt,idt) dates(1)=iyt*10000+imt*100+idt else dates(1)=iy*10000+im*100+id endif c c convert hrs into hhmmss c times(1)=times(1)*10000 times(2)=times(2)*10000 c c construct the datafile names for NOAA - 15, 16, 18 etc.. c cdata15='noaa15.out' cdata16='noaa16.out' cdata18='noaa18.out' c c construct output files for this portion of the loop c coutp15='noaa15.ret' coutp16='noaa16.ret' coutp18='noaa18.ret' c c construct log file names for this date c clogf15='noaa15.log' clogf16='noaa16.log' clogf18='noaa18.log' c c open data(old),output(new),and log (new) files c inquire(file=cdata15,exist=open15) if(open15)then open(15,file=cdata15,status='old') open(25,file=coutp15,status='replace') open(35,file=clogf15,status='replace') endif inquire(file=cdata16,exist=open16) if(open16)then open(16,file=cdata16,status='old') open(26,file=coutp16,status='replace') open(36,file=clogf16,status='replace') endif inquire(file=cdata18,exist=open18) if(open18)then open(18,file=cdata18,status='old') open(28,file=coutp18,status='replace') open(38,file=clogf18,status='replace') endif c c read data NOAA 15 c if(open15)then m=1 115 read(15,75,end=915) idate15(m),itime15(m),rlat15(m), & rlon15(m),isat15(m),isen15(m),iscan15(m),ls15(m), & satzen15(m),solzen15(m),ielev15(m),isatz15(m), & (tb15(m,n),n=1,15) 75 format(i8,1x,i6,1x,2(f8.3,1x),1x,i2,2x,i3,2x,i2,2x,i1,2x, & f8.3,1x,f8.3,1x,i5,2x,i6,2x,15(f6.2,2x)) m=m+1 goto 115 c 915 continue if (m.eq.1)then write(35,835) cdata15,coutp15 835 format('Error encountered when reading file ',a10,/, & 'The file ',a10,' will be contain COORDINATE and ', & 'TIMES data, but no data') endif write(35,935) cdata15,m-1 935 format('data file ',a10,' read with ',i8,' data points') num15=m-1 closest=9999.9 c c loop through the data to find the closest satellite point to the storm c center. c c sattime=times(2) satdate=iy*10000+im*100+id satjday=ijul strmlatc=lats(2) strmlonc=lons(2) c do i=1,num15 c c Calculate time difference for the extrapolation/interpolation c iy2=dates(2)/10000 im2=(dates(2)-iy2*10000)/100 id2=(dates(2)-(iy2*10000+im2*100)) hh2=times(2)/10000 mm2=(times(2)-hh2*10000)/100 ss2=(times(2)-(hh2*10000+mm2*100)) jd2=md2jul(id2,im2,iy2) iy1=idate15(i)/10000 im1=(idate15(i)-iy1*10000)/100 id1=(idate15(i)-(iy1*10000+im1*100)) jd1=md2jul(id1,im1,iy1) hh1=itime15(i)/10000 mm1=(itime15(i)-hh1*10000)/100 ss1=(itime15(i)-(hh1*10000+mm1*100)) tdiff=(jd2-jd1)*24.+(hh2-hh1)+(mm2-mm1)/60.+ . (ss2-ss1)/3600. c c find the storm center c strmlat=lats(2) + tdiff*((lats(1)-lats(2))/12.) strmlon=lons(2) + tdiff*((lons(1)-lons(2))/12.) c c calculate distance from strom center to this point at the point relative c time. call distk(strmlon,strmlat,rlon15(i),rlat15(i),dx,dy,rad) c c save the closest point c if (rad.lt.closest)then closest=rad sattime=itime15(i) satdate=idate15(i) satjday=jd1 strmlatc=strmlat strmlonc=strmlon endif enddo write(35,945)strmlatc,strmlonc,satdate,sattime,closest 945 format('storm center located is at ',2f8.2,' at ',2i10,/, . 'The closest point to this center is ',f8.1,' km') c c one the satellite relative storm center is found, find the center of c the amsu swath that corresponds to the time of the observation. c satclat1=0.0 satclon1=0.0 satclat2=0.0 satclon2=0.0 do i=1,num15 if (itime15(i).eq.sattime)then if (iscan15(i).eq.15)then satclat1=rlat15(i) satclon1=rlon15(i) endif if (iscan15(i).eq.16)then satclat2=rlat15(i) satclon2=rlon15(i) endif endif enddo satclat=(satclat1+satclat2)/2. satclon=(satclon1+satclon2)/2. call distk(strmlon,strmlat,satclon,satclat,dx,dy,rad) dist2cen=rad write(35,955)satclat,satclon,sattime,dist2cen 955 format('satellite center is located at ',2f8.2,' at ',i6,/, . 'This point is ',f8.1,' km from the closest point') c c writing coordinate and times information to output file and log file c write(35,965) 965 format('writing COORDINATE and TIMES information') write(25, 30)inum,iy,im,id,it1,it2,times(2)/10000 . ,lats(2),lons(2),lats(1),lons(1),dir0,spd0,irmw, . vtmax,vtmax12,stnum,stname write(25,975)inum,satdate/10000,satjday,sattime, . satclat,satclon,strmlatc,strmlonc 975 format(i2.2,1x,i4.4,i3.3,1x,i6.6,4f9.2) write(25,985)num15 985 format(i6) c c call temperature retrieval c if (num15.ne.0) write(35,995) 995 format('starting NOAA 15 temperature retieval') c do i=1, num15 iflag=0 do j=1,15 dummy(j)=tb15(i,j) if (j.ne.11.and.j.ne.14.and.dummy(j).lt.0.0)iflag=1 enddo call amsua_ret_n15(nup,dummy,bt,iscan15(i),ls15(i),bl, & tret,tpw,clw) do j=1,40 if (iflag.eq.0)then temp15(i,j)=tret(j) else temp15(i,j)=-99.00 endif enddo if (iflag.eq.0)then tpw15(i)=tpw clw15(i)=clw else tpw15(i)=-99.00 clw15(i)=-99.00 endif c c write NOAA 15 ascii data with temperature retreivals, total c precipitable water, and cloud liquid water tacked on the end. c write(25,1575) idate15(i),itime15(i),rlat15(i), & rlon15(i),isat15(i),isen15(i),iscan15(i),ls15(i), & satzen15(i),solzen15(i),ielev15(i),isatz15(i), & (tb15(i,n),n=1,15),(temp15(i,n),n=1,40),tpw15(i), & clw15(i) 1575 format(i8,1x,i6,1x,2(f8.3,1x),1x,i2,2x,i3,2x,i2,2x,i1,2x, & f8.3,1x,f8.3,1x,i5,2x,i6,2x,57(f6.2,2x)) enddo c c close NOAA 15 output file for coordinate file inum if (num15.ne.0)write(35,990) 990 format('Normal completion of NOAA 15 temperature Retrieval') 1015 continue call flush(15) close(15) call flush(25) close(25) call flush(35) close(35) endif c c read data NOAA 16 c if(open16)then m=1 116 read(16,75,end=916) idate16(m),itime16(m),rlat16(m), & rlon16(m),isat16(m),isen16(m),iscan16(m),ls16(m), & satzen16(m),solzen16(m),ielev16(m),isatz16(m), & (tb16(m,n),n=1,15) m=m+1 goto 116 916 continue if(m.eq.1) then write(36,835) cdata16,coutp16 endif write(36,935) cdata16,m-1 num16=m-1 closest=9999.9 c c loop through the data to find the closest satellite point to the storm c center. c c sattime=times(2) satdate=iy*10000+im*100+id satjday=ijul strmlatc=lats(2) strmlonc=lons(2) c do i=1,num16 c c Calculate time difference for the extrapolation/interpolation c iy2=dates(2)/10000 im2=(dates(2)-iy2*10000)/100 id2=(dates(2)-(iy2*10000+im2*100)) hh2=times(2)/10000 mm2=(times(2)-hh2*10000)/100 ss2=(times(2)-(hh2*10000+mm2*100)) jd2=md2jul(id2,im2,iy2) iy1=idate16(i)/10000 im1=(idate16(i)-iy1*10000)/100 id1=(idate16(i)-(iy1*10000+im1*100)) jd1=md2jul(id1,im1,iy1) hh1=itime16(i)/10000 mm1=(itime16(i)-hh1*10000)/100 ss1=(itime16(i)-(hh1*10000+mm1*100)) tdiff=(jd2-jd1)*24.+(hh2-hh1)+(mm2-mm1)/60.+ . (ss2-ss1)/3600. c c find the storm center c strmlat=lats(2) + tdiff*((lats(1)-lats(2))/12.) strmlon=lons(2) + tdiff*((lons(1)-lons(2))/12.) c c calculate distance from strom center to this point at the point relative c time. call distk(strmlon,strmlat,rlon16(i),rlat16(i),dx,dy,rad) c c save the closest point c if (rad.lt.closest)then closest=rad sattime=itime16(i) satdate=idate16(i) satjday=jd1 strmlatc=strmlat strmlonc=strmlon endif enddo write(36,945)strmlatc,strmlonc,satdate,sattime,closest c c one the satellite relative storm center is found, find the center of c the amsu swath that corresponds to the time of the observation. c satclat1=0.0 satclon1=0.0 satclat2=0.0 satclon2=0.0 do i=1,num16 if (itime16(i).eq.sattime)then if (iscan16(i).eq.16)then satclat1=rlat16(i) satclon1=rlon16(i) endif if (iscan16(i).eq.16)then satclat2=rlat16(i) satclon2=rlon16(i) endif endif enddo satclat=(satclat1+satclat2)/2. satclon=(satclon1+satclon2)/2. call distk(strmlon,strmlat,satclon,satclat,dx,dy,rad) dist2cen=rad write(36,955)satclat,satclon,sattime,dist2cen c c writing COORDINATE, and timesinformation to output file and log file c write(36,965) write(26, 30)inum,iy,im,id,it1,it2,times(2)/10000 . ,lats(2),lons(2),lats(1),lons(1),dir0,spd0,irmw, . vtmax,vtmax12,stnum,stname write(26,975)inum,satdate/10000,satjday,sattime, . satclat,satclon,strmlatc,strmlonc write(26,985)num16 c c call temperature retrieval c if(num16.ne.0)write(36,996) 996 format('starting NOAA 16 temperature retieval') do i=1, num16 iflag=0 do j=1,15 dummy(j)=tb16(i,j) if (dummy(j).lt.0.0)iflag=1 enddo call amsua_ret_n16(nup,dummy,bt,iscan16(i),ls16(i),bl, & tret,tpw,clw) do j=1,40 if (iflag.eq.0)then temp16(i,j)=tret(j) else temp16(i,j)=-99.00 endif enddo if (iflag.eq.0)then tpw16(i)=tpw clw16(i)=clw else tpw16(i)=-99.00 clw16(i)=-99.00 endif c c write NOAA 16 ascii data with temperature retreivals, total c precipitable water, and cloud liquid water tacked on the end. c write(26,1575) idate16(i),itime16(i),rlat16(i), & rlon16(i),isat16(i),isen16(i),iscan16(i),ls16(i), & satzen16(i),solzen16(i),ielev16(i),isatz16(i), & (tb16(i,n),n=1,15),(temp16(i,n),n=1,40),tpw16(i), & clw16(i) enddo c c close NOAA 16 output file for coordinate file inum if(num16.ne.0)write(36,991) 991 format('Normal completion of NOAA 16 temperature Retrieval') 1016 continue call flush(16) close(16) call flush(26) close(26) call flush(36) close(36) endif c c read data NOAA 18 c if(open18)then m=1 118 read(18,75,end=918) idate18(m),itime18(m),rlat18(m), & rlon18(m),isat18(m),isen18(m),iscan18(m),ls18(m), & satzen18(m),solzen18(m),ielev18(m),isatz18(m), & (tb18(m,n),n=1,15) m=m+1 goto 118 918 continue if(m.eq.1) then write(38,835) cdata18,coutp18 endif write(38,935) cdata18,m-1 num18=m-1 closest=9999.9 c c loop through the data to find the closest satellite point to the storm c center. c c c sattime=times(2) satdate=iy*10000+im*100+id satjday=ijul strmlatc=lats(2) strmlonc=lons(2) c do i=1,num18 c c Calculate time difference for the extrapolation/interpolation c iy2=dates(2)/10000 im2=(dates(2)-iy2*10000)/100 id2=(dates(2)-(iy2*10000+im2*100)) hh2=times(2)/10000 mm2=(times(2)-hh2*10000)/100 ss2=(times(2)-(hh2*10000+mm2*100)) jd2=md2jul(id2,im2,iy2) iy1=idate18(i)/10000 im1=(idate18(i)-iy1*10000)/100 id1=(idate18(i)-(iy1*10000+im1*100)) jd1=md2jul(id1,im1,iy1) hh1=itime18(i)/10000 mm1=(itime18(i)-hh1*10000)/100 ss1=(itime18(i)-(hh1*10000+mm1*100)) tdiff=(jd2-jd1)*24.+(hh2-hh1)+(mm2-mm1)/60.+ . (ss2-ss1)/3600. c c find the storm center c strmlat=lats(2) + tdiff*((lats(1)-lats(2))/12.) strmlon=lons(2) + tdiff*((lons(1)-lons(2))/12.) c c calculate distance from strom center to this point at the point relative c time. call distk(strmlon,strmlat,rlon18(i),rlat18(i),dx,dy,rad) c c save the closest point c if (rad.lt.closest)then closest=rad sattime=itime18(i) satdate=idate18(i) satjday=jd1 strmlatc=strmlat strmlonc=strmlon endif enddo write(38,945)strmlatc,strmlonc,satdate,sattime,closest c c one the satellite relative storm center is found, find the center of c the amsu swath that corresponds to the time of the observation. c satclat1=0.0 satclon1=0.0 satclat2=0.0 satclon2=0.0 do i=1,num18 if (itime18(i).eq.sattime)then if (iscan18(i).eq.18)then satclat1=rlat18(i) satclon1=rlon18(i) endif if (iscan18(i).eq.18)then satclat2=rlat18(i) satclon2=rlon18(i) endif endif enddo satclat=(satclat1+satclat2)/2. satclon=(satclon1+satclon2)/2. call distk(strmlon,strmlat,satclon,satclat,dx,dy,rad) dist2cen=rad write(38,955)satclat,satclon,sattime,dist2cen c c writing COORDINATE, TIMES, number of points information to output file c write(38,965) write(28, 30)inum,iy,im,id,it1,it2,times(2)/10000 . ,lats(2),lons(2),lats(1),lons(1),dir0,spd0,irmw, . vtmax,vtmax12,stnum,stname write(28,975)inum,satdate/10000,satjday,sattime, . satclat,satclon,strmlatc,strmlonc write(28,985)num18 c c call temperature retrieval c if(num18.ne.0)write(38,997) 997 format('starting NOAA 18 temperature retieval') do i=1, num18 iflag=0 do j=1,15 dummy(j)=tb18(i,j) if (j.ne.7.and.dummy(j).lt.0.0)iflag=1 enddo call amsua_ret_n18(nup,dummy,bt,iscan18(i),ls18(i),bl, & tret,tpw,clw) do j=1,40 if (iflag.eq.0)then temp18(i,j)=tret(j) else temp18(i,j)=-99.00 endif enddo if (iflag.eq.0)then tpw18(i)=tpw clw18(i)=clw else tpw18(i)=-99.00 clw18(i)=-99.00 endif c c write NOAA 18 ascii data with temperature retreivals, total c precipitable water, and cloud liquid water tacked on the end. c write(28,1575) idate18(i),itime18(i),rlat18(i), & rlon18(i),isat18(i),isen18(i),iscan18(i),ls18(i), & satzen18(i),solzen18(i),ielev18(i),isatz18(i), & (tb18(i,n),n=1,15),(temp18(i,n),n=1,40),tpw18(i), & clw18(i) enddo c c close NOAA 18 output file for coordinate file inum if(num18.ne.0)write(38,992) 992 format('Normal completion of NOAA 18 temperature Retrieval') 1018 continue call flush(18) close(18) call flush(28) close(28) call flush(38) close(38) endif goto 999 888 open (28,file='ret.log') write(28,889)coorfile close(28) 889 format('END OF ',a22,' REACHED PREMATURELY EXICUTION STOPPED') 999 stop end