program ilsin c This program makes the input file for the lsdiag.f c program. The input file is for the current storm c with stormcard and track files istormcard.dat and c itrack.dat. The itrack.dat file should have the format c of the new ATCF. c c Modified Jul 2002 (MD) to add Ocean Heat Content line c Modified Aug 2000 (MD) to make input file for IR program on the IBM c Modified May 2005 (MD) to switch from 12 to 6 hour intervals c Modified Aug 2006 (MD) to fix dateline crossing problem c Modified Apr 2009 (MD) to add proxy OHC line for use when RHCN is missing c and add east Pacific OHC calculation c Modified May 2012 (MD) to accomodate NCODA OHC files c parameter (mft=20) dimension rlat(-2:mft),rlon(-2:mft) dimension ilat(-2:mft),ilon(-2:mft) dimension inlat(0:mft),inlon(0:mft),intime(-2:mft) dimension sst(0:mft), rsst(0:mft), ohc(0:mft), dist(0:mft) dimension isst(0:mft),irsst(0:mft),iohc(0:mft),idist(0:mft) dimension ipohc(0:mft) dimension iadum99(0:mft) c character *1 cbasin,satlab1,ccom,cbla,iewl character *6 cbasin6 character *4 sname4,vlab,satlab character *10 sname character *25 fnis,fnit,fnls,fnir c data irdelt /9999/ data inlat /21*0/ data inlon /21*0/ c data intime /-24,-12,0,12,24,36,48,60,72,84,96,108,120/ data intime /-12,-6,0, 6,12,18,24,30,36, 42, 48, 54, 60, + 66,72,78,84,90,96,102,108,114,120/ data sst /21*999.9/ data rlat /23*999.9/ data rlon /23*999.9/ data iadum99 /21*9999/ data idum99 /9999/ data idum00 /0/ data dtr /0.0174533/ c data luis,luit,luls,luir /30,31,32,33/ c rmiss = 999.9 imiss = 9999 c c Open the input and output files fnis = 'istormcard.dat' fnit = 'itrack.dat' fnls = 'lsdiag.inp' fnir = 'INFILE.in' c open(unit=luis,file=fnis,form='formatted',status='old',err=900) open(unit=luit,file=fnit,form='formatted',status='old',err=900) open(unit=luls,file=fnls,form='formatted',status='new',err=900) c c Read the stormcard file read(luis,100) isyr,ismon,isday read(luis,102) istime read(luis,103) ilat0,ilatm12 read(luis,103) ilon0,ilonm12 read(luis,104) read(luis,104) read(luis,105) ivmx0,ivmxm12 read(luis,106) ihead read(luis,106) ispeed read(luis,107) sname read(luis,108) cbasin6 100 format(3(i2)) 102 format(i2) 103 format(i4,1x,i4) 104 format(a72) 105 format(i3,2x,i3) 106 format(i3) 107 format(1x,a10) 108 format(a6) cbasin(1:1) = cbasin6(1:1) c c Read the forecast positions and put them in their proper place idelt = 6 1000 continue read(luit,110,end=1001) intimet,inlatt,inlont,iewl 110 format(29x,i4,1x,i4,2x,i5,a1) indx = intimet/6 if (indx .le. mft .and. indx .ge. 0) then c Check longitude for deg E if (iewl .eq. 'E') then inlont = 1800 + (1800-inlont) if (inlont .gt. 3000) inlont=inlont-3600 endif inlat(indx) = inlatt inlon(indx) = inlont endif go to 1000 1001 continue c do 10 k=1,mft rlat(k) = 0.1*float(inlat(k)) rlon(k) = 0.1*float(inlon(k)) 10 continue c c Check for missing lat,lon do 12 k=1,mft if (rlat(k) .le. 0.0 .or. rlon(k) .le. 0.0) then rlat(k) = rmiss rlon(k) = rmiss endif 12 continue c c Get t=0 lat,lon values from the stormcard input. rlat(0) = 0.1*float(ilat0) rlon(0) = 0.1*float(ilon0) c c Convert storm speed and heading into 12 hr lat/lon displacements headr = dtr*float(ihead) dislat = 12.0*float(ispeed)/60.0 dlat = dislat*cos(headr) dlon = -dislat*sin(headr)/cos(rlat(0)*dtr) c c Put t=-12 h lat,lon from stormcard in rlat,rlon arrays. c If the t=-12 h lat,lon are missing, estimate them from c the t=0 position and storm motion vector. im12 = -12/idelt if (ilatm12 .gt. 0 .and. ilonm12 .gt. 0) then rlat(im12) = 0.1*float(ilatm12) rlon(im12) = 0.1*float(ilonm12) else rlat(im12) = rlat(0) - dlat rlon(im12) = rlon(0) - dlon endif c c Linearly interpolate to fill in missing lat/lons c Set max number of times to search for a good value ismax=3 do i=-1,mft if (rlat(i) .gt. 900.0 .or. rlon(i) .gt. 900.0) then istr = i-ismax iend = i+ismax if (istr .lt. -2) istr= -2 if (iend .gt. mft) iend=mft c c Search for nearest good value before current time tbefore = 0.0 do k=i-1,istr,-1 if (rlat(k) .lt. rmiss .and. + rlon(k) .lt. rmiss ) then tbefore = float(k-i) rlatb = rlat(k) rlonb = rlon(k) go to 5000 endif enddo 5000 continue c c Search for nearest good value after current time tafter = 0.0 do k=i+1,iend if (rlat(k) .lt. rmiss .and. + rlon(k) .lt. rmiss ) then tafter = float(k-i) rlata = rlat(k) rlona = rlon(k) go to 5001 endif enddo 5001 continue c if (tbefore .lt. 0.0 .and. + tafter .gt. 0.0 ) then c c Replace missing value with linear interpolation deltat = tafter-tbefore wtb = tafter/deltat wta = -tbefore/deltat rlat(i) = wtb*rlatb + wta*rlata rlon(i) = wtb*rlonb + wta*rlona endif endif enddo c c Calculate the sea surface temperature, ocean heat content and distance c to land along the storm track do 15 i=0,mft if (rlat(i) .lt. 500.0 .and. rlon(i) .lt. 500.0) then rlatt = rlat(i) rlont = -rlon(i) c call aland(rlont,rlatt,distt) if (distt .gt. 5000.0) distt=5000.0 if (distt .lt. -999.0) distt=-999.0 c if (cbasin .eq. 'A') then call asst(rlont,rlatt,ismon,isday,sstt,dmlt) else call esst(rlont,rlatt,ismon,isday,sstt,dmlt) endif if (sstt .gt. 50.0) sstt=999.9 c call rysst(rlatt,rlont,rsstt,iryr,irmon,irday,irtime,ierr) if (ierr .ne. 0) rsstt=999.9 c c After the first call to rysst, calculate the time lag between c the stormcard date and the SST data file date. if (i .eq. 0 .and. ierr .eq. 0) then call tdiff(isyr,ismon,isday,istime, + iryr,irmon,irday,irtime,irdelt) endif c if (cbasin .eq. 'A' .or. cbasin .eq. 'E' + .or. cbasin .eq. 'C') then call ryohc(rlatt,rlont,ohct,ioyear,iomon,ioday,iotime, + ierr) c print *, "OHC-ryohc", rlatt,rlont,ohct,ierr if (ierr .ne. 0) ohct=9999. if (distt .lt. -20.0) ohct=9999. c c After the first call to ryohc, calculate the time lag between c the stormcard date and the OHC data file date. if (i .eq. 0) then if (ierr .eq. 0) then ioyr = ioyear - 100*(ioyear/100) call tdiff(isyr,ismon,isday,istime, + ioyr,iomon,ioday,iotime,iodelt) c else iodelt=9999 endif endif else ohct=9999. iodelt=9999 endif c c Calculate the proxy OHC data if (isyr .le. 50) then iyr4 = 2000 + isyr else iyr4 = 1900 + isyr endif c rlontw = -1.0*rlont call ohcparam(iyr4,ismon,isday,rlatt,rlontw, + rsstt,rmiss,pohc,ierr) c if (ierr .ne. 0) then ipohc(i) = imiss else ipohc(i) = ifix(pohc+0.5) endif if (ipohc(i) .gt. 200 .or. ipohc(i) .lt. 0) ipohc(i) = imiss c sst(i) = sstt rsst(i) = rsstt ohc(i) = ohct dist(i) = distt else sst(i) = 999.9 rsst(i) = 999.9 ohc(i) = 9999. dist(i) = 9999. ipohc(i)= imiss endif isst(i) = ifix(10.0*sst(i)) irsst(i) = ifix(10.0*rsst(i)) iohc(i) = ifix(ohc(i)) idist(i) = ifix(dist(i)) 15 continue c do 20 k=-2,mft ilat(k) = ifix(0.5 + 10.0*rlat(k)) ilon(k) = ifix(0.5 + 10.0*rlon(k)) 20 continue c c Write the lsdiag.inp file c sname4 = sname(1:4) sname4 = cbasin6(1:4) vlab = 'HEAD' write(luls,200) sname4,isyr,ismon,isday,istime,ivmx0,vlab 200 format(1x,a4,1x,3(i2.2),1x,i2.2,2x,i3,96x,a4) c vlab = 'TIME' write(luls,201) (intime(k),k=-2,mft),vlab 201 format(23(1x,i4),1x,a4) c if (ivmx0 .gt. 19 .and. ivmxm12 .gt. 19) then idv12 = ivmx0-ivmxm12 else idv12 = 0 endif c vlab='DELV' write(luls,205) idv12,idum99,(iadum99(k),k=0,mft),vlab 205 format(23(1x,i4),1x,a4) c vlab='LAT ' write(luls,215) (ilat(k),k=-2,mft),vlab vlab='LON ' write(luls,215) (ilon(k),k=-2,mft),vlab 215 format(23(1x,i4),1x,a4) c vlab='CSST' write(luls,220) (isst(k) ,k=0,mft),vlab vlab='DTL ' write(luls,220) (idist(k),k=0,mft),vlab 220 format(10x,21(1x,i4),1x,a4) c vlab='RSST' if (irdelt .gt. 9999) irdelt=9999 write(luls,221) (irsst(k),k=0,mft),vlab,irdelt 221 format(10x,21(1x,i4),1x,a4,1x,i4) c c Check RSST and RHCN values for consistency do k=0,mft write(6,*) k,irsst(k) if (irsst(k) .gt. 260 .and. irsst(k) .lt. 350) then if (iohc(k) .eq. 0) iohc(k) = 9999 endif enddo c vlab='RHCN' if (iodelt .gt. 9999) iodelt=9999 write(luls,221) (iohc(k),k=0,mft),vlab,iodelt c vlab='PHCN' write(luls,220) (ipohc(k),k=0,mft),vlab c vlab='LAST' write(luls,220) (iadum99(k),k=0,mft),vlab c c Make the input file for getting the IR profiles on the IBM open(unit=luir,file=fnir,form='formatted',status='replace', + err=900) c c Write the comment line if (cbasin .eq. 'A') then satlab='EAST' ist1 = 15 ist2 = 45 elseif (cbasin .eq. 'E') then satlab='WEST' ist1 = 0 ist2 = 30 else satlab='WEST' ist1 = 0 ist2 = 30 endif satlab1(1:1) = satlab(1:1) c c Get storm position now and three hours previous slatn = rlat(0) slonn = -rlon(0) c slatp = 0.75*rlat(0) + 0.25*rlat(-1) slonp = 0.75*rlon(0) + 0.25*rlon(-1) slonp = -slonp c c Calculate the hour 1-5 hours ago istime1 = istime-1 istime2 = istime-2 istime3 = istime-3 istime4 = istime-4 istime5 = istime-5 c if (istime1 .lt. 0) istime1 = istime1 + 24 if (istime2 .lt. 0) istime2 = istime2 + 24 if (istime3 .lt. 0) istime3 = istime3 + 24 if (istime4 .lt. 0) istime4 = istime4 + 24 if (istime5 .lt. 0) istime5 = istime5 + 24 c c Write the header line on the file for the IR data write(luir,300,err=900) cbasin6,istime,sname 300 format(a6,1x,i2.2,1x,a10) c c Write the satellite name and generic name for IR c profile output files write(luir,310) satlab 310 format(a4,/,'IRRP') c c Write the current time, storm position and imagery choices ccom = ',' cbla = ' ' c write(luir,320) istime,slatn,slonn 320 format(i2.2,'Z',7x,f5.1,4x,f6.1) if (cbasin .eq. 'A') then write(luir,330) satlab1,istime1,ist1,ccom, + satlab1,istime1,ist2,ccom, + satlab1,istime2,ist2,ccom, + satlab1,istime ,ist1,cbla else write(luir,330) satlab1,istime1,ist2,ccom, + satlab1,istime1,ist1,ccom, + satlab1,istime2,ist2,ccom, + satlab1,istime ,ist1,cbla endif 330 format(4(a1,i2.2,i2.2,'ZI4',a1)) c c Write the previous time, storm position and imagery choices (t=-3 hr) write(luir,320) istime3,slatp,slonp if (cbasin .eq. 'A') then write(luir,330) satlab1,istime3,ist1,ccom, + satlab1,istime4,ist2,ccom, + satlab1,istime3,ist2,ccom, + satlab1,istime4,ist1,cbla else write(luir,330) satlab1,istime3,ist1,ccom, + satlab1,istime4,ist2,ccom, + satlab1,istime3,ist2,ccom, + satlab1,istime4,ist1,cbla endif c close(luir) c 900 stop end