program mkirprof
c
c     Version 2.0.3
c
c     This program makes the azimuthally averaged IR profiles
c     for the SHIPS and related models. The following input
c     files are needed:
c        
c     Input:
c        istormcard.dat  - Contains basic storm information
c        irdump1.dat     - The IR data closest to t = 0   hr
c        irdump2.dat     - The IR data closest to t =-1.5 hr
c        irdump3.dat     - The IR data closest to t =-3   hr
c
c     Output:
c        IRRP1.dat - The IR profile closest to t = 0 hr
c        IRRP3.dat - The IR profile closest to t= -1.5 hr
c        IRRP3.dat - The IR profile closest to t= -3 hr
c        IRRP1.inf - Basic satellite info for profile 1
c        IRRP2.inf - Basic satellite info for profile 2
c        IRRP3.inf - Basic satellite info for profile 3
c
c     Modified May 21, 2013 for WCOSS (MD)
c     Modified May 02, 2016 to add third radial profile (ver 2.0.0)
c                           IRPC processing also added. 
c     Modfifed May 11, 2016 to correct bug to create all 
c                           IRPC lines instead of just one (ver 2.0.1)
c     Modified May 15, 2018 corrected IRPC rotation (declared ispd/idir as real)
c     Modified Jun 27, 2019 modified btmin/btmax for QC (SS)
c

      parameter (mx=1024,my=1024)
      dimension bt(mx,my),rlat(mx,my),rlon(mx,my)
c
      parameter (mxr=1024,mxpc=100)
      dimension radkm(mxr),btmean(mxr),btstd(mxr),ppcount(mxr,mxpc)
c
      character*3,dimension(12) :: mname
      data mname /'JAN','FEB','MAR','APR','MAY','JUN',
     +            'JUL','AUG','SEP','OCT','NOV','DEC'/
c
      real ispd,idir
      
      character *6 atcfid6
      character *11 sname
      character *11 irinp1,irinp2,irinp3
      character *9  irdat1,irdat2,irdat3
      character *9  irinf1,irinf2,irinf3
c
c     Variables for IRPC processing
      character *24 fnirpc
      logical firstsat
c
c     I/O variables
      data luinp,lustm,ludat,luinf /10,11,12,13/
c
c     Specify input and output file names
      irinp1 = 'irdump1.dat'
      irinp2 = 'irdump2.dat'
      irinp3 = 'irdump3.dat'
c
      irdat1 = 'IRRP1.dat'
      irdat2 = 'IRRP2.dat'
      irdat3 = 'IRRP3.dat'
c
      irinf1 = 'IRRP1.inf'
      irinf2 = 'IRRP2.inf'
      irinf3 = 'IRRP3.inf'
c
c     Open and read the stormcard info
      open(file='istormcard.dat',unit=lustm,form='formatted',
     +     status='old',err=900)
c
      read(lustm,201,err=900,end=900) iyrs2,imons,idays
  201 format(3(i2))
      read(lustm,202,err=900,end=900) ihrs
  202 format(i2)
      read(lustm,203,err=900,end=900) islat00,islatm12
  203 format(i4,1x,i4)
      read(lustm,204,err=900,end=900) islon00,islonm12
  204 format(i4,1x,i4)
      read(lustm,*) idum1
      read(lustm,*) idum2
      read(lustm,*) ivmax00,ivmaxm12
      read(lustm,*) idir
      read(lustm,*) ispd
      read(lustm,206) sname
  206 format(a11)
      read(lustm,208) atcfid6
  208 format(a6)
c
      close(lustm)
c
      if (iyrs2 .le. 50) then 
         iyears = iyrs2 + 2000
      else
         iyears = iyrs2 + 1900
      endif
c
      slat00  = 0.1*float(islat00)
      slatm12 = 0.1*float(islatm12)
      slon00  = 0.1*float(islon00)
      slonm12 = 0.1*float(islonm12)
c
      if (atcfid6(1:2) .eq. 'AL' .or.
     +    atcfid6(1:2) .eq. 'EP' .or.
     +    atcfid6(1:2) .eq. 'CP'      ) then
         slon00 = 360.0  - slon00
         slonm12 = 360.0 - slonm12
      endif
c
      t00 =   0.0
      tm12= -12.0
c
      adlats = abs(slat00-slatm12)
      adlons = abs(slon00-slonm12)
c
c     Variables for IRPC routine
      rstride=1.5
      firstsat=.true.
      write(fnirpc,333) iyrs2,imons,idays,ihrs,atcfid6
  333 format(4(i2.2),a6,'_IRPC0.dat')
      if (atcfid6(1:2) .eq. 'AL') then
         offsethr = 0.0
      else
         offsethr = -.25
      endif
c
c     Read the first IR data file
      call irundump(irinp1,luinp,mx,my,ifound,nelex,nliny,
     +              iyear,imon,iday,ihr,imin,isec,bt,rlon,rlat)
c
      if (imon .lt. 1 .or. imon .gt. 12) go to 800
c
      if (ifound .eq. 1) then
c
c        Interpolate storm center to the time of image 1
         if (adlats .lt. 20.0 .and. adlons .lt. 20.0) then
c
            call tdiff(iyear ,imon ,iday ,ihr ,
     +                 iyears,imons,idays,ihrs, idelt)
            timage = float(idelt) + float(imin)/60.0
c
            if (timage .gt. -13.0 .and. timage .lt. 13.0) then
               w00  = (timage-tm12)/(t00-tm12)
               wm12 = 1.0 - w00
c
               slati = w00*slat00 + wm12*slatm12
               sloni = w00*slon00 + wm12*slonm12
            else
               slati = slat00
               sloni = slon00
            endif
         else
            slati = slat00
            sloni = slon00
          endif
c 
c         Make radial profile 1
          call irradprof(sloni,slati,bt,rlon,rlat,mx,my,nelex,nliny,
     +                   radkm,btmean,btstd,ppcount,mxr,nxr,mxpc,npc)
c
c         Write .dat for profile 1
          open(file=irdat1,unit=ludat,form='formatted',
     +         status='replace',err=800)
c
          do k=1,nxr
             write(ludat,250) radkm(k),btmean(k),btstd(k),
     +                        (ppcount(k,m),m=1,npc)
  250        format(f6.1,20(f7.2))
          enddo
c
c         Write .inf for profile 1
          open(file=irinf1,unit=luinf,form='formatted',
     +         status='replace',err=800)
c
          write(luinf,300) slati,sloni,nelex,nliny
  300     format('Center LAT, Center LON,  im_x, im_y',/,
     +            f7.2,f8.2,i4,i4,/,
     +            'GEO IR')
c
          write(luinf,310) iday,mname(imon),iyear,ihr,imin
  310     format('IDAY, MONTH, YEAR, TIME',/,
     +                1x,i2.2,2x,a3,2x,i4,1x,i2.2,i2.2)
c
          write(luinf,320) rlat(    1,1),rlat(nelex,    1),
     +                     rlat(1,nliny),rlat(nelex,nliny)
  320     format('Corner Lat',/,4(f8.2))
c
          write(luinf,330) rlon(    1,1),rlon(nelex,    1),
     +                     rlon(1,nliny),rlon(nelex,nliny)
  330     format('Corner Lon',/,4(f8.2))
c
c         Create IRPC data for first image
          nl = 1
          numsat = 1
          timagep = 0.25 + offsethr
          call mkirpc(irinp1,fnirpc,nl,numsat,timagep,rstride,
     +                       ispd,idir,slati,sloni,firstsat)
          firstsat=.false.
      endif
c
  800 continue
      close(ludat)
      close(luinf)
c
c     Read the second IR data file
      call irundump(irinp2,luinp,mx,my,ifound,nelex,nliny,
     +              iyear,imon,iday,ihr,imin,isec,bt,rlon,rlat)
c
      if (imon .lt. 1 .or. imon .gt. 12) go to 900
c
      if (ifound .eq. 1) then
c
c        Interpolate storm center to the time of image 2
         if (adlats .lt. 20.0 .and. adlons .lt. 20.0) then
c
            call tdiff(iyear ,imon ,iday ,ihr ,
     +                 iyears,imons,idays,ihrs, idelt)
            timage = float(idelt) + float(imin)/60.0
c
            if (timage .gt. -13.0 .and. timage .lt. 13.0) then
               w00  = (timage-tm12)/(t00-tm12)
               wm12 = 1.0 - w00
c
               slati = w00*slat00 + wm12*slatm12
               sloni = w00*slon00 + wm12*slonm12
            else
               slati = slat00
               sloni = slon00
            endif
         else
            slati = slat00
            sloni = slon00
          endif
c 
c         Make radial profile 2
          call irradprof(sloni,slati,bt,rlon,rlat,mx,my,nelex,nliny,
     +                   radkm,btmean,btstd,ppcount,mxr,nxr,mxpc,npc)
c
c         Write .dat for profile 2
          open(file=irdat2,unit=ludat,form='formatted',
     +         status='replace',err=900)
c
          do k=1,nxr
             write(ludat,250) radkm(k),btmean(k),btstd(k),
     +                        (ppcount(k,m),m=1,npc)
          enddo
          close(ludat)
c
c         Write .inf for profile 2
          open(file=irinf2,unit=luinf,form='formatted',
     +         status='replace',err=900)
c
          write(luinf,300) slati,sloni,nelex,nliny
          write(luinf,310) iday,mname(imon),iyear,ihr,imin
          write(luinf,320) rlat(    1,1),rlat(nelex,    1),
     +                     rlat(1,nliny),rlat(nelex,nliny)
          write(luinf,330) rlon(    1,1),rlon(nelex,    1),
     +                     rlon(1,nliny),rlon(nelex,nliny)
c
c         Create IRPC data for second image
          nl     = 2
          numsat = 2
          timagep = -1.25 + offsethr
          call mkirpc(irinp2,fnirpc,nl,numsat,timagep,rstride,
     +                       ispd,idir,slati,sloni,firstsat)
          firstsat=.false.
      endif
c
  900 continue
      close(ludat)
      close(luinf)
c
c     Read the third IR data file
      call irundump(irinp3,luinp,mx,my,ifound,nelex,nliny,
     +              iyear,imon,iday,ihr,imin,isec,bt,rlon,rlat)
c
      if (imon .lt. 1 .or. imon .gt. 12) go to 950
c
      if (ifound .eq. 1) then
c
c        Interpolate storm center to the time of image 2
         if (adlats .lt. 20.0 .and. adlons .lt. 20.0) then
c
            call tdiff(iyear ,imon ,iday ,ihr ,
     +                 iyears,imons,idays,ihrs, idelt)
            timage = float(idelt) + float(imin)/60.0
c
            if (timage .gt. -13.0 .and. timage .lt. 13.0) then
               w00  = (timage-tm12)/(t00-tm12)
               wm12 = 1.0 - w00
c
               slati = w00*slat00 + wm12*slatm12
               sloni = w00*slon00 + wm12*slonm12
            else
               slati = slat00
               sloni = slon00
            endif
         else
            slati = slat00
            sloni = slon00
          endif
c 
c         Make radial profile 3
          call irradprof(sloni,slati,bt,rlon,rlat,mx,my,nelex,nliny,
     +                   radkm,btmean,btstd,ppcount,mxr,nxr,mxpc,npc)
c
c         Write .dat for profile 3
          open(file=irdat3,unit=ludat,form='formatted',
     +         status='replace',err=900)
c
          do k=1,nxr
             write(ludat,250) radkm(k),btmean(k),btstd(k),
     +                        (ppcount(k,m),m=1,npc)
          enddo
          close(ludat)
c
c         Write .inf for profile 3
          open(file=irinf3,unit=luinf,form='formatted',
     +         status='replace',err=900)
c
          write(luinf,300) slati,sloni,nelex,nliny
          write(luinf,310) iday,mname(imon),iyear,ihr,imin
          write(luinf,320) rlat(    1,1),rlat(nelex,    1),
     +                     rlat(1,nliny),rlat(nelex,nliny)
          write(luinf,330) rlon(    1,1),rlon(nelex,    1),
     +                     rlon(1,nliny),rlon(nelex,nliny)
c
c         Create IRPC data for third image
          nl     = 3
          numsat = 3
          timagep = -2.75 + offsethr
          call mkirpc(irinp3,fnirpc,nl,numsat,timagep,rstride,
     +                       ispd,idir,slati,sloni,firstsat)
          firstsat=.false.
      endif
c
  950 continue
      close(ludat)
      close(luinf)
c
      stop
      end
      subroutine irundump(fninp,luinp,mx,my,ifound,nelex,nliny,
     +                    iyear,imon,iday,ihr,imin,isec,bt,rlon,rlat)
c     This routine reads the lat,lon and IR brightness information
c     from an ascii file created by routine irdump.
c
      dimension bt(mx,my),rlon(mx,my),rlat(mx,my)
      character *11 fninp 
c
      open(file=fninp,unit=luinp,form='formatted',status='old',
     +     err=900)
c
c     Read the header
      read(luinp,100,err=900,end=900) iyear,imon,iday,ihr,imin,isec,
     +                                nelex,nliny 
  100 format(i4,i2,i2,1x,3(i2),1x,i6,1x,i6)
c
c     Read the BT
      read(luinp,110,end=900,err=900) ((  bt(i,j),i=1,nelex),j=1,nliny)
      read(luinp,110,end=900,err=900) ((rlon(i,j),i=1,nelex),j=1,nliny)
      read(luinp,110,end=900,err=900) ((rlat(i,j),i=1,nelex),j=1,nliny)
  110 format(20(f9.3,1x))
c
c     Normal return
      close(luinp)
      ifound = 1
      return
c
  900 continue
c     File not found or error reading file
      close(luinp)
      ifound = 0
c
      end
      subroutine irradprof(sloni,slati,bt,rlon,rlat,mx,my,nelex,nliny,
     +                     radkm,btmean,btstd,ppcount,mxr,nxr,mxpc,npc)
c
c     This routine calculates radial profiles of azimuthally average brightness
c     IR brightness temperature, std dev of BT, and percent pixel counts
c     colder than various thresholds.
c
      dimension bt(mx,my),rlon(mx,my),rlat(mx,my)
      dimension radkm(mxr),btmean(mxr),btstd(mxr),ppcount(mxr,mxpc)
c
c     Local variables
      parameter (mxrl=1024,mxpcl=100)
      dimension tthresh(mxpcl),rcount(mxrl)
      parameter (mxl=1024,myl=1024)
      dimension radxy(mxl,myl)
c
c     Set default value for missing IR calculations
      rmiss = -1.0
c
c     Specify range acceptable brightness temperatures (K)
c      btmin = 140.0
c      btmax = 340.0
      btmin = 183.0 
      btmax = 313.0
c
c     Specify max radial grid size (km)
      rgmax = 600.0
c
c     Specify BT thresholds for pixel count variables
      npc = 8
      tthresh(1) = 273.0
      tthresh(2) = 263.0
      tthresh(3) = 253.0
      tthresh(4) = 243.0
      tthresh(5) = 233.0
      tthresh(6) = 223.0
      tthresh(7) = 213.0
      tthresh(8) = 203.0
c
c     Specify radial grid (no. of points calcualted later depending on 
c                          image size and storm position)
      drkm = 4.0
      r0   = 2.0
c     
c     Check to make sure drkm is not too small. If so, make it bigger
      dycheck = 111.1*0.5*
     +          ( abs(rlat(    1,    2) - rlat(    1,      1)) +
     +            abs(rlat(nelex,nliny) - rlat(nelex,nliny-1)) )
c
      if (dycheck .ge. 6.0) then
         drkm = 10.0
         r0   = 5.0 
      endif
c
c     Calculate the radius from the storm center to each pixel
      do i=1,nelex
      do j=1,nliny
         call distk(sloni,slati,rlon(i,j),rlat(i,j),dx,dy,radxy(i,j))
cx       possibly two different longitude definitions, take care of it here
         if (radxy(i,j) .gt. 1500.0) then
           call distk
     &     ( (360.0-sloni),slati,rlon(i,j),rlat(i,j),dx,dy,radxy(i,j))
         endif 
      enddo
      enddo
c
c     Calculate distance from storm to closest edge of the domain
      redge = 1.0e+10
      do i=1,nelex
         if (redge .gt. radxy(i,    1)) redge=radxy(i,    1)
         if (redge .gt. radxy(i,nliny)) redge=radxy(i,nliny)
      enddo
c
      do j=1,nliny
         if (redge .gt. radxy(    1,j)) redge=radxy(1,    j)
         if (redge .gt. radxy(nelex,j)) redge=radxy(nelex,j)
      enddo
c
      redge = 100.0*float(ifix(redge/100.0))
      if (redge .gt. rgmax) redge = rgmax
c
c     Calculate center point of radial grid intervals
      nxr = ifix(0.0001 + redge/drkm)
c
      do k=1,nxr
         radkm(k) = r0 + drkm*float(k-1)
      enddo
c
c     Initialize counters and calculated variables
      do k=1,nxr
         btmean(k) = 0.0
         btstd(k)  = 0.0
         rcount(k) = 0.0
         do m=1,npc
            ppcount(k,m) = 0.0
         enddo
      enddo
c
c     Calculate variables
      do i=1,nelex
      do j=1,nliny
c        Check for bt out of range
         if (bt(i,j) .lt. btmin .or. bt(i,j) .gt. btmax) go to 1000
c
         indexr = 1 + ifix(radxy(i,j)/drkm)
         if (indexr .gt. nxr) go to 1000
c
         rcount(indexr) = rcount(indexr) + 1.0         
         btmean(indexr) = btmean(indexr) + bt(i,j)
         btstd(indexr)  = btstd(indexr)  + bt(i,j)**2
c
         do m=1,npc
            if (bt(i,j) .le. tthresh(m)) ppcount(indexr,m) = 
     +                                   ppcount(indexr,m) + 1.0
         enddo
 1000    continue
      enddo
      enddo
c
      do k=1,nxr
         if (rcount(k) .ge. 2) then
            btmean(k) = btmean(k)/rcount(k)
c           btstd(k)  = sqrt( (btstd(k) - rcount(k)*(btmean(k)**2)) /
c    +                         (rcount(k)-1.0)) 
c
            btemp  =  (btstd(k) - rcount(k)*(btmean(k)**2)) /
     +                           (rcount(k)-1.0)
c
c           Check to make sure sqrt argument in stdev calculation
c           does not go negative due to round off error. If it does, set
c           btstd to very small value because problem only occurs when std is
c           close to zero. 
c
            if (btemp .gt. 0.0) then
                btstd(k) = sqrt(btemp)
            else
                btstd(k) = 0.1
            endif
c
            do m=1,npc
               ppcount(k,m) = 100.0*(ppcount(k,m)/rcount(k))
            enddo
         else
            btmean(k) = rmiss
            btstd(k)  = rmiss
         endif
      enddo
c
      return
      end