subroutine da_read_obs_bufrssmis (obstype,iv,infile) !--------------------------------------------------------------------------- ! Purpose: read in NCEP bufr SSM/IS data to innovation structure ! ! METHOD: use F90 sequential data structure to avoid reading file twice ! 1. read file radiance data in sequential data structure ! 2. do gross QC check ! 3. assign sequential data structure to innovation structure ! and deallocate sequential data structure !--------------------------------------------------------------------------- use da_control implicit none character(5) , intent (in) :: obstype ! ssmis character(20), intent (in) :: infile ! ssmis.bufr type (iv_type), intent (inout) :: iv #ifdef BUFR integer(i_kind), parameter :: bufsat_dmsp16 = 249 ! DMSP16 BUFR identifier integer(i_kind), parameter :: bufsat_dmsp17 = 285 ! DMSP17 BUFR identifier integer(i_kind), parameter :: bufsat_dmsp18 = 286 ! DMSP18 BUFR identifier integer(i_kind), parameter :: n1bhdr = 15 integer(i_kind), parameter :: maxchanl = 24 real(r_kind), parameter :: tbmin = 70.0_r_kind real(r_kind), parameter :: tbmax = 320.0_r_kind character(80) :: hdr1b !data hdr1b /'SAID FOVN YEAR MNTH DAYS HOUR MINU SECO CLAT CLON SLNM ORBN SELV SURF RAINF'/ data hdr1b /'SAID FOVN YEAR MNTH DAYS HOUR MINU SECO CLAT CLON SLNM ORBN SELV SFLG RFLAG'/ character(10) :: date character(8) :: subset, subfgn character(20) :: filename logical :: outside, outside_all integer(i_kind) :: iost, inst, lnbufr, ifgat integer(i_kind) :: num_bufr(7),numbufr,ibufr integer(i_kind) :: ihh, i, j, n, k, slnm, ifov, idd, ireadmg, ireadsb integer(i_kind) :: iret, idate, im, iy integer(i_kind) :: jc, incangl, bch, landsea_mask, rain_flag integer(i_kind) :: platform_id, satellite_id, sensor_id, nchan, num_ssmis_file integer(i_kind) :: num_ssmis_local, num_ssmis_global, num_ssmis_used, num_ssmis_thinned integer(i_kind) :: num_ssmis_used_tmp real(r_double), dimension(2,maxchanl) :: bufrtbb real(r_double), dimension(n1bhdr) :: bfr1bhdr ! pixel information integer(i_kind) :: year,month,day,hour,minute,second ! observation time real(kind=8) :: obs_time real(r_double), allocatable :: tb_inv(:) ! bright temperatures type (datalink_type), pointer :: head, p, current, prev type(info_type) :: info type(model_loc_type) :: loc ! thinning variables integer(i_kind) :: itt,itx,iobs,iout real(r_kind) :: terrain,timedif,crit,dist real(r_kind) :: dlon_earth,dlat_earth logical :: iuse real, allocatable :: in(:), out(:) logical :: found, head_found integer(i_kind), allocatable :: ptotal(:), nread(:) call da_trace_entry("da_read_obs_bufrssmis") allocate(nread(1:rtminit_nsensor)) allocate(ptotal(0:num_fgat_time)) nread(1:rtminit_nsensor) = 0 ptotal(0:num_fgat_time) = 0 platform_id = 2 ! for DMSP series sensor_id = 10 ! for SSMIS nchan = nchan_ssmis allocate (tb_inv(nchan)) num_ssmis_file = 0 num_ssmis_local = 0 num_ssmis_global = 0 num_ssmis_used = 0 num_ssmis_thinned = 0 iobs = 0 ! for thinning, argument is inout ! 0.0 Open unit to satellite bufr file and read file header !-------------------------------------------------------------- ! call da_get_unit(lnbufr) num_bufr(:)=0 numbufr=0 if (num_fgat_time>1) then do i=1,7 call da_get_unit(lnbufr) write(filename,fmt='(A,2I1,A)') trim(infile),0,i,'.bufr' open(unit = lnbufr, FILE = trim(filename),iostat = iost, form = 'unformatted', STATUS = 'OLD') if (iost == 0) then numbufr=numbufr+1 num_bufr(numbufr)=i else close (lnbufr) end if call da_free_unit(lnbufr) end do else numbufr=1 end if if (numbufr==0) numbufr=1 bufrfile: do ibufr=1,numbufr if (num_fgat_time==1) then filename=trim(infile)//'.bufr' else if ((numbufr ==1) .and. (num_bufr(ibufr) == 0)) then filename=trim(infile)//'.bufr' else write(filename,fmt='(A,2I1,A)') trim(infile),0,num_bufr(ibufr),'.bufr' end if end if lnbufr=98 open(unit=lnbufr,file=trim(filename),form='unformatted', & iostat = iost, status = 'old') if (iost /= 0) then call da_warning(__FILE__,__LINE__, & (/"Cannot open file "//infile/)) return end if call openbf(lnbufr,'IN',lnbufr) call datelen(10) call readmg(lnbufr,subset,idate,iret) iy=0 im=0 idd=0 ihh=0 write(unit=date,fmt='( i10)') idate read(unit=date,fmt='(i4,3i2)') iy,im,idd,ihh write(unit=stdout,fmt=*) & 'Bufr file date is ',iy,im,idd,ihh,infile ! Loop to read bufr file and assign information to a sequential structure !------------------------------------------------------------------------- if ( ibufr == 1 ) then allocate (head) nullify ( head % next ) p => head end if ! Set various variables depending on type of data to be read !subfgn = 'NC003003' subfgn = 'NC021201' incangl = 53.2_r_kind subset_loop: do while (ireadmg(lnbufr,subset,idate)==0) read_loop: do while (ireadsb(lnbufr)==0 .and. subset==subfgn) num_ssmis_file = num_ssmis_file + 1 ! 1.0 Read header record and data record call ufbint(lnbufr,bfr1bhdr,n1bhdr,1,iret,hdr1b) call ufbrep(lnbufr,bufrtbb,2,maxchanl,iret,"CHNM TMBR" ) ! check if observation outside range ! 2.0 Extract observation location and other required information ! QC1: judge if data is in the domain, read next record if not !------------------------------------------------------------------------ info%lat = bfr1bhdr(bufr_lat) info%lon = bfr1bhdr(bufr_lon) call da_llxy (info, loc, outside, outside_all) if (outside_all) cycle ! 3.0 Extract other information info%elv = 0.0 landsea_mask = nint(bfr1bhdr(bufr_landsea_mask)) ! ssmis surface flag ! 0:land, 2:near coast, 3:ice, ! 4:possible ice, 5:ocean, 6:coast ! RTTOV surftype: 0:land, 1:sea, 2:sea ice if ( landsea_mask == 5 ) then landsea_mask = 1 else if ( landsea_mask == 2 .or. landsea_mask == 6 ) then landsea_mask = 0 else if ( landsea_mask == 3 .or. landsea_mask == 4 ) then landsea_mask = 2 end if rain_flag = nint(bfr1bhdr(15)) ! 0:no rain, 1:rain !------------------------------------------------------ ! 3.1 Extract satellite id and scan position. if (nint(bfr1bhdr(bufr_satellite_id)) == bufsat_dmsp16) then satellite_id = 16 else if (nint(bfr1bhdr(bufr_satellite_id)) == bufsat_dmsp17) then satellite_id = 17 else if (nint(bfr1bhdr(bufr_satellite_id)) == bufsat_dmsp18) then satellite_id = 18 end if ! 3.3 Find wrfvar instrument index from RTTOV instrument triplet ! go to next data if id is not in the lists inst = 0 do i = 1, rtminit_nsensor if (platform_id == rtminit_platform(i) & .and. satellite_id == rtminit_satid(i) & .and. sensor_id == rtminit_sensor(i)) then inst = i exit end if end do if (inst == 0) cycle read_loop ! 3.1 Extract scan number and scan position. slnm = nint(bfr1bhdr(11)) ifov = nint(bfr1bhdr(bufr_ifov)) ! 3.2 Extract date information. year = bfr1bhdr(bufr_year) month = bfr1bhdr(bufr_month) day = bfr1bhdr(bufr_day) hour = bfr1bhdr(bufr_hour) minute = bfr1bhdr(bufr_minute) second = bfr1bhdr(bufr_second) write(unit=info%date_char, fmt='(i4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)') & year, '-', month, '-', day, '_', hour, ':', minute, ':', second ! QC3: time consistency check with the background date if (year <= 99) then if (year < 78) then year = year + 2000 else year = year + 1900 end if end if call da_get_julian_time(year,month,day,hour,minute,obs_time) if (obs_time < time_slots(0) .or. & obs_time >= time_slots(num_fgat_time)) cycle read_loop ! 3.2.1 determine FGAT index ifgat do ifgat=1,num_fgat_time if (obs_time >= time_slots(ifgat-1) .and. & obs_time < time_slots(ifgat)) exit end do num_ssmis_global = num_ssmis_global + 1 ptotal(ifgat) = ptotal(ifgat) + 1 if (outside) cycle ! No good for this PE num_ssmis_local = num_ssmis_local + 1 ! Make Thinning ! Map obs to thinning grid !------------------------------------------------------------------- if (thinning) then dlat_earth = info%lat dlon_earth = info%lon if (dlon_earth=r360) dlon_earth = dlon_earth-r360 dlat_earth = dlat_earth*deg2rad dlon_earth = dlon_earth*deg2rad timedif = 0. !2.0_r_kind*abs(tdiff) ! range: 0 to 6 terrain = 0.01_r_kind*abs(bfr1bhdr(13)) crit = 1. !0.01_r_kind+terrain + timedif !+ 10._r_kind*float(iskip) call map2grids(inst,ifgat,dlat_earth,dlon_earth,crit,iobs,itx,1,itt,iout,iuse) if (.not. iuse) then num_ssmis_thinned=num_ssmis_thinned+1 cycle end if end if num_ssmis_used = num_ssmis_used + 1 nread(inst) = nread(inst) + 1 if (num_ssmis_used > max_ssmis_input) then write(unit=message(1),fmt='(A,I10,A)') & "Max number of ssmis",max_ssmis_input," reached" call da_warning(__FILE__,__LINE__,message(1:1)) num_ssmis_used = num_ssmis_used - 1 exit read_loop end if ! 3.4 extract satellite and solar angle ! 3.5 extract surface information ! 3.6 extract channel bright temperature tb_inv(1:nchan) = missing_r do jc = 1, nchan bch = nint(bufrtbb(1,jc)) !ch index from bufr tb_inv(jc) = bufrtbb(2,jc) if (tb_inv(jc) < tbmin .or. tb_inv(jc) > tbmax .or. bch /= jc) then tb_inv(jc) = missing_r end if end do if ( maxval(tb_inv(:)) > missing_r ) then ! 4.0 assign information to sequential radiance structure !-------------------------------------------------------------------------- allocate (p % tb_inv (1:nchan)) p%info = info p%loc = loc p%landsea_mask = landsea_mask p%scanline = slnm p%scanpos = ifov p%satzen = incangl p%satazi = 0.0 ! dummy value p%solzen = 0.0 ! dummy value p%tb_inv(1:nchan) = tb_inv(1:nchan) p%sensor_index = inst p%ifgat = ifgat p%rain_flag = rain_flag allocate (p%next) ! add next data p => p%next nullify (p%next) else num_ssmis_local = num_ssmis_local - 1 num_ssmis_used = num_ssmis_used - 1 end if end do read_loop end do subset_loop call closbf(lnbufr) close(lnbufr) end do bufrfile if (thinning .and. num_ssmis_global > 0 ) then #ifdef DM_PARALLEL ! Get minimum crit and associated processor index. j = 0 do ifgat = 1, num_fgat_time do n = 1, iv%num_inst j = j + thinning_grid(n,ifgat)%itxmax end do end do allocate ( in (j) ) allocate ( out (j) ) j = 0 do ifgat = 1, num_fgat_time do n = 1, iv%num_inst do i = 1, thinning_grid(n,ifgat)%itxmax j = j + 1 in(j) = thinning_grid(n,ifgat)%score_crit(i) end do end do end do call mpi_reduce(in, out, j, true_mpi_real, mpi_min, root, comm, ierr) call wrf_dm_bcast_real (out, j) j = 0 do ifgat = 1, num_fgat_time do n = 1, iv%num_inst do i = 1, thinning_grid(n,ifgat)%itxmax j = j + 1 if ( ABS(out(j)-thinning_grid(n,ifgat)%score_crit(i)) > 1.0E-10 ) thinning_grid(n,ifgat)%ibest_obs(i) = 0 end do end do end do deallocate( in ) deallocate( out ) #endif ! Delete the nodes which being thinning out p => head prev => head head_found = .false. num_ssmis_used_tmp = num_ssmis_used do j = 1, num_ssmis_used_tmp n = p%sensor_index ifgat = p%ifgat found = .false. do i = 1, thinning_grid(n,ifgat)%itxmax if ( thinning_grid(n,ifgat)%ibest_obs(i) == j .and. thinning_grid(n,ifgat)%score_crit(i) < 9.99e6_r_kind ) then found = .true. exit end if end do ! free current data if ( .not. found ) then current => p p => p%next if ( head_found ) then prev%next => p else head => p prev => p end if deallocate ( current % tb_inv ) deallocate ( current ) num_ssmis_thinned = num_ssmis_thinned + 1 num_ssmis_used = num_ssmis_used - 1 nread(n) = nread(n) - 1 continue end if if ( found .and. head_found ) then prev => p p => p%next continue end if if ( found .and. .not. head_found ) then head_found = .true. head => p prev => p p => p%next end if end do end if ! End of thinning iv%total_rad_pixel = iv%total_rad_pixel + num_ssmis_used iv%total_rad_channel = iv%total_rad_channel + num_ssmis_used*nchan iv%info(radiance)%nlocal = iv%info(radiance)%nlocal + num_ssmis_used iv%info(radiance)%ntotal = iv%info(radiance)%ntotal + num_ssmis_global do i = 1, num_fgat_time ptotal(i) = ptotal(i) + ptotal(i-1) iv%info(radiance)%ptotal(i) = iv%info(radiance)%ptotal(i) + ptotal(i) end do if ( iv%info(radiance)%ptotal(num_fgat_time) /= iv%info(radiance)%ntotal ) then write(unit=message(1),fmt='(A,I10,A,I10)') & "Number of ntotal:",iv%info(radiance)%ntotal," is different from the sum of ptotal:", iv%info(radiance)%ptotal(num_fgat_time) call da_warning(__FILE__,__LINE__,message(1:1)) endif write(unit=stdout,fmt='(a)') 'num_ssmis_file, num_ssmis_global, num_ssmis_local, num_ssmis_used, num_ssmis_thinned' write(stdout,*) num_ssmis_file, num_ssmis_global, num_ssmis_local, num_ssmis_used, num_ssmis_thinned deallocate(tb_inv) ! 5.0 allocate innovation radiance structure !---------------------------------------------------------------- do i = 1, iv%num_inst if (nread(i) < 1) cycle iv%instid(i)%num_rad = nread(i) iv%instid(i)%info%nlocal = nread(i) write(UNIT=stdout,FMT='(a,i3,2x,a,3x,i10)') & 'Allocating space for radiance innov structure', & i, iv%instid(i)%rttovid_string, iv%instid(i)%num_rad call da_allocate_rad_iv (i, nchan, iv) end do ! 6.0 assign sequential structure to innovation structure !------------------------------------------------------------- nread(1:rtminit_nsensor) = 0 p => head do n = 1, num_ssmis_used i = p%sensor_index nread(i) = nread(i) + 1 call da_initialize_rad_iv (i, nread(i), iv, p) iv%instid(i)%rain_flag(n) = p%rain_flag current => p p => p%next ! free current data deallocate ( current % tb_inv ) deallocate ( current ) end do deallocate ( p ) deallocate (nread) deallocate (ptotal) call closbf(lnbufr) close(lnbufr) ! call da_free_unit(lnbufr) call da_trace_exit("da_read_obs_bufrssmis") #else call da_error(__FILE__,__LINE__,(/"Needs to be compiled with a BUFR library"/)) #endif end subroutine da_read_obs_bufrssmis