subroutine zalcals(t,ts,ps,p,z,lulog)
c This routine calculates the height z at all AMSU levels,
c and the surface pressure ps at the domain interior points.
c
c This version assumes z at the top level has been specified
c and is used as an upper boundary condition
c
parameter (nx=61,ny=61)
parameter (np=23)
c
dimension t(nx,ny,np),z(nx,ny,np)
dimension ts(nx,ny),ps(nx,ny)
dimension p(np)
c
c Integrate z down to lowest AMSU level
do 25 j=1,ny
do 25 i=1,nx
do 27 k=2,np
t1 = t(i,j,k-1)
t2 = t(i,j,k)
p1 = p(k-1)
p2 = p(k)
call tkness(p1,p2,t1,t2,delz)
z(i,j,k) = z(i,j,k-1) - delz
27 continue
25 continue
c
c Calculate surface pressure by integrating from lowest
c AMSU level to the surface
do 30 j=1,ny
do 30 i=1,nx
t1 = t(i,j,np)
z1 = z(i,j,np)
p1 = p(np)
t2 = ts(i,j)
z2 = 0.0
call p2cal(z1,z2,t1,t2,p1,p2)
ps(i,j) = p2
30 continue
c
return
end
subroutine zalcal(t,ts,ps,p,z,lulog)
c This routine calculates the height z at all AMSU levels,
c and the surface pressure ps at the domain interior points.
c
parameter (nx=61,ny=61)
parameter (np=23)
c
dimension t(nx,ny,np),z(nx,ny,np)
dimension ts(nx,ny),ps(nx,ny)
dimension p(np)
c
c Local work arrays
dimension ta1(nx,ny),ta2(nx,ny)
c
c Calculate height of lowest AMSU level at boundary points
do 10 i=1,nx
c Bottom edge of domain
j=1
t1 = ts(i,j)
t2 = t(i,j,np)
p1 = ps(i,j)
p2 = p(np)
call tkness(p1,p2,t1,t2,delz)
z(i,j,np) = delz
c
c Top edge of domain
j=ny
t1 = ts(i,j)
t2 = t(i,j,np)
p1 = ps(i,j)
p2 = p(np)
call tkness(p1,p2,t1,t2,delz)
z(i,j,np) = delz
10 continue
c
do 12 j=2,ny-1
c Left edge of domain
i=1
t1 = ts(i,j)
t2 = t(i,j,np)
p1 = ps(i,j)
p2 = p(np)
call tkness(p1,p2,t1,t2,delz)
z(i,j,np) = delz
c
c Right edge of domain
i=nx
t1 = ts(i,j)
t2 = t(i,j,np)
p1 = ps(i,j)
p2 = p(np)
call tkness(p1,p2,t1,t2,delz)
z(i,j,np) = delz
12 continue
c
c Calculate z at the rest of the AMSU levels at the boundary points
do 15 k=np-1,1,-1
do 17 i=1,nx
c Bottom edge of domain
j=1
t1 = t(i,j,k+1)
t2 = t(i,j,k )
p1 = p(k+1)
p2 = p(k)
call tkness(p1,p2,t1,t2,delz)
z(i,j,k) = z(i,j,k+1) + delz
c
c Top edge of domain
j=ny
t1 = t(i,j,k+1)
t2 = t(i,j,k )
p1 = p(k+1)
p2 = p(k)
call tkness(p1,p2,t1,t2,delz)
z(i,j,k) = z(i,j,k+1) + delz
if (k .eq. 1) then
endif
17 continue
c
do 19 j=2,ny-1
c Left edge of domain
i=1
t1 = t(i,j,k+1)
t2 = t(i,j,k )
p1 = p(k+1)
p2 = p(k)
call tkness(p1,p2,t1,t2,delz)
z(i,j,k) = z(i,j,k+1) + delz
c
c Right edge of domain
i=nx
t1 = t(i,j,k+1)
t2 = t(i,j,k )
p1 = p(k+1)
p2 = p(k)
call tkness(p1,p2,t1,t2,delz)
z(i,j,k) = z(i,j,k+1) + delz
if (k .eq. 1) then
endif
19 continue
15 continue
c
c Interpolate z (top pressure level) at the boundary points
c to the domain interior using a Laplacian filter
do 20 j=1,ny
do 20 i=1,nx
ta1(i,j) = 0.0
ta2(i,j) = z(i,j,1)
20 continue
c
call pson(ta1,ta2,z(1,1,1),ierr)
if (ierr .ne. 0) then
write(lulog,960)
960 format(/,' pson routine did not converge in z calculation')
stop
endif
c
c Integrate z down to lowest AMSU level at interior points
do 25 j=2,ny-1
do 25 i=2,nx-1
do 27 k=2,np
t1 = t(i,j,k-1)
t2 = t(i,j,k)
p1 = p(k-1)
p2 = p(k)
call tkness(p1,p2,t1,t2,delz)
z(i,j,k) = z(i,j,k-1) - delz
27 continue
25 continue
c
c Calculate surface pressure by integrating from lowest
c AMSU level to the surface
do 30 j=2,ny-1
do 30 i=2,nx-1
t1 = t(i,j,np)
z1 = z(i,j,np)
p1 = p(np)
t2 = ts(i,j)
z2 = 0.0
call p2cal(z1,z2,t1,t2,p1,p2)
ps(i,j) = p2
30 continue
c
return
end
subroutine ncepget(pn,un,vn,zn,tn,tn1000,rlond,rlatd,nx,ny,npn,
+ dlona,dlata,dayx,utcx,lundat,lulog,ierr)
c This routine gets u,v,z and t from the packed ncep analysis file
c
parameter (ixmax=181,iymax=91)
c
dimension pn(npn)
dimension un(nx,ny,npn),vn(nx,ny,npn),zn(nx,ny,npn),tn(nx,ny,npn)
dimension rlond(nx),rlatd(ny)
c
dimension tn1000(ixmax,iymax)
c
c Local variables
c
c Data for unpacking
parameter(imax=20000)
dimension tra(imax)
character *1 type
character *2 code(imax)
c
dimension t2d(ixmax,iymax)
c
c Array for checking that all requested NCEP levels were found
parameter (nvar=4)
dimension ncheck(100,nvar)
character *1 ccheck(nvar)
c
data ccheck /'U','V','Z','T'/
c
common /ncepfg/ rlonln,rlatbn,dlonn,dlatn,nlonn,nlatn
c
c Set interp=1 to allow option of interpolating between
c levels of NCEP data. There must be data above and below
c the levels requiring interpolation.
interp=1
c
c Set checking array to zero
do 5 k=1,npn
do 5 i=1,nvar
ncheck(k,i) = 0
5 continue
c
lu = lundat
linen = 0
c
c Read main header on packed data file
read(lu,200,err=901,end=901)
+ wx,dayx,utcx,rlatmn,rlatmx,rlonmn,rlonmx,dlatn,dlonn
200 format(1x,f3.0,f7.0,f5.0,4f8.3,1x,2f4.2)
linen = linen + 1
c
c Calculate lat and lon intervals and number of points on NCEP
c data file
epsil = 0.0001
nlat1 = 1 + ifix((rlatmx-rlatmn)/dlatn + epsil)
nlon1 = 1 + ifix((rlonmx-rlonmn)/dlonn + epsil)
ipts = nlat1*nlon1
rlatb = rlatmn
rlatt = rlatmx
rlonl = -rlonmx
rlonr = -rlonmn
c
rlonln = rlonl
rlatbn = rlatb
nlonn = nlon1
nlatn = nlat1
c
write(lulog,300) rlonl,rlonr,dlonn,rlatb,rlatt,dlatn,
+ nlon1,nlat1,ipts
300 format(/,' Data read from NCEP file',/,
+ ' rlonl=',f6.1,' rlonr=',f6.1,' dlon=',f4.2,/,
+ ' rlatb=',f6.1,' rlatt=',f6.1,' dlat=',f4.2,/,
+ ' nlon1=',i6,' nlat1=',i6,' ipts=',i6)
c
c Check array size
if (ipts .gt. imax) then
ierr = 2
return
endif
c
if (nlon1 .gt. ixmax .or. nlat1 .gt. iymax) then
ierr = 2
return
endif
c
c Read the rest of the data file
nrow = 1 + (ipts-1)/36
c
500 continue
read(lu,202,err=901,end=600) type,ptmp,bsub,smpy
202 format(1x,a1,1x,f6.1,2(1x,g15.9))
linen = linen + 1
c
c write(lulog,320) type,ptmp,linen,bsub,smpy
c 320 format(1x,a1,' p=',f6.1,' read from line ',i5,
c + ' of ncep file, b,s: ',2(e10.3))
c
c Read packed data
do 10 n=1,nrow
is = 1 + (n-1)*36
ie = is + 35
read(lu,204,end=901,err=901) (code(i),i=is,ie)
linen = linen + 1
204 format(36(a2))
10 continue
c
c Search for required NCEP pressure level and data type
ii = 0
do 11 i=1,nvar
if (type .eq. ccheck(i)) then
ii = i
endif
11 continue
c
kk = 0
do 12 k=1,npn
if (ptmp .eq. pn(k)) then
kk = k
endif
12 continue
c
if (kk .eq. 0 .or. ii .eq. 0) go to 500
c
c This variable is need, so update checking array and unpack it.
ncheck(kk,ii) = 1
c
c Unpack current variable
do 15 i=1,ipts
idec = idecod(code(i))
tra(i) = float(idec)*smpy - bsub
15 continue
c
c Add standard height to height perturbation
if (type .eq. 'Z') then
call stndz(ptmp,zstd,tstd,thstd)
do 16 i=1,ipts
tra(i) = tra(i) + zstd
16 continue
endif
c
c Put data in 2-D lon/lat array
do 20 j = 1,nlat1
do 20 i = 1,nlon1
ij = i + (j-1)*nlon1
t2d(i,j) = tra(ij)
20 continue
c
alonl = rlond(1)
if (rlonr .eq. 360.0 .and. alonl .lt. 0.0) then
alonl = alonl + 360.0
endif
c
alatb = rlatd(1)
izp = 0
c
c write(6,990) rlonl,rlatb,dlonn,dlatn,ixmax,iymax,nlon1,nlat1
c 990 format(/,' rlonl=',f6.1,' rlatb=',f6.1,' dlonn=',f6.1,' dlatn=',
c + f6.1,/,' ixmax=',i4,' iymax=',i4,' nlon1=',i4,
c + ' nlat1=',i4)
c
c write(6,991) alonl,alatb,dlona,dlata,nx,ny
c 991 format(' alonl=',f6.1,' alatb=',f6.1,' dlona=',f6.1,' dlata=',
c + f6.1,/,' nx=',i4,' ny=',i4)
c
c Interpolate data from NCEP to AMSU grid points
if (type .eq. 'U') then
call llintp( t2d,rlonl,rlatb,dlonn,dlatn,ixmax,iymax,
+ nlon1,nlat1,
+ un(1,1,kk),alonl,alatb,dlona,dlata, nx, ny,
+ nx, ny,
+ izp,lerr)
else if (type .eq. 'V') then
call llintp( t2d,rlonl,rlatb,dlonn,dlatn,ixmax,iymax,
+ nlon1,nlat1,
+ vn(1,1,kk),alonl,alatb,dlona,dlata, nx, ny,
+ nx, ny,
+ izp,lerr)
else if (type .eq. 'T') then
call llintp( t2d,rlonl,rlatb,dlonn,dlatn,ixmax,iymax,
+ nlon1,nlat1,
+ tn(1,1,kk),alonl,alatb,dlona,dlata, nx, ny,
+ nx, ny,
+ izp,lerr)
c
if (ptmp .eq. 1000.0) then
c Save 1000 mb temperature on for tcorsv routine
do j=1,nlat1
do i=1,nlon1
tn1000(i,j) = t2d(i,j)
enddo
enddo
endif
c
else if (type .eq. 'Z') then
call llintp( t2d,rlonl,rlatb,dlonn,dlatn,ixmax,iymax,
+ nlon1,nlat1,
+ zn(1,1,kk),alonl,alatb,dlona,dlata, nx, ny,
+ nx, ny,
+ izp,lerr)
else
ierr = 5
return
endif
c
go to 500
600 continue
c
if (interp .eq. 1) then
c Check to see if data at missing pressure levels can
c be filled by interpolation. Top and bottom levels can
c not be interpolated.
do k=2,npn-1
do m=1,nvar
if (ncheck(k,m) .eq. 0 .and. ncheck(k+1,m) .eq. 1
+ .and. ncheck(k-1,m) .eq. 1) then
c
ncheck(k,m) = 1
c
pk = pn(k)
pkp = pn(k+1)
pkm = pn(k-1)
c
wtm = (pkp-pk)/(pkp-pkm)
wtp = (pk-pkm)/(pkp-pkm)
c
if (m .eq. 1) then
do j=1,ny
do i=1,nx
un(i,j,k) = wtm*un(i,j,k-1) + wtp*un(i,j,k+1)
enddo
enddo
elseif (m .eq. 2) then
do j=1,ny
do i=1,nx
vn(i,j,k) = wtm*vn(i,j,k-1) + wtp*vn(i,j,k+1)
enddo
enddo
elseif (m .eq. 3) then
do j=1,ny
do i=1,nx
zn(i,j,k) = wtm*zn(i,j,k-1) + wtp*zn(i,j,k+1)
enddo
enddo
elseif (m .eq. 4) then
do j=1,ny
do i=1,nx
tn(i,j,k) = wtm*tn(i,j,k-1) + wtp*tn(i,j,k+1)
enddo
enddo
endif
endif
enddo
enddo
endif
c Check to make sure all levels have been found
do 90 k=1,npn
do 90 i=1,nvar
if (ncheck(k,i) .eq. 0) then
write(lulog,950) ccheck(i),pn(k)
950 format(/,' Required NCEP variable not found, type= ',a1,
+ ' P=',f6.1)
ierr = 4
endif
90 continue
c
write(lulog,330) nvar,npn
330 format(/,'All ',i1,' NCEP variables found at all ',i2,
+ ' requested levels')
c
ierr = 0
return
c
901 continue
ierr = 1
return
c
end
subroutine aprint(a,nx,ny,plev,iscale,inc,lulog,label)
c This routine prints the 2-d array a.
c If iscale=1, the array is scaled to a number between 0 and 100.
c If iscale<0, the array is scaled by abs(iscale)
c
dimension a(nx,ny)
character *20 label
c
slope = 1.0
offset = 0.0
c
c Find max and min of a
aamax = -1.0e+20
aamin = 1.0e+20
do 10 j=1,ny
do 10 i=1,nx
if (a(i,j) .gt. aamax) aamax=a(i,j)
if (a(i,j) .lt. aamin) aamin=a(i,j)
10 continue
c
if (aamax .ne. aamin .and. iscale .eq. 1) then
offset = aamin
slope = 100.0/(aamax-aamin)
endif
c
if (iscale .lt. 0) then
offset = 0.0
slope = 1.0/float(-iscale)
endif
c
write(lulog,200) label,plev/100.0,aamax,aamin,iscale
200 format(/,1x,a20,' p(mb)=',f6.1,' max,min= ',e11.4,1x,e11.4,
+ ' iscale=',i6)
c
istart = 1
iend = nx
c
if (nx .gt. 26) then
icut = 1 + (nx-27)/2
istart = istart + icut
iend = iend - icut
endif
c
if (inc .gt. 0) then
do 20 j=1,ny,inc
write(lulog,210) j,( (slope*(a(i,j)-offset)),i=istart,iend )
210 format(1x,i2,1x,26(f6.1),/,4x,26(f6.1))
20 continue
elseif (inc .lt. 0) then
do 30 j=ny,1,inc
write(lulog,210) j,( (slope*(a(i,j)-offset)),i=istart,iend)
30 continue
endif
c
write(lulog,220) (i,i=istart,iend)
220 format(6x,26(i2,4x),/,6x,26(i2,4x))
c
return
end
subroutine distk(rlon1,rlat1,rlon2,rlat2,dx,dy,rad)
c This routine calculates the distance in km (rad) between the
c points (rlon1,rlat1) and (rlon2,rlat2) using an approximate
c formula. The lon and lat are in deg E and N. The east and
c north components of the distance (dx,dy) are also calculated.
c
common /cons/ pi,g,rd,dtr,erad,erot
c
dtk = 111.1
cfac = cos(0.5*dtr*(rlat1+rlat2))
c
dx = dtk*(rlon2-rlon1)*cfac
dy = dtk*(rlat2-rlat1)
rad = sqrt(dx*dx + dy*dy)
c
return
end
subroutine wpof(u,v,t,z,ps,iyr,imon,iday,itime,luout)
c This routine writes u,v,t,z and ps to a packed output file
c
parameter (nx=61,ny=61)
parameter (npn=12)
dimension u(nx,ny,npn),v(nx,ny,npn),t(nx,ny,npn),z(nx,ny,npn)
dimension ps(nx,ny)
c
c NCEP pressure levels
dimension pn(npn)
c
c Lat/Lon arrays
dimension rlond(nx),rlatd(ny)
dimension rlonr(nx),rlatr(ny)
dimension sinlat(ny),coslat(ny),tanlat(ny)
c
c Local variables
parameter (imax=20000)
dimension tra(imax)
character *1 type
character *2 code(imax)
c
common /log/ lulog
common /ncepll/ rlond,rlonr,rlatd,rlatr,dlon,dlat,
+ dlonr,dlatr,sinlat,coslat,tanlat
common /ncepp/ pn
c
c Write header line on the file
wx = 12.
rdate = float(10000*iyr + 100*imon + iday)
rtime = float(itime)
ft = 0.0
c
write(luout,200) wx,rdate,rtime,rlatd( 1), rlatd(ny),
+ -rlond(nx),-rlond( 1),dlat,dlon,ft
200 format(1x,f3.0,f7.0,f5.0,4f8.3,2f4.2,f6.0)
c
c Calculate number of grid points
ipts = nx*ny
nrow = 1 + (ipts-1)/36
c
if (ipts .gt. imax) then
write(lulog,900)
900 format(/,' AMSU analysis domain too large, increase imax ',
+ /,' in routine wpof. ')
stop
endif
c
c Fill in code array with 1s before packing
do 10 i=1,ipts+36
code(i) = '11'
10 continue
c
epsil = 0.001
c Pack sea-level pressure (in mb)
type = 'S'
c
icount = 0
do 15 j=1,ny
do 15 i=1,nx
icount = icount + 1
tra(icount) = ps(i,j)/100.0
15 continue
c
call maxmin(tra,1,ipts,tmax,tmin)
call tstcod(tra,1,ipts,tmax,tmin,bsub,smpy,code)
c
plevx = 1070.0
write(luout,210) type,plevx,bsub,smpy
210 format(1x,a1,1x,f6.1,2(1x,g15.9))
c
do 20 n=1,nrow
is = 1 + (n-1)*36
ie = is+35
write(luout,220) (code(i),i=is,ie)
220 format(36(a2))
20 continue
c
c Pack the rest of the variables
do 25 k=1,npn
plevx = pn(k)/100.0
call stndz(plevx,zstd,tstd,thstd)
c
type = 'U'
icount = 0
do 31 j=1,ny
do 31 i=1,nx
icount = icount + 1
tra(icount) = u(i,j,k)
31 continue
c
call maxmin(tra,1,ipts,tmax,tmin)
call tstcod(tra,1,ipts,tmax,tmin,bsub,smpy,code)
c
write(luout,210) type,plevx,bsub,smpy
c
do 32 n=1,nrow
is = 1 + (n-1)*36
ie = is+35
write(luout,220) (code(i),i=is,ie)
32 continue
c
type = 'V'
icount = 0
do 41 j=1,ny
do 41 i=1,nx
icount = icount + 1
tra(icount) = v(i,j,k)
41 continue
c
call maxmin(tra,1,ipts,tmax,tmin)
call tstcod(tra,1,ipts,tmax,tmin,bsub,smpy,code)
c
write(luout,210) type,plevx,bsub,smpy
c
do 42 n=1,nrow
is = 1 + (n-1)*36
ie = is+35
write(luout,220) (code(i),i=is,ie)
42 continue
c
type = 'Z'
icount = 0
do 51 j=1,ny
do 51 i=1,nx
icount = icount + 1
tra(icount) = z(i,j,k) - zstd
51 continue
c
call maxmin(tra,1,ipts,tmax,tmin)
call tstcod(tra,1,ipts,tmax,tmin,bsub,smpy,code)
c
write(luout,210) type,plevx,bsub,smpy
c
do 52 n=1,nrow
is = 1 + (n-1)*36
ie = is+35
write(luout,220) (code(i),i=is,ie)
52 continue
c
type = 'T'
icount = 0
do 61 j=1,ny
do 61 i=1,nx
icount = icount + 1
tra(icount) = t(i,j,k)
61 continue
c
call maxmin(tra,1,ipts,tmax,tmin)
call tstcod(tra,1,ipts,tmax,tmin,bsub,smpy,code)
c
write(luout,210) type,plevx,bsub,smpy
c
do 62 n=1,nrow
is = 1 + (n-1)*36
ie = is+35
write(luout,220) (code(i),i=is,ie)
62 continue
25 continue
c
return
end
subroutine tcors(mxas,nxas,t,clw,aslon,aslat,
+ sslon,sslat,works,dt,dtred,tcrmax)
c This routine reduces the temperature anomaly below the
c threshold dt by a factor dtred. This version uses the
c pre-analyzed swath temperature data.
c
dimension t(mxas),clw(mxas)
dimension aslon(mxas),aslat(mxas),works(mxas)
c
c Calculate radius of each swath point
do 5 i=1,nxas
call distk(sslon,sslat,aslon(i),aslat(i),dx,dy,rad)
works(i) = rad
5 continue
c
c Find mean temperature
tbar = 0.0
rpts = 0.0
do 10 i=1,nxas
if (works(i) .gt. tcrmax) go to 10
tbar = tbar + t(i)
rpts = rpts + 1.0
10 continue
c
tbar = tbar/rpts
c
do 20 i=1,nxas
if (works(i) .gt. tcrmax) go to 20
dtt = t(i)-tbar
if (dtt .lt. dt) then
t(i) = tbar + dt + (dtt-dt)*dtred
endif
20 continue
c
return
end
subroutine tcorsv(mxas,nxas,mpas,tas,tn1000,pamsu,
+ clwas,aslon,aslat,clwth1,clwth2,kp1,kp2,lulog)
c This routine adjusts the temperature profile t to a constant
c lapse rate profile in regions where clw > clwth1. For clw > clwth2,
c The original profile is completely replaced by the constant lapse
c rate profile between pressure levels kp1 and kp2. The constant
c lapse rate is determined from the amsu temperature at pressure
c level kp1 and the NCEP analysis temperature at 1000 mb.
c
parameter (ixmax=181,iymax=91)
dimension tas(mxas,mpas),clwas(mxas),aslon(mxas),aslat(mxas)
dimension tn1000(ixmax,iymax)
dimension pamsu(mpas)
c
common /cons/ pi,g,rd,dtr,erad,erot
common /ncepfg/ rlonln,rlatbn,dlonn,dlatn,nlonn,nlatn
c
c Assign pressure levels for constant lapse rate atmosphere
p1 = 1000.0e+2
p2 = 100.0*pamsu(kp1)
c
c Search for points that require adjustment
icnt = 0
do 10 i=1,nxas
if (clwas(i) .gt. clwth1) then
c
c Calculate temperature profile weights
clwt = clwas(i)
if (clwt .gt. clwth2) then
wtam = 0.0
wtcl = 1.0
else
wtam = (clwth2-clwt)/(clwth2-clwth1)
wtcl = (clwt-clwth1)/(clwth2-clwth1)
endif
c
c Calculate parameters for the constant lapse rate
c atmosphere
t2 = tas(i,kp1)
izp = 0
call llintsp(tn1000,rlonln,rlatbn,dlonn,dlatn,ixmax,iymax,
+ nlonn,nlatn,t1,aslon(i),aslat(i),izp,ierr)
gmma = (g/rd)*alog(t2/t1)/alog(p1/p2)
aa = -rd*gmma/g
c
icnt = icnt + 1
write(lulog,300) icnt,aslat(i),aslon(i),clwt,gmma
300 format(i4,' lat,lon: ',f6.2,1x,f7.2,' clw: ',f5.2,
+ ' gmma:',e11.4)
c
do 20 k=kp1,kp2
tcl = t1*(100.0*pamsu(k)/p1)**aa
tas(i,k) = wtam*tas(i,k) + wtcl*tcl
20 continue
endif
10 continue
return
end
subroutine tcorsr(mxas,nxas,t,clw)
c This routine reduces the temperature anomaly by using the
c slope relationship between cloud liquid water and temperature
c anomalies (base temperature comes from regions with no clw).
c This program also adjusts the temperatures over the land areas
c where clw is not available. In this case, if negative temperature
c anomalies exist over land they are replaced with a temperature .2
c degrees colder than the mean + twenty percent of the original
c temperature anomaly.
c This version uses the pre-analyzed swath temperature data.
c
dimension t(mxas),clw(mxas)
c
c Find mean temperature
tbar = 0.0
npts=0
rpts = float(nxas)
do 10 i=1,nxas
if (clw(i) .le. 0.001 .and. clw(i) .gt. -1.0) then
tbar = tbar + t(i)
npts=npts+1
endif
10 continue
c
tbar = tbar/float(npts) ! mean
sxy=0.0 ! Sum of x*y
sxx=0.0 ! Sum of x^2
c
do i=1,nxas
if (clw(i) .gt. 0.0) then
sxx=sxx+clw(i)**2
sxy=sxy+clw(i)*(t(i)-tbar)
endif
enddo
c
slope=sxy/sxx
c write(6,*) 'tbar,npts,slope ',tbar,npts,slope
c
c over land temperature fix
do 20 i=1,nxas
if (slope.lt.0 .and. clw(i) .ge. 0.0) then
t(i) = t(i)-slope*clw(i)
elseif(clw(i).lt.-50.and.(t(i)-tbar).lt.-0.2)then
t(i) = tbar+(-0.2+0.2*(t(i)-tbar))
endif
20 continue
c
return
end
subroutine tcor(nx,ny,t,dt,dtred)
c This routine reduces the temperature anomaly below the
c threshold dt by a factor dtred. This version uses the analyzed
c temperature data.
c
dimension t(nx,ny)
c
c Find mean temperature
tbar = 0.0
rpts = float(nx*ny)
do 10 j=1,ny
do 10 i=1,nx
tbar = tbar + t(i,j)
10 continue
c
tbar = tbar/rpts
c
do 20 j=1,ny
do 20 i=1,nx
dtt = t(i,j)-tbar
if (dtt .lt. dt) then
t(i,j) = tbar + dt + (dtt-dt)*dtred
endif
20 continue
c
return
end
subroutine tcorclw(mxas,nxas,amkl,tas,clwas,k,lulog)
c This routine corrects for CLW attenuation using temp and CLW
c on the swath. Regression slopes at each pressure level
c were derived by comparing tdev vs CLW for 64 noland storms.
c
dimension oldtas(mxas)
dimension tas(mxas)
dimension clwas(mxas)
c
c Read in initial temp values
do 69 i=1,nxas
oldtas(i) = tas(i)
69 continue
c
c Values greater than clwthres are corrected
clwthres = 0.3
do 21 i=1,nxas
if (clwas(i) .gt. clwthres) then
tas(i) = oldtas(i) + (amkl*clwas(i))
endif
21 continue
return
end
subroutine tcorice(k,nx,ny,t,clwxy,
& told,tnew,diff,lulog,luice,
& atcfid,imon,iday,itime)
c This routine reduces the temperature anomaly via Ben Linstid's
c algorithm. It corrects for ice crystal attenuation where
c temp < (tmean - thrvar); uses Poisson's (actually Laplace's b/c equ=0)
c equation to replace bad temps with ave of nearest 2 to 4 neighbors.
c
dimension t(nx,ny)
dimension clwxy(nx,ny)
dimension told(nx,ny)
dimension tnew(nx,ny)
dimension diff(nx,ny)
c
c Specify parameters
clwaa=0.2
thrvar = 0.5
c tol = 0.001
tol = 0.005
smoothmx = 1000.0
c
c Find mean temperature of points with clwxy < clwaa to establish thresh
tmean = 0.0
count = 0.0
do 22 j=1,ny
do 22 i=1,nx
told(i,j) = t(i,j)
if (clwxy(i,j) .lt. clwaa) then
tmean = tmean + t(i,j)
count = count + 1.0
endif
22 continue
c
tmean = tmean/count
thresh = tmean - thrvar
c
c Count number of points to be corrected
iice = 0
do j=1,ny
do i=1,nx
if (t(i,j) .lt. thresh) iice=iice+1
enddo
enddo
c
c Correct where tprime is less than some threshold using
c Poisson/Laplace's equation of averages. Exit subroutine if
c temperature don't converge before mxsmooth number of iterations.
smoothct = 0.0
99 continue
diffmax = 0.0
smoothct = smoothct + 1.0
if (smoothct .eq. smoothmx) then
write(luice,969) 'NO temp convergence ',k,atcfid,imon,iday,itime
969 format(a20,i2,1x,a6,1x,3(i2.2))
goto 101
endif
c
do 23 j=1,ny
do 23 i=1,nx
if (t(i,j) .lt. thresh) then
if (i.eq.1 .and. j.eq.1) then
tnew(i,j) = (told(i,j+1) + told(i+1,j))/2.0
else if (i.eq.nx .and. j.eq.1) then
tnew(i,j) = (told(i-1,j) + told(i,j+1))/2.0
else if (i.eq.1 .and. j.eq.ny) then
tnew(i,j) = (told(i,j-1) + told(i+1,j))/2.0
else if (i.eq.nx .and. j.eq.ny) then
tnew(i,j) = (told(i-1,j) + told(i,j-1))/2.0
else if (i.eq.1 .and. j.gt.1 .and. j.lt.ny) then
tnew(i,j) = (told(i,j+1) + told(i+1,j) + told(i,j-1))/3.0
else if(i.eq.nx .and. j.gt.1 .and. j.lt.ny) then
tnew(i,j) = (told(i,j+1) + told(i-1,j) + told(i,j-1))/3.0
else if(j.eq.1 .and. i.gt.1 .and. i.lt.nx) then
tnew(i,j) = (told(i+1,j) + told(i,j+1) + told(i-1,j))/3.0
else if(j.eq.ny .and. i.gt.1 .and. i.lt.nx) then
tnew(i,j) = (told(i+1,j) + told(i,j-1) + told(i-1,j))/3.0
else
tnew(i,j) = ((told(i+1,j) + told(i-1,j) + told(i,j+1) +
& told(i,j-1))/4.0)
endif
else
tnew(i,j) = told(i,j)
endif
c
diff(i,j) = abs(tnew(i,j) - told(i,j))
if (diff(i,j) .gt. diffmax) then
diffmax = diff(i,j)
endif
23 continue
c
if (diffmax .gt. tol) then
do 28 j=1,ny
do 28 i=1,nx
told(i,j) = tnew(i,j)
28 continue
goto 99
endif
c
c Copy corrected temperature data into old array to return
do 29 j=1,ny
do 29 i=1,nx
t(i,j) = tnew(i,j)
29 continue
101 continue
c
icount = ifix(count)
write(lulog,400) iice,nx*ny
400 format(/,' Ice correction completed, ',i5,' of ',i5,
+ ' points adjusted')
c
c Output number of smoothing iterations
970 format('Ice correction used ',i6.1, ' smoothing iterations
& for level ',i2)
return
end
subroutine fadj(f,fa)
c This routine adjusts the Coriolis parameter f near the
c equator to keep it from going to zero
c
fmin = 1.0e-6
absf = abs(f)
sgnf = f/absf
c
if (absf .lt. fmin) then
fa = fmin*sgnf
else
fa = f
endif
c
return
end
subroutine wloc(sslat,sslon,aslat,aslon,luloc,np,
+ jyr,jmon,jday,jtimeh,jtimem,atcfid,sname)
c This routine writes the storm center position (sslat,sslon)
c and AMSU data locations to a file for later plotting.
c
dimension aslat(np),aslon(np)
character *6 atcfid
character *9 sname
c
c Write label
write(luloc,210) jyr,jmon,jday,jtimeh,jtimem,atcfid,sname
210 format('AMSU SWATH ',i4.4,1x,i2.2,i2.2,1x,i2.2,i2.2,' UTC ',
+ a6,1x,a9)
c
write(luloc,200) sslon,sslat
200 format(f8.2,f7.2)
c
do 10 i=1,np
write(luloc,200) aslon(i),aslat(i)
10 continue
c
return
end
subroutine barr(fk,flat,flon,nk,efld,exfac,rinf,
+ rlat,rlon,fr,rr,nr,ierr)
c This routine performs a single pass Barnes analysis
c of the unevenly values fk to give fr on an evenly
c spaced radial grid.
c
c Input:
c fk - function values for the analysis
c flat - latitudes (deg N positive) of fk
c flon - longitudes (deg E positive) of fk
c nk - number of fk points
c efld - e-folding radius (km) for Barnes analysis
c exfac- expansion factor to increasing efld as a function of radius
c rinf - Influence radius (km) for Barnes analysis
c rlat - latitude of center of radial grid (deg N)
c rlon - longitude of center of radial grid (deg E)
c rr - radial grid points (m)
c nr - Number of radial grid points
c ierr - error flat (0=normal completion)
c (1=error during analysis)
c
c Output
c fr - Analysis of fk on analysis grid
c
dimension fk(nk),flat(nk),flon(nk)
dimension rr(nr),fr(nr)
c
c Initialize fr to zero
ierr=0
do 10 i=1,nr
fr(i) = 0.0
10 continue
c
rmax = rr(nr)/1000.0
c
c Perform analysis
do 20 i=1,nr
wtsum = 0.0
rkm = rr(i)/1000.0
efldt = efld*(1.0 + exfac*rkm/rmax)
do 30 k=1,nk
call distk(rlon,rlat,flon(k),flat(k),dx,dy,rad)
rad = abs(rad-rkm)
if (rad .le. rinf) then
dnorm = rad/efldt
wt = exp(-dnorm*dnorm)
wtsum = wtsum + wt
fr(i) = fr(i) + wt*fk(k)
endif
30 continue
c
if (wtsum .gt. 0.0) then
fr(i) = fr(i)/wtsum
else
do 40 ii=1,nr
fr(ii) = 0.0
40 continue
ierr=1
return
endif
20 continue
c
return
end
subroutine barxy(fk,flat,flon,nk,efld,rinf,ispf,
+ rlat,rlon,fxy,nx,ny,ierr)
c This routine performs a Barnes analysis
c of the unevenly values fk to give fxy on an evenly
c spaced lat/lon grid. If ispf=1, a second pass is performed.
c
c Input:
c fk - function values for the analysis
c flat - latitudes (deg N positive) of fk
c flon - longitudes (deg E positive) of fk
c nk - number of fk points
c efld - e-folding radius (km) for Barnes analysis
c rinf - Influence radius (km) for Barnes analysis
c rlat - 1-D array of latitudes of analysis grid
c rlon - 1-D array of longitudes of analysis grid
c ispf - flag for second Barnes pass (ispf=1)
c nx - Number of longitudes on analysis grid
c ny - Number of latitudes on analysis grid
c ierr - error flat (0=normal completion)
c (1=error during analysis)
c
c Output
c fxy - Analysis of fk on analysis grid
c
dimension fk(nk),flat(nk),flon(nk)
dimension rlon(nx),rlat(ny),fxy(nx,ny)
c
c Work arrays
parameter (mxas=7000)
dimension fki(mxas),iflag(mxas)
c
c Initialize fxy to zero
ierr=0
do 10 j=1,ny
do 10 i=1,nx
fxy(i,j) = 0.0
10 continue
c
c Perform analysis
do 20 j=1,ny
do 20 i=1,nx
wtsum = 0.0
do 30 k=1,nk
call distk(rlon(i),rlat(j),flon(k),flat(k),dx,dy,rad)
if (rad .le. rinf) then
dnorm = rad/efld
wt = exp(-dnorm*dnorm)
wtsum = wtsum + wt
fxy(i,j) = fxy(i,j) + wt*fk(k)
endif
30 continue
c
if (wtsum .gt. 0.0) then
fxy(i,j) = fxy(i,j)/wtsum
else
do 40 jj=1,ny
do 40 ii=1,nx
fxy(ii,jj) = 0.0
40 continue
ierr=1
return
endif
20 continue
c
if (ispf .le. 0) return
c
c Perform second pass of Barnes analysis
c
c Interpolate analysis variable back to observation locations
c and calcualte deviations from data input values
slat1 = rlat(1)
slon1 = rlon(1)
dlat1 = rlat(2)-rlat(1)
dlon1 = rlon(2)-rlon(1)
izp = 0
c
do k=1,nk
call llintsp(fxy,slon1,slat1,dlon1,dlat1,nx,ny,nx,ny,
+ fsp,flon(k),flat(k),izp,ierrt)
c
if (ierrt .eq. 0) then
fki(k) = fk(k)-fsp
iflag(k) = 0
else
fki(k) = 0.0
iflag(k) = 1
endif
enddo
c
c Perform analysis of deviations
do 60 j=1,ny
do 60 i=1,nx
fxyt = 0.0
wtsum = 0.0
do 70 k=1,nk
call distk(rlon(i),rlat(j),flon(k),flat(k),dx,dy,rad)
if (rad .le. rinf .and. iflag(k) .eq. 0) then
dnorm = rad/efld
wt = exp(-dnorm*dnorm)
wtsum = wtsum + wt
fxyt = fxyt + wt*fki(k)
endif
70 continue
c
if (wtsum .gt. 0.0) then
fxyt = fxyt/wtsum
else
fxyt = 0.0
endif
c
fxy(i,j) = fxy(i,j) + fxyt
60 continue
c
return
end
subroutine altonl(z,t,p,za,ta,pn,ts,ps,nx,ny,np,npn,lulog)
c This routine extracts z,t at NCEP pressure levels
c from za,ta at AMSU pressure levels. Interpolation
c is used if necessary.
c
dimension z(nx,ny,np),t(nx,ny,np),p(np)
dimension za(nx,ny,npn),ta(nx,ny,npn),pn(npn)
dimension ts(nx,ny),ps(nx,ny)
c
c Local variable
dimension iamsuf(100)
c
c Initialize flag variable to zero
do 10 k=1,npn
iamsuf(k) = 0
10 continue
c
c Extract AMSU t,z at NCEP levels for output
do 20 k=1,npn
do 30 kk=1,np
if (p(kk) .eq. pn(k)) then
iamsuf(k) = 1
c
do 40 j=1,ny
do 40 i=1,nx
za(i,j,k) = z(i,j,kk)
ta(i,j,k) = t(i,j,kk)
40 continue
go to 20
endif
30 continue
20 continue
c
c Interpolate between amsu levels for t,z at ncep levels, if necessary
do 41 k=1,npn
if (iamsuf(k) .eq. 0) then
c Find amsu level just above current ncep level
do 42 kk=1,np
if (p(kk) .lt. pn(k)) then
kup = kk
pup = p(kk)
endif
42 continue
c
if (kup .eq. np) then
plo = 1070.0e+2
else
plo = p(kup+1)
endif
c
write(lulog,200) pn(k)/100.0,pup/100.0,plo/100.0
200 format(/,' Interpolate T,Z at p=',f5.0,
+ ' from available AMSU levels ',f5.0,1x,f5.0)
c
do 43 j=1,ny
do 43 i=1,nx
pt = pn(k)
c
p2 = p(kup)
t2 = t(i,j,kup)
z2 = z(i,j,kup)
c
if (kup .eq. np) then
p1 = ps(i,j)
t1 = ts(i,j)
z1 = 0.0
else
p1 = p(kup+1)
t1 = t(i,j,kup+1)
z1 = z(i,j,kup+1)
endif
c
call ztint(p1,p2,z1,z2,t1,t2,pt,zt,tt)
ta(i,j,k) = tt
za(i,j,k) = zt
43 continue
endif
41 continue
c
return
end
subroutine rzwrit(f,scale,nr,dr,nz,dz,slat,slon,atcfid,sname,
+ jyr,jmon,jday,jtimeh,jtimem,luout,lulog,label)
c This routine writes f to a file for later plotting
c
dimension f(nr,nz)
character *6 atcfid
character *10 sname
character *20 label
c
c Local array
parameter (mxd=5000)
dimension dum(mxd)
c
c Make sure dummy array is large enough
npts = nr*nz
if (mxd .lt. npts) then
write(lulog,900) mxd
900 format(/,'Dimension mxd too small in routine rzwrit: ',i4)
stop
endif
c
c Check scale factor
if (scale .le. 0.0) scale = 1.0
c
c Write header lines
write(luout,200) label,jyr,jmon,jday,jtimeh,jtimem,atcfid,sname
200 format(a20,1x,i4.4,1x,i2.2,i2.2,1x,i2.2,i2.2,1x,a6,1x,a10)
write(luout,210) slat,slon,dr/1000.,nr,dz/1000.,nz
210 format(f5.1,1x,f6.1,1x,f6.2,1x,i3,1x,f6.2,1x,i3)
c
c Put f in 1-d array and scale for printing
k = 0
do 10 j=1,nz
do 10 i=1,nr
k = k + 1
dum(k) = scale*f(i,j)
10 continue
c
c Write dum to file
write(luout,220) (dum(k),k=1,npts)
220 format(8(f9.2,1x))
c
return
end
subroutine xywrit(f,scale,nx,rlonl,rlonr,ny,rlatb,rlatt,
+ slat,slon,atcfid,sname,
+ jyr,jmon,jday,jtimeh,jtimem,luout,lulog,label)
c This routine writes f to a file for later plotting
c
dimension f(nx,ny)
character *6 atcfid
character *10 sname
character *20 label
c
c Local array
parameter (mxd=5000)
dimension dum(mxd)
c
c Make sure dummy array is large enough
npts = nx*ny
if (mxd .lt. npts) then
write(lulog,900) mxd
900 format(/,'Dimension mxd too small in routine xywrit: ',i4)
stop
endif
c
c Check scale factor
if (scale .le. 0.0) scale = 1.0
c
c Write header lines
write(luout,200) label,jyr,jmon,jday,jtimeh,jtimem,slat,slon,
+ atcfid,sname
200 format(a20,1x,i4.4,1x,i2.2,i2.2,1x,i2.2,i2.2,1x,f5.1,1x,f6.1,
+ 1x,a6,1x,a10)
write(luout,210) rlonl,rlonr,nx,rlatb,rlatt,ny
210 format(f6.1,1x,f6.1,1x,i3,1x,f6.1,1x,f6.1,1x,i3)
c
c Put f in 1-d array and scale for printing
k = 0
do 10 j=1,ny
do 10 i=1,nx
k = k + 1
dum(k) = scale*f(i,j)
10 continue
c
c Write dum to file
write(luout,220) (dum(k),k=1,npts)
220 format(8(f9.2,1x))
c
return
end
subroutine tcwrit(p,ps,ts,t,clw,nx,ny,np,lutcp,label,
+ x1,x2,y1,y2)
c
c This routine writes the surface P,T, the temperature at the AMSU
c pressure levels and cloud liquid water to a file
c
dimension p(np)
dimension ps(nx,ny),ts(nx,ny)
dimension t(nx,ny,np),clw(nx,ny)
character *20 label
c
c Write file header
write(lutcp,200) label,nx,x1,x2,ny,y1,y2
200 format(' STORMID ',a13,1x,' x: ',i2,1x,f6.1,1x,f6.1,
+ ' y: ',i2,1x,f6.1,1x,f6.1)
c
c Write the cloud liquid water
write(lutcp,201)
201 format(' CLOUD LIQUID WATER')
write(lutcp,210) ((clw(i,j),i=1,nx),j=1,ny)
210 format(10(f6.2,1x))
c
c Write the surface pressure
write(lutcp,202)
202 format(' SURFACE PRESSURE (MB)')
write(lutcp,211) ((ps(i,j)/100.,i=1,nx),j=1,ny)
211 format(10(f6.1,1x))
c
c Write the surface temperature
write(lutcp,203)
203 format(' SURFACE TEMPERATURE (K)')
write(lutcp,210) ((ts(i,j),i=1,nx),j=1,ny)
c
c Write the temperatures at each pressure level
do 10 k=1,np
write(lutcp,204) p(k)/100.
204 format(' TEMPERATURE (K) AT P= ',f5.0,' MB')
write(lutcp,210) ((t(i,j,k),i=1,nx),j=1,ny)
10 continue
c
return
end
subroutine tcrwrit(prz,rhorz,trz,vrz,nr,nz,sslat,sslon,lutcrp,
+ label,coord,times,r1,r2,z1,z2)
c
c This routine writes the p,rho,t and v as a function of r,z
c to a file
c
dimension prz(nr,nz),rhorz(nr,nz),trz(nr,nz),vrz(nr,nz)
character *23 label
c
character *90 coord,times
c
c Write file header
r1k = r1/1000.
r2k = r2/1000.
z1k = z1/1000.
z2k = z2/1000.
c
write(lutcrp,199) coord
write(lutcrp,199) times
199 format(a90)
c
write(lutcrp,200) label,nr,r1k,r2k,nz,z1k,z2k,sslat,sslon
200 format(a23,' r(km): ',i2,1x,f6.1,1x,f6.1,
+ ' z(km): ',i2,1x,f6.1,1x,f6.1,
+ ' lat=',f6.2,' lon=',f7.2)
c
c Write the pressure
write(lutcrp,201)
201 format(' PRESSURE (MB)')
write(lutcrp,211) ((prz(i,j)/100.,i=1,nr),j=1,nz)
211 format(10(f6.1,1x))
c
c Write the density
write(lutcrp,202)
202 format(' DENSITY (KG/M3)')
write(lutcrp,212) ((rhorz(i,j),i=1,nr),j=1,nz)
212 format(10(f6.4,1x))
c
c Write the surface temperature
write(lutcrp,203)
203 format(' TEMPERATURE (K)')
write(lutcrp,213) ((trz(i,j),i=1,nr),j=1,nz)
213 format(10(f6.2,1x))
c
c Write the gradient wind
write(lutcrp,204)
204 format(' GRADIENT WIND (M/S)')
write(lutcrp,214) ((vrz(i,j),i=1,nr),j=1,nz)
214 format(10(f6.2,1x))
c
return
end
subroutine spata(f1,f2,p,nx,ny,np,ibrem,lulog,label)
c This routine calculates area averages of f1 and f2.
c If ibrem=1, a constant value to the first field (f1)
c is added so that the average difference between the fields
c is zero, each each pressure level.
c
dimension f1(nx,ny,np),f2(nx,ny,np),p(np)
character *20 label
c
write(lulog,200) label
200 format(/,'Domain average f1,f2: ',a20)
c
rpts = float(nx*ny)
do 10 k=1,np
f1m = 0.0
f2m = 0.0
do 20 j=1,ny
do 20 i=1,nx
f1m = f1m + f1(i,j,k)
f2m = f2m + f2(i,j,k)
20 continue
f1m = f1m/rpts
f2m = f2m/rpts
c
pmb = p(k)/100.0
c
write(lulog,210) f1m,f2m,pmb
210 format(' f1=',f8.1,' f2=',f8.1,' p=',f5.0)
c
if (ibrem .eq. 1) then
do 30 j=1,ny
do 30 i=1,nx
f1(i,j,k) = f1(i,j,k) + (f2m-f1m)
30 continue
endif
10 continue
c
if (ibrem .eq. 1) then
write(lulog,220)
220 format(' f1 adjusted to make the mean equal to the f2 mean')
endif
c
return
end
subroutine spata1(f1,f2,nx,ny,ibrem,lulog,label)
c This routine calculates area averages of f1 and f2.
c If ibrem=1, a constant value to the first field (f1)
c is added so that the average difference between the fields
c is zero.
c
c This routine is similar to spata, but for a f1 and f2 defined at
c a single pressure level
c
dimension f1(nx,ny),f2(nx,ny)
character *20 label
c
write(lulog,200) label
200 format(/,'Domain average f1,f2: ',a20)
c
rpts = float(nx*ny)
f1m = 0.0
f2m = 0.0
do 20 j=1,ny
do 20 i=1,nx
f1m = f1m + f1(i,j)
f2m = f2m + f2(i,j)
20 continue
f1m = f1m/rpts
f2m = f2m/rpts
c
write(lulog,210) f1m,f2m
210 format(' f1=',f8.1,' f2=',f8.1)
c
if (ibrem .eq. 1) then
do 30 j=1,ny
do 30 i=1,nx
f1(i,j) = f1(i,j) + (f2m-f1m)
30 continue
endif
c
if (ibrem .eq. 1) then
write(lulog,220)
220 format(' f1 adjusted so to make the mean equal to the f2 mean')
endif
c
return
end
subroutine swata(f1,p,mx,np,nx,lulog,label)
c This routine calculates the average of f1 at the swath points
c
dimension f1(mx,np),p(np)
character *20 label
c
write(lulog,200) label
200 format(/,'Swath average f1: ',a20)
c
rpts = float(nx)
do 10 k=1,np
f1m = 0.0
do 20 i=1,nx
f1m = f1m + f1(i,k)
20 continue
f1m = f1m/rpts
c
pmb = p(k)
c
write(lulog,210) f1m,pmb
210 format(' f1=',f6.1,' p=',f6.1)
10 continue
c
return
end
subroutine qualcon(tas,clwas,tpwas,aslat,aslon,pamsu,dumas,
+ dumasp,idumas,mxas,mpas,nxas,np,npst,lulog)
c This routine quality controls the AMSU temperatures by comparing
c them with standard atmosphere temperatures. If the temperature
c difference between the AMSU and standard atmosphere exceeds a
c specified amount (devm) at any of the levels to be analyzed,
c the temperatures for the entire column, the clw and tpw are
c eliminated. The varible containing the number of AMSU swath points
c (nxas) is reduced to account for any deleted data points.
c
dimension tas(mxas,mpas),clwas(mxas),tpwas(mxas)
dimension aslat(mxas),aslon(mxas),pamsu(mpas)
dimension dumas(mxas),dumasp(mxas,mpas),idumas(mxas)
c
c Local arrays
dimension p(200),st(200)
c
c Specify maximum allowable temperature deviation (C)
devm = 40.0
c
write(lulog,200) devm
200 format(/,' Quality control AMSU temperatures with devm= ',f5.1)
c
c Calculate standard atmosphere temperatures
c and extract amsu level pressures
do 10 k=1,np
kk = npst+k-1
pmb = pamsu(kk)
p(k) = pmb
c
if (pmb .lt. 10.0) then
write(lulog,210) pmb
210 format(' Quality control not designed for p<10 mb, p=',f6.1)
stop
endif
c
c Find standard atmosphere temperature
call stndz(pmb,zstd,tstd,thetas)
st(k) = tstd
c
c write(lulog,220) pmb,tstd,zstd/1000.0
c 220 format(1x,' p (mb)=',f6.1,' tstd (K)=',f6.1,1x,
c + ' zstd (km)=',f5.1)
10 continue
c
c Start loop over swath points
do 20 i=1,nxas
c Initialize error flag to zero
idumas(i) = 0
c
c Check all pressure levels at current location
do 30 k=1,np
kk = npst+k-1
dev = abs(tas(i,kk)-st(k))
if (dev .gt. devm) then
idumas(i) = 1
write(lulog,230) aslat(i),aslon(i),p(k),tas(i,kk)
230 format(' Bad temp at lat/lon: ',f7.2,1x,f7.2,
+ ' p=',f6.1,1x,' t=',f8.1)
go to 1000
endif
30 continue
c
1000 continue
20 continue
c
c iland=1
c distth=0.0
c
c if (iland .eq. 0) then
c Eliminate land points
c do 35 i=1,nxas
c write(6,*) aslon(i),aslat(i)
c call aland(aslon(i),aslat(i),dist)
c write(6,*) dist
c if (dist .le. distth) idumas(i) = 1
c 35 continue
c endif
c
c Count the number of bad data points
nbad = 0
do 40 i=1,nxas
if (idumas(i) .eq. 1) then
nbad = nbad + 1
endif
40 continue
c
if (nbad .le. 0) go to 2000
c
nxasq = nxas - nbad
c
c Eliminate bad data points
c ** Temperatures
ii = 0
do 50 i=1,nxas
if (idumas(i) .eq. 0) then
ii = ii + 1
do 55 k=1,np
kk = npst+k-1
dumasp(ii,kk) = tas(i,kk)
55 continue
endif
50 continue
c
do 60 i=1,nxasq
do 60 k=1,np
kk = npst+k-1
tas(i,kk) = dumasp(i,kk)
60 continue
c
c ** clw,tpw,lat,lon
ii = 0
do 70 i=1,nxas
if (idumas(i) .eq. 0) then
ii = ii + 1
dumasp(ii,1) = clwas(i)
dumasp(ii,2) = tpwas(i)
dumasp(ii,3) = aslat(i)
dumasp(ii,4) = aslon(i)
endif
70 continue
c
do 80 i=1,nxasq
clwas(i) = dumasp(i,1)
tpwas(i) = dumasp(i,2)
aslat(i) = dumasp(i,3)
aslon(i) = dumasp(i,4)
80 continue
c
nxas = nxasq
c
2000 continue
write(lulog,240) nbad
240 format(/,i4,' data points eliminated by quality control')
c
return
end
subroutine delim(tas,clwas,tpwas,aslat,aslon,
+ rlatmn,rlatmx,rlatbf,rlonmn,rlonmx,rlonbf,
+ dumas,dumasp,idumas,
+ mxas,mpas,nxas,np,npst,lulog)
c
c This routine eliminates data that is outside of the analysis
c domain. Data is retained in a buffer region outside the domain
c defined by rlonbf and rlatbf.
c
dimension tas(mxas,mpas),clwas(mxas),tpwas(mxas)
dimension aslat(mxas),aslon(mxas),pamsu(mpas)
dimension dumas(mxas),dumasp(mxas,mpas),idumas(mxas)
c
c Initialize elimination flag to zero
do 10 i=1,nxas
idumas(i) = 0
10 continue
c
c Flag the points outside of the domain
rltn = rlatmn-rlatbf
rltx = rlatmx+rlatbf
rlnn = rlonmn-rlonbf
rlnx = rlonmx+rlonbf
c
do 15 i=1,nxas
if (aslat(i) .lt. rltn .or. aslat(i) .gt. rltx) idumas(i)=1
if (aslon(i) .lt. rlnn .or. aslon(i) .gt. rlnx) idumas(i)=1
15 continue
c
c Count the number of excluded data points
nex = 0
do 20 i=1,nxas
if (idumas(i) .eq. 1) then
nex = nex + 1
endif
20 continue
c
if (nex .le. 0) go to 2000
c
nxasq = nxas - nex
c
c Eliminate data points
c ** Temperatures
ii = 0
do 30 i=1,nxas
if (idumas(i) .eq. 0) then
ii = ii + 1
do 35 k=1,np
kk = npst+k-1
dumasp(ii,kk) = tas(i,kk)
35 continue
endif
30 continue
c
do 40 i=1,nxasq
do 40 k=1,np
kk = npst+k-1
tas(i,kk) = dumasp(i,kk)
40 continue
c
c ** clw,tpw,lat,lon
ii = 0
do 50 i=1,nxas
if (idumas(i) .eq. 0) then
ii = ii + 1
dumasp(ii,1) = clwas(i)
dumasp(ii,2) = tpwas(i)
dumasp(ii,3) = aslat(i)
dumasp(ii,4) = aslon(i)
endif
50 continue
c
do 60 i=1,nxasq
clwas(i) = dumasp(i,1)
tpwas(i) = dumasp(i,2)
aslat(i) = dumasp(i,3)
aslon(i) = dumasp(i,4)
60 continue
c
2000 continue
write(lulog,240) nex,nxas
240 format(/,i4,' of ',i4,
+ ' data points eliminated by routine delim')
c
nxas = nxasq
c
return
end
subroutine tdevtest(mxas,nxas,tas,clwas,lutdev,atcfid,
& imon,iday,itime)
c This routine outputs uncorrected temp deviations (from average at
c each pressure level) vs CLW, where tdev(i)=tas(i)-tmean
dimension tas(mxas)
dimension clwas(mxas)
dimension tdev(mxas)
c
c Calculate average of temps at this pressure level
c Initialize tdev array in this loop also
tmean = 0.0
count = 0.0
do 24 i=1,nxas
tmean = tmean + tas(i)
count = count + 1.0
tdev(i) = 0.0
24 continue
tmean=tmean/count
c
do 26 i=1,nxas
tdev(i) = tas(i)-tmean
if (clwas(i) .lt. 0) then
clwas(i) = 0
endif
write(lutdev,971)atcfid,imon,iday,itime,tdev(i),clwas(i)
971 format(a6,3(i2.2),1x,f8.3,1x,f6.2)
26 continue
end
subroutine uvmac(u,v,t,jyr,jlday,jtime,intv)
c This routine prints u and v to an ASCII file that can be used
c to plot data in McIDAS. Every intv point is printed.
c
parameter (nx=61,ny=61)
parameter (npn=12)
c
dimension u(nx,ny,npn),v(nx,ny,npn),t(nx,ny,npn)
c
dimension rlond(nx),rlatd(ny)
dimension rlonr(nx),rlatr(ny)
dimension sinlat(ny),coslat(ny),tanlat(ny)
dimension pn(npn)
c
common /ncepll/ rlond,rlonr,rlatd,rlatr,dlon,dlat,
+ dlonr,dlatr,sinlat,coslat,tanlat
common /ncepp/ pn
c
lumac = 66
c
c Set imac = 1 for original format for J. Dostalek
c or imac = 2 for modified format for J. Evans
imac = 2
c
open(file='uvamsu.dat',unit=lumac,form='formatted',
+ status='unknown')
c
if (jyr .lt. 2000) then
jyr2 = jyr-1900
else
jyr2 = jyr-2000
endif
c
if (imac .eq. 1) then
do 10 j=1,ny,intv
do 10 i=1,nx,intv
do 10 k=1,npn
pmb = pn(k)/100.0
utem = u(i,j,k)
vtem = v(i,j,k)
call ctor(utem,vtem,spd,dir)
dir = 270.-dir
if (dir .lt. 0.0) dir = dir+360.0
c
write(lumac,200) jyr2,jlday,jtime,rlatd(j),-rlond(i),
+ pmb,u(i,j,k),v(i,j,k),spd,dir
200 format(1x,i2.2,i3,1x,i6,7(1x,f6.1))
10 continue
else
call jdayi(jlday,jyr,jmon,jday)
jtime4 = jtime/100
c
do 20 j=1,ny,intv
do 20 i=1,nx,intv
do 20 k=1,npn
pmb = pn(k)/100.0
c
write(lumac,210) jyr,jmon,jday,jtime4,rlatd(j),rlond(i),
+ pmb,u(i,j,k),v(i,j,k),t(i,j,k)
210 format(1x,i4,1x,i2.2,i2.2,1x,i4.4,6(1x,f6.1))
20 continue
endif
c
close(lumac)
return
end
subroutine AMSUfdeck(cbasin,istnum,cstype,
& iyr,mm,dd,hh,min,rlat,rlon,
& mslp,ivmx,ir34,ir50,ir64,irmw,isdist,lufix)
c
c This subroutine creates an ATCF fdeck entry for an AMSU fix.
c
c
c Written by
c John Knaff
c CIRA
c July 2005
c
c Last Modified: August 1, 2005
c
c Language: FORTRAN 90
c
c***********************************************************************
c F-DECK structures
c***********************************************************************
c atcf fix formats as of 4/02/2002
c
c The first 32 words are common followed by specific structure for
c 1)Subjective Dvorak (DVTS)
c 2)Objective Dvorak (DVTO)
c 3)Microwave (SSMI,ERS2,SEAW,TRMM,AMSU,QSCT,ALTI,ADOS)
c 4)Radar (RDRT,RDRD)
c 5)Aircraft (AIRC)
c 6)Dropsonde (DRPS)
c 7)Analysis (ANAL)
c
c This subroutine used structures for (3) above.
c
c**********************************************************************
c
c COMMON FEATURES TO ALL FIXES
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
c
c All character fix record (first 32 words)
c
type cfx_rec
SEQUENCE
character basin*2
character cynum*2
character dtg*12
character fform*3
character ftype*4
character cip*10
character rflag*1
character latns*5
character lonew*6
character hob*5
character cconf*1
character vmax*3
character vconf*1
character mslp*4
character pconf*1
character pder*4
character rad*3
character wcode*4
character rad1*4
character rad2*4
character rad3*4
character rad4*4
character radm1*1
character radm2*1
character radm3*1
character radm4*1
character rconf*1
character rmw*3
character eyed*3
character sreg*1
character fsite*5
character finit*3
end type cfx_rec
c
c Mixed structure... causes a warning message when compiled..
c
c Common part of the fix record (32+2 words )
c
type fx_rec
SEQUENCE
character basin*2
integer cynum
integer dtg
integer fform
character ftype*4
character cip*10
character rflag*1
integer lat
character ns*1
integer lon
character ew*1
integer hob
integer cconf
integer vmax
integer vconf
integer mslp
integer pconf
character pder*4
integer rad
character wcode*4
integer rad1
integer rad2
integer rad3
integer rad4
character radm1*1
character radm2*1
character radm3*1
character radm4*1
integer rconf
integer rmw
integer eyed
character sreg*1
character fsite*5
character finit*3
end type fx_rec
c
c MICROWAVE FIX SPECIFIC FIX STRUCTURE
c
c Character only type
c
type cmifix
sequence
character rflag*1
character rrate*3
character process*6
character waveh*2
character temp*4
character rawslp*4
character retslp*4
character maxsea*3
character stype*6
character rad*3
character wcode*4
character rad1*4
character rad2*4
character rad3*4
character rad4*4
character rad5*4
character rad6*4
character rad7*4
character rad8*4
character radm1*1
character radm2*1
character radm3*1
character radm4*1
character radm5*1
character radm6*1
character radm7*1
character radm8*1
character rconf*1
character comments*52
end type cmifix
c
c Mixed structure type ... causes warnings when compiled
c
type mifix
sequence
character rflag*1
integer rrate
character process*6
integer waveh
integer temp
integer rawslp
integer retslp
integer maxsea
character stype*6
integer rad
character wcode*4
integer rad1
integer rad2
integer rad3
integer rad4
integer rad5
integer rad6
integer rad7
integer rad8
character radm1*1
character radm2*1
character radm3*1
character radm4*1
character radm5*1
character radm6*1
character radm7*1
character radm8*1
integer rconf
character comments*52
end type mifix
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
c End F90 Structures for the Fixes
c
c************************************************************************
c Program defined variables
c************************************************************************
integer iyr,mm,dd,hh,min,sdist,istnum,mslp,ivmx
integer ir34(4),ir50(4),ir64(4)
integer irmw,isdist,lufix
real rlat,rlon
character*4 clat
character*5 clon
character*31 fixname
character*2 cbasin
character*6 cstype
logical lr34,lr50,lr64
type (cfx_rec)fix ! all characters
type (fx_rec) cfx ! mix of types
type (mifix) mfx ! mix of types
type (cmifix) mfix ! all characters
c
c Mark...here's your data statement...preceeding the executables
c
c**********************************************************************
c Data Statements
c**********************************************************************
c Note: The data statements below were converted to assignment statements
c
c data lr34,lr50,lr64/.false.,.false.,.false./
c
c**********************************************************************
c End Variables
c**********************************************************************
c
lr34=.false.
lr50=.false.
lr64=.false.
c
c Set some constants for AMSU retrievals
c
fix%fform=' 30'
fix%ftype='AMSU' ! HARDWIRED TO AMSU
fix%cip=' IP' ! HARDWIRED Intensity and Pressure estimate
fix%rflag=' '
fix%hob=' '
fix%pder='MEAS'
fix%rad=' '
fix%wcode=' '
fix%rad1=' '
fix%rad2=' '
fix%rad3=' '
fix%rad4=' '
fix%radm1=' '
fix%radm2=' '
fix%radm3=' '
fix%radm4=' '
fix%eyed=' '
fix%fsite=' CIRA'
fix%finit='JAK'
fix%cconf='1'
c
mfix%rflag=' '
mfix%rrate=' '
mfix%process=' '
mfix%waveh=' '
mfix%temp=' '
mfix%rawslp=' '
mfix%maxsea=' '
mfix%stype=cstype
mfix%rad=' '
mfix%wcode=' '
mfix%rad1=' '
mfix%rad2=' '
mfix%rad3=' '
mfix%rad4=' '
mfix%rad5=' '
mfix%rad6=' '
mfix%rad7=' '
mfix%rad8=' '
mfix%radm1=' '
mfix%radm2=' '
mfix%radm3=' '
mfix%radm4=' '
mfix%radm5=' '
mfix%radm6=' '
mfix%radm7=' '
mfix%radm8=' '
mfix%comments='storm center extrapolated from t=-12 and t=0 '
. //'adeck'
c
fix%basin=cbasin
write(fix%cynum,15)istnum
15 format(i2.2)
c read(*,15,end=999)fix%basin,fix%cynum
c 15 format(/,19x,a2,a2,/,/,/)
c read(*,30,end=999) iyr,mm,dd,hh,min
c 30 format(23x,i4,1x,i2,i2,1x,i2,i2,/)
write(fix%dtg,30)iyr,mm,dd,hh,min
30 format(i4.4,i2.2,i2.2,i2.2,i2.2)
write(fix%mslp,45)mslp
45 format(i4)
mfix%retslp=fix%mslp
write(fix%vmax,60)ivmx
60 format(i3)
do i=1,4
if(ir34(i).gt.0)lr34=.true.
if(ir50(i).gt.0)lr50=.true.
if(ir64(i).gt.0)lr64=.true.
enddo
write(fix%rmw,75)irmw
75 format(i3)
c
c determine confidence given the distance from the satellite swath center
c (i.e., isdist) passed to the subroutine.
c
cfx%cconf=2
if(isdist.le.300)then
fix%vconf='1'
fix%pconf='1'
fix%rconf='1'
mfix%rconf='1'
else
fix%vconf='2'
fix%pconf='2'
fix%rconf='2'
mfix%rconf='2'
endif
c
c Latitude and Longitude
cfx%lat=nint(rlat*100)
cfx%lon=nint(rlon*100)
if(cfx%lat.ge.0)cfx%ns='N'
if(cfx%lat.lt.0)cfx%ns='S'
if(cfx%lon.ge.0)cfx%ew='E'
if(cfx%lon.lt.0)cfx%ew='W'
if(cfx%lon.gt.18000) then
cfx%lon=cfx%lon-36000.
cfx%ew='W'
endif
write(fix%latns,151)abs(cfx%lat),cfx%ns
151 format(i4,a1)
write(fix%lonew,152)abs(cfx%lon),cfx%ew
152 format(i5,a1)
c
c determine sub-basin
c
if (fix%basin.eq.'WP')then
fix%sreg=fix%basin(1:1)
elseif (fix%basin.eq.'AL')then
fix%sreg='L'
elseif (fix%basin.eq.'CP')then
fix%sreg=fix%basin(1:1)
elseif (fix%basin.eq.'EP')then
fix%sreg=fix%basin(1:1)
elseif (fix%basin.eq.'IO')then
if (rlon .lt. 80.)fix%sreg='A'
if (rlon .ge. 80.)fix%sreg='B'
elseif (fix%basin.eq.'SH')then
if (rlon .lt. 135.0 .and. rlon .ge. 0.)fix%sreg='S'
if (rlon .ge. 135.0 .or. rlon .lt. 0.)fix%sreg='P'
endif
c
c
c write(fixname,200)fix%dtg,fix%cynum,fix%sreg
c200 format('CIRA_',a12,'_',a2,a1,'_FIX')
c open(unit=lufix,file=fixname,status='replace')
c
c write the center/intensity/pressure fix
if(.not.lr34)then
write(lufix,165)fix,mfix
endif
c
c if 34 knot winds exist print the intensity/pressure/Radii fix
c for 34 knot winds
c
if(lr34)then
fix%cip=' IPR'
fix%rad=' 34'
fix%wcode=' NEQ'
mfix%rad=' 34'
mfix%wcode=' NEQ'
write(fix%rad1,175)ir34(1)
write(fix%rad2,175)ir34(2)
write(fix%rad3,175)ir34(3)
write(fix%rad4,175)ir34(4)
175 format(i4)
write(mfix%rad1,175)ir34(1)
write(mfix%rad2,175)ir34(2)
write(mfix%rad3,175)ir34(3)
write(mfix%rad4,175)ir34(4)
write(lufix,165)fix,mfix
endif
c
c if 50 knot winds exist print the Radii fix
c for 50 knot winds
c
if(lr50)then
fix%cip=' R'
fix%rad=' 50'
fix%wcode=' NEQ'
mfix%rad=' 50'
mfix%wcode=' NEQ'
write(fix%rad1,175)ir50(1)
write(fix%rad2,175)ir50(2)
write(fix%rad3,175)ir50(3)
write(fix%rad4,175)ir50(4)
write(mfix%rad1,175)ir50(1)
write(mfix%rad2,175)ir50(2)
write(mfix%rad3,175)ir50(3)
write(mfix%rad4,175)ir50(4)
write(lufix,165)fix,mfix
endif
c
c if 64 knot winds exist print the Radii fix
c for 64 knot winds
c
if(lr64)then
fix%cip=' R'
fix%rad=' 64'
fix%wcode=' NEQ'
mfix%rad=' 64'
mfix%wcode=' NEQ'
write(fix%rad1,175)ir64(1)
write(fix%rad2,175)ir64(2)
write(fix%rad3,175)ir64(3)
write(fix%rad4,175)ir64(4)
write(mfix%rad1,175)ir64(1)
write(mfix%rad2,175)ir64(2)
write(mfix%rad3,175)ir64(3)
write(mfix%rad4,175)ir64(4)
write(lufix,165)fix,mfix
endif
165 format(32(a,', '),28(a,', '),a)
999 close(lufix)
stop
end