c
      subroutine rapidloge_2018( perri, ilat, iepss, ienss, irhlo,
     +                  irhmd, id200, ishdc, ishgc, vsst,
     +                  vmx, sir5ri, sir6ri, sir10ri, sir12ri,
     +                  itpw1, itpw4, itpw6, icflx, probrilog)
c
c     subroutine rapidloge.f
c
c     v2.2.2
c
c     Updated: 21 Feb 2018 for 2018 season
c     Updated: 11 Apr 2020 for 7day input (SS)
c
c     Authors: 
c          Christopher M. Rozoff (chris.rozoff@ssec.wisc.edu)
c          James P. Kossin (kossin@ssec.wisc.edu)
c
c     Purpose: This subroutine is used to compute the probability of RI 
c          for East Pacific tropical cyclones via a logistic regression-based 
c          model (Rozoff and Kossin; Wea. Forecasting). This model 
c          provides estimates of the probability of RI for the 
c          20 kt / 12 h, 25 kt / 24 h, 30 kt / 24 h, 35 kt / 24 h, 
c          40 kt / 24 h, 45 kt / 36 h, 55 kt / 48 h RI,
c          and 65 kt / 72 h RI thresholds.
c
c     Notes: This subroutine assumes that missing input values
c         are equal to  9999 (e.g., perri=9999 if missing). Also, 
c         missing output values will be coded as 999.
c
c     subroutine inputs:
c
c     perri   : 12 h intensity change observed for the preceding 12 h (kt)
c     latriX  : latitude averaged from t=0 to t=X h
c     epssriX : The average theta_e difference between a parcel lifted from
c             the surface and its saturated environment (r = 200-800 km). Only
c             positive differences are included in the average.
c     enssriX : Vertical average of the negative differences between the 
c             equivalent potential temperature of a parcel lifted from the surface 
c             and the environmental saturation equivalent potential temperature 
c             (r = 200–800 km) averaged from T=0 to T=X h
c     rhloriX : 850-700 hPa relative humidity (%) vs time (200-800 km)
c     d200riX : 200-hPa divergence vs time (r=0-1000 km)
c     shdcriX : 850-200 mb vertical shear (kt) computed from 0-500 km radius a
c             after the vortex has been removed and averaged from T=0 to T=X h
c     shgcriX : 850-200-hPa shear magnitude with vortex removed and averaged from
c             0-500 km relative to 850 hPa vortex center
c     icflx   : Array containing Boundary-layer latent heat flux (watts/m2) difference 
c             predictor (q10layr-q10)*Vmx0 where q10 is obtained by bring air down 
c             from 1000 mb and then allowing air to cool as it reaches 
c             95% RH in the core. q10layr is obtained using the same 
c              methodology except using RH= 1/3* (RH_sfc + Rh_850-700 mb 
c              + RH_700-500 mb) (time-averaged).
c     potriX  : Potential intensity (kt) averaged from T=0 to T=X h.
c             The MPI that is used when computing the potential intensity 
c             is determined using Joe Cione's inner-core cooled SST (number 3). 
c             This MPI also accounts for storm speed based upon the algorithm 
c             developed by Schwerdt et al. (1979) and is capped at 165 kt.
c     sir5ri   : Stan. Dev. of GOES BT (deg C*10), r=100-300 km
c     sir6ri   : Percent area r=50-200 km of GOES ch 4 BT < -10 C 
c     sir10ri  : Percent area r=50-200 km of GOES ch 4 BT < -50 C
c     sir12ri  : max BT from 0 to 30 km radius (deg C*10)
c     tpw1   : TPW averaged over r = 0-200 km averaged from t=0 to t=X h.
c     tpw4X  : TPW standard deviation 200-400 km averaged from
c	      T=0 to T=X h.
c     tpw6X  : TPW standard deviation 400-600 km averaged from
c	      T=0 to T=X h.
c
c     subroutine output:
c 
c     probrilog : the forecasted probabilities of RI for each of the 8
c          RI thresholds
c
c     other relevant variables: 
c     Mft     : number of predictor times 
c     Nindx   : number of RI index thresholds computed
c     Nthrss  : number of RI predictors
c     Nlt     : number of lead times
c     xfeat   : feature vector containing predictors
c
c ------------------------------------------------------------------------------
c
      parameter (Mft = 28)
c
      parameter (Nindx = 8)
c
      parameter (Nthrss = 18)
c
      parameter (Nlt = 5)
c
c input
c
      dimension ishdc(0:Mft), ishgc(0:Mft), vsst(0:Mft)
      dimension ienss(0:Mft), iepss(0:Mft), irhlo(0:Mft)
	  dimension id200(0:Mft), ilat(-2:Mft), icflx(0:Mft)
	  dimension itpw1(0:Mft), itpw4(0:Mft), itpw6(0:Mft)
c
c  logistic regression coefficients for each RI threshold
c
      dimension coef(Nindx, Nthrss)
c
c internal data
c
      integer ideltat, itimeinc
      parameter (ideltat = 12, itimeinc = 6)
c
c 20 kt / 12-h RI threshold
      data (coef(1,i), i = 1, Nthrss)
     * / 1.2209804e+00,  8.7730916e-02, -8.9767385e-02, 
     *   0.0000000e+00, -2.6443750e-01, -4.6571404e-02, 
     *   0.0000000e+00, -1.6549010e-01,  0.0000000e+00, 
     *   0.0000000e+00,  1.9663133e-02, -8.9449259e-02,
     *   0.0000000e+00,  2.1433632e-02, -2.1009249e-02,
     *   0.0000000e+00,  0.0000000e+00,  0.0000000e+00/
c 25 kt / 24-h RI threshold
      data (coef(2,i), i = 1, Nthrss)
     * /-6.2039298e+00,  5.7316404e-02, -1.4268996e-01, 
     *   0.0000000e+00, -1.6300751e-01,  0.0000000e+00, 
     *   0.0000000e+00, -1.6672966e-01,  0.0000000e+00, 
     *   0.0000000e+00,  3.4943671e-02, -8.5928661e-02,
     *   0.0000000e+00,  1.3954956e-02, -2.4680896e-02,
     *   7.2557844e-03,  0.0000000e+00,  0.0000000e+00/
c 30 kt / 24-h RI threshold
      data (coef(3,i), i = 1, Nthrss)
     * /-4.5976194e+00,  7.3382353e-02, -1.2573638e-01, 
     *   1.0852531e-01,  0.0000000e+00,  0.0000000e+00, 
     *   5.7261973e-03, -2.1673989e-01,  0.0000000e+00, 
     *   2.4025490e-03,  4.2444222e-02, -7.9752456e-02,
     *   0.0000000e+00,  1.4537975e-02, -2.5783472e-02,
     *   0.0000000e+00,  0.0000000e+00,  0.0000000e+00/
c 35 kt / 24-h RI threshold
      data (coef(4,i), i = 1, Nthrss)
     * /-4.2832624e+00,  8.2532131e-02, -1.0842655e-01, 
     *   0.0000000e+00, -3.4210924e-01,  0.0000000e+00, 
     *   5.8364603e-03, -2.2501453e-01,  0.0000000e+00, 
     *   3.1856627e-03,  4.4006943e-02, -9.9575180e-02,
     *   0.0000000e+00,  1.3698791e-02, -2.8795436e-02,
     *   0.0000000e+00,  0.0000000e+00,  0.0000000e+00/
c 40 kt / 24-h RI threshold
      data (coef(5,i), i = 1, Nthrss)
     * /-4.7425951e+00,  8.2843230e-02, -1.5498203e-01, 
     *   2.4106060e-01,  0.0000000e+00,  0.0000000e+00, 
     *   4.4153476e-03,  0.0000000e+00, -2.0254987e-01, 
     *   4.6496317e-03,  4.4450559e-02, -1.0451886e-01,
     *   0.0000000e+00,  1.5184583e-02, -3.4670601e-02,
     *   0.0000000e+00,  0.0000000e+00,  0.0000000e+00/
c 45 kt / 36-h RI threshold
      data (coef(6,i), i = 1, Nthrss)
     * /-7.4998748e+00,  8.1362542e-02,  0.0000000e+00, 
     *   0.0000000e+00, -4.4633696e-01,  0.0000000e+00, 
     *   1.0427654e-02, -2.2537346e-01,  0.0000000e+00, 
     *   5.6076425e-03,  7.5096203e-02, -9.4472311e-02,
     *   0.0000000e+00,  1.0757435e-02, -2.0531136e-02,
     *   0.0000000e+00, -3.1914773e-02,  0.0000000e+00/
c 55 kt / 48-h RI threshold
      data (coef(7,i), i = 1, Nthrss)
     * /-1.1680254e+01,  8.1628555e-02,  0.0000000e+00, 
     *   0.0000000e+00, -5.4295407e-01,  0.0000000e+00, 
     *   1.2961255e-02,  0.0000000e+00, -1.6499823e-01, 
     *   7.6544523e-03,  9.0549360e-02,  0.0000000e+00,
     *   4.5568716e-02,  0.0000000e+00,  0.0000000e+00,
     *   0.0000000e+00, -5.5813880e-02,  0.0000000e+00/
c 65 kt / 72-h RI threshold
	  data (coef(8,i), i = 1, Nthrss)
     * /-2.6811208e+00,  8.7205182e-02, -2.6445919e-01, 
     *   0.0000000e+00,  0.0000000e+00, -6.1664861e-02, 
     *   0.0000000e+00,  0.0000000e+00, -1.1502962e-01, 
     *   0.0000000e+00,  8.6852624e-02,  0.0000000e+00,
     *   3.4802905e-02,  0.0000000e+00,  0.0000000e+00,
     *   0.0000000e+00,  0.0000000e+00, -3.6794201e-02/
c
c output probabilities
c
      dimension probrilog(Nindx)
c
      dimension xfeat(Nthrss)
c
c  internal variables
c
      integer i, j
      real sumExp
      integer nlat, nshdc, nshgc, nrhlo
      integer npot, nepss, nenss, nd200, ncflx
	  integer ntpw1, ntpw4, ntpw6
      real sumlat, sumshdc, sumshgc, sumepss, sumrhlo
      real sumpot, sumenss, sumd200, sumcflx
	  real sumtpw1, sumtpw4, sumtpw6
c
      dimension rlatrilt(nlt), shdcrilt(nlt), shgcrilt(nlt)
      dimension potrilt(nlt), epssrilt(nlt), enssrilt(nlt)
	  dimension rhlorilt(nlt), d200rilt(nlt), cflxrilt(nlt)
	  dimension tpw1rilt(nlt), tpw4rilt(nlt), tpw6rilt(nlt)
c
c ------------------------------------------------------------------------------
c
c Make missing values consistent with subroutine here.
      if (int(sir5ri) .eq. 999) sir5ri = 9999.
      if (int(sir6ri) .eq. 999) sir6ri = 9999.
      if (int(sir10ri) .eq. 999) sir10ri = 9999.
      if (int(sir12ri) .eq. 999) sir12ri = 9999.  
c
c Compute mean predictors values for each of the forecast
c lead times (12, 24, 36, 48 h, and 72 h)
c
      do j = 1, nlt
          nft = (j * ideltat) / itimeinc
		  if (j .eq. 5) nft = nft + 2
		  nlat = 0
          nshdc = 0
		  nshgc = 0
          npot = 0
		  nepss = 0
          nenss = 0 
		  nd200 = 0
		  nrhlo = 0
		  ncflx = 0
		  ntpw1 = 0
		  ntpw6 = 0
		  ntpw4 = 0
c
          sumlat = 0.
          sumshdc = 0.
		  sumshgc = 0.
          sumpot = 0.
		  sumepss = 0.
          sumenss = 0.
		  sumd200 = 0.
		  sumrhlo = 0.
		  sumcflx = 0.
		  sumtpw1 = 0.
		  sumtpw6 = 0.
		  sumtpw4 = 0.
c
          do i = 0, nft
c
            if (ilat(i) .lt. 9999) then
              nlat = nlat + 1
              sumlat = sumlat + 0.1 * float(ilat(i))
            endif
c
            if (ishdc(i) .lt. 9999) then
              nshdc = nshdc + 1
              sumshdc = sumshdc + 0.1 * float(ishdc(i))
            endif
c
            if (ishgc(i) .lt. 9999) then
              nshgc = nshgc + 1
              sumshgc = sumshgc + 0.1 * float(ishgc(i))
            endif
c
            if (vsst(i) .lt. 200 .and. vmx .lt. 200) then
              npot = npot + 1
              sumpot = sumpot + (vsst(i) - vmx)
            endif
c
            if (iepss(i) .lt. 9999) then
              nepss = nepss + 1
              sumepss = sumepss + 0.1 * float(iepss(i))
            endif
c
            if (ienss(i) .lt. 9999) then
              nenss = nenss + 1
              sumenss = sumenss + 0.1 * float(ienss(i))
            endif
c
            if (irhlo(i) .lt. 9999) then
              nrhlo = nrhlo + 1
              sumrhlo = sumrhlo + float(irhlo(i))
            endif
c
            if (id200(i) .lt. 9999) then
              nd200 = nd200 + 1
              sumd200 = sumd200 + float(id200(i))
            endif
c
            if (icflx(i) .lt. 9999) then
              ncflx = ncflx + 1
              sumcflx = sumcflx + float(icflx(i))
            endif
c
            if (itpw1(i) .lt. 9999) then
              ntpw1 = ntpw1 + 1
              sumtpw1 = sumtpw1 + float(itpw1(i))
            endif
c
            if (itpw4(i) .lt. 9999) then
              ntpw4 = ntpw4 + 1
              sumtpw4 = sumtpw4 + float(itpw4(i))
            endif
c
            if (itpw6(i) .lt. 9999) then
              ntpw6 = ntpw6 + 1
              sumtpw6 = sumtpw6 + float(itpw6(i))
            endif
c
          enddo
c
          nmin = nft / 2 + 1
c
          if (nlat .ge. nmin) then
            rlatrilt(j) = sumlat / float(nlat)
          else
            rlatrilt(j) = 9999.
          endif
c
          if (nshdc .ge. nmin) then
            shdcrilt(j) = sumshdc / float(nshdc)
          else
            shdcrilt(j) = 9999.
          endif
c
          if (nshgc .ge. nmin) then
            shgcrilt(j) = sumshgc / float(nshgc)
          else
            shgcrilt(j) = 9999.
          endif
c
          if (npot .ge. nmin) then
            potrilt(j) = sumpot / float(npot)
          else
            potrilt(j) = 9999.
          endif
c
          if (nepss .ge. nmin) then
            epssrilt(j) = sumepss / float(nepss)
          else
            epssrilt(j) = 9999.
          endif
c
          if (nenss .ge. nmin) then
            enssrilt(j) = sumenss / float(nenss)
          else
            enssrilt(j) = 9999.
          endif
c
          if (nrhlo .ge. nmin) then
            rhlorilt(j) = sumrhlo / float(nrhlo)
          else
            rhlorilt(j) = 9999.
          endif
c
          if (nd200 .ge. nmin) then
            d200rilt(j) = sumd200 / float(nd200)
          else
            d200rilt(j) = 9999.
          endif
c
          if (ncflx .ge. nmin) then
            cflxrilt(j) = sumcflx / float(ncflx)
          else
            cflxrilt(j) = 9999.
          endif
c
          if (ntpw1 .ge. nmin) then
            tpw1rilt(j) = sumtpw1 / float(ntpw1)
          else
            tpw1rilt(j) = 9999.
          endif
c
          if (ntpw4 .ge. nmin) then
            tpw4rilt(j) = sumtpw4 / float(ntpw4)
          else
            tpw4rilt(j) = 9999.
          endif
c
          if (ntpw6 .ge. nmin) then
            tpw6rilt(j) = sumtpw6 / float(ntpw6)
          else
            tpw6rilt(j) = 9999.
          endif
c
      enddo
c
      do j = 1, Nindx
c
c Initialize feature vectors depending on RI threshold
c
          xfeat(1) = 1.0
          xfeat(2) = perri
c
c RI = 20 kt / 12 h
c
          if (j .eq. 1) then
c
			xfeat(3) = rlatrilt(1)
			xfeat(4) = 0.0
			xfeat(5) = enssrilt(1)
			xfeat(6) = rhlorilt(1)
			xfeat(7) = 0.0
			xfeat(8) = shdcrilt(1)
			xfeat(9) = 0.0
			xfeat(10) = 0.0
			xfeat(11) = potrilt(1)
			xfeat(12) = sir5ri
			xfeat(13) = 0.0
			xfeat(14) = sir10ri
			xfeat(15) = sir12ri
			xfeat(16) = 0.0
			xfeat(17) = 0.0
			xfeat(18) = 0.0
c
c RI = 25 kt / 24 h
c
          elseif (j .eq. 2) then
c
  			xfeat(3) = rlatrilt(2)
  			xfeat(4) = 0.0
  			xfeat(5) = enssrilt(2)
  			xfeat(6) = 0.0
  			xfeat(7) = 0.0
  			xfeat(8) = shdcrilt(2)
  			xfeat(9) = 0.0
  			xfeat(10) = 0.0
  			xfeat(11) = potrilt(2)
  			xfeat(12) = sir5ri
  			xfeat(13) = 0.0
  			xfeat(14) = sir10ri
  			xfeat(15) = sir12ri
  			xfeat(16) = tpw1rilt(2)
  			xfeat(17) = 0.0
  			xfeat(18) = 0.0
c
c RI = 30 kt / 24 h
c
          elseif (j .eq. 3) then
c
  			xfeat(3) = rlatrilt(2)
  			xfeat(4) = epssrilt(2)
  			xfeat(5) = 0.0
  			xfeat(6) = 0.0
  			xfeat(7) = d200rilt(2)
  			xfeat(8) = shdcrilt(2)
  			xfeat(9) = 0.0
  			xfeat(10) = cflxrilt(2)
  			xfeat(11) = potrilt(2)
  			xfeat(12) = sir5ri
  			xfeat(13) = 0.0
  			xfeat(14) = sir10ri
  			xfeat(15) = sir12ri
  			xfeat(16) = 0.0
  			xfeat(17) = 0.0
  			xfeat(18) = 0.0
c
c RI = 35 kt / 24 h
c
          elseif (j .eq. 4) then
c
  			xfeat(3) = rlatrilt(2)
  			xfeat(4) = 0.0
  			xfeat(5) = enssrilt(2)
  			xfeat(6) = 0.0
  			xfeat(7) = d200rilt(2)
  			xfeat(8) = shdcrilt(2)
  			xfeat(9) = 0.0
  			xfeat(10) = cflxrilt(2)
  			xfeat(11) = potrilt(2)
  			xfeat(12) = sir5ri
  			xfeat(13) = 0.0
  			xfeat(14) = sir10ri
  			xfeat(15) = sir12ri
  			xfeat(16) = 0.0
  			xfeat(17) = 0.0
  			xfeat(18) = 0.0
c
c RI = 40 kt / 24 h
c
          elseif (j .eq. 5) then
c
  			xfeat(3) = rlatrilt(2)
  			xfeat(4) = epssrilt(2)
  			xfeat(5) = 0.0
  			xfeat(6) = 0.0
  			xfeat(7) = d200rilt(2)
  			xfeat(8) = 0.0
  			xfeat(9) = shgcrilt(2)
  			xfeat(10) = cflxrilt(2)
  			xfeat(11) = potrilt(2)
  			xfeat(12) = sir5ri
  			xfeat(13) = 0.0
  			xfeat(14) = sir10ri
  			xfeat(15) = sir12ri
  			xfeat(16) = 0.0
  			xfeat(17) = 0.0
  			xfeat(18) = 0.0
c
c RI = 45 kt / 36 h
c
          elseif (j .eq. 6) then
c
  			xfeat(3) = 0.0
  			xfeat(4) = 0.0
  			xfeat(5) = enssrilt(3)
  			xfeat(6) = 0.0
  			xfeat(7) = d200rilt(3)
  			xfeat(8) = shdcrilt(3)
  			xfeat(9) = 0.0
  			xfeat(10) = cflxrilt(3)
  			xfeat(11) = potrilt(3)
  			xfeat(12) = sir5ri
  			xfeat(13) = 0.0
  			xfeat(14) = sir10ri
  			xfeat(15) = sir12ri
  			xfeat(16) = 0.0
  			xfeat(17) = tpw4rilt(3)
  			xfeat(18) = 0.0
c
c RI = 55 kt / 48 h
c
          elseif (j .eq. 7) then
c
  			xfeat(3) = 0.0
  			xfeat(4) = 0.0
  			xfeat(5) = enssrilt(4)
  			xfeat(6) = 0.0
  			xfeat(7) = d200rilt(4)
  			xfeat(8) = 0.0
  			xfeat(9) = shgcrilt(4)
  			xfeat(10) = cflxrilt(4)
  			xfeat(11) = potrilt(4)
  			xfeat(12) = 0.0
  			xfeat(13) = sir6ri
  			xfeat(14) = 0.0
  			xfeat(15) = 0.0 
  			xfeat(16) = 0.0
  			xfeat(17) = tpw4rilt(4)
  			xfeat(18) = 0.0
c
c RI = 65 kt / 72 h
c
          elseif (j .eq. 8) then
c
  			xfeat(3) = rlatrilt(5)
  			xfeat(4) = 0.0
  			xfeat(5) = 0.0
  			xfeat(6) = rhlorilt(5)
  			xfeat(7) = 0.0
  			xfeat(8) = 0.0
  			xfeat(9) = shgcrilt(5)
  			xfeat(10) = 0.0
  			xfeat(11) = potrilt(5)
  			xfeat(12) = 0.0
  			xfeat(13) = sir6ri
  			xfeat(14) = 0.0
  			xfeat(15) = 0.0
  			xfeat(16) = 0.0
  			xfeat(17) = 0.0
  			xfeat(18) = tpw6rilt(5)
c			
          endif
c
          sumExp = 0.0
c
c Add together coefficients for the logistic regression equation
c
          do i = 1, Nthrss
               if (int(xfeat(i)) .ne. 9999 .and. 
     +           int(sumExp) .ne. 9999) then
                   sumExp = sumExp + xfeat(i) * coef(j, i)
               else
                   sumExp = 9999.
               endif
          enddo
c
c Compute the logistic regression based probability
c
          if (int(sumExp) .ne. 9999) then
            probrilog(j) = 100. / (1. + exp( -sumExp) )
          else
            probrilog(j) = 999.
          endif
      enddo
c
c Make missing values consistent with iships.
      if (int(sir5ri) .eq. 9999) sir5ri = 999.
      if (int(sir6ri) .eq. 9999) sir6ri = 999.
      if (int(sir10ri) .eq. 9999) sir10ri = 999.
      if (int(sir12ri) .eq. 9999) sir12ri = 999.	
c
      return
      end
c