decode_airs.f90
References to this file elsewhere.
1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 ! This program decodes a swath of AIRS L2 standard retrievals, computing RH
3 ! from q and writing the resulting soundings into a little_r formatted file
4 ! named soundings.little_r. The code is based on that from the example reader
5 ! programs supplied through the AIRS web site.
6 !
7 ! Michael Duda -- 26 May 2006
8 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
9 program decode_airs
10
11 use read_airs
12
13 implicit none
14
15 real, external :: calc_rh
16 integer, external :: iargc
17 logical, external :: granule_out_of_timewin
18
19 character (len=69), parameter :: rpt_format = '(2f20.5,4a40,f20.5,5i10,3L10,2i10,a20,13(f13.5,i7))'
20 character (len=15), parameter :: meas_format = '(10(f13.5, i7))'
21 character (len=5), parameter :: end_format = '(3i7)'
22 integer, parameter :: LESS = -1
23 integer, parameter :: EQUAL = 0
24 integer, parameter :: MORE = 1
25
26 ! Data stored in the header record
27 real :: latitude, longitude, elevation
28 real :: slp, psfc, refp, grndt, sst, x1, x2, x3, x4, x5, x6, x7, x8
29 integer :: islp, ipsfc, n_valid, n_err, n_warning, iseq, n_dups
30 integer :: isut, julian, irefp, igrndt, isst
31 integer :: ix1, ix2, ix3, ix4, ix5, ix6, ix7, ix8
32 character (len=20) :: date_char
33 character (len=40) :: id, name, platform, source
34 logical :: is_sound, is_bogus, is_discard
35
36 ! Data stored in the data record
37 real :: p, z, t, dewpt, spd, dir, u, v, rh, thk
38 integer :: ip, iz, it, idewpt, ispd, idir, iu, iv, irh, ithk
39
40 ! Data stored in the end record
41 integer :: num_valid, num_err, num_warn
42
43 type (airs_ret_gran_t) :: ret ! structure with entire granule
44 integer :: layer ! index for cloud or atmospheric layers
45 integer :: nargs ! number of arguments
46 integer :: track ! along-track index (scan number)
47 integer :: xtrack ! across-track index (FOV number)
48 character (len=256) :: arg ! buffer for argument
49 character (len=256) :: file_name ! name of AIRS Level 2 file to read
50 character (len=11) :: ref_time
51
52 character (len=14) :: min_time, max_time
53 namelist /time_window/ min_time, max_time
54
55 integer :: delta_time
56 integer :: istatus
57 integer :: nvalid
58 integer :: fn
59 integer :: cmp_min, cmp_max
60 integer :: mm, ss
61 integer :: qual ! quality. 0 is best, then 1. 2 is bad.
62 integer :: i, j ! indices 1..3 for 3x3 of IR FOVs per retrieval
63 real :: temp ! temperature or -9999 if not good enough
64 real :: h2o ! H2O MMR or -9999 if not good enough
65
66 nargs = iargc()
67
68 if (nargs <= 0) then
69 write(6,*) ' '
70 write(6,*) 'Usage: decode_airs.exe filename ...'
71 write(6,*) ' where filename is the name of an HDFEOS file containing'
72 write(6,*) ' a swath of L2 AIRS retreivals.'
73 write(6,*) ' '
74 stop
75 end if
76
77 min_time = '00000000000000'
78 max_time = '99999999999999'
79
80 open(11,file='time_window.nl',status='old',iostat=istatus)
81 if (istatus == 0) then
82 read(11,time_window)
83 end if
84
85 ! The time for each FOV is given as the number of seconds
86 ! since 1993 Jan 1 00Z
87 ref_time = '1993010100 '
88
89 open(10, file='soundings.little_r', form='formatted', status='replace')
90
91 do fn = 1, nargs
92 call getarg (fn, file_name)
93
94 ! Read all retrievals from file_name and place them in ret
95 call airs_ret_rdr(file_name, ret)
96
97 if (granule_out_of_timewin(ret%Time, min_time, max_time)) cycle
98
99 do track = 1, AIRS_RET_GEOTRACK
100 do xtrack = 1, AIRS_RET_GEOXTRACK
101 nvalid = 0
102
103 delta_time = nint(ret%Time(xtrack,track))/3600
104 date_char(1:6) = ' '
105 call geth_newdate(ref_time, delta_time, date_char(7:20))
106 mm = mod(nint(ret%Time(xtrack,track)),3600)/60
107 ss = mod(nint(ret%Time(xtrack,track)),60)
108 write(date_char(17:20),'(i2.2,i2.2)') mm, ss
109
110 if (ret%Time(xtrack,track) < 0.) then
111 cmp_min = LESS
112 else
113 call cmp_datestr(date_char(7:20), min_time, cmp_min)
114 call cmp_datestr(date_char(7:20), max_time, cmp_max)
115 end if
116
117
118 if (cmp_min /= LESS .and. cmp_max /= MORE) then
119
120 write (6,'(a,f8.4,a5,f9.4,a3)') 'Location: ', &
121 ret%Latitude(xtrack,track),'lat, ', &
122 ret%Longitude(xtrack,track),'lon'
123 write(6,'(a)') 'Time: '//date_char(7:20)
124 write(6,'(a,i9)') 'RetQAFlag: ', ret%RetQAFlag(xtrack,track)
125 write(6,*) ' '
126
127 !
128 ! First loop through all levels and find out how many valid obs there are
129 !
130 do layer=AIRS_RET_STDPRESSURELAY,1,-1
131
132 ! Which Qual flag to check depends on layer
133 if (ret%PressStd(layer) < ret%Press_mid_top_bndry(xtrack,track) ) then
134 qual = ret%Qual_Temp_Profile_Top(xtrack,track)
135
136 else if (ret%PressStd(layer) < ret%Press_bot_mid_bndry(xtrack,track) ) then
137 qual = ret%Qual_Temp_Profile_Mid(xtrack,track)
138
139 else
140 qual = ret%Qual_Temp_Profile_Bot(xtrack,track)
141
142 endif
143
144 ! Temperature. If quality is bad (2) then put out flag value of -9999.0
145 if (qual < 2 .and. layer > ret%nSurfStd(xtrack,track)) then
146 temp = ret%TAirStd(layer,xtrack,track)
147 h2o = ret%H2OMMRStd(layer,xtrack,track)
148
149 ! See p.3 of Level 2 Product Levels and Layers
150 ! else if (qual < 2 .and. abs(ret%PressStd(layer) - ret%PSurfStd(xtrack,track)) < 5.) then
151 ! temp = ret%TAirStd(layer,xtrack,track)
152 ! h2o = ret%H2OMMRStd(layer,xtrack,track)
153
154 else
155 temp = -9999.0
156 h2o = -9999.0
157
158 end if
159
160 if (temp /= -9999.) then
161 nvalid = nvalid + 1
162 if (h2o /= -9999.) nvalid = nvalid + 1
163 end if
164
165 end do
166
167
168 !
169 ! If there are valid obs to be reported, write them out to little_r format
170 !
171 if (nvalid > 0) then
172
173 latitude = ret%Latitude(xtrack,track)
174 longitude = ret%Longitude(xtrack,track)
175 id = 'AIRS L2 Std Retrieval '
176 name = 'AIRSRET '
177 platform = 'FM-133 '
178 source = 'GES DAAC '
179 elevation = -888888.
180 n_valid = nvalid
181 n_err = 0
182 n_warning = 0
183 iseq = 0
184 n_dups = 0
185 is_sound = .true.
186 is_bogus = .false.
187 is_discard = .false.
188 isut = -888888
189 julian = -888888
190 slp = -888888.
191 islp = -88
192 refp = -888888.
193 irefp = -88
194 grndt = -888888.
195 igrndt = -88
196 sst = -888888.
197 isst = -88
198 psfc = -888888.
199 ipsfc = -88
200 x1 = -888888.
201 ix1 = -88
202 x2 = -888888.
203 ix2 = -88
204 x3 = -888888.
205 ix3 = -88
206 x4 = -888888.
207 ix4 = -88
208 x5 = -888888.
209 ix5 = -88
210 x6 = -888888.
211 ix6 = -88
212 x7 = -888888.
213 ix7 = -88
214 x8 = -888888.
215 ix8 = -88
216
217 write(10, fmt=rpt_format) latitude, longitude, id, name, platform, &
218 source, elevation, n_valid, n_err, n_warning, iseq, n_dups, &
219 is_sound, is_bogus, is_discard, isut, julian, date_char, &
220 slp, islp, refp, irefp, grndt, igrndt, sst, isst, psfc, ipsfc, &
221 x1, ix1, x2, ix2, x3, ix3, x4, ix4, x5, ix5, x6, ix6, x7, ix7, &
222 x8, ix8
223
224 do layer=AIRS_RET_STDPRESSURELAY,1,-1
225
226 ! Which Qual flag to check depends on layer
227 if (ret%PressStd(layer) < ret%Press_mid_top_bndry(xtrack,track) ) then
228 qual = ret%Qual_Temp_Profile_Top(xtrack,track)
229
230 else if (ret%PressStd(layer) < ret%Press_bot_mid_bndry(xtrack,track) ) then
231 qual = ret%Qual_Temp_Profile_Mid(xtrack,track)
232
233 else
234 qual = ret%Qual_Temp_Profile_Bot(xtrack,track)
235
236 endif
237
238 ! Temperature. If quality is bad (2) then put out flag value of -888888.0
239 if (qual < 2 .and. layer > ret%nSurfStd(xtrack,track)) then
240 temp = ret%TAirStd(layer,xtrack,track)
241 if (temp == -9999.) temp = -888888.
242 h2o = ret%H2OMMRStd(layer,xtrack,track)
243 if (h2o == -9999.) h2o = -888888.
244
245 else
246 temp = -888888.0
247 h2o = -888888.0
248
249 end if
250
251 if (temp /= -888888.) then
252
253 p = ret%pressStd(layer)*100.
254 ip = 0
255 z = ret%GP_Height(layer,xtrack,track)
256 iz = 0
257 t = temp
258 it = 0
259 dewpt = -888888.
260 idewpt = 0
261 spd = -888888.
262 ispd = 0
263 dir = -888888.
264 idir = 0
265 u = -888888.
266 iu = 0
267 v = -888888.
268 iv = 0
269 if (h2o == -888888. .or. t == -888888.) then
270 rh = -888888.
271 else
272 rh = calc_rh(t, h2o/1000., p)
273 if (rh < 0.) rh = 0.
274 end if
275 irh = 0
276 thk = -888888.
277 ithk = 0
278 write(10, fmt=meas_format) p, ip, z, iz, t, it, dewpt, idewpt, &
279 spd, ispd, dir, idir, u, iu, v, iv, rh, irh, thk, ithk
280 end if
281 end do
282
283 p = -777777.
284 ip = 0
285 z = -777777.
286 iz = 0
287 t = -888888.
288 it = 0
289 dewpt = -888888.
290 idewpt = 0
291 spd = -888888.
292 ispd = 0
293 dir = -888888.
294 idir = 0
295 u = -888888.
296 iu = 0
297 v = -888888.
298 iv = 0
299 rh = -888888.
300 irh = 0
301 thk = -888888.
302 ithk = 0
303
304 write(10, fmt=meas_format) p, ip, z, iz, t, it, dewpt, idewpt, &
305 spd, ispd, dir, idir, u, iu, v, iv, rh, irh, thk, ithk
306
307 num_valid = nvalid
308 num_err = 0
309 num_warn = 0
310 write(10, fmt=end_format) num_valid, num_err, num_warn
311
312 end if
313
314 ! Need to output quality info so users can tell if we have a good
315 ! retrieval with no clouds or a bad retrieval. Both cases will put
316 ! all -9999s below.
317 ! write(6,*) '# Cloud quality:'
318 ! if (ret%Qual_Cloud_OLR(xtrack,track) .eq. 0) then
319 ! write(6,*) ' 0 Very Good'
320 ! else if (ret%Qual_Cloud_OLR(xtrack,track) .eq. 1) then
321 ! write(6,*) ' 1 Good'
322 ! else
323 ! write(6,*) ' 2 Do Not Use'
324 ! end if
325
326 end if ! Profile falls within time window
327
328 end do
329 end do
330
331 end do ! Loop over filenames
332
333 close(10)
334
335 stop
336
337 end program decode_airs
338
339 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
340 ! Given two date strings of the form YYYYMMDDHHmmss, icmp is set to
341 ! 0 if the two dates are the same,
342 ! -1 if date1 < date2, and
343 ! +1 if date1 > date2.
344 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
345 subroutine cmp_datestr(date1, date2, icmp)
346
347 implicit none
348
349 ! Arguments
350 character (len=14), intent(in) :: date1, date2
351 integer, intent(out) :: icmp
352
353 ! Local variables
354 integer :: yyyy1, yyyy2
355 integer :: mm1, mm2
356 integer :: dd1, dd2
357 integer :: hh1, hh2
358 integer :: min1, min2
359 integer :: sec1, sec2
360
361 read(date1(1:4), '(i4)') yyyy1
362 read(date1(5:6), '(i2)') mm1
363 read(date1(7:8), '(i2)') dd1
364 read(date1(9:10), '(i2)') hh1
365 read(date1(11:12), '(i2)') min1
366 read(date1(13:14), '(i2)') sec1
367
368 read(date2(1:4), '(i4)') yyyy2
369 read(date2(5:6), '(i2)') mm2
370 read(date2(7:8), '(i2)') dd2
371 read(date2(9:10), '(i2)') hh2
372 read(date2(11:12), '(i2)') min2
373 read(date2(13:14), '(i2)') sec2
374
375 if (yyyy1 > yyyy2) then
376 icmp = 1
377 return
378 else if (yyyy1 < yyyy2) then
379 icmp = -1
380 return
381 end if
382
383 if (mm1 > mm2) then
384 icmp = 1
385 return
386 else if (mm1 < mm2) then
387 icmp = -1
388 return
389 end if
390
391 if (dd1 > dd2) then
392 icmp = 1
393 return
394 else if (dd1 < dd2) then
395 icmp = -1
396 return
397 end if
398
399 if (hh1 > hh2) then
400 icmp = 1
401 return
402 else if (hh1 < hh2) then
403 icmp = -1
404 return
405 end if
406
407 if (min1 > min2) then
408 icmp = 1
409 return
410 else if (min1 < min2) then
411 icmp = -1
412 return
413 end if
414
415 if (sec1 > sec2) then
416 icmp = 1
417 return
418 else if (sec1 < sec2) then
419 icmp = -1
420 return
421 end if
422
423 icmp = 0
424 return
425
426 end subroutine cmp_datestr
427
428
429 function granule_out_of_timewin(time_arr, min_time, max_time)
430
431 use read_airs
432
433 implicit none
434
435 ! Parameters
436 integer, parameter :: LESS = -1
437 integer, parameter :: EQUAL = 0
438 integer, parameter :: MORE = 1
439
440 ! Arguments
441 double precision, dimension(AIRS_RET_GEOXTRACK,AIRS_RET_GEOTRACK), intent(in) :: time_arr
442 character (len=14), intent(in) :: min_time, max_time
443
444 ! Local variables
445 integer :: i, j
446 integer :: delta_time, cmp_min, cmp_max, mm, ss
447 character (len=11) :: ref_time
448 character (len=14) :: date_char
449
450 ! Return value
451 logical :: granule_out_of_timewin
452
453 granule_out_of_timewin = .true.
454 ref_time = '1993010100 '
455
456 do i=1,AIRS_RET_GEOXTRACK,AIRS_RET_GEOXTRACK-1
457 do j=1,AIRS_RET_GEOTRACK,AIRS_RET_GEOTRACK-1
458
459 delta_time = nint(time_arr(i,j))/3600
460 call geth_newdate(ref_time, delta_time, date_char)
461 mm = mod(nint(time_arr(i,j)),3600)/60
462 ss = mod(nint(time_arr(i,j)),60)
463 write(date_char(11:14),'(i2.2,i2.2)') mm, ss
464
465 if (time_arr(i,j) < 0.) then
466 cmp_min = LESS
467 else
468 call cmp_datestr(date_char, min_time, cmp_min)
469 call cmp_datestr(date_char, max_time, cmp_max)
470 end if
471
472 if (cmp_min /= LESS .and. cmp_max /= MORE) then
473 granule_out_of_timewin = .false.
474 return
475 end if
476
477 end do
478 end do
479
480 end function granule_out_of_timewin