| blob | 8bb20c18e516280465978a9b93661ff25aaefc6b |
1 subroutine da_read_obs_hdf5amsr2 (iv, infile_tb,infile_clw)
2 !--------------------------------------------------------
3 ! Purpose: read in JAXA AMSR2 Level-1R data in HDF5 format
4 ! and form innovation structure
5 !
6 ! METHOD: use F90 sequantial data structure to avoid read the file twice
7 ! 1. read file radiance data in sequential data structure
8 ! 2. do gross QC check
9 ! 3. assign sequential data structure to innovation structure
10 ! and deallocate sequential data structure
11 !
12 ! HISTORY: 2013/10/22 - Creation Syed RH Rizvi, NCAR/NESL/MMM/DAS
13 ! 2014 Modification Chun Yang
14 !------------------------------------------------------------------------------
16 implicit none
18 character(len=*), intent(in) :: infile_tb, infile_clw
19 type(iv_type), intent(inout) :: iv
21 #if defined(HDF5)
22 ! fixed parameter values
23 integer,parameter::max_scan=2200 ! Maximum allowed NumberOfScans
24 integer,parameter::ovr=20 ! Number of OverlapScans
25 integer,parameter::hi_rez_fov=486 ! high resolution pixel width
26 integer,parameter::lo_rez_fov=243 ! low resolution pixel width
27 integer,parameter::time_dims=6 ! Time dimension
28 integer,parameter::nfile_max = 8 ! each hdf file contains ~50min of data
29 ! at most 8 files for a 6-h time window
30 ! interface variable
31 integer iret ! return status
32 integer(HID_T) fhnd1, fhnd2, fhnd3 ! file handle
33 integer(HID_T) ahnd1, ahnd2, ahnd3 ! attribute handle
34 integer(HID_T) dhnd1, dhnd2, dhnd3 ! dataset handle
35 integer(HID_T) shnd1, shnd2, shnd3 ! dataspace handle
36 integer(HSIZE_T) sz1(3) ! array size 1
37 integer(HSIZE_T) sz2(3) ! array size 2
39 integer(4) :: nscan ! NumberOfScans
40 real(4) :: sca ! ScaleFactor
41 real(8) :: r8d1(max_scan) ! array buffer for real(8) D1
42 integer(4) :: i4d2hi(hi_rez_fov,max_scan) ! array buffer for integer(4) D2 high
44 type AM2_COMMON_SCANTIME
45 real(8) tai93sec
46 integer(2) year
47 integer(2) month
48 integer(2) day
49 integer(2) hour
50 integer(2) minute
51 integer(2) second
52 integer(2) ms
53 integer(2) reserve
54 endtype
56 ! array data
57 type(AM2_COMMON_SCANTIME) st(max_scan) ! scantime
58 real(4) :: lat89ar(hi_rez_fov,max_scan) ! lat for 89a altitude revised
59 real(4) :: latlr(lo_rez_fov,max_scan) ! lat for low resolution
60 real(4) :: lon89ar(hi_rez_fov,max_scan) ! lon for 89a altitude revised
61 real(4) :: lonlr(lo_rez_fov,max_scan) ! lon for low resolution
63 real(4) :: tb89ah(hi_rez_fov,max_scan) ! tb for 89ah
64 real(4) :: tb89av(hi_rez_fov,max_scan) ! tb for 89av
66 integer(4) :: loflo(lo_rez_fov,max_scan*4) ! land ocean flag for low
67 integer(1) :: lof06(lo_rez_fov,max_scan) ! land ocean flag for 06
68 integer(1) :: lof10(lo_rez_fov,max_scan) ! land ocean flag for 10
69 integer(1) :: lof23(lo_rez_fov,max_scan) ! land ocean flag for 23
70 integer(1) :: lof36(lo_rez_fov,max_scan) ! land ocean flag for 36
72 integer(4) :: lofhi(hi_rez_fov,max_scan*2) ! land ocean flag for high
73 integer(1) :: lof89a(hi_rez_fov,max_scan) ! land ocean flag for 89a
74 integer(1) :: lof89b(hi_rez_fov,max_scan) ! land ocean flag for 89b
76 real(4) :: ear_in(lo_rez_fov,max_scan) ! earth incidence
77 real(4) :: ear_az(lo_rez_fov,max_scan) ! earth azimuth
78 real(4) :: clw(lo_rez_fov,max_scan) ! obs retrieved cloud liquid water
80 real(4) :: sun_az(lo_rez_fov,max_scan) ! sun_azimuth
81 real(4) :: sun_el(lo_rez_fov,max_scan) ! sun_elevation
82 real(4) :: sun_zen(lo_rez_fov,max_scan) ! sun_zenith
84 real(r_kind) :: R90 = 90.0_r_kind
85 real(r_kind),parameter :: tbmin = 50._r_kind
86 real(r_kind),parameter :: tbmax = 550._r_kind
88 real(4) :: sat_zenith(lo_rez_fov,max_scan)
89 real(4) :: sat_azimuth(lo_rez_fov,max_scan)
91 real(kind=8) :: obs_time
92 type (datalink_type),pointer :: head, p, current, prev
93 type(info_type) :: info
94 type(model_loc_type) :: loc
96 integer(i_kind) :: idate5(6)
98 integer(i_kind) :: inst,platform_id,satellite_id,sensor_id
99 real(r_kind) :: tb, crit
100 integer(i_kind) :: ifgat, iout, iobs
101 logical :: outside, outside_all, iuse
103 integer :: i,j,k,l,m,n, ifile, landsea_mask
104 logical :: found, head_found, head_allocated
106 ! Other work variables
107 real(r_kind) :: dlon_earth,dlat_earth
108 integer(i_kind) :: num_amsr2_local, num_amsr2_global, num_amsr2_used, num_amsr2_thinned
109 integer(i_kind) :: num_amsr2_used_tmp, num_amsr2_file
110 integer(i_kind) :: num_amsr2_local_local, num_amsr2_global_local, num_amsr2_file_local
111 integer(i_kind) :: itx, itt
112 character(80) :: filename1, filename2
113 integer :: nchan,ifov,iscan,ichannels
114 integer :: nfile
115 character(80) :: fname_tb(nfile_max)
116 character(80) :: fname_clw(nfile_max)
117 logical :: fexist, got_clw_file
119 ! Allocatable arrays
120 integer(i_kind),allocatable :: ptotal(:)
121 real,allocatable :: in(:), out(:)
122 real(r_kind),allocatable :: data_all(:)
124 real,allocatable :: obstime(:,:)
126 real(r_kind) :: sun_zenith, sun_azimuth
128 integer,parameter :: num_low_freq_chan=12
129 real(4) :: tb_low(lo_rez_fov,max_scan,num_low_freq_chan)
130 character(len=80) tb_low_name(num_low_freq_chan)
131 data tb_low_name/'Brightness Temperature (res06,6.9GHz,V)','Brightness Temperature (res06,6.9GHz,H)',&
132 'Brightness Temperature (res06,7.3GHz,V)','Brightness Temperature (res06,7.3GHz,H)',&
133 'Brightness Temperature (res10,10.7GHz,V)','Brightness Temperature (res10,10.7GHz,H)',&
134 'Brightness Temperature (res23,18.7GHz,V)','Brightness Temperature (res23,18.7GHz,H)',&
135 'Brightness Temperature (res23,23.8GHz,V)','Brightness Temperature (res23,23.8GHz,H)',&
136 'Brightness Temperature (res36,36.5GHz,V)','Brightness Temperature (res36,36.5GHz,H)'/
138 if (trace_use) call da_trace_entry("da_read_obs_hdf5amsr2")
140 ! 0.0 Initialize variables
141 !-----------------------------------
142 head_allocated = .false.
143 platform_id = 29 ! Table-2 Col 1 corresponding to 'gcom-w'
144 satellite_id = 1 ! Table-2 Col 3
145 sensor_id = 63 ! Table-3 Col 2 corresponding to 'amsr2'
147 allocate(ptotal(0:num_fgat_time))
148 ptotal(0:num_fgat_time) = 0
149 iobs = 0 ! for thinning, argument is inout
150 num_amsr2_file = 0
151 num_amsr2_local = 0
152 num_amsr2_global = 0
153 num_amsr2_used = 0
154 num_amsr2_thinned = 0
156 do i = 1, rtminit_nsensor
157 if (platform_id == rtminit_platform(i) &
158 .and. satellite_id == rtminit_satid(i) &
159 .and. sensor_id == rtminit_sensor(i)) then
160 inst = i
161 exit
162 end if
163 end do
164 if (inst == 0) then
165 call da_warning(__FILE__,__LINE__, &
166 (/"The combination of Satellite_Id and Sensor_Id for AMSR-2 is not found"/))
167 if (trace_use) call da_trace_exit("da_read_obs_hdf5amsr2")
168 return
169 end if
171 ! Initialize HDF5 library and Fortran90 interface
172 call H5open_f(iret)
173 if(iret.lt.0)then
174 call da_warning(__FILE__,__LINE__, &
175 (/"Problems in Initializing HDF5 library. Can not read AMSR-2 HDF5 data. "/))
176 if (trace_use) call da_trace_exit("da_read_obs_hdf5amsr2")
177 return
178 endif
180 nchan = iv%instid(inst)%nchan
181 write(unit=stdout,fmt=*)'AMSR2 nchan: ',nchan
182 allocate(data_all(1:nchan))
184 ! 1.0 Assign file names and prepare to read amsr2 files
185 !-------------------------------------------------------------------------
186 nfile = 0 !initialize
187 fname_tb(:) = '' !initialize
188 ! first check if L1SGRTBR.h5 is available
189 filename1 = trim(infile_tb)//'.h5'
190 filename2 = trim(infile_clw)//'.h5'
191 inquire (file=filename1, exist=fexist)
192 if ( fexist ) then
193 nfile = 1
194 fname_tb(nfile) = filename1
195 fname_clw(nfile) = filename2
196 else
197 ! check if L1SGRTBR-0x.h5 is available for multiple input files
198 ! here 0x is the input file sequence number
199 ! do not confuse it with fgat time slot index
200 do i = 1, nfile_max
201 write(filename1,fmt='(A,A,I2.2,A)') trim(infile_tb),'-',i,'.h5'
202 write(filename2,fmt='(A,A,I2.2,A)') trim(infile_clw),'-',i,'.h5'
203 inquire (file=filename1, exist=fexist)
204 if ( fexist ) then
205 nfile = nfile + 1
206 fname_tb(nfile) = filename1
207 fname_clw(nfile) = filename2
208 else
209 exit
210 end if
211 end do
212 end if
214 if ( nfile == 0 ) then
215 call da_warning(__FILE__,__LINE__, &
216 (/"No valid AMSR-2 L1SGRTBR.h5 or L1SGRTBR-01.h5 file found."/))
217 if (trace_use) call da_trace_exit("da_read_obs_hdf5amsr2")
218 return
219 end if
221 ! Check to see if leap second file exists for graceful failure
222 inquire( file='leapsec.dat', exist=fexist )
223 if (.not. fexist) call da_error(__FILE__,__LINE__, &
224 (/'Can not find leapsec.dat for AMSR2 data: copy or link from WRFDA/var/run'/))
226 infile_loop: do ifile = 1, nfile
227 num_amsr2_file_local = 0
228 num_amsr2_local_local = 0
229 num_amsr2_global_local = 0
231 ! open HDF5 file for read
232 call H5Fopen_f(fname_tb(ifile),H5F_ACC_RDONLY_F,fhnd1,iret,H5P_DEFAULT_F)
233 if(iret.lt.0)then
234 call da_warning(__FILE__,__LINE__, &
235 (/"Cannot open HDF5 file "//trim(fname_tb(ifile))/))
236 cycle infile_loop
237 endif
238 got_clw_file = .false.
239 call H5Fopen_f(fname_clw(ifile),H5F_ACC_RDONLY_F,fhnd2,iret,H5P_DEFAULT_F)
240 if ( iret == 0 ) then
241 got_clw_file = .true.
242 endif
243 ! to do: when got_clw_file=.true., need to check GranuleID for consistency
244 ! betweee tb and clw files
246 ! calculate NumberOfScans from array size and OverlapScans
247 call H5Dopen_f(fhnd1,'Scan Time',dhnd1,iret)
248 call H5Dget_space_f(dhnd1,shnd1,iret)
249 call H5Sget_simple_extent_dims_f(shnd1,sz1,sz2,iret)
250 if(iret.lt.0)then
251 call da_warning(__FILE__,__LINE__, &
252 (/"HDF5 read problem for: Scan Time"/))
253 endif
254 call H5Sclose_f(shnd1,iret)
255 call H5Dclose_f(dhnd1,iret)
257 nscan=sz1(1)-ovr*2
259 write(unit=stdout,fmt=*)'NumberOfScans(RETRIEVE BY ARRAY SIZE): ',nscan
260 write(unit=stdout,fmt=*)'OverlapScans(FIXED VALUE): ',ovr
262 ! check limit
263 if(nscan.gt.max_scan)then
264 write(unit=stdout,fmt=*)'limit of NumberOfScans = ',max_scan
265 call da_warning(__FILE__,__LINE__, &
266 (/"HDF5 lmit error for: max_scan"/))
267 endif
269 ! read array: scantime
270 ! read
271 call H5Dopen_f(fhnd1,'Scan Time',dhnd1,iret)
272 sz1(1)=max_scan
273 call H5Dread_f(dhnd1, &
274 H5T_NATIVE_DOUBLE,r8d1,sz1,iret,H5S_ALL_F,H5S_ALL_F)
275 if(iret.lt.0)then
276 call da_warning(__FILE__,__LINE__, &
277 (/"HDF5 read error for: Scan Time"/))
278 endif
279 call H5Dclose_f(dhnd1,iret)
280 ! cutoff overlap
281 do j=1,nscan
282 r8d1(j)=r8d1(j+ovr)
283 enddo
284 do j=nscan+1,max_scan
285 r8d1(j)=0
286 enddo
287 ! convert
288 call amsr2time(nscan,r8d1,st)
289 ! sample display
290 allocate (obstime(1:time_dims,1:nscan)) ! year, month, day, hour, min, sec
291 do j = 1, nscan
292 obstime(1,j) = st(j)%year
293 obstime(2,j) = st(j)%month
294 obstime(3,j) = st(j)%day
295 obstime(4,j) = st(j)%hour
296 obstime(5,j) = st(j)%minute
297 obstime(6,j) = st(j)%second
298 end do
299 write(unit=stdout,fmt=*)'time(scan=1) year: ',st(1)%year,' month:',st(1)%month,' day: ',st(1)%day,&
300 ' hour: ',st(1)%hour,' minute: ',st(1)%minute,' second: ',st(1)%second
302 ! read array: latlon for 89a altitude revised
303 ! read lat
304 call H5Dopen_f(fhnd1, &
305 'Latitude of Observation Point for 89A',dhnd1,iret)
306 sz1(1)=max_scan
307 sz1(2)=hi_rez_fov
308 call H5Dread_f(dhnd1, &
309 H5T_NATIVE_REAL,lat89ar,sz1,iret,H5S_ALL_F,H5S_ALL_F)
310 if(iret.lt.0)then
311 call da_warning(__FILE__,__LINE__, &
312 (/"HDF5 read error for: Latitude of Observation Point for 89A"/))
313 endif
314 call H5Dclose_f(dhnd1,iret)
316 ! read lon
317 call H5Dopen_f(fhnd1, &
318 'Longitude of Observation Point for 89A',dhnd1,iret)
319 sz1(1)=max_scan
320 sz1(2)=hi_rez_fov
321 call H5Dread_f(dhnd1, &
322 H5T_NATIVE_REAL,lon89ar,sz1,iret,H5S_ALL_F,H5S_ALL_F)
323 if(iret.lt.0)then
324 call da_warning(__FILE__,__LINE__, &
325 (/"HDF5 read error for: Longitude of Observation Point for 89A"/))
326 call da_trace_exit("da_read_obs_hdf5amsr2")
327 endif
328 call H5Dclose_f(dhnd1,iret)
330 ! cutoff overlap
331 do j=1,nscan
332 lat89ar(:,j)=lat89ar(:,j+ovr)
333 lon89ar(:,j)=lon89ar(:,j+ovr)
334 enddo
335 do j=nscan+1,max_scan
336 lat89ar(:,j)=0
337 lon89ar(:,j)=0
338 enddo
339 ! sample display
340 !write(unit=stdout,fmt=*)'latlon89ar(pixel=1,scan=1): ',lat89ar(1,1),lon89ar(1,1)
342 ! read array: latlon for low resolution
343 do j=1,nscan
344 do i=1,lo_rez_fov
345 latlr(i,j)=lat89ar(i*2-1,j)
346 lonlr(i,j)=lon89ar(i*2-1,j)
347 enddo
348 enddo
349 ! sample display
350 !write(unit=stdout,fmt=*)&
351 ! 'latlonlr(pixel=1,scan=1): ',latlr(1,1),lonlr(1,1)
353 ! read array: tb for low frequency channels
354 do k=1,num_low_freq_chan
355 call H5Dopen_f(fhnd1,tb_low_name(k),dhnd1,iret)
356 call H5Aopen_f(dhnd1,'SCALE FACTOR',ahnd1,iret)
357 call H5Aread_f(ahnd1,H5T_NATIVE_REAL,sca,sz1,iret)
358 call H5Aclose_f(ahnd1,iret)
359 sz1(1)=max_scan
360 sz1(2)=lo_rez_fov
361 call H5Dread_f(dhnd1, &
362 H5T_NATIVE_REAL,tb_low(:,:,k),sz1,iret,H5S_ALL_F,H5S_ALL_F)
363 if(iret.lt.0)then
364 call da_warning(__FILE__,__LINE__, &
365 (/"HDF5 read error for: Brightness Temperature"/))
366 endif
367 call H5Dclose_f(dhnd1,iret)
368 ! cutoff overlap & convert to unsignd & change scale
369 do j=1,nscan
370 do i=1,lo_rez_fov
371 tb_low(i,j,k)=tb_low(i,j+ovr,k)
372 if(tb_low(i,j,k).lt.65534)tb_low(i,j,k)=tb_low(i,j,k)*sca
373 enddo
374 enddo
375 do j=nscan+1,max_scan
376 tb_low(:,j,:)=0
377 enddo
378 ! sample display
379 if (print_detail_rad) then
380 write(unit=message(1),fmt='(A,I6,A,F10.4)')&
381 'tb_low(pixel=1,scan=1,chan=',k,'): ',tb_low(1,1,k)
382 call da_message(message(1:1))
383 endif
384 enddo
386 ! read array: tb for 89ah
387 ! read
388 call H5Dopen_f(fhnd1, &
389 'Brightness Temperature (original,89GHz-A,H)',dhnd1,iret)
390 call H5Aopen_f(dhnd1,'SCALE FACTOR',ahnd1,iret) ! get scale
391 call H5Aread_f(ahnd1,H5T_NATIVE_REAL,sca,sz1,iret)
392 call H5Aclose_f(ahnd1,iret)
393 sz1(1)=max_scan
394 sz1(2)=hi_rez_fov
395 call H5Dread_f(dhnd1, &
396 H5T_NATIVE_REAL,tb89ah,sz1,iret,H5S_ALL_F,H5S_ALL_F)
397 if(iret.lt.0)then
398 call da_warning(__FILE__,__LINE__, &
399 (/"HDF5 read error for: Brightness Temperature (original,89GHz-A,H)"/))
400 endif
401 call H5Dclose_f(dhnd1,iret)
402 ! cutoff overlap & convert to unsignd & change scale
403 do j=1,nscan
404 do i=1,hi_rez_fov
405 tb89ah(i,j)=tb89ah(i,j+ovr)
406 if(tb89ah(i,j).lt.65534)tb89ah(i,j)=tb89ah(i,j)*sca
407 enddo
408 enddo
409 do j=nscan+1,max_scan
410 tb89ah(:,j)=0
411 enddo
412 ! sample display
413 if (print_detail_rad) then
414 write(unit=message(1),fmt='(A,F10.4)')&
415 'tb89ah(pixel=1,scan=1): ',tb89ah(1,1)
416 call da_message(message(1:1))
417 endif
419 ! read array: tb for 89av
420 ! read
421 call H5Dopen_f(fhnd1, &
422 'Brightness Temperature (original,89GHz-A,V)',dhnd1,iret)
423 call H5Aopen_f(dhnd1,'SCALE FACTOR',ahnd1,iret) ! get scale
424 call H5Aread_f(ahnd1,H5T_NATIVE_REAL,sca,sz1,iret)
425 call H5Aclose_f(ahnd1,iret)
426 sz1(1)=max_scan
427 sz1(2)=hi_rez_fov
428 call H5Dread_f(dhnd1, &
429 H5T_NATIVE_REAL,tb89av,sz1,iret,H5S_ALL_F,H5S_ALL_F)
430 if(iret.lt.0)then
431 call da_warning(__FILE__,__LINE__, &
432 (/"HDF5 read error for: Brightness Temperature (original,89GHz-A,V)"/))
433 endif
434 call H5Dclose_f(dhnd1,iret)
435 ! cutoff overlap & convert to unsignd & change scale
436 do j=1,nscan
437 do i=1,hi_rez_fov
438 tb89av(i,j)=tb89av(i,j+ovr)
439 if(tb89av(i,j).lt.65534)tb89av(i,j)=tb89av(i,j)*sca
440 enddo
441 enddo
442 do j=nscan+1,max_scan
443 tb89av(:,j)=0
444 enddo
445 ! sample display
446 if (print_detail_rad) then
447 write(unit=message(1),fmt='(A,F10.4)')&
448 'tb89av(pixel=1,scan=1): ',tb89av(1,1)
449 call da_message(message(1:1))
450 endif
452 ! read array: land ocean flag for low
453 ! read
454 call H5Dopen_f(fhnd1, &
455 'Land_Ocean Flag 6 to 36',dhnd1,iret)
456 sz1(1)=max_scan*6
457 sz1(2)=lo_rez_fov
458 call H5Dread_f(dhnd1, &
459 H5T_NATIVE_INTEGER,loflo,sz1,iret,H5S_ALL_F,H5S_ALL_F)
460 if(iret.lt.0)then
461 call da_warning(__FILE__,__LINE__, &
462 (/"HDF5 read error for: Land_Ocean Flag 6 to 36"/))
463 endif
464 call H5Dclose_f(dhnd1,iret)
465 ! separate
466 do j=1,nscan+ovr*2
467 do i=1,lo_rez_fov
468 lof06(i,j)=loflo(i,(nscan+ovr*2)*0+j)
469 lof10(i,j)=loflo(i,(nscan+ovr*2)*1+j)
470 lof23(i,j)=loflo(i,(nscan+ovr*2)*2+j)
471 lof36(i,j)=loflo(i,(nscan+ovr*2)*3+j)
472 enddo
473 enddo
474 ! cutoff overlap
475 do j=1,nscan
476 do i=1,lo_rez_fov
477 lof06(i,j)=lof06(i,j+ovr)
478 lof10(i,j)=lof10(i,j+ovr)
479 lof23(i,j)=lof23(i,j+ovr)
480 lof36(i,j)=lof36(i,j+ovr)
481 enddo
482 enddo
483 do j=nscan+1,max_scan
484 lof06(:,j)=0
485 lof10(:,j)=0
486 lof23(:,j)=0
487 lof36(:,j)=0
488 enddo
489 ! sample display
490 !write(unit=stdout,fmt=*)'lof06(pixel=1,scan=1): ',lof06(1,1)
491 !write(unit=stdout,fmt=*)'lof10(pixel=1,scan=1): ',lof10(1,1)
492 !write(unit=stdout,fmt=*)'lof23(pixel=1,scan=1): ',lof23(1,1)
493 !write(unit=stdout,fmt=*)'lof36(pixel=1,scan=1): ',lof36(1,1)
495 ! read array: land ocean flag for high
496 ! read
497 call H5Dopen_f(fhnd1, &
498 'Land_Ocean Flag 89',dhnd1,iret)
499 sz1(1)=max_scan*2
500 sz1(2)=hi_rez_fov
501 call H5Dread_f(dhnd1, &
502 H5T_NATIVE_INTEGER,lofhi,sz1,iret,H5S_ALL_F,H5S_ALL_F)
503 if(iret.lt.0)then
504 call da_warning(__FILE__,__LINE__, &
505 (/"HDF5 read error for: Land_Ocean Flag 89"/))
506 endif
507 call H5Dclose_f(dhnd1,iret)
508 ! separate
509 do j=1,nscan+ovr*2
510 do i=1,hi_rez_fov
511 lof89a(i,j)=lofhi(i,(nscan+ovr*2)*0+j)
512 lof89b(i,j)=lofhi(i,(nscan+ovr*2)*1+j)
513 enddo
514 enddo
515 do j=1,nscan
516 do i=1,hi_rez_fov
517 lof89a(i,j)=lof89a(i,j+ovr)
518 lof89b(i,j)=lof89b(i,j+ovr)
519 enddo
520 enddo
521 do j=nscan+1,max_scan
522 lof89a(:,j)=0
523 lof89b(:,j)=0
524 enddo
525 ! sample display
526 !write(unit=stdout,fmt=*)'lof89a(pixel=1,scan=1): ',lof89a(1,1)
527 !write(unit=stdout,fmt=*)'lof89b(pixel=1,scan=1): ',lof89b(1,1)
529 ! read array: earth incidence
530 ! read
531 call H5Dopen_f(fhnd1, &
532 'Earth Incidence',dhnd1,iret)
533 call H5Aopen_f(dhnd1,'SCALE FACTOR',ahnd1,iret) ! get scale
534 call H5Aread_f(ahnd1,H5T_NATIVE_REAL,sca,sz1,iret)
535 call H5Aclose_f(ahnd1,iret)
536 sz1(1)=max_scan
537 sz1(2)=lo_rez_fov
538 call H5Dread_f(dhnd1, &
539 H5T_NATIVE_REAL,ear_in,sz1,iret,H5S_ALL_F,H5S_ALL_F)
540 if(iret.lt.0)then
541 call da_warning(__FILE__,__LINE__, &
542 (/"HDF5 read error for: Earth Incidence"/))
543 endif
544 call H5Dclose_f(dhnd1,iret)
545 ! cutoff overlap & change scale
546 do j=1,nscan
547 do i=1,lo_rez_fov
548 ear_in(i,j)=ear_in(i,j+ovr)
549 if(ear_in(i,j).gt.-32767)ear_in(i,j)=ear_in(i,j)*sca
550 enddo
551 enddo
552 do j=nscan+1,max_scan
553 ear_in(:,j)=0
554 enddo
555 ! sample display
556 !write(unit=stdout,fmt=*)'ear_in(pixel=1,scan=1): ',ear_in(1,1)
558 ! read array: earth azimuth
559 ! read
560 call H5Dopen_f(fhnd1, &
561 'Earth Azimuth',dhnd1,iret)
562 call H5Aopen_f(dhnd1,'SCALE FACTOR',ahnd1,iret) ! get scale
563 call H5Aread_f(ahnd1,H5T_NATIVE_REAL,sca,sz1,iret)
564 call H5Aclose_f(ahnd1,iret)
565 sz1(1)=max_scan
566 sz1(2)=lo_rez_fov
567 call H5Dread_f(dhnd1, &
568 H5T_NATIVE_REAL,ear_az,sz1,iret,H5S_ALL_F,H5S_ALL_F)
569 if(iret.lt.0)then
570 call da_warning(__FILE__,__LINE__, &
571 (/"HDF5 read error for: Earth Azimuth"/))
572 endif
573 call H5Dclose_f(dhnd1,iret)
574 ! cutoff overlap & change scale
575 do j=1,nscan
576 do i=1,lo_rez_fov
577 ear_az(i,j)=ear_az(i,j+ovr)
578 if(ear_az(i,j).gt.-32767)ear_az(i,j)=ear_az(i,j)*sca
579 enddo
580 enddo
581 do j=nscan+1,max_scan
582 ear_az(:,j)=0
583 enddo
584 ! sample display
585 !write(unit=stdout,fmt=*)'ear_az(pixel=1,scan=1): ',ear_az(1,1)
587 ! read array: sun azimuth
588 ! read
589 call H5Dopen_f(fhnd1, &
590 'Sun Azimuth',dhnd1,iret)
591 call H5Aopen_f(dhnd1,'SCALE FACTOR',ahnd1,iret) ! get scale
592 call H5Aread_f(ahnd1,H5T_NATIVE_REAL,sca,sz1,iret)
593 call H5Aclose_f(ahnd1,iret)
594 sz1(1)=max_scan
595 sz1(2)=lo_rez_fov
596 call H5Dread_f(dhnd1, &
597 H5T_NATIVE_REAL,sun_az,sz1,iret,H5S_ALL_F,H5S_ALL_F)
598 if(iret.lt.0)then
599 call da_warning(__FILE__,__LINE__, &
600 (/"HDF5 read error for: Sun Azimuth"/))
601 endif
602 call H5Dclose_f(dhnd1,iret)
603 ! cutoff overlap & change scale
604 do j=1,nscan
605 do i=1,lo_rez_fov
606 sun_az(i,j)=sun_az(i,j+ovr)
607 if(sun_az(i,j).gt.-32767)sun_az(i,j)=sun_az(i,j)*sca
608 enddo
609 enddo
610 do j=nscan+1,max_scan
611 sun_az(:,j)=0
612 enddo
613 ! sample display
614 !write(unit=stdout,fmt=*)'sun_az(pixel=1,scan=1): ',sun_az(1,1)
616 ! read array: sun elevation
617 ! read
618 call H5Dopen_f(fhnd1, &
619 'Sun Elevation',dhnd1,iret)
620 call H5Aopen_f(dhnd1,'SCALE FACTOR',ahnd1,iret) ! get scale
621 call H5Aread_f(ahnd1,H5T_NATIVE_REAL,sca,sz1,iret)
622 call H5Aclose_f(ahnd1,iret)
623 sz1(1)=max_scan
624 sz1(2)=lo_rez_fov
625 call H5Dread_f(dhnd1, &
626 H5T_NATIVE_REAL,sun_el,sz1,iret,H5S_ALL_F,H5S_ALL_F)
627 if(iret.lt.0)then
628 call da_warning(__FILE__,__LINE__, &
629 (/"HDF5 read error for: Sun Elevation"/))
630 endif
631 call H5Dclose_f(dhnd1,iret)
632 ! cutoff overlap & change scale
633 do j=1,nscan
634 do i=1,lo_rez_fov
635 sun_el(i,j)=sun_el(i,j+ovr)
636 if(sun_el(i,j).gt.-32767)sun_el(i,j)=sun_el(i,j)*sca
637 enddo
638 enddo
639 do j=nscan+1,max_scan
640 sun_el(:,j)=0
641 enddo
642 ! sample display
643 !write(unit=stdout,fmt=*)'sun_el(pixel=1,scan=1): ',sun_el(1,1)
644 sun_zen(:,:)=R90-sun_el(:,:)
645 sat_zenith(:,:)=ear_in(:,:)
646 sat_azimuth(:,:)=ear_az(:,:)
648 ! close file and HDF5
649 call H5Fclose_f(fhnd1,iret)
651 if ( got_clw_file ) then
652 ! read CLW from infile_clw:
653 call H5Dopen_f(fhnd2,'Geophysical Data',dhnd2,iret)
654 call H5Aopen_f(dhnd2,'SCALE FACTOR',ahnd2,iret)
655 call H5Aread_f(ahnd2,H5T_NATIVE_REAL,sca,sz1,iret)
656 call H5Aclose_f(ahnd2,iret)
657 sz1(1)=max_scan
658 sz1(2)=lo_rez_fov
659 call H5Dread_f(dhnd2, &
660 H5T_NATIVE_REAL,clw,sz1,iret,H5S_ALL_F,H5S_ALL_F)
661 if(iret.lt.0)then
662 call da_warning(__FILE__,__LINE__, &
663 (/"HDF5 read error for: CLW data"/))
664 endif
665 call H5Dclose_f(dhnd2,iret)
666 ! change scale
667 do j=1,nscan
668 do i=1,lo_rez_fov
669 if(clw(i,j).gt.-32767)clw(i,j)=clw(i,j)*sca
670 enddo
671 enddo
672 do j=nscan+1,max_scan
673 clw(:,j)=0
674 enddo
675 ! sample display
676 !write(unit=stdout,fmt=*)'clw(pixel=1,scan=1): ',clw(1,1)
677 ! close file and HDF5
678 call H5Fclose_f(fhnd2,iret)
679 end if
681 ! 2.0 Loop to read hdf file and assign information to a sequential structure
682 !-------------------------------------------------------------------------
684 ! Allocate arrays to hold data
685 if ( .not. head_allocated ) then
686 allocate (head)
687 nullify ( head % next )
688 p => head
689 head_allocated = .true.
690 end if
691 ! start scan_loop
692 scan_loop: do iscan=1, nscan
693 do i = 1, 6
694 idate5(i)=obstime(i, iscan)
695 end do
696 call da_get_julian_time(idate5(1),idate5(2),idate5(3),idate5(4),idate5(5),obs_time)
697 if ( obs_time < time_slots(0) .or. &
698 obs_time >= time_slots(num_fgat_time) ) cycle scan_loop
699 do ifgat=1,num_fgat_time
700 if ( obs_time >= time_slots(ifgat-1) .and. &
701 obs_time < time_slots(ifgat) ) exit
702 end do
704 ! start fov_loop
705 fov_loop: do ifov=1, lo_rez_fov
706 num_amsr2_file = num_amsr2_file + 1
707 num_amsr2_file_local = num_amsr2_file_local + 1
708 info%lat = latlr(ifov,iscan)
709 info%lon = lonlr(ifov,iscan)
711 call da_llxy (info, loc, outside, outside_all)
712 if (outside_all) cycle fov_loop
714 num_amsr2_global = num_amsr2_global + 1
715 num_amsr2_global_local = num_amsr2_global_local + 1
716 ptotal(ifgat) = ptotal(ifgat) + 1
717 if (outside) cycle fov_loop ! No good for this PE
718 ! Discard data over Land (landmask =0 -->Land =1 -->Sea)
719 landsea_mask = 0
720 if(lof06(ifov,iscan) < 1.0 .and. lof10(ifov,iscan) < 1.0 .and. &
721 lof23(ifov,iscan) < 1.0 .and. lof36(ifov,iscan) < 1.0 .and. &
722 lof89a(2*ifov-1,iscan) < 1.0 ) landsea_mask = 1
723 if( landsea_mask == 0 ) cycle fov_loop
725 num_amsr2_local = num_amsr2_local + 1
726 num_amsr2_local_local = num_amsr2_local_local + 1
727 write(unit=info%date_char, &
728 fmt='(i4.4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)') &
729 idate5(1), '-', idate5(2), '-', idate5(3), '_', idate5(4), &
730 ':', idate5(5), ':', idate5(6)
731 info%elv = 0.0
733 ! 3.0 Make Thinning
734 ! Map obs to thinning grid
735 !-------------------------------------------------------------------
736 if (thinning) then
737 dlat_earth = info%lat !degree
738 dlon_earth = info%lon
739 if (dlon_earth<zero) dlon_earth = dlon_earth+r360
740 if (dlon_earth>=r360) dlon_earth = dlon_earth-r360
741 dlat_earth = dlat_earth*deg2rad !radian
742 dlon_earth = dlon_earth*deg2rad
743 crit = 1.
744 call map2grids(inst,ifgat,dlat_earth,dlon_earth,crit,iobs,itx,1,itt,iout,iuse)
745 if (.not. iuse) then
746 num_amsr2_thinned = num_amsr2_thinned+1
747 cycle fov_loop
748 end if
749 end if
751 num_amsr2_used = num_amsr2_used + 1
752 data_all = missing_r
754 do k=1,num_low_freq_chan
755 tb = tb_low(ifov,iscan,k)
756 if( tb < tbmin .or. tb > tbmax ) tb = missing_r
757 data_all(k)= tb
758 enddo
760 tb = tb89av(2*ifov-1,iscan)
761 if( tb < tbmin .or. tb > tbmax ) tb = missing_r
762 data_all(13)= tb
764 tb = tb89ah(2*ifov-1,iscan)
765 if( tb < tbmin .or. tb > tbmax ) tb = missing_r
766 data_all(14)= tb
768 ! 4.0 assign information to sequential radiance structure
769 !--------------------------------------------------------------------------
770 allocate ( p % tb_inv (1:nchan ))
771 p%info = info
772 p%loc = loc
773 p%landsea_mask = landsea_mask
774 p%scanpos = ifov
775 p%satzen = sat_zenith(ifov,iscan)
776 p%satazi = sat_azimuth(ifov,iscan)
777 p%solzen = sun_zen(ifov,iscan)
778 p%solazi = sun_az(ifov,iscan)
779 p%clw = clw(ifov,iscan)
780 p%tb_inv(1:nchan) = data_all(1:nchan)
781 p%sensor_index = inst
782 p%ifgat = ifgat
784 allocate (p%next) ! add next data
785 p => p%next
786 nullify (p%next)
787 end do fov_loop
788 end do scan_loop
790 ! Dellocate arrays
791 deallocate (obstime)
793 write(stdout,fmt='(3a,i7)') ' In file: ',trim(fname_tb(ifile)),' got num_amsr2_file : ',num_amsr2_file_local
794 write(stdout,fmt='(3a,i7)') ' In file: ',trim(fname_tb(ifile)),' got num_amsr2_global : ',num_amsr2_global_local
795 write(stdout,fmt='(3a,i7)') ' In file: ',trim(fname_tb(ifile)),' got num_amsr2_local : ',num_amsr2_local_local
796 end do infile_loop
798 call H5close_f(iret)
800 deallocate(data_all) ! Deallocate data arrays
802 if (thinning .and. num_amsr2_global > 0 ) then
803 #ifdef DM_PARALLEL
804 ! Get minimum crit and associated processor index.
805 j = 0
806 do ifgat = 1, num_fgat_time
807 j = j + thinning_grid(inst,ifgat)%itxmax
808 end do
810 allocate ( in (j) )
811 allocate ( out (j) )
812 j = 0
813 do ifgat = 1, num_fgat_time
814 do i = 1, thinning_grid(inst,ifgat)%itxmax
815 j = j + 1
816 in(j) = thinning_grid(inst,ifgat)%score_crit(i)
817 end do
818 end do
819 call mpi_reduce(in, out, j, true_mpi_real, mpi_min, root, comm, ierr)
821 call wrf_dm_bcast_real (out, j)
823 j = 0
824 do ifgat = 1, num_fgat_time
825 do i = 1, thinning_grid(inst,ifgat)%itxmax
826 j = j + 1
827 if ( ABS(out(j)-thinning_grid(inst,ifgat)%score_crit(i)) > 1.0E-10 ) &
828 thinning_grid(inst,ifgat)%ibest_obs(i) = 0
829 end do
830 end do
832 deallocate( in )
833 deallocate( out )
835 #endif
837 ! Delete the nodes which being thinning out
838 p => head
839 prev => head
840 head_found = .false.
841 num_amsr2_used_tmp = num_amsr2_used
842 do j = 1, num_amsr2_used_tmp
843 n = p%sensor_index
844 ifgat = p%ifgat
845 found = .false.
847 do i = 1, thinning_grid(n,ifgat)%itxmax
848 if ( thinning_grid(n,ifgat)%ibest_obs(i) == j .and. thinning_grid(n,ifgat)%score_crit(i) < 9.99e6_r_kind ) then
849 found = .true.
850 exit
851 end if
852 end do
854 ! free current data
855 if ( .not. found ) then
856 current => p
857 p => p%next
858 if ( head_found ) then
859 prev%next => p
860 else
861 head => p
862 prev => p
863 end if
864 deallocate ( current % tb_inv )
865 deallocate ( current )
866 num_amsr2_thinned = num_amsr2_thinned + 1
867 num_amsr2_used = num_amsr2_used - 1
868 continue
869 end if
871 if ( found .and. head_found ) then
872 prev => p
873 p => p%next
874 continue
875 end if
877 if ( found .and. .not. head_found ) then
878 head_found = .true.
879 head => p
880 prev => p
881 p => p%next
882 end if
884 end do
886 end if ! End of thinning
888 iv%total_rad_pixel = iv%total_rad_pixel + num_amsr2_used
889 iv%total_rad_channel = iv%total_rad_channel + num_amsr2_used*nchan
891 iv%info(radiance)%nlocal = iv%info(radiance)%nlocal + num_amsr2_used
892 iv%info(radiance)%ntotal = iv%info(radiance)%ntotal + num_amsr2_global
894 do i = 1, num_fgat_time
895 ptotal(i) = ptotal(i) + ptotal(i-1)
896 iv%info(radiance)%ptotal(i) = iv%info(radiance)%ptotal(i) + ptotal(i)
897 end do
898 if ( iv%info(radiance)%ptotal(num_fgat_time) /= iv%info(radiance)%ntotal ) then
899 write(unit=message(1),fmt='(A,I10,A,I10)') &
900 "Number of ntotal:",iv%info(radiance)%ntotal," is different from the sum of ptotal:", iv%info(radiance)%ptotal(num_fgat_time)
901 call da_warning(__FILE__,__LINE__,message(1:1))
902 endif
904 write(unit=stdout,fmt='(a)') 'AMSR2 data counts: '
905 write(stdout,fmt='(a,i7)') ' In file: ',num_amsr2_file
906 write(stdout,fmt='(a,i7)') ' Global : ',num_amsr2_global
907 write(stdout,fmt='(a,i7)') ' Local : ',num_amsr2_local
908 write(stdout,fmt='(a,i7)') ' Used : ',num_amsr2_used
909 write(stdout,fmt='(a,i7)') ' Thinned: ',num_amsr2_thinned
911 ! 5.0 allocate innovation radiance structure
912 !----------------------------------------------------------------
914 if (num_amsr2_used > 0) then
915 iv%instid(inst)%num_rad = num_amsr2_used
916 iv%instid(inst)%info%nlocal = num_amsr2_used
917 write(UNIT=stdout,FMT='(a,i3,2x,a,3x,i10)') &
918 'Allocating space for radiance innov structure', &
919 inst, iv%instid(inst)%rttovid_string, iv%instid(inst)%num_rad
920 call da_allocate_rad_iv (inst, nchan, iv)
921 end if
923 ! 6.0 assign sequential structure to innovation structure
924 !-------------------------------------------------------------
925 p => head
927 do n = 1, num_amsr2_used
928 i = p%sensor_index
929 call da_initialize_rad_iv (i, n, iv, p)
930 current => p
931 p => p%next
932 ! free current data
933 deallocate ( current % tb_inv )
934 deallocate ( current )
935 end do
936 deallocate ( p )
937 deallocate (ptotal)
939 if (trace_use) call da_trace_exit("da_read_obs_hdf5amsr2")
940 #else
941 call da_error(__FILE__,__LINE__,(/"Needs to be compiled with HDF5 library"/))
942 #endif
943 end subroutine da_read_obs_hdf5amsr2