setupt includes forward operators for virtual, sensible temperature, T2m and PBL pseudo surface observations
General workflow:
QC methods from setupt generally can be divided into methods that use only observation data (+ some configs):
and methods that require state information:
Note: gross checks potentially can use ensemble spread information.
type(t_obs_space) :: obs type(t_obs_operator) :: obs_op type(t_linear_obs_operator) :: lin_obs_op type(t_obs_error) :: obserr type(obs_data) :: hx type(fields_at_locations) :: int_state call obs%create(???) ! reading observations call obserr%create(t_obs_err_config, obs) ! initializing observation error (TODO: needs yobs on input, ???) call obserr%quality_control() ! QC that uses only observation data: handle multiple reports; observations out of time window (3, 4) call int_state%create(obs, geom, vars) ! interpolate fields to observation locations? (TODO: need fields on input!) call obs_op%create(t_hx_config) ! initializing observation operator call lin_obs_op%create(obs_op, int_state, bias_pred, t_linhx_config) ! initializing linear observation operator hx = obs_op%hoper(int_state, bias_pred) ! calculate innovations; do surface temp correction and bias correction for aircraft data (1, 5, 6) call obserr%quality_control(hx) ! QC that uses state info: buddy check, gross check (2, 7, 8?) |
AOD bufr created from HDF files, different HDF formats for MODIS, VIIRS. Also ABI GOES-R possible to implement
General workflow
subroutine setuprad(lunin,mype,aivals,stats,nchanl,nreal,nobs,& obstype,isis,is,rad_diagsave,init_pass,last_pass)
!$$$ subprogram documentation block
! . . . .
! subprogram: setuprad compute rhs of oi equation for radiances
! prgmmr: derber org: np23 date: 1995-07-06
!
! abstract: read in data, first guess, and obtain rhs of oi equation
! for radiances.
!
! program history log:
!
! input argument list:
! lunin - unit from which to read radiance (brightness temperature, tb) obs
! mype - mpi task id
! nchanl - number of channels per obs
! nreal - number of pieces of non-tb information per obs
! nobs - number of tb observations to process
! obstype - type of tb observation
! isis - sensor/instrument/satellite id ex.amsua_n15
! is - integer counter for number of observation types to process
! rad_diagsave - logical to switch on diagnostic output (.false.=no output)
! channelinfo - structure containing satellite sensor information
!
! output argument list:
! aivals - array holding sums for various statistics as a function of obs type
! stats - array holding sums for various statistics as a function of channel
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use mpeu_util, only: die,perr,getindex
use kinds, only: r_kind,r_single,i_kind
use crtm_spccoeff, only: sc
use radinfo, only: nuchan,tlapmean,predx,cbias,ermax_rad,tzr_qc,&
npred,jpch_rad,varch,varch_cld,iuse_rad,icld_det,nusis,fbias,retrieval,b_rad,pg_rad,&
air_rad,ang_rad,adp_anglebc,angord,ssmis_precond,emiss_bc,upd_pred, &
passive_bc,ostats,rstats,newpc4pred,radjacnames,radjacindxs,nsigradjac
use gsi_nstcouplermod, only: nstinfo
use read_diag, only: get_radiag,ireal_radiag,ipchan_radiag
use guess_grids, only: sfcmod_gfs,sfcmod_mm5,comp_fact10
use m_prad, only: radheadm
use m_obsdiags, only: radhead
use obsmod, only: ianldate,ndat,mype_diaghdr,nchan_total, &
dplat,dtbduv_on,&
i_rad_ob_type,obsdiags,obsptr,lobsdiagsave,nobskeep,lobsdiag_allocated,&
dirname,time_offset,lwrite_predterms,lwrite_peakwt,reduce_diag
use m_obsNode, only: obsNode
use m_radNode, only: radNode, radNode_typecast
use m_obsLList, only: obsLList_appendNode
use m_obsLList, only: obsLList_tailNode
use obsmod, only: obs_diag,luse_obsdiag,dval_use
use gsi_4dvar, only: nobs_bins,hr_obsbin,l4dvar
use gridmod, only: nsig,regional,get_ij
use satthin, only: super_val1
use constants, only: quarter,half,tiny_r_kind,zero,one,deg2rad,rad2deg,one_tenth, &
two,three,cg_term,wgtlim,r100,r10,r0_01,r_missing
use jfunc, only: jiter,miter,jiterstart
use sst_retrieval, only: setup_sst_retrieval,avhrr_sst_retrieval,&
finish_sst_retrieval,spline_cub
use m_dtime, only: dtime_setup, dtime_check, dtime_show
use crtm_interface, only: init_crtm,call_crtm,destroy_crtm,sensorindex,surface, &
itime,ilon,ilat,ilzen_ang,ilazi_ang,iscan_ang,iscan_pos,iszen_ang,isazi_ang, &
ifrac_sea,ifrac_lnd,ifrac_ice,ifrac_sno,itsavg, &
izz,idomsfc,isfcr,iff10,ilone,ilate, &
isst_hires,isst_navy,idata_type,iclr_sky,itref,idtw,idtc,itz_tr
use clw_mod, only: calc_clw, ret_amsua
use qcmod, only: qc_ssmi,qc_seviri,qc_ssu,qc_avhrr,qc_goesimg,qc_msu,qc_irsnd,qc_amsua,qc_mhs,qc_atms
use qcmod, only: igood_qc,ifail_gross_qc,ifail_interchan_qc,ifail_crtm_qc,ifail_satinfo_qc,qc_noirjaco3,ifail_cloud_qc
use qcmod, only: qc_gmi,qc_saphir,qc_amsr2
use qcmod, only: setup_tzr_qc,ifail_scanedge_qc,ifail_outside_range
use gsi_metguess_mod, only: gsi_metguess_get
use oneobmod, only: lsingleradob,obchan,oblat,oblon,oneob_type
use radinfo, only: radinfo_adjust_jacobian,radinfo_get_rsqrtinv
implicit none
! Declare passed variables
logical ,intent(in ) :: rad_diagsave
character(10) ,intent(in ) :: obstype
character(20) ,intent(in ) :: isis
integer(i_kind) ,intent(in ) :: lunin,mype,nchanl,nreal,nobs,is
real(r_kind),dimension(40,ndat) ,intent(inout) :: aivals
real(r_kind),dimension(7,jpch_rad),intent(inout) :: stats
logical ,intent(in ) :: init_pass,last_pass ! state of "setup" processing
|
Initialization:
Initialize variables and constants.
Initialize logical flags for satellite platform
Determine whether or not cloud radiance: lcw4crtm
Parameters for the observation error model when lcw4crtm=true
Initialize channel related information: predchan; iuse_rad;passive_bc;...
Set error instrument channels
Logic to turn off print of reading coefficients if not first interation or not mype_diaghdr or not init_pass
Initialize radiative transfer and pointers to values in data_s
call init_crtm(init_pass,iwrmype,mype,nchanl,isis,obstype)
Get indexes of variables in jacobian to handle exceptions down below
Initialize ozone jacobian flags to .false. (retain ozone jacobian)
If SSM/I, check for non-use of 85GHz channel, for QC workaround set no85GHz true if any 85GHz is not used, and other freq channel is used no85GHz = .false.
Setup for diagnostic
read(lunin) data_s,luse,ioid
! These variables are initialized in init_crtm ! isatid = 1 ! index of satellite id ! itime = 2 ! index of analysis relative obs time ! ilon = 3 ! index of grid relative obs location (x) ! ilat = 4 ! index of grid relative obs location (y) ! ilzen_ang = 5 ! index of local (satellite) zenith angle (radians) ! ilazi_ang = 6 ! index of local (satellite) azimuth angle (radians) ! iscan_ang = 7 ! index of scan (look) angle (radians) ! iscan_pos = 8 ! index of integer scan position ! iszen_ang = 9 ! index of solar zenith angle (degrees) ! isazi_ang = 10 ! index of solar azimuth angle (degrees) ! ifrac_sea = 11 ! index of ocean percentage ! ifrac_lnd = 12 ! index of land percentage ! ifrac_ice = 13 ! index of ice percentage ! ifrac_sno = 14 ! index of snow percentage ! its_sea = 15 ! index of ocean temperature ! its_lnd = 16 ! index of land temperature ! its_ice = 17 ! index of ice temperature ! its_sno = 18 ! index of snow temperature ! itsavg = 19 ! index of average temperature ! ivty = 20 ! index of vegetation type ! ivfr = 21 ! index of vegetation fraction ! isty = 22 ! index of soil type ! istp = 23 ! index of soil temperature ! ism = 24 ! index of soil moisture ! isn = 25 ! index of snow depth ! izz = 26 ! index of surface height ! idomsfc = 27 ! index of dominate surface type ! isfcr = 28 ! index of surface roughness ! iff10 = 29 ! index of ten meter wind factor ! ilone = 30 ! index of earth relative longitude (degrees) ! ilate = 31 ! index of earth relative latitude (degrees) ! itref = 34/36 ! index of foundation temperature: Tr ! idtw = 35/37 ! index of diurnal warming: d(Tw) at depth zob ! idtc = 36/38 ! index of sub-layer cooling: d(Tc) at depth zob ! itz_tr = 37/39 ! index of d(Tz)/d(Tr) ! Initialize sensor specific array pointers ! if (goes_img) then ! iclr_sky = 7 ! index of clear sky amount ! elseif (avhrr_navy) then ! isst_navy = 7 ! index of navy sst (K) retrieval ! idata_type = 30 ! index of data type (151=day, 152=night) ! isst_hires = 31 ! index of interpolated hires sst (K) ! elseif (avhrr) then ! iclavr = 32 ! index CLAVR cloud flag with AVHRR data ! isst_hires = 33 ! index of interpolated hires sst (K) ! elseif (seviri) then ! iclr_sky = 7 ! index of clear sky amount ! endif |
PROCESSING OF SATELLITE DATA: Loop over data in this block
do n = 1,nobs
Extract analysis relative observation time.
If desired recompute 10meter wind factor
Set land/sea, snow, ice percentages and flags (no time interpolation)
Count data of different surface types
Set relative weight value
Interpolate model fields to observation location, call crtm and create jacobians. Output both tsim and tsim_clr for allsky
radiance operator
1 2 3 4 5 6 7 8 9 10 11 12 13 | if (lcw4crtm) then call call_crtm(obstype,dtime,data_s(:,n),nchanl,nreal,ich, & tvp,qvp,clw_guess,prsltmp,prsitmp, & trop5,tzbgr,dtsavg,sfc_speed, & tsim,emissivity,ptau5,ts,emissivity_k, & temp,wmix,jacobian,error_status,tsim_clr=tsim_clr)else call call_crtm(obstype,dtime,data_s(:,n),nchanl,nreal,ich, & tvp,qvp,clw_guess,prsltmp,prsitmp, & trop5,tzbgr,dtsavg,sfc_speed, & tsim,emissivity,ptau5,ts,emissivity_k, & temp,wmix,jacobian,error_status)endif |
If the CRTM returns an error flag, do not assimilate any channels for this ob and set the QC flag to ifail_crtm_qc.
For SST retrieval, use interpolated NCEP SST analysis
If using adaptive angle dependent bias correction, update the predictors for this part of bias correction. The AMSUA cloud liquid water algorithm uses total angle dependent bias correction for channels 1 and 2
Compute microwave cloud liquid water or graupel water path for bias correction and QC.
COMPUTE AND APPLY BIAS CORRECTION TO SIMULATED VALUES
Calculate cloud effect for QC
QC OBSERVATIONS BASED ON VARIOUS CRITERIA: Separate blocks for various instruments.
! ---------- IR -------------------
! QC HIRS/2, GOES, HIRS/3 and AIRS sounder data
!
ObsQCs: if (hirs .or. goessndr .or. airs .or. iasi .or. cris) then
frac_sea=data_s(ifrac_sea,n)
! NOTE: The qc in qc_irsnd uses the inverse squared obs error.
! The loop below loads array varinv_use accounting for whether the
! cloud detection flag is set. Array
! varinv_use is then used in the qc calculations.
! For the case when all channels of a sensor are passive, all
! channels with iuse_rad=-1 or 0 are used in cloud detection.
do i=1,nchanl
m=ich(i)
if (varinv(i) < tiny_r_kind) then
varinv_use(i) = zero
else
if ((icld_det(m)>0)) then
varinv_use(i) = varinv(i)
else
varinv_use(i) = zero
end if
end if
end do
call qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse(n),goessndr, &
cris,zsges,cenlat,frac_sea,pangs,trop5,zasat,tzbgr,tsavg5,tbc,tb_obs,tnoise, &
wavenumber,ptau5,prsltmp,tvp,temp,wmix,emissivity_k,ts, &
id_qc,aivals,errf,varinv,varinv_use,cld,cldp,kmax,zero_irjaco3_pole(n))
! --------- MSU -------------------
! QC MSU data
else if (msu) then
call qc_msu(nchanl,is,ndat,nsig,sea,land,ice,snow,luse(n), &
zsges,cenlat,tbc,ptau5,emissivity_k,ts,id_qc,aivals,errf,varinv)
! ---------- AMSU-A -------------------
! QC AMSU-A data
else if (amsua) then
if (adp_anglebc) then
tb_obsbc1=tb_obs(1)-cbias(nadir,ich(1))-predx(1,ich(1))
else
tb_obsbc1=tb_obs(1)-cbias(nadir,ich(1))
end if
cldeff_obs5=cldeff_obs(5) ! observed cloud effect for channel 5
call qc_amsua(nchanl,is,ndat,nsig,npred,sea,land,ice,snow,mixed,luse(n), &
zsges,cenlat,tb_obsbc1,cosza,clw,tbc,ptau5,emissivity_k,ts, &
pred,predchan,id_qc,aivals,errf,errf0,clwp_amsua,varinv,cldeff_obs5,factch6, &
cld_rbc_idx,sfc_speed,error0,clw_guess_retrieval,scatp)
! If cloud impacted channels not used turn off predictor
do i=1,nchanl
if ( (i <= 5 .or. i == 15) .and. (varinv(i)<1.e-9_r_kind) ) then
pred(3,i) = zero
end if
end do
! ---------- AMSU-B -------------------
! QC AMSU-B and MHS data
else if (amsub .or. hsb .or. mhs) then
call qc_mhs(nchanl,ndat,nsig,is,sea,land,ice,snow,mhs,luse(n), &
zsges,tbc,tb_obs,ptau5,emissivity_k,ts, &
id_qc,aivals,errf,varinv,clw,tpwc)
! ---------- ATMS -------------------
! QC ATMS data
else if (atms) then
if (adp_anglebc) then
tb_obsbc1=tb_obs(1)-cbias(nadir,ich(1))-predx(1,ich(1))
cldeff_obs5=cldeff_obs(6) ! observed cloud effect for ATMS channel 6
else
tb_obsbc1=tb_obs(1)-cbias(nadir,ich(1))
end if
call qc_atms(nchanl,is,ndat,nsig,npred,sea,land,ice,snow,mixed,luse(n), &
zsges,cenlat,tb_obsbc1,cosza,clw,tbc,ptau5,emissivity_k,ts, &
pred,predchan,id_qc,aivals,errf,errf0,clwp_amsua,varinv,cldeff_obs5,factch6, &
cld_rbc_idx,sfc_speed,error0,clw_guess_retrieval,scatp)
! ---------- GOES imager --------------
! GOES imager Q C
!
else if(goes_img)then
cld = data_s(iclr_sky,n)
do i = 1,nchanl
tb_obs_sdv(i) = data_s(i+29,n)
end do
call qc_goesimg(nchanl,is,ndat,nsig,ich,dplat(is),sea,land,ice,snow,luse(n), &
zsges,cld,tzbgr,tb_obs,tb_obs_sdv,tbc,tnoise,temp,wmix,emissivity_k,ts,id_qc,aivals,errf,varinv)
! ---------- SEVIRI -------------------
! SEVIRI Q C
else if (seviri) then
cld = 100-data_s(iclr_sky,n)
call qc_seviri(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse(n), &
zsges,tzbgr,tbc,tnoise,temp,wmix,emissivity_k,ts,id_qc,aivals,errf,varinv)
!
! ---------- AVRHRR --------------
! NAVY AVRHRR Q C
else if (avhrr_navy .or. avhrr) then
frac_sea=data_s(ifrac_sea,n)
! NOTE: The qc in qc_avhrr uses the inverse squared obs error.
! The loop below loads array varinv_use accounting for whether the
! cloud detection flag is set. Array
! varinv_use is then used in the qc calculations.
! For the case when all channels of a sensor are passive, all
! channels with iuse_rad=-1 or 0 are used in cloud detection.
do i=1,nchanl
m=ich(i)
if (varinv(i) < tiny_r_kind) then
varinv_use(i) = zero
else
if ((icld_det(m)>0)) then
varinv_use(i) = varinv(i)
else
varinv_use(i) = zero
end if
end if
end do
call qc_avhrr(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse(n), &
zsges,cenlat,frac_sea,pangs,trop5,tzbgr,tsavg5,tbc,tb_obs,tnoise, &
wavenumber,ptau5,prsltmp,tvp,temp,wmix,emissivity_k,ts, &
id_qc,aivals,errf,varinv,varinv_use,cld,cldp)
! ---------- SSM/I , SSMIS, AMSRE -------------------
! SSM/I, SSMIS, & AMSRE Q C
else if( ssmi .or. amsre .or. ssmis )then
frac_sea=data_s(ifrac_sea,n)
if(amsre)then
bearaz= (270._r_kind-data_s(ilazi_ang,n))*deg2rad
sun_zenith=data_s(iszen_ang,n)*deg2rad
sun_azimuth=(r90-data_s(isazi_ang,n))*deg2rad
sgagl = acos(coscon * cos( bearaz ) * cos( sun_zenith ) * cos( sun_azimuth ) + &
coscon * sin( bearaz ) * cos( sun_zenith ) * sin( sun_azimuth ) + &
sincon * sin( sun_zenith )) * rad2deg
end if
call qc_ssmi(nchanl,nsig,ich, &
zsges,luse(n),sea,mixed, &
temp,wmix,ts,emissivity_k,ierrret,kraintype,tpwc,clw,sgagl,tzbgr, &
tbc,tbcnob,tsim,tnoise,ssmi,amsre_low,amsre_mid,amsre_hig,ssmis, &
varinv,errf,aivals(1,is),id_qc)
! ---------- AMSR2 -------------------
! AMSR2 Q C
else if (amsr2) then
sun_azimuth=data_s(isazi_ang,n)
sun_zenith=data_s(iszen_ang,n)
call qc_amsr2(nchanl,zsges,luse(n),sea, &
kraintype,clw_obs,tsavg5,tb_obs,sun_azimuth,sun_zenith,amsr2,varinv,aivals(1,is),id_qc)
! ---------- GMI -------------------
! GMI Q C
else if (gmi) then
call qc_gmi(nchanl,zsges,luse(n),sea,cenlat, &
kraintype,clw_obs,tsavg5,tb_obs,gmi,varinv,aivals(1,is),id_qc)
! ---------- SAPHIR -----------------
! SAPHIR Q C
else if (saphir) then
call qc_saphir(nchanl,zsges,luse(n),sea, &
kraintype,varinv,aivals(1,is),id_qc)
! ---------- SSU -------------------
! SSU Q C
elseif (ssu) then
call qc_ssu(nchanl,is,ndat,nsig,sea,land,ice,snow,luse(n), &
zsges,cenlat,tb_obs,ptau5,emissivity_k,ts,id_qc,aivals,errf,varinv)
end if ObsQCs
! Done with sensor qc blocks. Now make final qc decisions.
! Apply gross check to observations. Toss obs failing test.
do i = 1,nchanl
if (varinv(i) > tiny_r_kind ) then
m=ich(i)
if(lcw4crtm .and. sea) then
if (i <= 3 .or. i==15) then
errf(i) = 3.00_r_kind*errf(i)
else if (i == 4) then
errf(i) = 3.00_r_kind*errf(i)
else if (i == 5) then
errf(i) = 3.00_r_kind*errf(i)
else
errf(i) = min(three*errf(i),ermax_rad(m))
endif
else if (ssmis) then
errf(i) = min(1.5_r_kind*errf(i),ermax_rad(m)) ! tighten up gross check for SSMIS
else if (gmi .or. saphir .or. amsr2) then
errf(i) = ermax_rad(m) ! use ermax for GMI, SAPHIR, and AMSR2 gross check
else
errf(i) = min(three*errf(i),ermax_rad(m))
endif
if (abs(tbc(i)) > errf(i)) then
! 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(id_qc(i) == igood_qc)id_qc(i)=ifail_gross_qc
varinv(i) = zero
if(luse(n))stats(2,m) = stats(2,m) + one
if(luse(n))aivals(7,is) = aivals(7,is) + one
end if
end if
end do
if(amsua .or. amsub .or. mhs .or. msu .or. hsb)then
if(amsua)nlev=6
if(amsub .or. mhs)nlev=5
if(hsb)nlev=4
if(msu)nlev=4
kval=0
do i=2,nlev
! do i=1,nlev
channel_passive=iuse_rad(ich(i))==-1 .or. iuse_rad(ich(i))==0
if (varinv(i)<tiny_r_kind .and. ((iuse_rad(ich(i))>=1) .or. &
(passive_bc .and. channel_passive))) then
kval=max(i-1,kval)
if(amsub .or. hsb .or. mhs)kval=nlev
if(amsua .and. i <= 3)kval = zero
end if
end do
if(kval > 0)then
do i=1,kval
varinv(i)=zero
if(id_qc(i) == igood_qc)id_qc(i)=ifail_interchan_qc
end do
if(amsua)then
varinv(15)=zero
if(id_qc(15) == igood_qc)id_qc(15)=ifail_interchan_qc
end if
end if
end if
|
Only process observations to be assimilated (fill in radhead, radtail)
At the end of analysis, prepare for bias correction for monitored channels.
Only "good monitoring" obs are included in J_passive calculation.
summation of observation number
Write diagnostics to output file: diagbuf; diagbufex; diagbufchan
! Write diagnostics to output file.
if (rad_diagsave .and. luse(n) .and. nchanl_diag > 0) then
diagbuf(1) = cenlat ! observation latitude (degrees)
diagbuf(2) = cenlon ! observation longitude (degrees)
diagbuf(3) = zsges ! model (guess) elevation at observation location
diagbuf(4) = dtime-time_offset ! observation time (hours relative to analysis time)
diagbuf(5) = data_s(iscan_pos,n) ! sensor scan position
diagbuf(6) = zasat*rad2deg ! satellite zenith angle (degrees)
diagbuf(7) = data_s(ilazi_ang,n) ! satellite azimuth angle (degrees)
diagbuf(8) = pangs ! solar zenith angle (degrees)
diagbuf(9) = data_s(isazi_ang,n) ! solar azimuth angle (degrees)
diagbuf(10) = sgagl ! sun glint angle (degrees) (sgagl)
diagbuf(11) = surface(1)%water_coverage ! fractional coverage by water
diagbuf(12) = surface(1)%land_coverage ! fractional coverage by land
diagbuf(13) = surface(1)%ice_coverage ! fractional coverage by ice
diagbuf(14) = surface(1)%snow_coverage ! fractional coverage by snow
if(.not. retrieval)then
diagbuf(15) = surface(1)%water_temperature ! surface temperature over water (K)
diagbuf(16) = surface(1)%land_temperature ! surface temperature over land (K)
diagbuf(17) = surface(1)%ice_temperature ! surface temperature over ice (K)
diagbuf(18) = surface(1)%snow_temperature ! surface temperature over snow (K)
diagbuf(19) = surface(1)%soil_temperature ! soil temperature (K)
if (gmi .or. saphir) then
diagbuf(20) = gwp ! graupel water path
else
diagbuf(20) = surface(1)%soil_moisture_content ! soil moisture
endif
diagbuf(21) = surface(1)%land_type ! surface land type
else
diagbuf(15) = tsavg5 ! SST first guess used for SST retrieval
diagbuf(16) = sstcu ! NCEP SST analysis at t
diagbuf(17) = sstph ! Physical SST retrieval
diagbuf(18) = sstnv ! Navy SST retrieval
diagbuf(19) = dta ! d(ta) corresponding to sstph
diagbuf(20) = dqa ! d(qa) corresponding to sstph
diagbuf(21) = dtp_avh ! data type
endif
if(lcw4crtm .and. sea) then
! diagbuf(22) = tpwc_amsua
diagbuf(22) = scat ! scattering index from AMSU-A
diagbuf(23) = clw_guess ! integrated CLWP (kg/m**2) from background
else
diagbuf(22) = surface(1)%vegetation_fraction ! vegetation fraction
diagbuf(23) = surface(1)%snow_depth ! snow depth
endif
diagbuf(24) = surface(1)%wind_speed ! surface wind speed (m/s)
! Note: The following quantities are not computed for all sensors
if (.not.microwave) then
diagbuf(25) = cld ! cloud fraction (%)
diagbuf(26) = cldp ! cloud top pressure (hPa)
else
if((lcw4crtm .and. sea) .or. gmi .or. amsr2) then
if (gmi .or. amsr2) then
diagbuf(25) = clw_obs ! clw (kg/m**2) from retrievals
else
diagbuf(25) = clwp_amsua ! cloud liquid water (kg/m**2)
endif
diagbuf(26) = clw_guess_retrieval ! retrieved CLWP (kg/m**2) from simulated BT
else
diagbuf(25) = clw ! cloud liquid water (kg/m**2)
diagbuf(26) = tpwc ! total column precip. water (km/m**2)
endif
endif
! For NST
if (nstinfo==0) then
diagbuf(27) = r_missing
diagbuf(28) = r_missing
diagbuf(29) = r_missing
diagbuf(30) = r_missing
else
diagbuf(27) = data_s(itref,n)
diagbuf(28) = data_s(idtw,n)
diagbuf(29) = data_s(idtc,n)
diagbuf(30) = data_s(itz_tr,n)
endif
if (lwrite_peakwt) then
do i=1,nchanl_diag
diagbufex(1,i)=weightmax(ich_diag(i)) ! press. at max of weighting fn (mb)
end do
if (goes_img) then
do i=1,nchanl_diag
diagbufex(2,i)=tb_obs_sdv(ich_diag(i))
end do
end if
else if (goes_img .and. .not.lwrite_peakwt) then
do i=1,nchanl_diag
diagbufex(1,i)=tb_obs_sdv(ich_diag(i))
end do
end if
do i=1,nchanl_diag
diagbufchan(1,i)=tb_obs(ich_diag(i)) ! observed brightness temperature (K)
diagbufchan(2,i)=tbc(ich_diag(i)) ! observed - simulated Tb with bias corrrection (K)
diagbufchan(3,i)=tbcnob(ich_diag(i)) ! observed - simulated Tb with no bias correction (K)
errinv = sqrt(varinv(ich_diag(i)))
diagbufchan(4,i)=errinv ! inverse observation error
useflag=one
if (iuse_rad(ich(ich_diag(i))) < 1) useflag=-one
diagbufchan(5,i)= id_qc(ich_diag(i))*useflag ! quality control mark or event indicator
if (lcw4crtm) then
diagbufchan(6,i)=error0(ich_diag(i))
else
diagbufchan(6,i)=emissivity(ich_diag(i)) ! surface emissivity
endif
diagbufchan(7,i)=tlapchn(ich_diag(i)) ! stability index
if (lcw4crtm) then
diagbufchan(8,i)=cld_rbc_idx(ich_diag(i)) ! indicator of cloudy consistency
else
diagbufchan(8,i)=ts(ich_diag(i)) ! d(Tb)/d(Ts)
end if
if (lwrite_predterms) then
predterms=zero
do j = 1,npred
predterms(j) = pred(j,ich_diag(i))
end do
predterms(npred+1) = cbias(nadir,ich(ich_diag(i)))
do j=1,npred+2
diagbufchan(ipchan_radiag+j,i)=predterms(j) ! Tb bias correction terms (K)
end do
else ! Default to write out predicted bias
do j=1,npred+2
diagbufchan(ipchan_radiag+j,i)=predbias(j,ich_diag(i)) ! Tb bias correction terms (K)
end do
end if
end do |
End of n-loop over obs