c
      subroutine rapidloga_2018( perri, vsstc, sst, vmx, inohc,
     +                 iu200, ienss, id200, ishgc, ishrg, sir2ri,
     +                 sbtri, sir5ri, sir10ri, sir12ri, sir13ri,
     +                 sir14ri, itpw4, itpw19, pc2ri, probrilog)
c
c     subroutine rapidloga.f
c
c     v2.2.2
c
c     Updated: 21 Feb 2018 for 2018 season (CR)
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 Atlantic 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, and,
c          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     rsstriX : Reynolds sea surface temperature (C) averaged from t=0 to 
c             t=X h
c     rnohciX  : NCODA Ocean heat content (kj/cm2) averaged from t=0 to t=X h
c     u200riX : 200-hPa zonal wind (kt) vs time (r=200-800 km) averaged from
c	      T=0 to T=X h.
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     d200riX : 200-hPa divergence vs time (r=0-1000 km)
c     shgcriX : 850-200-hPa shear magnitude with vortex removed and averaged from
c             0-500 km relative to 850 hPa vortex center
c     shrgriX : Generalized 850-200 mb shear magintude (kt*10) vs time (takes
c             into account all levels) averaged from t=0 to t=X h.
c     sir2ri  : Average GOES ch 4 BT, r = 0-200 km
c     sbtri  : Stan. Dev. of GOES BT (deg C*10), r=0-200 km
c     sir5ri  : Stan. Dev. of GOES BT (deg C*10), r=100-300 km
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     sir13ri : avg BT from 0 to 30 km radius (deg C*10)
c     sir14ri : radius of maximum BT (km)
c     tpw4X  : TPW standard deviation 200-400 km averaged from
c	      T=0 to T=X h.
c     tpw19X : %TPW less than 45 mm, r=0=500 km averaged from
c	      T=0 to T=X h.
c     pc2ri   : IR PC 2
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
c     subroutine output:
c 
c     probrilog : the forecasted probabilities of RI for each of the 7
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 = 20)
c
      parameter (Nlt = 5)
c
c input
c
      dimension vsstc(0:Mft), sst(0:Mft), inohc(0:Mft), 
     +          iu200(0:Mft), ienss(0:Mft), id200(0:Mft),
     +          ishgc(0:Mft), ishrg(0:Mft), itpw4(0:Mft),
     +			itpw19(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)
     * /-3.4605983e+00,  6.1652376e-02,  0.0000000e+00,
     *   5.9769695e-03,  0.0000000e+00,  0.0000000e+00,
     *   0.0000000e+00,  0.0000000e+00, -4.2456800e-02,
     *   0.0000000e+00,  0.0000000e+00, -5.1923660e-02,
     *   1.3634771e-02, -1.0714040e-02,  0.0000000e+00,
     *   0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
     *  -3.0342484e-03,  1.2117134e-02/
c 25 kt / 24-h RI threshold
      data (coef(2,i), i = 1, Nthrss)
     * /-4.5404326e+00,  6.1904925e-02,  1.1090245e-01,
     *   0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
     *   5.0707953e-03, -8.5959640e-02,  0.0000000e+00,
     *   0.0000000e+00,  0.0000000e+00, -4.7119729e-02,
     *   0.0000000e+00, -1.8056111e-02,  0.0000000e+00,
     *   0.0000000e+00, -1.4151865e-02,  0.0000000e+00,
     *  -2.8349312e-03,  2.1098805e-02/
c 30 kt / 24-h RI threshold
      data (coef(3,i), i = 1, Nthrss)
     * /-6.7835244e+00,  6.7081963e-02,  1.6525489e-01,
     *   0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
     *   0.0000000e+00, -8.9281624e-02,  0.0000000e+00,
     *   0.0000000e+00,  0.0000000e+00, -6.3079049e-02,
     *   0.0000000e+00, -2.0169315e-02,  0.0000000e+00,
     *   0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
     *  -3.1024811e-03,  2.0980751e-02/
c 35 kt / 24-h RI threshold
      data (coef(4,i), i = 1, Nthrss)
     * /-9.9937755e+00,  6.5041929e-02,  2.5014488e-01,
     *   0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
     *   0.0000000e+00,  0.0000000e+00, -8.5334801e-02,
     *   0.0000000e+00,  0.0000000e+00, -8.8898976e-02,
     *   0.0000000e+00, -2.8334458e-02,  0.0000000e+00,
     *   0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
     *   0.0000000e+00,  2.0857206e-02/
c 40 kt / 24-h RI threshold
      data (coef(5,i), i = 1, Nthrss)
     * /-1.1610919e+01,  4.8080750e-02,  3.1785327e-01,
     *   0.0000000e+00, -2.2936840e-02,  0.0000000e+00,
     *   0.0000000e+00,  0.0000000e+00, -7.3979077e-02,
     *   0.0000000e+00,  0.0000000e+00, -8.4686754e-02,
     *   0.0000000e+00, -3.3872699e-02,  0.0000000e+00,
     *   0.0000000e+00, -2.5009715e-02,  0.0000000e+00,
     *   0.0000000e+00,  1.4848489e-02/
c 45 kt / 36-h RI threshold
      data (coef(6,i), i = 1, Nthrss)	  
     * /-2.4074890e+00,  5.9129735e-02,  0.0000000e+00,
     *   0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
     *   0.0000000e+00,  0.0000000e+00, -1.1713109e-01,
     *   0.0000000e+00, -4.1537587e-02,  0.0000000e+00,
     *   0.0000000e+00, -2.3313761e-02,  0.0000000e+00,
     *   0.0000000e+00, -2.4790999e-02,  0.0000000e+00,
     *   0.0000000e+00,  2.7081385e-02/
c 55 kt / 48-h RI threshold
      data (coef(7,i), i = 1, Nthrss)
     * /-2.4748353e+00,  4.4528069e-02,  0.0000000e+00,
     *   0.0000000e+00,  0.0000000e+00, -1.5194748e-01,
     *   0.0000000e+00,  0.0000000e+00, -1.4132625e-01,
     *  -2.9610428e-02,  0.0000000e+00,  0.0000000e+00,
     *   0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
     *   1.0837965e-02, -4.3072008e-02,  0.0000000e+00,
     *   0.0000000e+00,  3.5455484e-02/
c 65 kt / 72-h RI threshold
      data (coef(8,i), i = 1, Nthrss)
     * / 3.6949588e+00,  5.2594102e-02, -1.6573095e-01,
     *   0.0000000e+00,  0.0000000e+00, -6.7651930e-01,
     *   0.0000000e+00,  0.0000000e+00, -1.5979605e-01,
     *   0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
     *   0.0000000e+00,  0.0000000e+00, -1.4703265e-02,
     *   0.0000000e+00, -6.4524510e-02,  0.0000000e+00,
     *   0.0000000e+00,  5.6296395e-02/
c
      dimension probrilog(Nindx)
c
      dimension xfeat(Nthrss)
c
c  internal variables
c
      integer i, j
      real sumExp
	  integer nsst, nnohc, nu200, nenss, nd200, nshgc,
     +        nshrg, ntpw4, ntpw19, npot
	  real sumsst, sumnohc, sumu200, sumenss,
     +     sumd200, sumshgc, sumshrg, sumtpw4,
     +	   sumtpw19
c
	  dimension sstrilt(nlt), rnohcrilt(nlt), u200rilt(nlt), 
     +          enssrilt(nlt), d200rilt(nlt), shgcrilt(nlt), 
     +          shrgrilt(nlt), tpw4rilt(nlt), tpw19rilt(nlt), 
     +          potrilt(nlt)
c
c ------------------------------------------------------------------------------
c
c Make missing values consistent with subroutine here.
      if (int(sir2ri) .eq. 999) sir2ri = 9999.
	  if (int(sbtri) .eq. 999) sbtri = 9999.
      if (int(sir5ri) .eq. 999) sir5ri = 9999.
      if (int(sir10ri) .eq. 999) sir10ri = 9999.
      if (int(sir12ri) .eq. 999) sir12ri = 9999.
	  if (int(sir13ri) .eq. 999) sir13ri = 9999.
      if (int(sir14ri) .eq. 999) sir14ri = 9999.
      if (int(pc2ri) .eq. 999) pc2ri = 9999.
c
c Compute mean predictors values for each of the forecast
c lead times (12, 24, 36, and 48 h)
c
      do j = 1, nlt
          nft = (j * ideltat) / itimeinc
          if (j .eq. 5) nft = nft + 2
		  nsst = 0
		  nnohc = 0
          nu200 = 0
          nenss = 0
          nd200 = 0
          nshgc = 0
          nshrs = 0
          nshrg = 0
		  ntpw4 = 0
		  ntpw19 = 0
          npot = 0
c
          sumsst = 0.
          sumnohc = 0.
          sumu200 = 0.
          sumenss = 0.
          sumd200 = 0.
          sumshgc = 0.
          sumshrs = 0.
          sumshrg = 0.
		  sumtpw4 = 0.
		  sumtpw19 = 0.
          sumpot = 0.
c
          do i = 0, nft
c
            if (vsstc(i) .lt. 200.0) then
              nsst = nsst + 1
              sumsst = sumsst + sst(i)
            endif
c
            if (inohc(i) .lt. 9999) then
              nnohc = nnohc + 1
              sumnohc = sumnohc + float(inohc(i))
            endif
c
            if (iu200(i) .lt. 9999) then
              nu200 = nu200 + 1
              sumu200 = sumu200 + 0.1 * float(iu200(i))
            endif
c
            if (ienss(i) .lt. 9999) then
              nenss = nenss + 1
              sumenss = sumenss + 0.1 * float(ienss(i))
            endif
c
            if (id200(i) .lt. 9999) then
              nd200 = nd200 + 1
              sumd200 = sumd200 + float(id200(i))
            endif
c
            if (ishgc(i) .lt. 9999) then
              nshgc = nshgc + 1
              sumshgc = sumshgc + 0.1 * float(ishgc(i))
            endif
c
            if (ishrg(i) .lt. 9999) then
              nshrg = nshrg + 1
              sumshrg = sumshrg + 0.1 * float(ishrg(i))
            endif
c
            if (vsstc(i) .lt. 200 .and. vmx .lt. 200) then
              npot = npot + 1
              sumpot = sumpot + (vsstc(i) - vmx)
            endif
c
            if (itpw4(i) .lt. 9999) then
              ntpw4 = ntpw4 + 1
              sumtpw4 = sumtpw4 + float(itpw4(i))
            endif
c
            if (itpw19(i) .lt. 9999) then
              ntpw19 = ntpw19 + 1
              sumtpw19 = sumtpw19 + float(itpw19(i))
            endif
c
          enddo
c
          nmin = nft / 2 + 1
c
          if (nsst .ge. nmin) then
            sstrilt(j) = sumsst / float(nsst)
          else
            sstrilt(j) = 9999.
          endif
c
          if (nnohc .ge. nmin) then
            rnohcrilt(j) = sumnohc / float(nnohc)
          else
            rnohcrilt(j) = 9999.
          endif
c
          if (nu200 .ge. nmin) then
            u200rilt(j) = sumu200 / float(nu200)
          else
            u200rilt(j) = 9999.
          endif
c
          if (nenss .ge. nmin) then
            enssrilt(j) = sumenss / float(nenss)
          else
            enssrilt(j) = 9999.
          endif
c
          if (nd200 .ge. nmin) then
            d200rilt(j) = sumd200 / float(nd200)
          else
            d200rilt(j) = 9999.
          endif
c
          if (nshgc .ge. nmin) then
            shgcrilt(j) = sumshgc / float(nshgc)
          else
            shgcrilt(j) = 9999.
          endif
c
          if (nshrg .ge. nmin) then
            shrgrilt(j) = sumshrg / float(nshrg)
          else
            shrgrilt(j) = 9999.
          endif
c
          if (ntpw4 .ge. nmin) then
            tpw4rilt(j) = sumtpw4 / float(ntpw4)
          else
            tpw4rilt(j) = 9999.
          endif
c
          if (ntpw19 .ge. nmin) then
            tpw19rilt(j) = sumtpw19 / float(ntpw19)
          else
            tpw19rilt(j) = 9999.
          endif
c
          if (npot .ge. nmin) then
            potrilt(j) = sumpot / float(npot)
          else
            potrilt(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) = 0.0
			xfeat(4) = rnohcrilt(1)
			xfeat(5) = 0.0
			xfeat(6) = 0.0
			xfeat(7) = 0.0
			xfeat(8) = 0.0
			xfeat(9) = shrgrilt(1)
			xfeat(10) = 0.0
			xfeat(11) = 0.0
			xfeat(12) = sir5ri
			xfeat(13) = sir10ri
			xfeat(14) = sir12ri
			xfeat(15) = 0.0
			xfeat(16) = 0.0
			xfeat(17) = 0.0
			xfeat(18) = 0.0
			xfeat(19) = pc2ri
			xfeat(20) = potrilt(1)
c
c RI = 25 kt / 24 h
c
          elseif (j .eq. 2) then
c
            xfeat(3) = sstrilt(2)
			xfeat(4) = 0.0
			xfeat(5) = 0.0
			xfeat(6) = 0.0
			xfeat(7) = d200rilt(2)
			xfeat(8) = shgcrilt(2)
			xfeat(9) = 0.0
			xfeat(10) = 0.0
			xfeat(11) = 0.0
			xfeat(12) = sir5ri
			xfeat(13) = 0.0
			xfeat(14) = sir12ri
			xfeat(15) = 0.0
			xfeat(16) = 0.0
			xfeat(17) = tpw4rilt(2)
			xfeat(18) = 0.0
			xfeat(19) = pc2ri
			xfeat(20) = potrilt(2)

c
c RI = 30 kt / 24 h
c
          elseif (j .eq. 3) then
c
            xfeat(3) = sstrilt(2)
            xfeat(4) = 0.0
            xfeat(5) = 0.0
            xfeat(6) = 0.0
            xfeat(7) = 0.0
            xfeat(8) = shgcrilt(2)
            xfeat(9) = 0.0
            xfeat(10) = 0.0
            xfeat(11) = 0.0
            xfeat(12) = sir5ri
            xfeat(13) = 0.0
            xfeat(14) = sir12ri
            xfeat(15) = 0.0
            xfeat(16) = 0.0
            xfeat(17) = 0.0
            xfeat(18) = 0.0
            xfeat(19) = pc2ri
            xfeat(20) = potrilt(2)
c
c RI = 35 kt / 24 h
c
          elseif (j .eq. 4) then
c
            xfeat(3) = sstrilt(2)
            xfeat(4) = 0.0
            xfeat(5) = 0.0
            xfeat(6) = 0.0
            xfeat(7) = 0.0
            xfeat(8) = 0.0
            xfeat(9) = shrgrilt(2)
            xfeat(10) = 0.0
            xfeat(11) = 0.0
            xfeat(12) = sir5ri
            xfeat(13) = 0.0
            xfeat(14) = sir12ri
            xfeat(15) = 0.0
            xfeat(16) = 0.0
            xfeat(17) = 0.0
            xfeat(18) = 0.0
            xfeat(19) = 0.0
            xfeat(20) = potrilt(2)
c
c RI = 40 kt / 24 h
c
          elseif (j .eq. 5) then
c
            xfeat(3) = sstrilt(2)
            xfeat(4) = 0.0
            xfeat(5) = u200rilt(2)
            xfeat(6) = 0.0
            xfeat(7) = 0.0
            xfeat(8) = 0.0
            xfeat(9) = shrgrilt(2)
            xfeat(10) = 0.0
            xfeat(11) = 0.0
            xfeat(12) = sir5ri
            xfeat(13) = 0.0
            xfeat(14) = sir12ri
            xfeat(15) = 0.0
            xfeat(16) = 0.0
            xfeat(17) = tpw4rilt(2)
            xfeat(18) = 0.0
            xfeat(19) = 0.0
            xfeat(20) = potrilt(2)
c
c RI = 45 kt / 36 h
c
          elseif (j .eq. 6) then
c
            xfeat(3) = 0.0
            xfeat(4) = 0.0
            xfeat(5) = 0.0
            xfeat(6) = 0.0
            xfeat(7) = 0.0
            xfeat(8) = 0.0
            xfeat(9) = shrgrilt(3)
            xfeat(10) = 0.0
            xfeat(11) = sbtri
            xfeat(12) = 0.0
            xfeat(13) = 0.0
            xfeat(14) = sir12ri
            xfeat(15) = 0.0
            xfeat(16) = 0.0
            xfeat(17) = tpw4rilt(3)
            xfeat(18) = 0.0
            xfeat(19) = 0.0
            xfeat(20) = potrilt(3)
c
c RI = 55 kt / 48 h
c
          elseif (j .eq. 7) then
c
            xfeat(3) = 0.0
            xfeat(4) = 0.0
            xfeat(5) = 0.0
            xfeat(6) = enssrilt(4)
            xfeat(7) = 0.0
            xfeat(8) = 0.0
            xfeat(9) = shrgrilt(4)
            xfeat(10) = sir2ri
            xfeat(11) = 0.0
            xfeat(12) = 0.0
            xfeat(13) = 0.0
            xfeat(14) = 0.0
            xfeat(15) = 0.0
            xfeat(16) = sir14ri
            xfeat(17) = tpw4rilt(4)
            xfeat(18) = 0.0
            xfeat(19) = 0.0  
            xfeat(20) = potrilt(4)
c
c RI = 65 kt / 72 h
c
          elseif (j .eq. 8) then
c
	        xfeat(3) = sstrilt(5)
            xfeat(4) = 0.0
            xfeat(5) = 0.0
            xfeat(6) = enssrilt(5)
            xfeat(7) = 0.0
            xfeat(8) = 0.0
            xfeat(9) = shrgrilt(5)
            xfeat(10) = 0.0
            xfeat(11) = 0.0
            xfeat(12) = 0.0
            xfeat(13) = 0.0
            xfeat(14) = 0.0
            xfeat(15) = sir13ri
            xfeat(16) = 0.0
            xfeat(17) = tpw4rilt(5)
            xfeat(18) = 0.0
            xfeat(19) = 0.0
            xfeat(20) = potrilt(5)
          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
      if (int(sir2ri) .eq. 9999) sir2ri = 999.
	  if (int(sbtri) .eq. 9999) sbtri = 999.
      if (int(sir5ri) .eq. 9999) sir5ri = 999.
      if (int(sir10ri) .eq. 9999) sir10ri = 999.
      if (int(sir12ri) .eq. 9999) sir12ri = 999.
	  if (int(sir13ri) .eq. 9999) sir13ri = 999.
      if (int(sir14ri) .eq. 9999) sir14ri = 999.
      if (int(pc2ri) .eq. 9999) pc2ri = 999.
c
      return
      end
c