subroutine setuprad(jiter,mype_rad,deltat,sigl,sigi,& 1,18
nsig,lat1,lon1,nsat1,nsatdat_s,lunin,lunout,varst,varprd,&
npred,mype,nsat,nlath,nlon,rlats,rlons,cbias,tlapmean,&
predx,ggrid_g31,ggrid_g32,&
fact10,sfct,isli,sno,veg_type,veg_frac,soil_type,soil_temp,soil_moi, &
diag_rad,idiag_rad,iout_rad,obstype_all,npes,&
jppf,jpch,jpchus,tropprs)
!$$$ subprogram documentation block
! . . . .
! subprogram: setuprad compute rhs of oi equation for radiances
! prgmmr: derber org: w/nmc23 date: 95-07-06
!
! abstract: read in data, first guess, and obtain rhs of oi equation
! for radiances.
!
! program history log:
! 95-07-06 derber
! 96-11 wu data from prepbufr file
! 96-12- mcnally -changes for diagnostic file and bugfix
! 98-04-30 weiyu yang mpi version
! 99-08-24 derber, j., treadon, r., yang, w., first frozen mpp version
!
! input argument list:
! jiter - outer iteration counter
! mype_rad - processor used for print out
! deltat - time interval of first guess fields
! sigl - sigma values at mid-point of each sigma layer
! sigi - sigma values at interfaces of sigma layers
! nsig - number of sigma levels
! lat1 - number of gaussian lats in the sub-domain array
! lon1 - number of gaussian longitudes in the sub-domain array
! nsat1 - number of satellite observations for every satellites
! in each pe subdomain
! nsatdat_s - total number of observations for all satellites
! lunin - unit to read observation distributed data
! lunout - unit to save satllite observation
! processed distributed data
! npred - number of predictors
! mype - pe number
! nsat - number of satellites
! nlath - number of gaussian lats in one hemisphere
! nlon - number of longitudes
! rlats - grid latitudes (radians)
! rlons - grid longitudes (radians)
! cbias - angle dependent bias (fixed from input file)
! tlapmean - mean lapse rate (fixed from input file)
! predx - current estimate of bias correction coefficients
! ggrid_g31- atmospheric state at analysis time
! ggrid_g32- atmospheric state at +/- 3 hrs of analysis time
! fact10 - 10 meter wind factor
! sfct - guess skin temperature
! isli - snow-land-ice mask
! sno - snow-ice mask
! veg_type - veg type
! veg_frac - veg fraction
! soil_type- soil type
! soil_temp- soil temp (first layer)
! soil_moi - soil moisture (first layer)
! gsstm - time index of the first guess fields
! diag_rad - file name for radiance diagnostic file
! iout_rad - output file for printed output
! obstype_all - obs. types for input files
! npes - number of pe's
! jppf - number of profiles processed at once
! jpch - total number of channels
! jpchus - maximum number of channels for any instruments
!
! output argument list:
! varst - variance for skin temperature
! varprd - variance for predictor coefficients
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use type_kinds
, only: fp_kind,single
use error_handler
use initialize
use spectral_coefficients
implicit none
!
! Declare include files
include 'mpi_inc.h'
include 'ijstart.h'
! Declare and set parameters
integer iint,ireal,ipchan
parameter (iint=5, ireal=13, ipchan=7)
!
! Declare variables
!
!
! Declare input/output arrays
real(fp_kind),dimension(nsig):: sigl
real(fp_kind),dimension(nsig+1):: sigi
real(fp_kind),dimension(2*nlath):: varst
real(fp_kind),dimension(jpch*npred):: varprd
real(fp_kind),dimension(2*nlath):: rlats
real(fp_kind),dimension(nlon):: rlons
real(fp_kind),dimension(90,jpch):: cbias
real(fp_kind),dimension(jpch)::tlapmean
real(fp_kind),dimension(jpch,npred):: predx
real(fp_kind),dimension(3,lat1+2,lon1+2,npes):: ggrid_g31
real(fp_kind),dimension(3,lat1+2,lon1+2,npes,2):: ggrid_g32
real(fp_kind),dimension(lat1+2,lon1+2):: fact10,sfct,sno,veg_type,soil_type,soil_temp
real(fp_kind),dimension(lat1+2,lon1+2):: soil_moi,veg_frac,tropprs
integer,dimension(lat1+2,lon1+2):: isli
character(12) diag_rad
logical idiag_rad
integer iout_rad,mype_rad
character(10),dimension(55)::obstype_all
! Declare variables
integer ll,ichan5,iyp,l,ixp,ix1,ix,iy,ichan8,ichan18,ichan4,ichan12,kcount
integer knchpf,m,mm,mype,nlon,nlath,nobs,nblktot,jc,nobs1,nblk,newchn,isz
integer itoss,nchan,nobs2,nsdata,lon1,nsatdat_s,jier,nsig,lunin,npes,jppf
integer jpch,lunout,npred,nsat,k,n1,isat,nele,j,i,is,nsigx1,jiter,lat1,jpchus
integer nstats,nx,mype1,icoast,ireduce,itopo,iobs1,iland,isnoice,ivarl,ich4
integer ich12,ich8,igross,ich2,ich3,j11,nn,iobs2,j12,j21,j22,iang,ich1,iice
integer jsat,icerr,isum,ich21,jo,lnbias,ich5,iemiss,jj,kk,icc,jjj,n,nlev,kval
integer nsig2,nsigx
real(fp_kind) vfact,efact,efactir,efacts,vfactir,efactmc,dtbf,efact21,vfactmc
real(fp_kind) efact22,ch18flg,ch12flg,demisf,fact,dtempf,ch8flg,ch8ch10,sph
real(fp_kind) cenlon,term,sum,tlap,ch6dbc,ch4dbc,dsval,ch15dbc,bcbtm1,emix
real(fp_kind) ch2dbc,ch1dbc,facth6,dtde3,de2,efact6,de3,dtde1,facth4,dtde2
real(fp_kind) de1,vfact22,ch4flg,vfact4,vfact21,ch5flg,vfacts,vfact3,vfact2
real(fp_kind) deltat,diff1,fact1,ch4diff,ch2diff,dsi,si,rsum,svar,cpen,stdev
real(fp_kind) ch5diff,fact3,drad,ch3diff,fact45,fact2,factch6,vfact7,cor
real(fp_kind) ch1diff,factch4,efact7,vfact6,d2,dg2rad,rad2dg,d1,dy,wgt10
real(fp_kind) wgt01,wgt11,dx1,dy1,wgt00,penalty_all,dx,dradnob,varrad,onpp
real(fp_kind) constoz,fifty,zero
! Declare local arrays
real(single) svar4,bcor4(2*(npred+1))
real(fp_kind),dimension(lat1+2,lon1+2):: w10fact
real(fp_kind),dimension(ireal,jppf):: diagbuf
real(fp_kind),dimension(npred+jpchus-1,jppf):: pred
real(fp_kind),dimension(jpchus,jppf):: varinv,errf,btm,tbc,tbcnob,tlapchn,emissav
real(fp_kind),dimension(npred+1,jpchus,jppf):: bcor
real(fp_kind),dimension(jpchus):: rhssat,var,tmpvar,tnoise
real(fp_kind),dimension(1000):: rcount,rcount_all,rpenal
real(fp_kind),dimension(nsig,jppf):: prsr5,temp5,wmix5,omix5,pin5,qvp,tvp,poz
real(fp_kind),dimension(jppf):: ts5,soilt5,vegf5,vegtype5,soilm5,soiltype5,snow5,trop5
! xpath: secant of viewing path angle at surface
! xpaths: secant of solar zenith angle at surface
real(fp_kind),dimension(jppf):: xpath,xpaths,dqval,pangs,panglr,cld,cosza
real(fp_kind),dimension(jppf):: wsp10,cenlat,slats,slons,stfact,pp,zz,sfchgt,zasat
real(fp_kind),dimension(jppf*jpchus):: pems5,ts,pems,ptbcl5
real(fp_kind),dimension(nsig,jpchus*jppf):: wmix,temp,omix,ptau5
real(fp_kind),dimension(3*nsig+1,jpchus):: htlto
real(fp_kind),dimension(15,nsat):: aivals,aivals1
real(fp_kind),dimension(6+2*(npred+1),jpch):: stats,stats1
real(fp_kind),dimension(nsig):: c1,c2,c3,c4
real(fp_kind),dimension(nsig):: siglr,sigir,ccsig
! varch: rms error for each channel
! polar: channel polarization (0=vertical, 1=horizontal, or
! 0 to 1=mix of V and H)
! nuchan: satellite channel numbers
! nusat: satellite numbers
! itype: 0=IR channel, 1=microwave channel
! wvnum: wavenumber in cm**-1
real(fp_kind),dimension(jpch):: varch,polar
integer,dimension(jppf*jpchus):: kchan,kochan,kprof
integer,dimension(jpchus,jppf):: indexn
integer,dimension(jppf):: isflg
integer,dimension(jppf):: lndsea,nadir,knchan
integer,dimension(iint,jppf):: idiagbuf
integer,dimension(55):: kidsat,kidsat1,nsat1
integer,dimension(jpch):: nuchan,nusat,itype,mchannel,ifactq
integer,dimension(0:jpch):: iuse,iouse
integer,dimension(jpchus):: icx,iochan
integer,dimension(jpchus):: ich
character(10),dimension(50)::sattype
character(10) obstype
logical hirs2,msu,goes,hirs3,amsua,amsub,airs,eos_amsua,hsb,goesimg
logical,dimension(jpchus,jppf):: kuse
logical,dimension(jppf):: keep,coast,land,ice
logical delt,first,diagsave
logical ifirst(1000)
real(fp_kind),allocatable,dimension(:,:)::data_s
real(fp_kind),allocatable,dimension(:)::weight_s
!
!
! Initialize variables and constants.
delt = deltat > 0.
first = jiter == 1
diagsave = first .and. idiag_rad
nsig2 = 2*nsig
nsigx = 3*nsig
nsigx1 = 3*nsig+1
nstats = 6+2*(npred+1)
dg2rad = atan(1.0)/45.0
rad2dg = 1.0/dg2rad
d1=.754
d2=-2.265
aivals = 0.0
stats = 0.0
nsatdat_s=0
nx=2*nlath; mype1=mype+1
kidsat=0
iuse(0)=-999
constoz = 604229. * 9.80665 * 21.4e-9
fifty=50.0
zero=0.0
do n=1,nsig
sigir(n)=10.*sigi(nsig-n+1)
siglr(n)=10.*sigl(nsig-n+1)
ccsig(n)=constoz/(sigi(n)-sigi(n+1))
end do
!
!
! Set skin temperature and bias correction standard deviations
varst = sqrt(1.0)
varprd = 1./sqrt(.1)
!
! Scale sigma layer 1 wind by 10m wind factor. Currently, only use
! 6hr guess winds. These winds are used in the emissivity calculation.
! At a later date, obtain 10m guess wind from time interpolated 3,6,9
! hour surface fields
do j = 1,lon1+2
do i = 1,lat1+2
!ibm* prefetch_for_store(w10fact(i,j))
w10fact(i,j)=fact10(i,j)*sqrt(ggrid_g31(1,i,j,1)**2+& ! u.
ggrid_g31(2,i,j,1)**2) ! v.
end do
end do
! Get satellite information
call satinfo(nusat,nuchan,itype,iuse(1),ifactq,mchannel,varch,polar,wavenumber, &
frequency,jpch,mype,mype_rad)
!
rewind lunin
if (diagsave) then
open(4,file=diag_rad,form='unformatted')
rewind 4
! Initialize/write parameters for satellite diagnostic file on
! first outer iteration.
if(mype==0) then
write(4) mype,iint,ireal,ipchan,jpchus
write(6,*)'SETUPRAD: write header record for mype=',&
mype,iint,ireal,ipchan,jpchus,' ',diag_rad
endif
endif
!
!
!****************************************************************************
! MAIN LOOP OVER SATELLITES (ONE SATELLITE PLATFORM PER INPUT FILE)
!
do is = 1,nsat+5
n1 = 0
obstype=' '; isat=0; nele=0; nobs=0
isz=is-5
if(is > 5)sattype(isz)=obstype_all(is)
if (nsat1(is) == 0) goto 135
!
read(lunin,end=135) obstype,isat,nele,nobs
!
! Initialize logical flags for satellite platform
hirs2 = obstype == 'hirs/2'
hirs3 = obstype == 'hirs/3'
msu = obstype == 'msu'
goes = obstype == 'goes'
amsua = obstype == 'amsua'
amsub = obstype == 'amsub'
airs = obstype == 'airs'
hsb = obstype == 'hsb'
eos_amsua = obstype == 'eos_amsua'
goesimg = obstype == 'goesimg'
if(.not. (amsua .or. amsub .or. msu .or.&
hirs2 .or. hirs3 .or. goes .or. &
airs .or. hsb .or. eos_amsua .or. &
goesimg))then
read(lunin)
go to 135
end if
if(obstype /= obstype_all(is) .or. is <= 5)then
write(6,*)' error in setuprad ',is,obstype,obstype_all(is)
call stop2(93)
end if
!
! Based on isat, set new satellite number
!
kidsat(isz)=isat
if(goes)kidsat(isz)=isat+50
if(msu)kidsat(isz)=isat+200
if(goesimg)kidsat(isz)=isat+250
if(amsua)kidsat(isz)=isat+300
if(eos_amsua)kidsat(isz)=isat+300
if(amsub)kidsat(isz)=isat+400
if(hsb)kidsat(isz)=isat+400
if(airs)kidsat(isz)=isat+500
! If AMSU-A data set channel indices for later use
if (amsua .or. eos_amsua) then
ich(1) = newchn(kidsat(isz),1,nusat,nuchan,jpch)
ich(2) = newchn(kidsat(isz),2,nusat,nuchan,jpch)
ich(3) = newchn(kidsat(isz),3,nusat,nuchan,jpch)
ich(4) = newchn(kidsat(isz),4,nusat,nuchan,jpch)
ich(5) = newchn(kidsat(isz),5,nusat,nuchan,jpch)
ich(15)= newchn(kidsat(isz),15,nusat,nuchan,jpch)
endif
nchan=nele-10
if (airs) then
write(6,*) ' airs data currently not set up'
read(lunin)
go to 135
end if
if (nchan > jpchus) write(6,*)'WARNING*** on mype,is=',mype,is,&
nchan,jpchus
tnoise = 1.e10
itoss = 1
do jc=1,nchan
!
! Load channel numbers into local array based on satellite type
iochan(jc)=newchn(kidsat(isz),jc,nusat,nuchan,jpch)
!
! Temporary fix for problems with NOAA 11 channel 9 coefficients
! Use the coefficients from NOAA 14, channel 9.
! if ((kidsat(isz) == 11) .and. (jc == 9) ) &
! iochan(jc) = newchn(14,9,nusat,nuchan,jpch)
!
! Set error instrument channels
iouse(jc)=iuse(iochan(jc))
if((iouse(jc) < 0 .and. .not. first) .or. iouse(jc) < -1) then
tnoise(jc)=1.e10
else
tnoise(jc) = varch(iochan(jc))
end if
if (tnoise(jc) < 1.e4) itoss = 0
end do
if (itoss == 1) then
if(mype == 0)write(6,*)'SETUPRAD: all obs variance > 1e4. do not use ',&
'data from satellite is=',isat,obstype
read(lunin)
goto 135
endif
idiagbuf= 0
diagbuf = 0.
pred=0.0
! Load data array for current satellite
allocate(data_s(nele,nsat1(is)))
allocate(weight_s(nsat1(is)))
read(lunin) data_s,weight_s
nblktot=(nsat1(is)-1)/jppf+1
!
!
!*****
! PROCESS DATA IN BLOCKS OF JPPF
do nblk = 1,nblktot
!
! determine block of data to use
nobs1 = 1 + (nblk-1)*jppf
nobs2 = min(nobs1+jppf-1,nsat1(is))
if (nobs2-nobs1+1 <= 0) then
write(6,*)'setuprad: nblk not consistent with nblktot ',nblk,nblktot
goto 134
end if
!
! Initialize variables/arrays.
tbc = 0.
tbcnob = 0.
bcor = 0.
emissav = 0.
pems5 = 1.
cld = 0.
! INITIAL PROCESSING OF SATELLITE DATA
!*****
!
! Process data
nsdata = 0
do n = nobs1,nobs2
nsdata=nsdata + 1
if (nsdata>jppf .or. nsdata>nsat1(is)) &
write(6,*)'setuprad: ***error*** nsdata,jppf,nsat1(is)=',&
nsdata,jppf,nsat1(is),is,n,nobs1,nobs2,' ',obstype,' ',&
kidsat(isz)
! reuse beginning of weight_s array.
weight_s(nsdata) = weight_s(n)
aivals(1,isz) = aivals(1,isz) + weight_s(nsdata)
! Extract analysis relative observation time.
stfact(nsdata) = data_s(2,n)/deltat
! temporary
if (airs .or. eos_amsua .or. hsb) stfact(nsdata) = 0.
! Extract lon and lat.
cenlat(nsdata) = data_s(3,n)*rad2dg
slats(nsdata) = data_s(3,n)
cenlon = data_s(4,n)*rad2dg
slons(nsdata) = data_s(4,n)
zasat(nsdata) = data_s(5,n)
cosza(nsdata) = cos(zasat(nsdata))
xpath(nsdata) = 1./cosza(nsdata) ! view angle path factor
if(goesimg)then
panglr(nsdata) = 0.
cld(nsdata) = data_s(6,n)
else
panglr(nsdata) = data_s(6,n)
end if
! Extract nadir (scan step position)
nadir(nsdata) = nint(data_s(7,n))
pangs(nsdata) = data_s(8,n)
if (pangs(nsdata) <= 89.) then
xpaths(nsdata) = 1./cos(pangs(nsdata)*dg2rad) ! solar path factor
else ! if sun near/below horizon,
xpaths(nsdata) = 1.0e6 ! set to very large value
endif
sfchgt(nsdata) = data_s(9,n)
! Set land/sea flag. (sea=0, land=1)
lndsea(nsdata)=nint(data_s(10,n))
! Load channel data into work array.
do i = 1,nchan
btm(i,nsdata) = data_s(i+10,n)
end do
! Prepare for application of bias correction to simulated values.
! Construct predictors for 1B radiance bias correction.
pred(1,nsdata) = 0.01
pred(2,nsdata) = .1*(xpath(nsdata)-1.)**2-.015
pred(3,nsdata) = 0.0
! Save data in diagnostic arrays.
idiagbuf(1,nsdata) = kidsat(isz)
idiagbuf(2,nsdata) = nchan
idiagbuf(5,nsdata) = nadir(nsdata)
diagbuf(1,nsdata) = cenlat(nsdata)
diagbuf(2,nsdata) = cenlon
diagbuf(3,nsdata) = rad2dg*zasat(nsdata)
diagbuf(5,nsdata) = pangs(nsdata)
diagbuf(8,nsdata) = cld(nsdata)
diagbuf(9,nsdata) = pred(2,nsdata)
diagbuf(12,nsdata) = weight_s(nsdata)
! End of loop over obs
end do
!
aivals(3,isz) = aivals(3,isz) + nsdata
!
!
call gdcrdp(slats,nsdata,rlats,nx)
call gdcrdp(slons,nsdata,rlons,nlon)
!
! Interpolate model guess to observation locations
!
call intrppx(nx,ggrid_g31,ggrid_g32,&
stfact,sfct,w10fact,veg_type,veg_frac,soil_type,soil_temp, &
soil_moi,sno,tropprs,tvp,qvp,pp,poz,zz,ts5,wsp10,vegtype5,vegf5,soiltype5, &
soilt5,soilm5,snow5,trop5,&
slats,slons,nlon,lat1,lon1,nsig,nsdata,delt,mype,npes)
if (goes .or. airs .or. hsb .or. eos_amsua .or. goesimg) then
do n = 1,nsdata
sfchgt(n) = zz(n)
end do
endif
!
coast = .false.
isflg = 0
keep = .true.
do n = 1,nsdata
!
! Do some preliminary qc of the data. Check for points covered
! with ice/snow, land points, and coastal points.
!
ix1 = slats(n)
iy = slons(n)
ix1 = max(1,min(ix1,nx))
ix = ix1- istart(mype1)+2
iy = iy - jstart(mype1)+2
if (iy<1) iy = iy+nlon
if (iy>lon1+1) iy = iy-nlon
ixp = ix+1
iyp = iy+1
if (ix1==nx) ixp = ix
if (lndsea(n) /= 1)then
if (isli(ix,iy ) == 1 .or. isli(ixp,iy ) == 1 .or.&
isli(ix,iyp) == 1 .or. isli(ixp,iyp) == 1) lndsea(n) = 1
end if
!
! Check for coastal (transition) points
if (isli(ix,iy) /= isli(ixp,iy) .or.&
isli(ix,iy) /= isli(ix ,iyp) .or.&
isli(ix,iy) /= isli(ixp,iyp)) then
coast(n) = .true.
if(isli(ix,iy ) == 1 .or. isli(ixp,iy ) == 1 .or. &
isli(ix,iyp) == 1 .or. isli(ixp,iyp) == 1) lndsea(n)=1
end if
!
! Check for snow of sea ice in surrounding grid box
if(isli(ix,iy ) == 2 .or. isli(ixp,iy ) == 2 .or. &
isli(ix,iyp) == 2 .or. isli(ixp,iyp) == 2) isflg(n)=1
if(isflg(n) /= 1)then
if(sno(ix ,iy ) > 1. .and. isli(ix ,iy ) /= 0) isflg(n)=1
if(sno(ixp,iy ) > 1. .and. isli(ixp,iy ) /= 0) isflg(n)=1
if(sno(ix ,iyp) > 1. .and. isli(ix ,iyp) /= 0) isflg(n)=1
if(sno(ixp,iyp) > 1. .and. isli(ixp,iyp) /= 0) isflg(n)=1
end if
land(n) = lndsea(n) == 1
ice(n) = isflg(n) == 1
!
! Initial QC of data. Set keep to false if we don't want the data.
! The qc tests screen data over land, snow/ice, and coastal points.
! Toss obs over land
if (land(n)) then
keep(n) = .false.
aivals(4,isz) = aivals(4,isz) + weight_s(n)
end if
if (ice(n)) then
!
! Toss obs over snow or ice
keep(n) = .false.
aivals(5,isz) = aivals(5,isz) + weight_s(n)
end if
if (coast(n)) then
!
! Toss obs over coastal points
keep(n) = .false.
aivals(6,isz) = aivals(6,isz) + weight_s(n)
end if
if (hirs2 .or. goes .or. hirs3) then
!
! GROSS QC Check. Compare HIRS and GOES channel 8 with sst. If the
! difference is too large, flag the observation as "bad".
ch8flg = btm(8,n) - ts5(n)
ch8ch10= btm(8,n) - btm(10,n)
if(kidsat(isz) == 11)ch8flg=ch8flg-2.19955+1.28578*ch8ch10
if(kidsat(isz) == 14)ch8flg=ch8flg-2.61089+1.32519*ch8ch10
if(kidsat(isz) == 15)ch8flg=ch8flg-1.85483+1.30573*ch8ch10
if(kidsat(isz) == 16)ch8flg=ch8flg-1.85483+1.30573*ch8ch10
if(kidsat(isz) == 58)ch8flg=ch8flg-2.09303+.277606*ch8ch10
if(kidsat(isz) == 60)ch8flg=ch8flg-4.22641+.331097*ch8ch10
if (abs(ch8flg) > 8.) then
keep(n) = .false.
aivals(10,isz) = aivals(10,isz) + weight_s(n)
endif
endif
!
! Now load/compute values for model levels
snow5(n)=min(snow5(n),fifty)
pp(n) = exp(pp(n))
do l = 1,nsig
pin5(l,n) = sigir(l)*pp(n)
prsr5(l,n) = siglr(l)*pp(n)
end do
do l = 1,nsig
temp5(l,n) = tvp(nsig-l+1,n)/(1.+.61*qvp(nsig-l+1,n))
sph = max(zero,qvp(nsig-l+1,n))
wmix5(l,n) = 1000.*sph/(1.-sph)
omix5(l,n) = max(zero,poz(nsig-l+1,n)*604229.)
end do
!
! Save data in diagnostic arrays
idiagbuf(3,n)= lndsea(n)
idiagbuf(4,n)= isflg(n)
diagbuf(4,n) = ts5(n)
diagbuf(6,n) = sfchgt(n)
diagbuf(7,n) = zz(n)
diagbuf(13,n)= trop5(n)
!
! End loop over obs
end do
!
!
!
!*****
! IDENTIFY CHANNELS TO PROCESS
!*****
!
! Load arrays used by forward model to indicate which
! profiles and which channels to compute radiances for.
!
! The tests below screen from the dataset those profiles/
! channels which are not used (as of 3 Dec 01).
!
kchan = 0
kochan = 0
kprof = 0
knchan = 0
errf = 0.0
varinv = 1.e-12
kuse = .false.
!
! 1B HIRS and GOES: retain all channels over water
! over land, retain channels 2-5, 12.
!
if (hirs2 .or. hirs3 .or. goes) then
do j = 1,nsdata
do jc = 1,nchan
kuse(jc,j)=(keep(j) .or. jc <= 5 &
.or. jc == 12) .and. tnoise(jc) < 1.e4
end do
end do
elseif(msu .or. amsua .or. eos_amsua .or. amsub .or. hsb .or. goesimg)then
do j = 1,nsdata
do jc = 1,nchan
!
! retain all data, provided assigned channel error
! is "small"
kuse(jc,j)= tnoise(jc) < 1.e4
end do
end do
endif
knchpf = 0
indexn = 0
do j = 1,nsdata
kcount = 0
do jc = 1,nchan
if(kuse(jc,j))then
kcount = kcount+1
knchpf = knchpf+1
kchan(knchpf) = iochan(jc)
kochan(knchpf) = jc
kprof(knchpf) = j
varinv(jc,j) = 1.0/tnoise(jc)**2
errf(jc,j) = tnoise(jc)
indexn(kcount,j) = knchpf
endif
end do
knchan(j) = kcount
end do
if (knchpf < 1) go to 130
!
! Compute surface emissivity for observation
!
call emiss(knchpf,kprof,kchan,kochan,wsp10,zasat,panglr,&
isflg,lndsea,emissav,pems5,ts5,soiltype5,soilt5, &
soilm5,vegtype5,vegf5,snow5,polar,frequency,mchannel, &
itype,nsdata,jppf,jpch,jpchus)
!
! call radiative transfer (forward) model
!
call rad_tran_k(pin5,prsr5,temp5,wmix5,omix5,ts5,&
pems5,xpath,xpaths,knchan,kchan,ptau5,ptbcl5, &
temp,wmix,omix,ts,pems, &
itype,nsig,jpch,knchpf,nsdata)
!
!*****
! COMPUTE AND APPLY BIAS CORRECTION TO SIMULATED VALUES
!*****
! Compute predictor for AMSU-A cloud liquid water bias correction.
do i=1,knchpf
n=kprof(i)
tbcnob(kochan(i),n) = ptbcl5(i)
tbc(kochan(i),n) = ptbcl5(i) + cbias(nadir(n),kchan(i))
end do
if (amsua .or. eos_amsua) then
do n=1,nsdata
dqval(n)=0.0
if (tbc(1,n)<=284. .and. tbc(2,n)<=284. .and. keep(n)) &
dqval(n)=d1*(tbc(1,n)-btm(1,n))/(285.-tbc(1,n))+ &
d2*(tbc(2,n)-btm(2,n))/(285.-tbc(2,n))
if (btm(1,n)<0. .or. btm(2,n)<0.) then
dqval(n)=0.0
varinv(1,n)=0.0
varinv(2,n)=0.0
varinv(3,n)=0.0
varinv(4,n)=0.0
varinv(5,n)=0.0
varinv(6,n)=0.0
varinv(7,n)=0.0
varinv(15,n)=0.0
endif
pred(3,n) = dqval(n)*cosza(n)*cosza(n)
end do
endif
! Perform bias correction
do i=1,knchpf
n=kprof(i)
m=kochan(i)
mm=kchan(i)
tlapchn(m,n)= .01*(ptau5(nsig-1,i)-ptau5(nsig,i))* &
(ts5(n)-temp5(nsig-1,n))
do l=2,nsig-1
tlapchn(m,n)=tlapchn(m,n)+&
.01*(ptau5(l-1,i)-ptau5(l,i))*(temp5(l+1,n)-temp5(l-1,n))
end do
tlap = tlapchn(m,n)-tlapmean(mm)
sum = tlap*(predx(mm,npred)+tlap*predx(mm,npred-1))
bcor(1,m,n) = cbias(nadir(n),mm) !global_satangbias
bcor(2,m,n) = tlap*predx(mm,npred) !tlap
bcor(3,m,n) = tlap*tlap*predx(mm,npred-1) !tlap*tlap
do j = 1,npred-2
term = predx(mm,j)*pred(j,n)
sum = sum + term
bcor(j+3,m,n) = term
end do
tbc(m,n) = tbc(m,n) + sum
end do
!
!
!
!******
! QC OBSERVATIONS BASED ON VARIOUS CRITERIA
! Separate blocks for various instruments.
!******
!
!
! QC HIRS/2, GOES, HIRS/3 data
!
if (hirs2 .or. goes .or. hirs3 .or. airs) then
!default channel numbers
ichan8 = 8
ichan18 = 18
ichan12 = 12
ichan4 = 4
ichan5 = 5
! use proxy channels for AIRS
!============================
if (airs) then
ichan8 = 8
ichan18 = 18
ichan12 = 12
ichan4 = 4
ichan5 = 5
endif
!
! Loop over observations.
do n = 1,nsdata
efact = 1.
vfact = 0.25
efactir = 1.
vfactir = 1.
efacts = 1.
vfacts = 1.
vfact2 = 1.0
vfact3 = 1.0
vfact4 = 1.0
!
! Apply window check. Toss obs if fail window test.
ch4flg = btm(ichan4,n) - tbc(ichan4,n)
ch5flg = btm(ichan5,n) - tbc(ichan5,n)
ch8flg = btm(ichan8,n) - tbc(ichan8,n)
ch12flg = btm(ichan12,n)- tbc(ichan12,n)
ch18flg = btm(ichan18,n)- tbc(ichan18,n)
demisf = .01
dtempf = .5
!
! Cloud test for channel 2
if (ch4flg < -20.0) then
vfact2 = 0.0
vfact3 = 0.0
vfact4 = 0.0
vfactir = 0.0
aivals(11,isz) = aivals(11,isz) + weight_s(n)
else if (ch4flg < -1.0) then
!
! Cloud test for channel 3
vfact3 = 0.0
vfact4 = 0.0
vfactir = 0.0
aivals(12,isz) = aivals(12,isz) + weight_s(n)
else if(abs(ch4flg) > 1.0 .or. abs(ch5flg) > 1.0)then
vfact4 = 0.0
vfactir = 0.0
aivals(13,isz) = aivals(13,isz) + weight_s(n)
else if (ch12flg < -10.0)then !ch12 sees cloud
vfactir = 0.0
vfact4 = 0.0
aivals(14,isz) = aivals(14,isz) + weight_s(n)
! Channel 8 cloud test
! + chan 18 cloud test?
else if (abs(ch8flg) > 1.0 .or. .not. keep(n))then
! .or. (pangs(n) > 85. .and. abs(ch18flg) > 2.0) )then
!test
!test
vfactir = 0.0
if(keep(n))aivals(15,isz) = aivals(15,isz) + weight_s(n)
end if
!
! Reduce qc bounds in tropics
if ( (cenlat(n) > -25.) .and.&
(cenlat(n) < 25.) ) then
aivals(7,isz) = aivals(7,isz) + weight_s(n)
efact = (abs(cenlat(n))*.5/25.+.5)*efact
endif
!
! Reduce weight given to obs for shortwave ir if
! solar zenith angle small
if (pangs(n) <= 60. ) then
vfacts = .5*vfacts
efacts = .5*efacts
endif
if(ice(n))then
demisf = .01
dtempf = 3.0
else if(land(n))then
demisf = .01
dtempf = 2.0
end if
!
!
! Reduce weight for obs over higher topography
if (sfchgt(n) > 2000.) then
fact = 2000./sfchgt(n)
efactir = fact*efactir
vfact4 = fact*vfact4
! if(sfchgt(n) > 2000)then
! vfact4 = 0.0
! end if
! write(6,*)'qc: height=',mype,n,sfchgt(n),fact,efactif
endif
!
! Reduce weight if model and obs topography too different
if (abs(zz(n)-sfchgt(n)) > 200.) then
fact = 200./abs(zz(n)-sfchgt(n))
efactir = fact*efactir
vfact4 = fact*vfact4
aivals(8,isz) = aivals(8,isz) + weight_s(n)
endif
!
!
!
! Generate q.c. bounds and modified variances.
do ll=1,knchan(n)
l = kochan(indexn(ll,n))
varinv(l,n)=vfact*varinv(l,n)
errf(l,n)=3.*efact*errf(l,n)
! Channel 2 test
if (l == 2) then
varinv(l,n) = vfact2*varinv(l,n)
!
! Channel 3 test
else if (l == 3) then
varinv(l,n) = vfact3*varinv(l,n)
! Channel 4 and Channel 12 test
else if (l == 4 .or. l == 12) then
varinv(l,n) = vfact4*varinv(l,n)
else if (l /= 1)then
!
! Channels 4 and above
errf(l,n) = efactir*errf(l,n)
varinv(l,n) = vfactir*varinv(l,n)
if (l >= 13) then
errf(l,n) = efacts*errf(l,n)
varinv(l,n) = vfacts*varinv(l,n)
end if
end if
dtbf = demisf*abs(pems(indexn(ll,n)))+dtempf*abs(ts(indexn(ll,n)))
term = dtbf*dtbf
if (varinv(l,n)>1.e-12 .and. term>1.e-12) &
varinv(l,n)=1./(1./varinv(l,n)+term)
! emix=abs(pems(indexn(ll,n)))
! if(emix > .3 .and. (land(n) .or. ice(n) &
! ) .and. varinv(l,n) > 1.e-9)then
! print *,' emiss ',l,ll,emix,lndsea(n),isflg(n)
! varinv(l,n)=0.0
! end if
end do
!
!
! End of loop over observations.
end do
!
! End of HIRS and GOES QC blocks
!
!
! QC MSU data
else if (msu) then
!
! Loop over observations.
do n = 1,nsdata
efact = 1.
vfact = 0.25
efactmc = 1.
efact21 = 1.
efact22 = 1.
vfactmc = 1.
vfact21 = 1.
vfact22 = 1.
!
! Reduce qc bounds in tropics
if ( (cenlat(n) > -25.) .and.&
(cenlat(n) < 25.) ) then
aivals(7,isz) = aivals(7,isz) + weight_s(n)
efact = (abs(cenlat(n))*.5/25.+.5)*efact
endif
demisf = .015
dtempf = 0.5
!
! Reduce q.c. bound over land
if (land(n)) then
! efact22 = 0.0
! vfact22 = 0.0
demisf = .03
dtempf = 2.5
end if
if (coast(n)) then
! efact22 = 0.0
! vfact22 = 0.0
demisf = .20
dtempf = 4.5
end if
!
! Over water apply window test to channel 2 using channel 1
if (.not. land(n) .or. coast(n)) then
if (abs(btm(1,n) - tbc(1,n)) > 5.0) then
efact22 = 0.0
vfact22 = 0.0
aivals(10,isz) = aivals(10,isz) + weight_s(n)
endif
endif
!
! Reduce q.c. bound for likely snow or ice
if (ice(n)) then
! efact22 = .5*efact22
! vfact22 = .5*vfact22
demisf = .05
dtempf = 3.0
end if
!
!
! Reduce q.c. bounds over higher topography
if (sfchgt(n) > 2000.) then
fact = 2000./sfchgt(n)
efactmc = fact*efactmc
vfactmc = fact*vfactmc
end if
!
! Reduce q.c. bounds if model and obs topography
! too different.
if (abs(zz(n)-sfchgt(n)) > 200.) then
fact = 200./abs(zz(n)-sfchgt(n))
efactmc = fact*efactmc
vfactmc = fact*vfactmc
aivals(8,isz) = aivals(8,isz) + weight_s(n)
end if
!
! Generate q.c. bounds and modified variances.
errf(1,n) = efactmc*efact21*errf(1,n)
errf(2,n) = efactmc*efact22*errf(2,n)
errf(3,n) = 2.0*efactmc*errf(3,n)
errf(4,n) = 2.0*errf(4,n)
varinv(1,n) = vfactmc*vfact21*varinv(1,n)
varinv(2,n) = vfactmc*vfact22*varinv(2,n)
varinv(3,n) = vfactmc*varinv(3,n)
varinv(4,n) = vfactmc*varinv(4,n)
do ll=1,knchan(n)
l = kochan(indexn(ll,n))
errf(l,n) = 3.*efact*errf(l,n)
varinv(l,n) = vfact*varinv(l,n)
emix=abs(pems(indexn(ll,n)))
dtbf = demisf*abs(pems(indexn(ll,n)))+dtempf*abs(ts(indexn(ll,n)))
term = dtbf*dtbf
if (varinv(l,n)>1.e-12 .and. term>1.e-12) &
varinv(l,n)=1./(1./varinv(l,n)+term)
if(l == 1)diagbuf(9,n) = dtbf
if(l == 2)diagbuf(10,n) = dtbf
if(l == 3)diagbuf(11,n) = dtbf
! if(emix > .3 .and. (land(n) .or. ice(n) &
! ) .and. varinv(l,n) > 1.e-9)then
! print *,' emiss msu',l,ll,emix,lndsea(n),isflg(n)
! varinv(l,n)=0.0
! end if
end do
!
! End of loop over observations.
end do
!
! End of MSU QC block
!
!
! QC AMSU-A data
else if (amsua .or. eos_amsua) then
! Loop over observations.
do n = 1,nsdata
if(knchan(n) > 0)then
bcbtm1=btm(1,n) -cbias(nadir(n),ich(1))
! bcbtm2=btm(2,n) -cbias(nadir(n),ich(2))
! bcbtm15=btm(15,n)-cbias(nadir(n),ich(15))
ch1dbc=btm(1,n)-tbc(1,n)
ch2dbc=btm(2,n)-tbc(2,n)
ch4dbc=btm(4,n)-tbc(4,n)
ch6dbc=btm(6,n)-tbc(6,n)
ch15dbc=btm(15,n)-tbc(15,n)
! sval=-113.2+(2.41-.0049*bcbtm1)*bcbtm1+.454*bcbtm2-bcbtm15
dsval=(2.41-.0098*bcbtm1)*ch1dbc+.454*ch2dbc-ch15dbc
! factch6x=((sval-5.)/10.)**2+(ch6dbc/.8)**2
factch6=(dsval/10.)**2+(ch6dbc/.8)**2
factch4=(cosza(n)*dqval(n)/.3)**2+(ch4dbc/1.8)**2
! QC based on ratio of obs-ges increment versus the sensitivity of
! the simulated brightness temperature to the surface emissivity
! Y2K hurricane season runs by QingFu Liu found the hurricane
! forecast tracks to be degraded without this QC.
!
dtde1 = pems(indexn(1,n))
de1 = 0.
if (dtde1 /= 0.) de1=abs((btm(1,n)-tbcnob(1,n))/dtde1)
dtde2 = pems(indexn(2,n))
de2 = 0.
if (dtde2 /= 0.) de2=abs((btm(2,n)-tbcnob(2,n))/dtde2)
dtde3 = pems(indexn(3,n))
de3 = 0.
if (dtde3 /= 0.) de3=abs((btm(3,n)-tbcnob(3,n))/dtde3)
diagbuf(8,n) = factch4
! diagbuf(9,n) = de1
! diagbuf(10,n) = de2
! diagbuf(11,n) = de3
efact = 1.0
efact6 = 1.0
efact7 = 1.0
efactmc = 1.0
vfact = 1.0
vfact6 = 1.0
vfact7 = 1.0
vfactmc = 1.0
demisf = .01
dtempf = .5
if (ice(n)) then
! Decrease likelyhood of using data over snow and ice
demisf = .05
dtempf = 3.0
else if(land(n))then
!
! Decrease likelyhood and weight of data over land
demisf = .02
dtempf = 2.0
end if
!
! Reduce qc bounds in tropics
if ( (cenlat(n) > -25.) .and.&
(cenlat(n) < 25.) ) then
aivals(7,isz) = aivals(7,isz) + weight_s(n)
efact = (abs(cenlat(n))*.25/25.+.75)*efact
endif
!
if(factch6 >= 1. .or. coast(n))then
efactmc=0.0
vfactmc=0.0
efact6=0.0
vfact6=0.0
if(.not. coast(n))aivals(11,isz) = aivals(11,isz) + weight_s(n)
else if (factch4 > 1. .or. de2 > .03 .or. &
de3 > .05 .or. de1 > .05) then
!
efactmc=0.0
vfactmc=0.0
if(.not. coast(n) .and. factch4 > 1.) &
aivals(10,isz) = aivals(10,isz) + weight_s(n)
else if (abs(zz(n)-sfchgt(n)) > 200.) then
!
! Reduce q.c. bounds if model and obs topography
! too different.
fact = 200./abs(zz(n)-sfchgt(n))
efactmc = fact*efactmc
vfactmc = fact*vfactmc
aivals(8,isz) = aivals(8,isz) + weight_s(n)
end if
if (sfchgt(n) > 2000.) then
!
! Reduce q.c. bounds over higher topography
fact = 2000./sfchgt(n)
efactmc = fact*efactmc
efact6 = fact*efact6
vfactmc = fact*vfactmc
vfact6 = fact*vfact6
if(sfchgt(n) > 4000.)then
fact = 4000./sfchgt(n)
efact7 = fact*efact7
vfact7 = fact*vfact7
end if
end if
!
! Generate q.c. bounds and modified variances.
do ll=1,knchan(n)
l = kochan(indexn(ll,n))
errf(l,n) = 3.*efact*errf(l,n)
varinv(l,n) = vfact*varinv(l,n)
if(l <= 5 .or. l == 15)then
! Adjust observation error based on magnitude of liquid water correction
! 0.2 is empirical factor
cor=0.2*(predx(ich(l),3)*pred(3,n))**2
varinv(l,n) = varinv(l,n)/(1.+varinv(l,n)*cor)
errf(l,n) = efactmc*errf(l,n)
varinv(l,n) = vfactmc*varinv(l,n)
else if(l == 6)then
errf(l,n)=efact6*errf(l,n)
varinv(l,n) = vfact6*varinv(l,n)
else if(l == 7)then
errf(l,n)=efact7*errf(l,n)
varinv(l,n) = vfact7*varinv(l,n)
end if
dtbf = demisf*abs(pems(indexn(ll,n)))+dtempf*abs(ts(indexn(ll,n)))
term = dtbf*dtbf
if (varinv(l,n)>1.e-12 .and. term>1.e-12) &
varinv(l,n)=1./(1./varinv(l,n)+term)
! if(l == 3)diagbuf(8,n) = dtbf
if(l == 4)diagbuf(9,n) = dtbf
if(l == 5)diagbuf(10,n) = dtbf
if(l == 6)diagbuf(11,n) = dtbf
if ( (l <= 5 .or. l == 15) .and. &
(varinv(l,n)<1.e-9) ) pred(3,n) = 0.0
end do
end if
!
! End of loop over observations.
end do
!
! End of AMSU-A QC block
!
! QC AMSU-B data
else if (amsub .or. hsb) then
!
! Loop over observations.
do n = 1,nsdata
if(knchan(n) > 0)then
ch1diff=btm(1,n)-tbc(1,n)
ch2diff=btm(2,n)-tbc(2,n)
efact = 1.0
efactmc = 1.0
vfact = 1.0
vfactmc = 1.0
if(land(n))then
dsi=0.85*ch1diff-ch2diff
si=42.72+0.85*btm(1,n)-btm(2,n)
demisf = .03
dtempf = 2.0
else
dsi=999.
if(btm(2,n) < 300.)then
dsi=0.13*(ch1diff-33.58*ch2diff/(300.-btm(2,n)))
end if
si=42.72+0.85*ch1diff-ch2diff
demisf = .015
dtempf = 0.5
end if
if(ice(n))then
demisf = .05
dtempf = 3.0
end if
if(coast(n))then
demisf = .25
dtempf = 5.0
end if
diff1=ch1diff-7.5*dsi
fact1=(diff1/10.)**2+(dsi/1.)**2
diagbuf(9,n) = diff1
diagbuf(10,n) = dsi
diagbuf(11,n) = fact1
!
! Reduce q.c. bounds over higher topography
if (sfchgt(n) > 2000.) then
fact = 2000./sfchgt(n)
efactmc = fact*efactmc
vfactmc = fact*vfactmc
end if
!
! Reduce q.c. bounds if model and obs topography
! too different.
if (abs(zz(n)-sfchgt(n)) > 200.) then
fact = 200./abs(zz(n)-sfchgt(n))
efactmc = fact*efactmc
vfactmc = fact*vfactmc
aivals(8,isz) = aivals(8,isz) + weight_s(n)
end if
if(fact1 > 1.)then
efactmc=0.0
vfactmc=0.0
if(coast(n))aivals(10,isz) = aivals(10,isz) + weight_s(n)
else
efact = (1.-fact1*fact1)*efact
vfact = (1.-fact1*fact1)*vfact
end if
!
! Generate q.c. bounds and modified variances.
do ll=1,knchan(n)
l = kochan(indexn(ll,n))
errf(l,n) = 3.*efact*efactmc*errf(l,n)
varinv(l,n) = vfact*vfactmc*varinv(l,n)
dtbf = demisf*abs(pems(indexn(ll,n)))+dtempf*abs(ts(indexn(ll,n)))
term = dtbf*dtbf
if (varinv(l,n)>1.e-12 .and. term>1.e-12) &
varinv(l,n)=1./(1./varinv(l,n)+term)
end do
!
end if
! End of loop over observations.
end do
!
! End of AMSU-B QC block
!
! GOES imager Q C
!
else if(goesimg)then
! Loop over observations.
do n = 1,nsdata
ch4diff = btm(4,n) - tbcnob(4,n)
ch5diff = btm(5,n) - tbcnob(5,n)
ch2diff = btm(2,n) - tbcnob(2,n)
ch3diff = btm(3,n) - tbcnob(3,n)
if(knchan(n) > 0)then
efact = 1.0
vfact = 1.0
fact45 = 1.0
fact2 = 1.0
fact3 = 1.0
demisf = .01
dtempf = 0.5
if(land(n))then
fact45=0.
fact2=0.
if(cld(n) <40.0 ) fact3=0.0
demisf = .01
dtempf = 2.0
end if
if(ice(n))then
fact45=0.
fact2=0.
fact3=0.0
demisf = .02
dtempf = 3.0
end if
if(coast(n))then
fact45=0.
fact2=0.
fact3=0.0
demisf = .02
dtempf = 5.0
end if
!
!
! Reduce weight for obs over higher topography
! if(rad2dg*zasat(n) >60.)then
! vfact=0.
! aivals(10,isz)= aivals(10,isz) + weight_s(n)
! end if
if (sfchgt(n) > 2000.) then
fact = 2000./sfchgt(n)
efact = fact*efact
vfact = fact*vfact
aivals(11,isz)= aivals(11,isz) + weight_s(n)
! if(sfchgt(n) > 2000)then
! vfact4 = 0.0
! end if
! write(6,*)'qc: height=',mype,n,sfchgt(n),fact,efactif
end if
if( ch4diff <-3.6 .or. ch5diff <-3.6 ) then
fact3=0.0
aivals(13,isz)=aivals(13,isz)+weight_s(n)
end if
if (ch4diff >= -3.6 .and. ch4diff <-3.0) then
if(ch3diff <0.0 ) then
fact3=0.0
aivals(13,isz)=aivals(13,isz)+weight_s(n)
end if
end if
if (ch5diff >= -3.6 .and. ch5diff <-3.0) then
if(ch3diff <0.0 ) then
fact3=0.0
aivals(13,isz)=aivals(13,isz)+weight_s(n)
end if
end if
! end of channel 3 quality control
if(ch4diff <-3.0 .or. ch4diff > 5.0 .or.ch5diff <-3.0 .or. ch5diff > 5.0)then
fact45 = 0.
aivals(14,isz)= aivals(14,isz) + weight_s(n)
end if
if(pangs(n) >= 90. )then
if(ch2diff < -2.5)then
fact2 = 0.
aivals(12,isz)= aivals(12,isz) + weight_s(n)
end if
else
fact2 = 0.
end if
!
! Generate q.c. bounds and modified variances.
do ll=1,knchan(n)
l = kochan(indexn(ll,n))
errf(l,n) = 3.*efact*errf(l,n)
! if(l == 3) write(99,*) 'GOESIMG:errf=',errf(l,n),efact,l
varinv(l,n) = vfact*varinv(l,n)
if(l == 4 .or. l == 5)varinv(l,n) = fact45*varinv(l,n)
if(l == 2)varinv(l,n) = fact2*varinv(l,n)
if(l == 3) varinv(l,n) = fact3*varinv(l,n)
dtbf = demisf*abs(pems(indexn(ll,n)))+dtempf*abs(ts(indexn(ll,n)))
term = dtbf*dtbf
if (varinv(l,n)>1.e-12 .and. term>1.e-12) &
varinv(l,n)=1./(1./varinv(l,n)+dtbf**2)
if(l == 3)diagbuf(7,n) = sum
if(l == 3)diagbuf(9,n) = cbias(nadir(n),65)
if(l == 3)diagbuf(10,n) = pems(indexn(ll,n))
if(l == 3)diagbuf(11,n) = ts(indexn(ll,n))
end do
end if
! End of loop over observations.
end do
end if
!
!
!
! Apply gross check to observations. Toss obs failing test.
do k = 1,knchpf
if (varinv(kochan(k),kprof(k)) > 1.e-6) then
drad = abs(btm(kochan(k),kprof(k)) - tbc(kochan(k),kprof(k)))
!
! If mean obs-ges difference around observations
! location is too large and difference at the
! observation location is similarly large, then
! toss the observation.
if (drad > errf(kochan(k),kprof(k)) ) then
varinv(kochan(k),kprof(k)) = 0.
stats(2,kchan(k)) = stats(2,kchan(k)) + weight_s(kprof(k))
aivals(9,isz) = aivals(9,isz) + weight_s(kprof(k))
end if
!
end if
end do
if(amsua .or. amsub .or. msu .or. hsb)then
if(amsua .or. eos_amsua)nlev=6
if(amsub)nlev=5
if(hsb)nlev=4
if(msu)nlev=4
do n = 1, nsdata
kval=0
do k=2,nlev
if(varinv(k,n) < 1.e-6)then
kval=k-1
if(amsub .or. hsb)kval=nlev
end if
end do
if(kval > 0)then
do k=1,kval
varinv(k,n)=0.
end do
if(amsua)varinv(15,n)=0.
end if
end do
end if
do k = 1,knchpf
n = kprof(k)
mm = kochan(k)
if (varinv(mm,n) > 1.e-6) then
!
! Accumulate data for diagnostic print.
drad = btm(mm,n) - tbc(mm,n)
dradnob = btm(mm,n) - tbcnob(mm,n)
varrad = drad*varinv(mm,n)
m = kchan(k)
if (iuse(m) == 1)aivals(2,isz) = aivals(2,isz) + drad*varrad*weight_s(n)
stats(1,m) = stats(1,m) + weight_s(n) !number of obs
stats(3,m) = stats(3,m) + drad*weight_s(n) !obs-mod(w_biascor)
stats(4,m) = stats(4,m) + drad*drad*weight_s(n) !(obs-mod(w_biascor))**2
stats(5,m) = stats(5,m) + drad*varrad*weight_s(n)!penalty contribution
stats(6,m) = stats(6,m) + dradnob* weight_s(n) !obs-mod(w/o_biascor)
do j=1,npred+1 !jj=7=global_satangbias.txt
!jj=8=tlap, jj=9=tlap*tlap
!jj=10=mean, jj=11=xpath, jj=12=clw
jj = j+6
jjj=jj+6 !jjj=jj+6 is square of jj term
term = bcor(j,mm,n)
stats(jj, m)= stats(jj, m) + term*weight_s(n)
stats(jjj,m)= stats(jjj,m) + term*term*weight_s(n)
end do
endif
end do
!
!
!******
! CONSTRUCT SENSITIVITY VECTORS. WRITE TO OUTPUT FILE.
!******
!
! Generate o3, t, q, and ts sensitivity arrays.
do n = 1,nsdata
icc = 0
onpp=1./pp(n)
do i=1,nsig
c1(i)=ccsig(i)*onpp
c2(i)=1./(1.+.61*qvp(i,n))
c3(i)=1000./(1.-qvp(i,n))**2
c4(i)=.61*tvp(i,n)*c2(i)*c2(i)
end do
do k = 1,knchan(n)
kk = indexn(k,n)
j = kochan(kk)
jo = kchan(kk)
! Only "good" obs are included in calculation.
if (varinv(j,n) > 1.e-6 .and. iuse(jo)==1 ) then
icc = icc+1
icx(icc) = jo
var(icc) = varinv(j,n)
rhssat(icc) = btm(j,n)-tbc(j,n)
pred(npred-2+icc,n)=tlapchn(j,n)-tlapmean(jo)
do i = 1,nsig
htlto(i,icc) = temp(nsig-i+1,kk)*c2(i)
htlto(nsig+i,icc) = c3(i)*wmix(nsig-i+1,kk) - c4(i)*temp(nsig-i+1,kk)
htlto(nsig2+i,icc) = omix(nsig-i+1,kk)*c1(i)
end do
htlto(nsigx+1,icc) = ts(kk)
! If we do not use channel 9, turn off ozone sensitivity in all channels.
if ((hirs2 .or. goes .or. hirs3) .and. (varinv(9,n) < 1.e-6))then
do i = 1,nsig
htlto(nsig2+i,icc) = 0.0
end do
endif
! Reduce amplitude of q sensitivity above tropopause
do i=1,nsig
if (pin5(i,n) < trop5(n)) then
term = (pin5(i,n)-trop5(n))/(trop5(n)-pin5(1,n))
htlto(nsig+i,icc) = exp(ifactq(jo)*term)*htlto(nsig+i,icc)
endif
end do
!
end if
! End loop over channels.
end do
!
! Load data into output arrays
if (icc>0) then
n1 = n1+1
iy = slats(n)
ix = slons(n)
dx = slons(n) - ix
dy = slats(n) - iy
iy = min(max(1,iy),nx)
ix = min(max(0,ix),nlon)
dx1 = 1.0-dx
dy1 = 1.0-dy
iy = max(1,min(iy,nx))
iy = iy-istart(mype1)+2
ix = ix-jstart(mype1)+2
if (ix < 1) ix = ix + nlon
if (ix > lon1 + 1) ix = ix - nlon
wgt00=dx1*dy1; wgt10=dx1*dy; wgt01=dx*dy1; wgt11=dx*dy
j11=iy+(ix-1)*(lat1+2); j12=j11+lat1+2; j21=j11+1; j22=j12+1
! Write data to output file
write(lunout) j11,j21,j12,j22,wgt00,wgt10,wgt01,wgt11,icc,weight_s(n)
write(lunout) (icx(m),m=1,icc),(var(m),m=1,icc),&
(pred(m,n),m=1,npred+icc-2),(rhssat(m),m=1,icc),&
((htlto(m,nn),m=1,nsigx1),nn=1,icc)
!
end if
if (diagsave) then
!
! Write diagnostics to output file. Only generate
! file on first outer iteration.
do nn=1,nchan
tmpvar(nn) = varinv(nn,n)
if (iouse(nn)<1) tmpvar(nn)=-tmpvar(nn)
end do
write(4) mype,(idiagbuf(nn,n),nn = 1,iint),&
(diagbuf(nn,n),nn = 1,ireal),&
(float(nn),btm(nn,n),&
tbc(nn,n),tbcnob(nn,n),tmpvar(nn),&
emissav(nn,n),tlapchn(nn,n),nn=1,nchan)
endif
!
! End of loop over profiles
end do
!
!
! End of nblk loop
130 continue
end do
!
!
! Jump here when there is no more data to process for current satellite
134 continue
deallocate(data_s)
deallocate(weight_s)
135 continue
nsatdat_s=nsatdat_s+n1
!
!
!
!*****
! END MAIN LOOP OVER SATELLITES
!*****
end do
! Close unit to diagnostic file.
if (diagsave)close(4)
!
! Collect statistics
call mpi_reduce(aivals,aivals1,15*nsat,mpi_rtype,mpi_sum,mype_rad,mpi_comm_world,ierror)
call mpi_reduce(stats,stats1,nstats*jpch,mpi_rtype,mpi_sum,mype_rad,mpi_comm_world,ierror)
call mpi_reduce(kidsat,kidsat1,nsat,mpi_integer,mpi_max,mype_rad,mpi_comm_world,ierror)
! Compute and print statistics for current satellite.
if(mype==mype_rad) then
penalty_all=0.0
do is=1,nsat
obstype=sattype(is)
!
! Initialize logical flags for satellite platform
hirs2 = obstype == 'hirs/2'
hirs3 = obstype == 'hirs/3'
msu = obstype == 'msu'
goes = obstype == 'goes'
amsua = obstype == 'amsua'
amsub = obstype == 'amsub'
airs = obstype == 'airs'
hsb = obstype == 'hsb'
eos_amsua = obstype == 'eos_amsua'
goesimg = obstype == 'goesimg'
iobs2 = nint(aivals1(1,is))
penalty_all=penalty_all+aivals1(2,is)
if(iobs2 > 0)then
iobs1 = nint(aivals1(3,is))
iland = nint(aivals1(4,is))
isnoice = nint(aivals1(5,is))
icoast = nint(aivals1(6,is))
ireduce = nint(aivals1(7,is))
itopo = nint(aivals1(8,is))
ivarl = nint(aivals1(9,is))
igross = nint(aivals1(10,is))
if(hirs2 .or. hirs3 .or. goes .or. airs)then
ich2 = nint(aivals1(11,is))
ich3 = nint(aivals1(12,is))
ich4 = nint(aivals1(13,is))
ich12 = nint(aivals1(14,is))
ich8 = nint(aivals1(15,is))
write(iout_rad,2000) 'sat','type','num','numw','ich2','ich3', &
'ich4','ich8','ich12'
write(iout_rad,2010) kidsat1(is),sattype(is),iobs1,iobs2,ich2,ich3, &
ich4,ich8,ich12
write(iout_rad,2012) 'penalty','iland','isnoice','ireduce','icoast', &
'itopo','ivarl','igross'
write(iout_rad,2011) aivals1(2,is),iland,isnoice,ireduce,icoast,itopo,ivarl,igross
else if(msu)then
ich21 = nint(aivals1(10,is))
write(iout_rad,2000) 'sat','type','num','numw','ich21'
write(iout_rad,2010) kidsat1(is),sattype(is),iobs1,iobs2,ich21
write(iout_rad,2012) 'penalty','iland','isnoice','ireduce','icoast', &
'itopo','ivarl'
write(iout_rad,2011) aivals1(2,is),iland,isnoice,ireduce,icoast,itopo,ivarl
else if(amsua .or. eos_amsua)then
iemiss = nint(aivals1(10,is))
iice = nint(aivals1(11,is))
write(iout_rad,2000) 'sat','type','num','numw','ch4','ice'
write(iout_rad,2010) kidsat1(is),sattype(is),iobs1,iobs2,iemiss,iice
write(iout_rad,2012) 'penalty','iland','isnoice','ireduce','icoast', &
'itopo','ivarl'
write(iout_rad,2011) aivals1(2,is),iland,isnoice,ireduce,icoast,itopo,ivarl
else if(amsub .or. hsb)then
ich1 = nint(aivals1(10,is))
write(iout_rad,2000) 'sat','type','num','numw','ich1'
write(iout_rad,2010) kidsat1(is),sattype(is),iobs1,iobs2,ich1
write(iout_rad,2012) 'penalty','iland','isnoice','ireduce','icoast', &
'itopo','ivarl'
write(iout_rad,2011) aivals1(2,is),iland,isnoice,ireduce,icoast,itopo,ivarl
else if(goesimg)then
iang = nint(aivals1(10,is))
ich1 = nint(aivals1(11,is))
ich2 = nint(aivals1(12,is))
ich5 = nint(aivals1(13,is))
write(iout_rad,2000) 'sat','type','num','numw','hi topo','ch45','ch2'
write(iout_rad,2010) kidsat1(is),sattype(is),iobs1,iobs2,ich1,ich2,ich5
write(iout_rad,2012) 'penalty','iland','isnoice','ireduce','icoast', &
'itopo','ivarl','iang'
write(iout_rad,2011) aivals1(2,is),iland,isnoice,ireduce,icoast,itopo,ivarl,iang
else if (airs) then
write(iout_rad,2990)
else
write(iout_rad,2999)
end if
write(iout_rad,'(/)')
end if
end do
write(iout_rad,*)'rad total penalty_all=',penalty_all
! Print counts, bias, rms, stndev as a function of channel.
rcount=0.0; rcount_all=0.0; rpenal=0.0
ifirst=.false.
if (first) then
lnbias=59
open(lnbias,file='satbias_pred',form='unformatted')
end if
do i = 1,jpch
isum = nint(stats1(1,i))
if (isum > 0) then
jsat = nusat(i)
svar=varch(i)
rcount_all(jsat) = rcount_all(jsat) + isum
if (iuse(i)==1) then
rcount(jsat) = rcount(jsat) + isum
rpenal(jsat) = rpenal(jsat) + stats1(5,i)
else
svar=-svar
end if
ifirst(jsat)=.true.
rsum = 1./float(isum)
icerr = nint(stats1(2,i))
do j=3,6 ! j=3=obs-mod(w_biascor)
! j=4=(obs-mod(w_biascor))**2
! j=5=penalty contribution
! j=6=obs-mod(w/o_biascor)
stats1(j,i) = stats1(j,i)*rsum
end do
stats1(4,i) = sqrt(stats1(4,i))
if (isum > 1) then
stdev = sqrt(stats1(4,i)*stats1(4,i)-stats1(3,i)*stats1(3,i))
else
stdev = 0.
end if
write(iout_rad,1102) i,nuchan(i),nusat(i),isum,icerr,svar,&
stats1(6,i),stats1(3,i),stats1(5,i),stats1(4,i),stdev
! Print average magnitude of correction for each bias correction term.
! j=7=satangbias.txt, j=8=tlap, j=9=tlap*tlap, j=10=mean, j=11=xpath, j=12=clw
! j=13 through j=18 are square of j=7 through j=12
svar4=svar
do j=7,nstats
bcor4(j-6)=stats1(j,i)
end do
if(first)write(lnbias) i,nuchan(i),nusat(i),isum,svar4,bcor4
else
svar4=-999.
if(first)write(lnbias) i,nuchan(i),nusat(i),isum,svar4,bcor4
endif
end do
if(first)close(lnbias)
! Write obs count to runtime output file
write(iout_rad,1109)
do i=1,jpch
jsat=nusat(i)
if (ifirst(jsat)) then
ifirst(jsat)=.false.
cpen=0.0
if (rcount(jsat)>0.) cpen=rpenal(jsat)/rcount(jsat)
write(iout_rad,1115) jiter,jsat,rcount_all(jsat),rcount(jsat),&
rpenal(jsat),cpen
endif
end do
2000 format(a7,2x,A4,6x,8(a7,1x))
2010 format(i7,1x,A10,1x,8(i7,1x))
2011 format(3x,f16.8,7(i7,1x))
2012 format(7x,A7,5x,7(a7,1x))
2990 format(' airs data not implemented at this time ')
2999 format(' Illegal satellite type ')
1102 format(1x,2i3,i4,2i7,1x,f8.3,1x,6(f11.7,1x))
1109 format(t5,'it',t13,'sat',t20,'num all',t31,'num use',t41,'penalty',t56,&
'cpen')
1115 format('o-g',1x,i2.2,1x,'rad',2x,i3,2x,f9.0,2x,f9.0,2x,g12.6,2x,g12.6)
! End of diagnostic print block.
close(iout_rad)
end if
! End of routine
return
end subroutine setuprad