IODA Bundle
radiance_mod.f90
Go to the documentation of this file.
1 module radiance_mod
2 
3 use kinds, only: r_kind,i_kind,r_double
4 use define_mod, only: missing_r, missing_i, nstring, ndatetime, &
5  ninst, inst_list, set_name_satellite, set_name_sensor, xdata, name_sen_info, &
6  nvar_info, name_var_info, type_var_info, nsen_info, type_sen_info, set_brit_obserr
7 use ufo_vars_mod, only: ufo_vars_getindex
8 use netcdf, only: nf90_float, nf90_int, nf90_char
9 
10 implicit none
11 private
12 public :: read_amsua_amsub_mhs
13 public :: read_airs_colocate_amsua
14 public :: sort_obs_radiance
15 
16 real(r_kind), parameter :: r8bfms = 9.0e08 ! threshold to check for BUFR missing value
17 
18 type datalink_radiance
19  character(len=nstring) :: msg_type ! BUFR message type name
20  integer(i_kind) :: satid ! satellite identifier
21  integer(i_kind) :: instid ! instrument identifier
22  character(len=nstring) :: inst ! instrument name eg. amsua-n15
23  character(len=ndatetime) :: datetime ! ccyy-mm-ddThh:mm:ssZ
24  integer(i_kind) :: nchan ! number of channels
25  real(r_kind) :: lat ! latitude in degree
26  real(r_kind) :: lon ! longitude in degree
27  real(r_kind) :: elv ! elevation in m
28  real(r_kind) :: dhr ! obs time minus analysis time in hour
29  integer(i_kind) :: inst_idx ! index of inst in inst_list
30  integer(i_kind) :: landsea
31  integer(i_kind) :: scanpos
32  integer(i_kind) :: scanline
33  real(r_kind) :: satzen
34  real(r_kind) :: satazi
35  real(r_kind) :: solzen
36  real(r_kind) :: solazi
37  real(r_kind), allocatable :: tb(:)
38  real(r_kind), allocatable :: ch(:)
39  type (datalink_radiance), pointer :: next ! pointer to next data
40 end type datalink_radiance
41 
42 type(datalink_radiance), pointer :: rhead=>null(), rlink=>null()
43 
44 contains
45 
46 !--------------------------------------------------------------
47 
48 subroutine read_amsua_amsub_mhs (filename, filedate)
49 
50 !| NC021023 | A61223 | MTYP 021-023 PROC AMSU-A 1B Tb (NOAA-15-19, METOP-1,2) |
51 !| NC021024 | A61224 | MTYP 021-024 PROCESSED AMSU-B 1B Tb (NOAA-15-17) |
52 !| NC021025 | A61225 | MTYP 021-025 PROCESSED HIRS-3 1B Tb (NOAA-15-17) |
53 !| NC021027 | A61234 | MTYP 021-027 PROCESSED MHS Tb (NOAA-18-19, METOP-1,2) |
54 !| NC021028 | A61245 | MTYP 021-028 PROC HIRS-4 1B Tb (NOAA-18-19, METOP-1,2) |
55 !| | |
56 !| NC021023 | YEAR MNTH DAYS HOUR MINU SECO 207002 CLAT CLON 207000 |
57 !| NC021023 | SAID SIID FOVN LSQL SAZA SOZA HOLS 202127 HMSL 202000 |
58 !| NC021023 | SOLAZI BEARAZ "BRITCSTC"15 |
59 !| | |
60 !| NC021024 | YEAR MNTH DAYS HOUR MINU SECO 207002 CLAT CLON 207000 |
61 !| NC021024 | SAID SIID FOVN LSQL SAZA SOZA HOLS 202127 HMSL 202000 |
62 !| NC021024 | SOLAZI BEARAZ "BRITCSTC"5 |
63 !| | |
64 !| NC021025 | YEAR MNTH DAYS HOUR MINU SECO 207002 CLAT CLON 207000 |
65 !| NC021025 | SAID SIID FOVN LSQL SAZA SOZA HOLS 202127 HMSL 202000 |
66 !| NC021025 | SOLAZI BEARAZ "BRIT"20 |
67 !| | |
68 !| NC021027 | YEAR MNTH DAYS HOUR MINU SECO 207002 CLAT CLON 207000 |
69 !| NC021027 | SAID SIID FOVN LSQL SAZA SOZA HOLS 202127 HMSL 202000 |
70 !| NC021027 | SOLAZI BEARAZ "BRITCSTC"5 |
71 !| | |
72 !| NC021028 | YEAR MNTH DAYS HOUR MINU SECO 207002 CLAT CLON 207000 |
73 !| NC021028 | SAID SIID FOVN LSQL SAZA SOZA HOLS 202127 HMSL 202000 |
74 !| NC021028 | SOLAZI BEARAZ "BRIT"20 |
75 !| | |
76 !| BRITCSTC | CHNM TMBR CSTC |
77 !| BRIT | CHNM TMBR |
78 
79  implicit none
80 
81  character (len=*), intent(in) :: filename
82  character (len=10), intent(out) :: filedate ! ccyymmddhh
83 
84 
85  integer(i_kind), parameter :: ntime = 6 ! number of data to read in timestr
86  integer(i_kind), parameter :: ninfo = 10 ! number of data to read in infostr
87  integer(i_kind), parameter :: nlalo = 2 ! number of data to read in lalostr
88  integer(i_kind), parameter :: nbrit = 2 ! number of data to read in britstr
89  integer(i_kind), parameter :: maxchan = 15 ! max nchan among amsua, amsub, mhs
90 
91  character(len=80) :: timestr, infostr, lalostr, britstr
92 
93  real(r_double), dimension(ntime) :: timedat
94  real(r_double), dimension(ninfo) :: infodat
95  real(r_double), dimension(nlalo) :: lalodat
96  real(r_double), dimension(2,maxchan) :: data1b8
97 
98  character(len=8) :: subset
99  character(len=10) :: cdate
100 
101  integer(i_kind) :: iunit, iost, iret, i
102  integer(i_kind) :: nchan
103  integer(i_kind) :: idate
104  integer(i_kind) :: num_report_infile
105  integer(i_kind) :: ireadmg, ireadsb
106 
107  integer(i_kind) :: iyear, imonth, iday, ihour, imin, isec
108  real(r_double) :: ref_time, obs_time
109 
110  write(*,*) '--- reading '//trim(filename)//' ---'
111 
112  timestr = 'YEAR MNTH DAYS HOUR MINU SECO'
113  infostr = 'SAID SIID FOVN LSQL SAZA SOZA HOLS HMSL SOLAZI BEARAZ'
114  lalostr = 'CLAT CLON'
115  britstr = 'CHNM TMBR'
116 
117  num_report_infile = 0
118 
119  iunit = 96
120 
121  ! open bufr file
122  open (unit=iunit, file=trim(filename), &
123  iostat=iost, form='unformatted', status='old')
124  if (iost /= 0) then
125  write(unit=*,fmt='(a,i5,a)') &
126  "Error",iost," opening BUFR obs file "//trim(filename)
127  return
128  end if
129 
130  call openbf(iunit,'IN',iunit)
131  call datelen(10)
132  call readmg(iunit,subset,idate,iret)
133 
134  if ( iret /= 0 ) then
135  write(unit=*,fmt='(A,I5,A)') &
136  "Error",iret," reading BUFR obs file "//trim(filename)
137  call closbf(iunit)
138  return
139  end if
140  rewind(iunit)
141 
142  write(unit=*,fmt='(1x,a,i10)') trim(filename)//' file date is: ', idate
143  write(unit=filedate, fmt='(i10)') idate
144  read (filedate(1:10),'(i4,3i2)') iyear, imonth, iday, ihour
145  call get_julian_time (iyear,imonth,iday,ihour,0,ref_time)
146 
147  if ( .not. associated(rhead) ) then
148  nullify ( rhead )
149  allocate ( rhead )
150  nullify ( rhead%next )
151  end if
152 
153  if ( .not. associated(rlink) ) then
154  rlink => rhead
155  else
156  allocate ( rlink%next )
157  rlink => rlink%next
158  nullify ( rlink%next )
159  end if
160 
161  msg_loop: do while (ireadmg(iunit,subset,idate)==0)
162 !print*,subset
163  subset_loop: do while (ireadsb(iunit)==0)
164 
165  num_report_infile = num_report_infile + 1
166 
167  call ufbint(iunit,timedat,ntime,1,iret,timestr)
168 
169  iyear = nint(timedat(1))
170  imonth = nint(timedat(2))
171  iday = nint(timedat(3))
172  ihour = nint(timedat(4))
173  imin = nint(timedat(5))
174  isec = min(59, nint(timedat(6))) ! raw BUFR data that has SECO = 60.0 SECOND
175  ! that was probably rounded from 59.x seconds
176  ! reset isec to 59 rather than advancing one minute
177  if ( iyear > 1900 .and. iyear < 3000 .and. &
178  imonth >= 1 .and. imonth <= 12 .and. &
179  iday >= 1 .and. iday <= 31 .and. &
180  ihour >= 0 .and. ihour < 24 .and. &
181  imin >= 0 .and. imin < 60 .and. &
182  isec >= 0 .and. isec < 60 ) then
183  write(unit=rlink%datetime, fmt='(i4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a)') &
184  iyear, '-', imonth, '-', iday, 'T', ihour, ':', imin, ':', isec, 'Z'
185  call get_julian_time (iyear,imonth,iday,ihour,imin,obs_time)
186  rlink%dhr = (obs_time + (isec/60.0) - ref_time)/60.0
187  else
188  cycle subset_loop
189  end if
190 
191  call ufbint(iunit,lalodat,nlalo,1,iret,lalostr)
192  if ( abs(lalodat(1)) > 90.0 .or. abs(lalodat(1)) > 360.0 ) cycle subset_loop
193 
194  call ufbint(iunit,infodat,ninfo,1,iret,infostr)
195  call ufbrep(iunit,data1b8,2,maxchan,nchan,britstr)
196 
197  rlink % nchan = nchan
198  if ( nchan > 0 ) then
199  allocate ( rlink % tb(nchan) ) ! brightness temperature
200  allocate ( rlink % ch(nchan) ) ! channel number
201  end if
202 
203  call fill_datalink(rlink, missing_r, missing_i)
204 
205  if ( lalodat(1) < r8bfms ) rlink % lat = lalodat(1)
206  if ( lalodat(2) < r8bfms ) rlink % lon = lalodat(2)
207 
208  rlink % satid = nint(infodat(1)) ! SAID satellite identifier
209  rlink % instid = nint(infodat(2)) ! SIID instrument identifier
210 
211  if ( infodat(3) < r8bfms ) rlink % scanpos = nint(infodat(3)) ! FOVN field of view number
212  if ( infodat(4) < r8bfms ) rlink % landsea = infodat(4) ! LSQL land sea qualifier 0:land, 1:sea, 2:coast
213  if ( infodat(5) < r8bfms ) rlink % satzen = infodat(5) ! SAZA satellite zenith angle (degree)
214  if ( infodat(10) < r8bfms ) rlink % satazi = infodat(10) ! BEARAZ satellite azimuth (degree true)
215  if ( infodat(6) < r8bfms ) rlink % solzen = infodat(6) ! SOZA solar zenith angle (degree)
216  if ( infodat(9) < r8bfms ) rlink % solazi = infodat(9) ! SOLAZI solar azimuth (degree true)
217  if ( infodat(7) < r8bfms ) rlink % elv = infodat(7) ! HOLS height of land surface (m)
218 
219  if ( nchan > 0 ) then
220  do i = 1, nchan
221  if ( data1b8(1,i) < r8bfms ) rlink % ch(i) = nint(data1b8(1,i))
222  if ( data1b8(2,i) < r8bfms ) rlink % tb(i) = data1b8(2,i)
223  end do
224  end if
225 
226  allocate ( rlink%next )
227  rlink => rlink%next
228  nullify ( rlink%next )
229 
230  end do subset_loop ! ireadsb
231  end do msg_loop ! ireadmg
232 
233  call closbf(iunit)
234  close(iunit)
235 
236  write(*,'(1x,a,a,a,i10)') 'num_report_infile ', trim(filename), ' : ', num_report_infile
237 
238 end subroutine read_amsua_amsub_mhs
239 
240 !--------------------------------------------------------------
241 
242 subroutine read_airs_colocate_amsua (filename, filedate)
243 
244 ! adapted from WRFDA/var/da/da_radiance/da_read_obs_bufrairs.inc
245 
246 !| NC021249 | A50243 | MTYP 021-249 EVERY FOV AIRS/AMSU-A/HSB 1B BTEMPS(AQUA) |
247 !| | |
248 !| NC021249 | SCO1C3IN |
249 !| | |
250 !| SCO1C3IN | SPITSEQN SITPSEQN (SCBTSEQN) "SVCASEQN"4 TOCC AMSUSPOT |
251 !| SCO1C3IN | "AMSUCHAN"15 HSBSPOT "HSBCHAN"5 |
252 !| | |
253 !| SPITSEQN | SAID ORBN 201133 SLNM 201000 201132 MJFC 201000 202126 |
254 !| SPITSEQN | SELV 202000 SOZA SOLAZI "INTMS"9 |
255 !| | |
256 !| SITPSEQN | SIID YEAR MNTH DAYS HOUR MINU 202131 201138 SECO 201000 |
257 !| SITPSEQN | 202000 CLATH CLONH SAZA BEARAZ FOVN |
258 !| | |
259 !| SCBTSEQN | 201134 CHNM 201000 LOGRCW ACQF TMBR |
260 !| | |
261 !| AMSUSPOT | SIID YEAR MNTH DAYS HOUR MINU 202131 201138 SECO 201000 |
262 !| AMSUSPOT | 202000 CLATH CLONH SAZA BEARAZ FOVN |
263 !| | |
264 !| AMSUCHAN | 201134 CHNM 201000 LOGRCW ACQF TMBR |
265 !| | |
266 
267  implicit none
268 
269  character (len=*), intent(in) :: filename
270  character (len=10), intent(out) :: filedate ! ccyymmddhh
271 
272  integer(i_kind),parameter :: n_maxchan = 281 ! max nchan of airs-281-subset and amsua-a
273 
274  ! variables for BUFR SPITSEQN
275  integer(i_kind),parameter :: n_satellitespot_list = 25
276  type satellitespot_list
277  sequence
278  real(r_double) :: said ! Satellite identifier
279  real(r_double) :: orbn ! Orbit number
280  real(r_double) :: slnm ! Scan line number
281  real(r_double) :: mjfc ! Major frame count
282  real(r_double) :: selv ! Height of station
283  real(r_double) :: soza ! Solar zenith angle
284  real(r_double) :: solazi ! Solar azimuth angle
285  real(r_double) :: intms(2,9) ! SATELLITE inSTRUMENT TEMPERATURES
286  end type satellitespot_list
287  real(r_double), dimension(1:N_satellitespot_LIST) :: satellitespot_list_array
288 
289  ! variables for BUFR SITPSEQN/AMSUSPOT
290  integer(i_kind),parameter :: n_sensorspot_list = 12
291  type sensorspot_list
292  sequence
293  real(r_double) :: siid ! Satellite instruments
294  real(r_double) :: year
295  real(r_double) :: mnth
296  real(r_double) :: days
297  real(r_double) :: hour
298  real(r_double) :: minu
299  real(r_double) :: seco
300  real(r_double) :: clath ! Latitude (high accuracy)
301  real(r_double) :: clonh ! Longitude (high accuracy)
302  real(r_double) :: saza ! Satellite zenith angle
303  real(r_double) :: bearaz ! Bearing or azimuth
304  real(r_double) :: fovn ! Field of view number
305  end type sensorspot_list
306  real(r_double), dimension(1:N_sensorspot_LIST) :: sensorspot_list_array
307 
308  ! variables for BUFR SCBTSEQN/AMSUCHAN
309  integer(i_kind),parameter :: n_sensorchan_list = 4
310  type sensorchan_list
311  sequence
312  real(r_double) :: chnm ! Channel number
313  real(r_double) :: logrcw ! Log-10 of temperature-radiance central wavenumber
314  real(r_double) :: acqf ! Channel quality flags for ATOVS
315  real(r_double) :: tmbr ! Brightness temperature
316  end type sensorchan_list
317  real(r_double), dimension(1:N_sensorchan_LIST,1:N_MAXCHAN) :: sensorchan_list_array
318 
319  ! Variables for BUFR data
320  type(satellitespot_list) :: satellitespot
321  type(sensorspot_list) :: sensorspot
322  type(sensorchan_list) :: sensorchan(n_maxchan)
323 
324  character(len=8) :: subset
325  character(len=8) :: spotname
326  character(len=8) :: channame
327  integer(i_kind) :: ireadmg, ireadsb
328  integer(i_kind) :: nchan
329  integer(i_kind) :: iret
330 
331  ! Work variables for time
332  integer(i_kind) :: idate
333  integer(i_kind) :: iyear, imonth, iday, ihour, imin, isec
334  real(r_double) :: ref_time, obs_time
335 
336  ! Other work variables
337  integer(i_kind) :: i, ich
338  integer(i_kind) :: iost, iunit
339  integer(i_kind) :: num_report_infile
340  logical :: decode_airs, decode_amsua
341 
342 
343  write(*,*) '--- reading '//trim(filename)//' ---'
344 
345  decode_amsua = .false.
346  decode_airs = .false.
347  if ( ufo_vars_getindex(inst_list, 'amsua_aqua') > 0 ) decode_amsua = .true.
348  if ( ufo_vars_getindex(inst_list, 'airs_aqua') > 0 ) decode_airs = .true.
349 
350  num_report_infile = 0
351 
352  iunit = 97
353 
354  ! open bufr file
355  open (unit=iunit, file=trim(filename), &
356  iostat=iost, form='unformatted', status='old')
357  if (iost /= 0) then
358  write(unit=*,fmt='(a,i5,a)') &
359  "Error",iost," opening BUFR obs file "//trim(filename)
360  return
361  end if
362 
363  call openbf(iunit,'IN',iunit)
364  call datelen(10)
365  call readmg(iunit,subset,idate,iret)
366 
367  if ( iret /= 0 ) then
368  write(unit=*,fmt='(A,I5,A)') &
369  "Error",iret," reading BUFR obs file "//trim(filename)
370  call closbf(iunit)
371  return
372  end if
373  rewind(iunit)
374 
375  write(unit=*,fmt='(1x,a,i10)') trim(filename)//' file date is: ', idate
376  write(unit=filedate, fmt='(i10)') idate
377  read (filedate(1:10),'(i4,3i2)') iyear, imonth, iday, ihour
378  call get_julian_time (iyear,imonth,iday,ihour,0,ref_time)
379 
380  if ( .not. associated(rhead) ) then
381  nullify ( rhead )
382  allocate ( rhead )
383  nullify ( rhead%next )
384  end if
385 
386  if (.not. associated(rlink)) then
387  rlink => rhead
388  else
389  allocate ( rlink%next )
390  rlink => rlink%next
391  nullify ( rlink%next )
392  end if
393 
394  msg_loop: do while ( ireadmg(iunit,subset,idate)==0 )
395 
396  subset_loop: do while ( ireadsb(iunit)==0 )
397 
398  num_report_infile = num_report_infile + 1
399 
400  ! Read SPITSEQN
401  call ufbseq(iunit,satellitespot_list_array,n_satellitespot_list,1,iret,'SPITSEQN')
402  satellitespot = satellitespot_list( &
403  satellitespot_list_array(1), &
404  satellitespot_list_array(2), &
405  satellitespot_list_array(3), &
406  satellitespot_list_array(4), &
407  satellitespot_list_array(5), &
408  satellitespot_list_array(6), &
409  satellitespot_list_array(7), &
410  reshape(satellitespot_list_array(8:25), (/2,9/)) )
411 
412  loop_sensor: do i = 1, 2
413  if ( i == 1 ) then
414  if ( decode_airs ) then
415  spotname = 'SITPSEQN'
416  channame = 'SCBTSEQN'
417  else
418  cycle loop_sensor
419  end if
420  else if ( i == 2 ) then
421  if ( decode_amsua ) then
422  spotname = 'AMSUSPOT'
423  channame = 'AMSUCHAN'
424  else
425  exit loop_sensor
426  end if
427  end if
428 
429  ! Read SITPSEQN / AMSUSPOT
430  call ufbseq(iunit,sensorspot_list_array,n_sensorspot_list,1,iret,spotname)
431 
432  sensorspot = sensorspot_list( sensorspot_list_array(1), &
433  sensorspot_list_array(2), &
434  sensorspot_list_array(3), &
435  sensorspot_list_array(4), &
436  sensorspot_list_array(5), &
437  sensorspot_list_array(6), &
438  sensorspot_list_array(7), &
439  sensorspot_list_array(8), &
440  sensorspot_list_array(9), &
441  sensorspot_list_array(10), &
442  sensorspot_list_array(11), &
443  sensorspot_list_array(12) )
444 
445  iyear = nint(sensorspot%year)
446  imonth = nint(sensorspot%mnth)
447  iday = nint(sensorspot%days)
448  ihour = nint(sensorspot%hour)
449  imin = nint(sensorspot%minu)
450  isec = nint(sensorspot%seco)
451  if ( iyear > 1900 .and. iyear < 3000 .and. &
452  imonth >= 1 .and. imonth <= 12 .and. &
453  iday >= 1 .and. iday <= 31 .and. &
454  ihour >= 0 .and. ihour < 24 .and. &
455  imin >= 0 .and. imin < 60 .and. &
456  isec >= 0 .and. isec < 60 ) then
457  write(unit=rlink%datetime, fmt='(i4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a)') &
458  iyear, '-', imonth, '-', iday, 'T', ihour, ':', imin, ':', isec, 'Z'
459  call get_julian_time (iyear,imonth,iday,ihour,imin,obs_time)
460  rlink%dhr = (obs_time + (isec/60.0) - ref_time)/60.0
461  else
462  cycle subset_loop
463  end if
464 
465  if ( abs(sensorspot%clath) > 90.0 .or. abs(sensorspot%clonh) > 360.0 ) cycle subset_loop
466 
467  ! Read SCBTSEQN or AMSUCHAN
468  call ufbseq(iunit,sensorchan_list_array,n_sensorchan_list,n_maxchan,nchan,channame)
469 
470  rlink % nchan = nchan
471  if ( nchan > 0 ) then
472  allocate ( rlink % tb(1:nchan) )
473  allocate ( rlink % ch(1:nchan) )
474  end if
475 
476  call fill_datalink(rlink, missing_r, missing_i)
477 
478  do ich = 1 , nchan
479  sensorchan(ich) = sensorchan_list( sensorchan_list_array(1,ich), &
480  sensorchan_list_array(2,ich), &
481  sensorchan_list_array(3,ich), &
482  sensorchan_list_array(4,ich) )
483  if ( sensorchan(ich) % tmbr < r8bfms ) rlink % tb(ich) = sensorchan(ich) % tmbr
484  if ( sensorchan(ich) % chnm < r8bfms ) rlink % ch(ich) = sensorchan(ich) % chnm
485  end do
486 
487  if ( sensorspot%clath < r8bfms ) rlink % lat = sensorspot%clath
488  if ( sensorspot%clonh < r8bfms ) rlink % lon = sensorspot%clonh
489 
490  if ( sensorspot%fovn < r8bfms ) rlink % scanpos = nint( sensorspot%fovn )
491  if ( sensorspot%saza < r8bfms ) rlink % satzen = sensorspot%saza
492  if ( satellitespot%slnm < r8bfms ) rlink % scanline = nint(satellitespot%slnm)
493  if ( sensorspot%bearaz < r8bfms ) rlink % satazi = sensorspot%bearaz
494  if ( satellitespot%soza < r8bfms ) rlink % solzen = satellitespot%soza
495  if ( satellitespot%solazi < r8bfms ) rlink % solazi = satellitespot%solazi
496 
497  rlink % satid = nint(satellitespot % said)
498  rlink % instid = nint(sensorspot % siid)
499 
500  allocate ( rlink%next )
501  rlink => rlink%next
502  nullify ( rlink%next )
503 
504  end do loop_sensor
505  end do subset_loop ! ireadsb
506  end do msg_loop ! ireadmg
507 
508  call closbf(iunit)
509  close(iunit)
510 
511  write(*,'(1x,a,a,a,i10)') 'num_report_infile ', trim(filename), ' : ', num_report_infile
512 
513 end subroutine read_airs_colocate_amsua
514 
515 !--------------------------------------------------------------
516 
517 subroutine sort_obs_radiance
518 
519  implicit none
520 
521  integer(i_kind) :: i, iv, k
522  integer(i_kind) :: ityp, irec
523  integer(i_kind), dimension(ninst) :: nrecs
524  integer(i_kind), dimension(ninst) :: nlocs
525  integer(i_kind), dimension(ninst) :: nvars
526  integer(i_kind), dimension(ninst) :: iloc
527  character(len=nstring) :: satellite
528  character(len=nstring) :: sensor
529 
530  nrecs(:) = 0
531  nlocs(:) = 0
532  nvars(:) = 0
533 
534  write(*,*) '--- sorting radiance obs...'
535 
536  ! set inst type from satellite id and sensor id
537  ! and count the numbers
538  rlink => rhead
539  set_inst_loop: do while ( associated(rlink) )
540 
541  if ( rlink%nchan < 1 ) then
542  rlink => rlink%next
543  cycle set_inst_loop
544  end if
545 
546  call set_name_satellite(rlink%satid, satellite)
547  call set_name_sensor (rlink%instid, sensor)
548  rlink % inst = trim(sensor)//'_'//trim(satellite)
549 !print*, rlink%inst, rlink % nchan
550 
551  ! find index of inst in inst_list
552  rlink%inst_idx = ufo_vars_getindex(inst_list, rlink%inst)
553  if ( rlink % inst_idx > 0 ) then
554  ! obtype assigned, advance ob counts
555  nrecs(rlink%inst_idx) = nrecs(rlink%inst_idx) + 1
556  nlocs(rlink%inst_idx) = nlocs(rlink%inst_idx) + 1
557  nvars(rlink%inst_idx) = rlink % nchan
558  end if
559 
560  rlink => rlink%next
561  end do set_inst_loop
562 
563  !write(*,*) 'num_report_decoded = ', sum(nrecs(:))
564  write(*,'(1x,20x,a10)') 'nlocs'
565  do i = 1, ninst
566  !write(*,'(1x,a20,2i10)') inst_list(i), nrecs(i), nlocs(i)
567  write(*,'(1x,a20,i10)') inst_list(i), nlocs(i)
568  end do
569 
570  ! allocate data arrays with the just counted numbers
571  allocate (xdata(ninst))
572  do i = 1, ninst
573  xdata(i) % nrecs = nrecs(i)
574  xdata(i) % nlocs = nlocs(i)
575  xdata(i) % nvars = nvars(i)
576 
577  if ( nlocs(i) > 0 ) then
578  allocate (xdata(i)%xinfo_float(nlocs(i), nvar_info))
579  allocate (xdata(i)%xinfo_int (nlocs(i), nvar_info))
580  allocate (xdata(i)%xinfo_char (nlocs(i), nvar_info))
581  allocate (xdata(i)%xseninfo_float(nlocs(i), nsen_info))
582  allocate (xdata(i)%xseninfo_int (nvars(i), nsen_info))
583  xdata(i)%xinfo_float (:,:) = missing_r
584  xdata(i)%xinfo_int (:,:) = missing_i
585  xdata(i)%xinfo_char (:,:) = ''
586  xdata(i)%xseninfo_float(:,:) = missing_r
587  xdata(i)%xseninfo_int (:,:) = missing_i
588  if ( nvars(i) > 0 ) then
589  allocate (xdata(i)%xfield(nlocs(i), nvars(i)))
590  xdata(i)%xfield(:,:)%val = missing_r
591  xdata(i)%xfield(:,:)%qm = missing_i
592  xdata(i)%xfield(:,:)%err = missing_r
593  allocate (xdata(i)%var_idx(nvars(i)))
594  do iv = 1, nvars(i)
595  xdata(i)%var_idx(iv) = iv
596  end do
597  end if
598  end if
599  end do
600 
601  ! transfer data from rlink to xdata
602 
603  iloc(:) = 0
604  irec = 0
605 
606  rlink => rhead
607  reports: do while ( associated(rlink) )
608  irec = irec + 1
609  ityp = rlink%inst_idx
610  if ( ityp < 0 ) then
611  rlink => rlink%next
612  cycle reports
613  end if
614  if ( rlink%nchan < 1 ) then
615  rlink => rlink%next
616  cycle reports
617  end if
618 
619  iloc(ityp) = iloc(ityp) + 1
620 
621  do i = 1, nvar_info
622  if ( type_var_info(i) == nf90_int ) then
623  if ( trim(name_var_info(i)) == 'record_number' ) then
624  xdata(ityp)%xinfo_int(iloc(ityp),i) = irec
625  end if
626  else if ( type_var_info(i) == nf90_float ) then
627  if ( name_var_info(i) == 'time' ) then
628  xdata(ityp)%xinfo_float(iloc(ityp),i) = rlink%dhr
629  else if ( trim(name_var_info(i)) == 'station_elevation' ) then
630  xdata(ityp)%xinfo_float(iloc(ityp),i) = rlink%elv
631  else if ( trim(name_var_info(i)) == 'latitude' ) then
632  xdata(ityp)%xinfo_float(iloc(ityp),i) = rlink%lat
633  else if ( trim(name_var_info(i)) == 'longitude' ) then
634  xdata(ityp)%xinfo_float(iloc(ityp),i) = rlink%lon
635  end if
636  else if ( type_var_info(i) == nf90_char ) then
637  if ( trim(name_var_info(i)) == 'datetime' ) then
638  xdata(ityp)%xinfo_char(iloc(ityp),i) = rlink%datetime
639  else if ( trim(name_var_info(i)) == 'station_id' ) then
640  xdata(ityp)%xinfo_char(iloc(ityp),i) = rlink%inst
641  end if
642  end if
643  end do
644 
645  do i = 1, nsen_info
646  if ( type_sen_info(i) == nf90_float ) then
647  if ( trim(name_sen_info(i)) == 'scan_position' ) then
648  xdata(ityp)%xseninfo_float(iloc(ityp),i) = rlink%scanpos
649  else if ( trim(name_sen_info(i)) == 'sensor_zenith_angle' ) then
650  xdata(ityp)%xseninfo_float(iloc(ityp),i) = rlink%satzen
651  else if ( trim(name_sen_info(i)) == 'sensor_azimuth_angle' ) then
652  xdata(ityp)%xseninfo_float(iloc(ityp),i) = rlink%satazi
653  else if ( trim(name_sen_info(i)) == 'solar_azimuth_angle' ) then
654  xdata(ityp)%xseninfo_float(iloc(ityp),i) = rlink%solzen
655  else if ( trim(name_sen_info(i)) == 'sensor_azimuth_angle' ) then
656  xdata(ityp)%xseninfo_float(iloc(ityp),i) = rlink%solazi
657  else if ( trim(name_sen_info(i)) == 'sensor_view_angle' ) then
658  call calc_sensor_view_angle(trim(rlink%inst), rlink%scanpos, xdata(ityp)%xseninfo_float(iloc(ityp),i))
659  end if
660 ! else if ( type_sen_info(i) == nf90_int ) then
661 ! else if ( type_sen_info(i) == nf90_char ) then
662  end if
663  end do
664 
665  iv = ufo_vars_getindex(name_sen_info, 'sensor_channel')
666  xdata(ityp)%xseninfo_int(:,iv) = rlink%ch(:)
667 
668  do i = 1, nvars(ityp)
669  xdata(ityp)%xfield(iloc(ityp),i)%val = rlink%tb(i)
670  !xdata(ityp)%xfield(iloc(ityp),i)%err = 1.0
671  call set_brit_obserr(trim(rlink%inst), i, xdata(ityp)%xfield(iloc(ityp),i)%err)
672  xdata(ityp)%xfield(iloc(ityp),i)%qm = 0
673  end do
674  rlink => rlink%next
675  end do reports
676 
677  ! done with rlink
678  ! release the linked list
679  rlink => rhead
680  do while ( associated(rlink) )
681  rhead => rlink%next
682  if ( allocated (rlink%tb) ) deallocate (rlink%tb)
683  if ( allocated (rlink%ch) ) deallocate (rlink%ch)
684  if ( associated (rlink) ) deallocate (rlink)
685  rlink => rhead
686  end do
687  nullify (rhead)
688 
689 end subroutine sort_obs_radiance
690 
691 !--------------------------------------------------------------
692 
693 subroutine fill_datalink (datalink, rfill, ifill)
694 
695  implicit none
696 
697  type (datalink_radiance), intent(inout) :: datalink
698  real(r_kind), intent(in) :: rfill ! fill value in real
699  integer(i_kind), intent(in) :: ifill ! fill value in integer
700 
701  integer(i_kind) :: i
702 
703  !datalink % datetime = ''
704  datalink % lat = rfill
705  datalink % lon = rfill
706  datalink % satid = ifill
707  datalink % instid = ifill
708  datalink % scanpos = rfill
709  datalink % landsea = rfill
710  datalink % satzen = rfill
711  datalink % satazi = rfill
712  datalink % solzen = rfill
713  datalink % solazi = rfill
714  datalink % elv = rfill
715  if ( allocated (datalink % tb) ) then
716  datalink % tb(:) = rfill
717  end if
718  if ( allocated (datalink % ch) ) then
719  datalink % ch(:) = ifill
720  end if
721 
722 end subroutine fill_datalink
723 
724 subroutine calc_sensor_view_angle(name_inst, ifov, view_angle)
725 
726 ! calculate sensor view angle from given scan position (field of view number)
727 
728  implicit none
729 
730  character(len=*), intent(in) :: name_inst ! instrument name eg. amsua_n15
731  integer(i_kind), intent(in) :: ifov ! field of view number
732  real(r_kind), intent(out) :: view_angle ! sensor view angle
733 
734  integer(i_kind) :: idx
735  real(r_kind) :: start, step
736 
737  view_angle = missing_r
738 
739  idx = index(name_inst, '_')
740  select case ( name_inst(1:idx-1) )
741  case ( 'amsua' )
742  start = -48.0_r_kind - 1.0_r_kind/3.0_r_kind
743  step = 3.0_r_kind + 1.0_r_kind/3.0_r_kind
744  case ( 'amsub' )
745  start = -48.95_r_kind
746  step = 1.1_r_kind
747  case ( 'mhs' )
748  start = -445.0_r_kind/9.0_r_kind
749  step = 10.0_r_kind/9.0_r_kind
750  case ( 'atms' )
751  start = -52.725_r_kind
752  step = 1.11_r_kind
753  case default
754  return
755  end select
756 
757  view_angle = start + float(ifov-1) * step
758 
759 end subroutine calc_sensor_view_angle
760 
761 subroutine get_julian_time(year,month,day,hour,minute,gstime)
762 
763 ! taken from WRFDA/var/da/da_tools/da_get_julian_time.inc
764 
765  implicit none
766 
767  integer(i_kind), intent(in) :: year
768  integer(i_kind), intent(in) :: month
769  integer(i_kind), intent(in) :: day
770  integer(i_kind), intent(in) :: hour
771  integer(i_kind), intent(in) :: minute
772  real(r_double), intent(out) :: gstime
773 
774  integer(i_kind) :: iw3jdn, ndays, nmind
775 
776  iw3jdn = day - 32075 &
777  + 1461 * (year + 4800 + (month - 14) / 12) / 4 &
778  + 367 * (month - 2 - (month - 14) / 12 * 12) / 12 &
779  - 3 * ((year + 4900 + (month - 14) / 12) / 100) / 4
780  ndays = iw3jdn - 2443510
781 
782  nmind = ndays*1440 + hour * 60 + minute
783  gstime = float(nmind)
784 
785 end subroutine get_julian_time
786 
787 end module radiance_mod
define_mod
Definition: define_mod.f90:1
define_mod::nstring
integer(i_kind), parameter nstring
Definition: define_mod.f90:15
define_mod::name_var_info
character(len=nstring), dimension(nvar_info) name_var_info
Definition: define_mod.f90:96
define_mod::type_sen_info
integer(i_kind), dimension(nsen_info) type_sen_info
Definition: define_mod.f90:142
define_mod::set_name_sensor
subroutine set_name_sensor(instid, sensor)
Definition: define_mod.f90:276
define_mod::ninst
integer(i_kind), parameter ninst
Definition: define_mod.f90:22
define_mod::nsen_info
integer(i_kind), parameter nsen_info
Definition: define_mod.f90:21
define_mod::missing_i
integer(i_kind), parameter missing_i
Definition: define_mod.f90:11
define_mod::ndatetime
integer(i_kind), parameter ndatetime
Definition: define_mod.f90:16
define_mod::name_sen_info
character(len=nstring), dimension(nsen_info) name_sen_info
Definition: define_mod.f90:132
define_mod::nvar_info
integer(i_kind), parameter nvar_info
Definition: define_mod.f90:20
define_mod::missing_r
real(r_kind), parameter missing_r
Definition: define_mod.f90:10
define_mod::inst_list
character(len=nstring), dimension(ninst) inst_list
Definition: define_mod.f90:70
define_mod::xdata
type(xdata_type), dimension(:), allocatable xdata
Definition: define_mod.f90:185
define_mod::set_name_satellite
subroutine set_name_satellite(satid, satellite)
Definition: define_mod.f90:245
define_mod::set_brit_obserr
subroutine set_brit_obserr(name_inst, ichan, obserr)
Definition: define_mod.f90:305
define_mod::type_var_info
integer(i_kind), dimension(nvar_info) type_var_info
Definition: define_mod.f90:108
kinds
Definition: kinds.f90:1
kinds::i_kind
integer, parameter, public i_kind
Definition: kinds.f90:71
kinds::r_kind
integer, parameter, public r_kind
Definition: kinds.f90:100
kinds::r_double
integer, parameter, public r_double
Definition: kinds.f90:80
radiance_mod
Definition: radiance_mod.f90:1
radiance_mod::fill_datalink
subroutine fill_datalink(datalink, rfill, ifill)
Definition: radiance_mod.f90:694
radiance_mod::rhead
type(datalink_radiance), pointer rhead
Definition: radiance_mod.f90:42
radiance_mod::sort_obs_radiance
subroutine, public sort_obs_radiance
Definition: radiance_mod.f90:518
radiance_mod::read_amsua_amsub_mhs
subroutine, public read_amsua_amsub_mhs(filename, filedate)
Definition: radiance_mod.f90:49
radiance_mod::read_airs_colocate_amsua
subroutine, public read_airs_colocate_amsua(filename, filedate)
Definition: radiance_mod.f90:243
radiance_mod::calc_sensor_view_angle
subroutine calc_sensor_view_angle(name_inst, ifov, view_angle)
Definition: radiance_mod.f90:725
radiance_mod::r8bfms
real(r_kind), parameter r8bfms
Definition: radiance_mod.f90:16
radiance_mod::rlink
type(datalink_radiance), pointer rlink
Definition: radiance_mod.f90:42
radiance_mod::get_julian_time
subroutine get_julian_time(year, month, day, hour, minute, gstime)
Definition: radiance_mod.f90:762
ufo_vars_mod
Definition: ufo_variables_mod.F90:8
ufo_vars_mod::ufo_vars_getindex
integer function, public ufo_vars_getindex(vars, varname)
Definition: ufo_variables_mod.F90:175
radiance_mod::datalink_radiance
Definition: radiance_mod.f90:18