da_read_obs.inc
References to this file elsewhere.
1 subroutine da_read_obs (ob, xp, filename)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Read the header of a GTS observation file
5 !
6 ! Modifications: 10/26/2004 Syed RH Rizvi
7 !
8 ! Fix Obs-Long as -180 if it is +180
9 !
10 ! Date: 03/04/2005 - Syed RH Rizvi
11 !
12 ! (a) AMVs from Geostationary and Polar orbiting satellite are
13 ! seperated & used as profile. Currently there is a provision
14 ! to use MODIS winds only.
15 ! (b) MODIS obs error are currently assigned through namelist
16 ! parameter (modis_cmv_error)
17 !
18 ! 03/30/2005 Syed RH Rizvi
19 ! For global option duplicate obs at East/West boundary
20 !
21 ! 08/30/2005 Syed RH Rizvi
22 ! Writing original errors & thicknesses
23 ! desired for outputting QC obs with NTMAX = 0
24 !
25 !---------------------------------------------------------------------------
26
27 implicit none
28
29 type (xpose_type), intent(in) :: xp ! Domain decomposition vars.
30 type (ob_type), intent(inout) :: ob
31
32 character(*), intent(in) :: filename
33
34 character (LEN = 10) :: fmt_name
35
36 character (LEN = 160) :: info_string
37
38 character (LEN = 160) :: fmt_info, &
39 fmt_loc, &
40 fmt_each
41
42 integer :: i, j, iost, nlevels, fm,iunit
43
44 type (multi_level_type) :: platform
45
46 logical :: outside
47 logical :: outside_all
48 integer :: surface_level
49 real :: height_error, u_comp, v_comp
50 integer :: total_obs, &
51 num_sound, &
52 num_synop, &
53 num_pilot, &
54 num_geoamv, &
55 num_polaramv, &
56 num_satem, &
57 num_airep, &
58 num_metar, &
59 num_ships, &
60 num_gpspw, &
61 num_gpsref, &
62 num_ssmi_retrieval, &
63 num_ssmi_tb , &
64 num_ssmt1, &
65 num_ssmt2, &
66 num_qscat, &
67 num_profiler, &
68 num_airsr,num_buoy, num_bogus
69
70 integer :: ndup, n, report
71
72 if (trace_use) call da_trace_entry("da_read_obs")
73
74 total_obs = ob%total_obs
75 num_sound = ob%num_sound
76 num_synop = ob%num_synop
77 num_pilot = ob%num_pilot
78 num_geoamv = ob%num_geoamv
79 num_polaramv = ob%num_polaramv
80 num_satem = ob%num_satem
81 num_airep = ob%num_airep
82 num_metar = ob%num_metar
83 num_ships = ob%num_ships
84 num_gpspw = ob%num_gpspw
85 num_gpsref = ob%num_gpsref
86 num_ssmi_retrieval = ob%num_ssmi_retrieval
87 num_ssmi_tb = ob%num_ssmi_tb
88 num_ssmt1 = ob%num_ssmt1
89 num_ssmt2 = ob%num_ssmt2
90 num_qscat = ob%num_qscat
91 num_profiler = ob%num_profiler
92 num_buoy = ob%num_buoy
93 num_bogus= ob%num_bogus
94 num_airsr= ob%num_airsr
95
96 ! open file
97 ! =========
98
99 call da_get_unit(iunit)
100 open(unit = iunit, &
101 FILE = trim(filename), &
102 FORM = 'FORMATTED', &
103 ACCESS = 'SEQUENTIAL', &
104 iostat = iost, &
105 STATUS = 'OLD')
106
107 if (iost /= 0) then
108 write(unit=message(1),fmt='(A,I5,A)') &
109 "Error",iost," opening gts obs file "//trim(filename)
110 call da_warning(__FILE__,__LINE__,message(1:1))
111 return
112 end if
113
114 ! read header
115 ! ===========
116
117 do
118 read (unit = iunit, fmt = '(A)', iostat = iost) info_string
119 if (iost /= 0) then
120 call da_warning(__FILE__,__LINE__, &
121 (/"Problem reading gts obs header, skipping file"/))
122 return
123 end if
124
125 if (info_string(1:6) == 'EACH ') exit
126 end do
127
128 ! read formats
129 ! ------------
130
131 read (iunit, fmt = '(A,1X,A)', iostat = iost) &
132 fmt_name, fmt_info, &
133 fmt_name, fmt_loc, &
134 fmt_name, fmt_each
135
136 if (iost /= 0) then
137 call da_warning(__FILE__,__LINE__, &
138 (/"Problem reading gts obs file, skipping"/))
139 return
140 end if
141
142 if (print_detail_obs) then
143 write(unit=stdout, fmt='(2a)') &
144 'fmt_info=', fmt_info, &
145 'fmt_loc =', fmt_loc, &
146 'fmt_each=', fmt_each
147 end if
148
149 ! skip line
150 ! ----------
151
152 read (iunit, fmt = '(a)') fmt_name
153
154 ! loop over records
155 ! -----------------
156
157 report = 0 ! report number in file
158
159 reports: &
160 do
161
162 report = report+1
163
164 ! read station general info
165 ! =============================
166
167 read (iunit, fmt = fmt_info, iostat = iost) &
168 platform%info%platform, &
169 platform%info%date_char, &
170 platform%info%name, &
171 platform%info%levels, &
172 platform%info%lat, &
173 platform%info%lon, &
174 platform%info%elv, &
175 platform%info%id
176
177 if (iost /= 0) then
178 ! JRB This is expected, but its unclear how we handle failure
179 ! here without assuming the fortran2003 convention on
180 ! error statuses
181 exit reports
182 end if
183
184 if (print_detail_obs) then
185 write(unit=stdout, fmt=fmt_info) &
186 platform%info%platform, &
187 platform%info%date_char, &
188 platform%info%name, &
189 platform%info%levels, &
190 platform%info%lat, &
191 platform%info%lon, &
192 platform%info%elv, &
193 platform%info%id
194 end if
195
196 if (platform%info%lon == 180.0 ) platform%info%lon =-180.000
197 ! Fix funny wind direction at Poles
198 if (platform%info%lat < -89.9999 .or. platform%info%lat > 89.9999) then
199 platform%info%lon = 0.0
200 end if
201
202 if (platform%info%platform(6:6) == ' ') then
203 read(platform%info%platform(4:5), '(I2)') fm
204 else
205 read(platform%info%platform(4:6), '(I3)') fm
206 end if
207
208 ! read model location
209 ! =========================
210
211 read (iunit, fmt = fmt_loc) platform%loc%slp, platform%loc%pw
212 total_obs = total_obs + 1
213
214 ! levels < 1 and not GPSPW, go back to reports
215
216 if ((platform%info%levels < 1) .AND. &
217 (index(platform%info%platform, 'GPSPW') <= 0)) then
218 cycle reports
219 end if
220
221 ! read each level
222 ! ---------------
223
224 do i = 1, platform%info%levels
225 platform%each (i) = each_level_type(missing_r, missing, -1.0, & ! height
226 field_type(missing_r, missing, missing_r), & ! u
227 field_type(missing_r, missing, missing_r), & ! v
228 field_type(missing_r, missing, missing_r), & ! p
229 field_type(missing_r, missing, missing_r), & ! t
230 field_type(missing_r, missing, missing_r), & ! q
231 field_type(missing_r, missing, missing_r), & ! rh
232 field_type(missing_r, missing, missing_r), & ! td
233 field_type(missing_r, missing, missing_r)) ! speed
234
235 read (unit = iunit, fmt = trim (fmt_each)) &
236 platform%each (i)%p, &
237 platform%each (i)%speed, &
238 ! Here the 'direction' is stored in platform%each (i)%v:
239 platform%each (i)%v, &
240 platform%each (i)%height, &
241 platform%each (i)%height_qc, &
242 height_error, &
243 platform%each (i)%t, &
244 platform%each (i)%td, &
245 platform%each (i)%rh
246
247 if (print_detail_obs) then
248 write(unit=stdout, fmt = '(a, i6)') 'i=', i
249 write(unit=stdout, fmt = trim (fmt_each)) &
250 platform%each (i)%p, &
251 platform%each (i)%speed, &
252 platform%each (i)%v, &
253 platform%each (i)%height, &
254 platform%each (i)%height_qc, &
255 height_error, &
256 platform%each (i)%t, &
257 platform%each (i)%td, &
258 platform%each (i)%rh
259 end if
260
261 ! Uncomment the following if errors in reported heights (test later):
262 ! if (fm == 97 .or. fm == 96 .or. fm == 42) then
263 ! platform%each (i)%height_qc = -5 ! Aircraft heights are not above surface.
264 ! end if
265
266 ! To convert the wind speed and direction to
267 ! the model wind components: u, v
268
269 if (platform%each (i)%speed%qc /= missing_data .and. &
270 platform%each (i)%v%qc /= missing_data) then
271 call da_ffdduv (platform%each (i)%speed%inv, &
272 platform%each (i)%v%inv, &
273 U_comp, V_comp, platform%info%lon, convert_fd2uv)
274 platform%each (i)%u = platform%each (i)%speed
275 platform%each (i)%v = platform%each (i)%speed
276 platform%each (i)%u%inv = U_comp
277 platform%each (i)%v%inv = V_comp
278 else
279 platform%each (i)%u%inv = missing_r
280 platform%each (i)%v%inv = missing_r
281 platform%each (i)%u%error = missing_r
282 platform%each (i)%v%error = missing_r
283 platform%each (i)%u%qc = missing_data
284 platform%each (i)%v%qc = missing_data
285 end if
286 end do
287
288 ! Restrict to a range of reports, useful for debugging
289
290 if (report < report_start) then
291 cycle
292 end if
293
294 if (report > report_end) then
295 exit
296 end if
297
298 call da_ll_to_xy (platform%info, platform%loc, &
299 xp, outside, outside_all)
300
301 if (outside_all) then
302 cycle reports
303 end if
304
305 if (print_detail_obs) then
306 ! Simplistic approach, could be improved to get it all done on PE 0
307 if (.NOT. outside) then
308 write(unit=stdout,fmt='(A,I4,A,F8.2,A,F8.2,A,I3,A,2F8.2)') &
309 "Report",report," at",platform%info%lon,"E",platform%info%lat, &
310 "N on processor", myproc,", position", platform%loc%x,platform%loc%y
311 end if
312 end if
313
314 nlevels = platform%info%levels
315 if (nlevels > max_ob_levels) then
316 nlevels = max_ob_levels
317 write(unit=message(1), fmt='(/4(2a,2x),a,2f8.2,2x,2(a,f9.2,2x),2(a,i4,2x)/)') &
318 'Station: ', trim(platform%info%name), &
319 'ID: ', trim(platform%info%id), &
320 'Platform: ', trim(platform%info%platform), &
321 'Date: ', trim(platform%info%date_char), &
322 'At Loc(lat, lon): ', platform%info%lat, platform%info%lon, &
323 'At elv: ', platform%info%elv, &
324 'with pstar: ', platform%info%pstar, &
325 'Has level: ', platform%info%levels, &
326 'which is great than max_ob_levels: ', max_ob_levels
327
328 write (unit=message(2), fmt = '(A,1X,A,1X,A,1X,I4,2f9.3,f9.1,1X,A)') &
329 platform%info%name, &
330 platform%info%platform, &
331 platform%info%date_char, &
332 platform%info%levels, &
333 platform%info%lat, &
334 platform%info%lon, &
335 platform%info%elv, &
336 platform%info%id
337 call da_warning(__FILE__,__LINE__,message(1:2))
338
339 platform%info%levels = nlevels
340 else if (nlevels < 1) then
341 ! Not GPSPW and GPSZD data:
342 if (fm /= 111 .and. fm /= 114) then
343 cycle reports
344 end if
345 end if
346
347 if (fm /= 86) call da_obs_proc_station(platform)
348
349 nlevels = platform%info%levels
350 ! Loop over duplicating obs for global
351 ndup = 1
352 if (global .and. &
353 (platform%loc%i < xp%ids .or. platform%loc%i >= xp%ide)) ndup= 2
354
355 ! It is possible that logic for counting obs is incorrect for the
356 ! global case with >1 MPI tasks due to obs duplication, halo, etc.
357 ! TBH: 20050913
358
359 if (.not.outside) then
360 if (print_detail_obs .and. ndup > 1) then
361 write(unit=stdout, fmt = fmt_info) &
362 platform%info%platform, &
363 platform%info%date_char, &
364 platform%info%name, &
365 platform%info%levels, &
366 platform%info%lat, &
367 platform%info%lon, &
368 platform%info%elv, &
369 platform%info%id
370
371 write(unit=stdout, fmt = '(a,2i5,4e20.10)') &
372 ' duplicating obs since loc% i,j,dx,dxm,dy & dym ', &
373 platform%loc%i, platform%loc%j, &
374 platform%loc%dx, platform%loc%dxm, &
375 platform%loc%dy, platform%loc%dym
376 end if
377 end if
378
379 dup_loop: do n = 1, ndup
380 select case(fm)
381
382 case (12) ;
383
384 if (n==1) ob%num_synop_glo = ob%num_synop_glo + 1
385 if (outside) then
386 cycle reports
387 end if
388
389 if (.not.use_SynopObs) cycle reports
390
391 num_synop = num_synop + 1
392 ! Track serial obs index for reassembly of obs during bit-for-bit
393 ! tests with different numbers of MPI tasks.
394 platform%loc%obs_global_index = ob%num_synop_glo
395
396 ob%synop(num_synop)%info = platform%info
397 ob%synop(num_synop)%loc = platform%loc
398
399 ob%synop(num_synop)%h = platform%each(1)%height
400 ob%synop(num_synop)%u = platform%each(1)%u
401 ob%synop(num_synop)%v = platform%each(1)%v
402 ob%synop(num_synop)%t = platform%each(1)%t
403 ob%synop(num_synop)%q = platform%each(1)%q
404 ob%synop(num_synop)%p = platform%each(1)%p
405
406 case (13, 17) ; ! ships
407
408 if (n == 1) ob%num_ships_glo = ob%num_ships_glo + 1
409 if (outside) then
410 cycle reports
411 end if
412
413 if (.not.use_ShipsObs) cycle reports
414
415 num_ships = num_ships + 1
416
417 ! Track serial obs index for reassembly of obs during bit-for-bit
418 ! tests with different numbers of MPI tasks.
419 platform%loc%obs_global_index = ob%num_ships_glo
420
421 ob%ships(num_ships)%info = platform%info
422 ob%ships(num_ships)%loc = platform%loc
423
424 ob%ships(num_ships)%h = platform%each(1)%height
425 ob%ships(num_ships)%u = platform%each(1)%u
426 ob%ships(num_ships)%v = platform%each(1)%v
427 ob%ships(num_ships)%t = platform%each(1)%t
428 ob%ships(num_ships)%p = platform%each(1)%p
429 ob%ships(num_ships)%q = platform%each(1)%q
430
431 case (15:16) ;
432
433 if (n==1) ob%num_metar_glo = ob%num_metar_glo + 1
434 if (outside) then
435 cycle reports
436 end if
437
438 if (.not.use_MetarObs) cycle reports
439
440 num_metar = num_metar + 1
441 ! Track serial obs index for reassembly of obs during bit-for-bit
442 ! tests with different numbers of MPI tasks.
443 platform%loc%obs_global_index = ob%num_metar_glo
444
445 ob%metar(num_metar)%info = platform%info
446 ob%metar(num_metar)%loc = platform%loc
447
448 ob%metar(num_metar)%h = platform%each(1)%height
449 ob%metar(num_metar)%u = platform%each(1)%u
450 ob%metar(num_metar)%v = platform%each(1)%v
451 ob%metar(num_metar)%t = platform%each(1)%t
452 ob%metar(num_metar)%p = platform%each(1)%p
453 ob%metar(num_metar)%q = platform%each(1)%q
454
455 case (32:34) ;
456
457 if (n==1) ob%num_pilot_glo = ob%num_pilot_glo + 1
458 if (outside) then
459 cycle reports
460 end if
461
462 if (.not.use_PilotObs) cycle reports
463
464 num_pilot = num_pilot + 1
465 ! Track serial obs index for reassembly of obs during bit-for-bit
466 ! tests with different numbers of MPI tasks.
467 platform%loc%obs_global_index = ob%num_pilot_glo
468
469 ob%pilot(num_pilot)%info = platform%info
470 ob%pilot(num_pilot)%loc = platform%loc
471
472 allocate (ob%pilot(num_pilot)%p (1:nlevels))
473 allocate (ob%pilot(num_pilot)%zk(1:nlevels))
474 allocate (ob%pilot(num_pilot)%u (1:nlevels))
475 allocate (ob%pilot(num_pilot)%v (1:nlevels))
476
477 do i = 1, nlevels
478 ob%pilot(num_pilot)%p(i) = platform%each(i)%p%inv
479 ob%pilot(num_pilot)%u(i) = platform%each(i)%u
480 ob%pilot(num_pilot)%v(i) = platform%each(i)%v
481 end do
482
483 case (35:38) ;
484
485 if (n==1) ob%num_sound_glo = ob%num_sound_glo + 1
486 if (outside) then
487 cycle reports
488 end if
489
490 if (.not.use_SoundObs) cycle reports
491
492 num_sound = num_sound + 1
493 ! Track serial obs index for reassembly of obs during bit-for-bit
494 ! tests with different numbers of MPI tasks.
495 platform%loc%obs_global_index = ob%num_sound_glo
496 platform%loc%obs_global_report = report
497
498 ob%sonde_sfc(num_sound)%info = platform%info
499 ob%sonde_sfc(num_sound)%loc = platform%loc
500
501 ! Search to see if we have surface obs.
502
503 surface_level = 0
504
505 do i = 1, nlevels
506 ! if (elevation and height are the same, it is surface.)
507 if (abs(platform%info%elv - platform%each(i)%height) < 0.1) then
508 surface_level = i
509
510 ! Save surface pressure.
511 ob%sonde_sfc(num_sound)%h = platform%each(i)%height
512 ob%sonde_sfc(num_sound)%u = platform%each(i)%u
513 ob%sonde_sfc(num_sound)%v = platform%each(i)%v
514 ob%sonde_sfc(num_sound)%t = platform%each(i)%t
515 ob%sonde_sfc(num_sound)%q = platform%each(i)%q
516 ob%sonde_sfc(num_sound)%p = platform%each(i)%p
517
518 exit
519 end if
520 end do
521
522 ! processing the sound_sfc data:
523
524 ob%sound(num_sound)%info = platform%info
525 ob%sound(num_sound)%loc = platform%loc
526
527 if (surface_level > 0) then
528 nlevels = nlevels - 1
529 else
530 ob%sonde_sfc(num_sound)%h = missing_r
531 ob%sonde_sfc(num_sound)%u%inv = missing_r
532 ob%sonde_sfc(num_sound)%u%qc = missing
533 ob%sonde_sfc(num_sound)%u%error = abs(missing_r)
534 ob%sonde_sfc(num_sound)%v = ob%sonde_sfc(num_sound)%u
535 ob%sonde_sfc(num_sound)%t = ob%sonde_sfc(num_sound)%u
536 ob%sonde_sfc(num_sound)%p = ob%sonde_sfc(num_sound)%u
537 ob%sonde_sfc(num_sound)%q = ob%sonde_sfc(num_sound)%u
538 end if
539
540 if (nlevels > 0) then
541
542 allocate (ob%sound(num_sound)%h (1:nlevels))
543 allocate (ob%sound(num_sound)%p (1:nlevels))
544 allocate (ob%sound(num_sound)%zk(1:nlevels))
545 allocate (ob%sound(num_sound)%u (1:nlevels))
546 allocate (ob%sound(num_sound)%v (1:nlevels))
547 allocate (ob%sound(num_sound)%t (1:nlevels))
548 allocate (ob%sound(num_sound)%q (1:nlevels))
549
550 j = 0
551 do i = 1, ob%sound(num_sound)%info%levels
552 if (i == surface_level) cycle
553
554 j=j+1
555
556 ob%sound(num_sound)%h(j) = platform%each(i)%height
557 ob%sound(num_sound)%p(j) = platform%each(i)%p%inv
558 ob%sound(num_sound)%u(j) = platform%each(i)%u
559 ob%sound(num_sound)%v(j) = platform%each(i)%v
560 ob%sound(num_sound)%t(j) = platform%each(i)%t
561 ob%sound(num_sound)%q(j) = platform%each(i)%q
562 end do
563 end if
564
565 ob%sound(num_sound)%info%levels = nlevels
566
567 case (86) ;
568
569 if (n==1) ob%num_satem_glo = ob%num_satem_glo + 1
570 if (outside) then
571 cycle reports
572 end if
573
574 if (.not.use_SatemObs) cycle reports
575
576 ! Reject cloudy Satem obs.
577
578 if (platform%loc%pw%inv > 10.) then
579 cycle reports
580 end if
581
582 num_satem = num_satem + 1
583 ! Track serial obs index for reassembly of obs during bit-for-bit
584 ! tests with different numbers of MPI tasks.
585 platform%loc%obs_global_index = ob%num_satem_glo
586
587 ob%satem(num_satem)%info = platform%info
588 ob%satem(num_satem)%loc = platform%loc
589
590 ! The Satem ref_p is put in the SLP position in OBSPROC data stream.
591
592 ob%satem(num_satem)%ref_p= platform%loc%slp%inv
593
594 allocate (ob%satem(num_satem)%p (1:nlevels))
595 allocate (ob%satem(num_satem)%thickness(1:nlevels))
596 allocate (ob%satem(num_satem)%org_thickness(1:nlevels))
597
598 ob%satem(num_satem)%p(1) = platform%each(1)%p%inv
599 ob%satem(num_satem)%thickness(1) = platform%each(1)%t
600 ! write original thickness errors for filtered obs
601 if (anal_type_qcobs) then
602 do i = 1, nlevels
603 ob%satem(num_satem)%org_thickness(i) = platform%each(i)%t
604 end do
605 end if
606
607 ! Splitting the reported Satem data into smaller layers.
608
609 do i = 2, nlevels
610 ob%satem(num_satem)%p(i) = platform%each(i)%p%inv
611 ob%satem(num_satem)%thickness(i) = platform%each(i)%t
612 if (platform%each(i)%t%qc /= missing_data .and. &
613 platform%each(i-1)%t%qc /= missing_data) then
614 ob%satem(num_satem)%thickness(i)%inv = &
615 platform%each(i)%t%inv - platform%each(i-1)%t%inv
616 else
617 ob%satem(num_satem)%thickness(i)%inv = missing_r
618 ob%satem(num_satem)%thickness(i)%qc = missing_data
619 end if
620 end do
621
622 ! Thickness error (m):
623
624 do i = nlevels, 2, -1
625 ob%satem(num_satem)%thickness(i)%error = &
626 sqrt(ob%satem(num_satem)%thickness(i )%error ** 2 + &
627 ob%satem(num_satem)%thickness(i-1)%error ** 2) / gravity
628 end do
629
630 ob%satem(num_satem)%thickness(1)%error = &
631 sqrt(ob%satem(num_satem)%thickness(1)%error ** 2 + &
632 platform%loc%pw%error ** 2) / gravity
633
634 ! Geostatinary ot Polar orbitting Satellite AMVs:
635
636 case (88) ;
637
638 if (index(platform%info%name, 'MODIS') > 0 .or. &
639 index(platform%info%name, 'modis') > 0) then
640
641 if (n==1) ob%num_polaramv_glo = ob%num_polaramv_glo + 1
642 if (outside) then
643 cycle reports
644 end if
645
646 if (.not.use_polaramvObs) cycle reports
647
648 num_polaramv = num_polaramv + 1
649 ! Track serial obs index for reassembly of obs during bit-for-bit
650 ! tests with different numbers of MPI tasks.
651 platform%loc%obs_global_index = ob%num_polaramv_glo
652
653 ob%polaramv(num_polaramv)%info = platform%info
654 ob%polaramv(num_polaramv)%loc = platform%loc
655
656 allocate (ob%polaramv(num_polaramv)%p (1:nlevels))
657 allocate (ob%polaramv(num_polaramv)%zk(1:nlevels))
658 allocate (ob%polaramv(num_polaramv)%u (1:nlevels))
659 allocate (ob%polaramv(num_polaramv)%v (1:nlevels))
660
661 do i = 1, nlevels
662 ob%polaramv(num_polaramv)%p(i) = platform%each(i)%p%inv
663 ob%polaramv(num_polaramv)%u(i) = platform%each(i)%u
664 ob%polaramv(num_polaramv)%v(i) = platform%each(i)%v
665 end do
666 else
667 if (n==1) ob%num_geoamv_glo = ob%num_geoamv_glo + 1
668 if (outside) then
669 cycle reports
670 end if
671
672 if (.not.use_geoamvObs) cycle reports
673
674 num_geoamv = num_geoamv + 1
675 ! Track serial obs index for reassembly of obs during bit-for-bit
676 ! tests with different numbers of MPI tasks.
677 platform%loc%obs_global_index = ob%num_geoamv_glo
678
679 ob%geoamv(num_geoamv)%info = platform%info
680 ob%geoamv(num_geoamv)%loc = platform%loc
681
682 allocate (ob%geoamv(num_geoamv)%p (1:nlevels))
683 allocate (ob%geoamv(num_geoamv)%zk(1:nlevels))
684 allocate (ob%geoamv(num_geoamv)%u (1:nlevels))
685 allocate (ob%geoamv(num_geoamv)%v (1:nlevels))
686
687 do i = 1, nlevels
688 ob%geoamv(num_geoamv)%p(i) = platform%each(i)%p%inv
689 ob%geoamv(num_geoamv)%u(i) = platform%each(i)%u
690 ob%geoamv(num_geoamv)%v(i) = platform%each(i)%v
691 end do
692 end if
693
694 case (42,96:97) ;
695
696 if (n==1) ob%num_airep_glo = ob%num_airep_glo + 1
697 if (outside) then
698 cycle reports
699 end if
700
701 if (.not.use_AirepObs) cycle reports
702
703 num_airep = num_airep + 1
704 ! Track serial obs index for reassembly of obs during bit-for-bit
705 ! tests with different numbers of MPI tasks.
706 platform%loc%obs_global_index = ob%num_airep_glo
707
708 ob%airep(num_airep)%info = platform%info
709 ob%airep(num_airep)%loc = platform%loc
710
711 allocate (ob%airep(num_airep)%h (1:nlevels))
712 allocate (ob%airep(num_airep)%p (1:nlevels))
713 allocate (ob%airep(num_airep)%zk (1:nlevels))
714 allocate (ob%airep(num_airep)%u (1:nlevels))
715 allocate (ob%airep(num_airep)%v (1:nlevels))
716 allocate (ob%airep(num_airep)%t (1:nlevels))
717
718 do i = 1, nlevels
719 ob%airep(num_airep)%h(i) = platform%each(i)%height
720 ob%airep(num_airep)%p(i) = platform%each(i)%p%inv
721 ob%airep(num_airep)%u(i) = platform%each(i)%u
722 ob%airep(num_airep)%v(i) = platform%each(i)%v
723 ob%airep(num_airep)%t(i) = platform%each(i)%t
724 end do
725
726 case (111, 114) ;
727
728 if (n==1) ob%num_gpspw_glo = ob%num_gpspw_glo + 1
729 if (outside) then
730 cycle reports
731 end if
732
733 if (.not.use_GpspwObs) cycle reports
734
735 num_gpspw = num_gpspw + 1
736 ! Track serial obs index for reassembly of obs during bit-for-bit
737 ! tests with different numbers of MPI tasks.
738 platform%loc%obs_global_index = ob%num_gpspw_glo
739
740 ob%gpspw(num_gpspw)%info = platform%info
741 ob%gpspw(num_gpspw)%loc = platform%loc
742 ob%gpspw(num_gpspw)%tpw = platform%loc%pw
743
744 case (116) ;
745
746 if (n==1) ob%num_gpsref_glo = ob%num_gpsref_glo + 1
747 if (outside) then
748 cycle reports
749 end if
750
751 if (.not.use_GpsrefObs) cycle reports
752
753 num_gpsref = num_gpsref + 1
754 ! Track serial obs index for reassembly of obs during bit-for-bit
755 ! tests with different numbers of MPI tasks.
756 platform%loc%obs_global_index = ob%num_gpsref_glo
757
758 ob%gpsref(num_gpsref)%info = platform%info
759 ob%gpsref(num_gpsref)%loc = platform%loc
760
761 allocate (ob%gpsref(num_gpsref)%h (1:nlevels))
762 allocate (ob%gpsref(num_gpsref)%zk(1:nlevels))
763 allocate (ob%gpsref(num_gpsref)%ref(1:nlevels))
764 allocate (ob%gpsref(num_gpsref)% p(1:nlevels))
765 allocate (ob%gpsref(num_gpsref)% t(1:nlevels))
766 allocate (ob%gpsref(num_gpsref)% q(1:nlevels))
767
768 do i = 1, nlevels
769 ob%gpsref(num_gpsref)%h(i) = platform%each(i)%height
770
771 ! In OBSPROC, use "td" field to store "gpsref"
772 ob%gpsref(num_gpsref)%ref(i) = platform%each(i)%td
773 ! Keep the retrieved p and t (and q) as "field_type":
774 ob%gpsref(num_gpsref)%p(i) = platform%each(i)%p
775 ob%gpsref(num_gpsref)%t(i) = platform%each(i)%t
776 ob%gpsref(num_gpsref)%q(i) = platform%each(i)%q
777 end do
778
779 case (121) ;
780 ! SSM/T1 temperatures:
781
782 if (n==1) ob%num_ssmt1_glo = ob%num_ssmt1_glo + 1
783 if (outside) then
784 cycle reports
785 end if
786
787 if (.not.use_ssmt1obs) cycle reports
788
789 num_ssmt1 = num_ssmt1 + 1
790 ! Track serial obs index for reassembly of obs during bit-for-bit
791 ! tests with different numbers of MPI tasks.
792 platform%loc%obs_global_index = ob%num_ssmt1_glo
793
794 ob%ssmt1(num_ssmt1)%info = platform%info
795 ob%ssmt1(num_ssmt1)%loc = platform%loc
796
797 allocate (ob%ssmt1(num_ssmt1)%h (1:nlevels))
798 allocate (ob%ssmt1(num_ssmt1)%p (1:nlevels))
799 allocate (ob%ssmt1(num_ssmt1)%t (1:nlevels))
800 allocate (ob%ssmt1(num_ssmt1)%zk(1:nlevels))
801
802 do i = 1, nlevels
803 ob%ssmt1(num_ssmt1)%h(i) = platform%each(i)%height
804 ob%ssmt1(num_ssmt1)%p(i) = platform%each(i)%p%inv
805 ob%ssmt1(num_ssmt1)%t(i) = platform%each(i)%t
806 end do
807
808 case (122) ;
809 ! SSM/T2 relative humidities:
810
811 if (n==1) ob%num_ssmt2_glo = ob%num_ssmt2_glo + 1
812 if (outside) then
813 cycle reports
814 end if
815
816 if (.not.use_ssmt2obs) cycle reports
817
818 num_ssmt2 = num_ssmt2 + 1
819 ! Track serial obs index for reassembly of obs during bit-for-bit
820 ! tests with different numbers of MPI tasks.
821 platform%loc%obs_global_index = ob%num_ssmt2_glo
822
823 ob%ssmt2(num_ssmt2)%info = platform%info
824 ob%ssmt2(num_ssmt2)%loc = platform%loc
825
826 allocate (ob%ssmt2(num_ssmt2)%h (1:nlevels))
827 allocate (ob%ssmt2(num_ssmt2)%p (1:nlevels))
828 allocate (ob%ssmt2(num_ssmt2)%zk(1:nlevels))
829 allocate (ob%ssmt2(num_ssmt2)%rh(1:nlevels))
830
831 do i = 1, nlevels
832 ob%ssmt2(num_ssmt2)% h(i) = platform%each(i)%height
833 ob%ssmt2(num_ssmt2)% p(i) = platform%each(i)%p%inv
834 ob%ssmt2(num_ssmt2)%rh(i) = platform%each(i)%rh
835 end do
836
837 case (281) ;
838 ! Scatterometer:
839
840 if (n==1) ob%num_qscat_glo = ob%num_qscat_glo + 1
841 if (outside) then
842 cycle reports
843 end if
844
845 if (.not.use_qscatobs) cycle reports
846
847 num_qscat = num_qscat + 1
848 ! Track serial obs index for reassembly of obs during bit-for-bit
849 ! tests with different numbers of MPI tasks.
850 platform%loc%obs_global_index = ob%num_qscat_glo
851
852 ob%qscat(num_qscat)%info = platform%info
853 ob%qscat(num_qscat)%loc = platform%loc
854
855 ob%qscat(num_qscat)%h = platform%each(1)%height
856 ob%qscat(num_qscat)%u = platform%each(1)%u
857 ob%qscat(num_qscat)%v = platform%each(1)%v
858
859 ! Impose minimum observation error = 1.0m/s for Quikscat data:
860 ob%qscat(num_qscat)%u%error = max(platform%each(1)%u%error,1.0)
861 ob%qscat(num_qscat)%v%error = max(platform%each(1)%v%error,1.0)
862
863 case (132) ; ! profiler
864
865 if (n==1) ob%num_profiler_glo = ob%num_profiler_glo + 1
866 if (outside) then
867 cycle reports
868 end if
869
870 if (.not.use_ProfilerObs) cycle reports
871
872 num_profiler = num_profiler + 1
873 ! Track serial obs index for reassembly of obs during bit-for-bit
874 ! tests with different numbers of MPI tasks.
875 platform%loc%obs_global_index = ob%num_profiler_glo
876
877 ob%profiler(num_profiler)%info = platform%info
878 ob%profiler(num_profiler)%loc = platform%loc
879
880 allocate (ob%profiler(num_profiler)%p (1:nlevels))
881 allocate (ob%profiler(num_profiler)%zk(1:nlevels))
882 allocate (ob%profiler(num_profiler)%u (1:nlevels))
883 allocate (ob%profiler(num_profiler)%v (1:nlevels))
884
885 do i = 1, nlevels
886 ob%profiler(num_profiler)%p(i) = platform%each(i)%p%inv
887 ob%profiler(num_profiler)%u(i) = platform%each(i)%u
888 ob%profiler(num_profiler)%v(i) = platform%each(i)%v
889 end do
890
891 case (135) ; ! Bogus
892
893 if (n==1) ob%num_bogus_glo = ob%num_bogus_glo + 1
894 if (outside) then
895 cycle reports
896 end if
897
898 if (.not.use_BogusObs) cycle reports
899
900 num_bogus = num_bogus + 1
901 ! Track serial obs index for reassembly of obs during bit-for-bit
902 ! tests with different numbers of MPI tasks.
903 platform%loc%obs_global_index = ob%num_bogus_glo
904
905 if (num_bogus > max_bogus_input) then
906 write(unit=message(1),fmt='(A,I6,A,I6)') &
907 'Bogus #=', num_bogus, ' > max_bogus_input=', max_bogus_input
908 call da_error(__FILE__,__LINE__,message(1:1))
909 end if
910
911 ob%bogus(num_bogus)%info = platform%info
912 ob%bogus(num_bogus)%loc = platform%loc
913
914 allocate (ob%bogus(num_bogus)%h (1:nlevels))
915 allocate (ob%bogus(num_bogus)%p (1:nlevels))
916 allocate (ob%bogus(num_bogus)%zk(1:nlevels))
917 allocate (ob%bogus(num_bogus)%u (1:nlevels))
918 allocate (ob%bogus(num_bogus)%v (1:nlevels))
919 allocate (ob%bogus(num_bogus)%t (1:nlevels))
920 allocate (ob%bogus(num_bogus)%q (1:nlevels))
921
922 do i = 1, nlevels
923 ob%bogus(num_bogus)%h(i) = platform%each(i)%height
924 ob%bogus(num_bogus)%p(i) = platform%each(i)%p%inv
925 ob%bogus(num_bogus)%u(i) = platform%each(i)%u
926 ob%bogus(num_bogus)%v(i) = platform%each(i)%v
927 ob%bogus(num_bogus)%t(i) = platform%each(i)%t
928 ob%bogus(num_bogus)%q(i) = platform%each(i)%q
929 end do
930
931 ob%bogus(num_bogus)%slp = platform%loc%slp
932
933 if (print_detail_obs) then
934 write(unit=stdout,fmt=*) 'nlevels=', nlevels
935 write(unit=stdout,fmt=*) 'ob % num_bogus,slp', num_bogus, &
936 ob % bogus (num_bogus) % slp
937 do i=1,nlevels
938 write(unit=stdout,fmt=*) 'ob % num_bogus, i ', num_bogus,i
939 write(unit=stdout,fmt=*) 'ob%bogus(ob%num_bogus)%u,v=', &
940 ob%bogus(num_bogus)%u(i),ob%bogus(num_bogus)%v(i)
941 write(unit=stdout,fmt=*) 'ob%bogus(ob%num_bogus)%t,q=', &
942 ob%bogus(num_bogus)%t(i),ob%bogus(num_bogus)%q(i)
943 end do
944 write(unit=stdout,fmt='(2(a,i4))') 'ob%num_bogus=',num_bogus, &
945 'nlevels=',nlevels
946 end if
947
948 case (18,19) ; ! buoy
949
950 if (n==1) ob%num_buoy_glo = ob%num_buoy_glo + 1
951 if (outside) then
952 cycle reports
953 end if
954
955 if (.not.use_BuoyObs) cycle reports
956
957 num_buoy = num_buoy + 1
958 ! Track serial obs index for reassembly of obs during bit-for-bit
959 ! tests with different numbers of MPI tasks.
960 platform%loc%obs_global_index = ob%num_buoy_glo
961
962 ob%buoy(num_buoy)%info = platform%info
963 ob%buoy(num_buoy)%loc = platform%loc
964
965 ob%buoy(num_buoy)%h = platform%each(1)%height
966 ob%buoy(num_buoy)%u = platform%each(1)%u
967 ob%buoy(num_buoy)%v = platform%each(1)%v
968 ob%buoy(num_buoy)%t = platform%each(1)%t
969 ob%buoy(num_buoy)%p = platform%each(1)%p
970 ob%buoy(num_buoy)%q = platform%each(1)%q
971
972 case (133) ; ! AIRS retrievals
973
974 if (n==1) ob%num_airsr_glo = ob%num_airsr_glo + 1
975 if (outside) then
976 cycle reports
977 end if
978
979 if (.not.use_AIRSRETObs) cycle reports
980
981 num_airsr = num_airsr + 1
982
983 platform%loc%obs_global_index = ob%num_airsr_glo
984
985 ob%airsr(num_airsr)%info = platform%info
986 ob%airsr(num_airsr)%loc = platform%loc
987
988
989 allocate (ob%airsr(num_airsr)%zk(1:nlevels))
990 allocate (ob%airsr(num_airsr)%h (1:nlevels))
991 allocate (ob%airsr(num_airsr)%p (1:nlevels))
992 allocate (ob%airsr(num_airsr)%t (1:nlevels))
993 allocate (ob%airsr(num_airsr)%q (1:nlevels))
994 do i = 1, nlevels
995 ob%airsr(num_airsr)%h(i) = platform%each(i)%height
996 ob%airsr(num_airsr)%p(i) = platform%each(i)%p%inv
997 ob%airsr(num_airsr)%t(i) = platform%each(i)%t
998 ob%airsr(num_airsr)%q(i) = platform%each(i)%q
999 end do
1000
1001 case default;
1002
1003 write(unit=message(1), fmt='(a)') 'unsaved obs found:'
1004 write(unit=message(2), fmt='(2a)') &
1005 'platform%info%platform=', platform%info%platform
1006 write(unit=message(3), fmt='(a, i3)') &
1007 'platform%info%levels=', platform%info%levels
1008 call da_warning(__FILE__,__LINE__,message(1:3))
1009
1010 end select
1011
1012 if (global .and. n < 2) then
1013 if (Testing_WRFVAR) exit dup_loop
1014 if (platform%loc % i >= xp % ide) then
1015 platform%loc%i = platform%loc % i - xp % ide
1016 else if (platform%loc % i < xp % ids) then
1017 platform%loc%i = platform%loc % i + xp % ide
1018 end if
1019
1020 platform%loc%proc_domain = .not. platform%loc%proc_domain
1021 end if
1022 end do dup_loop
1023 end do reports
1024
1025 ob%num_sound = num_sound
1026 ob%num_synop = num_synop
1027 ob%num_pilot = num_pilot
1028 ob%num_satem = num_satem
1029 ob%num_geoamv = num_geoamv
1030 ob%num_polaramv = num_polaramv
1031 ob%num_airep = num_airep
1032 ob%num_gpspw = num_gpspw
1033 ob%num_gpsref = num_gpsref
1034 ob%num_metar = num_metar
1035 ob%num_ships = num_ships
1036 ob%num_qscat = num_qscat
1037 ob%num_buoy = num_buoy
1038 ob%num_profiler = num_profiler
1039 ob%num_bogus = num_bogus
1040 ob%num_airsr = num_airsr
1041
1042 ! WHY
1043 ! ob%num_ssmt1 = num_ssmt1
1044 ! ob%num_ssmt2 = num_ssmt2
1045 ! ob%num_ssmi_tb = num_ssmi_tb
1046 ! ob%num_ssmi_retrieval = num_ssmi_retrieval
1047
1048 if (print_detail_obs) then
1049 write(unit=stdout, fmt='(/a,i6/)') 'ob%current_ob_time=', ob%current_ob_time
1050
1051 write(unit=stdout, fmt='(5x,a,i6,a)') &
1052 'Read: ', num_sound, ' SOUND reports,', &
1053 'Read: ', num_synop, ' SYNOP reports,', &
1054 'Read: ', num_pilot, ' PILOT reports,', &
1055 'Read: ', num_satem, ' SATEM reports,', &
1056 'Read: ', num_geoamv, ' Geo AMV reports,', &
1057 'Read: ', num_polaramv, ' Polar AMV reports,', &
1058 'Read: ', num_airep, ' AIREP reports,', &
1059 'Read: ', num_gpspw , ' GPSPW/GPSZD reports,', &
1060 'Read: ', num_gpsref, ' GPSRF reports,', &
1061 'Read: ', num_metar, ' METAR reports,', &
1062 'Read: ', num_ships , ' SHIP reports,', &
1063 'Read: ', num_ssmi_retrieval , ' SSMI_RETRIEVAL reports,', &
1064 'Read: ', num_ssmi_tb, ' SSMI_TB reports,', &
1065 'Read: ', num_ssmt1, ' SSMT1 reports,', &
1066 'Read: ', num_ssmt2, ' SSMT2 reports,', &
1067 'Read: ', num_qscat, ' QSCAT reports,', &
1068 'Read: ', num_profiler, ' Profiler reports,', &
1069 'Read: ', num_buoy, ' Buoy reports,', &
1070 'Read: ', num_bogus, ' Bogus reports,', &
1071 'Read: ', num_airsr, ' AIRSRET reports,', &
1072 'Read: ', total_obs, ' Total reports.', &
1073 'There are ', total_obs - num_sound - num_synop &
1074 - num_pilot - num_satem &
1075 - num_geoamv - num_polaramv - num_airep &
1076 - num_metar - num_ships &
1077 - num_ssmi_retrieval &
1078 - num_ssmi_tb - num_ssmt1 - num_ssmt2 &
1079 - num_gpspw - num_gpsref - num_qscat &
1080 - num_profiler - num_buoy &
1081 - num_bogus - num_airsr , &
1082 ' reports unsaved.'
1083 write(unit=stdout,fmt=*) ' '
1084 end if
1085
1086 if ((ob%ob_numb(ob%current_ob_time)%sound /= ob%num_sound) .or. &
1087 (ob%ob_numb(ob%current_ob_time)%synop /= ob%num_synop) .or. &
1088 (ob%ob_numb(ob%current_ob_time)%pilot /= ob%num_pilot) .or. &
1089 (ob%ob_numb(ob%current_ob_time)%satem /= ob%num_satem) .or. &
1090 (ob%ob_numb(ob%current_ob_time)%geoamv /= ob%num_geoamv) .or. &
1091 (ob%ob_numb(ob%current_ob_time)%polaramv /= ob%num_polaramv) .or. &
1092 (ob%ob_numb(ob%current_ob_time)%airep /= ob%num_airep) .or. &
1093 (ob%ob_numb(ob%current_ob_time)%gpspw /= ob%num_gpspw) .or. &
1094 (ob%ob_numb(ob%current_ob_time)%gpsref /= ob%num_gpsref) .or. &
1095 (ob%ob_numb(ob%current_ob_time)%metar /= ob%num_metar) .or. &
1096 (ob%ob_numb(ob%current_ob_time)%ships /= ob%num_ships) .or. &
1097 (ob%ob_numb(ob%current_ob_time)%qscat /= ob%num_qscat) .or. &
1098 (ob%ob_numb(ob%current_ob_time)%buoy /= ob%num_buoy) .or. &
1099 (ob%ob_numb(ob%current_ob_time)%bogus /= ob%num_bogus) .or. &
1100 (ob%ob_numb(ob%current_ob_time)%ssmt2 /= ob%num_ssmt1) .or. &
1101 (ob%ob_numb(ob%current_ob_time)%ssmt2 /= ob%num_ssmt2) .or. &
1102 (ob%ob_numb(ob%current_ob_time)%airsr /= ob%num_airsr) .or. &
1103 (ob%ob_numb(ob%current_ob_time)%profiler /= ob%num_profiler)) then
1104
1105 write(unit=stderr, fmt='(a, i6, 2x, a, i6)') &
1106 'ob%ob_numb(ob%current_ob_time)%sound=', &
1107 ob%ob_numb(ob%current_ob_time)%sound, &
1108 'ob%num_sound=', ob%num_sound, &
1109 'ob%ob_numb(ob%current_ob_time)%synop=', &
1110 ob%ob_numb(ob%current_ob_time)%synop, &
1111 'ob%num_synop=', ob%num_synop, &
1112 'ob%ob_numb(ob%current_ob_time)%pilot=', &
1113 ob%ob_numb(ob%current_ob_time)%pilot, &
1114 'ob%num_pilot=', ob%num_pilot, &
1115 'ob%ob_numb(ob%current_ob_time)%satem=', &
1116 ob%ob_numb(ob%current_ob_time)%satem, &
1117 'ob%num_satem=', ob%num_satem, &
1118 'ob%ob_numb(ob%current_ob_time)%geoamv=', &
1119 ob%ob_numb(ob%current_ob_time)%geoamv, &
1120 'ob%num_geoamv=', ob%num_geoamv, &
1121 'ob%ob_numb(ob%current_ob_time)%polaramv=', &
1122 ob%ob_numb(ob%current_ob_time)%polaramv, &
1123 'ob%num_polaramv=', ob%num_polaramv, &
1124 'ob%ob_numb(ob%current_ob_time)%airep=', &
1125 ob%ob_numb(ob%current_ob_time)%airep, &
1126 'ob%num_airep=', ob%num_airep, &
1127 'ob%ob_numb(ob%current_ob_time)%gpspw=', &
1128 ob%ob_numb(ob%current_ob_time)%gpspw, &
1129 'ob%num_gpspw=', ob%num_gpspw, &
1130 'ob%ob_numb(ob%current_ob_time)%gpsref=', &
1131 ob%ob_numb(ob%current_ob_time)%gpsref,&
1132 'ob%num_gpsref=', ob%num_gpsref, &
1133 'ob%ob_numb(ob%current_ob_time)%metar=', &
1134 ob%ob_numb(ob%current_ob_time)%metar, &
1135 'ob%num_metar=', ob%num_metar, &
1136 'ob%ob_numb(ob%current_ob_time)%ships=', &
1137 ob%ob_numb(ob%current_ob_time)%ships, &
1138 'ob%num_ships=', ob%num_ships, &
1139 'ob%ob_numb(ob%current_ob_time)%qscat=', &
1140 ob%ob_numb(ob%current_ob_time)%qscat, &
1141 'ob%num_qscat=', ob%num_qscat, &
1142 'ob%ob_numb(ob%current_ob_time)%buoy =', &
1143 ob%ob_numb(ob%current_ob_time)%buoy , &
1144 'ob%num_buoy =', ob%num_buoy , &
1145 'ob%ob_numb(ob%current_ob_time)%bogus=', &
1146 ob%ob_numb(ob%current_ob_time)%bogus, &
1147 'ob%num_bogus=', ob%num_bogus, &
1148 'ob%ob_numb(ob%current_ob_time)%ssmt1=', &
1149 ob%ob_numb(ob%current_ob_time)%ssmt1, &
1150 'ob%num_ssmt1=', ob%num_ssmt1, &
1151 'ob%ob_numb(ob%current_ob_time)%ssmt2=', &
1152 ob%ob_numb(ob%current_ob_time)%ssmt2, &
1153 'ob%num_ssmt2=', ob%num_ssmt2, &
1154 'ob%ob_numb(ob%current_ob_time)%airsr=', &
1155 ob%ob_numb(ob%current_ob_time)%airsr, &
1156 'ob%num_airsr=', ob%num_airsr, &
1157 'ob%ob_numb(ob%current_ob_time)%profiler=', &
1158 ob%ob_numb(ob%current_ob_time)%profiler, &
1159 'ob%num_profiler=', ob%num_profiler
1160
1161 call da_error(__FILE__,__LINE__,(/"Obs mismatch"/))
1162 end if
1163
1164 ! WHY
1165 ! if ((ob%ob_numb(ob%current_ob_time)%ssmi_tb /= ob%num_ssmi_tb) .or. &
1166 ! (ob%ob_numb(ob%current_ob_time)%ssmi_retrieval /= ob%num_ssmi_retrieval)) then
1167 ! write(unit=message(1), fmt='(a, i6, 2x, a, i6)') &
1168 ! 'ob%ob_numb(ob%current_ob_time)%ssmi_tb =', ob%ob_numb(ob%current_ob_time)%ssmi_tb, &
1169 ! 'ob%num_ssmi_tb =', ob%num_ssmi_tb, &
1170 ! 'ob%ob_numb(ob%current_ob_time)%ssmi_retrieval=', ob%ob_numb(ob%current_ob_time)%ssmi_retrieval, &
1171 ! 'ob%num_ssmi_retrieval=', ob%num_ssmi_retrieval
1172 ! call da_error(__FILE__,__LINE,message(1:1))
1173 ! end if
1174
1175 close(iunit)
1176 call da_free_unit(iunit)
1177
1178 if (trace_use) call da_trace_exit("da_read_obs")
1179
1180 end subroutine da_read_obs
1181
1182