program oparet
c This program performs retrievals of temperature, heights,
c surface pressure and winds using the AMSU temperatures at
c the satellite footprints.
c
c This is the updated version where all satellite and tropical cyclone
c information is in a single file.
c
c File Updates:
c New version created 5/15/2003 M. DeMaria
c ATCF fix file option added 9/26/2005 M. DeMaria and J. Knaff
c DATELINE ISSUE Resolved 2/20/2006 J. Knaff
c
c
c Input files:
c temp_ret.dat - File with storm input parameters
c and AMSU-A radiances and temp/CLW retrievals
c (combination of temp, times and coordinate
c files, in the old version of oparet)
c
c AVN.DAT - Packed ASCII NCEP analysis file
c
c Output files:
c oparet.log - ASCII file with basic information about the
c retrieval calcuations
c Lbnnyy_mmddtt.DATss - loc: ASCII file containing lat/lons of AMSU input data
c
c Rbnnyy_mmddtt.DATss - rza: ASCII file containing V,T,P,rho as function of
c (r,z) from AMSU retrieval
c Abnnyy_mmddtt.DATss - xya: ASCII file containing U,V,T,Z,Ps,rho,CLW
c as function of (x,y,P) from AMSU retrieval
c
c Sbnnyy_mmddtt.DATss - sta: ASCII file containing statistical parameters for
c the TC estimation algorithm
c Fbnnyy_mmddtt.DATss - fix: ASCII file containing estimated TC intensity/size
c
c Gbnnyy_mmddtt.DATss - afx: ATCF formatted fix file
c
c Abnnyy_XmmDD_PACK.DAT - apk: Packed ASCII file containing U,V,T,Z,Ps
c as function of (x,y,P) from the AMSU retrieval
c Nbnnyy_XmmDD_PACK.DAT - npk: Packed ASCII file containing U,V,T,Z,Ps
c as function of (x,y,P) from NCEP fields
c interpolated to the AMSU analysis grid
c
c All of the output file names can be replaced by generic file
c names by setting the parameter igenfn=1. The files will be of the
c for NOAAss.olg, .loc, .rza, etc according the to above list.
c
c Note: b = basin ID (A=Atlantic, E=East Pacific, W=West Pacific)
c nn = ATCF storm number
c yy = year
c mm = month
c dd = day
c tt = time in UTC of ATCF record used for storm center
c ss = satellite number (15=NOAA15, etc)
c
c X = X if 00 or 12 UTC analysis
c = Y if 06 or 18 UTC analysis
c DD = dd if 00 or 12 UTC analysis
c = dd+50 if 06 or 18 UTC analysis
c
c Specify the lon/lat dimensions of the AMSU analysis domain
parameter (nx=61,ny=61)
c
c Specify number of AMSU pressure levels for calculations
c and first pressure (out of 40, see array pamsu) to use
parameter (np=23,npst=16)
c
c Specify number of NCEP pressure levels
parameter (npn=12)
c
c Specify maximum dimensions of NCEP grid arrays
parameter(ixmax=181,iymax=91)
c
c Specify max number of AMSU swath points and pressure levels
parameter (mxas=7000,mpas=40)
c
c Specify number of radial and height points for gradient wind
c calculations
parameter (nr=31,nz=21)
c
c Arrays for AMSU variables on analysis grid
dimension p(np)
dimension z(nx,ny,np),t(nx,ny,np)
dimension us(nx,ny),vs(nx,ny),ps(nx,ny),ts(nx,ny),rhos(nx,ny)
dimension clwxy(nx,ny)
c
c Array for AMSU swath variables
dimension tas(mxas,mpas),aslat(mxas),aslon(mxas)
dimension dumasp(mxas,mpas),dumas(mxas),idumas(mxas)
dimension clwas(mxas),tpwas(mxas),pamsu(mpas)
c
c Arrays for AMSU output (at NCEP pressure levels)
dimension ua(nx,ny,npn),va(nx,ny,npn),za(nx,ny,npn),ta(nx,ny,npn)
c
c Arrays for NCEP analysis variables
dimension pn(npn)
dimension un(nx,ny,npn),vn(nx,ny,npn),zn(nx,ny,npn),tn(nx,ny,npn)
dimension psn(nx,ny),tsn(nx,ny)
c
c Array for NCEP 1000 mb T for clw temperature adjustment routine
dimension tn1000(ixmax,iymax)
c
c Arrays for gradient wind calculations
dimension rr(nr),trp(nr,np),psr(nr),tsr(nr)
dimension zrp(nr,np),vrp(nr,np),zz(nz)
dimension trz(nr,nz),prz(nr,nz),rhorz(nr,nz),vrz(nr,nz)
dimension tarz(nr,nz)
c
c Lat/Lon arrays
dimension rlond(nx),rlatd(ny)
dimension rlonr(nx),rlatr(ny)
dimension sinlat(ny),coslat(ny),tanlat(ny)
c
c Temporary arrays for radial Barnes analysis
parameter (mxb=nx*ny)
dimension temlat(mxb),temlon(mxb),ftem(mxb)
c
c Array for predictors for the stats file
parameter (npred=18)
dimension preds(npred)
c
c Arrays for statistical wind radii prediction
dimension pr34(6),pr50(6),pr64(6)
dimension ir34(4),ir50(4),ir64(4)
c
c Work arrays
dimension work1(nx,ny),work2(nx,ny),work3(nx,ny)
dimension work4(nx,ny),work5(nx,ny),work6(nx,ny)
dimension worky(ny)
dimension works(mxas)
c
c Arrays for Cartesian grid calculations
dimension x(nx),y(ny),fy(ny),ibwin(npn)
c
character *1 rawdid
character *2 cnsat,cbasin
character *3 lab3
character *6 atcfid,noaass,cstype
character *9 sname
character *14 ctimdat
character *20 label
character *23 labelr
c
c Variables for input file names
character *25 fntinp,fnninp
c
c Variables for output file names
character *25 fnlog,fnloc
character *25 fnnpk,fnapk
character *25 fnrza,fnxya
character *25 fnsta,fnfix,fnafx
c
c Variables for output file headers
character *90 coord,times
c
c Variable for temporary file names
character *25 fntemp
c
common /cons/ pi,g,rd,dtr,erad,erot
common /cgrid/ x,y,fy,dx,dy,f0,beta
common /log/ lulog
common /ncepfg/ rlonln,rlatbn,dlonn,dlatn,nlonn,nlatn
common /ncepll/ rlond,rlonr,rlatd,rlatr,dlon,dlat,
+ dlonr,dlatr,sinlat,coslat,tanlat
common /ncepp/ pn
common /sinfo1/ atcfid,sname
common /sinfo2/ iyr,imon,iday,itime4,
+ jyr,jmon,jday,jtime4,
+ ivmax,spd,head,
+ slat00,slat12,sslat,
+ slon00,slon12,sslon,
+ rmp,xp,pr34,pr50,pr64
c
c Unit numbers for input files
data lutinp,luninp
+ /21,22/
c
c Unit numbers for output files
data luloc,lulog,lunpk,luapk,lurza,lusta,lufix,luafx,luxya
+ /31,32,33,34,35,36,37,38,39/
c
c Constants
data pi, dtr, g, rd, erad, erot
+ /3.14159, 0.0174533, 9.81, 287.0, 6371.0e+3, 7.292e-5/
c
c Specify input file with storm parameters and AMSU temperatures
data fntinp /'temp_ret.dat'/
c
c NCEP input data file name
data fnninp /'AVN.DAT'/
c
c Output log file name
data fnlog /'oparet.log'/
c
c Total AMSU pressure levels
data pamsu / 0.1, 0.2, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 7.0,
+ 10., 15., 20., 25., 30., 50., 60., 70., 85.,100.,
+ 115.,135.,150.,200.,250.,300.,350.,400.,430.,475.,
+ 500.,570.,620.,670.,700.,780.,850.,920.,950.,1000./
c
c NCEP pressure levels
data pn / 100., 150., 200., 250., 300., 400.,
+ 500., 600., 700., 850., 925.,1000./
c
c Flag array for calculating balanced winds at each NCEP pressure level.
c The calculation is performed (skipped) if ibwin=1 (0).
data ibwin / 1, 1, 1, 1, 1, 1,
+ 1, 1, 1, 1, 1, 1/
c
c Open log file
open(unit=lulog,file=fnlog,form='formatted',
+ status='unknown')
c
c **** Begin Parameter Specification ****
c
c ++ 3-D analysis grid parameters
c
c Specify lat/lon spacing of AMSU analysis grid
dlat = 0.20
dlon = 0.20
c
c Set idex=1 to exclude data outside of analysis domain. Also
c set lat/lon buffers for exclusion
idex = 1
rlonbf =3.0
rlatbf =3.0
c
c Specify maximum distance (km) the storm can be from the swath center
c to perform an analysis
threshr =700.0
c
c ++ Gradient wind parameters
c
c Set iaxsym =1 for axisymmetric analysis
iaxsym=1
c
c Specify the radial and height spacing (m) for gradient wind calculations
dr = 20.0e+3
dz = 1.0e+3
c
c ++ 3-D wind parameters
c Set iwinda = 0 to use NCEP x,y winds
c = 1 to use linear balance equation for x,y winds
c = 2 to use nonlinear balance equation for x,y winds
c (iterative solution)
c = 3 to use nonlinear balance equation for x,y winds
c (variational solution, u,v form)
c = 4 to use nonlinear balance equation for x,y winds
c (variational solution, psi form)
iwinda=3
c
c If iwinda = 2, 3 or 4 , select method for first guess u,v:
c Set ifgnbe = 0 to use NCEP winds as first guess
c Set ifgnbe = 1 to use linear balance winds at first guess
c Set ifgnbe = 2 to use zero curvature winds as first guess
ifgnbe=0
c
c ++ Set flags for output files (0=don't write file, 1=write file)
c
c Set igenfn=1 to use generic output file names, otherwise igenfn=0
igenfn=1
c
c fix file (requires iaxsym=1, ifixall=1 writes fix file even if TC
c analysis is not performed)
ipfix = 1
ifixall = 1
c
c loc file
iploc = 1
c
c rza file (requires iaxsym=1)
iprza=1
c
c xya file
ipxya = 1
c
c sta file (requires iaxsym=1 and ivrapd=1)
ipsta = 1
c
c apk file
ipapk = 0
c
c npk file
ipnpk = 0
c
c ++ Barnes analysis parameters
c
c Specify e-folding and influence radii (km) for 2-D Barnes analysis
c of AMSU temperatures.
c Rule of thumb: rinfxy = efldxy*5, but rinfxy must be at least 400
c efldxy=40.0
c rinfxy=200.0
c efldxy=100.0
efldxy=100.0
rinfxy=6.0*efldxy
c
c Set ispf=1 for a second pass of the x,y Barnes analyses or
c ispf=0 for single pass
ispf=1
c
c Specify e-folding and influence radii (km) for 1-D Barnes analysis
c of AMSU temperatures. Also specify expansion factor (exfac) for
c increasing e-folding radius as a function of radius.
efldr=20.0
rinfr=100.0
c efldr=120.0
c rinfr=480.0
exfac=0.0
c
c Set iradxy=1 to perform radial Barnes temperature analysis using
c gridded temperature from x,y analysis instead of swath data
iradxy=1
c
c Set iefa = 1 to adjust e-folding radii based upon horizontal
c resolution, else set iefa=0. Also, specify reduction
c factor for 2dx wave, to choose e-folding radius.
iefa = 0
c tdxfxy = 0.05
c tdxfr = 0.20
tdxfxy = 0.04
tdxfr = 0.16
c
c ++ CLW attenuation and ice scattering corretion parameters
c
c Set itcor = 0-10 to correct cold anomalies for water
c vapor attenuation.
c If itcor=0 no correction is applied
c If itcor=1 temp analysis is corrected, using a
c cold anomaly threshold
c If itcor=2 original swath data is corrected using a
c cold anomaly threshold
c If itcor=3 original swath data is corrected using a
c linear regression technique
c If itcor=4 original swath data is corrected using the
c linear regression method (3), followed
c by the cold anomaly threshold method (2).
c If itcor=5 orginal swath data is adjusted in regions
c with high liquid water content to a
c constant lapse rate sounding
c If itcor=6 method 3 is applied, followed by method 5
c If itcor=7 method 3, 5 and 2 are applied
c If itcor=8 temp analysis corrected; quadratically
c smoothed m(k) from Tdev vs CLW regression
c slopes; m(k)=ap^2+bp+c where p=pressure in Pa
c and b=c=0 according to initial conditions.
c If itcor=9 temp analysis corrected, for ice crystals
c via B. Linstid's algorithm
c If itcor=10 temp analysis correction via method (8) and (9)
c
c ncemx,ncsmx = first,last AMSU level (out of analysis levels)
c to apply correction (x=2 applies to method 1 or 2,
c x=3 applies to method 3,
c x=5 applies to method 5,
c x=8 applies to method 8 or 10,
c x=9 applied to method 9 or 10)
c dt = anomaly threshold for applying correction
c (for itcor=1 or 2)
c dtred = Factor to reduce fraction of cold anomaly below dt
c (for itcor=1 or 2)
c tcrmax = Max radius (km) from storm for applying method 3
c
itcor=10
c
c Specify variables for method 1 and 2
ncsm2 = np-5
ncem2 = np
dt = -0.5
dtred = 0.10
tcrmax = 900.0
c
c Specify variables for method 3
ncsm3 = 11
ncem3 = np
c
ncs2 = ncsm2 + npst - 1
nce2 = ncem2 + npst - 1
ncs3 = ncsm3 + npst - 1
nce3 = ncem3 + npst - 1
c
c Specify variables for method 5
ncsm5 = 13
ncem5 = np
clwth1 = 1.00
clwth2 = 2.00
c
ncs5 = ncsm5 + npst - 1
nce5 = ncem5 + npst - 1
c
c Specify variables for method 8
c Variables for Ben Linstid's original correction
c ncsm8 = 10
c ncem8 = np
c amkx = 5.78
c xlambda = 1.55
c
c Variables for Tdev vs CLW regression correction on grid
c The first amkx value is with no core removed from the grid data,
c and the second is with a 2 x 2 degree core removed (comment one out).
c ncsm8 = 12
c ncem8 = np
c amkx = 4.27435
c amkx = 4.04594
c xlambda = 1.0
c
c Variables for Tdev vs CLW regression correction on swath
ncsm8 = 12 + npst - 1
ncem8 = np + npst - 1
amkx = 4.57166
xlambda = 1.0
c
c Specify variables for method 9
ncsm9 = 12
ncem9 = np
c
c ++ NCEP analysis parameters
c
c Set ismooth to the number of times to smooth the NCEP
c fields after interpolating them to the AMSU analysis grid
ismooth=3
c
c Set intop=1 to use NCEP height field at top of AMSU domain
c as boundary condition. This is only possible if AMSU domain
c top is at 100 mb.
intop=0
c
c Set ibrem=1 to correct amsu T,Z and surface P to make the mean
c equal to that in the NCEP analysis
ibrem = 0
c
c **** End Parameter Specification ****
c
c **** Begin Storm/AMSU Data Input File Read and Related Processing ****
c
c Open the storm parameter and temperature retrieval input file
fntemp=fntinp
open(unit=lutinp,file=fntinp,form='formatted',
+ status='old',err=900)
c
c Read the storm/retrieval file
call srread(lutinp,jyr,jlday,jtime,swlat,swlon,
+ iyr,imon,iday,itime,slat00,slon00,slat12,slon12,
+ idir,ispeed,ivmax,atcfid,sname,nsat,
+ nxas,aslat,aslon,tas,clwas,tpwas,
+ coord,times,istat)
close(lutinp)
c
c J. Knaff -- DATELINE ISSUE resolved 2/20/2006
c This fix allows for the use of data across the dateline when storm
c is within 10 degrees of the dateline by using a 0 to 360 longitude
c convention in that region.
c Fix
if (abs(slon00).gt.170.0 .or. abs(slon12).gt.170.0)then
if(slon00 .lt. 0.0 ) slon00 = slon00 + 360.0
if(slon12 .lt. 0.0 ) slon12 = slon12 + 360.0
if(swlon .lt. 0.0 ) swlon = swlon + 360.0
do ikn=1,nxas
if(aslon(ikn) .lt. 0.0) aslon(ikn) = aslon(ikn) + 360.0
enddo
endif
c
c
if (istat .ne. 0) go to 900
c
call ucase(atcfid,6)
cbasin=atcfid(1:2)
iyr2 = iyr - 100*(iyr/100)
c
write(cnsat,309) nsat
309 format(i2.2)
noaass(1:4) = 'NOAA'
noaass(5:6) = cnsat
c
call jdayi(jlday,jyr,jmon,jday)
c
c Extract AMSU pressures for calculations and convert to Pa
do 5 k=1,np
p(k) = 100.0*pamsu(npst+k-1)
5 continue
c
c Estimate storm center at AMSU swath center time
jtimeh = jtime/10000
jtimem = (jtime - 10000*jtimeh)/100
jtime4 = 100*jtimeh + jtimem
c
call tdiff(jyr,jmon,jday,jtimeh,iyr,imon,iday,itime,idelt)
adelt = float(idelt) + float(jtimem)/60.0
itime4 = 100*itime
c
if (abs(slat12) .gt. 0.0 .and.
+ abs(slon12) .gt. 0.0 ) then
w0 = 1.0 + adelt/12
wm = -adelt/12
else
w0 = 1.0
wm = 0.0
endif
c
sslat = w0*slat00 + wm*slat12
sslon = w0*slon00 + wm*slon12
c
c Calculate distance from storm to swath center
call distk(swlon,swlat,sslon,sslat,dxk,dyk,rad)
radsts=rad
c
c Estimate AMSU resolution
call amsures(rad,ares)
c
if (iefa .eq. 1 .and. ares .gt. 0.0) then
c Adjust e-folding radii based upon horizontal resolution of the data
cofr = 2.0*sqrt(alog(1./tdxfr))/pi
cofxy = 2.0*sqrt(alog(1./tdxfxy))/pi
efldxy = float(ifix(cofxy*ares))
efldr = float(ifix(cofr *ares))
c
if (efldxy .lt. 35.0) efldxy = 35.0
if (efldxy .gt. 150.0) efldxy = 150.0
c
if (efldr .lt. 35.0) efldr = 35.0
if (efldr .gt. 150.0) efldr = 150.0
c
rinfxy = 5.0*efldxy
rinfr = 5.0*efldr
if (rinfxy .lt. 400.0) rinfxy = 400.0
if (rinfr .lt. 400.0) rinfr = 400.0
endif
c
c Write storm parameter info to log file
write(lulog,200) atcfid,sname,iyr,imon,iday,itime,
+ slat00,slon00,ivmax,slat12,slon12
200 format(/,' Start program oparet. ',
+ /,' ATCF INPUT: ',a6,1x,a9,1x,i4.4,1x,i2.2,i2.2,1x,
+ i2.2,' UTC',
+ /,' Storm lat/lon (t= 0 hr): ',f5.1,1x,f6.1,1x,
+ ' Vmax: ',i3,
+ /,' Storm lat/lon (t=-12 hr): ',f5.1,1x,f6.1)
c
write(lulog,202) np,p(1)/100.,p(np)/100.
202 format(/,i2,' AMSU levels from ',f5.0,' to ',f5.0,' included')
c
write(lulog,204) efldxy,rinfxy
204 format(/,' Cartesian Barnes analysis with e-fold, inf. radii: ',
+ f5.0,1x,f5.0)
c
if (iaxsym .eq. 1) then
write(lulog,206) efldr,rinfr,exfac
206 format(' Radial Barnes analysis with e-fold, inf. radii: ',
+ f5.0,1x,f5.0,' exfac=',f4.1)
endif
c
if (itcor .gt. 0) then
write(lulog,207) itcor,p(ncsm2)/100.,p(ncem2)/100.,
+ p(ncsm3)/100.,p(ncem3)/100.,
+ p(ncsm5)/100.,p(ncem5)/100.,
+ clwth1,clwth2
207 format(/,' Cold anomalies adjusted using method ',i2,
+ /,' pm2= ',f5.0,' to ',f5.0,
+ /,' pm3= ',f5.0,' to ',f5.0,
+ /,' pm5= ',f5.0,' to ',f5.0,
+ /,' clwth1=',f5.2,' clwth2=',f5.2)
else
write(lulog,*) ' itcor=0, no correction applied'
endif
c
write(lulog,210) swlat,swlon,jyr,jmon,jday,jtime
210 format(/,'AMSU data swath center: ',f6.2,1x,f7.2,' at ',
+ i4,1x,i2.2,i2.2,1x,i6)
c
write(lulog,212) adelt,sslat,sslon,rad,ares,nxas,nsat
212 format('Swath data is ',f5.1,' hr from analysis time',/,
+ 'Storm center at swath time:',f6.2,1x,f7.2,/,
+ 'Distance (km) from storm to swath center: ',f6.0,/,
+ 'Approximate data spacing at storm center: ',f6.0,/,
+ 'AMSU data points read from input, n=: ',i5,/,
+ 'Data from NOAA',i2.2)
c
c Check to make sure storm is not too close to the swath edge
inhce = 0
if (radsts .gt. threshr) then
write(lulog,905) radsts,threshr
905 format(/,' Storm too far from swath center, r (km)=',f6.0,
+ ' max allowable r (km) = ',f5.0)
inhce = 1
c
if (ifixall .eq. 1) go to 5000
stop
endif
c
c **** End Storm/AMSU Data Input File Read and Related Processing ****
c
c **** Begin Calculation of Analysis Grid Parameters ****
c Center the domain on the estimated storm position
c at the swath time
islat = ifix(sslat+0.5)
islon = ifix(sslon-0.5)
c
rlonl = float(islon) - dlon*float(nx/2)
rlatb = float(islat) - dlat*float(ny/2)
c
c Calculate lat and lon in deg and radians, and related variables
dlonr = dlon*dtr
dlatr = dlat*dtr
c
do 10 i=1,nx
rlond(i) = rlonl + dlon*float(i-1)
rlonr(i) = dtr*rlond(i)
10 continue
c
do 12 j=1,ny
rlatd(j) = rlatb + dlat*float(j-1)
rlatr(j) = dtr*rlatd(j)
c
sinlat(j) = sin(rlatr(j))
coslat(j) = cos(rlatr(j))
tanlat(j) = tan(rlatr(j))
12 continue
c
c Calcuate Cartesian grid variables
reflat = 0.5*(rlatd(1)+rlatd(ny))
reflon = 0.5*(rlond(1)+rlond(nx))
c
do j=1,ny
y(j) = erad*dtr*(rlatd(j)-reflat)
enddo
c
crl = cos(dtr*reflat)
do i=1,nx
x(i) = erad*dtr*(rlond(i)-reflon)*crl
enddo
c
dx = x(2)-x(1)
dy = y(2)-y(1)
f0 = 2.0*erot*sin(dtr*reflat)
beta = 2.0*erot*cos(dtr*reflat)/erad
c
do j=1,ny
fy(j) = f0 + beta*y(j)
enddo
c
c Calculate radial and height grids
do 13 i=1,nr
rr(i) = dr*float(i-1)
13 continue
c
do 14 m=1,nz
zz(m) = dz*float(m-1)
14 continue
c
c Write grid info to the log file
write(lulog,230) rlond(1),rlond(nx),dlon,
+ rlatd(1),rlatd(ny),dlat
230 format(/,' Longitude domain: ',f6.1,1x,f6.1,' dlon=',f5.2,/,
+ ' Latitude domain: ',f6.1,1x,f6.1,' dlat=',f5.2)
c
if (iaxsym .eq. 1) then
write(lulog,232) nr,dr/1000.0,rr(1)/1000.0,rr(nr)/1000.0
232 format(i3,' radial grid points, dr=',f5.1,1x,
+ ' rmin=',f5.1,' rmax=',f6.1)
endif
c
write(lulog,234) nz,dz/1000.0,zz(1)/1000.0,zz(nz)/1000.0
234 format(i3,' z grid points, dz=',f5.1,1x,
+ ' zmin=',f5.1,' zmax=',f6.1)
c
write(lulog,239) x(1),x(nx),dx,y(1),y(ny),dy,f0,beta
239 format(/,' Cartesian grid created for x,y balanced winds: ',
+ /,' x domain: ',e11.4,' to ',e11.4,' dx=',e11.4,
+ /,' y domain: ',e11.4,' to ',e11.4,' dy=',e11.4,
+ /,' f0=',e11.4,' beta=',e11.4)
c
if (iwinda .eq. 4) then
write(lulog,244) ifgnbe
244 format(/,' AMSU x,y winds from nonlinear balance equation',
+ /,' Variational solution (psi form) with ifgnbe= ',i1)
elseif (iwinda .eq. 3) then
write(lulog,243) ifgnbe
243 format(/,' AMSU x,y winds from nonlinear balance equation',
+ /,' Variational solution (u,v form) with ifgnbe= ',i1)
elseif (iwinda .eq. 2) then
write(lulog,242) ifgnbe
242 format(/,' AMSU x,y winds from nonlinear balance equation',
+ /,' Iterative solution with ifgnbe= ',i1)
elseif (iwinda .eq. 1) then
write(lulog,241)
241 format(/,' AMSU x,y winds from linear balance equation')
else
write(lulog,240)
240 format(/,' AMSU x,y winds from NCEP analysis')
endif
c
c **** End Calculation of Analysis Grid Parameters ****
c
c **** Begin NCEP Input File Read and Related Processing ****
c
c Open NCEP data input file
fntemp=fnninp
open(unit=luninp,file=fnninp,form='formatted',status='old',
+ err=900)
c
c Get NCEP analysis variables
call ncepget(pn,un,vn,zn,tn,tn1000,rlond,rlatd,nx,ny,npn,
+ dlon,dlat,dayx,utcx,luninp,lulog,ierr)
close(luninp)
c
if (ierr .ne. 0) then
write(lulog,910) ierr
910 format(/,' Halting due to error in routine ncepget, ierr=',
+ i2)
endif
c
if (ismooth .gt. 0) then
do n=1,ismooth
do k=1,npn
call smooth(un(1,1,k),work1,nx,ny,nx,ny,0)
call smooth(vn(1,1,k),work1,nx,ny,nx,ny,0)
call smooth(tn(1,1,k),work1,nx,ny,nx,ny,0)
call smooth(zn(1,1,k),work1,nx,ny,nx,ny,0)
enddo
enddo
endif
c
c Convert NCEP date/time varibles to integers
idayx = ifix(dayx)
iutcx = ifix(utcx)
c
nyrt = idayx/10000
nmon = (idayx - 10000*nyrt)/100
nday = (idayx - 10000*nyrt - 100*nmon)
ntime = iutcx/100
c
if (nyrt .lt. 50) then
nyr = nyrt + 2000
else
nyr = nyrt + 1900
endif
c
c Calculate time difference between NCEP analysis to AMSU swath data
call tdiff(jyr,jmon,jday,jtimeh,nyr,nmon,nday,ntime,idelt)
adelt = float(idelt) + float(jtimem)/60.0
c
write(lulog,216) fnninp
216 format(/,'NCEP data input file: ',a25)
c
write(lulog,236) nyr,nmon,nday,ntime,adelt
236 format(/,'NCEP analysis date/time: ',i4,1x,2(i2.2),1x,i2,/,
+ 'NCEP analysis is ',f5.1,' hrs from AMSU swath')
c
if (abs(adelt) .gt. 36.0) then
write(lulog,*) ' NCEP analysis is too old'
stop
endif
c
c Convert NCEP pressures to Pa
do 15 k=1,npn
pn(k) = 100.0*pn(k)
15 continue
c
c Use data at lowest two NCEP levels to get surface temperature
c and pressure
do 17 j=1,ny
do 17 i=1,nx
z1 = zn(i,j,npn)
z2 = zn(i,j,npn-1)
t1 = tn(i,j,npn)
t2 = tn(i,j,npn-1)
p1 = pn(npn)
call pstcal(z1,z2,t1,t2,p1,pstem,tstem)
c
psn(i,j) = pstem
tsn(i,j) = tstem
17 continue
c
c **** End NCEP Input File Read and Related Processing ****
c
c ***** Begin Output File Naming/Opening ****
c
c ++ Name for loc file (AMSU footprint locations)
fnloc = 'LX0000_000000.DAT'
fnloc( 2: 2) = atcfid(1:1)
fnloc( 3: 6) = atcfid(3:6)
fnloc(18:19) = cnsat
write(fnloc(8:13),226) imon,iday,itime
226 format(3(i2.2))
c
if (igenfn .eq. 1) then
fnloc=' '
fnloc(1:6) = noaass
fnloc(7:10) = '.LOC'
endif
c
c ++ Name for rza file (2-D AMSU retrieval output file)
fnrza = 'RX0000_000000.DAT'
fnrza( 2: 2) = atcfid(1:1)
fnrza( 3: 6) = atcfid(3:6)
fnrza(18:19) = cnsat
write(fnrza(8:13),226) imon,iday,itime
c
if (igenfn .eq. 1) then
fnrza=' '
fnrza(1:6) = noaass
fnrza(7:10) = '.RZA'
endif
c
c ++ Name for xya file (3-D AMSU retrieval output file)
fnxya = 'AX0000_000000.DAT'
fnxya( 2: 2) = atcfid(1:1)
fnxya( 3: 6) = atcfid(3:6)
fnxya(18:19) = cnsat
write(fnxya(8:13),226) imon,iday,itime
c
if (igenfn .eq. 1) then
fnxya=' '
fnxya(1:6) = noaass
fnxya(7:10) = '.XYA'
endif
c
c ++ Name for sta file (TC statistics file)
fnsta = 'SX0000_000000.DAT'
fnsta( 2: 2) = atcfid(1:1)
fnsta( 3: 6) = atcfid(3:6)
fnsta(18:19) = cnsat
write(fnsta(8:13),226) imon,iday,itime
c
if (igenfn .eq. 1) then
fnsta=' '
fnsta(1:6) = noaass
fnsta(7:10) = '.STA'
endif
c
c ++ Create file name for fix file
fnfix = 'FX0000_000000.DAT'
fnfix( 2: 2) = atcfid(1:1)
fnfix( 3: 6) = atcfid(3:6)
fnfix(18:19) = cnsat
write(fnfix(8:13),226) imon,iday,itime
c
if (igenfn .eq. 1) then
fnfix=' '
fnfix(1:6) = noaass
fnfix(7:10) = '.FIX'
endif
c
c ++ Create file name for afx file
fnafx = 'GX0000_000000.DAT'
fnafx( 2: 2) = atcfid(1:1)
fnafx( 3: 6) = atcfid(3:6)
fnafx(18:19) = cnsat
write(fnafx(8:13),226) imon,iday,itime
c
if (igenfn .eq. 1) then
fnafx=' '
fnafx(1:6) = noaass
fnafx(7:10) = '.AFX'
endif
c
c ++ Name for npk file (packed NCEP output file)
if (itime .eq. 6 .or. itime .eq. 18) then
rawdid = 'Y'
else
rawdid = 'X'
endif
c
if (itime .ge. 12) then
idayr = iday+50
else
idayr = iday
endif
c
fnnpk = ' '
fnnpk(1:1) = 'N'
fnnpk(2:2) = atcfid(1:1)
fnnpk(3:6) = atcfid(3:6)
fnnpk(7:21) = '_X0000_PACK.DAT'
fnnpk(22:23)= cnsat
write(fnnpk(8:12),300) rawdid,imon,idayr
300 format(a1,i2.2,i2.2)
c
if (igenfn .eq. 1) then
fnnpk=' '
fnnpk(1:6) = noaass
fnnpk(7:10) = '.NPK'
endif
c
c ++ Name for apk file (packed 3-D AMSU retrieval output file)
c on lat,lon,P grid
fnapk=fnnpk
fnapk(1:1) = 'A'
c
if (igenfn .eq. 1) then
fnapk=' '
fnapk(1:6) = noaass
fnapk(7:10) = '.APK'
endif
c
c ++ Write file names to log file
lab3 = 'loc'
fntemp = fnloc
iprt = iploc
write(lulog,218) lab3,fntemp,iprt
c
lab3 = 'rza'
fntemp = fnrza
iprt = iprza
write(lulog,218) lab3,fntemp,iprt
c
lab3 = 'xya'
fntemp = fnxya
iprt = ipxya
write(lulog,218) lab3,fntemp,iprt
c
lab3 = 'sta'
fntemp = fnsta
iprt = ipsta
write(lulog,218) lab3,fntemp,iprt
c
lab3 = 'fix'
fntemp = fnfix
iprt = ipfix
write(lulog,218) lab3,fntemp,iprt
c
lab3 = 'afx'
fntemp = fnafx
iprt = ipfix
write(lulog,218) lab3,fntemp,iprt
c
lab3 = 'apk'
fntemp = fnapk
iprt = ipapk
write(lulog,218) lab3,fntemp,iprt
c
lab3 = 'npk'
fntemp = fnnpk
iprt = ipnpk
write(lulog,218) lab3,fntemp,iprt
218 format(a3,' data output file: ',a25,' write flag=',i1)
c
c ++ Open the output files, if necessary
c Note: Fix file is opened later if needed, depending
c on the outcome of the TC retrieval
c
if (iploc .eq. 1) then
fntemp=fnloc
open(unit=luloc,file=fnloc,form='formatted',status='unknown',
+ err=900)
endif
c
if (iprza .eq. 1 .and. iaxsym .eq. 1) then
fntemp = fnrza
open(file=fnrza,unit=lurza,form='formatted',status='unknown',
+ err=900)
endif
c
if (ipxya .eq. 1) then
fntemp=fnxya
open(unit=luxya,file=fnxya,form='formatted',status='unknown',
+ err=900)
endif
c
if (ipsta .eq. 1 .and. iaxsym .eq. 1) then
fntemp = fnsta
open(file=fnsta,unit=lusta,form='formatted',status='unknown',
+ err=900)
endif
c
if (ipapk .eq. 1) then
fntemp=fnapk
open(unit=luapk,file=fnapk,form='formatted',status='unknown',
+ err=900)
endif
c
if (ipnpk .eq. 1) then
fntemp=fnnpk
open(unit=lunpk,file=fnnpk,form='formatted',status='unknown',
+ err=900)
endif
c
c ++ Write NCEP variables at AMSU analysis grid points to packed output file
if (ipnpk .eq. 1) then
call wpof(un,vn,tn,zn,psn,iyr2,imon,iday,itime,lunpk)
close(lunpk)
endif
c
c ++ Write amsu data locations to a file
if (iploc .eq. 1) then
call wloc(sslat,sslon,aslat,aslon,luloc,nxas,
+ jyr,jmon,jday,jtimeh,jtimem,atcfid,sname)
close(luloc)
endif
c
c **** End Output File Naming/Opening ****
c
c **** Begin Processing of AMSU temperature/CLW data ****
if (idex .eq. 1) then
c Eliminate data outside analysis domain
call delim(tas,clwas,tpwas,aslat,aslon,
+ rlatd(1),rlatd(ny),rlatbf,
+ rlond(1),rlond(nx),rlonbf,
+ dumas,dumasp,idumas,
+ mxas,mpas,nxas,np,npst,lulog)
endif
c
c ++ Quality control the AMSU temperatures
call qualcon(tas,clwas,tpwas,aslat,aslon,pamsu,dumas,
+ dumasp,idumas,mxas,mpas,nxas,np,npst,lulog)
c
c Write AMSU swath info to the log file
write(lulog,238) nxas,aslat( 1),aslon( 1),
+ aslat(nxas),aslon(nxas)
238 format(/,' Final AMSU swath data includes ',i4,' points',
+ /,' First lat/lon: ',f5.1,1x,f6.1,/,
+ ' Last lat/lon: ',f5.1,1x,f6.1,/)
c
c ++ Correct swath temperatures for cold anomalies if necessary
if (itcor .eq. 3 .or. itcor .eq. 4 .or. itcor .eq. 6 .or.
& itcor .eq. 7) then
write(lulog,*) ' itcor=3 correction applied'
do 23 k=ncs3,nce3
call tcorsr(mxas,nxas,tas(1,k),clwas)
23 continue
endif
c
if (itcor .eq. 5 .or. itcor .eq. 6 .or. itcor .eq. 7) then
write(lulog,*) ' itcor=5 correction applied'
call tcorsv(mxas,nxas,mpas,tas,tn1000,pamsu,
+ clwas,aslon,aslat,clwth1,clwth2,ncs5,nce5,lulog)
endif
c
if (itcor .eq. 2 .or. itcor .eq. 4 .or. itcor .eq. 7) then
write(lulog,*) ' itcor=2 correction applied'
do 22 k=ncs2,nce2
call tcors(mxas,nxas,tas(1,k),clwas,aslon,aslat,
+ sslon,sslat,works,dt,dtred,tcrmax)
22 continue
endif
c
if (itcor .eq. 8 .or. itcor .eq. 10) then
write(lulog,*) ' itcor=8 correction applied'
do 26 k=ncsm8,ncem8
amkl=((amkx*xlambda)/((pamsu(ncem8)-pamsu(ncsm8))**2))*
& (pamsu(k)-pamsu(ncsm8))**2
write(lulog,460) pamsu(k),amkl
460 format(' itcor=8 correction for p=',f6.1,1x,' amkl=',e11.4)
call tcorclw(mxas,nxas,amkl,tas(1,k),clwas,k,lulog)
26 continue
endif
c
c ++ Interpolate temperatures to analysis grid
do 20 k=1,np
do 25 i=1,nxas
dumas(i) = tas(i,npst+k-1)
25 continue
c
call barxy(dumas,aslat,aslon,nxas,efldxy,rinfxy,ispf,
+ rlatd,rlond,t(1,1,k),nx,ny,ierr)
if (ierr .ne. 0) then
write(lulog,940)
940 format(/,'Error in barxy')
stop
endif
20 continue
c
c Interpolate cloud liquid water to analysis grid
do i=1,nxas
if (clwas(i) .gt. 0.0) then
dumas(i) = clwas(i)
else
dumas(i) = 0.0
endif
enddo
c
call barxy(dumas,aslat,aslon,nxas,efldxy,rinfxy,ispf,
+ rlatd,rlond,clwxy,nx,ny,ierr)
if (ierr .ne. 0) then
write(lulog,940)
stop
endif
c
c Correct analysis temperatures for cold anomalies if necessary
if (itcor .eq. 1) then
write(lulog,*) ' itcor=1 correction applied'
do 21 k=ncsm2,ncem2
call tcor(nx,ny,t(1,1,k),dt,dtred)
21 continue
endif
c
c If ice correction doesn't converge, the corresponding pressure
c level, storm, mon, day, and time will be output to the log file
if (itcor .eq. 9 .or. itcor .eq. 10) then
write(lulog,*) ' itcor=9 correction applied'
do 27 k=ncsm9,ncem9
write(lulog,470) p(k)/100.0
470 format(' itcor=9 correction for p=',f6.1)
call tcorice(k,nx,ny,t(1,1,k),clwxy,
+ work1,work2,work3,lulog,lulog,
+ atcfid,imon,iday,itime)
27 continue
write(lulog,*) 'Finished ice corrections.'
endif
c
itest=0
if (itest .eq. 1) then
write (lulog,*) ' itest=1, so program is stopping!'
stop
endif
c
c **** End Processing of AMSU temperature/CLW data ****
c
c **** Begin hydrostatic integration and wind retrieval on 3-D grid ****
c
c Specify surface temperature and pressure of AMSU analysis
c from NCEP analysis
c (All pressures except lateral boundary values will be recomputed)
do 30 j=1,ny
do 30 i=1,nx
ts(i,j) = tsn(i,j)
ps(i,j) = psn(i,j)
30 continue
c
c Check for special case with top AMSU and NCEP levels of 100 mb.
c In this case, the NCEP height field can be used as an upper
c boundary condition.
ispec=0
if (p(1) .eq. 10000.0 .and. pn(1) .eq. 10000.0) then
ispec=1
endif
c
if (intop .eq. 1 .and. ispec .ne. 1) then
write(lulog,920)
920 format(/,' AMSU domain top must be at 100 mb for intop=1')
stop
endif
c
if (intop .eq. 1) then
c Use NCEP height field for upper boundary condition
do 31 j=1,ny
do 31 i=1,nx
z(i,j,1) = zn(i,j,1)
31 continue
c
c Calculate height field at AMSU levels and surface pressure
c for special case
call zalcals(t,ts,ps,p,z,lulog)
else
c Calculate height field at AMSU levels and surface pressure
call zalcal(t,ts,ps,p,z,lulog)
endif
c
c Calculate surface density
do 35 j=1,ny
do 35 i=1,nx
rhos(i,j) = ps(i,j)/(rd*ts(i,j))
35 continue
c
c Extract AMSU t,z at NCEP levels for output
call altonl(z,t,p,za,ta,pn,ts,ps,nx,ny,np,npn,lulog)
c
c Make first guess for AMSU retrieved winds at NCEP pressure levels.
c Use the lowest available NCEP pressure level for the surface winds.
do 40 k=1,npn
do 40 j=1,ny
do 40 i=1,nx
ua(i,j,k) = un(i,j,k)
va(i,j,k) = vn(i,j,k)
40 continue
c
do 45 j=1,ny
do 45 i=1,nx
us(i,j) = un(i,j,npn)
vs(i,j) = vn(i,j,npn)
45 continue
c
if (iwinda .eq. 1 .or. (iwinda .gt. 1 .and. ifgnbe .eq. 1)) then
c Calculate winds from heights using linear balance equation
do k=1,npn
do j=1,ny
do i=1,nx
work1(i,j) = g*za(i,j,k)
work2(i,j) = ua(i,j,k)
work3(i,j) = va(i,j,k)
enddo
enddo
c
if (ibwin(k) .eq. 1) then
call lbe(work1,work2,work3,nx,ny)
endif
c
if (iwinda .eq. 1) then
do j=1,ny
do i=1,nx
ua(i,j,k) = work2(i,j)
va(i,j,k) = work3(i,j)
enddo
enddo
else
do j=2,ny-1
do i=2,nx-1
ua(i,j,k) = work2(i,j)
va(i,j,k) = work3(i,j)
enddo
enddo
endif
enddo
endif
c
if (iwinda .gt. 1) then
c
do 46 k=1,npn
if (ibwin(k) .ne. 1) go to 46
c
do 47 j=1,ny
do 47 i=1,nx
work1(i,j) = g*za(i,j,k)
work2(i,j) = ua(i,j,k)
work3(i,j) = va(i,j,k)
47 continue
c
if (ifgnbe .eq. 2) then
c Calculate zero curvature first guess
do j=1,ny
worky(j) = 0.0
enddo
c
do j=1,ny
do i=1,nx
work4(i,j) = 0.0
work5(i,j) = ua(i,j,k)
work6(i,j) = va(i,j,k)
enddo
enddo
c
call zinter(work5,0.0,nx,ny)
call zinter(work6,0.0,nx,ny)
call psonxy(work4,worky,work2,work5,ierr)
call psonxy(work4,worky,work3,work6,ierr)
endif
c
if (iwinda .eq. 2) then
call nbei(work1,work2,work3,nx,ny)
elseif (iwinda .eq. 3) then
write(6,*) 'call nbev for p=',pn(k)
inon=1
c call nbei(work1,work2,work3,nx,ny)
call nbev(work1,work2,work3,inon,nx,ny)
elseif (iwinda .eq. 4) then
write(6,*) 'call nbevs for p=',pn(k)
inon=1
call nbevs(work1,work2,work3,inon,nx,ny)
endif
c
do 48 j=1,ny
do 48 i=1,nx
ua(i,j,k) = work2(i,j)
va(i,j,k) = work3(i,j)
48 continue
46 continue
endif
c
itest=0
if (itest .eq. 1) then
write(lulog,879)
879 format(/,' u,v near domain center')
do j=23,19,-1
write(lulog,880) (ua(i,j,8),i=19,23),(va(i,j,8),i=19,23)
880 format(1x,5(f7.1),2x,5(f7.1))
enddo
endif
c
c Calculate spatial averages of ta, tn
label = ' AMSU T, NCEP T'
call spata(ta,tn,pn,nx,ny,npn,ibrem,lulog,label)
c
c Calculate spatial averages of za, zn
label = ' AMSU Z, NCEP Z'
call spata(za,zn,pn,nx,ny,npn,ibrem,lulog,label)
c
c Calculate spatial averages of ps,psn
label = ' AMSU PS, NCEP PS'
call spata1(ps,psn,nx,ny,ibrem,lulog,label)
c
c Calculate spatial averages of tas
label = ' AMSU swath T'
call swata(tas,pamsu,mxas,mpas,nxas,lulog,label)
c
if (ipapk .eq. 1) then
c Write AMSU fields to packed file
call wpof(ua,va,ta,za,ps,iyr2,imon,iday,itime,luapk)
close(luapk)
endif
c
if (ipxya .eq. 1) then
c Write fields (u,v,t,z,ps,clw) to ASCII file
c
c Main header
c write(luxya,176) atcfid,slat00,slon00,iyr,imon,iday,
c + itime,sname,
c + sslat,sslon,jyr,jmon,jday,jtimeh,jtimem
c 176 format(' Storm ID: ',a6,f6.2,1x,f7.2,
c + ' at ',i4,1x,2(i2.2),1x,i2.2,1x,a10,/,
c + ' Swath information: ',6x,f6.2,1x,f7.2,
c + ' at ',i4,1x,2(i2.2),1x,i2.2,i2.2)
write(luxya,176) coord
write(luxya,176) times
176 format(a90)
c
do 77 i=1,npn
c Write pressure level to file
write(luxya,78) pn(i),nx,ny
78 format('Pressure level (Pa)=',f8.1,1x,
+ 'nlat=',i2,1x,'nlon=',i2)
c
c Write header line to file
write(luxya,79)
79 format(4x,'Lat',5x,'Lon',11x,'U',11x,'V',11x,'T',11x,'Z')
c
do 177 j=1,ny
do 177 k=1,nx
write(luxya,178) rlatd(j),rlond(k),ua(k,j,i),va(k,j,i),
+ ta(k,j,i),za(k,j,i)
178 format(f7.2,1x,f7.2,4(1x,f11.2))
177 continue
77 continue
c
c Write information about surface variables
write(luxya,80) ny,ny
80 format('Surface pressure level',1x,'nlat=',i2,1x,'nlon=',i2)
c
c Write header line for surface pressure level
write(luxya,81)
81 format(4x,'Lat',5x,'Lon',11x,'U',11x,'V',11x,'T',11x,'P',
+ 10x,'CLW')
c
do 179 j=1,ny
do 179 k=1,nx
write(luxya,180) rlatd(j),rlond(k),us(k,j),vs(k,j),
+ ts(k,j),ps(k,j),clwxy(k,j)
180 format(f7.2,1x,f7.2,5(1x,f11.2))
179 continue
c
close(luxya)
endif
c
c Find location of minimum surface pressure in x,y analysis
pminxy = 1.0e+10
ibuf=10
jbuf=10
im = 0
jm = 0
do i=ibuf,nx-ibuf
do j=jbuf,ny-jbuf
if (ps(i,j) .lt. pminxy) then
rlatpm = rlatd(j)
rlonpm = rlond(i)
im = i
jm = j
pminxy = ps(i,j)
endif
enddo
enddo
pminxy = pminxy/100.0
c
c Refine the min pressure location
pmaxl = -1.0e+10
do i=-1,1
do j=-1,1
ii = im + i
jj = jm + j
if (ps(ii,jj) .gt. pmaxl) pmaxl = ps(ii,jj)
enddo
enddo
c
wtsum = 0.0
rlatpma = 0.0
rlonpma = 0.0
do i=-1,1
do j=-1,1
ii = im + i
jj = jm + j
wt = (pmaxl-ps(ii,jj))
wtsum = wtsum + wt
rlatpma = rlatpma + wt*rlatd(jj)
rlonpma = rlonpma + wt*rlond(ii)
enddo
enddo
c
if (wtsum .gt. 0.0) then
rlatpma = rlatpma/wtsum
rlonpma = rlonpma/wtsum
else
rlatpma = 0.0
rlonpma = 0.0
endif
c
c **** End hydrostatic integration and wind retrieval on 3-D grid ****
c
if (iaxsym .ne. 1) go to 5000
c
c **** Begin axisymmetric retrievals ****
c
c sslat = rlatpma
c sslon = rlonpma
c
c Interpolate AMSU swath temperatures to radial grid
if (iradxy .eq. 1) then
kk=0
do 57 j=1,ny
do 57 i=1,nx
kk = kk + 1
temlat(kk) = rlatd(j)
temlon(kk) = rlond(i)
57 continue
c
do 58 k=1,np
kk=0
do 59 j=1,ny
do 59 i=1,nx
kk = kk + 1
ftem(kk) = t(i,j,k)
59 continue
call barr(ftem,temlat,temlon,kk,efldr,exfac,rinfr,
+ sslat,sslon,trp(1,k),rr,nr,ierr)
58 continue
c
else
do 50 k=1,np
do 51 i=1,nxas
dumas(i) = tas(i,npst+k-1)
51 continue
c
call barr(dumas,aslat,aslon,nxas,efldr,exfac,rinfr,
+ sslat,sslon,trp(1,k),rr,nr,ierr)
c
if (ierr .ne. 0) then
write(lulog,950)
950 format(/,'Error in barr')
stop
endif
50 continue
endif
c
c Interpolate surface temperature and pressure to radial grid
c for first guess
kk=0
do 55 j=1,ny
do 55 i=1,nx
kk = kk + 1
ftem(kk) = tsn(i,j)
temlat(kk) = rlatd(j)
temlon(kk) = rlond(i)
55 continue
c
call barr(ftem,temlat,temlon,kk,efldr,exfac,rinfr,
+ sslat,sslon,tsr,rr,nr,ierr)
c
if (ierr .ne. 0) then
write(lulog,950)
stop
endif
c
kk=0
do 56 j=1,ny
do 56 i=1,nx
kk = kk + 1
ftem(kk) = psn(i,j)
temlat(kk) = rlatd(j)
temlon(kk) = rlond(i)
56 continue
c
call barr(ftem,temlat,temlon,kk,efldr,exfac,rinfr,
+ sslat,sslon,psr,rr,nr,ierr)
c
if (ierr .ne. 0) then
write(lulog,950)
stop
endif
c
c Calculate gradient wind
call bwgcal(rr,p,trp,psr,tsr,dz,nr,np,nz,sslat,
+ zrp,vrp,zz,trz,prz,rhorz,vrz)
c
c Print gradient winds to log file
write(lulog,284)
284 format(/,'Tangential wind profiles')
c
do 90 k=1,nr
write(lulog,285) rr(k)/1000.,prz(k,1)/100.,
+ zz( 1)/1000.,vrz(k, 1),
+ zz( 6)/1000.,vrz(k, 6),
+ zz(11)/1000.,vrz(k,11)
285 format(' r=',f5.0,1x,' ps=',f6.1,3(' z=',f5.1,' v=',f5.1))
90 continue
c
c Calculate temperature anomalies
do k=1,nz
do i=1,nr
tarz(i,k) = trz(i,k)-trz(nr,k)
enddo
enddo
c
if (iprza .eq. 1 .and. iaxsym .eq. 1) then
c Write all r,z variables to a file
labelr = ' '
write(labelr,701) atcfid,imon,iday,itime,
+ jmon,jday,jtimeh,jtimem
701 format(a6,1x,3(i2.2),1x,4(i2.2))
call tcrwrit(prz,rhorz,trz,vrz,nr,nz,sslat,sslon,lurza,
+ labelr,coord,times,rr(1),rr(nr),zz(1),zz(nz))
endif
c
if (ipsta .eq. 1 .and. iaxsym .eq. 1) then
c Write basic statistics of the analysis to a file
write(lusta,700) atcfid,slat00,slon00,iyr,imon,iday,
+ itime,sname,
+ sslat,sslon,jyr,jmon,jday,jtimeh,jtimem
700 format(' Analysis statistics for ',a6,f6.2,1x,f7.2,
+ ' at ',i4,1x,2(i2.2),1x,i2.2,1x,a10,/,
+ ' Swath information: ',6x,f6.2,1x,f7.2,
+ ' at ',i4,1x,2(i2.2),1x,i2.2,i2.2)
c
c Evaluate parameters from analyses
c
c Find minimum pressure and pressure drop
pmin = prz(1,1)/100.0
delp = (prz(nr,1)-prz(1,1))/100.0
c
c Find pressure drop at 3 km
delp3 = (prz(nr,4)-prz(1,4))/100.0
c
c Find max temperature anomalies at 3 km or above
tmax = -999.0
ztmax= -99.0
do k=4,nz
if (tarz(1,k) .gt. tmax) then
tmax = tarz(1,k)
ztmax= zz(k)/1000.0
endif
enddo
c
c Find max winds (first local maximum starting at r=0)
thlm = 5.0
vmx0 = -99.
vmx3 = -99.
rmx0 = -99.
rmx3 = -99.
c
isz0 = 1
isz3 = 1
do i=1,nr-1
dvmax0 = vmx0-vrz(i,1)
if (isz0 .eq. 1) then
if (dvmax0 .gt. thlm) isz0=0
endif
c
if (vrz(i,1) .gt. vmx0 .and. isz0 .eq. 1) then
vmx0 = vrz(i,1)
rmx0 = rr(i)/1000.0
endif
c
dvmax3 = vmx0-vrz(i,4)
if (isz3 .eq. 1) then
if (dvmax3 .gt. thlm) isz3=0
endif
c
if (vrz(i,4) .gt. vmx3 .and. isz3 .eq. 1) then
vmx3 = vrz(i,4)
rmx3 = rr(i)/1000.0
endif
enddo
vmx0 = 1.944*vmx0
vmx3 = 1.944*vmx3
c
c Find wind radii
v15 = 15.0/1.944
v25 = 25.0/1.944
v35 = 35.0/1.944
c
r035 = -99.0
r025 = -99.0
r015 = -99.0
c
r335 = -99.0
r325 = -99.0
r315 = -99.0
c
if (vmx0 .gt. 35.0) then
rmx = 1000.0*rmx0
do i=1,nr
if (vrz(i,1) .lt. v35 .and. rr(i) .gt. rmx) then
r035 = rr(i)/1000.0
go to 1001
endif
enddo
1001 endif
c
if (vmx0 .gt. 25.0) then
rmx = 1000.0*rmx0
do i=1,nr
if (vrz(i,1) .lt. v25 .and. rr(i) .gt. rmx) then
r025 = rr(i)/1000.0
go to 1002
endif
enddo
1002 endif
c
if (vmx0 .gt. 15.0) then
rmx = 1000.0*rmx0
do i=1,nr
if (vrz(i,1) .lt. v15 .and. rr(i) .gt. rmx) then
r015 = rr(i)/1000.0
go to 1003
endif
enddo
r015 = rr(nr)/1000.0
1003 endif
c
if (vmx3 .gt. 35.0) then
rmx = 1000.0*rmx3
do i=1,nr
if (vrz(i,4) .lt. v35 .and. rr(i) .gt. rmx) then
r335 = rr(i)/1000.0
go to 4001
endif
enddo
4001 endif
c
if (vmx3 .gt. 25.0) then
rmx = 1000.0*rmx0
do i=1,nr
if (vrz(i,4) .lt. v25 .and. rr(i) .gt. rmx) then
r325 = rr(i)/1000.0
go to 4002
endif
enddo
4002 endif
c
if (vmx3 .gt. 15.0) then
rmx = 1000.0*rmx0
do i=1,nr
if (vrz(i,4) .lt. v15 .and. rr(i) .gt. rmx) then
r315 = rr(i)/1000.0
go to 4003
endif
enddo
r315 = rr(nr)/1000.0
4003 endif
c
c Find mean tangential wind at 0,3 and 5 km
c for r=0 to 250 km and r=250 to 500 km
rthi = 250.0e+3
rtho = 500.0e+3
c
vbi0 = 0.0
vbi3 = 0.0
vbi5 = 0.0
vbo0 = 0.0
vbo3 = 0.0
vbo5 = 0.0
c
cbi = 0.0
cbo = 0.0
c
do i=1,nr
if (rr(i) .le. rthi) then
cbi = cbi + 1.0
vbi0 = vbi0 + vrz(i,1)
vbi3 = vbi3 + vrz(i,4)
vbi5 = vbi5 + vrz(i,6)
endif
c
if (rr(i) .gt. rthi .and. rr(i) .le. rtho) then
cbo = cbo + 1.0
vbo0 = vbo0 + vrz(i,1)
vbo3 = vbo3 + vrz(i,4)
vbo5 = vbo5 + vrz(i,6)
endif
enddo
c
vbi0 = vbi0/cbi
vbi3 = vbi3/cbi
vbi5 = vbi5/cbi
c
vbo0 = vbo0/cbo
vbo3 = vbo3/cbo
vbo5 = vbo5/cbo
c
c Average cloud liquid water near storm center
c and find fractional area within radat km covered
c by cthresh mm of CLW
clwavg = 0.0
cnt = 0.0
radmin = 100.0
c
radat = 300.0
cthresh = 0.5
pcount = 0.0
pthresh = 0.0
c
do j=1,ny
do i=1,nx
call distk(sslon,sslat,rlond(i),rlatd(j),dx,dy,rad)
c
if (rad .lt. radmin .and. clwxy(i,j) .ge. 0.0) then
clwavg = clwavg + clwxy(i,j)
cnt = cnt + 1.0
endif
c
if (rad .lt. radat) then
pcount = pcount + 1.0
if (clwxy(i,j) .gt. cthresh) then
pthresh = pthresh + 1.0
endif
endif
enddo
enddo
c
if (cnt .gt. 0.0) then
clwavg = clwavg/cnt
else
clwavg = 0.0
endif
c
if (pcount .gt. 0.0) then
pcrat = 100.0*(pthresh/pcount)
else
pcrat = 0.0
endif
c
write(lusta,710) pmin,delp,delp3,tmax,ztmax,ares,
+ vmx0,rmx0,vmx3,rmx3
710 format(' Min. Pressure= ',f6.1,/,
+ ' Sfc. P drop = ',f6.1,/,
+ ' 3 km P drop = ',f6.1,/,
+ ' Max T anomaly= ',f6.1,/,
+ ' z(max T) = ',f6.1,/,
+ ' Swath spacing= ',f6.1,/,
+ ' Vmx(z=0) = ',f6.1,/,
+ ' r(vmx0) = ',f6.1,/,
+ ' Vmx(z=3) = ',f6.1,/,
+ ' r(vmx3) = ',f6.1)
c
write(lusta,720) vbi0,vbi3,vbi5,vbo0,vbo3,vbo5,clwavg,
+ pcrat
720 format(' vbi0 = ',f6.1,/,
+ ' vbi3 = ',f6.1,/,
+ ' vbi5 = ',f6.1,/,
+ ' vbo0 = ',f6.1,/,
+ ' vbo3 = ',f6.1,/,
+ ' vbo5 = ',f6.1,/,
+ ' Average CLW = ',f6.2,/,
+ ' CLW percent = ',f6.2)
c
write(lusta,730) pminxy,rlatpm,rlonpm,rlatpma,rlonpma
730 format(' pminxy = ',f6.1,/,
+ ' lat of pminxy= ',f6.1,/,
+ ' lon of pminxy= ',f6.1,/,
+ ' adj lat = ', f7.2,/,
+ ' adj lon = ', f7.2)
c
if (ipsta .eq. 1) then
c Estimate max wind and azimuthal mean wind radii
c using statistical relationships
c Fill predictor array for new form of vradp routine
preds( 1) = pmin
preds( 2) = delp
preds( 3) = delp3
preds( 4) = tmax
preds( 5) = ztmax
preds( 6) = ares
preds( 7) = vmx0
preds( 8) = rmx0
preds( 9) = vmx3
preds(10) = rmx3
preds(11) = vbi0
preds(12) = vbi3
preds(13) = vbi5
preds(14) = vbo0
preds(15) = vbo3
preds(16) = vbo5
preds(17) = clwavg
preds(18) = pcrat
c
spd = float(ispeed)
head = float(idir)
vmaxop = float(ivmax)
ias = 0
c
call vradp(preds,sslat,sslon,spd,head,vmaxop,ias,
+ pvmx,ppmin,rmp,xp,pr34,pr50,pr64)
c
write(lusta,735) pvmx,ppmin,
+ pr34(5),pr50(5),pr64(5),
+ pr34(6),pr50(6),pr64(6),
+ vmaxop,spd,head,
+ rmp,xp,(pr34(kk),kk=1,4),
+ (pr50(kk),kk=1,4),
+ (pr64(kk),kk=1,4)
735 format(/,' Statistical Intensity/Radii Estimates',
+ /,1x,f5.0,1x,f5.0,7x,' AMSU Max wind (kt), Min P (hPa) ',
+ /,1x,3(f5.0,1x), ' Mean 34,50,64 kt wind radius (nm)',
+ /,1x,3(f5.0,1x), ' Fit 34,50,64 kt wind radius (nm)',
+ /,1x,3(f5.0,1x), ' ATCF max wind, speed, heading',
+ /,1x,f5.0,1x,f5.2,7x,' rm (nm) and x from vortex fit',
+ /,1x,4(f5.0,1x),' 34 kt radii NE,SE,SW,NW',
+ /,1x,4(f5.0,1x),' 50 kt radii NE,SE,SW,NW',
+ /,1x,4(f5.0,1x),' 64 kt radii NE,SE,SW,NW')
endif
endif
c
5000 continue
c
if (ipfix .eq. 1 .and. iaxsym .eq. 1) then
c Write information to NHC fix file
ippmin = ifix(ppmin+0.5)
ipvmax = ifix(pvmx +0.5)
c
if (ifixall .eq. 1 .and. inhce .ne. 0) then
c No analysis was performed, and fix file name not created,
c so create it now
fnfix = 'FX0000_000000.DAT'
fnfix( 2: 2) = atcfid(1:1)
fnfix( 3: 6) = atcfid(3:6)
fnfix(18:19) = cnsat
write(fnfix(8:13),226) imon,iday,itime
c
if (igenfn .eq. 1) then
fnfix=' '
fnfix(1:6) = noaass
fnfix(7:10) = '.FIX'
endif
endif
c
fntemp=fnfix
open(file=fnfix,unit=lufix,form='formatted',status='unknown',
+ err=900)
c See DATELINE ISSUE. - Knaff 2/20/2006
if(sslon .gt. 180.0) sslon = sslon - 360.0
if(slon00 .gt. 180.0) slon00 = slon00 - 360.0
if(slon12 .gt. 180.0) slon12 = slon12 - 360.0
call fixprt(lufix,threshr,radsts,inhce,ippmin,ipvmax,cnsat)
c
if (inhce .eq. 0) then
c Open and write fix file in ATCF format
c
fntemp=fnafx
open(file=fnafx,unit=luafx,form='formatted',
+ status='unknown',err=900)
c
read(atcfid(3:4),295) istnum
295 format(i2)
irmp = ifix(0.5 + rmp)
isdist = ifix(rad)
cstype(1:4)='NOAA'
cstype(5:6)=cnsat
c
do i=1,4
ir34(i) = ifix(0.5 + pr34(i))
ir50(i) = ifix(0.5 + pr50(i))
ir64(i) = ifix(0.5 + pr64(i))
enddo
c
call AMSUfdeck(cbasin,istnum,cstype,
+ jyr,jmon,jday,jtimeh,jtimem,sslat,sslon,
+ ippmin,ipvmax,ir34,ir50,ir64,
+ irmp,isdist,luafx)
c
endif
endif
c
write(lulog,290)
290 format(/,' Program oparet completed normally.')
c
stop
c
c Error processing
900 continue
write(lulog,990) fntemp
990 format(/,' Error opening file ',a25)
c
stop
end
subroutine fixprt(lufix,thresh,rad,ierr,ippmin,ipvmax,cnsat)
+
c This routine writes the fix file
c
dimension pr34(6),pr50(6),pr64(6)
character *2 cnsat
character *3 cdow
character *6 atcfid
character *9 sname
c
dimension ipr34(6),ipr50(6),ipr64(6)
c
common /sinfo1/ atcfid,sname
common /sinfo2/ iyr,imon,iday,itime4,
+ jyr,jmon,jday,jtime4,
+ ivmax,spd,head,
+ slat00,slat12,sslat,
+ slon00,slon12,sslon,
+ rmp,xp,pr34,pr50,pr64
c
c Convert output variables to integers
do k=1,6
ipr34(k) = ifix(0.5 + pr34(k))
ipr50(k) = ifix(0.5 + pr50(k))
ipr64(k) = ifix(0.5 + pr64(k))
enddo
ispd = ifix(0.5 + spd)
ihead= ifix(0.5 + head)
irmp = ifix(0.5 + rmp)
irad = ifix(0.5 + rad)
c
c Get current date and time
call cdt(kyr,kmon,kday,ktime4,cdow)
c
c Calculate time difference between AMSU and ATCF date/times
itime2 = itime4/100
jtime2 = jtime4/100
call tdiff(jyr,jmon,jday,jtime2,iyr,imon,iday,itime2,idtsa)
c
if (ierr .eq. 0) then
write(lufix,200) cnsat
200 format(' ***************************************',
+ '*******************************',
+ /,' CIRA/NESDIS Experimental AMSU-A TC',
+ ' Intensity/Size Estimation - NOAA',a2,/)
c
write(lufix,210) atcfid(1:4),iyr,sname
210 format(' Tropical Cyclone ',1x,a4,i4.4,2x,a9)
c
write(lufix,220) kyr,kmon,kday,ktime4,
+ iyr,imon,iday,itime4,
+ jyr,jmon,jday,jtime4
220 format(' Current date/time: ',i4,1x,2(i2.2),1x,i4.4,' UTC',
+ /,' ATCF file date/time: ',i4,1x,2(i2.2),1x,i4.4,' UTC',
+ /,/,' AMSU swath date/time: ',i4,1x,2(i2.2),1x,i4.4,' UTC')
c
write(lufix,300) ippmin,ipvmax
300 format(/,' Minimum Sea-Level Pressure: ',i5,' hPa',
+ /,' Maximum Surface Winds: ',i5,' kt',/)
c
write(lufix,310) (ipr34(k),k=1,4)
310 format(' 34 kt wind radii (NE,SE,SW,NW): ',4(i4),' nmi')
write(lufix,312) (ipr50(k),k=1,4)
312 format(' 50 kt wind radii (NE,SE,SW,NW): ',4(i4),' nmi')
write(lufix,314) (ipr64(k),k=1,4)
314 format(' 64 kt wind radii (NE,SE,SW,NW): ',4(i4),' nmi')
c
write(lufix,320) irmp
320 format(/,' AMSU-retrieved max wind radius: ',i4,' nmi')
c
write(lufix,330) irad
330 format(/,' Storm center is ',i4,' km from AMSU swath center',
+ /,' 0-300 km is optimal ',
+ /,' 300-600 km is adequate',
+ /,' >600 km is marginal')
c
write(lufix,340) idtsa
340 format(/,' AMSU data is ',i4,' hr from time of ATCF input',/)
c
write(lufix,201)
201 format(' ***************************************',
+ '*******************************')
c
write(lufix,400) atcfid(1:4),iyr,imon,iday,itime4
400 format(' ATCF File Input: ',
+ /,1x,a4,i4.4,1x,2(i2.2),1x,i4.4,' UTC')
c
write(lufix,410) slat00,slon00,slat12,slon12,
+ idtsa,sslat,sslon
410 format(/,' Storm lat,lon (t = 0 hr): ',f6.2,1x,f7.2,
+ /,' Storm lat,lon (t = -12 hr): ',f6.2,1x,f7.2,
+ /,' Storm lat,lon (t =',i4,' hr): ',f6.2,1x,f7.2,
+ ' (AMSU swath time)')
c
write(lufix,420) ivmax,ihead,ispd
420 format(/,' Storm max winds (ATCF): ',i3 ,' kt',
+ /,' Storm heading: ',i3.3,' deg',
+ /,' Storm translation speed: ',i3 ,' kt',
+ //,' Note: AMSU wind radii provided for all wind thresholds',
+ /,' up to the ATCF max winds. Thus, AMSU wind radii ',
+ /,' may be provided for thresholds that exceed the ',
+ /,' AMSU max wind estimate. ')
c
write(lufix,201)
return
elseif (ierr .eq. 1) then
write(lufix,200) cnsat
write(lufix,210) atcfid(1:4),iyr,sname
write(lufix,220) kyr,kmon,kday,ktime4,
+ iyr,imon,iday,itime4,
+ jyr,jmon,jday,jtime4
c
write(lufix,901) rad,thresh
901 format(/,' INTENSITY/SIZE ESTIMATION FAILED',
+ /,' Storm too far from center of AMSU data swath',
+ /,' Observed Distance: ',f5.0,' km',
+ /,' Max allowable distance: ',f5.0,' km')
c
write(lufix,201)
stop
elseif (ierr .eq. 2) then
write(lufix,200) cnsat
write(lufix,210) atcfid(1:4),iyr,sname
write(lufix,220) kyr,kmon,kday,ktime4,
+ iyr,imon,iday,itime4,
+ jyr,jmon,jday,jtime4
c
write(lufix,902)
902 format(/,' INTENSITY/SIZE ESTIMATION FAILED',
+ /,' Error Processing AMSU data',
+ /,' (Usually caused by data that is not available yet)')
c
write(lufix,201)
stop
endif
c
return
end
subroutine cdt(kyr,kmon,kday,ktime,cdow)
c This routine gets the current dat and time (UTC)
c from the file tdate.dat that is created by the mdasmu.ksh
c script using the HP-UX command date -u. The output of this
c command is of the form
c
c Fri May 17 16:00:42 GMT 2002
c
dimension cmonl(12)
c
character *3 cdow,cmon,cmonl
c
kyr = 0
kmon= 0
kday= 0
ktime=0
cdow='XXX'
c
cmonl( 1) = 'Jan'
cmonl( 2) = 'Feb'
cmonl( 3) = 'Mar'
cmonl( 4) = 'Apr'
cmonl( 5) = 'May'
cmonl( 6) = 'Jun'
cmonl( 7) = 'Jul'
cmonl( 8) = 'Aug'
cmonl( 9) = 'Sep'
cmonl(10) = 'Oct'
cmonl(11) = 'Nov'
cmonl(12) = 'Dec'
c
open(unit=65,file='tdate.dat',form='formatted',status='old',
+ err=900)
c
read(65,100,err=900) cdow,cmon,kday,khr,kmin,kyr
100 format(a3,1x,a3,1x,i2,1x,i2,1x,i2,8x,i4)
close(65)
c
do k=1,12
if (cmon .eq. cmonl(k)) kmon = k
enddo
if (kmon .eq. 0) go to 900
c
ktime = 100*khr + kmin
c
return
c
900 continue
c
kyr = 0
kmon= 0
kday= 0
ktime=0
cdow='XXX'
c
return
end