subroutine read_ieeetovs(ii,iadtime,ndata,ndatppe,mype,step,start,& 1,6
nlath,nlon,rlats,rlons,obs_amount,obs_boundary,aload,cosfilter,&
ithin,jsatid,gstime,infile,hirs2,hirs3,msu,amsua,amsub,lunout)
!$$$ subprogram documentation block
! subprogram: read_ieeetovs read in tovs data and reformat it
!
! called by read_obs
!
! input argument list:
!
! iadtime - analysis date and time array
! infile - input file array
! jsatid - satellite id
! nlath - number of gaussian lats in one hemisphere
! nlon - num of longitudes
! rlats - grid latitudes (radians)
! rlons - grid longitudes (radians)
!
!
!
! output argument list:
! ndata - number of observation records
! obs_amount - total load for observation computations in pcgsoi
!
!
!
!$$$
use type_kinds
, only: fp_kind,single
implicit none
integer maxchanl,maxblock
parameter (maxchanl=100)
parameter (maxblock=200)
logical hirs2,hirs3,msu,amsua,amsub,cosfilter
character(10) date
character(10),dimension(ndatppe) :: infile
character(8) subset,subfgn
integer nlath,ich2,nrealx,nchanlx,nreal,nchanl,n1data,n1xdata,inadir
integer ilat,ich15,ils,isza,ihgt,ilon,ilza,nthin,nread,lunout,ndatppe
integer ii,mype,nlon,ich8,nmind,isc,nrec,ich10,ith,k,jlon,lonp1,ifov
integer nn,lunin,nlath2,ithx,ithx2,i,iscanpos,lon1,latp1,jlat,lat1
integer,dimension(ndatppe)::ndata,ithin,jsatid
integer,dimension(5):: iadtime
integer idate5(5)
real(fp_kind) gstime,start,aix,rato,dg2rad,sstime,dlat,tdiff,alat,ainc
real(fp_kind) step,sstx,ch8ch10,crit,each,alns,ch8flg,panglr,dlon
real(fp_kind) sch8flg,ch8max
real(fp_kind) rlats(2*nlath),rlons(nlon)
real(fp_kind),dimension(2*nlath,nlon):: obs_amount
real(fp_kind),dimension(2*nlath,nlon,4):: obs_boundary
real(fp_kind),dimension(ndatppe):: aload
real(fp_kind),allocatable,dimension(:,:):: sli,sst
real(fp_kind),dimension(maxchanl,maxblock):: data1b
real(single),dimension(maxchanl):: data1b4
data lunin / 10 /
!**************************************************************************
! Initialize variables
nlath2 = 2*nlath
dg2rad = atan(1.0)/45.0
nread = 0
nthin = 0
ich8 = 10+8
ich10 = 10+10
ich2 = 10+2
ich15 = 10+15
! Open units to satellite data i/o files
open(lunin,file=infile(ii),form='unformatted')
! Read sst and land/sea mask for hirs data
if(hirs2 .or. hirs3)then
allocate(sst(2*nlath,nlon),sli(2*nlath,nlon))
call rdsfull(sst,nlath,nlon,3,mype)
call rdsfull(sli,nlath,nlon,13,mype)
end if
! Set indices for lon,lat location in data array
ils = 7
inadir = 8
ilat = 9
ilon = 10
ilza = 11
isza = 12
ihgt = 13
! Read input file header. Check for bad data file.
read(lunin,end=120,err=120) nreal,nchanl
nrealx=nreal
nchanlx=nchanl
nreal=10
if(hirs2 .or. hirs3)nchanl=19
n1data=nrealx+nchanlx
write(6,*)'READ_IEEETOVS: read 1B header record for ',infile(ii),lunin,&
' with ',nreal,nchanl,n1data
if (abs(n1data)>maxchanl) then
write(6,*)'READ_IEEETOVS: *** WARNING: BAD 1B FILE***'
write(6,*)' iscr,infile=',lunin,infile(ii)
write(6,*)' SKIP PROCESSING OF THIS 1B FILE'
goto 130
end if
write(lunout) nreal,nchanl
n1xdata=nreal+nchanl
! Read and optionally thin data.
aix=0.
if(hirs2 .or. msu) rato=1.1363987
do nrec=1,1000000
ith=0
do while (aix .lt. float(ithin(ii)) .and. ith<maxblock)
read(lunin,end=180,err=120) (data1b4(k),k=1,n1data)
ith=ith+1
ainc=1.
alat=abs(data1b4(ilat))
if(alat .gt. 30. .and. cosfilter)ainc=cos(dg2rad*(alat-30.))
aix=aix+ainc
idate5(1) = data1b4(3) !year
idate5(2) = data1b4(4) !month
idate5(3) = data1b4(5) !day
idate5(4) = data1b4(6)/3600. !hour
idate5(5) = (data1b4(6)-idate5(4)*3600)/60 !minute
isc = data1b4(6)-idate5(4)*3600-idate5(5)*60 !second
call w3fs21(idate5,nmind)
sstime=float(nmind) + isc/60.
tdiff=sstime-gstime
if(tdiff > 180.)then
tdiff = 180.
else if(tdiff < -180.)then
tdiff = -180.
end if
dlat=data1b4(ilat)
dlon=data1b4(ilon)
if(dlon<0.)dlon=dlon+360.
dlat=dlat*dg2rad
dlon=dlon*dg2rad
ifov = nint(data1b4(inadir))
panglr=(start+(ifov-1)*step)*dg2rad
data1b(1,ith) = data1b4(1) ! satellite ID
data1b(2,ith) = tdiff ! time
data1b(3,ith) = dlat ! lat
data1b(4,ith) = dlon ! long
if(hirs2 .or. msu)then
data1b(5,ith) = asin(rato*sin(panglr))
else
data1b(5,ith) = real(data1b4(ilza))*dg2rad ! local zenith angle
if(amsua .and. ifov .le. 15) data1b(5,ith)=-data1b(5,ith)
if(amsub .and. ifov .le. 45) data1b(5,ith)=-data1b(5,ith)
if(hirs3 .and. ifov .le. 28) data1b(5,ith)=-data1b(5,ith)
end if
data1b(6,ith) = panglr ! look angle
data1b(7,ith) = data1b4(inadir)+.001 ! scan position
data1b(8,ith) = data1b4(isza) ! solar zenith angle
data1b(9,ith) = data1b4(ihgt) ! station height
data1b(10,ith)= data1b4(ils)+.001 ! land sea mask
do nn=1,nchanl
data1b(nn+10,ith) = data1b4(nn+nrealx)
end do
end do
aix=aix-float(ithin(ii))
180 continue
nread = nread + ith
if(ith .eq. 0)go to 120
! Set scan position for MSU data
iscanpos = 0
sch8flg=0.
if (msu) then
iscanpos = nint(data1b(7,ith))
else
if(hirs2 .or. hirs3)then
sch8flg=999999.
ch8max=-100.
ithx=ith
ithx2=ith
do i=1,ith
dlat=data1b(3,i)
dlon=data1b(4,i)
call gdcrdp(dlat,1,rlats,nlath2)
call gdcrdp(dlon,1,rlons,nlon)
jlat=nint(dlat)
jlon=nint(dlon)
jlat=max(1,min(jlat,nlath2))
jlon=max(1,min(jlon,nlon+1))
if (jlon==nlon+1) jlon=1
sstx=sst(jlat,jlon)
lat1=int(dlat); lon1=int(dlon)
lat1=min(max(1,lat1),nlath2); lon1=min(max(0,lon1),nlon)
if(lon1==0) lon1=nlon
latp1=min(nlath2,lat1+1); lonp1=lon1+1
if(lonp1==nlon+1) lonp1=1
alns=sli(lat1,lon1)+sli(lat1,lonp1)+sli(latp1,lon1)+sli(latp1,lonp1)
ch8flg = data1b(ich8,i)-sstx
ch8ch10 = data1b(ich8,i)-data1b(ich10,i)
if(jsatid(ii) .eq. 11)ch8flg=ch8flg-2.19955+1.28578*ch8ch10
if(jsatid(ii) .eq. 14)ch8flg=ch8flg-2.61089+1.32519*ch8ch10
if(jsatid(ii) .eq. 15)ch8flg=ch8flg-1.85483+1.30573*ch8ch10
if(jsatid(ii) .eq. 16)ch8flg=ch8flg-1.85483+1.30573*ch8ch10
crit=abs(ch8flg)+abs(data1b(9,i))+20.*alns
if(crit .lt. sch8flg )then
ithx=i
sch8flg=crit
end if
if(data1b(ich8,i) .gt. ch8max)then
ithx2=i
ch8max=data1b(ich8,i)
end if
end do
if(sch8flg .lt. 20.)then
ith=ithx
else
ith=ithx2
end if
else
if(amsua .or. amsub)then
ithx=ith
ith=ithx
end if
end if
end if
! Remove MSU scan positions 1 and 11.
! These are the outermost MSU scan positions.
if ( iscanpos/=1 .and. iscanpos/=11) then
! Increment data counter. Write out data.
ndata(ii) = ndata(ii) + 1
write(lunout) (data1b(k,ith),k=1,n1xdata)
! Compute observational load on analysis grid.
dlat=data1b(3,ith); dlon=data1b(4,ith)
each=aload(ii)
if(sch8flg .gt. 8.)each=.1*each
call sumload(dlat,dlon,each,rlats,rlons,nlath,nlon,&
lat1,lon1,latp1,lonp1,obs_amount,obs_boundary)
endif
! End of loop over input data records
end do
! Jump here when read EOF
120 continue
! Write data counts to stdout
write(6,8000) infile(ii),lunin,mype,nread,ithin(ii),ndata(ii)
8000 format(' READ_IEEETOVS: read EOF for infile=',a10,&
' unit=',i3,' mype=',i3,&
' nread=',i10,' ithin=',i2,' ndata=',i10)
! Jump here when a problem is detected reading 1b header record.
130 continue
if(hirs2 .or. hirs3) deallocate(sst,sli)
return
end subroutine read_ieeetovs