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