program ilsin c This program makes the input file for the lsdiag.f program. c The input file is for the current storm with c stormcard and track files called istormcard.dat and c ntrack.dat. The ntrack.dat file should have the format c of the ATCF. c c All of the SST and OHC processing is also included in this routine. c c Ver 2.5.1 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 Dec 2011 (MD) to add west Pacific option c Modified May 2012 (MD) to accomodate NCODA OHC files c Modified May 2013 (MD) for new input required by WCOSS IR readers c Modified May 2013 (MD) to add climo OHC, D20 and D26 variables, c and replace a,e,wsst routines with global rsstclim routine c Modified Aug 2013 (MD) for global land routine c Modified Sep 2013 (MD) to set sign of lat/lon based on new stormcard input variables c instead of basin ID. Climo OHC, D20 and D26 also added. c Modified Aug 2014 (MD) for WCOSS implementation c Modified May 2016 (MD) for ver2.0.0 to add t=-1.5 hr IR input file name c and head/speed to file for creating IR predictors c Modified Nov 2016 (MD) to use ntrack.dat instead of itrack.dat c (part of clean up for WCOSS Cray conversion). c ilsin.log file added. (Ver 2.1.0) c Modified Feb 2018 (MD) to add DSST,DSTA processing (Ver 2.2.0 ) c phcn variable removed. c Remaining ifix functions changed to nint. c Modified May 2018 (SNS) removed RHCN from lsdiag output (new: NOHC) c added 16 more climo ocean variables c changed satellite logic - based on lon rather than basin c Modified Jul 2018 (SNS) added imiss values to ncoda variables for missing lat/lons c Modified Dec 2018 (MD) for ver2.3.0. c 1)Call to ryohc subroutine removed c (old OHC files processing). c 2)PHCN removed. c 3)Time lag for new OHC files added c 4)Check to set new OHC data to zero for c inland points removed c Modified Jan 2019 (GC) for ver2.3.1. c - added option to use updated ncodaclimo2 c - that can process both 0.25 deg and 0.5 c - deg climatology c Modified Apr 2019 (GC) for ver2.4.0 c - removed all ncodaclimo2 and eohc c processing. use eohcaddd instead. c Modified May 2019 (SS) for ver2.5.0 c - added GOES-17 processing, including new variable (loncutoff) c set at top to define IR cutoff c Modified Apr 2020 (SS) for 7-day forecasts c Modified May 2020 (SS) for setting lon near PACK boundaries to missing values c Modified Jan 2023 (KM) removed go to statements c Modified Feb 2023 (MD) to skip call to rysst for 2023 and later c c Use statement for new ncodaclimo2 routine c use ncoda_climo2,only: n_data_arr_max,ncodaclimo2 c parameter (mft=28) 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 dsst(0:mft),dssta(0:mft) dimension idsst(0:mft),idssta(0:mft) dimension id20c(-2:mft),id26c(-2:mft),iohcc(-2:mft),ipohc(0:mft) dimension iadum99(0:mft) c real rlonta c character *1 latNS,lonEW,lonlab character *1 cbasin,satlab1,ccom,cbla,iewl,insl character *6 cbasin6 character *8 cbasin8 character *4 sname4,vlab character *20 satlab character *10 sname character *25 fnis,fnit,fnls,fnir,fnlg character *55 hfile c c data irdelt /9999/ data inlat /29*0/ data inlon /29*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, + 126,132,138,144,150,156,162,168/ data sst /29*999.9/ data rlat /31*999.9/ data rlon /31*999.9/ data iadum99 /29*9999/ data idum99 /9999/ data idum00 /0/ data dtr /0.0174533/ c data luis,luit,luls,luir,lulg /30,31,32,33,34/ c rmiss = 999.9 imiss = 9999 c c Set cutoff longitude for G16 v G17 IR c loncutoff must include (-) if degW loncutoff = -130.0 loncut = loncutoff + 360.0 c c Open the input and output files fnis = 'istormcard.dat' fnit = 'ntrack.dat' fnls = 'lsdiag.inp' fnir = 'INFILE.in' fnlg = 'ilsin.log' c open(unit=lulg,file=fnlg,form='formatted',status='replace') 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='replace', + err=900) c c Read the stormcard file read(luis,100) isyr,ismon,isday read(luis,102) istime read(luis,103) ilat0,ilatm12,latNS read(luis,103) ilon0,ilonm12,lonEW 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,1x,a1) 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 Calculate 4 digit year if (isyr .le. 50) then iyr4 = 2000 + isyr else iyr4 = 1900 + isyr endif c c Make lat negative for SH if (latNS .eq. 'S') then ilat0 = -iabs(ilat0) ilatm12 = -iabs(ilatm12) endif c c Calculate analysis time storm location, with lon in 0-360 convention alat = float(ilat0)/10.0 alon = float(ilon0)/10.0 c c Create 8 character ATCF ID cbasin8(1:4) = cbasin6(1:4) cbasin8(7:8) = cbasin6(5:6) c read(cbasin6(5:6),150) iyrstm 150 format(i2) if (iyrstm .ge. 50) then cbasin8(5:6) = '19' iyear4 = iyrstm + 1900 else cbasin8(5:6) = '20' iyear4 = iyrstm + 2000 endif c write(lulg,500) iyear4,ismon,isday,istime,cbasin8 500 format('Start ilsin for first part of SHIPS diagnostic file',/, + 'Initialization date/time: ',i4,1x,i2.2,i2.2, + 1x,i2.2,' UTC for ',a8) c write(lulg,504) ivmx0,alat,latNS,alon,lonEW 504 format(5x,'t= 0 vmax, lat, lon: ',i3,1x,f5.1,1x,a1,1x, + f6.1,1x,a1) c if (lonEW .eq. 'W') then alon = 360.0-alon endif c c Determine if longitude is deg W negative (iloncon=1) or deg E positive (iloncon=2) c based on the initial position. Deg W neg is for storms in the AL, EP or CP domains even c if the storm originated in a different basin. Similarly, Deg E pos is for storms in the c WP, IO or SH even if the storm originated in a different basin. if (alat .lt. 0.0) then iloncon = 2 else if (alon .ge. 180.0 .or. alon .lt. 15.0) then iloncon = 1 else iloncon = 2 endif endif c c Read the forecast positions and put them in their proper place. c Use 0-360 convention for longitude for now. idelt = 6 do read(luit,110,end=1001) intimet,inlatt,insl,inlont,iewl 110 format(29x,i4,1x,i4,a1,1x,i5,a1) indx = intimet/6 if (indx .le. mft .and. indx .ge. 0) then c Check longitude for deg W if (iewl .eq. 'W') then inlont = 3600 - inlont endif c c Check latitude for deg S if (insl .eq. 'S') then inlatt = -1*iabs(inlatt) endif c inlat(indx) = inlatt inlon(indx) = inlont endif enddo 1001 continue c c Set the basin flag by the agency flag and the analysis lat/lon c ibasin = 1, 2, 2, 3, 4, 5 for c AL, EP, CP, WP, IO, SH c if (alat .lt. 0.0) then ibasin = 5 else if (alon .ge. 15.0 .and. alon .lt. 100.0) then ibasin = 4 elseif (alon .ge. 100.0 .and. alon .lt. 180.0) then ibasin = 3 else ibasin = 2 if (cbasin .eq. 'A') ibasin = 1 endif endif c do k=1,mft rlat(k) = 0.1*float(inlat(k)) rlon(k) = 0.1*float(inlon(k)) if (iloncon .eq. 1) then c For deg W negative convention, remove track points near c right and left boundaries of lsdiag.f computational domain c for typical A*PACK files if (rlon(k) .ge. 15.0 .and. rlon(k) .le. 165.0) then rlat(k) = rmiss rlon(k) = rmiss endif endif enddo c c Check for missing lat,lon do k=1,mft ilatt = nint(rlat(k)) if (ilatt .eq. 0) then rlat(k) = rmiss rlon(k) = rmiss endif enddo c c Get t=0 lat,lon values from the stormcard input. rlat(0) = alat rlon(0) = alon c c Convert long to deg W negative for NHC/CPHC domain forecasts. c Skip small area east of Greenwich Meridian to avoid discontinuity. if (iloncon .eq. 1) then do k=0,mft if (rlon(k) .gt. 15.0 .and. rlon(k) .lt. rmiss) then rlon(k) = rlon(k)-360.0 endif enddo endif 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 Calculate t=-12 h lat,lon from the t=0 position and storm motion vector. im12 = -12/idelt rlat(im12) = rlat(0) - dlat rlon(im12) = rlon(0) - dlon 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) exit endif enddo 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) exit endif enddo 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 Initialization for new ncodaclimo2 routine ioper = 1 c Calculate the sea surface temperature, ocean heat content and distance c to land along the storm track ifirstw=1 ifirstd=1 ifirsto=1 do i=0,mft if (rlat(i) .lt. 500.0 .and. rlon(i) .lt. 500.0) then rlatt = rlat(i) rlont = rlon(i) c c Longitude in 0 to 360 convention needed for some routines rlonto = rlont if (rlonto .lt. 0.0) rlonto = rlonto + 360.0 c c ++ Distance to land call gdland_table(rlonto,rlatt,distt) if (distt .gt. 5000.0) distt=5000.0 if (distt .lt. -999.0) distt=-999.0 c c ++ Climo SST call rsstclim(isyr,ismon,isday,rlatt,rlonto,rmiss, + sstt,ierr) if (sstt .gt. 50.0 .or. ierr .ne. 0) sstt=999.9 write(lulg,880) sstt,ierr 880 format('call to rsstclim, sstt=',f8.3,1x,' ierr=',i3) c c ++ Climo D20, D26, OHC call ncodaclim(rlatt,rlonto,isyr,ismon,isday, + ohctem,d20tem,d26tem,ssttem,ierr) c write(6,992) i,ssttem,ohctem,d20tem,d26tem c 992 format('i,ssttem,ohc,d20,d26: ',i2,4(f8.1)) if (ierr .eq. 0) then if (ohctem .gt. -100.0) then if (ohctem .le. 0.0) ohctem = 0.0 iohcc(i) = nint(ohctem) else iohcc(i) = imiss endif c if (d20tem .ge. 0.0) then id20c(i) = nint(d20tem) else if (ssttem .ge. 0.0 .and. ssttem .lt. 20.0) then id20c(i) = 0 else id20c(i) = imiss endif endif c if (d26tem .ge. 0.0) then id26c(i) = nint(d26tem) else if (ssttem .ge. 0.0 .and. ssttem .lt. 26.0) then id26c(i) = 0 else id26c(i) = imiss endif endif else id20c(i) = imiss id26c(i) = imiss iohcc(i) = imiss endif c ++ Observed weekly SST c Note: Weekly SST no longer available after Nov 2022 c so only call rysst in 2022 or earlier ierr=1 if (iyr4 .le. 2022) then call rysst(rlatt,rlont,rsstt, + iryr,irmon,irday,irtime,ifirstw,ierr) endif 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) irdelt = irdelt/24 endif c c ++ Observed daily SST ierr=1 distavg = 50.0 call rydsst(rlatt,rlont,dsstt,dsstta,distavg, + iryrd,irmond,irdayd,irtimed,ifirstd,ierr) if (ierr .ne. 0) then rssttd =999.9 rssttda=999.9 endif c c After the first call to rydsst, calculate the time lag between c the stormcard date and the DSST data file date. if (i .eq. 0 .and. ierr .eq. 0) then call tdiff(isyr ,ismon ,isday ,istime, + iryrd,irmond,irdayd,irtimed,irdeltd) irdeltd = irdeltd/24 endif c c ++ Observed OHC c (old formatted file, commented out in ver2.3.0) c c ierr=1 c call ryohc(rlatt,rlont,ohct,ioyear,iomon,ioday,iotime, c + ifirsto,ierr) c print *, "OHC-ryohc", rlatt,rlont,ohct,ierr c if (ierr .ne. 0) ohct=9999. c 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. c if (i .eq. 0) then c if (ierr .eq. 0) then c ioyr = ioyear - 100*(ioyear/100) c call tdiff(isyr,ismon,isday,istime, c + ioyr,iomon,ioday,iotime,iodelt) c iodelt = iodelt/24 c c else c iodelt=9999 c endif c endif c c ++ OHC predicted from climo and observed SST (AL, EP, CP, WP only) c (Commented out in ver2.3.0) c c Calculate the proxy OHC data c call ohcparam2(ibasin,iyr4,ismon,isday,rlatt,rlonto, c + rsstt,rmiss,pohc,ierr) c c if (ierr .ne. 0) then c ipohc(i) = imiss c else c ipohc(i) = nint(pohc) c endif c sst(i) = sstt rsst(i) = rsstt dsst(i) = dsstt dssta(i) = dsstta ohc(i) = ohct dist(i) = distt else sst(i) = 999.9 rsst(i) = 999.9 dsst(i) = 999.9 dssta(i) = 999.9 ohc(i) = 9999. dist(i) = 9999. ipohc(i) = imiss iohcc(i) = imiss id20c(i) = imiss id26c(i) = imiss endif c c ++ Convert to scaled integers for lsdiag.dat file isst(i) = nint(10.0*sst(i)) irsst(i) = nint(10.0*rsst(i)) idsst(i) = nint(10.0*dsst(i)) idssta(i) = nint(10.0*dssta(i)) iohc(i) = nint(ohc(i)) idist(i) = nint(dist(i)) enddo c do k=-2,mft if (rlat(k) .lt. rmiss .and. abs(rlon(k)) .lt. rmiss) then ilat(k) = nint(10.0*rlat(k)) ilon(k) = nint(10.0*rlon(k)) if (iloncon .eq. 1) ilon(k) = -ilon(k) else ilat(k) = imiss ilon(k) = imiss endif enddo c c Write the lsdiag.inp file c sname4 = sname(1:4) sname4 = cbasin6(1:4) c vlab = 'HEAD' write(luls,200) sname4,isyr,ismon,isday,istime,ivmx0, + rlat(0),rlon(0),idum99,cbasin8,vlab 200 format(1x,a4,1x,3(i2.2),1x,i2.2,2x,i3, + 1x,f6.1,1x,f6.1,1x,i4,1x,a8,108x,a4) c vlab = 'TIME' write(luls,201) (intime(k),k=-2,mft),vlab 201 format(31(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(31(1x,i4),1x,a4) c vlab='LAT ' write(luls,215) (ilat(k),k=-2,mft),vlab 215 format(31(1x,i4),1x,a4) vlab='LON ' if (iloncon .eq. 1) then lonlab = 'W' else lonlab = 'E' endif write(luls,216) (ilon(k),k=-2,mft),vlab,lonlab 216 format(31(1x,i4),1x,a4,4x,a1) 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,29(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,29(1x,i4),1x,a4,1x,i4) c vlab='DSST' if (irdeltd .gt. 9999) irdelt=9999 write(luls,221) (idsst(k),k=0,mft),vlab,irdeltd c vlab='DSTA' write(luls,221) (idssta(k),k=0,mft),vlab,irdeltd c c Check RSST and RHCN values for consistenty c do k=0,mft c write(6,*) k,irsst(k) c if (irsst(k) .gt. 260 .and. irsst(k) .lt. 350) then c if (iohc(k) .eq. 0) iohc(k) = 9999 c endif c enddo c vlab='CD20' write(luls,220) (id20c(k),k=0,mft),vlab c vlab='CD26' write(luls,220) (id26c(k),k=0,mft),vlab c vlab='COHC' write(luls,220) (iohcc(k),k=0,mft),vlab c c Also write iohc to RHCN for now so 2017 SHIPS model can be run c vlab='RHCN' c write(luls,221) (iohc(k),k=0,mft),vlab,iodelt c c vlab='PHCN' c write(luls,220) (ipohc(k),k=0,mft),vlab c c vlab='LAST' write(luls,220) (iadum99(k),k=0,mft),vlab c c Make the input file for getting the IR profiles on WCOSS open(unit=luir,file=fnir,form='formatted',status='replace', + err=900) c c Write the storm position, where lon is in the 0 to 360 convention. slatn = rlat(0) if (iloncon .eq. 1) then slonn = 360.0+rlon(0) else slonn = rlon(0) endif c write(luir,300) slatn,slonn,ihead,ispeed 300 format(f5.1,f6.1,1x,i4,1x,i4) c c Determine satellite data file prefix based on longitude c East of 105W uses GOES-16, West of 105W uses GOES-15 if (slonn .ge. loncut) then satlab='GE/CH13_IR_FDR' irmin = 0 elseif (slonn .lt. loncut) then c satlab='GW/GWIR' satlab='GW/CH13_IR_FDR' irmin = 0 endif irmin1 = irmin + 30 c ihra = -2 call tadd(iyr4 ,ismon, isday ,istime ,ihra, + iyr4m1,ismonm1,isdaym1,istimem1) c ihra = -3 call tadd(iyr4 ,ismon, isday ,istime ,ihra, + iyr4m3,ismonm3,isdaym3,istimem3) c c if (satlab .eq. "GW/GWIR") then cc Write the name of the IR file at t= 0 hr c write(luir,310,err=900) satlab,iyr4,ismon,isday,istime,irmin cc cc Write the name of the IR file at t= -1.5 hr c write(luir,310,err=900) satlab,iyr4m1,ismonm1,isdaym1, c + istimem1,irmin1 cc cc Write the name of the IR file at t=-3 hr c write(luir,310,err=900) satlab,iyr4m3,ismonm3,isdaym3, c + istimem3,irmin c 310 format(a7,'_',i4.4,i2.2,i2.2,'_',i2.2,i2.2) c endif c if (satlab .eq. "GE/CH13_IR_FDR") then c Write the name of the IR file at t= 0 hr write(luir,311,err=900) satlab,iyr4,ismon,isday,istime,irmin c c Write the name of the IR file at t= -1.5 hr write(luir,311,err=900) satlab,iyr4m1,ismonm1,isdaym1, + istimem1,irmin1 c c Write the name of the IR file at t=-3 hr write(luir,311,err=900) satlab,iyr4m3,ismonm3,isdaym3, + istimem3,irmin 311 format(a14,'_',i4.4,i2.2,i2.2,'_',i2.2,i2.2) c endif c close(luir) c write(lulg,650) 650 format(/,'ilsin completed normally') c stop c c Error processing 900 write(lulg,950) 950 format(/,'I/O error in ilsin') stop c end