subroutine ircal4(slat,vmax,itime,rad,vrad,mxr,mxv,idat,ierr)
c******
c PREAMBLE:
c
c NOTE:: storm latitude and vmax are added to the list of arguments
c
c This routine calculates parameters from the IR radial profiles
c and puts them in an integer array idat. The azimuthally
c averaged BT is converted to deg C. Additional parameters that
c are related to TC size are also calulated. TC size parameters are
c documented in Knaff et al. (2014, J. Clim), and Knaff et al. (2014,
c an extended abstract for the 31st AMS conference on Hurricanes
c and tropical meteorology)
c******
c USES:
c These are external
c SUBROUTINE r5wrapper -f90, calls the module r5utils
c These are contained in r5utils
c SUBROUTINE mk_12_1d_ir_PCs -f90, calculates the Principle Components
c (PCs) of radial Tb profiles
c REAL FUNCTION V500viaTb_sat -f90, calculates V500 (mean Vt @ 500 km)
c from latitude, PCs of the radial Tb profiles
c REAL FUNCTION R5_850 -f90, calculates R5 (TC size) from V500
c REAL FUNCTION r5mean -f90, calculates the R5 climatology as a
c function of intensity (vmax)
c
c******
c OUTPUT:
c
c ierr, INTEGER, 0 -no errors, 1- errors were encountered
c
c idat, INTEGER ARRAY (0:20):
c The calculated variables include the following:
c idat( 0) = Age of the GOES imagery relative to the storm case time (10*hr)
c idat( 1) = 0 to 200 scaled km radially averaged TB (10*deg C)
c idat( 2) = 0 to 200 scaled km radially averaged TB std dev (10*deg C)
c idat( 3) = 100 to 300 scaled km radially averaged TB (10*deg C)
c idat( 4) = 100 to 300 scaled km radially avareged TB std dev (10*deg C)
c idat( 5) = percent area from r=50 to 200 scaled km with TB < -10 C
c idat( 6) = same as (5) with TB < -20 C
c idat( 7) = same as (5) with TB < -30 C
c idat( 8) = same as (5) with TB < -40 C
c idat( 9) = same as (5) with TB < -50 C
c idat(10) = same as (5) with TB < -60 C
c idat(11) = max TB from 0 to 30 km radius (10*deg C)
c idat(12) = avg TB from 0 to 30 km radius (10*deg C)
c idat(13) = radius of max TB (km)
c idat(14) = min TB from 20 to 120 km radius (10*deg C)
c idat(15) = avg TB from 20 to 120 km radius (10*deg C)
c idat(16) = radius of min TB (km)
c idat(17) = IR-based estimate of V500
c idat(18) = IR-based estimate of R5
c idat(19) = IR-based size scaling factor - multiply by the real radius
c
c Note: TB = azimuthally averaged GOES channel 4 brightness temperature
c as a function of radius.
c******
c HISTORY:
c
c Modified 10/18/2013 to scale the averaging areas by an IR-based
c TC size estimate. This added one passed variable
c and several additional subroutines are called
c when compared to the ircal3 - JAK
c
c USES these external routines
c
c SUBROUTINE mk_12_1d_ir_PCs
c
c FUNCTIONS V500viaTb_sat, R5_850, r5mean
c
c Modified 05/22/2014 additional comments added to aid others - JAK
c Modified 11/21/2014 additional comments added, subroutines and function
c concatonatated to the end of this subroutine - JAK
c Modified 01/10/2017 adjusted idat assignment loop to 19 (instead of 16),
c fixed assignment of rproft last point - KDM
c
c******
c LAST MODIFIED: 01 January 2017
c
c*****************************************************************************
use r5utils
dimension rad(mxr),vrad(mxr,mxv)
dimension idat(0:19)
dimension rproft(0:150)
dimension pc(12)
c
c functions
c real V500viaTb_sat, R5_850, r5mean
c call r5wrapper
c
c Set min and max radii for BT and STD variables
c (Average BT and STD both calculated for both sets of radial intervals)
rmin1 = 0.0
rmax1 = 200.0
rmin2 = 100.0
rmax2 = 300.0
c
c Set min,max radius for pixel count variables
rmin = 50.0
rmax = 200.0
c
c Set flag for use of size adjustments to radii (1 = use TC size scaling)
isizerad = 0
c
tk = 273.15
imiss = 9999
c
idat( 0) = itime
do k=1,19
idat(k) = imiss
enddo
c
c Quality control the brightness temperature data
c
c Set min and max allowable Tb (deg C)
tminqc = -99.9
tmaxqc = 45.0
c
c Convert to K
tminqc = tminqc + tk
tmaxqc = tmaxqc + tk
c
c Check Tb. If any are bad skip this case
do i=1,mxr
t = vrad(i,1)
if (t .gt. 0.0) then
if (t .gt. tmaxqc .or. t .lt. tminqc) then
idat(0) = imiss
ierr=1
return
endif
endif
enddo
c
c Estimate the V500 and TC Size from the IR profile
do i=1,mxr
rproft(i-1)=vrad(i,1)
enddo
c account for the last point in the expected radial profile
rproft(mxr)=vrad(mxr-1,1)
npc=12
c
c 2-d principle components of the IR radial profile are calculated here
CALL mk_12_1d_ir_PCs(rproft,mxr,PC,npc,istat)
if (istat .ne. 0) then
print*,'pc calculations not made'
goto 900
endif
c
c size parameters are calculated here
V500=V500viaTb_sat(slat,PC(1),PC(2),PC(3))
R5=R5_850(V500)
size_scale=R5/r5mean(vmax)
c
c Adjust radii by TC size scaling if flag is set
if (isizerad .ne. 0) then
c
c Set min and max radii for BT and STD variables
c (Average BT and STD both calculated for both sets of radial intervals)
rmin1 = 0.0 * size_scale
rmax1 = 200.0 * size_scale
rmin2 = 100.0 * size_scale
rmax2 = 300.0 * size_scale
c
c Set min,max radius for pixel count variables
rmin = 50.0 * size_scale
rmax = 200.0 * size_scale
endif
c
c Calculate radial average t and s
w1 = 0.0
w2 = 0.0
c
t1 = 0.0
t2 = 0.0
s1 = 0.0
s2 = 0.0
c
do 15 i=1,mxr
t = vrad(i,1)
s = vrad(i,2)
r = rad(i)
c
if (t .lt. 0.0 .or. s .lt. 0.0) go to 15
c
if (r .gt. rmin1 .and. r .le. rmax1) then
t1 = t1 + t
s1 = s1 + s
w1 = w1 + 1.0
endif
c
if (r .gt. rmin2 .and. r .le. rmax2) then
t2 = t2 + t
s2 = s2 + s
w2 = w2 + 1.0
endif
15 continue
c
if (w1 .le. 0.0 .or. w2 .le. 0.0) go to 900
t1 = t1/w1 - tk
t2 = t2/w2 - tk
s1 = s1/w1
s2 = s2/w2
c
c Calculate pixel count variables
p1 = 0.0
p2 = 0.0
p3 = 0.0
p4 = 0.0
p5 = 0.0
p6 = 0.0
c
w1 = 0.0
w2 = 0.0
w3 = 0.0
w4 = 0.0
w5 = 0.0
w6 = 0.0
c
do i=1,mxr
r = rad(i)
p10 = vrad(i,4)
p20 = vrad(i,5)
p30 = vrad(i,6)
p40 = vrad(i,7)
p50 = vrad(i,8)
p60 = vrad(i,9)
c
if (p10 .ge. 0.0 .and. r .gt. rmin .and. r .le. rmax) then
w1 = w1 + r
p1 = p1 + r*p10
endif
c
if (p20 .ge. 0.0 .and. r .gt. rmin .and. r .le. rmax) then
w2 = w2 + r
p2 = p2 + r*p20
endif
c
if (p30 .ge. 0.0 .and. r .gt. rmin .and. r .le. rmax) then
w3 = w3 + r
p3 = p3 + r*p30
endif
c
if (p40 .ge. 0.0 .and. r .gt. rmin .and. r .le. rmax) then
w4 = w4 + r
p4 = p4 + r*p40
endif
c
if (p50 .ge. 0.0 .and. r .gt. rmin .and. r .le. rmax) then
w5 = w5 + r
p5 = p5 + r*p50
endif
c
if (p60 .ge. 0.0 .and. r .gt. rmin .and. r .le. rmax) then
w6 = w6 + r
p6 = p6 + r*p60
endif
enddo
c
if (w1 .le. 0.0 .or. w2 .le. 0.0 .or. w3 .le. 0.0) go to 900
if (w4 .le. 0.0 .or. w5 .le. 0.0 .or. w6 .le. 0.0) go to 900
c
p1 = p1/w1
p2 = p2/w2
p3 = p3/w3
p4 = p4/w4
p5 = p5/w5
p6 = p6/w6
c
c write(6,*) rmin,tmin,smin,tinn
c
c Find eye and cold ring variables
temax = -9999.0
teavg = 0.0
remax = 0.0
nre = 0
c
trmin = 9999.0
travg = 0.0
rrmin = 0.0
nrr = 0
c
r1e = 0.0
r2e = 30.0
r1r = 20.0
r2r = 120.0
c
do 35 i=1,mxr
t = vrad(i,1)
r = rad(i)
if (r .le. r1e) go to 35
if (r .gt. r2e) go to 35
if (t .lt. 0.0) go to 35
c
nre = nre + 1
teavg = teavg + t
c
c Find max eye temperature
if (t .gt. temax) then
temax = t
remax = r
endif
35 continue
if (nre .lt. 1) go to 900
teavg = (teavg/float(nre)) - tk
temax = temax - tk
c
do 45 i=1,mxr
t = vrad(i,1)
r = rad(i)
if (r .le. r1r) go to 45
if (r .gt. r2r) go to 45
if (t .lt. 0.0) go to 45
c
nrr = nrr + 1
travg = travg + t
c
c Find min ring temperature
if (t .lt. trmin) then
trmin = t
rrmin = r
endif
45 continue
if (nrr .lt. 1) go to 900
travg = (travg/float(nrr)) - tk
trmin = trmin - tk
c
idat(1) = int(10.0*t1)
idat(2) = int(10.0*s1)
idat(3) = int(10.0*t2)
idat(4) = int(10.0*s2)
idat(5) = int(p1)
idat(6) = int(p2)
idat(7) = int(p3)
idat(8) = int(p4)
idat(9) = int(p5)
idat(10) = int(p6)
idat(11) = int(10.0*temax)
idat(12) = int(10.0*teavg)
idat(13) = int(remax)
idat(14) = int(10.0*trmin)
idat(15) = int(10.0*travg)
idat(16) = int(rrmin)
idat(17) = int(V500*10)
idat(18) = int(R5*10)
idat(19) = int(size_scale*100)
c
ierr = 0
return
c
900 write(6,*) 'fatal error in routine ircal4'
ierr = 1
return
c
end