c subroutine rapidbaye_2018(perri, ilat, inohc, iepos, iepss, + irhhi, id200, ishdc, ishgc, sbtri, sir5ri, + sir7ri, sir9ri, sir10ri, sir11ri, sir12ri, + sir14ri, vsst, sst, vmx, itpw6, probribay) c c subroutine rapidbaye.f c c v2.2.2 c c Updated: 22 Feb 2018 for 2018 season (CR) c 04 Mar 2018 to put ioper parameter back in for WCOSS (MD) c 08 Apr 2020 removed ioper flag (SS) c 11 Apr 2020 updated 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 Bayesian-based c model (Rozoff and Kossin 2011; Wea. Forecasting). c This model provides estimates of the probability of RI c for the 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 rlatriX : latitude averaged from t=0 to t=X h c nohcriX : NCODA heat content (kj/cm2) averaged from t=0 to t=X h c eposriX : The average theta_e difference between a parcel lifted from c the surface and its environment (200-800 km average) vs. c time. Only positive differences are included in average 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 rhhiriX : 500-250 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 tpw6X : TPW standard deviation 400-600 km averaged from c T=0 to T=X h. c sbtri : std. dev. of the 0-200 km GOES channel 4 BT at t=0 h c sir5ri : std. dev. of the 100-300 km GOES channel 4 BT at t=0 h c sir7ri : % area r=50-200 km of GOES ch 4 BT < -20 C c sir9ri : % area r=50-200 km of GOES ch 4 BT < -40 C c sir10ri : % area r=50-200 km of GOES ch 4 BT < -50 C c sir11ri : % area r=50-200 km of GOES ch 4 BT < -60 C c sir12ri : max BT from 0 to 30 km radius (deg C*10) c sir14ri : radius of maximum BT (km) 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 probribay : 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 x1X : non-RI bins c x2X : RI bins c p1X : non-RI pdf c p2X : RI pdf c pclim : climatological probability of RI for each RI defn c c ------------------------------------------------------------------------------ c parameter (Mft = 28) c parameter (Nindx = 8) c dimension Nthrss(Nindx) data Nthrss/7, 8, 9, 9, 9, 10, 10, 9/ c parameter (Nlt = 5) c dimension xfeat(10) c real x11(7, 100), x21(7, 100) real p11(7, 100), p21(7, 100) real x12(8, 100), x22(8, 100) real p12(8, 100), p22(8, 100) real x13(9, 100), x23(9, 100) real p13(9, 100), p23(9, 100) real x14(9, 100), x24(9, 100) real p14(9, 100), p24(9, 100) real x15(9, 100), x25(9, 100) real p15(9, 100), p25(9, 100) real x16(10, 100), x26(10, 100) real p16(10, 100), p26(10, 100) real x17(10, 100), x27(10, 100) real p17(10, 100), p27(10, 100) real x18(9, 100), x28(9, 100) real p18(9, 100), p28(9, 100) c c Climatological probabilities of RI for each RI defn c real pclim(Nindx) data pclim /0.0653,0.1258,0.0863,0.0620,0.0426, + 0.0656,0.0573,0.0409/ c c input c dimension ilat(-2:Mft), id200(0:Mft), ishdc(0:Mft), + ishgc(0:Mft), vsst(0:Mft), sst(0:Mft), iepos(0:Mft), + iepss(0:Mft), inohc(0:Mft), irhhi(0:Mft), itpw6(0:Mft) c c output probabilities c dimension probribay(Nindx) c c internal variables c integer i, j, iskip integer nlat, nd200, nshdc, nshgc, npot, ntpw6, + nepss, nnohc, nepos, nrhhi real sumlat, sumd200, sumshdc, sumpot, sumtpw6, + sumepss, sumnohc, sumepos, sumrhhi c dimension d200rilt(nlt), shdcrilt(nlt), shgcrilt(nlt), + potrilt(nlt), tpw6rilt(nlt), rlatrilt(nlt), epssrilt(nlt), + rnohcrilt(nlt), eposrilt(nlt), rhhirilt(nlt) c c c PDF file name variables character *256 fnpdf, coef_location c c internal data c integer ideltat, itimeinc parameter (ideltat = 12, itimeinc = 6) c c ------------------------------------------------------------------------------ c c Make missing values consistent with subroutine here. if (int(sbtri) .eq. 999) sbtri = 9999. if (int(sir5ri) .eq. 999) sir5ri = 9999. if (int(sir7ri) .eq. 999) sir7ri = 9999. if (int(sir9ri) .eq. 999) sir9ri = 9999. if (int(sir10ri) .eq. 999) sir10ri = 9999. if (int(sir11ri) .eq. 999) sir11ri = 9999. if (int(sir12ri) .eq. 999) sir12ri = 9999. if (int(sir14ri) .eq. 999) sir14ri = 9999. c c Compute mean predictors values for each of the forecast c lead times (12, 24, 36, 48 h, 72 h) c do j = 1, nlt nft = (j * ideltat) / itimeinc if (j .eq. 5) nft = nft + 2 nlat = 0 nd200 = 0 nshdc = 0 nshgc = 0 npot = 0 nepss = 0 nnohc = 0 nepos = 0 nrhhi = 0 ntpw6 = 0 c sumlat = 0. sumd200 = 0. sumshdc = 0. sumshgc = 0. sumpot = 0. sumepss = 0. sumnohc = 0. sumepos = 0. sumrhhi = 0. sumtpw6 = 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 (id200(i) .lt. 9999) then nd200 = nd200 + 1 sumd200 = sumd200 + float(id200(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 (inohc(i) .lt. 9999) then nnohc = nnohc + 1 sumnohc = sumnohc + float(inohc(i)) endif c if (iepos(i) .lt. 9999) then nepos = nepos + 1 sumepos = sumepos + 0.1 * float(iepos(i)) endif c if (irhhi(i) .lt. 9999) then nrhhi = nrhhi + 1 sumrhhi = sumrhhi + float(irhhi(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 (nd200 .ge. nmin) then d200rilt(j) = sumd200 / float(nd200) else d200rilt(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 (nnohc .ge. nmin) then rnohcrilt(j) = sumnohc / float(nnohc) else rnohcrilt(j) = 9999. endif c if (nepos .ge. nmin) then eposrilt(j) = sumepos / float(nepos) else eposrilt(j) = 9999. endif c if (nrhhi .ge. nmin) then rhhirilt(j) = sumrhhi / float(nrhhi) else rhhirilt(j) = 9999. endif c if (ntpw6 .ge. nmin) then tpw6rilt(j) = sumtpw6 / float(ntpw6) else tpw6rilt(j) = 9999. endif c enddo c c Open file of empirical PDFS and read in data c call getenv( "SHIPS_COEF",coef_location ) fnpdf = trim( coef_location ) //'pdfs_RI_epac_2018.dat' c open(unit = 10, file = fnpdf, status = 'old') read(10,*) x11 read(10,*) x21 read(10,*) p11 read(10,*) p21 read(10,*) x12 read(10,*) x22 read(10,*) p12 read(10,*) p22 read(10,*) x13 read(10,*) x23 read(10,*) p13 read(10,*) p23 read(10,*) x14 read(10,*) x24 read(10,*) p14 read(10,*) p24 read(10,*) x15 read(10,*) x25 read(10,*) p15 read(10,*) p25 read(10,*) x16 read(10,*) x26 read(10,*) p16 read(10,*) p26 read(10,*) x17 read(10,*) x27 read(10,*) p17 read(10,*) p27 read(10,*) x18 read(10,*) x28 read(10,*) p18 read(10,*) p28 close(10) c do j = 1, Nindx c c RI = 20 kt / 12 h c if (j .eq. 1) then xfeat(1) = perri xfeat(2) = rlatrilt(1) xfeat(3) = rnohcrilt(1) xfeat(4) = d200rilt(1) xfeat(5) = shdcrilt(1) xfeat(6) = sir10ri xfeat(7) = potrilt(1) c c Call primary model function c iskip = 0 do i = 1, Nthrss(j) if (int(xfeat(i)) .eq. 9999) then iskip = 1 endif enddo if (iskip .eq. 0) then probribay(j) = postp2rie(xfeat, Nthrss(j), + pclim(j), x11, x21, p11, p21) else probribay(j) = 999. endif c c RI = 25 kt / 24 h c elseif (j .eq. 2) then xfeat(1) = perri xfeat(2) = rlatrilt(2) xfeat(3) = rnohcrilt(2) xfeat(4) = epssrilt(2) xfeat(5) = shdcrilt(2) xfeat(6) = sir5ri xfeat(7) = sir10ri xfeat(8) = potrilt(2) iskip = 0 do i = 1, Nthrss(j) if (int(xfeat(i)) .eq. 9999) then iskip = 1 endif enddo if (iskip .eq. 0) then probribay(j) = postp2rie(xfeat, Nthrss(j), + pclim(j), x12, x22, p12, p22) else probribay(j) = 999. endif c c RI = 30 kt / 24 h c elseif (j .eq. 3) then xfeat(1) = perri xfeat(2) = rlatrilt(2) xfeat(3) = rnohcrilt(2) xfeat(4) = epssrilt(2) xfeat(5) = d200rilt(2) xfeat(6) = shdcrilt(2) xfeat(7) = sir5ri xfeat(8) = sir10ri xfeat(9) = potrilt(2) iskip = 0 do i = 1, Nthrss(j) if (int(xfeat(i)) .eq. 9999) then iskip = 1 endif enddo if (iskip .eq. 0) then probribay(j) = postp2rie(xfeat, Nthrss(j), + pclim(j), x13, x23, p13, p23) else probribay(j) = 999. endif c c RI = 35 kt / 24 h c elseif (j .eq. 4) then xfeat(1) = perri xfeat(2) = rlatrilt(2) xfeat(3) = rnohcrilt(2) xfeat(4) = epssrilt(2) xfeat(5) = d200rilt(2) xfeat(6) = shdcrilt(2) xfeat(7) = sir5ri xfeat(8) = sir10ri xfeat(9) = potrilt(2) iskip = 0 do i = 1, Nthrss(j) if (int(xfeat(i)) .eq. 9999) then iskip = 1 endif enddo if (iskip .eq. 0) then probribay(j) = postp2rie(xfeat, Nthrss(j), + pclim(j), x14, x24, p14, p24) else probribay(j) = 999. endif c c RI = 40 kt / 24 h c elseif (j .eq. 5) then xfeat(1) = perri xfeat(2) = rlatrilt(2) xfeat(3) = rnohcrilt(2) xfeat(4) = epssrilt(2) xfeat(5) = shgcrilt(2) xfeat(6) = sir5ri xfeat(7) = sir11ri xfeat(8) = sir14ri xfeat(9) = potrilt(2) iskip = 0 do i = 1, Nthrss(j) if (int(xfeat(i)) .eq. 9999) then iskip = 1 endif enddo if (iskip .eq. 0) then probribay(j) = postp2rie(xfeat, Nthrss(j), + pclim(j), x15, x25, p15, p25) else probribay(j) = 999. endif c c RI = 45 kt / 36 h c elseif (j .eq. 6) then xfeat(1) = perri xfeat(2) = rlatrilt(3) xfeat(3) = rnohcrilt(3) xfeat(4) = epssrilt(3) xfeat(5) = rhhirilt(3) xfeat(6) = shdcrilt(3) xfeat(7) = sir5ri xfeat(8) = sir9ri xfeat(9) = sir12ri xfeat(10) = potrilt(3) iskip = 0 do i = 1, Nthrss(j) if (int(xfeat(i)) .eq. 9999) then iskip = 1 endif enddo if (iskip .eq. 0) then probribay(j) = postp2rie(xfeat, Nthrss(j), + pclim(j), x16, x26, p16, p26) else probribay(j) = 999. endif c c RI = 55 kt / 48 h c elseif (j .eq. 7) then xfeat(1) = perri xfeat(2) = rlatrilt(4) xfeat(3) = rnohcrilt(4) xfeat(4) = epssrilt(4) xfeat(5) = rhhirilt(4) xfeat(6) = shdcrilt(4) xfeat(7) = sbtri xfeat(8) = sir10ri xfeat(9) = sir12ri xfeat(10) = potrilt(4) iskip = 0 do i = 1, Nthrss(j) if (int(xfeat(i)) .eq. 9999) then iskip = 1 endif enddo if (iskip .eq. 0) then probribay(j) = postp2rie(xfeat, Nthrss(j), + pclim(j), x17, x27, p17, p27) else probribay(j) = 999. endif c c RI = 65 kt / 72 h c elseif (j .eq. 8) then xfeat(1) = perri xfeat(2) = rnohcrilt(5) xfeat(3) = eposrilt(5) xfeat(4) = rhhirilt(5) xfeat(5) = d200rilt(5) xfeat(6) = shgcrilt(5) xfeat(7) = tpw6rilt(5) xfeat(8) = sir7ri xfeat(9) = potrilt(5) iskip = 0 do i = 1, Nthrss(j) if (int(xfeat(i)) .eq. 9999) then iskip = 1 endif enddo if (iskip .eq. 0) then probribay(j) = postp2rie(xfeat, Nthrss(j), + pclim(j), x18, x28, p18, p28) else probribay(j) = 999. endif endif enddo c c Make missing values consistent with iships. if (int(sbtri) .eq. 9999) sbtri = 999. if (int(sir5ri) .eq. 9999) sir5ri = 999. if (int(sir7ri) .eq. 9999) sir7ri = 999. if (int(sir9ri) .eq. 9999) sir9ri = 999. if (int(sir10ri) .eq. 9999) sir10ri = 999. if (int(sir11ri) .eq. 9999) sir11ri = 999. if (int(sir12ri) .eq. 9999) sir12ri = 999. if (int(sir14ri) .eq. 9999) sir14ri = 999. c return end c c This is the naive Bayesian probabilistic model of c Rozoff and Kossin (WAF; 2011) for the NATL c function postp2rie(xfeat, nfeat, prior2, x1, x2, p1, p2) c c Christopher Rozoff (CIMSS), 2012 c (Code originally written by James P. Kossin) c real xfeat(10) integer j real x1(nfeat, 100), x2(nfeat, 100) real p1(nfeat, 100), p2(nfeat, 100) real x1s(100), x2s(100), p1s(100), p2s(100) real prior2 c prior1 = 1. - prior2 ! prior probability of non-RI c prod1 = 1. prod2 = 1. do ifeat = 1, nfeat do ix = 1, 100 x1s(ix) = x1(ifeat, ix) x2s(ix) = x2(ifeat, ix) p1s(ix) = p1(ifeat, ix) p2s(ix) = p2(ifeat, ix) enddo c c Call the straddle subroutine inside of PrSEFoNe.f c call straddle(xfeat(ifeat), x1s, ib1, ib2) dx1 = abs(x1s(ib1) - x1s(ib2)) c c Integrate under PDF (trapezoidal) c prob = 0.5 * (p1s(ib1) + p1s(ib2)) * dx1 c c naive assumption c prod1 = prod1 * prob c call straddle(xfeat(ifeat), x2s, ib1, ib2) dx2 = abs(x2s(ib1) - x2s(ib2)) prob = 0.5 * (p2s(ib1) + p2s(ib2)) * dx2 prod2 = prod2 * prob c enddo c den = prior1 * prod1 + prior2 * prod2 postp2rie = (prior2 / den) * prod2 * 100 c return end