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