subroutine ric2o(mft,ibasin,dvm12,vm00,vmpi,ishdc,id200,igoes, + irhcn,irhlo,itwac,imiss,rmiss, + sname,natcf8,ismon,isday,isyr,istime, + luout,probric2,potmin_final,ierr) c The routine uses a discrimant analysis to determine the probability c of average intensification or rapid intensification. This code was developed c Aug 2013 as a simplification of the ric3 code that also includes a c rapid weakening class. c c This version is for operational runs for global tropical cyclones c c Version 3.0.3 c c Modified 7/26/2014 - MD ver 2.1.0 to skip writing if luout < 0 for ver2.1 c Modified 5/04/2015 - MD ver 2.2.0 to add option to set procric2 to zero if c if potential intensity is too low to support RI c Modified 4/11/2017 - MD ver 3.0.0 to add two new predictors c (vmax and vmax**2). Also, set potmin for c 48 hr RII for 2017 season c Modified 3/23/2019 - GC ver 3.0.1 to rename potmin to c potmin_final, and make potmin_final c a calling argument c Modified 4/22/2019 - GC ver 3.0.2 divide Vmax*82 predictor by c vmax2_scale to match what is used in c rii2.f and make Vmax**2 values similar c to Vmax c Modified 4/4/2020 - SS removed ioper flag, increased mftl for 7day c c Passed variables dimension vmpi(0:mft),ishdc(0:mft),id200(0:mft),igoes(0:mft) dimension irhcn(0:mft),irhlo(0:mft),itwac(0:mft) character *8 natcf8 character *10 sname c c Local variables parameter (mvda=50,mgda=10,mftl=28) dimension dfcoef(0:mvda,mgda,0:mftl) dimension dmean(0:mvda,0:mftl),dstd(0:mvda,0:mftl) dimension dfinp(0:mvda),dfinpc(0:mvda),dfinpnl(0:mvda) dimension dfval(mgda),dfvalc(mgda) character *6 dlab(0:mvda) character *25 plab(0:mvda) c dimension shdc(0:mftl),d200(0:mftl),rhcn(0:mftl) dimension rhlo(0:mftl),twac(0:mftl) c character *7 glab character *30 cbasin character *19 fncoef,fnprob character *256 fncol character *256 coef_location c c Variables for converting df values to probabilities parameter (mx=101,mxg=2) dimension iptime(0:mftl),nx(0:mft) dimension dx(0:mftl) dimension xbar(0:mftl),xstd(0:mftl),x1(0:mftl),x2(0:mftl) dimension dfprob(0:mx,mxg,0:mftl) dimension gprob(mxg), gprobc(mxg), gprobnl(mxg) dimension igprob(mxg),igprobc(mxg),igprobnl(mxg) dimension gcont(0:mvda) c vmax**2 predictor will be divided by vmax2_scale to c match what is used by rii2.f and normalize Vmax**2 values to be c more similar to Vmax real vmax2_scale parameter(vmax2_scale = 70.) c ierr=0 c c Set output probability to missing for now probric2 = 999. c c Set minimum required potential intensity. For 2017 NHC model c 55 kt change in 48 hr is used for RI definition. The threshold c is set to 5 kt less that than. c potmin = 50.0 c ccccc potmin is replaced by potmin_final c ccccc potmin_final is set in the calling subroutine c c Specify predictor labels for output file plab( 0) = 'Climatology ' plab( 1) = 'SST Potential ' plab( 2) = '850-200 hPa Shear ' plab( 3) = '200 hPa Divergence ' plab( 4) = 'Persistence ' plab( 5) = 'Sat. IR Cold Pixels ' plab( 6) = 'Sat. IR asymmetry ' plab( 7) = 'Ocean Heat Content ' plab( 8) = '850-700 hPa RH ' plab( 9) = 'Model Tang Wind Tendency ' plab(10) = 'vmax ' plab(11) = 'vmax*vmax ' c Convert input to real as needed do k=0,mft if (ishdc(k) .lt. imiss) then shdc(k) = 0.1*float(ishdc(k)) else shdc(k) = rmiss endif c if (id200(k) .lt. imiss) then d200(k) = float(id200(k)) else d200(k) = rmiss endif c if (irhcn(k) .lt. imiss) then rhcn(k) = float(irhcn(k)) else rhcn(k) = rmiss endif c if (irhlo(k) .lt. imiss) then rhlo(k) = float(irhlo(k)) else rhlo(k) = rmiss endif c if (itwac(k) .lt. imiss) then twac(k) = 0.1*float(itwac(k)) else twac(k) = rmiss endif enddo c c Fill in missing values of twac call patch0(twac,rmiss,mft) c c +++ Start Probability calculation c c Open and read the discriminant analysis coefficients for the c algorithm with lightning lucof = 26 if (ibasin .eq. 1) then fncoef = 'rii2o_al_coef.dat' fnprob = 'rii2o_al_prob.dat' elseif (ibasin .eq. 2) then fncoef = 'rii2o_ep_coef.dat' fnprob = 'rii2o_ep_prob.dat' elseif (ibasin .eq. 3) then fncoef = 'rii2o_wp_coef.dat' fnprob = 'rii2o_wp_prob.dat' elseif (ibasin .eq. 4) then fncoef = 'rii2o_io_coef.dat' fnprob = 'rii2o_io_prob.dat' elseif (ibasin .eq. 5) then fncoef = 'rii2o_sh_coef.dat' fnprob = 'rii2o_sh_prob.dat' else ierr = 1 endif c c get SHIPS_COEF env variable call getenv( "SHIPS_COEF",coef_location ) fncol=trim( coef_location )//fncoef c open(file=fncol,unit=lucof,form='formatted', + status='old',err=900) c c Read the header line read(lucof,*,err=900,end=900) it1,it2,idelt,nvar,ngrp nft = (it2-it1)/idelt c c Read the coefficients, predictor labels, sample means and standard deviations do k=0,nft do j=0,nvar read(lucof,*,err=900,end=900) idum1,(dfcoef(j,m,k),m=1,ngrp), + itdum,dlab(j),dmean(j,k), + dstd(j,k) enddo enddo c close(lucof) c c Open and read the file for converting discriminant function values c to probabilities c get SHIPS_COEF env variable call getenv( "SHIPS_COEF",coef_location ) fncol=trim( coef_location )//fnprob c open(file=fncol,unit=lucof,form='formatted', + status='old',err=900) c c Read the main header line read(lucof,*,err=900,end=900) it1,it2,idelt,nvar,ngrp nft = (it2-it1)/idelt c c Read the probabilities do k=0,nft c Read the header line for the current time read(lucof,*,err=900,end=900) iptime(k), + xbar(k),xstd(k), + nx(k),x1(k),x2(k),dx(k) c do m=1,ngrp read(lucof,400,err=900,end=900) glab 400 format(a7) do i=0,nx(k) read(lucof,410,err=900,end=900) dfprob(i,m,k) 410 format(f6.1) enddo enddo enddo c close(lucof) c c Start the main time loop c c Specify the time increment (hr) of the input data idelti = 6 itinc = idelt/idelti c do k=0,nft ktime1 = k*idelt ktime2 = ktime1 + idelt k1 = ktime1/idelti k2 = ktime2/idelti c nmissing = 0 c c Assemble input data c c 0 Constant term dfinp(0) = 1.0 c c 1 Potential intensity call rtavg(vmpi,mft,k1,k2,rmiss,vmpiavg) if (vmpiavg .ge. rmiss) then nmissing = nmissing + 1 else dfinp( 1) = vmpiavg - vm00 endif c c Check for minimum required potential intensity if (dfinp(1) .lt. potmin_final) then ipotflag = 1 else ipotflag = 0 endif c c 2 SHDC call rtavg(shdc,mft,k1,k2,rmiss,shdcavg) if (shdcavg .ge. rmiss) then nmissing = nmissing + 1 else dfinp( 2) = shdcavg endif c c 3 D200 call rtavg(d200,mft,k1,k2,rmiss,d200avg) if (d200avg .ge. rmiss) then nmissing = nmissing + 1 else dfinp( 3) = d200avg endif c c 4 Persistence dfinp( 4) = dvm12 c c 5 PC30 if (igoes(7) .eq. imiss) then nmissing = nmissing + 1 else dfinp( 5) = float(igoes(7)) endif c c 6 TBSTD if (igoes(4) .eq. imiss) then nmissing = nmissing + 1 else dfinp( 6) = 0.1*float(igoes(4)) endif c c 7 RHCN call rtavg(rhcn,mft,k1,k2,rmiss,rhcnavg) if (rhcnavg .ge. rmiss) then nmissing = nmissing + 1 else dfinp( 7) = rhcnavg endif c c 8 RHLO call rtavg(rhlo,mft,k1,k2,rmiss,rhloavg) if (rhloavg .ge. rmiss) then nmissing = nmissing + 1 else dfinp( 8) = rhloavg endif c c 9 TWAT if (twac(k2) .lt. rmiss .and. twac(k1) .lt. rmiss) then twat = twac(k2)-twac(k1) dfinp( 9) = twat else nmissing = nmissing + 1 endif c c 10,11 VMAX and VMAX**2 if (vm00 .lt. rmiss) then dfinp(10) = vm00 c divide vmax predcitor by vmax2_scale to match what is c used by rii2.f and mormalize Vmax**2 to have values c similar to Vmax dfinp(11) = (vm00**2)/vmax2_scale else nmissing = nmissing + 1 endif c if (nmissing .eq. 0) then c Calculate discriminant function values for each group do m=1,ngrp dfval(m) = dfcoef(0,m,k) do j=1,nvar dfval(m) = dfval(m) + dfinp(j)*dfcoef(j,m,k) enddo enddo else do m=1,ngrp dfval(m) = rmiss enddo endif c c Calculate x coordinate to convert df to probabilities x = dfval(2) - dfval(1) xnorm = (x - xbar(k))/xstd(k) c call dftop1(xnorm,gprob,dfprob(0,1,k),mxg, + mx,nx(k),dx(k),x1(k),x2(k)) c if (nmissing .ne. 0 .or. ierr .ne. 0) then do m=1,mxg gprob(m) = -99.9 enddo endif c c Calculate RI probability with climo values of input parameters do j=1,nvar dfinpc(j) = dmean(j,k) enddo c do m=1,ngrp dfvalc(m) = dfcoef(0,m,k) do j=1,nvar dfvalc(m) = dfvalc(m) + dfinpc(j)*dfcoef(j,m,k) enddo enddo c x = dfvalc(2) - dfvalc(1) xnorm = (x - xbar(k))/xstd(k) c call dftop1(xnorm,gprobc,dfprob(0,1,k),mxg, + mx,nx(k),dx(k),x1(k),x2(k)) enddo c if (gprob(mxg) .ge. 0.0 .and. gprob(mxg) .le. 100.0) then probric2 = gprob(mxg) else probric2 = 999.0 endif c c Set probri to zero if there is not enough potential intensity to c support RI if (ipotflag .eq. 1) probric2 = 0.0 c c +++ End Primary RI probability calculation c c +++ Start output file section c c Skip writing if luout < 0 if (luout .lt. 0) return c c Write RI probabilities c do m=1,mxg igprob(m) = nint(gprob(m)) enddo c if (ibasin .eq. 1) then cbasin = 'Atlantic Algorithm ' elseif (ibasin .eq. 2) then cbasin = 'East/Central Pacific Algorithm' elseif (ibasin .eq. 3) then cbasin = 'West Pacific Algorithm ' elseif (ibasin .eq. 4) then cbasin = 'N. Indiain Ocean Algorithm ' elseif (ibasin .eq. 5) then cbasin = 'S. Hemisphere Algorithm ' endif c write(luout,600) cbasin,sname,natcf8,ismon,isday,isyr,itime 600 format(//,10x,'** RI INDEX ',a30,' for cyclone ', + a10,' **', + /,10x,'** ATCFID=',a10,' Initialized ', + i2.2,'/',i2.2,'/',i4,1x,i2.2, + ' UTC **') c write(luout,601) igprob(mxg) 601 format(/,'Probability of Rapid Intensification= ',i5,'%', + ' (55 kt or more max wind increase in the next 48 hr)') c c ++ Probability contribution table c gcont(0) = gprobc(mxg) ptotal = gprob(mxg)-gprobc(mxg) do j=1,nvar gcont(j) = (dfcoef(j,2,0)-dfcoef(j,1,0))*(dfinp(j)-dmean(j,0)) enddo c gtotal = 0.0 do j=1,nvar gtotal = gtotal + gcont(j) enddo c do j=1,nvar gcont(j) = ptotal*gcont(j)/gtotal enddo c write(luout,610) 610 format(/,'Predictor Contributions to RI Probability',/, + 'Predictor Name Normalized Value Prob Contribtuion') c do j=0,nvar if (j .eq. 0) then pnorm = 0.0 else pnorm = (dfinp(j)-dmean(j,0))/dstd(j,0) endif c igcont = nint(gcont(j)) write(luout,615) plab(j),pnorm,igcont 615 format(a25,1x,f6.1,7x,i5,'%') enddo c write(luout,620) 620 format(/,'Note: Normalized Predictor Values are in std deviation', + ' units relative to sample mean values.',/, + ' (e.g., Value of 1.5 is 1.5 std dev above the mean)') c c Normal return close(lucof) return c c Error processing 900 continue ierr=1 write(6,*) ' Error in subroutine ric2' close(lucof) return c end subroutine rtavg(ft,mft,k1,k2,rmiss,ftavg) c This routine averaged ftime from k1 to k2, excluding missing values. c If all values of ftime are missing, ftavg is set to rmiss. c dimension ft(0:mft) c rcount = 0.0 ftavg = 0.0 do k=k1,k2 c write(6,*) 'k,ft(k)',ft(k) if (ft(k) .lt. rmiss) then ftavg = ftavg + ft(k) rcount = rcount + 1.0 endif enddo c if (rcount .ge. 1.0) then ftavg = ftavg/rcount else ftavg = rmiss endif c c write(6,*) 'rcount,ftavg',rcount,ftavg c return end