subroutine read_ssmis(mype,val_ssmis,ithin,isfcalc,rmesh,jsatid,gstime,&
infile,lunout,obstype,nread,ndata,nodata,twind,sis,&
mype_root,mype_sub,npe_sub,mpi_comm_sub,nobs, &
nrec_start,dval_use)
! subprogram: read_ssmis read ssmis data
! prgmmr: okamoto org: np23 date: 2005-01-05
!
! abstract: This routine reads BUFR format SSM/IS radiance
! (brightness temperature) files. Optionally, the data
! are thinned to a specified resolution using simple
! quality control checks.
! 1) obs time check |obs-anal|<time_window;
! 2) remove overlap orbit;
! 3) climate check reject for tb<tbmin or tb>tbmax
!
! When running the gsi in regional mode, the code only
! retains those observations that fall within the regional
! domain
!
! program history log:
! 2005-01-05 okamoto
! 2005-10-07 Xu & Pawlak - modify the code related to ityp determination to
! use routine deter_ityp, added values for constants
! rlndsea for four ssmis instruments, fixed indentation
! 2005-10-10 treadon - replace deter_ityp with deter_sfc, modify rlndsea to be
! consistent with other read_* routines
! 2005-10-17 treadon - add grid and earth relative obs location to output file
! 2005-10-18 treadon - remove array obs_load and call to sumload
! 2005-11-29 parrish - modify getsfc to work for different regional options
! 2005-12-15 treadon - patch to constrain ssmi_img scan positions to be in 1-90 range
! 2006-02-01 parrish - remove getsfc (different version called now in read_obs)
! 2006-02-03 derber - modify for new obs control and obs count
! 2006-04-27 derber - some efficiency modifications
! 2006-05-19 eliu - add logic to reset relative weight when all channels not used
! 2006-07-28 derber - add solar and satellite azimuth angles remove isflg from output
! 2006-08-25 treadon - replace serial bufr i/o with parallel bufr i/o (mpi_io)
! 2007-01-24 kazumori- modify to read UKMO preprocessed SSMIS data
! 2008-04-08 Yan - fix bug in calculation of ifov for UPP SSMIS data
! 2007-03-01 tremolet - measure time from beginning of assimilation window
! 2008-05-27 safford - rm unused vars and uses
! 2009-01-09 gayno - new option to calculate surface fields within FOV
! (when isfcalc flag set to one)
! 2009-04-18 woollen - improve mpi_io interface with bufrlib routines
! 2009-04-21 derber - add ithin to call to makegrids
! 2011-04-08 li - (1) use nst_gsi, nstinfo, fac_dtl, fac_tsl and add NSST vars
! (2) get zob, tz_tr (call skindepth and cal_tztr)
! (3) interpolate NSST Variables to Obs. location (call deter_nst)
! (4) add more elements (nstinfo) in data array
! 2011-08-01 lueken - added module use deter_sfc_mod
! 2011-09-02 gayno - add processing of future satellites for FOV-based
! surface field calculation and improved its error handling
! (isfcalc=1)
! 2011-12-10 eliu - add handling and call to zensun to calculate solar zenith/azimuth angles;
! fix UTC hour used by zensun
! 2012-01-10 eliu - add handling to do spatial averaging (noise reduction) for
! observed brightness temperatures
! 2012-03-05 akella - nst now controlled via coupler
! 2012-07-10 sienkiewicz add control for choosing noise reduction method 0=no smoothing
! 2013-01-26 parrish - change from grdcrd to grdcrd1 (to allow successful debug compile on WCOSS)
! 2013-01-26 parrish - WCOSS debug compile error--change mype from intent(inout) to intent(in)
! 2014-12-03 derber remove unused variables
! 2015-02-23 Rancic/Thomas - add thin4d to time window logical
!
! input argument list:
! mype - mpi task id
! val_ssmis- weighting factor applied to super obs
! ithin - flag to thin data
! isfcalc - flag to specify method to calculate sfc fields within FOV
! rmesh - thinning mesh size (km)
! jsatid - satellite to read ex.15
! gstime - analysis time in minutes from reference date
! infile - unit from which to read BUFR data
! lunout - unit to which to write data for further processing
! obstype - observation type to process
! twind - input group time window (hours)
! sis - satellite/instrument/sensor indicator
! mype_root - "root" task for sub-communicator
! mype_sub - mpi task id within sub-communicator
! npe_sub - number of data read tasks
! mpi_comm_sub - sub-communicator for data read
! nrec_start - first subset with useful information
!
! output argument list:
! nread - number of BUFR MI 1b observations read
! ndata - number of BUFR MI 1b profiles retained for further processing
! nodata - number of BUFR MI 1b observations retained for further processing
! nobs - array of observations on each subdomain for each processor
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!$$$ end documentation block
use kinds, only: r_kind,r_double,i_kind
use satthin, only: super_val,itxmax,makegrids,map2tgrid,destroygrids, &
checkob,finalcheck,score_crit
use radinfo, only: ssmis_method
use radinfo, only: iuse_rad,newchn,nusis,jpch_rad,&
use_edges,radedge1,radedge2
use gridmod, only: diagnostic_reg,regional,rlats,rlons,nlat,nlon,&
tll2xy,txy2ll
use constants, only: deg2rad,rad2deg,zero,half,one,two,four,r60inv
use gsi_4dvar, only: l4dvar,l4densvar,iwinbgn,winlen,thin4d
use calc_fov_conical, only: instrument_init
use deter_sfc_mod, only: deter_sfc,deter_sfc_fov
use gsi_nstcouplermod, only: nst_gsi,nstinfo
use gsi_nstcouplermod, only: gsi_nstcoupler_skindepth, gsi_nstcoupler_deter
use ssmis_spatial_average_mod, only : ssmis_spatial_average
use m_sortind
use mpimod, only: npe
! use radiance_mod, only: rad_obs_type
implicit none
! Declare passed variables
character(len=*),intent(in ) :: infile,obstype,jsatid
character(len=20),intent(in ) :: sis
real(r_kind) ,intent(in ) :: rmesh,gstime,twind
real(r_kind) ,intent(inout) :: val_ssmis
integer(i_kind) ,intent(in ) :: mype,nrec_start
integer(i_kind) ,intent(inout) :: lunout,ithin,isfcalc
integer(i_kind),dimension(npe) ,intent(inout) :: nobs
integer(i_kind) ,intent(inout) :: nread
integer(i_kind) ,intent(inout) :: ndata,nodata
integer(i_kind) ,intent(in ) :: mype_root
integer(i_kind) ,intent(in ) :: mype_sub
integer(i_kind) ,intent(in ) :: npe_sub
integer(i_kind) ,intent(in ) :: mpi_comm_sub
logical ,intent(in ) :: dval_use
! Declare local variables
character(7),parameter :: fov_flag="conical"
integer(i_kind),parameter :: maxchanl = 24
integer(i_kind),parameter :: mxscen_img = 180 !img
integer(i_kind),parameter :: mxscen_env = 90 !env
integer(i_kind),parameter :: mxscen_las = 60 !las
integer(i_kind),parameter :: mxscen_uas = 30 !uas
integer(i_kind),parameter :: maxobs = 800000
real(r_kind),parameter :: r360=360.0_r_kind
real(r_kind),parameter :: tbmin=70.0_r_kind
real(r_kind),parameter :: tbmax=320.0_r_kind
real(r_kind),parameter :: one_minute=0.01666667_r_kind
logical :: do_noise_reduction
logical :: ssmis_las,ssmis_uas,ssmis_img,ssmis_env,ssmis
logical :: outside,iuse,assim,valid
character(len=8):: subset
integer(i_kind) :: i,k,ifovoff,ntest
integer(i_kind) :: nlv,idate,nchanl,nreal
integer(i_kind) :: n,ireadsb,ireadmg,irec
integer(i_kind) :: nmind,itx,nele,itt
integer(i_kind) :: iskip
integer(i_kind) :: lnbufr,isflg,idomsfc(1)
integer(i_kind) :: ilat,ilon
integer(i_kind) :: nscan,jc,bufsat,incangl,said
integer(i_kind) :: nfov_bad
integer(i_kind) :: ichan, instr
integer(i_kind) :: radedge_min, radedge_max
integer(i_kind) :: iobs,num_obs,method,iret
integer(i_kind) :: irain
integer(i_kind) :: doy,mon,m
integer(i_kind) :: ibfms,maxinfo
! integer(i_kind),pointer :: ifov,iscan,iorbn,inode
integer(i_kind),pointer :: ifov,inode
integer(i_kind),allocatable :: sorted_index(:)
integer(i_kind),allocatable,target :: ifov_save(:)
! integer(i_kind),allocatable,target :: iscan_save(:)
! integer(i_kind),allocatable,target :: iorbn_save(:)
integer(i_kind),allocatable,target :: inode_save(:)
integer(i_kind),dimension(12):: mlen,mday
integer(i_kind),dimension(5) :: iobsdate
integer(i_kind),allocatable :: nrec(:)
real(r_kind) :: sfcr,r07
! real(r_kind) :: pred
real(r_kind) :: tdiff,timedif,dist1
! real(r_kind) :: step,start
real(r_kind) :: tsavg,vty,vfr,sty,stp,sm,sn,zz,ff10
real(r_kind) :: zob,tref,dtw,dtc,tz_tr
real(r_kind) :: disterr,disterrmax,cdist,dlon00,dlat00
real(r_kind) :: fovn,sscan,orbn,rainf
! real(r_kind) :: sort_time1, sort_time2
real(r_kind) :: flgch
real(r_kind) :: clat,clon
real(r_kind) :: dlat,dlon
real(r_kind) :: dlon_earth_deg,dlat_earth_deg,expansion,sat_aziang
real(r_kind) :: utc_hour,sun_zenith,sun_azimuth
real(r_double),dimension(7) :: bufrinit
real(r_double),dimension(3,5) :: bufrymd
real(r_double),dimension(2,2) :: bufrhm
real(r_double),dimension(2,29) :: bufrloc
real(r_double),dimension(1,maxchanl):: bufrtbb
real(r_double) :: rnode
real(r_kind),dimension(0:3) :: sfcpct
real(r_kind),dimension(0:4) :: rlndsea
real(r_kind),dimension(0:3) :: ts
real(r_kind),pointer :: bt_in(:)
real(r_kind),pointer :: crit1,rsat,t4dv,solzen,solazi,dlon_earth,dlat_earth,satazi,lza
real(r_kind),allocatable,target :: rsat_save(:)
real(r_kind),allocatable,target :: t4dv_save(:)
real(r_kind),allocatable,target :: dlon_earth_save(:)
real(r_kind),allocatable,target :: dlat_earth_save(:)
real(r_kind),allocatable,target :: crit1_save(:)
real(r_kind),allocatable,target :: lza_save(:)
real(r_kind),allocatable,target :: satazi_save(:)
real(r_kind),allocatable,target :: solzen_save(:)
real(r_kind),allocatable,target :: solazi_save(:)
real(r_kind),allocatable,target :: bt_save(:,:)
real(r_kind),allocatable :: relative_time_in_seconds(:)
real(r_kind),allocatable :: data_all(:,:)
! For solar zenith/azimuth angles calculation
data mlen/31,28,31,30,31,30, &
31,31,30,31,30,31/
!----------------------------------------------------------------------
! Initialize variables
maxinfo = 31
m = 0
do mon=1,12
mday(mon) = m
m = m + mlen(mon)
end do
do_noise_reduction = .true.
if (ssmis_method .eq. 0) do_noise_reduction = .false.
nchanl = maxchanl
disterrmax = zero
ntest = 0
ndata = 0
nodata = 0
nread = 0
nfov_bad = 0
ilon = 3
ilat = 4
lnbufr = 15
r07 = 0.7_r_kind * deg2rad
if (nst_gsi > 0 ) then
call gsi_nstcoupler_skindepth(obstype, zob) ! get penetration depth (zob) for the obstype
endif
! If all channels of a given sensor are set to monitor or not
! assimilate mode (iuse_rad<1), reset relative weight to zero.
! We do not want such observations affecting the relative
! weighting between observations within a given thinning group.
assim=.false.
search: do i=1,jpch_rad
if ((nusis(i)==sis) .and. (iuse_rad(i)>0)) then
assim=.true.
exit search
endif
end do search
if (.not.assim) val_ssmis=zero
! Make thinning grids
call makegrids(rmesh,ithin)
! Set various variables depending on type of data to be read
ssmis_uas= obstype == 'ssmis_uas'
ssmis_las= obstype == 'ssmis_las'
ssmis_img= obstype == 'ssmis_img'
ssmis_env= obstype == 'ssmis_env'
ssmis = obstype == 'ssmis'
! Define WMO satellite number
bufsat = 0
if (trim(sis) == 'ssmis_f16') bufsat=249
if (trim(sis) == 'ssmis_f17') bufsat=285
if (trim(sis) == 'ssmis_f18') bufsat=286
if (trim(sis) == 'ssmis_f19') bufsat=287
write(6,*) 'READ_SSMIS: reading bufsat = ', bufsat, trim(sis)
! Humidity imager:180
if(ssmis)then
nscan = mxscen_las
ifovoff = 270
incangl = 53.0_r_kind
else if(ssmis_img) then
nscan = mxscen_img
ifovoff = 0
incangl = 53.0_r_kind
! env:90
else if(ssmis_env) then
nscan = mxscen_env
ifovoff = 180
incangl = 53.1_r_kind
! las:60
else if(ssmis_las) then
nscan = mxscen_las
ifovoff = 270
incangl = 53.0_r_kind
! uas:30
else if(ssmis_uas) then
nscan = mxscen_uas
ifovoff = 330
incangl = 52.4_r_kind
end if
! Initialize variables for use by FOV-based surface code
if (isfcalc == 1) then
instr=25 ! circular fov, use as default
if (trim(jsatid) == 'f16') instr=26
if (trim(jsatid) == 'f17') instr=27
if (trim(jsatid) == 'f18') instr=28
if (trim(jsatid) == 'f19') instr=29
! right now, all ssmis data is mapped to a common fov -
! that of the las channels.
ichan = 1
expansion = 2.9_r_kind
sat_aziang = 90.0_r_kind ! 'fill' value; need to get this from file
call instrument_init(instr, jsatid, expansion, valid)
if (.not. valid) then
if (assim) then
write(6,*)'READ_SSMIS: ***ERROR*** IN SETUP OF FOV-SFC CODE. STOP'
call stop2(71)
else
isfcalc = 0
write(6,*)'READ_SSMIS: ***ERROR*** IN SETUP OF FOV-SFC CODE'
endif
endif
endif
radedge_min = 0
radedge_max = 1000
do i=1,jpch_rad
if (trim(nusis(i))==trim(sis)) then
! step = radstep(i)
! start = radstart(i)
if (radedge1(i)/=-1 .or. radedge2(i)/=-1) then
radedge_min=radedge1(i)
radedge_max=radedge2(i)
end if
exit
endif
end do
rlndsea(0) = zero
rlndsea(1) = 15._r_kind
rlndsea(2) = 10._r_kind
rlndsea(3) = 15._r_kind
rlndsea(4) = 100._r_kind
! Allocate arrays for BUFR I/O
allocate(ifov_save(maxobs))
! allocate(iscan_save(maxobs))
! allocate(iorbn_save(maxobs))
allocate(inode_save(maxobs))
allocate(rsat_save(maxobs))
allocate(t4dv_save(maxobs))
allocate(dlon_earth_save(maxobs))
allocate(dlat_earth_save(maxobs))
allocate(crit1_save(maxobs))
allocate(lza_save(maxobs))
allocate(satazi_save(maxobs))
allocate(solzen_save(maxobs))
allocate(solazi_save(maxobs))
allocate(bt_save(maxchanl,maxobs))
inode_save = 0
! Read in data from bufr into arrays first
! Open unit to satellite bufr file
iobs=1
call closbf(lnbufr)
open(lnbufr,file=trim(infile),form='unformatted',status='old',err=500)
call openbf(lnbufr,'IN',lnbufr)
call datelen(10)
! Loop to read bufr file
irec=0
read_subset: do while(ireadmg(lnbufr,subset,idate)>=0 .and. iobs < maxobs)
irec = irec + 1
if(irec < nrec_start) cycle read_subset
read_loop: do while(ireadsb(lnbufr)==0 .and. iobs < maxobs)
rsat => rsat_save(iobs)
t4dv => t4dv_save(iobs)
dlon_earth => dlon_earth_save(iobs)
dlat_earth => dlat_earth_save(iobs)
crit1 => crit1_save(iobs)
ifov => ifov_save(iobs)
! iscan => iscan_save(iobs)
! iorbn => iorbn_save(iobs)
inode => inode_save(iobs)
lza => lza_save(iobs)
satazi => satazi_save(iobs)
solzen => solzen_save(iobs)
solazi => solazi_save(iobs)
! BUFR read 1/3
call ufbint(lnbufr,bufrinit,7,1,nlv, &
"SAID SECO SLNM FOVN SFLG RFLAG ORBN" )
! Extract satellite id. If not the one we want, read next record
said = nint( bufrinit(1))
if( said /= bufsat) cycle read_subset
rainf = bufrinit(6)
irain = nint(rainf)
! Rain check (-1=indeterminate 0=no rain 1=rain)
if(irain == 1 .or. irain < 0) cycle read_loop ! rain check
rsat=bufsat
fovn = bufrinit(4)
sscan = bufrinit(3)
orbn = bufrinit(7)
ifov = nint(fovn)
! iscan = nint(sscan)
! iorbn = nint(orbn)
! if not doing noise reduction, try reading node information
if ( .not. do_noise_reduction ) then
inode = 1000
call ufbint(lnbufr,rnode,1,1,nlv, 'STKO')
if (ibfms(rnode) == 0) then
if (rnode == 1.) inode = -1
if (rnode == 0.) inode = 1
end if
end if
! Check date/time
! BUFR read 2/3 --- read in observation date/time
call ufbrep(lnbufr,bufrymd,3,5,nlv,"YEAR MNTH DAYS" )
call ufbrep(lnbufr,bufrhm, 2,2,nlv,"HOUR MINU" )
! Calc obs seqential time If time outside window, skip this obs
iobsdate(1:3) = bufrymd(1:3,1) !year,month,day for scan start time kozo
iobsdate(4:5) = bufrhm(1:2,1) !hour,min for scan start time kozo
call w3fs21(iobsdate,nmind)
t4dv=(real(nmind-iwinbgn,r_kind) + real(bufrinit(2),r_kind)*r60inv)*r60inv
tdiff=t4dv+(iwinbgn-gstime)*r60inv
if (l4dvar.or.l4densvar) then
if (t4dv<zero .OR. t4dv>winlen) cycle read_loop
else
if(abs(tdiff) > twind+one_minute) cycle read_loop
endif
if (thin4d) then
! Give score based on time in the window
! crit1 = 0.01_r_kind+ flgch
crit1 = zero
else
timedif = 6.0_r_kind*abs(tdiff) ! range: 0 to 18
! crit1 = 0.01_r_kind+timedif + flgch
crit1 = timedif
endif
! Extract obs location, TBB, other information
! BUFR read 3/3 --- read in observation lat/lon
call ufbrep(lnbufr,bufrloc, 2,29, nlv,"CLAT CLON" )
! call ufbrep(lnbufr,bufrinfo, 1,3, nlv,"SELV" )
! call ufbrep(lnbufr,bufrlleaa,2,28, nlv,"RAIA BEARAZ" )
dlat_earth = bufrloc(1,1) !degrees
dlon_earth = bufrloc(2,1) !degrees
if(abs(dlat_earth)>90.0_r_kind .or. abs(dlon_earth)>r360) cycle read_loop
if(dlon_earth < zero) dlon_earth = dlon_earth+r360
if(dlon_earth == r360) dlon_earth = dlon_earth-r360
lza = incangl ! conical scanning instrument has fixed local zenith angle
satazi = 0 ! currently missing from bufr
! Calculate solar zenith/azimuth angle
clat = dlat_earth
clon = dlon_earth
if(clon > 180_r_kind) clon = clon-360.0_r_kind
utc_hour = real(iobsdate(4),r_kind)+real(iobsdate(5),r_kind)*r60inv+real(bufrinit(2),r_kind)*r60inv*r60inv
! utc_hour = real(iobsdate(4),r_kind) !orig
doy = mday( int(iobsdate(2)) ) + int(iobsdate(3))
if ((mod( int(iobsdate(1)),4)==0).and.( int(iobsdate(2)) > 2)) then
doy = doy + 1
end if
call zensun(doy,utc_hour,clat,clon,sun_zenith,sun_azimuth) ! output solar zenith angles are between -90 and 90
sun_zenith = 90.-sun_zenith ! make sure solar zenith angles are between 0 and 180
solzen = sun_zenith
solazi = sun_azimuth
! Check Tb
! Transfer observed brightness temperature to work array.
call ufbrep(lnbufr,bufrtbb, 1,maxchanl,nlv,"TMBR" )
bt_save(1:maxchanl,iobs) = bufrtbb(1,1:maxchanl)
iobs=iobs+1
end do read_loop
end do read_subset
call closbf(lnbufr)
num_obs = iobs-1
500 continue
write(*,*) 'READ_SSMIS: num_obs = ', num_obs, num_obs*nchanl
if (num_obs <= 0 ) then
write(*,*) 'READ_SSMIS: No ', trim(sis), &
' data read in at mype mype_sub ', mype, mype_sub
return
endif
if (do_noise_reduction) then
! call cpu_time(sort_time1)
write(*,*) 'READ_SSMIS: num_obs = ', num_obs, num_obs*nchanl
! Sort time in ascending order and get sorted index
! relative_time_in_seconds referenced at the beginning of the assimilation window
allocate(relative_time_in_seconds(num_obs))
allocate(sorted_index(num_obs))
relative_time_in_seconds = 3600.0_r_kind*t4dv_save(1:num_obs)
sorted_index = sortind(relative_time_in_seconds)
! Sort data according to observation time in ascending order
relative_time_in_seconds(1:num_obs) = relative_time_in_seconds(sorted_index)
rsat_save(1:num_obs) = rsat_save(sorted_index)
t4dv_save(1:num_obs) = t4dv_save(sorted_index)
dlon_earth_save(1:num_obs) = dlon_earth_save(sorted_index)
dlat_earth_save(1:num_obs) = dlat_earth_save(sorted_index)
crit1_save(1:num_obs) = crit1_save(sorted_index)
ifov_save(1:num_obs) = ifov_save(sorted_index)
! iscan_save(1:num_obs) = iscan_save(sorted_index)
! iorbn_save(1:num_obs) = iorbn_save(sorted_index)
lza_save(1:num_obs) = lza_save(sorted_index)
satazi_save(1:num_obs) = satazi_save(sorted_index)
solzen_save(1:num_obs) = solzen_save(sorted_index)
solazi_save(1:num_obs) = solazi_save(sorted_index)
bt_save(:,1:num_obs) = bt_save(:,sorted_index)
! call cpu_time(sort_time2)
! write(*,*)'READ_SSMIS: cpu_time (sorting) ', sort_time2-sort_time1
! write(*,*)'READ_SSMIS: min/max time ', minval(relative_time_in_seconds(1:num_obs)), &
! maxval(relative_time_in_seconds(1:num_obs))
! write(*,*)'READ_SSMIS: min/max lat ', minval(dlat_earth_save(1:num_obs)), &
! maxval(dlat_earth_save(1:num_obs))
! write(*,*)'READ_SSMIS: min/max lon ', minval(dlon_earth_save(1:num_obs)), &
! maxval(dlon_earth_save(1:num_obs))
! write(*,*)'READ_SSMIS: min/max iscan_save ', minval(iscan_save(1:num_obs)), &
! maxval(iscan_save(1:num_obs))
! write(*,*)'READ_SSMIS: min/max ifov_save ', minval(ifov_save(1:num_obs)), &
! maxval(ifov_save(1:num_obs))
! write(*,*)'READ_SSMIS: min/max bt_save ', minval(bt_save(:,1:num_obs)), &
! maxval(bt_save(:,1:num_obs))
!========================================================================================================================
! Do SSMIS spatial averaging
! method=1 --- simply averaging over circular domains centered on each field of view
! method=2 --- similar to 1 (for testing only)
! method=3 --- AAPP package
method = ssmis_method
write(*,*) 'READ_SSMIS: Calling ssmis_spatial_average, method =', method
call ssmis_spatial_average(bufsat,method,num_obs,nchanl, &
ifov_save,inode_save,relative_time_in_seconds, &
dlat_earth_save,dlon_earth_save, &
bt_save(1:nchanl,1:num_obs),iret) ! inout
if (iret /= 0) then
write(*,*) 'Error calling ssmis_spatial_average from READ_SSMIS'
return
endif
if (num_obs > 0) then
deallocate(sorted_index)
deallocate(relative_time_in_seconds)
endif
endif ! do_noise_reduction
!========================================================================================================================
! Complete thinning and QC steps for SSMIS
! Write header record to scratch file. Also allocate array
! to hold all data for given satellite
if(dval_use) maxinfo = maxinfo+2
nreal = maxinfo + nstinfo
nele = nreal + nchanl
allocate(data_all(nele,itxmax),nrec(itxmax))
nrec=999999
obsloop: do iobs = 1, num_obs
rsat => rsat_save(iobs)
t4dv => t4dv_save(iobs)
dlon_earth => dlon_earth_save(iobs)
dlat_earth => dlat_earth_save(iobs)
crit1 => crit1_save(iobs)
ifov => ifov_save(iobs)
inode => inode_save(iobs)
lza => lza_save(iobs)
satazi => satazi_save(iobs)
solzen => solzen_save(iobs)
solazi => solazi_save(iobs)
bt_in => bt_save(1:nchanl,iobs)
if (inode == 0) cycle obsloop ! this indicate duplicated data
if (.not. use_edges .and. (ifov < radedge_min .or. ifov > radedge_max)) &
cycle obsloop
dlat_earth_deg = dlat_earth
dlon_earth_deg = dlon_earth
dlat_earth = dlat_earth*deg2rad
dlon_earth = dlon_earth*deg2rad
! Regional case
if(regional)then
call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside)
if(diagnostic_reg) then
call txy2ll(dlon,dlat,dlon00,dlat00)
ntest=ntest+1
cdist=sin(dlat_earth)*sin(dlat00)+cos(dlat_earth)*cos(dlat00)* &
(sin(dlon_earth)*sin(dlon00)+cos(dlon_earth)*cos(dlon00))
cdist=max(-one,min(cdist,one))
disterr=acos(cdist)*rad2deg
disterrmax=max(disterrmax,disterr)
end if
! Check to see if in domain
if(outside) cycle obsloop
! Global case
else
dlat=dlat_earth
dlon=dlon_earth
call grdcrd1(dlat,rlats,nlat,1)
call grdcrd1(dlon,rlons,nlon,1)
endif
! Check time window
if (l4dvar.or.l4densvar) then
if (t4dv<zero .OR. t4dv>winlen) cycle obsloop
else
tdiff=t4dv+(iwinbgn-gstime)*r60inv
if(abs(tdiff) > twind) cycle ObsLoop
endif
! Map obs to thinning grid
call map2tgrid(dlat_earth,dlon_earth,dist1,crit1,itx,ithin,itt,iuse,sis)
if(.not. iuse)cycle obsloop
nread=nread+nchanl
! Transfer observed brightness temperature to work array.
! If any temperature exceeds limits, reset observation
! to "bad" value
iskip=0
do jc=1,nchanl
if(bt_in(jc)<tbmin .or. bt_in(jc)>tbmax) then
iskip = iskip + 1
endif
end do
if(iskip>=nchanl) cycle obsloop !if all ch for any position are bad, skip
flgch = iskip*two !used for thinning priority range 0-14
crit1 = crit1 + flgch
call checkob(dist1,crit1,itx,iuse)
if (.not. iuse) cycle obsloop
! Locate the observation on the analysis grid. Get sst and land/sea/ice
! mask.
! isflg - surface flag
! 0 sea
! 1 land
! 2 sea ice
! 3 snow
! 4 mixed
! sfcpct(0:3)- percentage of 4 surface types
! (0) - sea percentage
! (1) - land percentage
! (2) - sea ice percentage
! (3) - snow percentage
! FOV-based surface code requires fov number; if out-of-range, then
! skip this obs.
if (isfcalc==1) then
call deter_sfc_fov(fov_flag,ifov,instr,ichan,sat_aziang,dlat_earth_deg,&
dlon_earth_deg,expansion,t4dv,isflg,idomsfc(1), &
sfcpct,vfr,sty,vty,stp,sm,ff10,sfcr,zz,sn,ts,tsavg)
else
call deter_sfc(dlat,dlon,dlat_earth,dlon_earth,t4dv,isflg,idomsfc(1),sfcpct, &
ts,tsavg,vty,vfr,sty,stp,sm,sn,zz,ff10,sfcr)
endif ! isfcalc==1
crit1 = crit1 + rlndsea(isflg)
! call checkob(dist1,crit1,itx,iuse)
! if(.not. iuse)cycle obsloop
! Set common predictor parameters
! pred = zero
! Compute "score" for observation. All scores>=0.0. Lowest score is "best"
! crit1 = crit1+pred
call finalcheck(dist1,crit1,itx,iuse)
if(.not. iuse)cycle obsloop
!
! interpolate NSST variables to Obs. location and get dtw, dtc, tz_tr
!
if ( nst_gsi > 0 ) then
tref = ts(0)
dtw = zero
dtc = zero
tz_tr = one
if ( sfcpct(0) > zero ) then
call gsi_nstcoupler_deter(dlat_earth,dlon_earth,t4dv,zob,tref,dtw,dtc,tz_tr)
endif
endif
! Load selected observation into data array
data_all( 1,itx)= rsat ! satellite id
data_all( 2,itx)= t4dv ! time diff between obs and anal (min)
data_all( 3,itx)= dlon ! grid relative longitude
data_all( 4,itx)= dlat ! grid relative latitude
data_all( 5,itx)= incangl*deg2rad ! local zenith angle (rad)
data_all( 6,itx)= inode ! local azimuth angle (missing) ** AS/DS node infomation for SSMIS
data_all( 7,itx)= zero ! look angle (rad)
data_all( 8,itx)= ifov ! FOV scan position
! data_all( 9,itx)= zero ! solar zenith angle (deg) : not used for MW-RT calc
! data_all(10,itx)= zero ! solar azimuth angle (deg) : not used for MW-RT calc
data_all( 9,itx)= solzen ! calculated solar zenith angle (deg)
data_all(10,itx)= solazi ! calculated solar azimuth angle (deg)
data_all(11,itx) = sfcpct(0) ! sea percentage of
data_all(12,itx) = sfcpct(1) ! land percentage
data_all(13,itx) = sfcpct(2) ! sea ice percentage
data_all(14,itx) = sfcpct(3) ! snow percentage
data_all(15,itx)= ts(0) ! ocean skin temperature
data_all(16,itx)= ts(1) ! land skin temperature
data_all(17,itx)= ts(2) ! ice skin temperature
data_all(18,itx)= ts(3) ! snow skin temperature
data_all(19,itx)= tsavg ! average skin temperature
data_all(20,itx)= vty ! vegetation type
data_all(21,itx)= vfr ! vegetation fraction
data_all(22,itx)= sty ! soil type
data_all(23,itx)= stp ! soil temperature
data_all(24,itx)= sm ! soil moisture
data_all(25,itx)= sn ! snow depth
data_all(26,itx)= zz ! surface height
data_all(27,itx)= idomsfc(1) + 0.001_r_kind ! dominate surface type
data_all(28,itx)= sfcr ! surface roughness
data_all(29,itx)= ff10 ! ten meter wind factor
data_all(30,itx)= dlon_earth_deg ! earth relative longitude (degrees)
data_all(31,itx)= dlat_earth_deg ! earth relative latitude (degrees)
if(dval_use)then
data_all(32,itx)=val_ssmis
data_all(33,itx)=itt
end if
if ( nst_gsi > 0 ) then
data_all(maxinfo+1,itx) = tref ! foundation temperature
data_all(maxinfo+2,itx) = dtw ! dt_warm at zob
data_all(maxinfo+3,itx) = dtc ! dt_cool at zob
data_all(maxinfo+4,itx) = tz_tr ! d(Tz)/d(Tr)
endif
do jc=1,nchanl
data_all(nreal+jc,itx) = bt_in(jc)
end do
nrec(itx)=iobs
end do obsloop
! Deallocate I/O arrays
deallocate(rsat_save)
deallocate(ifov_save)
! deallocate(iscan_save)
deallocate(inode_save)
! deallocate(iorbn_save)
deallocate(t4dv_save)
deallocate(dlon_earth_save)
deallocate(dlat_earth_save)
deallocate(crit1_save)
deallocate(lza_save)
deallocate(satazi_save)
deallocate(solzen_save)
deallocate(solazi_save)
deallocate(bt_save)
! If multiple tasks read input bufr file, allow each tasks to write out
! information it retained and then let single task merge files together
call combine_radobs(mype_sub,mype_root,npe_sub,mpi_comm_sub,&
nele,itxmax,nread,ndata,data_all,score_crit,nrec)
! Allow single task to check for bad obs, update superobs sum,
! and write out data to scratch file for further processing.
if (mype_sub==mype_root.and.ndata>0) then
! Identify "bad" observation (unreasonable brightness temperatures).
! Update superobs sum according to observation location
do n=1,ndata
do i=1,nchanl
if(data_all(i+nreal,n) > tbmin .and. &
data_all(i+nreal,n) < tbmax)nodata=nodata+1
end do
end do
if(dval_use .and. assim)then
do n=1,ndata
itt=nint(data_all(33,n))
super_val(itt)=super_val(itt)+val_ssmis
end do
end if
! Write final set of "best" observations to output file
call count_obs(ndata,nele,ilat,ilon,data_all,nobs)
write(lunout) obstype,sis,nreal,nchanl,ilat,ilon
write(lunout) ((data_all(k,n),k=1,nele),n=1,ndata)
endif
! Deallocate local arrays
deallocate(data_all,nrec)
! Deallocate satthin arrays
1000 continue
call destroygrids
if(diagnostic_reg .and. ntest>0 .and. mype_sub==mype_root) &
write(6,*)'READ_SSMIS: mype,ntest,disterrmax=',&
mype,ntest,disterrmax
if (nfov_bad>0) &
write(6,*)'READ_SSMIS(',obstype,'): found ',nfov_bad,' questionable fov'
! End of routine
return
end subroutine read_ssmis