| blob | 7aa45a3d4a75db55f121b6705c10540958be58fa |
1 subroutine da_read_obs_bufriasi (obstype,iv,infile)
2 !$$$ subprogram documentation block
3 ! . . . .
4 ! subprogram: read_iasi read bufr format iasi data
6 implicit none
9 character(9) , intent (in) :: obstype
10 character(100) , intent (in) :: infile
11 type (iv_type) ,intent (inout) :: iv
12 #ifdef BUFR
13 real(kind=8) :: obs_time
14 type (datalink_type), pointer :: head, p, current, prev
15 type(info_type) :: info
16 type(model_loc_type) :: loc
17 type(model_loc_type) :: loc_fov
19 ! Subroutine constants
20 real(r_kind) :: ten = 10.0_r_kind
21 real(r_kind) :: r90 = 90.0_r_kind
22 real(r_kind) :: r360 = 360.0_r_kind
24 ! Number of channels for sensors in BUFR
25 integer(i_kind),parameter :: nchan = 616 !--- 616 subset ch out of 8078 ch for AIRS
26 integer(i_kind),parameter :: n_totchan = 616
27 integer(i_kind),parameter :: maxinfo = 33
28 integer(i_kind) :: inst,platform_id,satellite_id,sensor_id
29 integer(i_kind), allocatable :: ptotal(:), nread(:)
30 real(r_kind) :: crit
31 integer(i_kind) :: ifgat, iout, iobs
32 logical :: outside, outside_all, iuse
34 ! BUFR functions
35 integer(i_kind) :: ireadsb,ireadmg
37 ! Variables for BUFR IO
38 real(r_double),dimension(5) :: linele
39 real(r_double),dimension(14) :: allspot
40 real(r_double),dimension(2,n_totchan) :: allchan
41 real(r_double),dimension(3,10) :: cscale
42 real(r_double),dimension(6) :: cloud_frac
43 integer :: numbufr,ibufr
44 logical :: found, head_found
45 real(r_kind) :: step, start,step_adjust
46 character(len=8) :: subset
47 character(len=4) :: senname
48 character(len=80) :: allspotlist
49 integer(i_kind) :: jstart
50 integer(i_kind) :: iret
51 character(10) :: date
53 ! Work variables for time
54 integer(i_kind) :: idate, im, iy, idd, ihh
55 integer(i_kind) :: idate5(6)
57 ! Other work variables
58 real(r_kind) :: piece
59 real(r_kind) :: dlon_earth,dlat_earth
60 real(r_kind) :: lza, lzaest,sat_height_ratio
61 real(r_kind) :: sat_zenang
62 real(r_kind) :: radi
63 real(r_kind),dimension(10) :: sscale
64 real(r_kind),dimension(n_totchan) :: temperature
65 real(r_kind),allocatable,dimension(:) :: data_all
66 real(r_kind) :: Effective_Temperature
68 logical :: iasi
69 integer(i_kind) :: nreal, ksatid
70 integer(i_kind) :: ifov, iscn
71 integer(i_kind) :: num_iasi_file
72 integer(i_kind) :: num_iasi_local, num_iasi_global, num_iasi_used, num_iasi_thinned
73 integer(i_kind) :: num_iasi_used_tmp
74 integer(i_kind) :: i, j, l, iskip, ifovn, bad_line
75 integer(i_kind) :: itx, k, nele, itt, n
76 integer(i_kind) :: iexponent
77 integer(i_kind) :: num_bufr(7)
78 integer :: iost, lnbufr
79 character(20) :: filename
80 real,allocatable :: in(:), out(:)
82 ! Set standard parameters
83 integer(i_kind),parameter :: ichan=-999 ! fov-based surface code is not channel specific for iasi
84 real(r_kind),parameter :: expansion=one ! exansion factor for fov-based location code. ! use one for ir sensors.
85 real(r_kind),parameter :: tbmin = 50._r_kind
86 real(r_kind),parameter :: tbmax = 550._r_kind
87 real(r_kind),parameter :: earth_radius = 6371000._r_kind
90 if (trace_use) call da_trace_entry("da_read_obs_bufriasi")
91 ! 0.0 Initialize variables
92 !-----------------------------------
93 platform_id = 10 ! metop series
94 sensor_id = 16 !iasi
95 nreal = maxinfo
96 iasi= obstype == 'iasi'
97 bad_line=-1
98 step = 3.334
99 start = -48.33
100 step_adjust = 0.625_r_kind
101 senname = 'IASI'
102 allspotlist= &
103 'SIID YEAR MNTH DAYS HOUR MINU SECO CLATH CLONH SAZA BEARAZ SOZA SOLAZI SAID'
104 num_bufr(:)=0
105 numbufr=0
106 allocate(nread(1:rtminit_nsensor))
107 allocate(ptotal(0:num_fgat_time))
108 nread(1:rtminit_nsensor) = 0
109 ptotal(0:num_fgat_time) = 0
110 iobs = 0 ! for thinning, argument is inout
111 num_iasi_file = 0
112 num_iasi_local = 0
113 num_iasi_global = 0
114 num_iasi_used = 0
115 num_iasi_thinned = 0
116 if (num_fgat_time>1) then
117 do i=1,7
118 call da_get_unit(lnbufr)
119 write(filename,fmt='(A,2I1,A)') trim(infile),0,i,'.bufr'
120 open(unit = lnbufr, FILE = trim(filename),iostat = iost, form = 'unformatted', STATUS = 'OLD')
121 if (iost == 0) then
122 numbufr=numbufr+1
123 num_bufr(numbufr)=i
124 else
125 close (lnbufr)
126 end if
127 call da_free_unit(lnbufr)
128 end do
129 else
130 numbufr=1
131 end if
133 if (numbufr==0) numbufr=1
135 bufrfile: do ibufr=1,numbufr
136 if (num_fgat_time==1) then
137 filename=trim(infile)//'.bufr'
138 else
139 if ((numbufr ==1) .and. (num_bufr(ibufr) == 0)) then
140 filename=trim(infile)//'.bufr'
141 else
142 write(filename,fmt='(A,2I1,A)') trim(infile),0,num_bufr(ibufr),'.bufr'
143 end if
144 end if
145 lnbufr = 95
146 open(unit=lnbufr,file=trim(filename),form='unformatted', &
147 iostat = iost, status = 'old')
148 if (iost /= 0) then
149 call da_warning(__FILE__,__LINE__, &
150 (/"Cannot open file "//infile/))
151 if (trace_use) call da_trace_exit("da_read_obs_bufriasi")
152 return
153 end if
155 ! Open BUFR table
156 call openbf(lnbufr,'IN',lnbufr)
157 call datelen(10)
158 call readmg(lnbufr,subset,idate,iret)
159 iy=0
160 im=0
161 idd=0
162 ihh=0
163 write(unit=date,fmt='( i10)') idate
164 read(unit=date,fmt='(i4,3i2)') iy,im,idd,ihh
165 write(unit=stdout,fmt='(a,4i4,2x,a)') &
166 'Bufr file date is ',iy,im,idd,ihh,trim(infile)
168 ! Loop to read bufr file and assign information to a sequential structure
169 !-------------------------------------------------------------------------
170 ! Allocate arrays to hold data
171 nele=nreal+nchan
172 allocate(data_all(nele))
173 if ( ibufr == 1 ) then
174 allocate (head)
175 nullify ( head % next )
176 p => head
177 end if
180 ! Big loop to read data file
182 do while(ireadmg(lnbufr,subset,idate)>=0)
184 read_loop: do while (ireadsb(lnbufr)==0)
185 num_iasi_file = num_iasi_file + 1
187 ! Read IASI FOV information
188 call ufbint(lnbufr,linele,5,1,iret,'FOVN SLNM QGFQ MJFC SELV')
189 if ( linele(3) /= zero) then
190 cycle read_loop
191 end if
192 if ( bad_line == nint(linele(2))) then
193 ! zenith angle/scan spot mismatch, reject entire line
194 cycle read_loop
195 else
196 bad_line = -1
197 end if
198 call ufbint(lnbufr,allspot,14,1,iret,allspotlist)
199 if(iret /= 1) cycle read_loop
201 ksatid = nint(allspot(14))
202 ! SAID 3 is metop-b (metop-1)
203 ! SAID 4 is metop-a (metop-2)
204 ! SAID 5 is metop-c (metop-3)
205 if ( ( ksatid > 5) .or. ( ksatid < 3) ) then
206 write(unit=message(1),fmt='(A,I6)') 'Unknown platform: ', ksatid
207 call da_warning(__FILE__,__LINE__,message(1:1))
208 end if
209 if ( ksatid == 3 ) then
210 satellite_id = 1
211 else if ( ksatid == 4 ) then
212 satellite_id = 2
213 else if ( ksatid == 5 ) then
214 satellite_id = 3
215 end if
217 ! Check observing position
218 info%lat = allspot(8) ! latitude
219 info%lon = allspot(9) ! longitude)
220 if( abs(info%lat) > r90 .or. abs(info%lon) > r360 .or. &
221 (abs(info%lat) == r90 .and. info%lon /= zero) )then
222 write(unit=stdout,fmt=*) &
223 'READ_IASI: ### ERROR IN READING ', senname, ' BUFR DATA:', &
224 ' STRANGE OBS POINT (LAT,LON):', info%lat, info%lon
225 cycle read_loop
226 end if
228 call da_llxy (info, loc, outside, outside_all)
229 if (outside_all) cycle
230 inst = 0
231 do i = 1, rtminit_nsensor
232 if (platform_id == rtminit_platform(i) &
233 .and. satellite_id == rtminit_satid(i) &
234 .and. sensor_id == rtminit_sensor(i)) then
235 inst = i
236 exit
237 end if
238 end do
239 if (inst == 0) cycle read_loop
241 ! Check obs time
242 idate5(1) = nint(allspot(2)) ! year
243 idate5(2) = nint(allspot(3)) ! month
244 idate5(3) = nint(allspot(4)) ! day
245 idate5(4) = nint(allspot(5)) ! hour
246 idate5(5) = nint(allspot(6)) ! minute
247 idate5(6) = nint(allspot(7)) ! second
249 if( idate5(1) < 1900 .or. idate5(1) > 3000 .or. &
250 idate5(2) < 1 .or. idate5(2) > 12 .or. &
251 idate5(3) < 1 .or. idate5(3) > 31 .or. &
252 idate5(4) <0 .or. idate5(4) > 24 .or. &
253 idate5(5) <0 .or. idate5(5) > 60 )then
255 write(6,*)'READ_IASI: ### ERROR IN READING ', 'IASI', ' BUFR DATA:', &
256 ' STRANGE OBS TIME (YMDHM):', idate5(1:5)
257 cycle read_loop
258 end if
259 call da_get_julian_time(idate5(1),idate5(2),idate5(3),idate5(4),idate5(5),obs_time)
260 if ( obs_time < time_slots(0) .or. &
261 obs_time >= time_slots(num_fgat_time) ) cycle read_loop
262 do ifgat=1,num_fgat_time
263 if ( obs_time >= time_slots(ifgat-1) .and. &
264 obs_time < time_slots(ifgat) ) exit
265 end do
266 num_iasi_global = num_iasi_global + 1
267 ptotal(ifgat) = ptotal(ifgat) + 1
269 ! Observational info
270 sat_zenang = allspot(10) ! satellite zenith angle
271 ifov = nint(linele(1)) ! field of view
273 ! IASI fov ranges from 1 to 120. Current angle dependent bias
274 ! correction has a maximum of 90 scan positions. Geometry
275 ! of IASI scan allows us to remap 1-120 to 1-60. Variable
276 ! ifovn below contains the remapped IASI fov. This value is
277 ! passed on to and used in setuprad
278 ifovn = (ifov-1)/2 + 1
279 iscn = nint(linele(2)) ! scan line
281 ! Check field of view (FOVN) and satellite zenith angle (SAZA)
282 if( ifov <= 0 .or. ifov > 120 .or. sat_zenang > 90._r_kind ) then
283 write(unit=stdout,fmt=*) &
284 'READ_IASI: ### ERROR IN READING ', senname, ' BUFR DATA:', &
285 ' STRANGE OBS INFO(FOVN,SLNM,SAZA):', ifov, iscn, allspot(10)
286 cycle read_loop
287 end if
288 if ( ifov <= 60 ) sat_zenang = -sat_zenang
290 ! Compare IASI satellite scan angle and zenith angle
291 piece = -step_adjust
292 if ( mod(ifovn,2) == 1) piece = step_adjust
293 lza = ((start + float((ifov-1)/4)*step) + piece)*deg2rad
294 sat_height_ratio = (earth_radius + linele(5))/earth_radius
295 lzaest = asin(sat_height_ratio*sin(lza))*rad2deg
296 if (abs(sat_zenang - lzaest) > one) then
297 bad_line = iscn
298 write(unit=stdout,fmt=*) 'abs(sat_zenang - lzaest) > one',abs(sat_zenang - lzaest)
299 cycle read_loop
300 end if
301 ! Clear Amount (percent clear)
302 call ufbrep(lnbufr,cloud_frac,1,6,iret,'FCPH')
303 if (outside) cycle ! No good for this PE
304 num_iasi_local = num_iasi_local + 1
306 write(unit=info%date_char, &
307 fmt='(i4.4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)') &
308 idate5(1), '-', idate5(2), '-', idate5(3), '_', idate5(4), &
309 ':', idate5(5), ':', idate5(6)
310 info%elv = 0.0 !aquaspot%selv
312 ! Make Thinning
313 ! Map obs to thinning grid
314 !-------------------------------------------------------------------
315 if (thinning) then
316 dlat_earth = info%lat
317 dlon_earth = info%lon
318 if (dlon_earth<zero) dlon_earth = dlon_earth+r360
319 if (dlon_earth>=r360) dlon_earth = dlon_earth-r360
320 dlat_earth = dlat_earth*deg2rad
321 dlon_earth = dlon_earth*deg2rad
322 crit = 1.
323 call map2grids(inst,ifgat,dlat_earth,dlon_earth,crit,iobs,itx,1,itt,iout,iuse)
324 if (.not. iuse) then
325 num_iasi_thinned=num_iasi_thinned+1
326 cycle
327 end if
328 end if
329 call ufbrep(lnbufr,cscale,3,10,iret,'STCH ENCH CHSF')
330 if(iret /= 10) then
331 write(unit=stdout,fmt=*) 'READ_IASI read scale error ',iret
332 cycle read_loop
333 end if
335 ! The scaling factors are as follows, cscale(1) is the start channel number,
336 ! cscale(2) is the end channel number,
337 ! cscale(3) is the exponent scaling factor
338 ! In our case (616 channels) there are 10 groups of cscale (dimension :: cscale(3,10))
339 ! The units are W/m2..... you need to convert to mW/m2.... (subtract 5 from cscale(3)
340 do i=1,10 ! convert exponent scale factor to int and change units
341 iexponent = -(nint(cscale(3,i)) - 5)
342 sscale(i)=ten**iexponent
343 end do
345 ! Read IASI channel number(CHNM) and radiance (SCRA)
347 call ufbint(lnbufr,allchan,2,n_totchan,iret,'SCRA CHNM')
348 if( iret /= n_totchan)then
349 write(unit=stdout,fmt=*) &
350 'READ_IASI: ### ERROR IN READING ', senname, ' BUFR DATA:', &
351 iret, ' CH DATA IS READ INSTEAD OF ',n_totchan
352 cycle read_loop
353 end if
354 num_iasi_used = num_iasi_used + 1
355 nread(inst) = nread(inst) + 1
356 iskip = 0
357 jstart=1
358 do i=1,n_totchan
359 ! check that channel number is within reason
360 if (( allchan(1,i) > zero .and. allchan(1,i) < 99999._r_kind) .and. & ! radiance bounds
361 (allchan(2,i) < 8462._r_kind .and. allchan(2,i) > zero )) then ! chan # bounds
362 ! radiance to BT calculation
363 radi = allchan(1,i)
364 scaleloop: do j=jstart,10
365 if(allchan(2,i) >= cscale(1,j) .and. allchan(2,i) <= cscale(2,j))then
366 radi = allchan(1,i)*sscale(j)
367 jstart=j
368 exit scaleloop
369 end if
370 end do scaleloop
371 if (rtm_option == rtm_option_crtm) then
372 #ifdef CRTM
373 call CRTM_Planck_Temperature(inst,i,radi,temperature(i))
374 #endif
375 else if (rtm_option == rtm_option_rttov) then
376 #ifdef RTTOV
377 Effective_Temperature = LOG( ( coefs(inst)%coef%planck1(i) / radi ) + 1 )
378 Effective_Temperature = coefs(inst)%coef%planck2(i)/Effective_Temperature
379 temperature(i) = ( Effective_Temperature - coefs(inst)%coef%ff_bco(i) ) / &
380 coefs(inst)%coef%ff_bcs(i)
381 #endif
382 end if
383 if(temperature(i) < tbmin .or. temperature(i) > tbmax ) then
384 temperature(i) = min(tbmax,max(zero,temperature(i)))
386 end if
387 else ! error with channel number or radiance
388 ! write(6,*)'READ_IASI: iasi chan error',i,allchan(1,i), allchan(2,i)
389 temperature(i) = min(tbmax,max(zero,temperature(i)))
391 end if
392 end do
393 if(iskip > 0)write(6,*) ' READ_IASI : iskip > 0 ',iskip
394 ! if( iskip >= 10 )cycle read_loop
395 data_all(5) = sat_zenang ! satellite zenith angle (deg)
396 data_all(6) = allspot(11) ! satellite azimuth angle (deg)
397 data_all(9) = allspot(12) ! solar zenith angle (deg)
398 data_all(10)= allspot(13) ! solar azimuth angle (deg)
399 do l=1,nchan
400 data_all(l+nreal) = temperature(l) ! brightness temerature
401 end do
403 ! 4.0 assign information to sequential radiance structure
404 !--------------------------------------------------------------------------
405 allocate ( p % tb_inv (1:nchan) )
406 p%info = info
407 p%loc = loc
408 p%landsea_mask = 1
409 p%scanpos = ifovn
410 p%satzen = data_all(5)
411 p%satazi = data_all(6) ! look angle (deg) ! airsspot%bearaz
412 p%solzen = data_all(9)
413 p%solazi = data_all(10)
414 p%tb_inv(1:nchan) = data_all(nreal+1:nreal+nchan)
415 p%sensor_index = inst
416 p%ifgat = ifgat
418 allocate (p%next) ! add next data
419 p => p%next
420 nullify (p%next)
421 end do read_loop
422 end do
423 call closbf(lnbufr)
425 !Deallocate temporary array for next bufrfile do loop
426 deallocate(data_all)
427 end do bufrfile
429 if (thinning .and. num_iasi_global > 0 ) then
431 #ifdef DM_PARALLEL
433 ! Get minimum crit and associated processor index.
434 j = 0
435 do ifgat = 1, num_fgat_time
436 do n = 1, iv%num_inst
437 j = j + thinning_grid(n,ifgat)%itxmax
438 end do
439 end do
441 allocate ( in (j) )
442 allocate ( out (j) )
443 j = 0
444 do ifgat = 1, num_fgat_time
445 do n = 1, iv%num_inst
446 do i = 1, thinning_grid(n,ifgat)%itxmax
447 j = j + 1
448 in(j) = thinning_grid(n,ifgat)%score_crit(i)
449 end do
450 end do
451 end do
452 call mpi_reduce(in, out, j, true_mpi_real, mpi_min, root, comm, ierr)
454 call wrf_dm_bcast_real (out, j)
456 j = 0
457 do ifgat = 1, num_fgat_time
458 do n = 1, iv%num_inst
459 do i = 1, thinning_grid(n,ifgat)%itxmax
460 j = j + 1
461 if ( ABS(out(j)-thinning_grid(n,ifgat)%score_crit(i)) > 1.0E-10 ) thinning_grid(n,ifgat)%ibest_obs(i) = 0
462 end do
463 end do
464 end do
466 deallocate( in )
467 deallocate( out )
469 #endif
471 ! Delete the nodes which being thinning out
472 p => head
473 prev => head
474 head_found = .false.
475 num_iasi_used_tmp = num_iasi_used
476 do j = 1, num_iasi_used_tmp
477 n = p%sensor_index
478 ifgat = p%ifgat
479 found = .false.
481 do i = 1, thinning_grid(n,ifgat)%itxmax
482 if ( thinning_grid(n,ifgat)%ibest_obs(i) == j .and. thinning_grid(n,ifgat)%score_crit(i) < 9.99e6_r_kind ) then
483 found = .true.
484 exit
485 end if
486 end do
488 ! free current data
489 if ( .not. found ) then
490 current => p
491 p => p%next
492 if ( head_found ) then
493 prev%next => p
494 else
495 head => p
496 prev => p
497 end if
498 deallocate ( current % tb_inv )
499 deallocate ( current )
500 num_iasi_thinned = num_iasi_thinned + 1
501 num_iasi_used = num_iasi_used - 1
502 nread(n) = nread(n) - 1
503 continue
504 end if
506 if ( found .and. head_found ) then
507 prev => p
508 p => p%next
509 continue
510 end if
512 if ( found .and. .not. head_found ) then
513 head_found = .true.
514 head => p
515 prev => p
516 p => p%next
517 end if
519 end do
521 end if ! End of thinning
523 iv%total_rad_pixel = iv%total_rad_pixel + num_iasi_used
524 iv%total_rad_channel = iv%total_rad_channel + num_iasi_used*nchan
526 iv%info(radiance)%nlocal = iv%info(radiance)%nlocal + num_iasi_used
527 iv%info(radiance)%ntotal = iv%info(radiance)%ntotal + num_iasi_global
529 do i = 1, num_fgat_time
530 ptotal(i) = ptotal(i) + ptotal(i-1)
531 iv%info(radiance)%ptotal(i) = iv%info(radiance)%ptotal(i) + ptotal(i)
532 end do
533 if ( iv%info(radiance)%ptotal(num_fgat_time) /= iv%info(radiance)%ntotal ) then
534 write(unit=message(1),fmt='(A,I10,A,I10)') &
535 "Number of ntotal:",iv%info(radiance)%ntotal," is different from the sum of ptotal:", iv%info(radiance)%ptotal(num_fgat_time)
536 call da_warning(__FILE__,__LINE__,message(1:1))
537 endif
539 write(unit=message(1),fmt='(a)') ' num_iasi_file num_iasi_global num_iasi_local num_iasi_used num_iasi_thinned'
540 write(unit=message(2),fmt='(5(6x,i10))') num_iasi_file, num_iasi_global, num_iasi_local, num_iasi_used, num_iasi_thinned
541 call da_message(message(1:2))
543 ! 5.0 allocate innovation radiance structure
544 !----------------------------------------------------------------
547 do i = 1, iv%num_inst
549 if (nread(i) < 1) cycle
550 iv%instid(i)%num_rad = nread(i)
551 iv%instid(i)%info%nlocal = nread(i)
552 write(UNIT=stdout,FMT='(a,i3,2x,a,3x,i10)') &
553 'Allocating space for radiance innov structure', &
554 i, iv%instid(i)%rttovid_string, iv%instid(i)%num_rad
555 call da_allocate_rad_iv (i, nchan, iv)
557 end do
559 ! 6.0 assign sequential structure to innovation structure
560 !-------------------------------------------------------------
561 nread(1:rtminit_nsensor) = 0
562 p => head
564 do n = 1, num_iasi_used
565 i = p%sensor_index
566 nread(i) = nread(i) + 1
568 call da_initialize_rad_iv (i, nread(i), iv, p)
570 current => p
571 p => p%next
572 ! free current data
573 deallocate ( current % tb_inv )
574 deallocate ( current )
575 end do
576 deallocate ( p )
577 deallocate (nread)
578 deallocate (ptotal)
580 call closbf(lnbufr)
581 close(lnbufr)
583 ! call da_free_unit(lnbufr)
585 if (trace_use) call da_trace_exit("da_read_obs_bufriasi")
586 #else
587 call da_error(__FILE__,__LINE__,(/"Needs to be compiled with a BUFR library"/))
588 #endif
591 end subroutine da_read_obs_bufriasi