subroutine setuppcp(jiter,deltat,sigl,sigi,& 1,7
deltim,dtphys,jcap,&
nsig,lat1,lon1,&
npcp1,npcpdat_s,&
lunin,lunout,varprdp,npredp,predxp,&
mype,nsat,&
nlath,nlon,rlats,rlons,rbs2,ggrid_g31,ggrid_g32,rdivs,&
rcwm,rtlat,rtlon,rqlat,rqlon,&
sfct,isli,sno,diag_pcp,idiag_pcp,mype_pcp,iout_pcp,&
obstype_all,npes,&
xkt2d,jppfp,jtype)
!$$$ subprogram documentation block
! . . . .
! subprogram: setuppcp 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
! 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
! nn_pcp - length of satellites data
! npcp1 - number of satellite observations for every satellites
! in each pe subdomain
! npcpdat_s - total number of observations for all satellites
! npredp - number of predictors
! mype - pe number
! npcp - number of satellites
! dtskin - skin temperature increment in each subdomain
! nlath - number of gaussian lats in one hemisphere
! nlon - number of longitudes
! rlats - grid latitudes (radians)
! rlons - grid longitudes (radians)
! rz_ges - guess mountains on grid
! rp_ges - guess log(sfcp) on grid
! ru_ges - guess u on grid
! rv_ges - guess v on grid
! rt_ges - guess t on grid
! rq_ges - guess specific humidity on grid
! roz_ges - guess ozone on grid in each pe subdomain
! sfct - guess skin temperature
! isli - snow-land-ice mask
! sno - snow-ice mask
!
! output argument list:
! varprdp - variance for predictor coefficients
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
!
! Include data type module
use type_kinds
, only: fp_kind,single
! Turn off implicit typing
implicit none
! Other include files
include 'mpi_inc.h'
include 'constant.h'
include 'ijstart.h'
! Declare and set parameters
integer iint,ireal,ncld
parameter (ncld=1)
parameter (iint=8, ireal=23)
!
! Declare variables
logical idiag_pcp,sea
logical ssmi,amsua,tmi
logical,dimension(jppfp):: land,coast,ice
character(12) diag_pcp
character(10),dimension(50)::sattype
character(10) obstype
character(10),dimension(55)::obstype_all
integer ii,i,j,k,m,n,nlath2,nsdata,is,nele,nobs1,nobs2,nobs
integer isat,isz,lunin,lunout,k1,ns2,nsig1,ipt,jppfp,jtype
integer lat1,lon1,nlath,nlon,nsig,jcap,nsphys,ix,im
integer mype,j11,j22,j12,j21,iy,iy1,ix1,ixp,iyp,npes
integer npredp,nblk,nblktot,iland
integer iout_pcp,mype_pcp,km,mype1
integer isum,icoast,icerr,kx
integer isnoice,itotal,ireduce,jiter
integer n1,npcpdat_s,nsat,kdt,iratio
integer ncloud,igradsml,igradbig
integer inumw,inum,iclip,ikeep,itoss
integer,dimension(55):: npcp1,npcp1d
integer,dimension(lat1+2,lon1+2):: isli
integer,dimension(iint,jppfp):: idiagbuf
integer,dimension(jppfp,4):: iloc,jloc,ijloc
integer,dimension(jppfp):: kbcon,jmin,ktcon,kuo,istype,&
lndsea,ksatid,isflg
integer,dimension(nsat):: kidsat,kidsat1
integer,dimension(jtype):: nupcp,itypep,iusep,ibias,ifutr
real(fp_kind) avg,sdv,rterm1,rterm2,rterm,zero,one,oneneg
real(fp_kind) pcpobs0,pcpnbc0,pcpbc0,pcpobs1,pcpnbc1,pcpbc1,error,a0,a1,simerr
real(fp_kind) pcpsas0,pcpsas1,errlog,rsum
real(fp_kind) rdocp,frain,dtp,dtf,dtphys,deltat,deltim,sum,rad2dg,halfpi
real(fp_kind) wgt00,wgt01,wgt10,wgt11,drad,vfact,efact,fhour,rtime
real(fp_kind) penalty_all,dx,dx1,dy,dy1,varrad,xlo,xhi
real(fp_kind) obserr_wgt,simerr_wgt,wgtij
real(fp_kind) rmiss,pterm,tf1,tf2,tf3,tfo
real(fp_kind) rmmhr,term
real(single),dimension(ireal,jppfp):: diagbuf
real(fp_kind),dimension(jtype):: varchp,vfutr
real(fp_kind),dimension(2*nlath):: rlats,rbs2
real(fp_kind),dimension(nlon):: rlons
real(fp_kind),dimension(jtype*npredp):: varprdp
real(fp_kind),dimension(nsig):: sigl,del,slk
real(fp_kind),dimension(nsig+1):: sigi
real(fp_kind),dimension(npredp,jppfp):: pred
real(fp_kind),dimension(jppfp):: plon0,plat0,rmask0,xkt2,psexp0,&
rbs0,ts0,cldwrk,rn1,pcpsas,rn4,pcpsas4,pcplrg4,pcpnbc,pcpbc,pcplrg,&
ad_rn1,var,obs,ges,varinv,errf,obserr,slats,slons,stfact,cenlat,&
satpcp,satclw,satcnv,satcli,pcpnum,pcpsdv
real(fp_kind),dimension(jppfp,4):: wgtloc
real(fp_kind),dimension(nsig):: dpcpdt,dpcpdq,dpcpdu,dpcpdv,dpcpddiv,dpcpdcwm
real(fp_kind),dimension(2):: dpcpmag
real(fp_kind),dimension(jtype,npredp):: predxp
real(fp_kind),dimension(36,nsat):: aivals,aivals1
real(fp_kind),dimension(nsat):: rpen,cpen
real(fp_kind),dimension(lat1+2,lon1+2):: sfct,sno,xkt2d
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,nsig,3):: rtlat,rtlon,rqlat,rqlon,rcwm,rdivs
real(fp_kind),allocatable,dimension(:):: weight_p
real(fp_kind),allocatable,dimension(:,:):: data_p,tlat0,tlon0,qlat0,qlon0,&
t0,q0,cwm0,u0,v0,div0,t1,q1,cwm1,u1,v1,ad_t0,ad_q0,ad_cwm0,&
ad_u0,ad_v0,ad_div0,ad_t1,ad_q1,ad_cwm1,ad_u1,ad_v1,ad_t4,&
ad_q4,ad_cwm4,ad_u4,ad_v4,ad_div4,tsas,qsas,tlrg,qlrg
data rmiss / -999. /
!
!
!
!*********************************************************************************
! ONE TIME, INITIAL SETUP PRIOR TO PROCESSING SATELLITE DATA
! Allocate arrays
allocate(tsas(nsig,jppfp),qsas(nsig,jppfp),tlrg(nsig,jppfp),qlrg(nsig,jppfp))
allocate(t0(nsig,jppfp),q0(nsig,jppfp),cwm0(nsig,jppfp),u0(nsig,jppfp),&
v0(nsig,jppfp),div0(nsig,jppfp),t1(nsig,jppfp),q1(nsig,jppfp),&
cwm1(nsig,jppfp),u1(nsig,jppfp),v1(nsig,jppfp),tlat0(nsig,jppfp),&
tlon0(nsig,jppfp),qlat0(nsig,jppfp),qlon0(nsig,jppfp))
allocate(ad_t0(nsig,jppfp),ad_q0(nsig,jppfp),ad_cwm0(nsig,jppfp),&
ad_u0(nsig,jppfp),ad_v0(nsig,jppfp),ad_div0(nsig,jppfp),&
ad_t1(nsig,jppfp),ad_q1(nsig,jppfp),ad_cwm1(nsig,jppfp),&
ad_u1(nsig,jppfp),ad_v1(nsig,jppfp),ad_t4(nsig,jppfp),&
ad_q4(nsig,jppfp),ad_cwm4(nsig,jppfp),ad_u4(nsig,jppfp),&
ad_v4(nsig,jppfp),ad_div4(nsig,jppfp))
!
! Initialize variables
ncloud = ncld
nsphys = max(int(2*deltim/dtphys+0.9999),1)
dtp = 2*deltim/nsphys
dtf = 0.5*dtp
frain = dtf/dtp
kdt = 1
fhour = kdt*deltim/3600.
rtime = 1./(3600.*fhour)
rmmhr = 1.e3*rtime * 3600.
nlath2=2*nlath
mype1=mype+1
nsig1=nsig+1; ns2=2*nsig+2
zero=0.0
one=1.0
oneneg=-1.0
call pcpinfo(mype,jtype,nupcp,itypep,iusep,ibias,ifutr,varchp,vfutr)
if (jiter==1 .and. mype==mype_pcp) then
write(iout_pcp,*)'deltim,dtphys =',deltim,dtphys
write(iout_pcp,*)'nsphys,dtp,dtf=',nsphys,dtp,dtf
write(iout_pcp,*)'frain,fhour =',frain,fhour
write(iout_pcp,*)'rtime,rmmhr =',rtime,rmmhr
write(iout_pcp,*)'kdt =',kdt
do j=1,jtype
write(iout_pcp,1005) nupcp(j),itypep(j),iusep(j),ibias(j),varchp(j)
end do
1000 format(1X,5I5,2E18.10)
1005 format(' ipcp = ',i4, ' itypep= ',i4, &
' iusep = ',i4, ' ibias = ',i4, &
' var = ',e10.4)
endif
!
!
! Open units. Set pointers.
rewind(lunin)
kidsat=0
if (jiter==1 .and. idiag_pcp) then
open(4,file=diag_pcp,form='unformatted')
rewind 4
endif
!
! Initialize/write parameters for satellite diagnostic file on
! first outer iteration.
if (jiter==1 .and. idiag_pcp) then
if(mype==0) then
write(4) mype,iint,ireal
write(6,*)'SETUPPCP: write header record for mype=',&
mype,iint,ireal,' ',diag_pcp
endif
endif
!
!
! Initialize variables.
rdocp = rd/cp
!
!
! Set other constants.
rad2dg = 45./atan(1.0)
halfpi = 0.5*acos(-1.)
!
!
! Set bias correction standard deviations
varprdp = 1./sqrt(.1)
aivals = 0.0; aivals1=0.0
npcp1d = 0; n1=0
!
!
!****************************************************************************
! MAIN LOOP OVER OBSERVATION PLATFORMS
!
!*****
! INITIALIZE VARIABLE AND DATA ARRAYS
!*****
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 (npcp1(is) == 0) goto 135
read(lunin,end=135) obstype,isat,nele,nobs
! Initialize logical flags for satellite platform
ssmi = obstype == 'pcp_ssm/i'
tmi = obstype == 'pcp_tmi'
amsua = obstype == 'pcp_amsua'
if (.not. (ssmi .or. tmi .or. amsua) ) then
read(lunin)
goto 135
endif
if (obstype /= obstype_all(is) .or. is <=5) then
write(6,*)' error in setuppcp ',is,obstype,obstype_all(is)
call stop2(93)
end if
! Save isat. Initialize arrays
kidsat(isz) = isat
idiagbuf= 0
diagbuf = 0.
!
! Load data array for current satellite
allocate(data_p(nele,npcp1(is)))
allocate(weight_p(npcp1(is)))
read(lunin) data_p,weight_p
!
!*****
! PROCESS DATA IN BLOCKS OF JPPFP
nblktot=(npcp1(is)-1)/jppfp+1
do nblk = 1,nblktot
!
! Determine block of data to use
nobs1 = 1 + (nblk-1)*jppfp
nobs2 = min(nobs1+jppfp-1,npcp1(is))
if (nobs2-nobs1+1 <= 0) then
write(6,*)'setuppcp: nblk not consistent with nblktot ',nblk,nblktot
goto 134
end if
!
! Initialize variables/arrays.
varinv = 0.
errf = 0.
satpcp = 0.
satcnv = 0.
satclw = 0.
satcli = 0.
rn1 = 0.
pcpsas = 0.
pcplrg = 0.
lndsea = 0
stfact = 0.
idiagbuf= 0
diagbuf = 0.
slats = 0.
slons = 0.
cenlat = 0.
ts0 = 0.
! INITIAL PROCESSING OF SATELLITE DATA
!*****
!
! Process data
nsdata=0
do n=nobs1,nobs2
nsdata = nsdata+1
if (nsdata>jppfp .or. nsdata>npcp1(is)) &
write(6,*)'setuppcp: ***error*** nsdata,jppfp,npcp1(is)=',&
nsdata,jppfp,npcp1(is),is,n,nobs1,nobs2,' ',obstype,' ',&
kidsat(isz)
! Reuse beginning of weight_p array
weight_p(nsdata) = weight_p(n)
aivals(1,isz) = aivals(1,isz) + weight_p(nsdata)
! Extract satellite id
ksatid(nsdata) = nint(data_p(1,n))
kx = -999
do i = 1,jtype
if (ksatid(nsdata).eq.nupcp(i)) kx = i
end do
if (kx<0) write(6,*) &
'***WARNING*** Problem with satellite id. ksatid,kx=',&
ksatid(nsdata),kx
istype(nsdata) = kx
!
! Set observation error
errf(nsdata) = varchp(kx)
varinv(nsdata) = 1./(varchp(kx)**2)
! Extract land/sea/ice flag (0=sea, 1=land, 2=ice)
lndsea(nsdata) = nint(data_p(5,n))
sea = lndsea(nsdata) == 0
!
! Extract obs lon and lat.
slons(nsdata) = data_p(4,n)
slats(nsdata) = data_p(3,n)
cenlat(nsdata) = slats(nsdata)*rad2dg
!
! Extract obs date/time.
stfact(nsdata) = data_p(2,n)/3.
!
! Extract observations
if (ssmi) then
satpcp(nsdata) = data_p(6,n)
pcpsdv(nsdata) = data_p(7,n)
pcpnum(nsdata) = data_p(8,n)
satcnv(nsdata) = rmiss
satclw(nsdata) = rmiss
satcli(nsdata) = rmiss
elseif (amsua) then
! nothing
elseif (tmi) then
satpcp(nsdata) = data_p(6,n)
satcnv(nsdata) = data_p(7,n)
satclw(nsdata) = data_p(8,n)
satcli(nsdata) = data_p(9,n)
pcpnum(nsdata) = data_p(10,n)
pcpsdv(nsdata) = rmiss
endif
!
! Adjust obs error based on precipitation rate.
!
! Bill Olson TMI error model.
if (ksatid(nsdata)==211) then
if (sea) then
a0=0.137; a1=0.118 !tmi,sea
else
a0=0.335; a1=0.045 !tmi, land
endif
! Bob Kuligowski SSMI error model
elseif (ksatid(nsdata)==264) then
if (sea) then
a0=0.3835; a1=0.1699 !ssmi, sea (n=50 bin) (scattering)
else
a0=0.3148; a1=0.1781 !ssmi, land (n=50 bin)
endif
else
a0=0.0; a1=0.0
endif
term = log(1.+satpcp(nsdata))
errlog = a0 + a1*term
obserr(nsdata) = errlog
error = varchp(kx) + obserr(nsdata)
errf(nsdata) = error
varinv(nsdata) = 1./(error*error)
!
! If negative precipitation, reset to 0.0
if (satpcp(nsdata) < 0.) satpcp(nsdata)=0.0
!
! Save data in diagnostic arrays.
idiagbuf(1,nsdata) = ksatid(nsdata)
idiagbuf(2,nsdata) = lndsea(nsdata)
diagbuf(1,nsdata) = slats(nsdata)*rad2dg
diagbuf(2,nsdata) = slons(nsdata)*rad2dg
diagbuf(3,nsdata) = stfact(nsdata)
diagbuf(4,nsdata) = pcpnum(nsdata)
diagbuf(5,nsdata) = satpcp(nsdata)
diagbuf(6,nsdata) = pcpsdv(nsdata)
diagbuf(7,nsdata) = satcnv(nsdata)
diagbuf(8,nsdata) = satclw(nsdata)
diagbuf(9,nsdata) = satcli(nsdata)
diagbuf(10,nsdata) = varinv(nsdata)
diagbuf(11,nsdata) = errf(nsdata)
! End of loop over obs.
end do
aivals(2,isz) = aivals(2,isz) + nsdata
!
!
!
!
!*****
! IDENTIFY OBSERVATIONS TO PROCESS
!*****
!
!
! Get skin temperature at observation location.
call gdcrdp(slats,nsdata,rlats,nlath2)
call gdcrdp(slons,nsdata,rlons,nlon)
call intrp2(nlath2,sfct,ts0,slats,slons,lat1,lon1,nlon,nsdata,mype)
!
! Do some prelimiary qc of the data. Check for points covered
! with ice/snow, land points, and coastal points.
!
do n = 1,nsdata
isflg(n)=0
icoast =0
land(n) =.false.
coast(n)=.false.
ice(n) =.false.
!
! Check for ocean temperatures less the freezing (sea ice)
! or probable snow over land
isflg(n)=0
if ( (ts0(n) < 272.16) .and. (abs(cenlat(n)) >= 40.) ) isflg(n)=1
!
! Determine local (i,j) indices of analysis gridpoints surrounding
! the observation. Also determine bilinear interpolation weights
! from analysis gridpoints to observation location.
! respect to the surrounding analysis grid box.
iy1 = slats(n)
ix1 = slons(n)
dx = slons(n)-ix1
dy = slats(n)-iy1
iy1=min(max(1,iy1),nlath2); ix1=min(max(0,ix1),nlon)
dx1 = 1.0-dx
dy1 = 1.0-dy
iy1 = max(1,min(iy1,nlath2))
iy = iy1-istart(mype1)+2
ix = ix1-jstart(mype1)+2
if (ix<1) then
ix1 = ix1+nlon
ix = ix1-jstart(mype1)+2
end if
if (ix>lon1+1) then
ix1 = ix1-nlon
ix = ix1-jstart(mype1)+2
end if
ixp = ix+1
iyp = iy+1
if (iy1==nlath2) iyp=iy
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
! Save (i,j) indices and interpolation weights.
iloc(n,1)=iy; iloc(n,2)=iy; iloc(n,3)=iyp; iloc(n,4)=iyp
jloc(n,1)=ix; jloc(n,2)=ixp; jloc(n,3)=ix; jloc(n,4)=ixp
wgtloc(n,1)=wgt00; wgtloc(n,2)=wgt10
wgtloc(n,3)=wgt01; wgtloc(n,4)=wgt11
ijloc(n,1)=j11; ijloc(n,2)=j21; ijloc(n,3)=j12; ijloc(n,4)=j22
! Check for coastal observations. Coastal is defined with
! respect to the surrounding analysis grid box.
if (lndsea(n) /= 1)then
if (isli(iy,ix) .eq. 1 .or. isli(iy,ixp) .eq. 1 .or.&
isli(iyp,ix) .eq. 1 .or. isli(iyp,ixp) .eq. 1)&
lndsea(n) = 1
end if
!
! Check for snow of sea ice in surrounding grid box
if(isflg(n) /= 1)then
if((sno(iy,ix) > 1. .and. isli(iy,ix) /= 0) .or.&
isli(iy,ix) .eq. 2 ) isflg(n)=1
if((sno(iy,ixp) > 1. .and. isli(iy,ixp) /= 0) .or.&
isli(iy,ixp) .eq. 2 ) isflg(n)=1
if((sno(iyp,ix) > 1. .and. isli(iyp,ix) /= 0) .or.&
isli(iyp,ix) .eq. 2 ) isflg(n)=1
if((sno(iyp,ixp) > 1. .and. isli(iyp,ixp) /= 0) .or.&
isli(iyp,ixp) .eq. 2 ) isflg(n)=1
end if
!
! Check for coastal (transition) points
if (isli(iy,ix) /= isli(iy,ixp) .or.&
isli(iy,ix) /= isli(iyp,ix) .or.&
isli(iy,ix) /= isli(iyp,ixp)) coast(n) = .true.
! Set logical flags
land(n) = lndsea(n) == 1
ice(n) = isflg(n) == 1
if (coast(n)) icoast = 1
! Save data in diagnostic arrays
idiagbuf(3,n)= lndsea(n)
idiagbuf(4,n)= isflg(n)
idiagbuf(5,n)= icoast
!
! End loop over obs
end do
!
!
!
!*****
! CALL Precipitation (FORWARD) MODEL
!*****
!
!
! Call precipitation model. Set ajm forcing
! for precipitation rate to unity.
! Initialize variables and arrays
im = nsdata
ix = jppfp
km = nsig
do k = 1,nsig
del(k) = sigi(k) - sigi(k+1)
slk(k) = sigl(k)**rdocp
end do
do i = 1,jppfp
do k = 1,nsig
ad_t0(k,i) = 0.0
ad_q0(k,i) = 0.0
ad_u0(k,i) = 0.0
ad_v0(k,i) = 0.0
ad_div0(k,i) = 0.0
ad_cwm0(k,i) = 0.0
end do
rn1(i) = 0.0
pcpsas(i) = 0.0
pcplrg(i) = 0.0
end do
! Loop over four surrounding analysis gridpoints
do ipt=1,4
! Loop over observation points.
do n=1,nsdata
! Set (i,j) index for ipt-th point for observation n.
i=iloc(n,ipt); j=jloc(n,ipt)
! Load arrays used by forward/adjoint model
xkt2(n) = xkt2d(i,j)
ii = slats(n)+0.5
rbs0(n) = rbs2(ii)
rmask0(n) = isli(i,j)
! Set weights for linear time interpolation.
if (deltat>0.0) then
tfo=stfact(n)
tf1=-min(max(oneneg,tfo),zero)
tf2=min(max(zero,1.0-abs(tfo)),one)
tf3=min(max(zero,tfo),one)
else
tf1=0.0; tf2=1.0; tf3=0.0
end if
! Perform interpolation.
pterm = tf2*ggrid_g31(3,i,j,nsig1) + tf1*ggrid_g32(3,i,j,nsig1,1) + &
tf3*ggrid_g32(3,i,j,nsig1,2)
!ibm* prefetch_for_store(psexp0(n),plon0(n),plat0(n))
psexp0(n) = exp(pterm)
plon0(n) = tf2*ggrid_g31(2,i,j,ns2) + tf1*ggrid_g32(2,i,j,ns2,1) + &
tf3*ggrid_g32(2,i,j,ns2,2)
plat0(n) = tf2*ggrid_g31(1,i,j,ns2) + tf1*ggrid_g32(1,i,j,ns2,1) + &
tf3*ggrid_g32(1,i,j,ns2,2)
!ibm* prefetch_for_store(k1)
k1=nsig1
do k=1,nsig
k1=k1+1
!ibm* prefetch_for_store(t0(k,n),q0(k,n))
t0(k,n) = tf2*ggrid_g31(3,i,j,k1) + tf1*ggrid_g32(3,i,j,k1,1) + &
tf3*ggrid_g32(3,i,j,k1,2)
q0(k,n) = tf2*ggrid_g31(1,i,j,k1) + tf1*ggrid_g32(1,i,j,k1,1) + &
tf3*ggrid_g32(1,i,j,k1,2)
end do
do k=1,nsig
!ibm* prefetch_for_store(u0(k,n),v0(k,n),div0(k,n),cwm0(k,n),tlon0(k,n),tlat0(k,n))
u0(k,n) = tf2*ggrid_g31(1,i,j,k) + tf1*ggrid_g32(1,i,j,k,1) + &
tf3*ggrid_g32(1,i,j,k,2)
v0(k,n) = tf2*ggrid_g31(2,i,j,k) + tf1*ggrid_g32(2,i,j,k,1) + &
tf3*ggrid_g32(2,i,j,k,2)
div0(k,n) = tf2*rdivs(i,j,k,1) + tf1*rdivs(i,j,k,2) + &
tf3*rdivs(i,j,k,3)
cwm0(k,n) = tf2*rcwm(i,j,k,1) + tf1*rcwm(i,j,k,2) + &
tf3*rcwm(i,j,k,3)
tlat0(k,n) = tf2*rtlat(i,j,k,1) + tf1*rtlat(i,j,k,2) + &
tf3*rtlat(i,j,k,3)
tlon0(k,n) = tf2*rtlon(i,j,k,1) + tf1*rtlon(i,j,k,2) + &
tf3*rtlon(i,j,k,3)
qlat0(k,n) = tf2*rqlat(i,j,k,1) + tf1*rqlat(i,j,k,2) + &
tf3*rqlat(i,j,k,3)
qlon0(k,n) = tf2*rqlon(i,j,k,1) + tf1*rqlon(i,j,k,2) + &
tf3*rqlon(i,j,k,3)
end do
end do
! Initialize arrays below prior to each call to forward/adjoint model
do i = 1,jppfp
do k = 1,nsig
ad_t4(k,i) = 0.0
ad_q4(k,i) = 0.0
ad_u4(k,i) = 0.0
ad_v4(k,i) = 0.0
ad_div4(k,i) = 0.0
ad_cwm4(k,i) = 0.0
ad_t1(k,i) = 0.0
ad_q1(k,i) = 0.0
ad_u1(k,i) = 0.0
ad_v1(k,i) = 0.0
ad_cwm1(k,i) = 0.0
t1(k,i) = 0.0
q1(k,i) = 0.0
u1(k,i) = 0.0
v1(k,i) = 0.0
cwm1(k,i) = 0.0
end do
rn4(i) = 0.0
ad_rn1(i) = 1.0
pcpsas4(i) = 0.0
pcplrg4(i) = 0.0
end do
! Call forward/adjoint model
call ajmpcp(im,ix,km,jcap,dtp,del,sigl,slk,rbs0, &
rmask0,nlon,xkt2,ncloud,frain,rmmhr,&
psexp0,plat0,plon0,&
t0,q0,u0,v0,div0,cwm0,tlat0,tlon0,qlat0,qlon0,&
tsas,qsas,pcpsas4,cldwrk,kbcon,ktcon,jmin,kuo, &
tlrg,qlrg,pcplrg4, &
t1,q1,cwm1,u1,v1,rn4,&
ad_t4,ad_q4,ad_cwm4,ad_u4,ad_v4,ad_div4,&
ad_t1,ad_q1,ad_cwm1,ad_u1,ad_v1,ad_rn1,mype)
! Based on bilinear interpolation weights, accumulate rain rates
! and sensitivity over four surrounding analysis gridpoints.
do n=1,nsdata
wgtij=wgtloc(n,ipt)
pcpsas(n) = pcpsas(n) + wgtij*pcpsas4(n)
pcplrg(n) = pcplrg(n) + wgtij*pcplrg4(n)
rn1(n) = rn1(n) + wgtij*rn4(n)
do k=1,nsig
ad_t0(k,n) = ad_t0(k,n) + wgtij*ad_t4(k,n)
ad_q0(k,n) = ad_q0(k,n) + wgtij*ad_q4(k,n)
ad_cwm0(k,n) = ad_cwm0(k,n) + wgtij*ad_cwm4(k,n)
ad_u0(k,n) = ad_u0(k,n) + wgtij*ad_u4(k,n)
ad_v0(k,n) = ad_v0(k,n) + wgtij*ad_v4(k,n)
ad_div0(k,n) = ad_div0(k,n) + wgtij*ad_div4(k,n)
end do
end do
! End of loop over analysis gridpoints
end do
! Load model simulated rain rates into work array.
! NOTE: simulated values converted to units of mm/hr
! in ajmpcp. obs already in mm/hr
do n = 1,nsdata
pcpnbc(n) = rn1(n)
end do
!!test
! Adjust error based on estimated forward model error.
do n=1,nsdata
if (ksatid(n)==211) then
if (land(n)) then
a0=0.1244544; a1=0.3883144 !tmi, land
else
a0=0.1764123; a1=0.3948465 !tmi, sea
endif
elseif (ksatid(n)==264) then
if (land(n)) then
a0=0.1282869; a1=0.4108055 !ssmi, land
else
a0=0.0917045; a1=0.4757942 !ssmi, sea
endif
else
a0=0.0; a1=0.0
endif
term = log(1.+pcpnbc(n))
simerr = a0 + a1*term
simerr = max(zero,simerr)
obserr_wgt = varchp(istype(n))
simerr_wgt = 1.-obserr_wgt
error = max(obserr(n),obserr_wgt*obserr(n) + simerr_wgt*simerr)
errf(n) = error
varinv(n)=1./(error*error)
end do
!!test
! Save cloud base, downdraft origin, cloud top, cloud work function
! for diagnostics
do n = 1,nsdata
idiagbuf(6,n) = kbcon(n)
idiagbuf(7,n) = ktcon(n)
idiagbuf(8,n) = jmin(n)
diagbuf(12,n) = rmiss
diagbuf(13,n) = rmiss
diagbuf(14,n) = rmiss
if (kbcon(n)>=1 .and. kbcon(n)<=nsig) diagbuf(12,n) = 10.*psexp0(n) * sigl(kbcon(n))
if (ktcon(n)>=1 .and. ktcon(n)<=nsig) diagbuf(13,n) = 10.*psexp0(n) * sigl(ktcon(n))
if (jmin(n)>=1 .and. jmin(n)<=nsig) diagbuf(14,n) = 10.*psexp0(n) * sigl(jmin(n))
diagbuf(15,n) = cldwrk(n)
diagbuf(23,n) = xkt2(n)
end do
! If magnitude of sensitivity vector is zero, set inverse observation
! error to zero since we don't have gradient information for this
! observation. For now, ommit divergence sensitivity since it is
! not used in pcgsoi. Likewise, if the gradient is too large, do not
! use this observation (this check rarely catches any points).
rterm1=1./float(nsig)
rterm2=1./float(nsig*(nsig-1))
do n = 1,nsdata
dpcpmag=0.0
do k=1,nsig
dpcpmag(1) = dpcpmag(1) + ad_t0(k,n)
dpcpmag(1) = dpcpmag(1) + ad_q0(k,n)
dpcpmag(1) = dpcpmag(1) + ad_u0(k,n)
dpcpmag(1) = dpcpmag(1) + ad_v0(k,n)
dpcpmag(1) = dpcpmag(1) + ad_cwm0(k,n)
dpcpmag(2) = dpcpmag(2) + ad_t0(k,n)*ad_t0(k,n)
dpcpmag(2) = dpcpmag(2) + ad_q0(k,n)*ad_q0(k,n)
dpcpmag(2) = dpcpmag(2) + ad_u0(k,n)*ad_u0(k,n)
dpcpmag(2) = dpcpmag(2) + ad_v0(k,n)*ad_v0(k,n)
dpcpmag(2) = dpcpmag(2) + ad_cwm0(k,n)*ad_cwm0(k,n)
end do
avg = dpcpmag(1)*rterm1
term= (nsig*dpcpmag(2)-dpcpmag(1)**2)*rterm2
sdv = 0.0
if (term>0) sdv=sqrt(term)
dpcpmag(2)=sqrt(dpcpmag(2))
if (dpcpmag(2)<1.e-10) then
aivals(26,isz) = aivals(26,isz) + weight_p(n)
varinv(n)=0.0
endif
if (dpcpmag(2)>1.e4) then
aivals(27,isz) = aivals(27,isz) + weight_p(n)
varinv(n)=0.0
elseif (abs(avg)>0.) then
if (abs(sdv/avg)>25.) then
aivals(28,isz) = aivals(28,isz) + weight_p(n)
varinv(n) = 0.0
endif
endif
if (varinv(n)>1.e-6) aivals(3,isz) = aivals(3,isz) + weight_p(n)
end do
!
!
!
!
!
!*****
! COMPUTE AND APPLY BIAS CORRECTION TO SIMULATED VALUES
!*****
!
! Prepare for application of bias correction to simulated values.
pred=0.0
!
! Construct predictors for bias correction.
do n = 1,nsdata
diagbuf(16,n) = pcpsas(n)
diagbuf(17,n) = pcpnbc(n)
pred(1,n) = 0.0
pred(2,n) = 0.0
pred(3,n) = 0.0
pred(4,n) = 0.0
pred(5,n) = 0.0
pred(6,n) = 0.0
end do
!
!
! Apply bias correction to simulated values.
do n = 1,nsdata
sum = 0.0
do j = 1,npredp
sum = sum + predxp(istype(n),j)*pred(j,n)* &
ibias(istype(n))
end do
pcpbc(n) = pcpnbc(n) + sum
diagbuf(18,n) = pcpbc(n)
end do
!
!
!
!
!***
! QC OBSERVATIONS BASED ON VARIOUS CRITERIA
!***
!
! QC SSMI, TMI, and AMSUA precipitation.
if (ssmi .or. tmi .or. amsua) then
!
! Loop over observations.
do n = 1,nsdata
efact = 1.
vfact = 0.25
!
! Reduce qc bounds in extratropics. Between 50
! and 80 degrees latitude linearly reduce error
! bound. Poleward of 80 degrees latitude do not use
if ( abs(cenlat(n)) > 50. ) then
aivals(4,isz) = aivals(4,isz) + weight_p(n)
term = (80.-abs(cenlat(n)))/30.
term = max(zero,term)
term = min(term,one)
if (abs(cenlat(n)) > 80.) term=0.0
efact = term*efact
endif
!
! If land obs, reduce error bounds and weight
if (land(n)) then
efact = 0.7*efact
vfact = 0.7*vfact
aivals(5,isz) = aivals(5,isz) + weight_p(n)
end if
!
! Toss obs if over snow or ice
if (ice(n)) then
efact = 0.0
vfact = 0.0
aivals(6,isz) = aivals(6,isz) + weight_p(n)
end if
!
! If coastal obs, reduce error bounds and weight
if (coast(n)) then
efact = 0.5*efact
vfact = 0.5*vfact
aivals(7,isz) = aivals(7,isz) + weight_p(n)
end if
!
! Reduce weight and tighten obs-ges bound as function of obs time.
term = cos(halfpi*stfact(n))
term = term**3
term = max(zero,min(term,one))
efact = term*efact
vfact = term*vfact
!
! Generate q.c. bounds and modified variances.
errf(n) = efact*errf(n)
varinv(n) = vfact*varinv(n)
diagbuf(19,n) = varinv(n)
!
! End of loop over observations.
end do
!
! End of SSMI, TMI, and AMSUA qc block
endif
!
!
! Apply gross check to good pcp obs.
do n = 1,nsdata
xlo = -3.*errf(n)
xhi = +3.*errf(n)
if (varinv(n) > 1.e-6) then
drad = log(1.+satpcp(n)) - log(1.+pcpbc(n))
if (drad < xlo .or. drad > xhi) then
varinv(n) = 0.0
aivals(12,isz) = aivals(12,isz) + weight_p(n)
endif
endif
diagbuf(20,n) = varinv(n)
diagbuf(21,n) = errf(n)
end do
!
!
!***
! CONSTRUCT SENSITIVITY VECTORS
!***
!
! Initialize dpcp/dT, dq, du, dv, dd sensitivity vectors.
do n = 1,nsdata
var(n) = 0.
obs(n) = 0.
ges(n) = 0.
end do
do k = 1,nsig
dpcpdt(k) = 0.
dpcpdq(k) = 0.
dpcpdu(k) = 0.
dpcpdv(k) = 0.
dpcpddiv(k) = 0.
dpcpdcwm(k) = 0.
end do
! Load information pcgsoi needs into output arrays.
do n = 1,nsdata
! Load values for "good" obs.
if (varinv(n) > 1.e-6) then
var(n) = varinv(n)
obs(n) = satpcp(n)
ges(n) = pcpbc(n)
! Load sensitivity vector.
do i = 1,nsig
dpcpdt(i) = ad_t0(i,n)
dpcpdq(i) = ad_q0(i,n)
dpcpdu(i) = ad_u0(i,n)
dpcpdv(i) = ad_v0(i,n)
dpcpddiv(i)= ad_div0(i,n)
dpcpdcwm(i)= ad_cwm0(i,n)
end do
end if
!
! Write diagnostics to output file.
if (jiter==1 .and. idiag_pcp) then
! if (varinv(n)>1.e-6) then
diagbuf(22,n) = weight_p(n)
write(4) mype,(idiagbuf(m,n),m=1,iint),&
(diagbuf(m,n),m=1,ireal)
! endif
endif
! Optionally, set parameter to turn off data.
if (iusep(istype(n))<1) then
do i=1,nsig
dpcpdt(i) = 0.0
dpcpdq(i) = 0.0
dpcpdu(i) = 0.0
dpcpdv(i) = 0.0
dpcpddiv(i)= 0.0
dpcpdcwm(i)= 0.0
end do
var(n) = 0.0
obs(n) = 0.0
ges(n) = 0.0
endif
! Write data needed by pcgsoi to file
if (varinv(n)>1.e-6) then
n1 = n1+1
write(lunout) (ijloc(n,m),m=1,4),(wgtloc(n,m),m=1,4),weight_p(n)
write(lunout) istype(n),var(n),(pred(m,n),m=1,npredp),&
obs(n),ges(n),&
(dpcpdt(m),m=1,nsig),(dpcpdq(m),m=1,nsig),&
(dpcpdu(m),m=1,nsig),(dpcpdv(m),m=1,nsig),&
(dpcpdcwm(m),m=1,nsig)
! For now, do not include divergence sensitivity
! (dpcpddiv(m,n),m=1,nsig)
! End of pcgsoi write block
end if
!
! End of loop over profiles
end do
!
!
!
!******
! DIAGNOSTIC CALCULATIONS.
!******
!
!
! Accumulate data for diagnostic print.
do n = 1,nsdata
if (varinv(n) > 1.e-6) then
drad = log(1.+satpcp(n)) - log(1.+pcpbc(n))
varrad = drad*varinv(n)
error = sqrt(1./varinv(n))
aivals(11,isz) = aivals(11,isz) + weight_p(n)
aivals(13,isz) = aivals(13,isz) + drad*weight_p(n) !bias
aivals(14,isz) = aivals(14,isz) + drad*drad*weight_p(n) !bias**2
aivals(15,isz) = aivals(15,isz) + drad*varrad*weight_p(n) !pen
aivals(16,isz) = aivals(16,isz) + error*weight_p(n) !obs error
aivals(21,isz) = aivals(21,isz) + satpcp(n)*weight_p(n)
aivals(22,isz) = aivals(22,isz) + pcpnbc(n)*weight_p(n)
aivals(23,isz) = aivals(23,isz) + pcpbc(n) *weight_p(n)
aivals(24,isz) = aivals(24,isz) + pcpsas(n)*weight_p(n)
aivals(25,isz) = aivals(25,isz) + weight_p(n)
else
aivals(31,isz) = aivals(31,isz) + satpcp(n)*weight_p(n)
aivals(32,isz) = aivals(32,isz) + pcpnbc(n)*weight_p(n)
aivals(33,isz) = aivals(33,isz) + pcpbc(n) *weight_p(n)
aivals(34,isz) = aivals(34,isz) + pcpsas(n)*weight_p(n)
aivals(35,isz) = aivals(35,isz) + weight_p(n)
endif
end do
!
! End of nblk loop
end do
!
!
! Jump here when there is no more data to process for current satellite
134 continue
deallocate(data_p)
deallocate(weight_p)
135 continue
npcp1d(is)=n1
!
!
!*****
! END MAIN LOOP OVER SATELLITES
!*****
end do
npcpdat_s=0
do k=1,nsat+5
npcpdat_s=npcpdat_s+npcp1d(k)
end do
! Close unit to diagnostic file.
if (jiter==1 .and. idiag_pcp) close(4)
!
! Collect statistics
call mpi_reduce(aivals,aivals1,36*nsat,mpi_rtype,mpi_sum,mype_pcp,mpi_comm_world,ierror)
call mpi_reduce(kidsat,kidsat1, nsat,mpi_integer,mpi_max,mype_pcp,mpi_comm_world,ierror)
! Compute and print statistics for current satellite.
if(mype==mype_pcp) then
penalty_all=0.0
do is=1,nsat
obstype=sattype(is)
if(obstype .eq. 'pcp_ssm/i')ssmi=.true.
if(obstype .eq. 'pcp_tmi')tmi=.true.
if(obstype .eq. 'pcp_amsua')amsua=.true.
inumw = nint(aivals1(1,is))
inum = nint(aivals1(2,is))
itotal = nint(aivals1(3,is))
ireduce = nint(aivals1(4,is))
iland = nint(aivals1(5,is))
isnoice = nint(aivals1(6,is))
icoast = nint(aivals1(7,is))
iclip = nint(aivals1(12,is))
igradsml = nint(aivals1(26,is))
igradbig = nint(aivals1(27,is))
iratio = nint(aivals1(28,is))
rterm=0.0
if (aivals1(25,is) > 0.) rterm=1./aivals1(25,is)
pcpobs0 = aivals1(21,is)*rterm
pcpnbc0 = aivals1(22,is)*rterm
pcpbc0 = aivals1(23,is)*rterm
pcpsas0 = aivals1(24,is)*rterm
ikeep = nint(aivals1(25,is))
rterm=0.0
if (aivals1(35,is) > 0.) rterm=1./aivals1(35,is)
pcpobs1 = aivals1(31,is)*rterm
pcpnbc1 = aivals1(32,is)*rterm
pcpbc1 = aivals1(33,is)*rterm
pcpsas1 = aivals1(34,is)*rterm
itoss = nint(aivals1(35,is))
! Print qc counts.
if (inum > 0) then
write(iout_pcp,*)' '
write(iout_pcp,2000) 'type',' ','num','numw',&
'itotal','ireduce','iland','isnoice','icoast',&
'grad<1e-10','grad>1e4','ratio>3'
write(iout_pcp,2010) obstype,inum,inumw,&
itotal,ireduce,iland,isnoice,icoast,&
igradsml,igradbig,iratio
write(iout_pcp,2012) obstype,aivals1(15,is)
write(iout_pcp,2000) 'qc-flag','satpcp ','pcpnbc ',&
'pcpsas ','pcplrg ','pcpbc ',' count '
write(iout_pcp,2013)'keep',pcpobs0,pcpnbc0,pcpsas0,&
pcpnbc0-pcpsas0,pcpbc0,ikeep
write(iout_pcp,2013)'toss',pcpobs1,pcpnbc1,pcpsas1,&
pcpnbc1-pcpsas1,pcpbc1,itoss
end if
2000 format(15(a10,1x))
2010 format(1x,A10,' counts ',2x,12(i10,1x))
2011 format(2012 format(1x,A10,' penalty ',g18.12)
2013 format(a10,1x,5(e10.4,1x),i6)
2014 format(
! Print penalties
penalty_all=penalty_all+aivals1(15,is)
end do
write(iout_pcp,*)' '
write(iout_pcp,*)'pcp total penalty_all=',penalty_all
! Print counts, bias, rms, stndev as a function of channel.
do is = 1,nsat
isum = nint(aivals1(11,is))
if (isum > 0) then
rpen(is) = aivals1(15,is)
rsum = 1./float(isum)
icerr = nint(aivals1(12,is))
do j=13,16
aivals1(j,is)=aivals1(j,is)*rsum
end do
aivals1(14,is) = sqrt(aivals1(14,is))
cpen(is) = aivals1(15,is)
if (isum > 1) then
sdv = sqrt(&
aivals1(14,is)*aivals1(14,is)- &
aivals1(13,is)*aivals1(13,is))
else
sdv = 0.0
endif
write(iout_pcp,1102) kidsat1(is),sattype(is),isum,icerr,&
aivals1(16,is),aivals1(13,is),aivals1(15,is),aivals1(14,is),sdv
1102 format(1x,i3,1x,a10,1x,2(i6,1x),6(g18.12,1x))
else
rpen(is) = 0.0
cpen(is) = 0.0
end if
end do
write(iout_pcp,1109)
do i=1,nsat
if (aivals1(11,i)>0.) &
write(iout_pcp,1115) jiter,kidsat1(i),aivals1(11,i),&
aivals1(13,i),aivals1(14,i),rpen(i),cpen(i)
end do
1109 format(t6,'it',t14,'sat',t21,'num use',t35,'bias',t45,'rms',t55,'penalty',t69,'cpen')
1115 format(t2,'o-g',1x,i2.2,1x,'pcp',2x,i3,2x,f9.0,2x,g10.4,2x,g10.4,2x,f9.4,2x,f12.6)
! End of diagnostic print block.
close(iout_pcp)
end if
!
! Deallocate arrays
deallocate(tsas,qsas,tlrg,qlrg)
deallocate(t0,q0,cwm0,u0,v0,div0,t1,q1,cwm1,u1,v1,tlat0,tlon0,qlat0,qlon0)
deallocate(ad_t0,ad_q0,ad_cwm0,ad_u0,ad_v0,ad_div0,ad_t1,ad_q1,ad_cwm1,&
ad_u1,ad_v1,ad_t4,ad_q4,ad_cwm4,ad_u4,ad_v4,ad_div4)
! End of routine
return
end subroutine setuppcp