| blob | fbcdebf8b4c07963dd9b30c7b7fa3ce07d3eccff |
1 subroutine da_read_obs_hdf5mwhs2 (iv, infile_tb)
2 !--------------------------------------------------------
3 ! Purpose: read in FY3C MWHS2 Level-1 data in HDF5 format
4 ! and form innovation structure
5 !
6 ! METHOD: use F90 sequential data structure to avoid reading 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: 2017/10/15 - Creation Wei sun
13 !------------------------------------------------------------------------------
15 implicit none
17 character(len=*), intent(in) :: infile_tb
18 type(iv_type), intent(inout) :: iv
20 #if defined(HDF5)
21 ! fixed parameter values
22 integer,parameter::max_scan=2500 ! Maximum allowed NumberOfScans
23 integer,parameter::ovr=0 ! Number of OverlapScans
24 integer,parameter::hi_rez_fov=486 ! high resolution pixel width
25 integer,parameter::lo_rez_fov=98 ! low resolution pixel width
26 integer,parameter::time_dims=6 ! Time dimension
27 integer,parameter::nfile_max = 4 ! each hdf file contains ~100min of data
28 ! at most 4 files for a 6-h time window
29 integer,parameter :: num_low_freq_chan=15
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) :: slp,itp ! 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 integer :: yr,mt,dy,hr,mn,sc
45 integer :: tess,tesp,allocstat
46 integer :: scnday(max_scan),scntim(max_scan)
47 real(8) :: ju
48 integer,allocatable :: date(:,:)
49 real(4),allocatable :: tbb(:,:,:)
51 ! array data
52 real(4) :: latlr(lo_rez_fov,max_scan) ! lat for low resolution
53 real(4) :: lonlr(lo_rez_fov,max_scan) ! lon for low resolution
55 real(4) :: tb89ah(hi_rez_fov,max_scan) ! tb for 89ah
56 real(4) :: tb89av(hi_rez_fov,max_scan) ! tb for 89av
58 ! real(4) :: tbb(98,2385,15)
60 integer(4) :: loflo(lo_rez_fov,max_scan) ! land ocean flag for low
61 integer(4) :: elev(lo_rez_fov,max_scan) ! elevation
63 integer :: satzen0(lo_rez_fov,max_scan) ! satellite zenith
64 integer :: satazi0(lo_rez_fov,max_scan) ! satellite azimuth
65 integer :: solzen0(lo_rez_fov,max_scan) ! sun zenith
66 integer :: solazi0(lo_rez_fov,max_scan) ! sun azimuth
68 real :: satzen(lo_rez_fov,max_scan) ! satellite zenith
69 real :: satazi(lo_rez_fov,max_scan) ! satellite azimuth
70 real :: solzen(lo_rez_fov,max_scan) ! sun zenith
71 real :: solazi(lo_rez_fov,max_scan) ! sun azimuth
74 real(r_kind) :: R90 = 90.0_r_kind
75 real(r_kind),parameter :: tbmin = 50._r_kind
76 real(r_kind),parameter :: tbmax = 550._r_kind
79 real(kind=8) :: obs_time
80 type (datalink_type),pointer :: head, p, current, prev
81 type(info_type) :: info
82 type(model_loc_type) :: loc
84 integer(i_kind) :: idate5(6)
86 integer(i_kind) :: inst,platform_id,satellite_id,sensor_id
87 real(r_kind) :: tb, crit
88 integer(i_kind) :: ifgat, iout, iobs
89 logical :: outside, outside_all, iuse
91 integer :: i,j,k,l,m,n, ifile, landsea_mask
92 logical :: found, head_found, head_allocated
94 ! Other work variables
95 real(r_kind) :: dlon_earth,dlat_earth
96 integer(i_kind) :: num_mwhs2_local, num_mwhs2_global, num_mwhs2_used, num_mwhs2_thinned
97 integer(i_kind) :: num_mwhs2_used_tmp, num_mwhs2_file
98 integer(i_kind) :: num_mwhs2_local_local, num_mwhs2_global_local, num_mwhs2_file_local
99 integer(i_kind) :: itx, itt
100 character(80) :: filename1, filename2
101 integer :: nchan,ifov,iscan,ichannels
102 integer :: nfile
103 character(80) :: fname_tb(nfile_max)
104 logical :: fexist
106 ! Allocatable arrays
107 integer(i_kind),allocatable :: ptotal(:)
108 real,allocatable :: in(:), out(:)
111 !real(kind=8) :: tb_low(lo_rez_fov,max_scan,num_low_freq_chan)
113 if (trace_use) call da_trace_entry("da_read_obs_hdf5mwhs2")
115 ! 0.0 Initialize variables
116 !-----------------------------------
117 head_allocated = .false.
118 platform_id = 23 ! Table-2 Col 1 corresponding to 'fy3'
119 satellite_id = 3 ! Table-2 Col 3
120 sensor_id = 73 ! Table-3 Col 2 corresponding to 'mwhs2'
122 allocate(ptotal(0:num_fgat_time))
123 ptotal(0:num_fgat_time) = 0
124 iobs = 0 ! for thinning, argument is inout
125 num_mwhs2_file = 0
126 num_mwhs2_local = 0
127 num_mwhs2_global = 0
128 num_mwhs2_used = 0
129 num_mwhs2_thinned = 0
131 do i = 1, rtminit_nsensor
132 if (platform_id == rtminit_platform(i) &
133 .and. satellite_id == rtminit_satid(i) &
134 .and. sensor_id == rtminit_sensor(i)) then
135 inst = i
136 exit
137 end if
138 end do
139 if (inst == 0) then
140 call da_warning(__FILE__,__LINE__, &
141 (/"The combination of Satellite_Id and Sensor_Id for MWHS-2 is not found"/))
142 if (trace_use) call da_trace_exit("da_read_obs_hdf5mwhs2")
143 return
144 end if
146 ! Initialize HDF5 library and Fortran90 interface
147 call H5open_f(iret)
148 if(iret.lt.0)then
149 call da_warning(__FILE__,__LINE__, &
150 (/"Problems in Initializing HDF5 library. Can not read MWHS-2 HDF data. "/))
151 if (trace_use) call da_trace_exit("da_read_obs_hdf5mwhs2")
152 return
153 endif
155 nchan = iv%instid(inst)%nchan
156 write(unit=stdout,fmt=*)'MWHS2 nchan: ',nchan
158 ! 1.0 Assign file names and prepare to read mwhs2 files
159 !-------------------------------------------------------------------------
160 nfile = 0 !initialize
161 fname_tb(:) = '' !initialize
162 ! first check if L1SGRTBR.h5 is available
163 filename1 = trim(infile_tb)//'.hdf'
164 inquire (file=filename1, exist=fexist)
166 if ( fexist ) then
167 nfile = 1
168 fname_tb(nfile) = filename1
169 else
170 ! check if L1SGRTBR-0x.h5 is available for multiple input files
171 ! here 0x is the input file sequence number
172 ! do not confuse it with fgat time slot index
173 do i = 1, nfile_max
174 write(filename1,fmt='(A,A,I2.2,A)') trim(infile_tb),'-',i,'.hdf'
175 inquire (file=filename1, exist=fexist)
176 if ( fexist ) then
177 nfile = nfile + 1
178 fname_tb(nfile) = filename1
179 else
180 exit
181 end if
182 end do
183 end if
185 if ( nfile == 0 ) then
186 call da_warning(__FILE__,__LINE__, &
187 (/"No valid MWHS-2 .hdf file found."/))
188 if (trace_use) call da_trace_exit("da_read_obs_hdf5mwhs2")
189 return
190 end if
192 !! ! Check to see if leap second file exists for graceful failure
193 !! inquire( file='leapsec.dat', exist=fexist )
194 !! if (.not. fexist) call da_error(__FILE__,__LINE__, &
195 !! (/'Can not find leapsec.dat for MWHS2 data: copy or link from WRFDA/var/run'/))
197 infile_loop: do ifile = 1, nfile
198 num_mwhs2_file_local = 0
199 num_mwhs2_local_local = 0
200 num_mwhs2_global_local = 0
202 ! open HDF5 file for read
203 call H5Fopen_f(fname_tb(ifile),H5F_ACC_RDONLY_F,fhnd1,iret,H5P_DEFAULT_F)
204 if(iret.lt.0)then
205 call da_warning(__FILE__,__LINE__, &
206 (/"Cannot open HDF5 file "//trim(fname_tb(ifile))/))
207 cycle infile_loop
208 endif
210 ! calculate NumberOfScans from array size and OverlapScans
211 call h5dopen_f(fhnd1,"/Geolocation/Latitude",dhnd1,iret)
212 call h5dget_space_f(dhnd1,shnd1,iret)
213 call h5sget_simple_extent_dims_f(shnd1,sz1,sz2,iret)
214 if(iret.lt.0)then
215 call da_warning(__FILE__,__LINE__, &
216 (/"HDF5 read problem for: Latitude"/))
217 endif
218 call H5Sclose_f(shnd1,iret)
219 call H5Dclose_f(dhnd1,iret)
220 nscan=sz1(2)
222 write(unit=stdout,fmt=*)'NumberOfScans(RETRIEVE BY ARRAY SIZE): ',nscan
223 write(unit=stdout,fmt=*)'OverlapScans(FIXED VALUE): ',ovr
225 ! check limit
226 if(nscan.gt.max_scan)then
227 write(unit=stdout,fmt=*)'limit of NumberOfScans = ',max_scan
228 call da_warning(__FILE__,__LINE__, &
229 (/"HDF5 lmit error for: max_scan"/))
230 endif
232 ! read array: scantime
233 ! read
234 sz1(1)=max_scan
235 call h5dopen_f(fhnd1,"/Geolocation/Scnlin_daycnt",dhnd1,iret)
236 call h5dread_f(dhnd1,H5T_NATIVE_INTEGER,scnday,sz1,iret)
237 if(iret.lt.0)then
238 call da_warning(__FILE__,__LINE__, &
239 (/"HDF5 read error for: Scan Day"/))
240 endif
241 call h5dclose_f(dhnd1, iret)
242 do j=nscan+1,max_scan
243 scnday(j)=0
244 end do
246 call h5dopen_f(fhnd1,"/Geolocation/Scnlin_mscnt",dhnd1,iret)
247 call h5dread_f(dhnd1,H5T_NATIVE_INTEGER,scntim,sz1,iret)
248 if(iret.lt.0)then
249 call da_warning(__FILE__,__LINE__, &
250 (/"HDF5 read error for: Scan Time"/))
251 endif
252 call h5dclose_f(dhnd1,iret)
254 do j=nscan+1,max_scan
255 scntim(j)=0
256 end do
258 allocate (date(nscan,6))
259 do j=1,nscan
260 ju=scnday(j)+dfloat(scntim(j))/86400000.
261 call mjd2cal(ju,yr,mt,dy,hr,mn,sc)
262 date(j,1)=yr
263 date(j,2)=mt
264 date(j,3)=dy
265 date(j,4)=hr
266 date(j,5)=mn
267 date(j,6)=sc
268 end do
270 tess = 2046
271 tesp = 47
272 ! sample display
273 write(unit=stdout,fmt=*)'time=',tesp,tess,date(tess,:)
276 ! read array: tb
277 ! read
279 ! sz1(1)=98
280 ! sz1(2)=2385
281 ! sz1(3)=15
283 sz1(1)=lo_rez_fov
284 sz1(2)=nscan
285 sz1(3)=num_low_freq_chan
288 ! allocate(tbb(lo_rez_fov,nscan,num_low_freq_chan))
289 allocate(tbb(sz1(1),sz1(2),sz1(3)))
291 print*,'newthree',lo_rez_fov,nscan,num_low_freq_chan
293 call h5dopen_f(fhnd1,"/Data/Earth_Obs_BT",dhnd1,iret)
294 call h5dread_f(dhnd1,H5T_NATIVE_REAL,tbb,sz1,iret)
296 if(iret.lt.0)then
297 call da_warning(__FILE__,__LINE__, &
298 (/"HDF5 read error for: Brightness Temperature"/))
299 endif
300 call h5dclose_f(dhnd1, iret)
302 ! cutoff overlap & convert to unsignd & change scale
303 !do j=nscan+1,max_scan
304 ! tb_low(:,j,:)=0
305 !enddo
306 ! sample display
307 ! if (print_detail_rad) then
308 ! write(unit=message(1),fmt='(A,F10.4)')&
309 ! 'tb(pixel=1,scan=1): ',tb_low(1,1,:)
310 ! call da_message(message(1:1))
311 ! endif
312 write(unit=stdout,fmt=*)'tb',tesp,tess,tbb(tesp,tess,:)
314 ! read array: latlon
315 ! read lat
316 sz1(1)=lo_rez_fov
317 sz1(2)=max_scan
318 call h5dopen_f(fhnd1,"/Geolocation/Latitude",dhnd1,iret)
319 call h5dread_f(dhnd1,H5T_NATIVE_REAL,latlr,sz1,iret)
320 if(iret.lt.0)then
321 call da_warning(__FILE__,__LINE__, &
322 (/"HDF5 read error for: Latitude"/))
323 endif
324 call h5dclose_f(dhnd1, iret)
326 ! read lon
327 sz1(1)=lo_rez_fov
328 sz1(2)=max_scan
329 call h5dopen_f(fhnd1,"/Geolocation/Longitude",dhnd1,iret)
330 call h5dread_f(dhnd1,H5T_NATIVE_REAL,lonlr,sz1,iret)
331 if(iret.lt.0)then
332 call da_warning(__FILE__,__LINE__, &
333 (/"HDF5 read error for: Longitude"/))
334 call da_trace_exit("da_read_obs_hdf5mwhs2")
335 endif
336 call h5dclose_f(dhnd1, iret)
338 ! cutoff overlap
339 do i=1,lo_rez_fov
340 do j=1,nscan
341 if(lonlr(i,j).lt.0.0) lonlr(i,j)=lonlr(i,j)+360.0
342 enddo
343 enddo
344 do j=nscan+1,max_scan
345 latlr(:,j)=0
346 lonlr(:,j)=0
347 enddo
348 ! sample display
350 write(unit=stdout,fmt=*)'lat',latlr(tesp,tess:tess+10)
351 write(unit=stdout,fmt=*)'lon',lonlr(tesp,tess:tess+10)
353 ! read array: elevation
354 ! read
355 sz1(1)=lo_rez_fov
356 sz1(2)=max_scan
357 call h5dopen_f(fhnd1,"/Data/DEM",dhnd1,iret)
358 call h5dread_f(dhnd1,H5T_NATIVE_INTEGER,elev,sz1,iret)
359 if(iret.lt.0)then
360 call da_warning(__FILE__,__LINE__, &
361 (/"HDF5 read error for: Elevation"/))
362 endif
363 call h5dclose_f(dhnd1,iret)
364 ! cutoff overlap
365 do j=nscan+1,max_scan
366 elev(:,j)=0
367 enddo
368 ! sample display
369 write(unit=stdout,fmt=*)'elev',elev(tesp,tess)
371 ! read array: land ocean flag for low
372 ! read
373 sz1(1)=lo_rez_fov
374 sz1(2)=max_scan
375 call h5dopen_f(fhnd1,"/Data/LandSeaMask",dhnd1,iret)
376 call h5dread_f(dhnd1,H5T_NATIVE_INTEGER,loflo,sz1,iret)
377 if(iret.lt.0)then
378 call da_warning(__FILE__,__LINE__, &
379 (/"HDF5 read error for: Land_Ocean Flag"/))
380 endif
381 call h5dclose_f(dhnd1,iret)
382 ! cutoff overlap
383 do j=nscan+1,max_scan
384 loflo(:,j)=0
385 enddo
386 ! sample display
387 write(unit=stdout,fmt=*)'lof',loflo(tesp,tess)
389 ! read array: satellite zenith
390 ! read
391 sz1(1)=lo_rez_fov
392 sz1(2)=max_scan
393 call h5dopen_f(fhnd1,"/Geolocation/SensorZenith",dhnd1,iret)
394 call h5dread_f(dhnd1,H5T_NATIVE_INTEGER,satzen0,sz1,iret)
395 call h5aopen_f(dhnd1,"Slope",ahnd1,iret)
396 call h5aread_f(ahnd1,H5T_NATIVE_REAL,slp,sz1,iret)
397 call h5aclose_f(ahnd1,iret)
398 call h5aopen_f(dhnd1,"Intercept",ahnd1,iret)
399 call h5aread_f(ahnd1,H5T_NATIVE_REAL,itp,sz1,iret)
400 call h5aclose_f(ahnd1,iret)
401 if(iret.lt.0)then
402 call da_warning(__FILE__,__LINE__, &
403 (/"HDF5 read error for: Satellite Zenith"/))
404 endif
405 call h5dclose_f(dhnd1,iret)
406 satzen=satzen0*slp+itp
407 ! cutoff overlap & change scale
408 do j=nscan+1,max_scan
409 satzen(:,j)=0
410 enddo
411 ! sample display
412 write(unit=stdout,fmt=*)'sat_zen',satzen(tesp,tess)
414 ! read array: satellite azimuth
415 ! read
416 sz1(1)=lo_rez_fov
417 sz1(2)=max_scan
418 call h5dopen_f(fhnd1,"/Geolocation/SensorAzimuth",dhnd1,iret)
419 call h5dread_f(dhnd1,H5T_NATIVE_INTEGER,satazi0,sz1,iret)
420 call h5aopen_f(dhnd1,"Slope",ahnd1,iret)
421 call h5aread_f(ahnd1,H5T_NATIVE_REAL,slp,sz1,iret)
422 call h5aclose_f(ahnd1,iret)
423 call h5aopen_f(dhnd1,"Intercept",ahnd1,iret)
424 call h5aread_f(ahnd1,H5T_NATIVE_REAL,itp,sz1,iret)
425 call h5aclose_f(ahnd1,iret)
426 if(iret.lt.0)then
427 call da_warning(__FILE__,__LINE__, &
428 (/"HDF5 read error for: Satellite Azimuth"/))
429 endif
430 call h5dclose_f(dhnd1,iret)
431 satazi=satazi0*slp+itp
432 ! cutoff overlap & change scale
434 do i=1,lo_rez_fov
435 do j=1,nscan
436 if(satazi(i,j).lt.0.0) satazi(i,j)=satazi(i,j)+360.0
437 enddo
438 enddo
439 do j=nscan+1,max_scan
440 satazi(:,j)=0
441 enddo
443 ! sample display
444 write(unit=stdout,fmt=*)'sat_azi',satazi(tesp,tess)
446 ! read array: solar zenith
447 ! read
448 sz1(1)=lo_rez_fov
449 sz1(2)=max_scan
450 call h5dopen_f(fhnd1,"/Geolocation/SolarZenith",dhnd1,iret)
451 call h5dread_f(dhnd1,H5T_NATIVE_INTEGER,solzen0,sz1,iret)
452 call h5aopen_f(dhnd1,"Slope",ahnd1,iret)
453 call h5aread_f(ahnd1,H5T_NATIVE_REAL,slp,sz1,iret)
454 call h5aclose_f(ahnd1,iret)
455 call h5aopen_f(dhnd1,"Intercept",ahnd1,iret)
456 call h5aread_f(ahnd1,H5T_NATIVE_REAL,itp,sz1,iret)
457 call h5aclose_f(ahnd1,iret)
458 if(iret.lt.0)then
459 call da_warning(__FILE__,__LINE__, &
460 (/"HDF5 read error for: Solar Zenith"/))
461 endif
462 call h5dclose_f(dhnd1,iret)
463 solzen=solzen0*slp+itp
464 ! cutoff overlap & change scale
465 do j=nscan+1,max_scan
466 solzen(:,j)=0
467 enddo
468 ! sample display
469 write(unit=stdout,fmt=*)'sol_zen',solzen(tesp,tess)
471 ! read array: solar azimuth
472 ! read
473 sz1(1)=lo_rez_fov
474 sz1(2)=max_scan
475 call h5dopen_f(fhnd1,"/Geolocation/SolarAzimuth",dhnd1,iret)
476 call h5dread_f(dhnd1,H5T_NATIVE_INTEGER,solazi0,sz1,iret)
477 call h5aopen_f(dhnd1,"Slope",ahnd1,iret)
478 call h5aread_f(ahnd1,H5T_NATIVE_REAL,slp,sz1,iret)
479 call h5aclose_f(ahnd1,iret)
480 call h5aopen_f(dhnd1,"Intercept",ahnd1,iret)
481 call h5aread_f(ahnd1,H5T_NATIVE_REAL,itp,sz1,iret)
482 call h5aclose_f(ahnd1,iret)
483 if(iret.lt.0)then
484 call da_warning(__FILE__,__LINE__, &
485 (/"HDF5 read error for: Solar Azimuth"/))
486 endif
487 call h5dclose_f(dhnd1,iret)
488 solazi=solazi0*slp+itp
490 ! cutoff overlap & change scale
491 do i=1,lo_rez_fov
492 do j=1,nscan
493 if(solazi(i,j).lt.0.0) solazi(i,j)=solazi(i,j)+360.0
494 enddo
495 enddo
496 do j=nscan+1,max_scan
497 solazi(:,j)=0
498 enddo
500 ! sample display
501 write(unit=stdout,fmt=*)'sol_azi',solazi(tesp,tess)
503 ! close file and HDF5
504 call h5fclose_f(fhnd1,iret)
506 ! 2.0 Loop to read hdf file and assign information to a sequential structure
507 !-------------------------------------------------------------------------
509 ! Allocate arrays to hold data
510 if ( .not. head_allocated ) then
511 allocate (head)
512 nullify ( head % next )
513 p => head
514 head_allocated = .true.
515 end if
516 ! start scan_loop
517 scan_loop: do iscan=1, nscan
518 do i = 1, 6
519 idate5(i)=date(iscan,i)
520 end do
521 call da_get_julian_time(idate5(1),idate5(2),idate5(3),idate5(4),idate5(5),obs_time)
522 if ( obs_time < time_slots(0) .or. &
523 obs_time >= time_slots(num_fgat_time) ) cycle scan_loop
524 do ifgat=1,num_fgat_time
525 if ( obs_time >= time_slots(ifgat-1) .and. &
526 obs_time < time_slots(ifgat) ) exit
527 end do
529 ! start fov_loop
530 fov_loop: do ifov=1, lo_rez_fov
531 num_mwhs2_file = num_mwhs2_file + 1
532 num_mwhs2_file_local = num_mwhs2_file_local + 1
533 info%lat = latlr(ifov,iscan)
534 info%lon = lonlr(ifov,iscan)
536 call da_llxy (info, loc, outside, outside_all)
537 if (outside_all) cycle fov_loop
539 num_mwhs2_global = num_mwhs2_global + 1
540 num_mwhs2_global_local = num_mwhs2_global_local + 1
541 ptotal(ifgat) = ptotal(ifgat) + 1
542 if (outside) cycle fov_loop ! No good for this PE
543 ! Discard data over Land (landmask =0 -->Land =1 -->Sea)
544 landsea_mask = 0
545 if(loflo(ifov,iscan).eq.3) landsea_mask = 1
546 !!! if( landsea_mask == 0 ) cycle fov_loop
548 num_mwhs2_local = num_mwhs2_local + 1
549 num_mwhs2_local_local = num_mwhs2_local_local + 1
550 write(unit=info%date_char, &
551 fmt='(i4.4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)') &
552 idate5(1), '-', idate5(2), '-', idate5(3), '_', idate5(4), &
553 ':', idate5(5), ':', idate5(6)
554 info%elv = elev(ifov,iscan)
556 ! 3.0 Make Thinning
557 ! Map obs to thinning grid
558 !-------------------------------------------------------------------
559 if (thinning) then
560 dlat_earth = info%lat !degree
561 dlon_earth = info%lon
562 if (dlon_earth<zero) dlon_earth = dlon_earth+r360
563 if (dlon_earth>=r360) dlon_earth = dlon_earth-r360
564 dlat_earth = dlat_earth*deg2rad !radian
565 dlon_earth = dlon_earth*deg2rad
566 crit = 1.
567 call map2grids(inst,ifgat,dlat_earth,dlon_earth,crit,iobs,itx,1,itt,iout,iuse)
568 if (.not. iuse) then
569 num_mwhs2_thinned = num_mwhs2_thinned+1
570 cycle fov_loop
571 end if
572 end if
574 num_mwhs2_used = num_mwhs2_used + 1
577 ! 4.0 assign information to sequential radiance structure
578 !--------------------------------------------------------------------------
579 allocate ( p % tb_inv (1:nchan ))
580 p%info = info
581 p%loc = loc
582 p%landsea_mask = landsea_mask
583 p%scanpos = ifov
584 p%satzen = satzen(ifov,iscan)
585 p%satazi = satazi(ifov,iscan)
586 p%solzen = solzen(ifov,iscan)
587 p%solazi = solazi(ifov,iscan)
588 p%tb_inv(1:nchan) = tbb(ifov,iscan,1:nchan)
589 p%sensor_index = inst
590 p%ifgat = ifgat
592 allocate (p%next) ! add next data
593 p => p%next
594 nullify (p%next)
595 end do fov_loop
596 end do scan_loop
598 ! Dellocate arrays
600 print*,'nscan=',nscan
601 print*,'nfov=',lo_rez_fov
602 write(stdout,fmt='(3a,i7)') ' In file: ',trim(fname_tb(ifile)),' got num_mwhs2_file : ',num_mwhs2_file_local
603 write(stdout,fmt='(3a,i7)') ' In file: ',trim(fname_tb(ifile)),' got num_mwhs2_global : ',num_mwhs2_global_local
604 write(stdout,fmt='(3a,i7)') ' In file: ',trim(fname_tb(ifile)),' got num_mwhs2_local : ',num_mwhs2_local_local
605 end do infile_loop
607 call H5close_f(iret)
610 if (thinning .and. num_mwhs2_global > 0 ) then
611 #ifdef DM_PARALLEL
612 ! Get minimum crit and associated processor index.
613 j = 0
614 do ifgat = 1, num_fgat_time
615 j = j + thinning_grid(inst,ifgat)%itxmax
616 end do
618 allocate ( in (j) )
619 allocate ( out (j) )
620 j = 0
621 do ifgat = 1, num_fgat_time
622 do i = 1, thinning_grid(inst,ifgat)%itxmax
623 j = j + 1
624 in(j) = thinning_grid(inst,ifgat)%score_crit(i)
625 end do
626 end do
627 call mpi_reduce(in, out, j, true_mpi_real, mpi_min, root, comm, ierr)
629 call wrf_dm_bcast_real (out, j)
631 j = 0
632 do ifgat = 1, num_fgat_time
633 do i = 1, thinning_grid(inst,ifgat)%itxmax
634 j = j + 1
635 if ( ABS(out(j)-thinning_grid(inst,ifgat)%score_crit(i)) > 1.0E-10 ) &
636 thinning_grid(inst,ifgat)%ibest_obs(i) = 0
637 end do
638 end do
640 deallocate( in )
641 deallocate( out )
643 #endif
645 ! Delete the nodes which being thinning out
646 p => head
647 prev => head
648 head_found = .false.
649 num_mwhs2_used_tmp = num_mwhs2_used
650 do j = 1, num_mwhs2_used_tmp
651 n = p%sensor_index
652 ifgat = p%ifgat
653 found = .false.
655 do i = 1, thinning_grid(n,ifgat)%itxmax
656 if ( thinning_grid(n,ifgat)%ibest_obs(i) == j .and. thinning_grid(n,ifgat)%score_crit(i) < 9.99e6_r_kind ) then
657 found = .true.
658 exit
659 end if
660 end do
662 ! free current data
663 if ( .not. found ) then
664 current => p
665 p => p%next
666 if ( head_found ) then
667 prev%next => p
668 else
669 head => p
670 prev => p
671 end if
672 deallocate ( current % tb_inv )
673 deallocate ( current )
674 num_mwhs2_thinned = num_mwhs2_thinned + 1
675 num_mwhs2_used = num_mwhs2_used - 1
676 continue
677 end if
679 if ( found .and. head_found ) then
680 prev => p
681 p => p%next
682 continue
683 end if
685 if ( found .and. .not. head_found ) then
686 head_found = .true.
687 head => p
688 prev => p
689 p => p%next
690 end if
692 end do
694 end if ! End of thinning
696 iv%total_rad_pixel = iv%total_rad_pixel + num_mwhs2_used
697 iv%total_rad_channel = iv%total_rad_channel + num_mwhs2_used*nchan
699 iv%info(radiance)%nlocal = iv%info(radiance)%nlocal + num_mwhs2_used
700 iv%info(radiance)%ntotal = iv%info(radiance)%ntotal + num_mwhs2_global
702 do i = 1, num_fgat_time
703 ptotal(i) = ptotal(i) + ptotal(i-1)
704 iv%info(radiance)%ptotal(i) = iv%info(radiance)%ptotal(i) + ptotal(i)
705 end do
706 if ( iv%info(radiance)%ptotal(num_fgat_time) /= iv%info(radiance)%ntotal ) then
707 write(unit=message(1),fmt='(A,I10,A,I10)') &
708 "Number of ntotal:",iv%info(radiance)%ntotal," is different from the sum of ptotal:", iv%info(radiance)%ptotal(num_fgat_time)
709 call da_warning(__FILE__,__LINE__,message(1:1))
710 endif
712 write(unit=stdout,fmt='(a)') 'MWHS2 data counts: '
713 write(stdout,fmt='(a,i7)') ' In file: ',num_mwhs2_file
714 write(stdout,fmt='(a,i7)') ' Global : ',num_mwhs2_global
715 write(stdout,fmt='(a,i7)') ' Local : ',num_mwhs2_local
716 write(stdout,fmt='(a,i7)') ' Used : ',num_mwhs2_used
717 write(stdout,fmt='(a,i7)') ' Thinned: ',num_mwhs2_thinned
719 ! 5.0 allocate innovation radiance structure
720 !----------------------------------------------------------------
722 if (num_mwhs2_used > 0) then
723 iv%instid(inst)%num_rad = num_mwhs2_used
724 iv%instid(inst)%info%nlocal = num_mwhs2_used
725 write(UNIT=stdout,FMT='(a,i3,2x,a,3x,i10)') &
726 'Allocating space for radiance innov structure', &
727 inst, iv%instid(inst)%rttovid_string, iv%instid(inst)%num_rad
728 call da_allocate_rad_iv (inst, nchan, iv)
729 end if
731 ! 6.0 assign sequential structure to innovation structure
732 !-------------------------------------------------------------
733 p => head
735 do n = 1, num_mwhs2_used
736 i = p%sensor_index
737 call da_initialize_rad_iv (i, n, iv, p)
738 current => p
739 p => p%next
740 ! free current data
741 deallocate ( current % tb_inv )
742 deallocate ( current )
743 end do
744 deallocate ( p )
745 deallocate (ptotal)
747 if (trace_use) call da_trace_exit("da_read_obs_hdf5mwhs2")
748 #else
749 call da_error(__FILE__,__LINE__,(/"Needs to be compiled with HDF5 library"/))
750 #endif
751 end subroutine da_read_obs_hdf5mwhs2
753 subroutine mjd2cal(ju,yr,mt,dy,hr,mn,sc)
754 implicit none
755 real(kind=8)::ju,j0
756 integer::yr,mt,dy,mn,yr0,sc
757 real(kind=4)::bc,dd,n1,n2,n3,sc0
758 integer::hr
760 ju=ju+2451545.0
762 if (ju.lt.1721423.5) then
763 bc=1
764 else
765 bc=0
766 end if
768 if (ju.lt.2299160.5) then
769 j0=floor(ju+0.5)
770 dd=ju+0.5-j0
771 else
772 n1=floor((ju-2342031.5)/36524.25/4)+1
773 n2=floor((ju-2378555.5)/36524.25/4)+1
774 n3=floor((ju-2415079.5)/36524.25/4)+1
775 j0=n1+n2+n3+ju+10
776 dd=j0+0.5-floor(j0+0.5)
777 j0=floor(j0+0.5)
778 end if
780 j0=j0+32083
781 yr0=ceiling(j0/365.25)-1
782 yr=yr0-4800
783 dy=j0-floor(yr0*365.25)
784 mt=floor((dy-0.6)/30.6)+3
785 dy=dy-nint((mt-3)*30.6)
787 if (mt.gt.12) then
788 mt=mt-12
789 yr=yr+1
790 end if
792 yr=yr-bc
793 sc0=nint(dd*86400)
796 hr=floor(sc0/3600)
797 sc0=sc0-hr*3600
798 mn=floor(sc0/60)
799 sc=int(sc0-mn*60)
801 end subroutine mjd2cal