c ++++BEGIN ANNULAR HURRICANE INDEX MODULE++++
c
c version 19.1.0
c
c
c This set of routines calculates the probability of a tropical
c cyclone becoming an annular cyclone in the next 24h. They were
c adapted to be cut and pasted into the SHIPS model and run in
c realtime.
c-----7---------------------------------------------------------------72
SUBROUTINE id_annular_op(luout,stname,tcfid,fn_ships,fn_ird1,
. fn_ird2,fn_iri1,fn_iri2,ilmon,ilday,
. iyr,iltime,dfval,ann_prob,ahi_sst)
c-----7---------------------------------------------------------------72
c id_annular_op.f determines the probability of a tropical
c cyclone becoming an annular cyclone within 24 h. It uses the most
c recent SHIPS run and requires at least 1 IR data file within 12h
c of the SHIPS data file, but ideally uses 2 IR data files.
c
c
c Operational version of id_annular.test3.f written by Thomas Cram,
c CIRA/CSU, February 2006
c
c Adapted from McIDAS GOES-IR AREA file reader program written
c by Jim Kossin, CIRA/CSU, Jan. 2002
c
c Input:
c luout = unit number for writing data output
c stname = storm name, 10 character string
c tcfid = ATCF storm id, 8 character string
c fn_ships = SHIPS data file name - eg. 'lsdiag.dat'
c fn_ird1 = IR data file name - eg. 'IRRP1.dat'
c fn_ird2 = IR data file name - eg. 'IRRP2.dat'
c fn_iri1 = IR information file name - eg. 'IRRP1.inf'
c fn_iri2 = IR information file name - eg. 'IRRP2.inf'
c ahi_sst = real sst array(0:mft) with AHI sst for all forecast times
c the SST is the same as used by SHIPS model, but J Cione
c cooling is NOT applied. The selections of SST from climo and
c observed,and different types of SST is done in iships.f
c
c Output:
c ilmon = storm month from SHIPS file, eg. '09' for Sep
c ilday = storm day from SHIPS file, eg. '23' for the 23rd
c iyr = storm year from SHIPS file, eg. '2005'
c iltime = storm time from SHIPS file, UTC, eg '0600'
c dfval = discriminant function value normalized to
c vary between 0 and 1 = annular hurr index
c ann_prob = probability of annular structure
c
c
c Should be machine independent (portable) assuming that:
c 1) The RECL specifiers in the OPEN statements refer to bytes,
c not words.
c 2) The default size of integers is 4 bytes.
c
c REQUIRES SUBROUTINES:
c read_ships
c ir_prof_avg2
c detend - removed (KM)
c read_radprof
c moment
c chartointmon
c ahindex
c tdiff
c
c Modified by Andrea Schumacher for operational use, CIRA, July 2006
c
c MODIFIED: February 16, 2007 - JAK
c 1) changed the number of disciminants used to 5
c 2) more changes to ahindex and block data.
c
c Modified by Stephanie Stevenson, April 2020
c Changed mft and lsdiag reading for 7day input
c
c-----7---------------------------------------------------------------72
parameter(maxavg=1000) ! max size of 1D radial profile arrays
parameter(radmax=600.) ! max radius of profile output files
parameter(maxrad=radmax/4) ! max number of elements in profiles
dimension raz(maxrad)
character *10 stname
character*9 stnm ! storm name
character*9 stnmabr ! abbreviated storm name
character*8 tcfid
character*2 bas ! ocean basin (AT,EP,WP,SP,CP,NI,SI)
real xmsng
character*3 irmonchar
character*40 errmsg0,errmsg1,errmsg2,errmsg3,errmsg4,errmsg5
c Averaged IR data
real prof_avg(maxrad),std_avg(maxrad)
integer irerr,irerr1,irerr2,irerrflg
integer ibuff
integer useir1,useir2
c Annular hurricane predictors
real std_cpix,eye_diff
real cpix,std,var,weye,shrd,shrs,shrg,u200,t200,vmax
! Other screening variables
real gen_sst
real mmaxval
c SHIPS data
parameter(mft=28)
integer ilyr,ilmon,ilday,iltime,ivmx,iyr
integer ishrd(0:mft),ishtd(0:mft),ishrs(0:mft),ishts(0:mft),
+ ishrg(0:mft),iu200(0:mft),it200(0:mft),
+ irefc(0:mft),ivmax(0:mft)
c this SST is passed from iships. It's the same SST as used by SHIPS
c model, with the selection feom climo, observed, and various available
c SST datasets. However, J Cione cooling is NOT used for AHI
real ahi_sst(0:mft)
integer endoffile
integer ludiag
character*4 sname4
character*34 fn_ships
character*34 fn_ird1
character*34 fn_ird2
character*34 fn_iri1
character*34 fn_iri2
integer iphold,irhold
real tbmin,tbmax
integer smiss ! SHIPS missing data value (=9999)
c Array counters
integer sample_count
integer luout
real dv,dfval,ann_prob
integer screen_fail
logical radexist
c-----7---------------------------------------------------------------72
! Standardized data
real weye_s, var_s, std_s, sst_s, u200_s
c-----7---------------------------------------------------------------72
! Statistical moments and discriminant weights (block data)
REAL sdev(5),means(5),dweights(5),div
common /MOMENTS/ sdev,means
common /DWEIGHTS/ dweights
common /DIVIDER/ div
c-----7---------------------------------------------------------------72
c actual lag in days of the observed SST
c irlag - WSST
c irdllag - daily
c iohclag - daily observed OHC
integer irlag,irdllag,iohclag,iobslag
c
c Set ishort=1 for short form of output for SHIPS text file
ishort=1
c Specify total number of screening variables for output
nscreen=7
c
ludiag=50
xmsng=-1.0
endoffile=0
sample_count=0
ibuff=0
! SHIPS missing data value
smiss=9999
! Discriminant value missing data
dmiss=9999.
ir1ex=0
ir2ex=0
irex=0
ir1tg=0
ir2tg=0
irtg=0
ivmx=0
dv=0.0
dfval=0.0
ann_prob=0.0
c write(errmsg0,631) 'NOTE:',', ANNULAR INDEX RAN NORMALLY'
errmsg0=' ANNULAR INDEX RAN NORMALLY'
write(errmsg1,632) 'ERR=1',', SHIPS FILE MISSING'
write(errmsg2,633) 'ERR=2',', BOTH IR FILES BAD OR MISSING'
write(errmsg3,634) 'ERR=3',', IR & SHIPS DATA > 12h APART'
write(errmsg4,635) 'ERR=4',', SHIPS DATA MISSING'
errmsg5='NOTE: 1 INSTEAD OF 2 GOES FILES USED'
c write(errmsg5,636)
c636 format(5x,'NOTE: 1 INSTEAD OF 2 GOES FILES USED ')
c637 format(/,/,/,/,a40)
637 format(' ##',4x,a40)
638 format(' ##',4x,a40)
639 format(' ##',4x,a40)
c639 format(/,13x,a40)
631 format(5x,a5,a28,2x)
632 format(5x,a5,a20,10x)
633 format(5x,a5,a30)
634 format(5x,a5,a29,1x)
635 format(5x,a5,a20,10x)
mmaxval=-3.105263
c-----7---------------------------------------------------------------72
c open data output file
c
c open(unit=luout,file='./annout.dat',status='replace')
c
c
c open warning message file
c
c open(unit=luwarn,file='./warning_msg2.out',status='unknown')
c
c
c determine endian-ness of host machine
c
c call detend(indwrd)
c-----7---------------------------------------------------------------72
c Verify that SHIPS data file exists. If not, fail. If so, read header.
c-----7---------------------------------------------------------------72
! Open SHIPS file on first pass
open(unit=ludiag,file=fn_ships,status='old',err=901)
!Read the date and time info from SHIPS data file
c read(ludiag,110,err=901,end=901) sname4,ilyr,ilmon,ilday,
c + iltime,ivmx,iheader
c 110 format(1x,a4,1x,3(i2.2),1x,i2.2,1x,i4,a110)
read(ludiag,110,err=901,end=901) sname4,ilyr,ilmon,ilday,
+ iltime,ivmx
110 format(1x,a4,1x,3(i2.2),1x,i2.2,1x,i4)
close(ludiag)
c-----7---------------------------------------------------------------72
c
c Write header to output file, first calculating 4-digit year:
c
c-----7---------------------------------------------------------------72
if (ilyr.lt.40) then
iyr=2000+ilyr
write(stnm,25) sname4(1:3)
write(stnmabr,26) ilyr,bas,sname4(1:3)
else
iyr=1900+ilyr
write(stnm,25) sname4(1:3)
write(stnmabr,28) ilyr,bas,sname4(1:3)
endif
25 format(a3)
26 format('20',i2.2,a2,a3)
28 format('19',i2.2,a2,a3)
c
write(luout,600) tcfid,stname,ilmon,ilday,
+ ilyr,iltime
600 format(/,' ##',9x,'ANNULAR HURRICANE INDEX (AHI) ',a8,1x,a10,
+ 1x,i2.2,'/',i2.2,'/',i2.2,2x,i2.2,' UTC ##')
c-----7---------------------------------------------------------------72
c Verify 1st .inf and .dat files both exist. If so, read header and
c determine if .dat file time is > 12h from SHIPS data file time.
c-----7---------------------------------------------------------------72
! Does IR radial profile file exist?
inquire(file=fn_ird1,exist=radexist)
ir1ex=0
if(radexist) then
open(unit=ludiag+2,file=fn_iri1,status='old',err=902)
read(unit=ludiag+2,fmt=111,err=902,end=902) irday,
+ irmonchar,iryr,irtime
111 format(/,/,/,/,1x,i2.2,2x,a3,4x,i2,1x,i2)
ir1ex=1
call chartointmon(irmonchar,irmon)
call tdiff(iryr,irmon,irday,irtime,ilyr,ilmon,ilday,iltime,
+ idelt)
if (idelt.gt.12.or.idelt.lt.-12) then
ir1tg=0
else
ir1tg=1
endif
902 endif
close(ludiag+2)
c-----7---------------------------------------------------------------72
c Verify 2nd .inf and .dat files both exist. If so, read header and
c determine if .dat file time is > 12h from SHIPS data file time.
c-----7---------------------------------------------------------------72
! Does IR radial profile file exist?
inquire(file=fn_ird2,exist=radexist)
ir2ex=0
if(radexist) then
open(unit=ludiag+2,file=fn_iri2,status='old',err=903)
read(unit=ludiag+2,fmt=112,err=903,end=903) irday,
+ irmonchar,iryr,irtime
112 format(/,/,/,/,1x,i2.2,2x,a3,4x,i2,1x,i2)
ir2ex=1
call chartointmon(irmonchar,irmon)
call tdiff(iryr,irmon,irday,irtime,ilyr,ilmon,ilday,iltime,
+ idelt)
if (idelt.gt.12.or.idelt.lt.-12) then
ir2tg=0
else
ir2tg=1
endif
903 endif
close(ludiag+2)
irex = ir1ex+ir2ex
irtg = ir1tg+ir2tg
! If both IR files are missing, display error in data output
! file and skip to the end (ie, calculation fails)
if (irex.eq.0) then
write(luout,638) errmsg2
return
endif
! If both IR files are more than 12h before or after the
! current SHIPS data, display error in data output file
! and skip to the end (ie, calcualtion fails)
if (irtg.eq.0) then
write(luout,638) errmsg3
return
endif
if ((ir1ex.eq.1).and.(ir1tg.eq.1)) then
useir1 = 1
else
useir1 = 0
endif
if ((ir2ex.eq.1).and.(ir2tg.eq.1)) then
useir2 = 1
else
useir2 = 0
endif
c-----7---------------------------------------------------------------72
c Open and Read in SHIPS data
c-----7---------------------------------------------------------------72
! Open SHIPS file
open(unit=ludiag,file=fn_ships,status='old',err=901)
5 call read_ships(ludiag,sname4,ilyr,ilmon,ilday,iltime,ivmx,
+ iu200,it200,
+ irefc,ishrd,ishtd,ishrs,
+ ishts,ishrg,endoffile)
c-----7---------------------------------------------------------------72
c Call routine 'ir_prof_avg2' to read in IR Tb radial profiles,
c compute the Tb standard deviations, and then compute average radial
c profiles of each.
c-----7---------------------------------------------------------------72
call ir_prof_avg2(fn_ird1,fn_ird2,fn_iri1,fn_iri2,
+ iyr,useir1,useir2,prof_avg,std_avg,
+ raz,ibuff,irerr)
if (sum(prof_avg).le.0) then
dv=dmiss
dfval=dmiss
ann_prob=0.0
write(luout,638) errmsg2
return
endif
c-----7---------------------------------------------------------------72
c Screen SHIPS data -- skip if environmental conditions would not
c support an annular hurricane
c-----7---------------------------------------------------------------72
screen_fail = 0
ahi = 0.0
c
vmxr = float(ivmx)
if (ivmx .lt. 85) then
dv=dmiss
dfval=dmiss
ann_prob=0.0
screen_fail=screen_fail+1
if (ishort .ne. 1) then
write(luout,652) vmxr
652 format(5x,'SCREENING STORM INTENSITY = ',f6.1,
+ ' kt > 84 kt? ---> FAILED')
endif
else
if (ishort .ne. 1) then
write(luout,653) vmxr
653 format(5x,'SCREENING STORM INTENSITY = ',f6.1,
+ ' kt > 84 kt? ---> PASSED')
endif
endif
rsstr = ahi_sst(0)
if (rsstr .lt. 24.3) then
dv=dmiss
dfval=dmiss
ann_prob=0.0
screen_fail=screen_fail+1
if (ishort .ne. 1) then
write(luout,654) rsstr
654 format(5x,'SCREENING SST = ',f6.1,
+ ' C > 24.3 C ? ---> FAILED')
endif
else
if (ishort .ne. 1) then
write(luout,655) rsstr
655 format(5x,'SCREENING SST = ',f6.1,
+ ' C > 24.3 C ? ---> PASSED')
endif
endif
if (rsstr .gt. 29.1) then
dv=dmiss
dfval=dmiss
ann_prob=0.0
screen_fail=screen_fail+1
if (ishort .ne. 1) then
write(luout,656) rsstr
656 format(5x,'SCREENING SST = ',f6.1,
+ ' C < 29.1 C ? ---> FAILED')
endif
else
if (ishort .ne. 1) then
write(luout,657) rsstr
657 format(5x,'SCREENING SST = ',f6.1,
+ ' C < 29.1 C ? ---> PASSED')
endif
endif
c
shrdr = 0.1*float(ishrd(0))
if (shrdr .gt. 22.0) then
dv=dmiss
dfval=dmiss
ann_prob=0.0
screen_fail=screen_fail+1
if (ishort .ne. 1) then
write(luout,658) shrdr
658 format(5x,'SCREENING VERT. SHR = ',f6.1,
+ ' kt < 22.0 kt? ---> FAILED')
endif
else
if (ishort .ne. 1) then
write(luout,659) shrdr
659 format(5x,'SCREENING VERT. SHR = ',f6.1,
+ ' kt < 22.0 kt? ---> PASSED')
endif
endif
u200r = 0.1*float(iu200(0))
if (u200r .lt. -23.0) then
dv=dmiss
dfval=dmiss
ann_prob=0.0
screen_fail=screen_fail+1
if (ishort .ne. 1) then
write(luout,662) u200r
662 format(5x,'SCREENING 200 hPa ZONAL WIND = ',f6.1,
+ ' kt > -23.0 kt? ---> FAILED')
endif
else
if (ishort .ne. 1) then
write(luout,663) u200r
663 format(5x,'SCREENING 200 hPa ZONAL WIND = ',f6.1,
+ ' kt > -23.0 kt? ---> PASSED')
endif
endif
if (u200r .gt. 3.0) then
dv=dmiss
dfval=dmiss
ann_prob=0.0
screen_fail=screen_fail+1
if (ishort .ne. 1) then
write(luout,664) u200r
664 format(5x,'SCREENING 200 hPa ZONAL WIND = ',f6.1,
+ ' kt < 3.0 kt? ---> FAILED')
endif
else
if (ishort .ne. 1) then
write(luout,665) u200r
665 format(5x,'SCREENING 200 hPa ZONAL WIND = ',f6.1,
+ ' kt < 3.0 kt? ---> PASSED')
endif
endif
refcr = 0.1*float(irefc(0))
if (refcr .lt. -9.0) then
dv=dmiss
dfval=dmiss
ann_prob=0.0
screen_fail=screen_fail+1
if (ishort .ne. 1) then
write(luout,666) refcr
666 format(5x,'SCREENING 200 hPa MOM FLUX CONV = ',f6.1,
+ ' m/s-day > -9.0 m/s-day? ---> FAILED')
endif
else
if (ishort .ne. 1) then
write(luout,667) refcr
667 format(5x,'SCREENING 200 hPa MOM FLUX CONV = ',f6.1,
+ ' m/s-day > -9.0 m/s-day? ---> PASSED')
endif
endif
if (refcr .gt. 11.0) then
dv=dmiss
dfval=dmiss
ann_prob=0.0
screen_fail=screen_fail+1
if (ishort .ne. 1) then
write(luout,668) refcr
668 format(5x,'SCREENING 200 hPa MOM FLUX CONV = ',f6.1,
+ ' m/s-day < 11.0 m/s-day? ---> FAILED')
endif
else
if (ishort .ne. 1) then
write(luout,669) refcr
669 format(5x,'SCREENING 200 hPa MOM FLUX CONV = ',f6.1,
+ ' m/s-day < 11.0 m/s-day? ---> PASSED')
endif
endif
702 format(5x,a27,i6,a13,13x,a11)
703 format(5x,a15,f6.1,a14,24x,a11)
704 format(5x,a21,f6.1,a16,16x,a11)
705 format(5x,a25,f6.1,a17,11x,a11)
706 format(5x,a26,f6.1,a25,2x,a11)
707 format(5x,a29,f6.1,a16,8x,a11)
712 format(19x,a35)
713 format(19x,a28)
714 format(19x,a44)
715 format(19x,a36)
716 format(19x,a41)
c-----7---------------------------------------------------------------72
c Skip if missing data
missdata_fail=0
if (ahi_sst(0).ge.dmiss-1) then
dv=dmiss
dfval=dmiss
ann_prob=0.0
missdata_fail=missdata_fail+1
endif
if (ishrd(0).eq.smiss) then
dv=dmiss
dfval=dmiss
ann_prob=0.0
missdata_fail=missdata_fail+1
endif
if (ishrs(0).eq.smiss) then
dv=dmiss
dfval=dmiss
missdata_fail=missdata_fail+1
endif
if (ishrg(0).eq.smiss) then
dv=dmiss
dfval=dmiss
ann_prob=0.0
missdata_fail=missdata_fail+1
endif
if (iu200(0).eq.smiss) then
dv=dmiss
dfval=dmiss
ann_prob=0.0
missdata_fail=missdata_fail+1
endif
if (it200(0).eq.smiss) then
dv=dmiss
dfval=dmiss
ann_prob=0.0
missdata_fail=missdata_fail+1
endif
if (irefc(0).eq.smiss) then
dv=dmiss
dfval=dmiss
ann_prob=0.0
missdata_fail=missdata_fail+1
endif
if (ivmx.eq.smiss) then
dv=dmiss
dfval=dmiss
ann_prob=0.0
missdata_fail=missdata_fail+1
endif
if (missdata_fail.gt.0) then
write(luout,638) errmsg4
return
endif
c-----7---------------------------------------------------------------72
c
c Find radius of coldest Tb pixel in averaged IR profile
c and grab standard deviation at that radius
c
tbmin=1.e3
tbmax=-1.
irhold=0
iphold=0
! Determine size of prof_avg (i.e. where prof_avg does not
! contain missing data)
do iprof=1,maxrad
if (prof_avg(iprof).ne.xmsng) then
iphold=iprof
endif
enddo
! Find radius of coldest pixel and std dev.
do iprof=1,maxrad
if ((prof_avg(iprof).lt.tbmin).and.
+ (prof_avg(iprof).ne.xmsng)) then
tbmin=prof_avg(iprof)
cpix=raz(iprof)
std_cpix=std_avg(iprof)
irhold=iprof
endif
enddo
c-----7---------------------------------------------------------------72
c Skip if missing data
if ((iphold.lt.2).or.(irhold.eq.0)) then
dv=dmiss
dfval=dmiss
ann_prob=0.0
return
endif
if (std_cpix.lt.0) then
dv=dmiss
dfval=dmiss
ann_prob=0.0
return
endif
c-----7---------------------------------------------------------------72
c Calculate variance of radial profile
c
call moment(prof_avg(1:iphold),iphold,ave_prof,adev_prof,
+ sdev_prof,var_prof,skew_prof,curt_prof)
c
c Calculate warm eye - cold ring difference
c
do iprof=1,irhold
if (prof_avg(iprof).gt.tbmax) then
tbmax=prof_avg(iprof)
eye_diff=tbmax-tbmin
endif
enddo
c-----7---------------------------------------------------------------72
c Screen IR data -- skip if cold eye exists
if (cpix .lt. 50.) then
dv=dmiss
dfval=dmiss
ann_prob=0.0
if (ishort .ne. 1) then
write(luout,670) cpix
670 format(5x,'SCREENING GOES RAD COLD BR TEMP = ',f6.1,
+ ' km > 50.0 km? ---> FAILED')
endif
screen_fail=screen_fail+1
else
if (ishort .ne. 1) then
write(luout,671) cpix
671 format(5x,'SCREENING GOES RAD COLD BR TEMP = ',f6.1,
+ ' km > 50.0 km? ---> PASSED')
endif
endif
if (eye_diff .lt. 15.) then
dv=dmiss
dfval=dmiss
ann_prob=0.0
if (ishort .ne. 1) then
write(luout,672) eye_diff
672 format(5x,'SCREENING GOES EYE-RING BR TEMP = ',f6.1,
+ ' C > 15.0 C ? ---> FAILED')
endif
screen_fail=screen_fail+1
else
if (ishort .ne. 1) then
write(luout,673) eye_diff
673 format(5x,'SCREENING GOES EYE-RING BR TEMP = ',f6.1,
+ ' C > 15.0 C ? ---> PASSED')
endif
endif
if (screen_fail.gt.0) then
if (ishort .eq. 1) then
nfail = screen_fail
npass = nscreen-nfail
write(luout,780) npass,nfail
780 format(' ## STORM NOT ANNULAR, SCREENING STEP FAILED,',
+ ' NPASS=',i1,' NFAIL=',i1,15x,
+ ' ##')
iahi = ifix(ahi+0.49)
write(luout,782) iahi
else
write(luout,723)
723 format(25x,'STORM NOT ANNULAR, FAILED SCREENING',/)
c
write(luout,722) '******************************',
+ '***********************'
write(luout,721) ahi
write(luout,722) '******************************',
+ '***********************'
endif
return
else
if (ishort .eq. 1) then
write(luout,781)
781 format(' ## PASSED SCREENING STEP, MIGHT BE ANNULAR,',
+ ' CALCULATE AHI FROM DISCRIMINANT ANALYSIS ##')
else
write(luout,724)
724 format(25x,'STORM MAY BE ANNULAR, PASSED SCREENING',/,
+ 23x,'CALCULATE AHI FROM DISCRIMINANT ANALYSIS',/)
endif
endif
701 format(25x,a31,/)
782 format(' ## AHI=',i3,' (AHI OF 100 IS BEST FIT TO ANN.',
+ ' STRUC., 1 IS MARGINAL, 0 IS NOT ANNULAR) ##')
721 format(15x,'* ANNULAR HURRICANE INDEX (AHI) VALUE = ',f5.0,
+ ' *',/,15x,
+ '* (AHI=100. IS BEST MATCH TO ANNULAR STRUCTURE) *',
+ /,15x,'* (AHI= 1. IS WORST MATCH TO ANNULAR STRUCTURE) *',
+ /,15x,'* (AHI= 0. FOR NO ANNULAR STRUCTURE) *')
722 format(15x,a30,a23)
708 format(5x,a33,f5.1,a13,8x,a11)
709 format(5x,a31,f5.1,a10,13x,a11,/)
c-----7---------------------------------------------------------------72
c Standardize data
gen_sst = ahi_sst(0)
u200 = real(iu200(0)) / 10.0
weye_s = (eye_diff - means(1)) / sdev(1)
var_s = (var_prof - means(2)) / sdev(2)
std_s = (std_cpix - means(3)) / sdev(3)
sst_s = (gen_sst - means(4)) / sdev(4)
u200_s = (u200 - means(5)) / sdev(5)
c Get annular hurricane index
call ahindex(weye_s,var_s,std_s,sst_s,u200_s ,dv,dfval,ann_prob)
c Calculate annular hurricane index (ahi) from the discriminant value
c
c Specific df scaling factors
c ahi = ahmin at dfmin
c ahi = ahmax at dfmax
c
c dfmin, dfmax based upon 1995-2006 developmental data and should
c result in a ~ 96% hit rate and ~4% false alarm rate - JAK
c
dfmin = -0.3
dfmax = 2.3
ahmin = 1.0
ahmax = 100.0
c
ahslope = (ahmax-ahmin)/(dfmax-dfmin)
ahyint = ahmin - ahslope*dfmin
c
ahi = dfval*ahslope + ahyint
c
if (ahi .lt. 1.0) ahi = 1.0
if (ahi .gt. 100.0) ahi = 100.0
sample_count=sample_count+1
c-----7---------------------------------------------------------------72
c Write data
if (ishort .eq. 1) then
iahi = ifix(ahi+0.49)
write(luout,782) iahi
else
write(luout,722) '******************************',
. '***********************'
write(luout,721) ahi
write(luout,722) '******************************',
. '***********************'
endif
c
c iwcon = 0
c if (iwcon .eq. 1) then
c write individual contributions to annular index
c write(luout,603) 'CONTRIBUTIONS TO ANNULAR INDEX'
c write(luout,604) 'Vertical Shear',' ',dweights(5)*shrd_s +
c + dweights(6)*shrg_s
c write(luout,605) 'Radius of Coldest Pixel',' ',
c + dweights(1)*cpix_s
c write(luout,606) 'GOES Core Symmetry Factor',' ',
c + dweights(3)*var_s + dweights(4)*std_s
c write(luout,607) 'GOES Eye Temperature',' ',dweights(2)*weye_s
c write(luout,608) '200 hPa Eddy Fluxes',' ',dweights(8)*refc_s
c write(luout,609) '200 hPa Zonal Wind',' ',dweights(7)*u200_s
c write(luout,610) 'Storm Intensity',' ',dweights(9)*vmax_s
c write(luout,621) '_','_____________________________',
c + '__________'
c write(luout,611) 'TOTAL',' ', dv
c write(luout,622) 'TOTAL > ',div, ' ---->',
c + ' ANNULAR STRUCTURE POSSIBLE'
c write(luout,623) 'TOTAL < ',div, ' ---->',
c + ' ANNULAR STRUCTURE VERY UNLIKELY'
c endif
602 format(a15,5x,a36,f4.1,a1,4x,a5)
603 format(/,a56/)
604 format(23x,a14,a16,f7.3)
605 format(23x,a23,a7,f7.3)
606 format(23x,a25,a5,f7.3)
607 format(23x,a20,a10,f7.3)
608 format(23x,a19,a11,f7.3)
609 format(23x,a18,a12,f7.3)
610 format(23x,a15,a15,f7.3)
611 format(23x,a5,a25,f7.3,/)
620 format(a15,a30,a25)
621 format(a24,a29,a7)
622 format(18x,a8,f3.1,a6,a27)
623 format(18x,a8,f3.1,a6,a33)
! If only 1 of the IR files is within 12h of SHIPS data,
! display a note in the data output file, proceed.
c if ((irtg.eq.1).or.(irerr.gt.0)) then
c write(luout,639) errmsg5
c else
c write(luout,639) errmsg0
c endif
return
901 write(luout,637) errmsg1
RETURN
END
c-----7---------------------------------------------------------------72
c-------------------------- SUBROUTINES ------------------------------
c-----7---------------------------------------------------------------72
SUBROUTINE ir_prof_avg2(fn_ird1,fn_ird2,fn_iri1,fn_iri2,
+ iyr,useir1,useir2,prof_avg,std_avg,
+ raz,ibuff,irerr)
c-----7---------------------------------------------------------------72
c Routine to read in 2 IR radial profile data and average
c
c
c INPUT
c fn_ird1 = filename of 1st IR data file
c fn_ird2 = filename of 2nd IR data file
c fn_iri1 = filename of 1st IR information file
c fn_iri2 = filename of 2nd IR information file
c iyr = 4-digit storm year
c useir1 = integer indicating validity of IRRP1.dat, 0=n, 1=y
c useir2 = integer indicating validity of IRRP2.dat, 0=n, 1=y
c
c OUTPUT
c prof_avg = time average of IR profile
c std_avg = time average of IR standard deviations
c raz = radii (km) of prof_avg and std_avg
c ibuff = (integer) If equal to zero, then disregard prof_avg
c and std_avg in main program.
c
c-----7---------------------------------------------------------------72
parameter(maxbyt=1000000) ! max file size (bytes; multiple of 4)
byte tot(maxbyt) ! total data file (bytes)
parameter(maxele=3000,maxlin=3000) ! max image size (elems,lines)
dimension ximage(maxele,maxlin)
parameter(maxavg=1000) ! max size of 1D radial profile arrays
parameter(radmax=600.) ! max radius of profile output files (km)
parameter(maxrad=radmax/4) ! max number of elements in profiles
dimension azav(maxrad),raz(maxrad),std(maxrad) ! averaging arrays
dimension specr(maxrad),speci(maxrad) ! spectral coefficient arrays
parameter(maxtim=1000) ! max number of data files
dimension timarray(maxtim)
parameter(maxwrd=maxbyt/4) ! max file size (words)
integer navi(maxwrd) ! navigation block (words)
integer area(64) ! area block (words)
equivalence (tot(1),navi(1)) ! fills navigation from total
character*80 fn ! area file name
character*80 fnout ! generic output file name
character*80 pthar ! path to area file
logical swap ! little or big endian query
logical fnxst ! file name query
integer LL(2) ! (lat,lon) <---> (y,x) flag
character*12 stnm ! storm name
character*9 stnmabr ! abbreviated storm name
character*2 bas ! ocean basin (AT,EP,WP,SP,CP,NI,SI)
character*11 lmstarr(maxtim) ! LMST in "MM/DD/hh/mm" format
parameter(pi=3.1415926535897932384626433832795) ! define pi
integer tdiffmin,tdiff,irerr,irerr1,irerr2
character*29 fnrad ! Radial Tb profile data file
integer lurad
logical radexist
real xmsng
integer ir1ex,ir2ex,irex
integer ir1tg,ir2tg,irtg
integer useir1,useir2
character*34 fn_ird1
character*34 fn_ird2
character*34 fn_iri1
character*34 fn_iri2
character*80 pthbt ! path to best track file
parameter(maxfix=200) ! max no of best track fixes (per storm)
dimension vnz(maxfix),pnz(maxfix)
dimension xlanz(maxfix),xlonz(maxfix)
dimension tnz(maxfix),sdla(maxfix),sdlo(maxfix)
integer im(maxfix/4),id(maxfix/4)
c Start and end time of IR time averaging period
real bt_start,bt_end
c percent pixels along circles
dimension pp00(maxrad),pp10(maxrad),pp20(maxrad),pp30(maxrad)
dimension pp40(maxrad),pp50(maxrad),pp60(maxrad),pp70(maxrad)
parameter(maxbuff=40)
real buff_azav(maxrad,maxbuff),buff_std(maxrad,maxbuff),
+ prof_avg(maxrad),std_avg(maxrad)
integer ibuff,ip,ib
integer iphold,irhold,iphold_min
c-----7---------------------------------------------------------------72
xmsng=-1.0
itim=0
ibuff=0
iphold_min=1000
irnum=0
lurad=46
c Zero out buffer and averaged IR arrays
do i=1,maxrad
do j=1,maxbuff
buff_azav(i,j)=0.0
buff_std(i,j)=0.0
enddo
prof_avg(i)=0.0
std_avg(i)=0.0
enddo
c-----7---------------------------------------------------------------72
c Read in 1st radial profile
c-----7---------------------------------------------------------------72
if (useir1.ne.0) then
fnrad=fn_ird1
irnum=irnum+1
open(unit=lurad,file=fnrad,status='old',err=446)
call read_radprof(lurad,raz,azav,std,pp50,irerr1)
close(lurad)
c-----7---------------------------------------------------------------72
c Find radius of coldest Tb pixel and grab standard deviation
c at that radius
c
tbmin=1.e3
irhold=0
iphold=0
do iprof=1,maxrad
if ((azav(iprof).ne.xmsng).and.(std(iprof).ne.xmsng)) then
iphold=iprof
endif
enddo
do iprof=1,maxrad
if ((azav(iprof).lt.tbmin).and.(azav(iprof).ne.xmsng)) then
tbmin=azav(iprof)
cpix=raz(iprof)
std_cpix=std(iprof)
irhold=iprof
endif
enddo
! Check data -- skip if missing data
if ((tbmin.eq.xmsng).or.(std_cpix.eq.xmsng)) then
tbmin=1.e3
else if ((irhold.ne.0).and.(iphold.ne.0)) then
if (iphold.lt.iphold_min) iphold_min=iphold
c-----7---------------------------------------------------------------72
c Store current IR data in buffer arrays for averaging
c
do ip=1,maxrad
buff_azav(ip,irnum)=azav(ip)
buff_std(ip,irnum)=std(ip)
enddo
endif
endif
c-----7---------------------------------------------------------------72
446 continue
c-----7---------------------------------------------------------------72
c-----7---------------------------------------------------------------72
c Read in 2nd radial profile
c-----7---------------------------------------------------------------72
if (useir2.ne.0) then
fnrad=fn_ird2
irnum=irnum+1
open(unit=lurad,file=fnrad,status='old',err=449)
call read_radprof(lurad,raz,azav,std,pp50,irerr2)
close(lurad)
c-----7---------------------------------------------------------------72
c Find radius of coldest Tb pixel and grab standard deviation
c at that radius
c
tbmin=1.e3
irhold=0
iphold=0
do iprof=1,maxrad
if ((azav(iprof).ne.xmsng).and.(std(iprof).ne.xmsng)) then
iphold=iprof
endif
enddo
do iprof=1,maxrad
if ((azav(iprof).lt.tbmin).and.(azav(iprof).ne.xmsng)) then
tbmin=azav(iprof)
cpix=raz(iprof)
std_cpix=std(iprof)
irhold=iprof
endif
enddo
! Check data -- skip if missing data
if ((tbmin.eq.xmsng).or.(std_cpix.eq.xmsng)) then
tbmin=1.e3
else if ((irhold.ne.0).and.(iphold.ne.0)) then
if (iphold.lt.iphold_min) iphold_min=iphold
c-----7---------------------------------------------------------------72
c Store current IR data in buffer arrays for averaging
c
do ip=1,maxrad
buff_azav(ip,irnum)=azav(ip)
buff_std(ip,irnum)=std(ip)
enddo
endif
endif
c-----7---------------------------------------------------------------72
449 continue
c-----7---------------------------------------------------------------72
if (irnum.ne.0) then
! Average IR data
do ip=1,iphold_min
prof_avg(ip)=sum(buff_azav(ip,1:irnum))/2
std_avg(ip)=sum(buff_std(ip,1:irnum))/2
enddo
! Fill out rest with missing value, if necessary
if (iphold_min.lt.maxrad) then
do ip=iphold_min+1,maxrad
prof_avg(ip)=xmsng
std_avg(ip)=xmsng
enddo
endif
! Zero out buffer arrays
do ip=1,maxrad
do ib=1,maxbuff
buff_azav(ip,ib)=0.0
buff_std(ip,ib)=0.0
enddo
enddo
endif
irerr = irerr1 + irerr2
ccccc output prof_avg and raz to file
RETURN
END
c-----7---------------------------------------------------------------72
c-----7---------------------------------------------------------------72
subroutine read_ships(luls,sname4,ilyr,ilmon,ilday,iltime,ivmx,
+ iu200,it200,
+ irefc,ishrd,ishtd,ishrs,
+ ishts,ishrg,endoffile)
c Thomas Cram, CIRA, November 2005
c
c Reads information from the SHIPS model reanalysis data files
c (lsdiaga_1982_2004_rean.dat, lsdiage_1982_2004_rean.dat)
c
c INPUT
c luls = SHIPS data file unit number
c
c OUTPUT
c sname4 = First 4 letters of storm name
c ilyr = 2-digit year
c ilmon = UTC month
c ilday = UTC day
c iltime = UTC hour
c ivmx = Maximum winds (kt)
c iu200 = 200 mb zonal wind (kt)
c it200 = 200 mb temperature (deg C)
c irefc = Relative eddy momentum flux convergence (m/s/day)
c ishrd = 850-200 mb shear magnitude (kt)
c ishtd = Heading (deg) of SHRD shear vector
c ishrs = 850-500 mb shear magnitude (kt)
c ishts = Heading (deg) of SHRS shear vector
c ishrg = Generalized 850-200 mb shear magnitude (kt; takes into
c account all levels)
c endoffile = Integer specifying if the end of the file has been
c reached (0 = no, 1 = yes)
c
c-----7---------------------------------------------------------------72
c
parameter (mft=28)
character*180 iline
character*110 iheader
integer ilyr,ilmon,ilday,iltime,ivmx
integer iper,itime(0:mft),iincv(0:mft),ilat(0:mft),
+ ilon(0:mft),id20c(0:mft),
+ id26c(0:mft),ihcon(0:mft)
integer idist(0:mft),iu200(0:mft),it200(0:mft),
+ ie000(0:mft),iepos(0:mft),ieneg(0:mft),iepss(0:mft),
+ ienss(0:mft),irhlo(0:mft),irhmd(0:mft),irhhi(0:mft),
+ ishrd(0:mft),ishtd(0:mft),ishrs(0:mft),ishts(0:mft),
+ ishrg(0:mft),ipslv(0:mft),iz850(0:mft),id200(0:mft),
+ irefc(0:mft),it000(0:mft),ir000(0:mft),iz000(0:mft),
+ igoes(0:mft),igoesm3(0:mft),ird20(0:mft),ird26(0:mft),
+ irhcn(0:mft)
integer endoffile
character*4 sname4
endoffile=0
c Read the SHIPS data file
read(luls,110,end=1201,err=1201) sname4,ilyr,ilmon,ilday,
+ iltime,ivmx,iheader
110 format(1x,a4,1x,3(i2.2),1x,i2.2,1x,i4,a110)
1100 continue
c Format read definitions
111 format(a180)
112 format(1x,31(i4,1x))
114 format(11x,i4)
116 format(11x,29(i4,1x))
117 format(11x,29(i4,1x),5x,i4)
c-----7---------------------------------------------------------------72
c Will read every line until 'LAST' is found or end of file
do
read(unit=luls,fmt=111,end=1200,err=1200) iline
select case (iline(157:160))
! TIME
case('TIME')
cc read(iline,112) (itime(k),k=-2,mft)
read(iline,116) (itime(k),k=0,mft)
! DELV
case('DELV')
read(iline,114) iper
! INCV
case('INCV')
cc read(iline,112) (iincv(k),k=-2,mft)
read(iline,116) (iincv(k),k=0,mft)
! LAT
case('LAT ')
cc read(iline,112) (ilat(k),k=-2,mft)
read(iline,116) (ilat(k),k=0,mft)
! LON
case('LON ')
cc read(iline,112) (ilon(k),k=-2,mft)
read(iline,116) (ilon(k),k=0,mft)
! DTL
case('DTL')
read(iline,116) (idist(k),k= 0,mft)
! D20C
case('D20C')
cc read(iline,112) (id20c(k),k=-2,mft)
read(iline,116) (id20c(k),k=0,mft)
! D26C
case('D26C')
cc read(iline,112) (id26c(k),k=-2,mft)
read(iline,116) (id26c(k),k=0,mft)
! HCON
case('HCON')
cc read(iline,112) (ihcon(k),k=-2,mft)
read(iline,116) (ihcon(k),k=0,mft)
! U200
case('U200')
read(iline,116) (iu200(k),k= 0,mft)
! T200
case('T200')
read(iline,116) (it200(k),k= 0,mft)
! E000
case('E000')
read(iline,116) (ie000(k),k= 0,mft)
! EPOS
case('EPOS')
read(iline,116) (iepos(k),k= 0,mft)
! ENEG
case('ENEG')
read(iline,116) (ieneg(k),k= 0,mft)
! EPSS
case('EPSS')
read(iline,116) (iepss(k),k= 0,mft)
! ENSS
case('ENSS')
read(iline,116) (ienss(k),k= 0,mft)
! RHLO
case('RHLO')
read(iline,116) (irhlo(k),k= 0,mft)
! RHMD
case('RHMD')
read(iline,116) (irhmd(k),k= 0,mft)
! RHHI
case('RHHI')
read(iline,116) (irhhi(k),k= 0,mft)
! SHRD
case('SHRD')
read(iline,116) (ishrd(k),k= 0,mft)
! SHTD
case('SHTD')
read(iline,116) (ishtd(k),k= 0,mft)
c Convert ishtd from heading to direction
do k=0,mft
if (ishtd(k) .lt. imiss) then
ishtd(k) = 180 + ishtd(k)
if (ishtd(k) .gt. 360) ishtd(k) = ishtd(k)-360
endif
enddo
! SHRS
case('SHRS')
read(iline,116) (ishrs(k),k= 0,mft)
! SHTS
case('SHTS')
read(iline,116) (ishts(k),k= 0,mft)
c Convert ishts from heading to direction
do k=0,mft
if (ishts(k) .lt. imiss) then
ishts(k) = 180 + ishts(k)
if (ishts(k) .gt. 360) ishts(k) = ishts(k)-360
endif
enddo
! SHRG
case('SHRG')
read(iline,116) (ishrg(k),k= 0,mft)
! PSLV
case('PSLV')
read(iline,116) (ipslv(k),k= 0,mft)
! Z850
case('Z850')
read(iline,116) (iz850(k),k= 0,mft)
! D200
case('D200')
read(iline,116) (id200(k),k= 0,mft)
! REFC
case('REFC')
read(iline,116) (irefc(k),k= 0,mft)
! T000
case('T000')
read(iline,116) (it000(k),k= 0,mft)
! R000
case('R000')
read(iline,116) (ir000(k),k= 0,mft)
! Z000
case('Z000')
read(iline,116) (iz000(k),k= 0,mft)
! IR00
case('IR00')
read(iline,116) (igoes(k),k= 0,mft)
! IRM3
case('IRM3')
read(iline,116) (igoesm3(k),k= 0,mft)
! RD20
case('RD20')
read(iline,116) (ird20(k),k= 0,mft)
! RD26
case('RD26')
read(iline,116) (ird26(k),k= 0,mft)
! NOHC
case('NOHC')
read(iline,117) (irhcn(k),k= 0,mft),iohclag
! LAST
case('LAST')
exit
!other lines
case default
continue
end select
cycle
c-----7---------------------------------------------------------------72
1200 continue
exit
c-----7---------------------------------------------------------------72
enddo
1201 continue
endoffile=1
RETURN
END
c-----7---------------------------------------------------------------72
c-----7---------------------------------------------------------------72
c
subroutine read_radprof(lurad,raz,azav,std,pp50,irerrflg)
c Thomas Cram, May 2005
c
c Reads information from the IR brightness temperature radial
c profile data.
c
c INPUT
c lurad = unit number for radial profile data file
c
c OUTPUT
c raz = radii (km)
c azav = Tb profile
c std = standard deviation
c pp50 = percentage of -50C pixels
c irerrflg = 1.0 if error reading ir radial profile file
c
c-----7---------------------------------------------------------------72
parameter(maxrad=150)
real raz(maxrad),azav(maxrad),std(maxrad),
+ pp00(maxrad),pp10(maxrad),pp20(maxrad),pp30(maxrad),
+ pp40(maxrad),pp50(maxrad),pp60(maxrad),pp70(maxrad)
integer irerrflg
irerrflg=0.0
radtot = 0.0
do iprof=1,maxrad
read(unit=lurad,fmt=100,end=567,err=567)
+ raz(iprof),azav(iprof),std(iprof),
+ pp00(iprof),pp10(iprof),pp20(iprof),
+ pp30(iprof),pp40(iprof),pp50(iprof),
+ pp60(iprof),pp70(iprof)
radtot = radtot + 1
enddo
100 format(f6.1,1x,10(f6.2,1x))
567 continue
if (radtot.lt.maxrad) then
irerrflg = 1.0
endif
return
end
c-----7---------------------------------------------------------------72
c-----7---------------------------------------------------------------72
SUBROUTINE moment(data1,n,ave,adev,sdev,var,skew,curt)
c Given an array of DATA1 of length N, this routine returns its mean AVE,
c average deviation ADEV, standard deviation SDEV, variance VAR,
c skewness SKEW, and kurtosis CURT
INTEGER n
REAL adev,ave,curt,sdev,skew,var,data1(n)
INTEGER j
REAL p,s,ep
c if(n.le.1)pause 'n must be at least 2 in moment'
if (n.le.1) then
print*, 'n must be at least 2 in moment'
stop
endif
c First pass to get the mean
s=0.
do j=1,n
s=s+data1(j)
enddo
ave=s/n
c Second pass to get the first (absolute), second, third, and fourth
c moments of the deviation from the mean
adev=0.
var=0.
skew=0.
curt=0.
ep=0.
do j=1,n
s=data1(j)-ave
ep=ep+s
adev=adev+abs(s)
p=s*s
var=var+p
p=p*s
skew=skew+p
p=p*s
curt=curt+p
enddo
adev=adev/n
var=(var-ep**2/n)/(n-1)
sdev=sqrt(var)
if(var.ne.0.)then
skew=skew/(n*sdev**3)
curt=curt/(n*var**2)-3.
else
c pause 'no skew or kurtosis when zero variance in moment'
write(6,*) 'no skew,kurtosis when zero variance in subr moment'
endif
return
END
c-----7---------------------------------------------------------------72
c-----7---------------------------------------------------------------72
BLOCK DATA
c***********************************************************************
c Contains standard deviations, means, discriminant parameter
c weights, and discriminant function dividing value for use in the
c annular hurricane linear discriminant analysis.
c
c Order of array elements:
c weye, var, std, sst, u200
c
c Last Modified: Feb 16, 2007: JAK
c
c 1) Reduced number of factors to 5 from 9
c 2) replaced the means and standard deviations to represent the screened
c sample in the 1995-2006 developmental data
c************************************************************************
REAL sdev(5),means(5),dweights(5),div
common /MOMENTS/ sdev,means
common /DWEIGHTS/ dweights
common /DIVIDER/ div
data means/ 56.73, 558.21, 4.23,
+ 27.68, -8.09 /
data sdev/ 21.61, 218.52, 2.45,
+ 1.04, 5.62 /
data dweights/ 0.80944, 0.61429, -0.44537,
+ -0.80173, -0.14644/
data div/0.75904/
END
c-----7---------------------------------------------------------------72
SUBROUTINE ahindex(weye,var,std,gen_sst,u200,dv,val,prob)
c-----7---------------------------------------------------------------72
c Subroutine to apply disciminant weights contained in the BLOCK
c DATA routine to create discrimnant fucntion value (val) and
c a probability associated with val.
c
c Last Modified: Feburary 16, 2007 - JAK
c 1) now uses 5 predictors
c 2) new fit to the probability estimation
c***********************************************************************
REAL weye,var,std,gen_sst,u200
REAL dv,val,prob
! Discriminant weights (block data)
common /DWEIGHTS/ dweights(5)
! Discriminant function divider (block data)
common /DIVIDER/ div
c-----7---------------------------------------------------------------72
dv=0.0
val=0.0
dv=dweights(1)*weye + dweights(2)*var +
+ dweights(3)*std + dweights(4)*gen_sst + dweights(5)*u200
val=dv-div
! Determine annular hurricane probability
!
! quadratic function fits the observe PDF from the 1995-2006
! dependent data
if (val.lt.-0.3)then
prob=0.0
else if (val .ge. -0.3 .and. val .lt. 1.929)then
prob=-0.0183 * val**3 + 0.0689 * val**2
1 +0.2483 * val + 0.1347
else if (val .ge. 1.929) then
prob=1.0
endif
if (prob.gt.1.0) prob=1.0
RETURN
END
c
c-----7---------------------------------------------------------------72
c-----7---------------------------------------------------------------72
c
subroutine chartointmon(charmon,intmon)
integer intmon
character*3 charmon
select case (charmon)
case ('JAN')
intmon=1
case ('FEB')
intmon=2
case ('MAR')
intmon=3
case ('APR')
intmon=4
case ('MAY')
intmon=5
case ('JUN')
intmon=6
case ('JUL')
intmon=7
case ('AUG')
intmon=8
case ('SEP')
intmon=9
case ('OCT')
intmon=10
case ('NOV')
intmon=11
case ('DEC')
intmon=12
end select
END
c ++++END ANNULAR HURRICANE INDEX MODULE++++