subroutine read_saphir(mype,val_tovs,ithin,isfcalc,& rmesh,jsatid,gstime,infile,lunout,obstype,& nread,ndata,nodata,twind,sis, & mype_root,mype_sub,npe_sub,mpi_comm_sub,nobs,dval_use) ! subprogram: read_saphir read bufr format saphir data ! prgmmr : ejones org: jcsda date: 2015-01-02 ! code copied from read_atms.f90 ! ! abstract: This routine reads BUFR format SAPHIR radiance (brightness ! temperature) data. ! ! When running the gsi in regional mode, the code only ! retains those observations that fall within the regional ! domain ! ! program history log: ! 2015-01-02 ejones - adapted from read_atms.f90 ! 2015-09-17 Thomas - add l4densvar and thin4d to data selection procedure ! 2016-03-09 ejones - update mnemonics for operational SAPHIR bufr ! 2016-04-01 ejones - add binning of fovs for scan angle bias correction ! 2016-07-25 ejones - remove binning of fovs ! 2016-10-05 acollard -Fix interaction with NSST and missing zenith angle issue. ! ! input argument list: ! mype - mpi task id ! val_tovs - weighting factor applied to super obs ! ithin - flag to thin data ! isfcalc - method to calculate surface fields within FOV ! when one, calculate accounting for size/shape of FOV. ! otherwise, use bilinear interpolation. ! rmesh - thinning mesh size (km) ! jsatid - satellite to read ! 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 - sensor/instrument/satellite 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 ! ! output argument list: ! nread - number of BUFR ATMS 1b observations read ! ndata - number of BUFR ATMS 1b profiles retained for further processing ! nodata - number of BUFR ATMS 1b observations retained for further processing ! nobs - array of observations on each subdomain for each processor ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: r_kind,r_double,i_kind use satthin, only: super_val,itxmax,makegrids,destroygrids,checkob, & finalcheck,map2tgrid,score_crit use radinfo, only: iuse_rad,nusis,jpch_rad, & use_edges,radedge1,radedge2,radstart,radstep use gridmod, only: diagnostic_reg,regional,nlat,nlon,tll2xy,txy2ll,rlats,rlons use constants, only: deg2rad,zero,one,two,three,rad2deg,r60inv use crtm_module, only : max_sensor_zenith_angle use calc_fov_crosstrk, only : instrument_init, fov_cleanup, fov_check use gsi_4dvar, only: l4dvar,iwinbgn,winlen,l4densvar,thin4d use deter_sfc_mod, only: deter_sfc_fov,deter_sfc use gsi_nstcouplermod, only: nst_gsi,nstinfo use gsi_nstcouplermod, only: gsi_nstcoupler_skindepth,gsi_nstcoupler_deter 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 integer(i_kind) ,intent(in ) :: mype,lunout,ithin integer(i_kind) ,intent(inout) :: isfcalc integer(i_kind) ,intent(inout) :: nread integer(i_kind) ,intent( out) :: ndata,nodata real(r_kind) ,intent(in ) :: rmesh,gstime,twind real(r_kind) ,intent(inout) :: val_tovs 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 integer(i_kind),dimension(npe) ,intent(inout) :: nobs logical ,intent(in ) :: dval_use ! Declare local parameters character(8),parameter:: fov_flag="crosstrk" integer(i_kind),parameter:: n1bhdr=12 integer(i_kind),parameter:: n2bhdr=4 integer(i_kind),parameter:: maxobs = 5000000 integer(i_kind),parameter:: max_chanl = 22 real(r_kind),parameter:: r360=360.0_r_kind real(r_kind),parameter:: tbmin=50.0_r_kind real(r_kind),parameter:: tbmax=550.0_r_kind ! The next two are one minute in hours real(r_kind),parameter:: one_minute=0.01666667_r_kind real(r_kind),parameter:: minus_one_minute=-0.01666667_r_kind real(r_kind),parameter:: rato=1.1363987_r_kind ! ratio of satellite height to ! distance from Earth's centre ! Declare local variables logical :: outside,iuse,assim,valid character(8) :: subset character(80) :: hdr1b,hdr2b integer(i_kind) :: ireadsb,ireadmg,irec integer(i_kind) :: i,j,k,ntest,iob integer(i_kind) :: iret,idate,nchanl,n,idomsfc(1) integer(i_kind) :: kidsat,maxinfo integer(i_kind) :: nmind,itx,nreal,nele,itt,num_obs integer(i_kind) :: iskip integer(i_kind) :: lnbufr,ksatid,isflg integer(i_kind) :: ilat,ilon,nadir integer(i_kind),dimension(5):: idate5 integer(i_kind) :: instr,ichan integer(i_kind) :: radedge_min, radedge_max integer(i_kind), POINTER :: ifov real(r_kind) :: sfcr real(r_kind) :: expansion real(r_kind),dimension(0:3):: sfcpct real(r_kind),dimension(0:3):: ts real(r_kind) :: tsavg,vty,vfr,sty,stp,sm,sn,zz,ff10 real(r_kind) :: zob,tref,dtw,dtc,tz_tr real(r_kind) :: dlat,dlon,tdiff,panglr real(r_kind) :: dlon_earth_deg,dlat_earth_deg real(r_kind) :: step,start,dist1 real(r_kind),dimension(0:4) :: rlndsea real(r_kind),allocatable,dimension(:,:):: data_all real(r_kind), POINTER :: bt_in(:), crit1,rsat, t4dv, solzen, solazi real(r_kind), POINTER :: dlon_earth,dlat_earth,satazi, lza integer(i_kind), ALLOCATABLE, TARGET :: ifov_save(:) 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(:,:) integer(i_kind),allocatable,dimension(:):: nrec real(r_double),allocatable,dimension(:) :: data1b8 real(r_double),dimension(n1bhdr):: bfr1bhdr real(r_double),dimension(n2bhdr):: bfr2bhdr real(r_kind) :: disterr,disterrmax,dlon00,dlat00 !************************************************************************** ! Initialize variables maxinfo=31 lnbufr = 15 disterrmax=zero ntest=0 ndata = 0 nodata = 0 nread = 0 ilon=3 ilat=4 if(nst_gsi>0) then call gsi_nstcoupler_skindepth(obstype,zob) endif ! Make thinning grids call makegrids(rmesh,ithin) ! Set nadir position nadir=65 ! Set various variables depending on type of data to be read if (obstype /= 'saphir') then write(*,*) 'READ_SAPHIR called for obstype '//obstype//': RETURNING' return end if if(jsatid == 'meghat') then kidsat=440 else write(*,*) 'READ_SAPHIR: Unrecognized value for jsatid '//jsatid//': RETURNING' return end if ! look up info in the scaninfo file 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 ! IFSCALC setup nchanl=6 if (dval_use) maxinfo = maxinfo+2 if (isfcalc==1) then instr=19 ! This section isn't really updated. expansion=2.9_r_kind ! use almost three for microwave sensors. endif ! Set rlndsea for types we would prefer selecting rlndsea(0) = zero rlndsea(1) = 15._r_kind rlndsea(2) = 10._r_kind rlndsea(3) = 15._r_kind rlndsea(4) = 100._r_kind ! 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_tovs=zero ! This may need updating, but ! Initialize variables for use by FOV-based surface code. if (isfcalc == 1) then call instrument_init(instr,jsatid,expansion,valid) if (.not. valid) then if (assim) then write(6,*)'READ_SAPHIR: ***ERROR*** IN SETUP OF FOV-SFC CODE. STOP' call stop2(71) else call fov_cleanup isfcalc = 0 write(6,*)'READ_SAPHIR: ***ERROR*** IN SETUP OF FOV-SFC CODE' endif endif endif ! Allocate arrays for BUFR I/O ALLOCATE(data1b8(nchanl)) ALLOCATE(rsat_save(maxobs)) ALLOCATE(t4dv_save(maxobs)) ALLOCATE(ifov_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(max_chanl,maxobs)) ! Reopen unit to satellite bufr file iob=1 call closbf(lnbufr) open(lnbufr,file=trim(infile),form='unformatted',status = 'old',err = 500) call openbf(lnbufr,'IN',lnbufr) hdr1b ='SAID FOVN YEAR MNTH DAYS HOUR MINU SECO CLATH CLONH' hdr2b ='IANG SOZA BEARAZ SOLAZI' ! mnemonics in non-operational SAPHIR bufr are a little different: ! hdr1b ='SAID FOVN YEAR MONTH DAY HOUR MINU SECO CLATH CLONH' ! hdr2b ='AGIND SOZA BEARAZ SOLAZI' ! AGIND instead of SAZA ! Loop to read bufr file irec=0 read_subset: do while(ireadmg(lnbufr,subset,idate)>=0 .AND. iob < maxobs) read_loop: do while (ireadsb(lnbufr)==0 .and. iob < maxobs) rsat => rsat_save(iob) t4dv => t4dv_save(iob) dlon_earth => dlon_earth_save(iob) dlat_earth => dlat_earth_save(iob) crit1 => crit1_save(iob) ifov => ifov_save(iob) lza => lza_save(iob) satazi => satazi_save(iob) solzen => solzen_save(iob) solazi => solazi_save(iob) call ufbint(lnbufr,bfr1bhdr,n1bhdr,1,iret,hdr1b) ! Extract satellite id. If not the one we want, read next record rsat=bfr1bhdr(1) ksatid=nint(bfr1bhdr(1)) if(ksatid /= kidsat) cycle read_subset ! Extract observation location and other required information if(abs(bfr1bhdr(9)) <= 90._r_kind .and. abs(bfr1bhdr(10)) <= r360)then dlat_earth = bfr1bhdr(9) dlon_earth = bfr1bhdr(10) else cycle read_loop end if if(dlon_earth=r360) dlon_earth = dlon_earth-r360 ! Extract date information. If time outside window, skip this obs idate5(1) = bfr1bhdr(3) !year idate5(2) = bfr1bhdr(4) !month idate5(3) = bfr1bhdr(5) !day idate5(4) = bfr1bhdr(6) !hour idate5(5) = bfr1bhdr(7) !minute call w3fs21(idate5,nmind) t4dv= (real((nmind-iwinbgn),r_kind) + bfr1bhdr(8)*r60inv)*r60inv ! add in seconds tdiff=t4dv+(iwinbgn-gstime)*r60inv if (l4dvar.or.l4densvar) then if (t4dvwinlen+one_minute) & cycle read_loop else if(abs(tdiff) > twind+one_minute) cycle read_loop endif if (thin4d) then crit1 = zero else crit1 = two*abs(tdiff) ! range: 0 to 6 endif call ufbint(lnbufr,bfr2bhdr,n2bhdr,1,iret,hdr2b) satazi=bfr2bhdr(3) if (abs(satazi) > r360) then satazi=zero endif ifov = nint(bfr1bhdr(2)) lza = bfr2bhdr(1)*deg2rad ! local zenith angle if(ifov <= 65) lza=-lza ! compute look angle (panglr) and check against max angle ! panglr=(start+float(ifov-1)*step)*deg2rad ! Use this calculation for now: step = .6660465 panglr = (42.96 - float(ifov-1)*step)*deg2rad if(abs(lza)*rad2deg > MAX_SENSOR_ZENITH_ANGLE) then write(6,*)'READ_SAPHIR WARNING lza error ',lza,panglr cycle read_loop end if ! solzen_save(iob)=bfr2bhdr(2) ! encoded as 0.1E+12 in BUFR file solzen_save(iob)=zero ! set to 0.0 to bypass CRTM check solazi_save(iob)=bfr2bhdr(4) ! Read data record. Increment data counter call ufbrep(lnbufr,data1b8,1,nchanl,iret,'TMBRST') ! non-operational SAPHIR bufr: ! call ufbrep(lnbufr,data1b8,1,nchanl,iret,'TMBR') bt_save(1:nchanl,iob) = data1b8(1:nchanl) iob=iob+1 end do read_loop end do read_subset call closbf(lnbufr) deallocate(data1b8) num_obs = iob-1 500 continue if (num_obs <= 0) then write(*,*) 'READ_SAPHIR: No SAPHIR Data were read in' return end if ! Allocate arrays to hold all data for given satellite nreal = maxinfo + nstinfo nele = nreal + nchanl allocate(data_all(nele,itxmax),nrec(itxmax)) nrec=999999 ObsLoop: do iob = 1, num_obs rsat => rsat_save(iob) t4dv => t4dv_save(iob) dlon_earth => dlon_earth_save(iob) dlat_earth => dlat_earth_save(iob) crit1 => crit1_save(iob) ifov => ifov_save(iob) lza => lza_save(iob) satazi => satazi_save(iob) solzen => solzen_save(iob) solazi => solazi_save(iob) bt_in => bt_save(1:nchanl,iob) 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 disterr=acos(sin(dlat_earth)*sin(dlat00)+cos(dlat_earth)*cos(dlat00)* & (sin(dlon_earth)*sin(dlon00)+cos(dlon_earth)*cos(dlon00)))*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 (t4dvwinlen) 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 ! ! Check FOV and scan-edge usage if (.not. use_edges .and. (ifov < radedge_min .OR. ifov > radedge_max )) & 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 j=1,nchanl if (bt_in(j) < tbmin .or. bt_in(j) > tbmax) then iskip = iskip + 1 endif end do if (iskip >= nchanl) cycle ObsLoop ! Determine surface properties based on ! sst and land/sea/ice mask ! ! isflg - surface flag ! 0 sea ! 1 land ! 2 sea ice ! 3 snow ! 4 mixed ! FOV-based surface code requires fov number. if out-of-range, then ! skip this ob. if (isfcalc == 1) then call fov_check(ifov,instr,ichan,valid) if (.not. valid) cycle ObsLoop ! When isfcalc is one, calculate surface fields based on size/shape of fov. ! Otherwise, use bilinear method. call deter_sfc_fov(fov_flag,ifov,instr,ichan,satazi,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) if(isflg/=0) cycle ObsLoop ! use data over water only endif crit1 = crit1 + rlndsea(isflg) + 10._r_kind*float(iskip) + 0.01_r_kind * abs(zz) call checkob(dist1,crit1,itx,iuse) if(.not. iuse)cycle ObsLoop ! Compute "score" for observation. All scores>=0.0. Lowest score is "best" ! crit1 = crit1+pred ! not using pred for now... -ej 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 ! Re-calculate look angle ! panglr=(start+float(ifov-1)*step)*deg2rad ! Use this calculation for now: step = .6660465 panglr = (42.96 - float(ifov-1)*step)*deg2rad ! Load selected observation into data array data_all(1 ,itx)= rsat ! satellite ID data_all(2 ,itx)= t4dv ! time data_all(3 ,itx)= dlon ! grid relative longitude data_all(4 ,itx)= dlat ! grid relative latitude data_all(5 ,itx)= lza ! local zenith angle data_all(6 ,itx)= satazi ! local azimuth angle data_all(7 ,itx)= panglr ! look angle data_all(8 ,itx)= ifov ! scan position data_all(9 ,itx)= solzen ! solar zenith angle data_all(10,itx)= solazi ! solar azimuth angle 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 (deg) data_all(31,itx) = dlat_earth_deg ! earth relative latitude (deg) if(dval_use) then data_all(32,itx)= val_tovs data_all(33,itx)= itt endif 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 i=1,nchanl data_all(i+nreal,itx)=bt_in(i) end do nrec(itx)=iob end do ObsLoop ! DEAllocate I/O arrays DEALLOCATE(rsat_save) DEALLOCATE(t4dv_save) DEALLOCATE(ifov_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) call combine_radobs(mype_sub,mype_root,npe_sub,mpi_comm_sub,& nele,itxmax,nread,ndata,data_all,score_crit,nrec) if(mype_sub==mype_root)then 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(maxinfo,n)) super_val(itt)=super_val(itt)+val_tovs end do endif ! 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) end if ! Deallocate local arrays deallocate(data_all,nrec) ! Deallocate satthin arrays call destroygrids ! Deallocate FOV surface code arrays and nullify pointers. if (isfcalc == 1) call fov_cleanup if(diagnostic_reg.and.ntest>0) write(6,*)'READ_SAPHIR: ',& 'mype,ntest,disterrmax=',mype,ntest,disterrmax ! End of routine return end subroutine read_saphir