da_tune_obs_desroziers.f90
References to this file elsewhere.
1 program da_tune_obs_desroziers
2 !---------------------------------------------------------------------
3 ! Abstract:
4 ! Purpose: Observation error tuning (Desroziers method)
5 ! Ref: QJRMS (2001), 127, pp. 1433-1452
6 ! Gerald Desroziers and Serguei Ivanov
7 !---------------------------------------------------------------------
8
9 implicit none
10
11 integer,parameter :: filename_len = 200
12 integer, parameter :: rand_unit = 45
13 integer, parameter :: yp_unit = 46
14 integer, parameter :: y_unit = 47
15 integer, parameter :: jo_unit = 48
16 integer, parameter :: in_unit = 49
17 integer, parameter :: ninst = 35 ! max sensor number
18 integer, Parameter :: nplatforms = 20
19 real, parameter :: missing_r = -888888.0
20 ! below copy from RTTOV
21 ! platform names
22 Character (len=8), Parameter :: platform_name(nplatforms) = &
23 (/ 'noaa ', 'dmsp ', 'meteosat', 'goes ', 'gms ', &
24 'fy2 ', 'trmm ', 'ers ', 'eos ', 'metop ', &
25 'envisat ', 'msg ', 'fy1 ', 'adeos ', 'mtsat ', &
26 'coriolis', 'npoess ', 'gifts ', 'xxxxxxxx', 'xxxxxxxx'/)
27
28 ! List of instruments !!!! HIRS is number 0
29 Character (len=8), Dimension(0:ninst-1) :: inst_name = &
30 & (/ 'hirs ', 'msu ', 'ssu ', 'amsua ', 'amsub ', &
31 & 'avhrr ', 'ssmi ', 'vtpr1 ', 'vtpr2 ', 'tmi ', &
32 & 'ssmis ', 'airs ', 'hsb ', 'modis ', 'atsr ', &
33 & 'mhs ', 'iasi ', 'amsr ', 'imager ', 'atms ', &
34 & 'mviri ', 'seviri ', 'imager ', 'sounder ', 'imager ', &
35 & 'vissr ', 'mvisr ', 'cris ', 'cmis ', 'viirs ', &
36 & 'windsat ', 'gifts ', 'xxxxxxxx', 'xxxxxxxx', 'xxxxxxxx' /)
37
38 integer :: n,ipixel
39
40 ! radiance namelist variables
41
42 integer :: rtminit_nsensor
43 integer, dimension(ninst) :: rtminit_platform, &
44 rtminit_satid, &
45 rtminit_sensor, &
46 rtminit_nchan
47 ! radiance relevant variables
48 character*1 :: str1,str2,str3
49 character*4 :: platform
50 character*5 :: sensor
51 integer :: platform_id,satellite_id,sensor_id,ichan
52
53 type field_type
54 real :: yp
55 real :: y
56 real :: error
57 real :: pert
58 end type field_type
59
60 type surfc_type
61 type (field_type) :: u
62 type (field_type) :: v
63 type (field_type) :: t
64 type (field_type) :: p
65 type (field_type) :: q
66 end type surfc_type
67
68 type geoamv_type
69 integer :: numlevs
70 type (field_type),pointer :: u(:)
71 type (field_type),pointer :: v(:)
72 end type geoamv_type
73
74 type polaramv_type
75 integer :: numlevs
76 type (field_type),pointer :: u(:)
77 type (field_type),pointer :: v(:)
78 end type polaramv_type
79
80 type gpspw_type
81 type (field_type) :: tpw
82 end type gpspw_type
83
84 type sound_type
85 integer :: numlevs
86 type (field_type), pointer :: u(:)
87 type (field_type), pointer :: v(:)
88 type (field_type), pointer :: t(:)
89 type (field_type), pointer :: q(:)
90 end type sound_type
91
92 type airsr_type
93 integer :: numlevs
94 type (field_type), pointer :: t(:)
95 type (field_type), pointer :: q(:)
96 end type airsr_type
97
98 type airep_type
99 integer :: numlevs
100 type (field_type), pointer :: u(:)
101 type (field_type), pointer :: v(:)
102 type (field_type), pointer :: t(:)
103 end type airep_type
104
105 type pilot_type
106 integer :: numlevs
107 type (field_type), pointer :: u(:)
108 type (field_type), pointer :: v(:)
109 end type pilot_type
110
111 type ssmir_type
112 type (field_type) :: speed
113 type (field_type) :: tpw
114 end type ssmir_type
115
116 type satem_type
117 integer :: numlevs
118 type (field_type), pointer :: thickness(:)
119 end type satem_type
120
121 type ssmt1_type
122 integer :: numlevs
123 type (field_type), pointer :: t(:)
124 end type ssmt1_type
125
126 type ssmt2_type
127 integer :: numlevs
128 type (field_type), pointer :: rh(:)
129 end type ssmt2_type
130
131 type qscat_type
132 type (field_type) :: u
133 type (field_type) :: v
134 end type qscat_type
135
136 type bogus_type
137 integer :: numlevs
138 type (field_type) :: slp
139 type (field_type), pointer :: u(:)
140 type (field_type), pointer :: v(:)
141 type (field_type), pointer :: t(:)
142 type (field_type), pointer :: q(:)
143 end type bogus_type
144
145 type pixel_type
146 type (field_type), pointer :: pixel(:) ! dimension: num_rad_tot(ichan)
147 end type pixel_type
148
149 type radiance_type
150 character*20 :: rttovid_string
151 integer :: nchan
152 integer, pointer :: num_rad_tot(:) ! dimension: nchan
153 real , pointer :: trace_rad (:), & ! dimension: nchan
154 jo_rad (:), & ! dimension: nchan
155 joa_rad (:), & ! dimension: nchan
156 factor_rad (:) ! dimension: nchan
157 type (pixel_type), pointer :: tb(:) ! dimension: nchan
158 end type radiance_type
159
160 type gpsref_type
161 integer :: numlevs
162 type (field_type), pointer :: ref(:)
163 end type gpsref_type
164
165
166 type iv_type
167 integer :: total_obs
168 integer :: num_synop, num_metar, num_ships, &
169 num_polaramv, num_geoamv,num_gpspw, num_sound, &
170 num_airep, num_pilot, num_ssmir, num_airsr, &
171 num_satem, num_ssmt1, num_ssmt2, num_gpsref, &
172 num_buoy, num_qscat, num_sonde_sfc, num_bogus, num_profiler
173 integer :: num_synop_tot, num_metar_tot, num_ships_tot, &
174 num_polaramv_tot, num_geoamv_tot, num_gpspw_tot, &
175 num_sound_tot, num_airsr_tot, &
176 num_airep_tot, num_pilot_tot, num_ssmir_tot, &
177 num_satem_tot, num_ssmt1_tot, num_ssmt2_tot, &
178 num_buoy_tot, num_qscat_tot, num_sonde_sfc_tot, num_bogus_tot, &
179 num_profiler_tot, num_gpsref_tot
180
181 real :: trace_total
182 real :: trace_synop, trace_metar, trace_ships, &
183 trace_polaramv, trace_geoamv, trace_gpspw, trace_sound, &
184 trace_airep, trace_pilot, trace_ssmir, trace_airsr, &
185 trace_satem, trace_ssmt1, trace_ssmt2, &
186 trace_buoy, trace_qscat, trace_sonde_sfc, trace_bogus, &
187 trace_profiler, trace_gpsref
188 real :: jo_synop_u, jo_synop_v, jo_synop_t, jo_synop_p, jo_synop_q
189 real :: jo_metar_u, jo_metar_v, jo_metar_t, jo_metar_p, jo_metar_q
190 real :: jo_ships_u, jo_ships_v, jo_ships_t, jo_ships_p, jo_ships_q
191 real :: jo_polaramv_u, jo_polaramv_v, jo_geoamv_u, jo_geoamv_v, jo_gpspw_tpw
192 real :: jo_sound_u, jo_sound_v, jo_sound_t, jo_sound_q
193 real :: jo_gpsref_ref
194 real :: jo_airsr_t, jo_airsr_q
195 real :: jo_airep_u, jo_airep_v, jo_airep_t
196 real :: jo_pilot_u, jo_pilot_v
197 real :: jo_ssmir_speed, jo_ssmir_tpw, jo_satem_thickness
198 real :: jo_ssmt1_t, jo_ssmt2_rh
199 real :: jo_buoy_u, jo_buoy_v, jo_buoy_t, jo_buoy_p, jo_buoy_q
200 real :: jo_bogus_u, jo_bogus_v, jo_bogus_t, jo_bogus_q, jo_bogus_slp
201 real :: jo_qscat_u, jo_qscat_v
202 real :: jo_sonde_sfc_u, jo_sonde_sfc_v, jo_sonde_sfc_t, jo_sonde_sfc_p, jo_sonde_sfc_q
203 real :: jo_profiler_u, jo_profiler_v
204
205 real :: joa_synop_u, joa_synop_v, joa_synop_t, joa_synop_p, joa_synop_q
206 real :: joa_metar_u, joa_metar_v, joa_metar_t, joa_metar_p, joa_metar_q
207 real :: joa_ships_u, joa_ships_v, joa_ships_t, joa_ships_p, joa_ships_q
208 real :: joa_polaramv_u, joa_polaramv_v, joa_geoamv_u, joa_geoamv_v
209 real :: joa_gpspw_tpw
210 real :: joa_sound_u, joa_sound_v, joa_sound_t, joa_sound_q
211 real :: joa_gpsref_ref
212 real :: joa_airsr_t, joa_airsr_q
213 real :: joa_airep_u, joa_airep_v, joa_airep_t
214 real :: joa_pilot_u, joa_pilot_v
215 real :: joa_ssmir_speed, joa_ssmir_tpw, joa_satem_thickness
216 real :: joa_ssmt1_t, joa_ssmt2_rh
217 real :: joa_buoy_u, joa_buoy_v, joa_buoy_t, joa_buoy_p, joa_buoy_q
218 real :: joa_bogus_u, joa_bogus_v, joa_bogus_t, joa_bogus_q, joa_bogus_slp
219 real :: joa_qscat_u, joa_qscat_v
220 real :: joa_sonde_sfc_u, joa_sonde_sfc_v, joa_sonde_sfc_t
221 real :: joa_sonde_sfc_p, joa_sonde_sfc_q
222 real :: joa_profiler_u, joa_profiler_v
223
224 real :: ef_synop_u, ef_synop_v, ef_synop_t, ef_synop_p, ef_synop_q
225 real :: ef_metar_u, ef_metar_v, ef_metar_t, ef_metar_p, ef_metar_q
226 real :: ef_ships_u, ef_ships_v, ef_ships_t, ef_ships_p, ef_ships_q
227 real :: ef_polaramv_u, ef_polaramv_v, ef_geoamv_u, ef_geoamv_v
228 real :: ef_gpspw_tpw
229 real :: ef_sound_u, ef_sound_v, ef_sound_t, ef_sound_q
230 real :: ef_gpsref_ref
231 real :: ef_airsr_t, ef_airsr_q
232 real :: ef_airep_u, ef_airep_v, ef_airep_t
233 real :: ef_pilot_u, ef_pilot_v
234 real :: ef_ssmir_speed, ef_ssmir_tpw, ef_satem_thickness
235 real :: ef_ssmt1_t, ef_ssmt2_rh
236 real :: ef_buoy_u, ef_buoy_v, ef_buoy_t, ef_buoy_p, ef_buoy_q
237 real :: ef_bogus_u, ef_bogus_v, ef_bogus_t, ef_bogus_q, ef_bogus_slp
238 real :: ef_qscat_u, ef_qscat_v
239 real :: ef_sonde_sfc_u, ef_sonde_sfc_v, ef_sonde_sfc_t, ef_sonde_sfc_p
240 real :: ef_sonde_sfc_q
241 real :: ef_profiler_u, ef_profiler_v
242
243 type (surfc_type), pointer :: synop(:)
244 type (surfc_type), pointer :: metar(:)
245 type (surfc_type), pointer :: ships(:)
246 type (polaramv_type), pointer :: polaramv(:)
247 type (geoamv_type), pointer :: geoamv(:)
248 type (gpspw_type), pointer :: gpspw(:)
249 type (sound_type), pointer :: sound(:)
250 type (gpsref_type), pointer :: gpsref(:)
251 type (airsr_type), pointer :: airsr(:)
252 type (airep_type), pointer :: airep(:)
253 type (pilot_type), pointer :: pilot(:)
254 type (ssmir_type), pointer :: ssmir(:)
255 type (satem_type), pointer :: satem(:)
256 type (ssmt1_type), pointer :: ssmt1(:)
257 type (ssmt2_type), pointer :: ssmt2(:)
258 type (surfc_type), pointer :: sonde_sfc(:)
259 type (surfc_type), pointer :: buoy(:)
260 type (qscat_type), pointer :: qscat(:)
261 type (pilot_type), pointer :: profiler(:)
262 type (bogus_type), pointer :: bogus(:)
263 type (radiance_type), pointer :: rad(:)
264
265 end type iv_type
266 type (iv_type) :: ob
267
268 !--------------------------------------------------------------------------
269 ! Initialise the counter
270 !--------------------------------------------------------------------------
271
272 ! ob % totla_obs = 0
273 ! ob % num_synop_tot = 0
274 ! ob % num_metar_tot = 0
275 ! ob % num_ships_tot = 0
276 ! ob % num_polaramv_tot = 0
277 ! ob % num_geoamv_tot = 0
278 ! ob % num_gpspw_tot = 0
279 ! ob % num_sound_tot = 0
280 ! ob % num_airsr_tot = 0
281 ! ob % num_airep_tot = 0
282 ! ob % num_pilot_tot = 0
283 ! ob % num_ssmir_tot = 0
284 ! ob % num_satem_tot = 0
285 ! ob % num_ssmt1_tot = 0
286 ! ob % num_ssmt2_tot = 0
287 ! ob % num_buoy_tot = 0
288 ! ob % num_sonde_sfc_tot = 0
289 ! ob % num_qscat_tot = 0
290 ! ob % num_profiler_tot = 0
291 ! ob% num_bogus_tot = 0
292
293 !--------------------------------------------------------------------------
294 ! [1.0] Count total number of observations and allocate arrays:
295 !--------------------------------------------------------------------------
296
297 call da_count_obs( y_unit, ob )
298
299 !--------------------------------------------------------------------------
300 ! [2.0] Read in UNperturbed y = H(x_inc) from each ob type:
301 !--------------------------------------------------------------------------
302
303 call da_read_y( y_unit, ob )
304
305 !--------------------------------------------------------------------------
306 ! [3.0] Read in perturbed yp = H(x_inc_p) from each ob type:
307 !--------------------------------------------------------------------------
308
309 call da_read_yp( yp_unit, ob )
310
311 !--------------------------------------------------------------------------
312 ! [4.0] Read in perturbations and errors from each ob type:
313 !--------------------------------------------------------------------------
314
315 call da_read_obs_rand( rand_unit, ob )
316
317 !--------------------------------------------------------------------------
318 ! [5.0] Calculate expected cost function J values:
319 !--------------------------------------------------------------------------
320
321 call da_calc_jo_expected( ob )
322
323 !--------------------------------------------------------------------------
324 ! [6.0] Read actual cost function J and error tuning factors used:
325 !--------------------------------------------------------------------------
326
327 call da_read_jo_actual( ob )
328
329 !--------------------------------------------------------------------------
330 ! [7.0] Calculate observation error tuning factors:
331 !--------------------------------------------------------------------------
332
333 call da_calc_new_factors( ob )
334
335 !--------------------------------------------------------------------------
336 ! [8.0] Calculate observation error tuning factors:
337 !--------------------------------------------------------------------------
338
339 call da_get_j( ob )
340
341 contains
342
343 subroutine da_count_obs( y_unit, ob )
344
345 implicit none
346
347 integer, intent(in) :: y_unit
348 type (iv_type), intent(inout) :: ob
349
350 character*20 :: ob_name, dummy
351 integer :: times, num_obs, k, num_levs
352
353 ! [1] Initialize ob numbers:
354
355 ob % num_synop = 0
356 ob % num_metar = 0
357 ob % num_ships = 0
358 ob % num_polaramv = 0
359 ob % num_geoamv = 0
360 ob % num_gpspw = 0
361 ob % num_sound = 0
362 ob % num_airsr = 0
363 ob % num_airep = 0
364 ob % num_pilot = 0
365 ob % num_ssmir = 0
366 ob % num_satem = 0
367 ob % num_ssmt1 = 0
368 ob % num_ssmt2 = 0
369 ob % num_sonde_sfc = 0
370 ob % num_buoy = 0
371 ob % num_profiler = 0
372 ob % num_bogus = 0
373 ob % num_qscat = 0
374 ob % num_gpsref = 0
375 ! [1.2] Initialize satellite instrument numbers
376 ! and channel numbers for each instrument
377
378 call read_namelist_radiance
379
380 times = 0
381 ! [2] Loop through input file to count number of obs:
382
383 do
384 read(y_unit,'(a20,i8)', end = 1000)ob_name, num_obs
385 if( num_obs > 0) times = times + 1
386 if ( index( ob_name,'synop') > 0 ) then
387 ob % num_synop = ob % num_synop + num_obs
388 elseif ( index( ob_name,'metar') > 0 ) then
389 ob % num_metar = ob % num_metar + num_obs
390 elseif ( index( ob_name,'ships') > 0 ) then
391 ob % num_ships = ob % num_ships + num_obs
392 elseif ( index( ob_name,'polaramv') > 0 ) then
393 ob % num_polaramv = ob % num_polaramv + num_obs
394 elseif ( index( ob_name,'geoamv') > 0 ) then
395 ob % num_geoamv = ob % num_geoamv + num_obs
396 elseif ( index( ob_name,'gpspw') > 0 ) then
397 ob % num_gpspw = ob % num_gpspw + num_obs
398 elseif ( index( ob_name,'sound') > 0 ) then
399 ob % num_sound = ob % num_sound + num_obs
400 elseif ( index( ob_name,'airsr') > 0 ) then
401 ob % num_airsr = ob % num_airsr + num_obs
402 elseif ( index( ob_name,'airep') > 0 ) then
403 ob % num_airep = ob % num_airep + num_obs
404 elseif ( index( ob_name,'pilot') > 0 ) then
405 ob % num_pilot = ob % num_pilot + num_obs
406 elseif ( index( ob_name,'ssmir') > 0 ) then
407 ob % num_ssmir = ob % num_ssmir + num_obs
408 elseif ( index( ob_name,'satem') > 0 ) then
409 ob % num_satem = ob % num_satem + num_obs
410 elseif ( index( ob_name,'ssmt1') > 0 ) then
411 ob % num_ssmt1 = ob % num_ssmt1 + num_obs
412 elseif ( index( ob_name,'ssmt2') > 0 ) then
413 ob % num_ssmt2 = ob % num_ssmt2 + num_obs
414 else if ( index( ob_name,'qscat') > 0 ) then
415 ob % num_qscat = ob % num_qscat + num_obs
416 else if ( index( ob_name,'sonde_sfc') > 0 ) then
417 ob % num_sonde_sfc = ob % num_sonde_sfc + num_obs
418 else if ( index( ob_name,'buoy') > 0 ) then
419 ob % num_buoy = ob % num_buoy + num_obs
420 else if ( index( ob_name,'profiler') > 0 ) then
421 ob % num_profiler = ob % num_profiler + num_obs
422 else if ( index( ob_name,'bogus') > 0 ) then
423 ob % num_bogus = ob % num_bogus + num_obs
424 elseif ( index( ob_name,'gpsref') > 0 ) then
425 ob % num_gpsref = ob % num_gpsref + num_obs
426 ! Radiance obs: consistent with RTTOV triplet and WRF-VAR
427 !--------------------------------------------------------------------
428 else if ( index( ob_name,'noaa') > 0 ) then
429 platform_id = 1
430 read (ob_name,'(a4,a1,i2,a1,a5,a1,i4)') &
431 platform, str1,satellite_id,str2,sensor,str3,ichan
432 if ( sensor == 'amsua' ) then
433 sensor_id = 3
434 else if ( sensor == 'amsub' ) then
435 sensor_id = 4
436 else
437 print*,' Unrecognized Sensor ',sensor
438 stop
439 end if
440 do n = 1, rtminit_nsensor
441 if ( platform_id == rtminit_platform(n) &
442 .and. satellite_id == rtminit_satid(n) &
443 .and. sensor_id == rtminit_sensor(n) ) then
444 ob%rad(n)%num_rad_tot(ichan) = ob%rad(n)%num_rad_tot(ichan) + &
445 num_obs
446 exit
447 end if
448 end do
449
450 ! WHY
451 ! read(y_unit,'(i8)')num_levs
452 ! read(y_unit,'(a20)')dummy
453 elseif ( index( ob_name,'*****') > 0 ) then
454 print*,' scaned for time ', times
455 exit
456 else
457 print*,' unknown obs type: ',trim(ob_name),' found on unit ',y_unit
458 end if
459
460 if ( index( ob_name,'noaa') > 0 ) then
461 do k = 1, num_obs
462 read(y_unit,'(a20)')dummy
463 end do
464 else
465 do n = 1, num_obs
466 read(y_unit,'(i8)')num_levs
467 do k = 1, num_levs
468 read(y_unit,'(a20)')dummy
469 end do
470 end do
471 end if
472 end do
473
474 1000 print*,' end of file reached on unit ',y_unit
475
476 ! [3] Allocate ob structures where obs exist:
477
478 if ( rtminit_nsensor > 0 ) then
479 do n=1,rtminit_nsensor
480 do ichan=1,ob%rad(n)%nchan
481 if ( ob%rad(n)%num_rad_tot(ichan) > 0 ) then
482 allocate(ob%rad(n)%tb(ichan)%pixel &
483 (1:ob % rad(n)%num_rad_tot(ichan)))
484 write(6,'(a,i5,a,i10)') &
485 ' Number of '//trim(ob%rad(n)%rttovid_string)//&
486 ' channel ', ichan , ' = ', &
487 ob%rad(n)%num_rad_tot(ichan)
488 end if
489 end do
490 end do
491 end if
492
493 if ( ob % num_synop > 0 ) then
494 allocate( ob % synop(1:ob % num_synop) )
495 write(6,'(a,i8)')' Number of synop obs = ', ob % num_synop
496 end if
497
498 if ( ob % num_metar > 0 ) then
499 allocate( ob % metar(1:ob % num_metar) )
500 write(6,'(a,i8)')' Number of metar obs = ', ob % num_metar
501 end if
502
503 if ( ob % num_ships > 0 ) then
504 allocate( ob % ships(1:ob % num_ships) )
505 write(6,'(a,i8)')' Number of ships obs = ', ob % num_ships
506 end if
507
508 if ( ob % num_polaramv > 0 ) then
509 allocate( ob % polaramv(1:ob % num_polaramv) )
510 write(6,'(a,i8)')' Number of Polar AMV obs = ', ob % num_polaramv
511 end if
512
513 if ( ob % num_geoamv > 0 ) then
514 allocate( ob % geoamv(1:ob % num_geoamv) )
515 write(6,'(a,i8)')' Number of Geo AMV obs = ', ob % num_geoamv
516 end if
517
518 if ( ob % num_gpspw > 0 ) then
519 allocate( ob % gpspw(1:ob % num_gpspw) )
520 write(6,'(a,i8)')' Number of gpspw obs = ', ob % num_gpspw
521 end if
522
523 if ( ob % num_sound > 0 ) then
524 allocate( ob % sound(1:ob % num_sound) )
525 write(6,'(a,i8)')' Number of sound obs = ', ob % num_sound
526 end if
527 if ( ob % num_airsr > 0 ) then
528 allocate( ob % airsr(1:ob % num_airsr) )
529 write(6,'(a,i8)')' Number of AIRS retrievals = ', ob % num_airsr
530 end if
531
532 if ( ob % num_airep > 0 ) then
533 allocate( ob % airep(1:ob % num_airep) )
534 write(6,'(a,i8)')' Number of airep obs = ', ob % num_airep
535 end if
536
537 if ( ob % num_pilot > 0 ) then
538 allocate( ob % pilot(1:ob % num_pilot) )
539 write(6,'(a,i8)')' Number of pilot obs = ', ob % num_pilot
540 end if
541
542 if ( ob % num_ssmir > 0 ) then
543 allocate( ob % ssmir(1:ob % num_ssmir) )
544 write(6,'(a,i8)')' Number of ssmir obs = ', ob % num_ssmir
545 end if
546
547 if ( ob % num_satem > 0 ) then
548 allocate( ob % satem(1:ob % num_satem) )
549 write(6,'(a,i8)')' Number of satem obs = ', ob % num_satem
550 end if
551
552 if ( ob % num_ssmt1 > 0 ) then
553 allocate( ob % ssmt1(1:ob % num_ssmt1) )
554 write(6,'(a,i8)')' Number of ssmt1 obs = ', ob % num_ssmt1
555 end if
556
557 if ( ob % num_ssmt2 > 0 ) then
558 allocate( ob % ssmt2(1:ob % num_ssmt2) )
559 write(6,'(a,i8)')' Number of ssmt2 obs = ', ob % num_ssmt2
560 end if
561 if ( ob % num_bogus > 0 ) then
562 allocate( ob % bogus(1:ob % num_bogus) )
563 write(6,'(a,i8)')' Number of bogus obs = ', ob % num_bogus
564 end if
565 if ( ob % num_buoy > 0 ) then
566 allocate( ob % buoy(1:ob % num_buoy) )
567 write(6,'(a,i8)')' Number of buoy obs = ', ob % num_buoy
568 end if
569 if ( ob % num_sonde_sfc > 0 ) then
570 allocate( ob % sonde_sfc(1:ob % num_sonde_sfc) )
571 write(6,'(a,i8)')' Number of sonde_sfc obs = ', ob % num_sonde_sfc
572 end if
573 if ( ob % num_qscat > 0 ) then
574 allocate( ob % qscat(1:ob % num_qscat) )
575 write(6,'(a,i8)')' Number of qscat obs = ', ob % num_qscat
576 end if
577 if ( ob % num_profiler> 0 ) then
578 allocate( ob % profiler(1:ob % num_profiler) )
579 write(6,'(a,i8)')' Number of Profiler obs = ', ob % num_profiler
580 end if
581 if ( ob % num_gpsref > 0 ) then
582 allocate( ob % gpsref(1:ob % num_gpsref) )
583 write(6,'(a,i8)')' Number of gpsref obs = ', ob % num_gpsref
584 end if
585
586 end subroutine da_count_obs
587
588 subroutine da_read_y( y_unit, ob )
589
590 implicit none
591
592 integer, intent(in) :: y_unit
593 type (iv_type), intent(inout) :: ob
594
595 character*20 :: ob_name
596 integer :: n, ndum, k, kdum, num_obs, num_levs
597 integer :: isynop, imetar, iships, ipolaramv, igeoamv, igpspw, isound, &
598 iairep, ipilot, issmir, isatem, issmt1, issmt2, iairsr, &
599 ibuoy, isonde_sfc, iqscat, ibogus, iprofiler, igpsref
600
601 rewind (y_unit)
602
603 isynop = 0
604 imetar = 0
605 iships = 0
606 ipolaramv = 0
607 igeoamv = 0
608 igpspw = 0
609 isound = 0
610 iairsr = 0
611 iairep = 0
612 ipilot = 0
613 issmir = 0
614 isatem = 0
615 issmt1 = 0
616 issmt2 = 0
617 ibuoy = 0
618 isonde_sfc = 0
619 iqscat = 0
620 ibogus = 0
621 iprofiler = 0
622 igpsref = 0
623
624 do n = 1,rtminit_nsensor
625 ob%rad(n)%num_rad_tot(:) = 0
626 end do
627
628 do
629 read(y_unit,'(a20,i8)', end = 1000)ob_name, num_obs
630
631 if ( index( ob_name,'synop') > 0 ) then
632 do n = 1, num_obs
633 isynop = isynop + 1
634 read(y_unit,'(i8)')num_levs
635 read(y_unit,'(2i8,7e15.7)')ndum, kdum, ob % synop(isynop) % u % y, &
636 ob % synop(isynop) % v % y, &
637 ob % synop(isynop) % t % y, &
638 ob % synop(isynop) % p % y, &
639 ob % synop(isynop) % q % y
640 end do
641 elseif ( index( ob_name,'metar') > 0 ) then
642 do n = 1, num_obs
643 imetar = imetar + 1
644 read(y_unit,'(i8)')num_levs
645 read(y_unit,'(2i8,7e15.7)')ndum, kdum, ob % metar(imetar) % u % y, &
646 ob % metar(imetar) % v % y, &
647 ob % metar(imetar) % t % y, &
648 ob % metar(imetar) % p % y, &
649 ob % metar(imetar) % q % y
650 end do
651 elseif ( index( ob_name,'ships') > 0 ) then
652 do n = 1, num_obs
653 iships = iships + 1
654 read(y_unit,'(i8)')num_levs
655 read(y_unit,'(2i8,7e15.7)')ndum, kdum, ob % ships(iships) % u % y, &
656 ob % ships(iships) % v % y, &
657 ob % ships(iships) % t % y, &
658 ob % ships(iships) % p % y, &
659 ob % ships(iships) % q % y
660 end do
661
662 elseif ( index( ob_name,'geoamv') > 0 ) then
663
664 do n = 1, num_obs
665 igeoamv = igeoamv + 1
666 read(y_unit,'(i8)')num_levs
667 ob % geoamv(igeoamv) % numlevs = num_levs
668 allocate( ob % geoamv(igeoamv) % u(1:num_levs) )
669 allocate( ob % geoamv(igeoamv) % v(1:num_levs) )
670 do k = 1, num_levs
671 read(y_unit,'(2i8,7e15.7)')ndum, kdum, &
672 ob % geoamv(igeoamv) % u(k) % y, &
673 ob % geoamv(igeoamv) % v(k) % y
674 end do
675 end do
676
677 elseif ( index( ob_name,'polaramv') > 0 ) then
678
679 do n = 1, num_obs
680 ipolaramv = ipolaramv + 1
681 read(y_unit,'(i8)')num_levs
682 ob % polaramv(ipolaramv) % numlevs = num_levs
683 allocate( ob % polaramv(ipolaramv) % u(1:num_levs) )
684 allocate( ob % polaramv(ipolaramv) % v(1:num_levs) )
685 do k = 1, num_levs
686 read(y_unit,'(2i8,7e15.7)')ndum, kdum, &
687 ob % polaramv(ipolaramv) % u(k) % y, &
688 ob % polaramv(ipolaramv) % v(k) % y
689 end do
690 end do
691
692 elseif ( index( ob_name,'gpspw') > 0 ) then
693
694 do n = 1, num_obs
695 igpspw = igpspw + 1
696 read(y_unit,'(i8)')num_levs
697 read(y_unit,'(2i8,7e15.7)')ndum, kdum, &
698 ob % gpspw(igpspw) % tpw % y
699 end do
700
701 elseif ( index( ob_name,'sound') > 0 ) then
702 do n = 1, num_obs
703 isound = isound + 1
704 read(y_unit,'(i8)')num_levs
705 ob % sound(isound) % numlevs = num_levs
706 allocate( ob % sound(isound) % u(1:num_levs) )
707 allocate( ob % sound(isound) % v(1:num_levs) )
708 allocate( ob % sound(isound) % t(1:num_levs) )
709 allocate( ob % sound(isound) % q(1:num_levs) )
710
711 do k = 1, num_levs
712 read(y_unit,'(2i8,7e15.7)')ndum, kdum, &
713 ob % sound(isound) % u(k) % y, &
714 ob % sound(isound) % v(k) % y, &
715 ob % sound(isound) % t(k) % y, &
716 ob % sound(isound) % q(k) % y
717 end do
718 end do
719
720 elseif ( index( ob_name,'airsr') > 0 ) then
721 do n = 1, num_obs
722 iairsr = iairsr + 1
723 read(y_unit,'(i8)')num_levs
724 ob % airsr(iairsr) % numlevs = num_levs
725 allocate( ob % airsr(iairsr) % t(1:num_levs) )
726 allocate( ob % airsr(iairsr) % q(1:num_levs) )
727
728 do k = 1, num_levs
729 read(y_unit,'(2i8,7e15.7)')ndum, kdum, &
730 ob % airsr(iairsr) % t(k) % y, &
731 ob % airsr(iairsr) % q(k) % y
732 end do
733 end do
734 elseif ( index( ob_name,'airep') > 0 ) then
735 do n = 1, num_obs
736 iairep = iairep + 1
737 read(y_unit,'(i8)')num_levs
738 ob % airep(iairep) % numlevs = num_levs
739 allocate( ob % airep(iairep) % u(1:num_levs) )
740 allocate( ob % airep(iairep) % v(1:num_levs) )
741 allocate( ob % airep(iairep) % t(1:num_levs) )
742
743 do k = 1, num_levs
744 read(y_unit,'(2i8,7e15.7)')ndum, kdum, &
745 ob % airep(iairep) % u(k) % y, &
746 ob % airep(iairep) % v(k) % y, &
747 ob % airep(iairep) % t(k) % y
748 end do
749 end do
750 elseif ( index( ob_name,'pilot') > 0 ) then
751 do n = 1, num_obs
752 ipilot = ipilot + 1
753 read(y_unit,'(i8)')num_levs
754 ob % pilot(ipilot) % numlevs = num_levs
755 allocate( ob % pilot(ipilot) % u(1:num_levs) )
756 allocate( ob % pilot(ipilot) % v(1:num_levs) )
757
758 do k = 1, num_levs
759 read(y_unit,'(2i8,7e15.7)')ndum, kdum, &
760 ob % pilot(ipilot) % u(k) % y, &
761 ob % pilot(ipilot) % v(k) % y
762 end do
763 end do
764 elseif ( index( ob_name,'ssmir') > 0 ) then
765 do n = 1, num_obs
766 issmir = issmir + 1
767 read(y_unit,'(i8)')num_levs
768 read(y_unit,'(2i8,7e15.7)')ndum, kdum, &
769 ob % ssmir(issmir) % speed % y, &
770 ob % ssmir(issmir) % tpw % y
771 end do
772 elseif ( index( ob_name,'satem') > 0 ) then
773 do n = 1, num_obs
774 isatem = isatem + 1
775 read(y_unit,'(i8)')num_levs
776 ob % satem(isatem) % numlevs = num_levs
777 allocate( ob % satem(isatem) % thickness(1:num_levs) )
778
779 do k = 1, num_levs
780 read(y_unit,'(2i8,7e15.7)')ndum, kdum, &
781 ob % satem(isatem) % thickness(k) % y
782 end do
783 end do
784 elseif ( index( ob_name,'ssmt1') > 0 ) then
785 do n = 1, num_obs
786 issmt1 = issmt1 + 1
787 read(y_unit,'(i8)')num_levs
788 ob % ssmt1(issmt1) % numlevs = num_levs
789 allocate( ob % ssmt1(issmt1) % t(1:num_levs) )
790
791 do k = 1, num_levs
792 read(y_unit,'(2i8,7e15.7)')ndum, kdum, &
793 ob % ssmt1(issmt1) % t(k) % y
794 end do
795 end do
796 elseif ( index( ob_name,'ssmt2') > 0 ) then
797 do n = 1, num_obs
798 issmt2 = issmt2 + 1
799 read(y_unit,'(i8)')num_levs
800 ob % ssmt2(issmt2) % numlevs = num_levs
801 allocate( ob % ssmt2(issmt2) % rh(1:num_levs) )
802
803 do k = 1, num_levs
804 read(y_unit,'(2i8,7e15.7)')ndum, kdum, &
805 ob % ssmt2(issmt2) % rh(k) % y
806 end do
807 end do
808 elseif ( index( ob_name,'buoy') > 0 ) then
809 do n = 1, num_obs
810 ibuoy = ibuoy + 1
811 read(y_unit,'(i8)')num_levs
812 read(y_unit,'(2i8,7e15.7)')ndum, kdum, ob % buoy(ibuoy) % u % y, &
813 ob % buoy(ibuoy) % v % y, &
814 ob % buoy(ibuoy) % t % y, &
815 ob % buoy(ibuoy) % p % y, &
816 ob % buoy(ibuoy) % q % y
817 end do
818 elseif ( index( ob_name,'sonde_sfc') > 0 ) then
819 do n = 1, num_obs
820 isonde_sfc = isonde_sfc + 1
821 read(y_unit,'(i8)')num_levs
822 read(y_unit,'(2i8,7e15.7)')ndum, kdum, ob % sonde_sfc(isonde_sfc) % u % y, &
823 ob % sonde_sfc(isonde_sfc) % v % y, &
824 ob % sonde_sfc(isonde_sfc) % t % y, &
825 ob % sonde_sfc(isonde_sfc) % p % y, &
826 ob % sonde_sfc(isonde_sfc) % q % y
827 end do
828
829 elseif ( index( ob_name,'qscat') > 0 ) then
830 do n = 1, num_obs
831 iqscat = iqscat + 1
832 read(y_unit,'(i8)')num_levs
833 read(y_unit,'(2i8,7e15.7)')ndum, kdum, ob % qscat(iqscat) % u % y, &
834 ob % qscat(iqscat) % v % y
835 end do
836 elseif ( index( ob_name,'profiler') > 0 ) then
837 do n = 1, num_obs
838 iprofiler = iprofiler + 1
839 read(y_unit,'(i8)')num_levs
840 ob % profiler(iprofiler) % numlevs = num_levs
841 allocate( ob % profiler(iprofiler) % u(1:num_levs) )
842 allocate( ob % profiler(iprofiler) % v(1:num_levs) )
843
844 do k = 1, num_levs
845 read(y_unit,'(2i8,7e15.7)')ndum, kdum, &
846 ob % profiler(iprofiler) % u(k) % y, &
847 ob % profiler(iprofiler) % v(k) % y
848 end do
849 end do
850 elseif ( index( ob_name,'bogus') > 0 ) then
851 do n = 1, num_obs
852 ibogus = ibogus + 1
853 read(y_unit,'(i8)')kdum
854 read(y_unit,'(2i8,e15.7)')ndum, kdum, &
855 ob % bogus(ibogus) % slp % y
856 read(y_unit,'(i8)')num_levs
857 ob % bogus(ibogus) % numlevs = num_levs
858 allocate( ob % bogus(ibogus) % u(1:num_levs) )
859 allocate( ob % bogus(ibogus) % v(1:num_levs) )
860 allocate( ob % bogus(ibogus) % t(1:num_levs) )
861 allocate( ob % bogus(ibogus) % q(1:num_levs) )
862
863 do k = 1, num_levs
864 read(y_unit,'(2i8,4e15.7)')ndum, kdum, &
865 ob % bogus(ibogus) % u(k) % y, &
866 ob % bogus(ibogus) % v(k) % y, &
867 ob % bogus(ibogus) % t(k) % y, &
868 ob % bogus(ibogus) % q(k) % y
869 end do
870 end do
871 elseif ( index( ob_name,'gpsref') > 0 ) then
872 do n = 1, num_obs
873 igpsref = igpsref + 1
874 read(y_unit,'(i8)')num_levs
875 ob % gpsref(igpsref) % numlevs = num_levs
876 allocate( ob % gpsref(igpsref) % ref(1:num_levs) )
877
878 do k = 1, num_levs
879 read(y_unit,'(2i8,7e15.7)')ndum, kdum, &
880 ob % gpsref(igpsref) % ref(k) % y
881 end do
882 end do
883 ! Radiance obs: consistent with RTTOV triplet and WRF-VAR
884 !--------------------------------------------------------------------
885 else if ( index( ob_name,'noaa') > 0 ) then
886 platform_id = 1
887 read (ob_name,'(a4,a1,i2,a1,a5,a1,i4)') &
888 platform, str1,satellite_id,str2,sensor,str3,ichan
889 if ( sensor == 'amsua' ) then
890 sensor_id = 3
891 else if ( sensor == 'amsub' ) then
892 sensor_id = 4
893 else
894 write(6,*) ' Unrecognized Sensor '
895 end if
896 do n = 1, rtminit_nsensor
897 if ( platform_id == rtminit_platform(n) &
898 .and. satellite_id == rtminit_satid(n) &
899 .and. sensor_id == rtminit_sensor(n) ) then
900 do k = 1,num_obs
901 ob%rad(n)%num_rad_tot(ichan) = ob%rad(n)%num_rad_tot(ichan) + 1
902 read(y_unit,'(2i8,e15.7)') ipixel,kdum, &
903 ob%rad(n)%tb(ichan)%pixel(ob%rad(n)%num_rad_tot(ichan))%y
904 end do
905 exit
906 end if
907 end do
908
909 elseif ( index( ob_name,'*****') > 0 ) then
910 exit
911 else
912 print*,' unknown obs type: ',trim(ob_name),' found on unit ',y_unit
913 end if
914 end do
915 1000 print*,' end of file reached on unit ',y_unit
916
917 end subroutine da_read_y
918
919
920
921 subroutine da_read_yp( yp_unit, ob )
922
923 implicit none
924
925 integer, intent(in) :: yp_unit
926 type (iv_type), intent(inout) :: ob
927
928 character*20 :: ob_name
929 integer :: n, ndum, k, kdum, num_obs, num_levs
930 integer :: isynop, imetar, iships, ipolaramv, igeoamv, igpspw, isound, &
931 iairep, ipilot, issmir, isatem, issmt1, issmt2, iairsr, &
932 ibuoy, isonde_sfc, iqscat, ibogus, iprofiler, igpsref
933
934 rewind (yp_unit)
935
936 isynop = 0
937 imetar = 0
938 iships = 0
939 ipolaramv = 0
940 igeoamv = 0
941 igpspw = 0
942 isound = 0
943 iairsr = 0
944 iairep = 0
945 ipilot = 0
946 issmir = 0
947 isatem = 0
948 issmt1 = 0
949 issmt2 = 0
950 ibuoy = 0
951 isonde_sfc = 0
952 iqscat = 0
953 ibogus = 0
954 iprofiler = 0
955 igpsref = 0
956
957 do n = 1,rtminit_nsensor
958 ob%rad(n)%num_rad_tot(:) = 0
959 end do
960
961 do
962 read(yp_unit,'(a20,i8)', end=1000)ob_name, num_obs
963
964 if ( index( ob_name,'synop') > 0 ) then
965 do n = 1, num_obs
966 isynop = isynop + 1
967 read(yp_unit,'(i8)')num_levs
968 read(yp_unit,'(2i8,7e15.7)') &
969 ndum, kdum, ob % synop(isynop) % u % yp, &
970 ob % synop(isynop) % v % yp, &
971 ob % synop(isynop) % t % yp, &
972 ob % synop(isynop) % p % yp, &
973 ob % synop(isynop) % q % yp
974 end do
975 elseif ( index( ob_name,'metar') > 0 ) then
976 do n = 1, num_obs
977 imetar = imetar + 1
978 read(yp_unit,'(i8)')num_levs
979 read(yp_unit,'(2i8,7e15.7)') &
980 ndum, kdum, ob % metar(imetar) % u % yp, &
981 ob % metar(imetar) % v % yp, &
982 ob % metar(imetar) % t % yp, &
983 ob % metar(imetar) % p % yp, &
984 ob % metar(imetar) % q % yp
985 end do
986 elseif ( index( ob_name,'ships') > 0 ) then
987 do n = 1, num_obs
988 iships = iships + 1
989 read(yp_unit,'(i8)')num_levs
990 read(yp_unit,'(2i8,7e15.7)')ndum, kdum, ob % ships(iships) % u % yp, &
991 ob % ships(iships) % v % yp, &
992 ob % ships(iships) % t % yp, &
993 ob % ships(iships) % p % yp, &
994 ob % ships(iships) % q % yp
995 end do
996
997 elseif ( index( ob_name,'geoamv') > 0 ) then
998 do n = 1, num_obs
999 igeoamv = igeoamv + 1
1000 read(yp_unit,'(i8)')num_levs
1001 do k = 1, num_levs
1002 read(yp_unit,'(2i8,7e15.7)')ndum, kdum, &
1003 ob % geoamv(igeoamv) % u(k) % yp, &
1004 ob % geoamv(igeoamv) % v(k) % yp
1005 end do
1006 end do
1007
1008 elseif ( index( ob_name,'polaramv') > 0 ) then
1009 do n = 1, num_obs
1010 ipolaramv = ipolaramv + 1
1011 read(yp_unit,'(i8)')num_levs
1012 do k = 1, num_levs
1013 read(yp_unit,'(2i8,7e15.7)')ndum, kdum, &
1014 ob % polaramv(ipolaramv) % u(k) % yp, &
1015 ob % polaramv(ipolaramv) % v(k) % yp
1016 end do
1017 end do
1018
1019 elseif ( index( ob_name,'gpspw') > 0 ) then
1020
1021 do n = 1, num_obs
1022 igpspw = igpspw + 1
1023 read(yp_unit,'(i8)')num_levs
1024 read(yp_unit,'(2i8,7e15.7)')ndum, kdum, &
1025 ob % gpspw(igpspw) % tpw % yp
1026 end do
1027
1028 elseif ( index( ob_name,'sound') > 0 ) then
1029 do n = 1, num_obs
1030 isound = isound + 1
1031 read(yp_unit,'(i8)')num_levs
1032 do k = 1, num_levs
1033 read(yp_unit,'(2i8,7e15.7)')ndum, kdum, &
1034 ob % sound(isound) % u(k) % yp, &
1035 ob % sound(isound) % v(k) % yp, &
1036 ob % sound(isound) % t(k) % yp, &
1037 ob % sound(isound) % q(k) % yp
1038 end do
1039 end do
1040
1041 elseif ( index( ob_name,'airsr') > 0 ) then
1042 do n = 1, num_obs
1043 iairsr = iairsr + 1
1044 read(yp_unit,'(i8)')num_levs
1045 do k = 1, num_levs
1046 read(yp_unit,'(2i8,7e15.7)')ndum, kdum, &
1047 ob % airsr(iairsr) % t(k) % yp, &
1048 ob % airsr(iairsr) % q(k) % yp
1049 end do
1050 end do
1051 elseif ( index( ob_name,'airep') > 0 ) then
1052 do n = 1, num_obs
1053 iairep = iairep + 1
1054 read(yp_unit,'(i8)')num_levs
1055 do k = 1, num_levs
1056 read(yp_unit,'(2i8,7e15.7)')ndum, kdum, &
1057 ob % airep(iairep) % u(k) % yp, &
1058 ob % airep(iairep) % v(k) % yp, &
1059 ob % airep(iairep) % t(k) % yp
1060 end do
1061 end do
1062 elseif ( index( ob_name,'pilot') > 0 ) then
1063 do n = 1, num_obs
1064 ipilot = ipilot + 1
1065 read(yp_unit,'(i8)')num_levs
1066 do k = 1, num_levs
1067 read(yp_unit,'(2i8,7e15.7)')ndum, kdum, &
1068 ob % pilot(ipilot) % u(k) % yp, &
1069 ob % pilot(ipilot) % v(k) % yp
1070 end do
1071 end do
1072 elseif ( index( ob_name,'ssmir') > 0 ) then
1073 do n = 1, num_obs
1074 issmir = issmir + 1
1075 read(yp_unit,'(i8)')num_levs
1076 read(yp_unit,'(2i8,7e15.7)')ndum, kdum, &
1077 ob % ssmir(issmir) % speed % yp, &
1078 ob % ssmir(issmir) % tpw % yp
1079 end do
1080 elseif ( index( ob_name,'satem') > 0 ) then
1081 do n = 1, num_obs
1082 isatem = isatem + 1
1083 read(yp_unit,'(i8)')num_levs
1084 do k = 1, num_levs
1085 read(yp_unit,'(2i8,7e15.7)')ndum, kdum, &
1086 ob % satem(isatem) % thickness(k) % yp
1087 end do
1088 end do
1089 elseif ( index( ob_name,'ssmt1') > 0 ) then
1090 do n = 1, num_obs
1091 issmt1 = issmt1 + 1
1092 read(yp_unit,'(i8)')num_levs
1093 do k = 1, num_levs
1094 read(yp_unit,'(2i8,7e15.7)')ndum, kdum, &
1095 ob % ssmt1(issmt1) % t(k) % yp
1096 end do
1097 end do
1098 elseif ( index( ob_name,'ssmt2') > 0 ) then
1099 do n = 1, num_obs
1100 issmt2 = issmt2 + 1
1101 read(yp_unit,'(i8)')num_levs
1102 do k = 1, num_levs
1103 read(yp_unit,'(2i8,7e15.7)')ndum, kdum, &
1104 ob % ssmt2(issmt2) % rh(k) % yp
1105 end do
1106 end do
1107 elseif ( index( ob_name,'buoy') > 0 ) then
1108 do n = 1, num_obs
1109 ibuoy = ibuoy + 1
1110 read(yp_unit,'(i8)')num_levs
1111 read(yp_unit,'(2i8,7e15.7)')ndum, kdum, ob % buoy(ibuoy) % u % yp, &
1112 ob % buoy(ibuoy) % v % yp, &
1113 ob % buoy(ibuoy) % t % yp, &
1114 ob % buoy(ibuoy) % p % yp, &
1115 ob % buoy(ibuoy) % q % yp
1116 end do
1117 elseif ( index( ob_name,'sonde_sfc') > 0 ) then
1118 do n = 1, num_obs
1119 isonde_sfc = isonde_sfc + 1
1120 read(yp_unit,'(i8)')num_levs
1121 read(yp_unit,'(2i8,7e15.7)')ndum, kdum, ob % sonde_sfc(isonde_sfc) % u % yp, &
1122 ob % sonde_sfc(isonde_sfc) % v % yp, &
1123 ob % sonde_sfc(isonde_sfc) % t % yp, &
1124 ob % sonde_sfc(isonde_sfc) % p % yp, &
1125 ob % sonde_sfc(isonde_sfc) % q % yp
1126 end do
1127
1128 elseif ( index( ob_name,'qscat') > 0 ) then
1129 do n = 1, num_obs
1130 iqscat = iqscat + 1
1131 read(yp_unit,'(i8)')num_levs
1132 read(yp_unit,'(2i8,7e15.7)')ndum, kdum, ob % qscat(iqscat) % u % yp, &
1133 ob % qscat(iqscat) % v % yp
1134 end do
1135 elseif ( index( ob_name,'profiler') > 0 ) then
1136 do n = 1, num_obs
1137 iprofiler = iprofiler + 1
1138 read(yp_unit,'(i8)')num_levs
1139 ob % profiler(iprofiler) % numlevs = num_levs
1140 do k = 1, num_levs
1141 read(yp_unit,'(2i8,7e15.7)')ndum, kdum, &
1142 ob % profiler(iprofiler) % u(k) % yp, &
1143 ob % profiler(iprofiler) % v(k) % yp
1144 end do
1145 end do
1146 elseif ( index( ob_name,'bogus') > 0 ) then
1147 do n = 1, num_obs
1148 ibogus = ibogus + 1
1149 read(yp_unit,'(i8)') kdum
1150 read(yp_unit,'(2i8,e15.7)')ndum, kdum, ob%bogus(ibogus)%slp%yp
1151 read(yp_unit,'(i8)')num_levs
1152 do k = 1, num_levs
1153 read(yp_unit,'(2i8,4e15.7)')ndum, kdum, &
1154 ob % bogus(ibogus) % u(k) % yp, &
1155 ob % bogus(ibogus) % v(k) % yp, &
1156 ob % bogus(ibogus) % t(k) % yp, &
1157 ob % bogus(ibogus) % q(k) % yp
1158 end do
1159 end do
1160 elseif ( index( ob_name,'gpsref') > 0 ) then
1161 do n = 1, num_obs
1162 igpsref = igpsref + 1
1163 read(yp_unit,'(i8)')num_levs
1164 do k = 1, num_levs
1165 read(yp_unit,'(2i8,7e15.7)')ndum, kdum, &
1166 ob % gpsref(igpsref) % ref(k) % yp
1167 end do
1168 end do
1169 ! Radiance obs: consistent with RTTOV triplet and WRF-VAR
1170 !--------------------------------------------------------------------
1171 elseif ( index( ob_name,'noaa') > 0 ) then
1172 platform_id = 1
1173 read (ob_name,'(a4,a1,i2,a1,a5,a1,i4)') &
1174 platform, str1,satellite_id,str2,sensor,str3,ichan
1175 if ( sensor == 'amsua' ) then
1176 sensor_id = 3
1177 else if ( sensor == 'amsub' ) then
1178 sensor_id = 4
1179 else
1180 write(6,*) ' Unrecognized Sensor '
1181 end if
1182 do n = 1, rtminit_nsensor
1183 if ( platform_id == rtminit_platform(n) &
1184 .and. satellite_id == rtminit_satid(n) &
1185 .and. sensor_id == rtminit_sensor(n) ) then
1186 do k = 1,num_obs
1187 ob%rad(n)%num_rad_tot(ichan) = ob%rad(n)%num_rad_tot(ichan)+1
1188 read(yp_unit,'(2i8,e15.7)') ipixel, kdum, &
1189 ob%rad(n)%tb(ichan)%pixel(ob%rad(n)%num_rad_tot(ichan))%yp
1190 end do
1191 exit
1192 end if
1193 end do
1194 elseif ( index( ob_name,'*****') > 0 ) then
1195 exit
1196 else
1197 print*,' unknown obs type: ',trim(ob_name),' found on unit ',yp_unit
1198 end if
1199 end do
1200 1000 print*,' end of file reached on unit ',yp_unit
1201
1202 end subroutine da_read_yp
1203
1204
1205
1206 subroutine da_read_obs_rand( rand_unit, ob )
1207
1208 implicit none
1209
1210 integer, intent(in) :: rand_unit
1211 type (iv_type), intent(inout) :: ob
1212
1213 character*20 :: ob_name
1214 integer :: n, ndum, k, kdum, num_obs, num_levs
1215 integer :: isynop, imetar, iships, ipolaramv, igeoamv, igpspw, isound, &
1216 iairep, ipilot, issmir, isatem, issmt1, issmt2, iairsr, &
1217 ibuoy, isonde_sfc, iqscat, ibogus, iprofiler, igpsref
1218
1219 isynop = 0
1220 imetar = 0
1221 iships = 0
1222 ipolaramv = 0
1223 igeoamv = 0
1224 igpspw = 0
1225 isound = 0
1226 iairsr = 0
1227 iairep = 0
1228 ipilot = 0
1229 issmir = 0
1230 isatem = 0
1231 issmt1 = 0
1232 issmt2 = 0
1233 ibuoy = 0
1234 isonde_sfc = 0
1235 iqscat = 0
1236 ibogus = 0
1237 iprofiler = 0
1238 igpsref = 0
1239 do n = 1,rtminit_nsensor
1240 ob%rad(n)%num_rad_tot(:) = 0
1241 end do
1242
1243 rewind( rand_unit )
1244
1245 do
1246 read(rand_unit,'(a20,i8)',end=1000)ob_name, num_obs
1247 if ( index( ob_name,'synop') > 0 ) then
1248 do n = 1, num_obs
1249 isynop = isynop + 1
1250 read(rand_unit,'(i8)')num_levs
1251 read(rand_unit,'(2i8,10e15.7)')ndum, kdum, &
1252 ob % synop(isynop) % u % error, ob % synop(isynop) % u % pert, &
1253 ob % synop(isynop) % v % error, ob % synop(isynop) % v % pert, &
1254 ob % synop(isynop) % t % error, ob % synop(isynop) % t % pert, &
1255 ob % synop(isynop) % p % error, ob % synop(isynop) % p % pert, &
1256 ob % synop(isynop) % q % error, ob % synop(isynop) % q % pert
1257 end do
1258 elseif ( index( ob_name,'metar') > 0 ) then
1259
1260 do n = 1, num_obs
1261 imetar = imetar + 1
1262 read(rand_unit,'(i8)')num_levs
1263 read(rand_unit,'(2i8,10e15.7)')ndum, kdum, &
1264 ob % metar(imetar) % u % error, ob % metar(imetar) % u % pert, &
1265 ob % metar(imetar) % v % error, ob % metar(imetar) % v % pert, &
1266 ob % metar(imetar) % t % error, ob % metar(imetar) % t % pert, &
1267 ob % metar(imetar) % p % error, ob % metar(imetar) % p % pert, &
1268 ob % metar(imetar) % q % error, ob % metar(imetar) % q % pert
1269 end do
1270
1271 elseif ( index( ob_name,'ships') > 0 ) then
1272
1273 do n = 1, num_obs
1274 iships = iships + 1
1275 read(rand_unit,'(i8)')num_levs
1276 read(rand_unit,'(2i8,10e15.7)')ndum, kdum, &
1277 ob % ships(iships) % u % error, ob % ships(iships) % u % pert, &
1278 ob % ships(iships) % v % error, ob % ships(iships) % v % pert, &
1279 ob % ships(iships) % t % error, ob % ships(iships) % t % pert, &
1280 ob % ships(iships) % p % error, ob % ships(iships) % p % pert, &
1281 ob % ships(iships) % q % error, ob % ships(iships) % q % pert
1282 end do
1283
1284 elseif ( index( ob_name,'geoamv') > 0 ) then
1285
1286 do n = 1, num_obs
1287 igeoamv = igeoamv + 1
1288 read(rand_unit,'(i8)')num_levs
1289 do k = 1, num_levs
1290 read(rand_unit,'(2i8,10e15.7)')ndum, kdum, &
1291 ob % geoamv(igeoamv) % u(k) % error, ob % geoamv(igeoamv) % u(k) % pert, &
1292 ob % geoamv(igeoamv) % v(k) % error, ob % geoamv(igeoamv) % v(k) % pert
1293 end do
1294 end do
1295
1296 elseif ( index( ob_name,'polaramv') > 0 ) then
1297
1298 do n = 1, num_obs
1299 ipolaramv = ipolaramv + 1
1300 read(rand_unit,'(i8)')num_levs
1301 do k = 1, num_levs
1302 read(rand_unit,'(2i8,10e15.7)')ndum, kdum, &
1303 ob % polaramv(ipolaramv) % u(k) % error, &
1304 ob % polaramv(ipolaramv) % u(k) % pert, &
1305 ob % polaramv(ipolaramv) % v(k) % error, &
1306 ob % polaramv(ipolaramv) % v(k) % pert
1307 end do
1308 end do
1309
1310 elseif ( index( ob_name,'gpspw') > 0 ) then
1311
1312 do n = 1, num_obs
1313 igpspw = igpspw + 1
1314 read(rand_unit,'(i8)')num_levs
1315 read(rand_unit,'(2i8,10e15.7)')ndum, kdum, &
1316 ob % gpspw(igpspw) % tpw % error, ob % gpspw(igpspw) % tpw % pert
1317 end do
1318
1319 elseif ( index( ob_name,'sound') > 0 ) then
1320
1321 do n = 1, num_obs
1322 isound = isound + 1
1323 read(rand_unit,'(i8)')num_levs
1324
1325 do k = 1, num_levs
1326 read(rand_unit,'(2i8,10e15.7)')ndum, kdum, &
1327 ob % sound(isound) % u(k) % error, &
1328 ob % sound(isound) % u(k) % pert, &
1329 ob % sound(isound) % v(k) % error, &
1330 ob % sound(isound) % v(k) % pert, &
1331 ob % sound(isound) % t(k) % error, &
1332 ob % sound(isound) % t(k) % pert, &
1333 ob % sound(isound) % q(k) % error, &
1334 ob % sound(isound) % q(k) % pert
1335 end do
1336 end do
1337
1338 elseif ( index( ob_name,'airsr') > 0 ) then
1339
1340 do n = 1, num_obs
1341 iairsr = iairsr + 1
1342 read(rand_unit,'(i8)')num_levs
1343
1344 do k = 1, num_levs
1345 read(rand_unit,'(2i8,10e15.7)')ndum, kdum, &
1346 ob % airsr(iairsr) % t(k) % error, &
1347 ob % airsr(iairsr) % t(k) % pert, &
1348 ob % airsr(iairsr) % q(k) % error, &
1349 ob % airsr(iairsr) % q(k) % pert
1350 end do
1351 end do
1352
1353 elseif ( index( ob_name,'airep') > 0 ) then
1354
1355 do n = 1, num_obs
1356 iairep = iairep + 1
1357 read(rand_unit,'(i8)')num_levs
1358
1359 do k = 1, num_levs
1360 read(rand_unit,'(2i8,10e15.7)')ndum, kdum, &
1361 ob % airep(iairep) % u(k) % error, &
1362 ob % airep(iairep) % u(k) % pert, &
1363 ob % airep(iairep) % v(k) % error, &
1364 ob % airep(iairep) % v(k) % pert, &
1365 ob % airep(iairep) % t(k) % error, &
1366 ob % airep(iairep) % t(k) % pert
1367 end do
1368 end do
1369
1370 elseif ( index( ob_name,'pilot') > 0 ) then
1371
1372 do n = 1, num_obs
1373 ipilot = ipilot + 1
1374 read(rand_unit,'(i8)')num_levs
1375
1376 do k = 1, num_levs
1377 read(rand_unit,'(2i8,10e15.7)')ndum, kdum, &
1378 ob % pilot(ipilot) % u(k) % error, &
1379 ob % pilot(ipilot) % u(k) % pert, &
1380 ob % pilot(ipilot) % v(k) % error, &
1381 ob % pilot(ipilot) % v(k) % pert
1382 end do
1383 end do
1384
1385 elseif ( index( ob_name,'ssmir') > 0 ) then
1386
1387 do n = 1, num_obs
1388 issmir = issmir + 1
1389 read(rand_unit,'(i8)')num_levs
1390 read(rand_unit,'(2i8,10e15.7)')ndum, kdum, &
1391 ob % ssmir(issmir) % speed % error, &
1392 ob % ssmir(issmir) % speed % pert, &
1393 ob % ssmir(issmir) % tpw % error, &
1394 ob % ssmir(issmir) % tpw % pert
1395 end do
1396
1397 elseif ( index( ob_name,'satem') > 0 ) then
1398
1399 do n = 1, num_obs
1400 isatem = isatem + 1
1401 read(rand_unit,'(i8)')num_levs
1402
1403 do k = 1, num_levs
1404 read(rand_unit,'(2i8,10e15.7)')ndum, kdum, &
1405 ob % satem(isatem) % thickness(k) % error, &
1406 ob % satem(isatem) % thickness(k) % pert
1407 end do
1408 end do
1409
1410 elseif ( index( ob_name,'ssmt1') > 0 ) then
1411
1412 do n = 1, num_obs
1413 issmt1 = issmt1 + 1
1414 read(rand_unit,'(i8)')num_levs
1415
1416 do k = 1, num_levs
1417 read(rand_unit,'(2i8,10e15.7)')ndum, kdum, &
1418 ob % ssmt1(issmt1) % t % error, ob % ssmt1(issmt1) % t % pert
1419 end do
1420 end do
1421
1422 elseif ( index( ob_name,'ssmt2') > 0 ) then
1423
1424 do n = 1, num_obs
1425 issmt2 = issmt2 + 1
1426 read(rand_unit,'(i8)')num_levs
1427
1428 do k = 1, num_levs
1429 read(rand_unit,'(2i8,10e15.7)')ndum, kdum, &
1430 ob % ssmt2(issmt2) % rh % error, &
1431 ob % ssmt2(issmt2) % rh % pert
1432 end do
1433 end do
1434 elseif ( index( ob_name,'buoy') > 0 ) then
1435 do n = 1, num_obs
1436 ibuoy = ibuoy + 1
1437 read(rand_unit,'(i8)')num_levs
1438 read(rand_unit,'(2i8,10e15.7)')ndum, kdum, &
1439 ob % buoy(ibuoy) % u % error, ob % buoy(ibuoy) % u % pert, &
1440 ob % buoy(ibuoy) % v % error, ob % buoy(ibuoy) % v % pert, &
1441 ob % buoy(ibuoy) % t % error, ob % buoy(ibuoy) % t % pert, &
1442 ob % buoy(ibuoy) % p % error, ob % buoy(ibuoy) % p % pert, &
1443 ob % buoy(ibuoy) % q % error, ob % buoy(ibuoy) % q % pert
1444 end do
1445 elseif ( index( ob_name,'sonde_sfc') > 0 ) then
1446 do n = 1, num_obs
1447 isonde_sfc = isonde_sfc + 1
1448 read(rand_unit,'(i8)')num_levs
1449 read(rand_unit,'(2i8,10e15.7)')ndum, kdum, &
1450 ob % sonde_sfc(isonde_sfc) % u % error, &
1451 ob % sonde_sfc(isonde_sfc) % u % pert, &
1452 ob % sonde_sfc(isonde_sfc) % v % error, &
1453 ob % sonde_sfc(isonde_sfc) % v % pert, &
1454 ob % sonde_sfc(isonde_sfc) % t % error, &
1455 ob % sonde_sfc(isonde_sfc) % t % pert, &
1456 ob % sonde_sfc(isonde_sfc) % p % error, &
1457 ob % sonde_sfc(isonde_sfc) % p % pert, &
1458 ob % sonde_sfc(isonde_sfc) % q % error, &
1459 ob % sonde_sfc(isonde_sfc) % q % pert
1460 end do
1461
1462 elseif ( index( ob_name,'qscat') > 0 ) then
1463 do n = 1, num_obs
1464 iqscat = iqscat + 1
1465 read(rand_unit,'(i8)')num_levs
1466 read(rand_unit,'(2i8,10e15.7)')ndum, kdum, &
1467 ob % qscat(iqscat) % u % error, ob % qscat(iqscat) % u % pert, &
1468 ob % qscat(iqscat) % v % error, ob % qscat(iqscat) % v % pert
1469 end do
1470 elseif ( index( ob_name,'profiler') > 0 ) then
1471 do n = 1, num_obs
1472 iprofiler = iprofiler + 1
1473 read(rand_unit,'(i8)')num_levs
1474 do k = 1, num_levs
1475 read(rand_unit,'(2i8,10e15.7)')ndum, kdum, &
1476 ob % profiler(iprofiler) % u(k) % error, &
1477 ob % profiler(iprofiler) % u(k) % pert, &
1478 ob % profiler(iprofiler) % v(k) % error, &
1479 ob % profiler(iprofiler) % v(k) % pert
1480 end do
1481 end do
1482 elseif ( index( ob_name,'bogus') > 0 ) then
1483 do n = 1, num_obs
1484 ibogus = ibogus + 1
1485 read(rand_unit,'(i8)')num_levs
1486 do k = 1, num_levs
1487 read(rand_unit,'(2i8,8e15.7)')ndum, kdum, &
1488 ob % bogus(ibogus) % u(k) % error, &
1489 ob % bogus(ibogus) % u(k) % pert, &
1490 ob % bogus(ibogus) % v(k) % error, &
1491 ob % bogus(ibogus) % v(k) % pert, &
1492 ob % bogus(ibogus) % t(k) % error, &
1493 ob % bogus(ibogus) % t(k) % pert, &
1494 ob % bogus(ibogus) % q(k) % error, &
1495 ob % bogus(ibogus) % q(k) % pert
1496 end do
1497 read(rand_unit,'(2i8,2e15.7)')ndum, kdum, &
1498 ob % bogus(ibogus) % slp % error, ob % bogus(ibogus) % slp % pert
1499 end do
1500 elseif ( index( ob_name,'gpsref') > 0 ) then
1501 do n = 1, num_obs
1502 igpsref = igpsref + 1
1503 read(rand_unit,'(i8)')num_levs
1504 do k = 1, num_levs
1505 read(rand_unit,'(2i8,7e15.7)')ndum, kdum, &
1506 ob % gpsref(igpsref) % ref(k) % error, &
1507 ob % gpsref(igpsref) % ref(k) % pert
1508 end do
1509 end do
1510
1511 ! Radiance obs: consistent with RTTOV triplet and WRF-VAR
1512 !--------------------------------------------------------------------
1513 elseif ( index( ob_name,'noaa') > 0 ) then
1514 platform_id = 1
1515 read (ob_name,'(a4,a1,i2,a1,a5,a1,i4)') &
1516 platform, str1,satellite_id,str2,sensor,str3,ichan
1517 if ( sensor == 'amsua' ) then
1518 sensor_id = 3
1519 else if ( sensor == 'amsub' ) then
1520 sensor_id = 4
1521 else
1522 write(6,*) ' Unrecognized Sensor '
1523 end if
1524 do n = 1, rtminit_nsensor
1525 if ( platform_id == rtminit_platform(n) &
1526 .and. satellite_id == rtminit_satid(n) &
1527 .and. sensor_id == rtminit_sensor(n) ) then
1528 do k = 1,num_obs
1529 ob%rad(n)%num_rad_tot(ichan) = ob%rad(n)%num_rad_tot(ichan) + 1
1530 read(rand_unit,'(2i8,f10.3,e15.7)') ipixel, kdum, &
1531 ob%rad(n)%tb(ichan)%pixel( &
1532 ob%rad(n)%num_rad_tot(ichan))%error, &
1533 ob%rad(n)%tb(ichan)%pixel( &
1534 ob%rad(n)%num_rad_tot(ichan))%pert
1535 end do
1536 exit
1537 end if
1538 end do
1539 elseif ( index( ob_name,'*****') > 0 ) then
1540 exit
1541 else
1542 print*,' unknown obs type: ',trim(ob_name),' found on unit ',rand_unit
1543 end if
1544 end do
1545 1000 print*,' end of file reached on unit ',rand_unit
1546
1547 end subroutine da_read_obs_rand
1548
1549 subroutine da_calc_jo_expected( ob )
1550
1551 implicit none
1552
1553 type (iv_type), intent(inout) :: ob
1554
1555 integer :: n, k
1556 integer :: count1, count2, count3, count4, count5
1557 real :: trace1, trace2, trace3, trace4, trace5
1558
1559 ob % trace_total = 0
1560
1561 ob % num_synop_tot = 0
1562 ob % trace_synop = 0
1563 if ( ob % num_synop > 0 ) then
1564 trace1 = 0.0; trace2 = 0.0; trace3 = 0.0; trace4 = 0.0; trace5 = 0.0
1565 count1 = 0; count2 = 0; count3 = 0; count4 = 0; count5 = 0
1566
1567 do n = 1, ob % num_synop
1568 call da_calc_trace_single( ob % synop(n) % u, count1, trace1 )
1569 call da_calc_trace_single( ob % synop(n) % v, count2, trace2 )
1570 call da_calc_trace_single( ob % synop(n) % t, count3, trace3 )
1571 call da_calc_trace_single( ob % synop(n) % p, count4, trace4 )
1572 call da_calc_trace_single( ob % synop(n) % q, count5, trace5 )
1573 end do
1574
1575 ob % jo_synop_u = 0.5 * ( count1 - trace1 )
1576 ob % jo_synop_v = 0.5 * ( count2 - trace2 )
1577 ob % jo_synop_t = 0.5 * ( count3 - trace3 )
1578 ob % jo_synop_p = 0.5 * ( count4 - trace4 )
1579 ob % jo_synop_q = 0.5 * ( count5 - trace5 )
1580 ob % trace_synop = trace1 + trace2 + trace3 + trace4 + trace5
1581
1582 if ( ob % trace_synop < 0.0 ) &
1583 write(6,'(a,f15.5)')' Warning: ob % trace_synop < 0 = ', ob % trace_synop
1584
1585 ob % trace_total = ob % trace_total + ob % trace_synop
1586 ob % num_synop_tot = count1 + count2 + count3 + count4 + count5
1587 end if
1588
1589 ob % num_metar_tot = 0
1590 ob % trace_metar = 0
1591 if ( ob % num_metar > 0 ) then
1592 trace1 = 0.0; trace2 = 0.0; trace3 = 0.0; trace4 = 0.0; trace5 = 0.0
1593 count1 = 0; count2 = 0; count3 = 0; count4 = 0; count5 = 0
1594
1595 do n = 1, ob % num_metar
1596 call da_calc_trace_single( ob % metar(n) % u, count1, trace1 )
1597 call da_calc_trace_single( ob % metar(n) % v, count2, trace2 )
1598 call da_calc_trace_single( ob % metar(n) % t, count3, trace3 )
1599 call da_calc_trace_single( ob % metar(n) % p, count4, trace4 )
1600 call da_calc_trace_single( ob % metar(n) % q, count5, trace5 )
1601 end do
1602
1603 ob % jo_metar_u = 0.5 * ( count1 - trace1 )
1604 ob % jo_metar_v = 0.5 * ( count2 - trace2 )
1605 ob % jo_metar_t = 0.5 * ( count3 - trace3 )
1606 ob % jo_metar_p = 0.5 * ( count4 - trace4 )
1607 ob % jo_metar_q = 0.5 * ( count5 - trace5 )
1608 ob % trace_metar = trace1 + trace2 + trace3 + trace4 + trace5
1609
1610 if ( ob % trace_metar < 0.0 ) &
1611 write(6,'(a,f15.5)') &
1612 ' Warning: ob % trace_metar < 0 = ', ob % trace_metar
1613
1614 ob % trace_total = ob % trace_total + ob % trace_metar
1615 ob % num_metar_tot = count1 + count2 + count3 + count4 + count5
1616
1617 end if
1618
1619 ob % num_ships_tot = 0
1620 ob % trace_ships = 0
1621 if ( ob % num_ships > 0 ) then
1622 trace1 = 0.0; trace2 = 0.0; trace3 = 0.0; trace4 = 0.0; trace5 = 0.0
1623 count1 = 0; count2 = 0; count3 = 0; count4 = 0; count5 = 0
1624
1625 do n = 1, ob % num_ships
1626 call da_calc_trace_single( ob % ships(n) % u, count1, trace1 )
1627 call da_calc_trace_single( ob % ships(n) % v, count2, trace2 )
1628 call da_calc_trace_single( ob % ships(n) % t, count3, trace3 )
1629 call da_calc_trace_single( ob % ships(n) % p, count4, trace4 )
1630 call da_calc_trace_single( ob % ships(n) % q, count5, trace5 )
1631 end do
1632
1633 ob % jo_ships_u = 0.5 * ( count1 - trace1 )
1634 ob % jo_ships_v = 0.5 * ( count2 - trace2 )
1635 ob % jo_ships_t = 0.5 * ( count3 - trace3 )
1636 ob % jo_ships_p = 0.5 * ( count4 - trace4 )
1637 ob % jo_ships_q = 0.5 * ( count5 - trace5 )
1638 ob % trace_ships = trace1 + trace2 + trace3 + trace4 + trace5
1639
1640 if ( ob % trace_ships < 0.0 ) &
1641 write(6,'(a,f15.5)') &
1642 ' Warning: ob % trace_ships < 0 = ', ob % trace_ships
1643 ob % trace_total = ob % trace_total + ob % trace_ships
1644 ob % num_ships_tot = count1 + count2 + count3 + count4 + count5
1645 end if
1646
1647 ob % num_geoamv_tot = 0
1648 ob % trace_geoamv = 0
1649 if ( ob % num_geoamv > 0 ) then
1650 trace1 = 0.0; trace2 = 0.0
1651 count1 = 0; count2 = 0
1652
1653 do n = 1, ob % num_geoamv
1654 do k = 1, ob % geoamv(n) % numlevs
1655 call da_calc_trace_single( ob % geoamv(n) % u(k), count1, trace1 )
1656 call da_calc_trace_single( ob % geoamv(n) % v(k), count2, trace2 )
1657 end do
1658 end do
1659
1660 ob % jo_geoamv_u = 0.5 * ( count1 - trace1 )
1661 ob % jo_geoamv_v = 0.5 * ( count2 - trace2 )
1662 ob % trace_geoamv = trace1 + trace2
1663
1664 if ( ob % trace_geoamv < 0.0 ) &
1665 write(6,'(a,f15.5)') &
1666 ' Warning: ob % trace_geoamv < 0 = ', ob % trace_geoamv
1667
1668 ob % trace_total = ob % trace_total + ob % trace_geoamv
1669 ob % num_geoamv_tot = count1 + count2
1670 end if
1671
1672 ob % num_polaramv_tot = 0
1673 ob % trace_polaramv = 0
1674 if ( ob % num_polaramv > 0 ) then
1675 trace1 = 0.0; trace2 = 0.0
1676 count1 = 0; count2 = 0
1677
1678 do n = 1, ob % num_polaramv
1679 do k = 1, ob % polaramv(n) % numlevs
1680 call da_calc_trace_single( ob % polaramv(n) % u(k), count1, trace1 )
1681 call da_calc_trace_single( ob % polaramv(n) % v(k), count2, trace2 )
1682 end do
1683 end do
1684
1685 ob % jo_polaramv_u = 0.5 * ( count1 - trace1 )
1686 ob % jo_polaramv_v = 0.5 * ( count2 - trace2 )
1687 ob % trace_polaramv = trace1 + trace2
1688
1689 if ( ob % trace_polaramv < 0.0 ) &
1690 write(6,'(a,f15.5)') &
1691 ' Warning: ob % trace_polaramv < 0 = ', ob % trace_polaramv
1692
1693 ob % trace_total = ob % trace_total + ob % trace_polaramv
1694 ob % num_polaramv_tot = count1 + count2
1695 end if
1696
1697 ob % num_gpspw_tot = 0
1698 ob % trace_gpspw = 0
1699 if ( ob % num_gpspw > 0 ) then
1700 trace1 = 0.0
1701 count1 = 0
1702
1703 do n = 1, ob % num_gpspw
1704 call da_calc_trace_single( ob % gpspw(n) % tpw, count1, trace1 )
1705 end do
1706
1707 ob % jo_gpspw_tpw = 0.5 * ( count1 - trace1 )
1708 ob % trace_gpspw = trace1
1709
1710 if ( ob % trace_gpspw < 0.0 ) &
1711 write(6,'(a,f15.5)') &
1712 ' Warning: ob % trace_gpspw < 0 = ', ob % trace_gpspw
1713
1714 ob % trace_total = ob % trace_total + ob % trace_gpspw
1715 ob % num_gpspw_tot = count1
1716 end if
1717
1718 ob % num_sound_tot = 0
1719 ob % trace_sound = 0
1720 if ( ob % num_sound > 0 ) then
1721 trace1 = 0.0; trace2 = 0.0; trace3 = 0.0; trace4 = 0.0
1722 count1 = 0; count2 = 0; count3 = 0; count4 = 0
1723
1724 do n = 1, ob % num_sound
1725 do k = 1, ob % sound(n) % numlevs
1726 call da_calc_trace_single( ob % sound(n) % u(k), count1, trace1 )
1727 call da_calc_trace_single( ob % sound(n) % v(k), count2, trace2 )
1728 call da_calc_trace_single( ob % sound(n) % t(k), count3, trace3 )
1729 call da_calc_trace_single( ob % sound(n) % q(k), count4, trace4 )
1730 end do
1731 end do
1732
1733 ob % jo_sound_u = 0.5 * ( count1 - trace1 )
1734 ob % jo_sound_v = 0.5 * ( count2 - trace2 )
1735 ob % jo_sound_t = 0.5 * ( count3 - trace3 )
1736 ob % jo_sound_q = 0.5 * ( count4 - trace4 )
1737 ob % trace_sound = trace1 + trace2 + trace3 + trace4
1738
1739 if ( ob % trace_sound < 0.0 ) &
1740 write(6,'(a,f15.5)') &
1741 ' Warning: ob % trace_sound < 0 = ', ob % trace_sound
1742
1743 ob % trace_total = ob % trace_total + ob % trace_sound
1744 ob % num_sound_tot = count1 + count2 + count3 + count4
1745
1746 end if
1747
1748 ob % num_airsr_tot = 0
1749 ob % trace_airsr = 0
1750 if ( ob % num_airsr > 0 ) then
1751 trace1 = 0.0; trace2 = 0.0
1752 count1 = 0; count2 = 0
1753
1754 do n = 1, ob % num_airsr
1755 do k = 1, ob % airsr(n) % numlevs
1756 call da_calc_trace_single( ob % airsr(n) % t(k), count1, trace1 )
1757 call da_calc_trace_single( ob % airsr(n) % q(k), count2, trace2 )
1758 end do
1759 end do
1760
1761 ob % jo_airsr_t = 0.5 * ( count1 - trace1 )
1762 ob % jo_airsr_q = 0.5 * ( count2 - trace2 )
1763 ob % trace_airsr = trace1 + trace2
1764
1765 if ( ob % trace_airsr < 0.0 ) &
1766 write(6,'(a,f15.5)') &
1767 ' Warning: ob % trace_airsr < 0 = ', ob % trace_airsr
1768
1769 ob % trace_total = ob % trace_total + ob % trace_airsr
1770 ob % num_airsr_tot = count1 + count2
1771 end if
1772
1773 ob % num_airep_tot = 0
1774 ob % trace_airep = 0
1775 if ( ob % num_airep > 0 ) then
1776 trace1 = 0.0; trace2 = 0.0; trace3 = 0.0; trace4 = 0.0
1777 count1 = 0; count2 = 0; count3 = 0; count4 = 0
1778
1779 do n = 1, ob % num_airep
1780 do k = 1, ob % airep(n) % numlevs
1781 call da_calc_trace_single( ob % airep(n) % u(k), count1, trace1 )
1782 call da_calc_trace_single( ob % airep(n) % v(k), count2, trace2 )
1783 call da_calc_trace_single( ob % airep(n) % t(k), count3, trace3 )
1784 end do
1785 end do
1786
1787 ob % jo_airep_u = 0.5 * ( count1 - trace1 )
1788 ob % jo_airep_v = 0.5 * ( count2 - trace2 )
1789 ob % jo_airep_t = 0.5 * ( count3 - trace3 )
1790 ob % trace_airep = trace1 + trace2 + trace3
1791
1792 if ( ob % trace_airep < 0.0 ) &
1793 write(6,'(a,f15.5)') &
1794 ' Warning: ob % trace_airep < 0 = ', ob % trace_airep
1795
1796 ob % trace_total = ob % trace_total + ob % trace_airep
1797 ob % num_airep_tot = count1 + count2 + count3
1798 end if
1799
1800 ob % num_pilot_tot = 0
1801 ob % trace_pilot = 0
1802 if ( ob % num_pilot > 0 ) then
1803 trace1 = 0.0; trace2 = 0.0
1804 count1 = 0; count2 = 0
1805
1806 do n = 1, ob % num_pilot
1807 do k = 1, ob % pilot(n) % numlevs
1808 call da_calc_trace_single( ob % pilot(n) % u(k), count1, trace1 )
1809 call da_calc_trace_single( ob % pilot(n) % v(k), count2, trace2 )
1810 end do
1811 end do
1812
1813 ob % jo_pilot_u = 0.5 * ( count1 - trace1 )
1814 ob % jo_pilot_v = 0.5 * ( count2 - trace2 )
1815 ob % trace_pilot = trace1 + trace2
1816
1817 if ( ob % trace_pilot < 0.0 ) &
1818 write(6,'(a,f15.5)') &
1819 ' Warning: ob % trace_pilot < 0 = ', ob % trace_pilot
1820
1821 ob % trace_total = ob % trace_total + ob % trace_pilot
1822 ob % num_pilot_tot = count1 + count2
1823
1824 end if
1825
1826 ob % num_ssmir_tot = 0
1827 ob % trace_ssmir = 0
1828 if ( ob % num_ssmir > 0 ) then
1829 trace1 = 0.0; trace2 = 0.0
1830 count1= 0; count2 = 0
1831
1832 do n = 1, ob % num_ssmir
1833 call da_calc_trace_single( ob % ssmir(n) % speed, count1, trace1 )
1834 call da_calc_trace_single( ob % ssmir(n) % tpw, count2, trace2 )
1835 end do
1836
1837 ob % jo_ssmir_speed = 0.5 * ( count1 - trace1 )
1838 ob % jo_ssmir_tpw = 0.5 * ( count2 - trace2 )
1839 ob % trace_ssmir = trace1 + trace2
1840
1841 if ( ob % trace_ssmir < 0.0 ) &
1842 write(6,'(a,f15.5)') &
1843 ' Warning: ob % trace_ssmir < 0 = ', ob % trace_ssmir
1844
1845 ob % trace_total = ob % trace_total + ob % trace_ssmir
1846 ob % num_ssmir_tot = count1 + count2
1847 end if
1848
1849 ob % num_satem_tot = 0
1850 ob % trace_satem = 0
1851 if ( ob % num_satem > 0 ) then
1852 trace1 = 0.0
1853 count1 = 0
1854
1855 do n = 1, ob % num_satem
1856 do k = 1, ob % satem(n) % numlevs
1857 call da_calc_trace_single( ob % satem(n) % thickness(k), count1, trace1 )
1858 end do
1859 end do
1860
1861 ob % jo_satem_thickness = 0.5 * ( count1 - trace1 )
1862 ob % trace_satem = trace1
1863
1864 if ( ob % trace_satem< 0.0 ) &
1865 write(6,'(a,f15.5)') &
1866 ' Warning: ob % trace_satem < 0 = ', ob % trace_satem
1867
1868 ob % trace_total = ob % trace_total + ob % trace_satem
1869 ob % num_satem_tot = count1
1870 end if
1871
1872 ob % num_ssmt1_tot = 0
1873 ob % trace_ssmt1 = 0
1874 if ( ob % num_ssmt1 > 0 ) then
1875 trace1 = 0.0
1876 count1 = 0
1877
1878 do n = 1, ob % num_ssmt1
1879 do k = 1, ob % ssmt1(n) % numlevs
1880 call da_calc_trace_single( ob % ssmt1(n) % t(k), count1, trace1 )
1881 end do
1882 end do
1883
1884 ob % jo_ssmt1_t = 0.5 * ( count1 - trace1 )
1885 ob % trace_ssmt1 = trace1
1886
1887 if ( ob % trace_ssmt1< 0.0 ) &
1888 write(6,'(a,f15.5)') &
1889 ' Warning: ob % trace_ssmt1 < 0 = ', ob % trace_ssmt1
1890
1891 ob % trace_total = ob % trace_total + ob % trace_ssmt1
1892 ob % num_ssmt1_tot = count1
1893 end if
1894
1895 ob % num_ssmt2_tot = 0
1896 ob % trace_ssmt2 = 0
1897 if ( ob % num_ssmt2 > 0 ) then
1898 trace1 = 0.0
1899 count1 = 0
1900
1901 do n = 1, ob % num_ssmt2
1902 do k = 1, ob % ssmt2(n) % numlevs
1903 call da_calc_trace_single( ob % ssmt2(n) % rh(k), count1, trace1 )
1904 end do
1905 end do
1906
1907 ob % jo_ssmt2_rh = 0.5 * ( count1 - trace1 )
1908 ob % trace_ssmt2 = trace1
1909
1910 if ( ob % trace_ssmt2< 0.0 ) &
1911 write(6,'(a,f15.5)') &
1912 ' Warning: ob % trace_ssmt2 < 0 = ', ob % trace_ssmt2
1913
1914 ob % trace_total = ob % trace_total + ob % trace_ssmt2
1915 ob % num_ssmt2_tot = count1
1916 end if
1917
1918 ob % num_buoy_tot = 0
1919 ob % trace_buoy = 0
1920 if ( ob % num_buoy > 0 ) then
1921 trace1 = 0.0; trace2 = 0.0; trace3 = 0.0; trace4 = 0.0; trace5 = 0.0
1922 count1 = 0; count2 = 0; count3 = 0; count4 = 0; count5 = 0
1923
1924 do n = 1, ob % num_buoy
1925 call da_calc_trace_single( ob % buoy (n) % u, count1, trace1 )
1926 call da_calc_trace_single( ob % buoy (n) % v, count2, trace2 )
1927 call da_calc_trace_single( ob % buoy (n) % t, count3, trace3 )
1928 call da_calc_trace_single( ob % buoy (n) % p, count4, trace4 )
1929 call da_calc_trace_single( ob % buoy (n) % q, count5, trace5 )
1930 end do
1931
1932 ob % jo_buoy_u = 0.5 * ( count1 - trace1 )
1933 ob % jo_buoy_v = 0.5 * ( count2 - trace2 )
1934 ob % jo_buoy_t = 0.5 * ( count3 - trace3 )
1935 ob % jo_buoy_p = 0.5 * ( count4 - trace4 )
1936 ob % jo_buoy_q = 0.5 * ( count5 - trace5 )
1937 ob % trace_buoy = trace1 + trace2 + trace3 + trace4 + trace5
1938
1939 if ( ob % trace_buoy < 0.0 ) &
1940 write(6,'(a,f15.5)') &
1941 ' Warning: ob % trace_buoy < 0 = ', ob % trace_buoy
1942
1943 ob % trace_total = ob % trace_total + ob % trace_buoy
1944 ob % num_buoy_tot = count1 + count2 + count3 + count4 + count5
1945 end if
1946 ob % num_sonde_sfc_tot = 0
1947 ob % trace_sonde_sfc = 0
1948 if ( ob % num_sonde_sfc > 0 ) then
1949 trace1 = 0.0; trace2 = 0.0; trace3 = 0.0; trace4 = 0.0; trace5 = 0.0
1950 count1 = 0; count2 = 0; count3 = 0; count4 = 0; count5 = 0
1951
1952 do n = 1, ob % num_sonde_sfc
1953 call da_calc_trace_single( ob % sonde_sfc (n) % u, count1, trace1 )
1954 call da_calc_trace_single( ob % sonde_sfc (n) % v, count2, trace2 )
1955 call da_calc_trace_single( ob % sonde_sfc (n) % t, count3, trace3 )
1956 call da_calc_trace_single( ob % sonde_sfc (n) % p, count4, trace4 )
1957 call da_calc_trace_single( ob % sonde_sfc (n) % q, count5, trace5 )
1958 end do
1959
1960 ob % jo_sonde_sfc_u = 0.5 * ( count1 - trace1 )
1961 ob % jo_sonde_sfc_v = 0.5 * ( count2 - trace2 )
1962 ob % jo_sonde_sfc_t = 0.5 * ( count3 - trace3 )
1963 ob % jo_sonde_sfc_p = 0.5 * ( count4 - trace4 )
1964 ob % jo_sonde_sfc_q = 0.5 * ( count5 - trace5 )
1965 ob % trace_sonde_sfc = trace1 + trace2 + trace3 + trace4 + trace5
1966 if ( ob % trace_sonde_sfc < 0.0 ) &
1967 write(6,'(a,f15.5)') &
1968 ' Warning: ob % trace_sonde_sfc < 0 = ', ob % trace_sonde_sfc
1969
1970 ob % trace_total = ob % trace_total + ob % trace_sonde_sfc
1971 ob % num_sonde_sfc_tot = count1 + count2 + count3 + count4 + count5
1972 end if
1973 ob % num_profiler_tot = 0
1974 ob % trace_profiler = 0
1975 if ( ob % num_profiler > 0 ) then
1976 trace1 = 0.0; trace2 = 0.0
1977 count1 = 0; count2 = 0
1978
1979 do n = 1, ob % num_profiler
1980 do k = 1, ob % profiler(n) % numlevs
1981 call da_calc_trace_single( ob % profiler(n) % u(k), count1, trace1 )
1982 call da_calc_trace_single( ob % profiler(n) % v(k), count2, trace2 )
1983 end do
1984 end do
1985
1986 ob % jo_profiler_u = 0.5 * ( count1 - trace1 )
1987 ob % jo_profiler_v = 0.5 * ( count2 - trace2 )
1988 ob % trace_profiler = trace1 + trace2
1989
1990 if ( ob % trace_profiler < 0.0 ) &
1991 write(6,'(a,f15.5)') &
1992 ' Warning: ob % trace_profiler < 0 = ', ob % trace_profiler
1993
1994 ob % trace_total = ob % trace_total + ob % trace_profiler
1995 ob % num_profiler_tot = count1 + count2
1996 end if
1997
1998 ob % num_qscat_tot = 0
1999 ob % trace_qscat = 0
2000 if ( ob % num_qscat > 0 ) then
2001 trace1 = 0.0; trace2 = 0.0
2002 count1 = 0; count2 = 0
2003
2004 do n = 1, ob % num_qscat
2005 call da_calc_trace_single( ob % qscat(n) % u, count1, trace1 )
2006 call da_calc_trace_single( ob % qscat(n) % v, count2, trace2 )
2007 end do
2008
2009 ob % jo_qscat_u = 0.5 * ( count1 - trace1 )
2010 ob % jo_qscat_v = 0.5 * ( count2 - trace2 )
2011 ob % trace_qscat = trace1 + trace2
2012
2013 if ( ob % trace_qscat < 0.0 ) &
2014 write(6,'(a,f15.5)') &
2015 ' Warning: ob % trace_qscat < 0 = ', ob % trace_qscat
2016
2017 ob % trace_total = ob % trace_total + ob % trace_qscat
2018 ob % num_qscat_tot = count1 + count2
2019
2020 end if
2021
2022 ob % num_bogus_tot = 0
2023 ob % trace_bogus = 0
2024 if ( ob % num_bogus > 0 ) then
2025 trace1 = 0.0; trace2 = 0.0; trace3 = 0.0; trace4 = 0.0; trace5 = 0.0
2026 count1 = 0; count2 = 0; count3 = 0; count4 = 0; count5 = 0
2027
2028 do n = 1, ob % num_bogus
2029 do k = 1, ob % bogus(n) % numlevs
2030 call da_calc_trace_single( ob % bogus(n) % u(k), count1, trace1 )
2031 call da_calc_trace_single( ob % bogus(n) % v(k), count2, trace2 )
2032 call da_calc_trace_single( ob % bogus(n) % t(k), count3, trace3 )
2033 call da_calc_trace_single( ob % bogus(n) % q(k), count4, trace4 )
2034 end do
2035 call da_calc_trace_single( ob % bogus(n) % slp, count5, trace5 )
2036 end do
2037
2038 ob % jo_bogus_u = 0.5 * ( count1 - trace1 )
2039 ob % jo_bogus_v = 0.5 * ( count2 - trace2 )
2040 ob % jo_bogus_t = 0.5 * ( count3 - trace3 )
2041 ob % jo_bogus_q = 0.5 * ( count4 - trace4 )
2042 ob % jo_bogus_slp = 0.5 * ( count5 - trace5 )
2043 ob % trace_bogus = trace1 + trace2 + trace3 + trace4 + trace5
2044
2045 if ( ob % trace_bogus < 0.0 ) &
2046 write(6,'(a,f15.5)') &
2047 ' Warning: ob % trace_bogus < 0 = ', ob % trace_bogus
2048
2049 ob % trace_total = ob % trace_total + ob % trace_bogus
2050 ob % num_bogus_tot = count1 + count2 + count3 + count4 + count5
2051 end if
2052
2053 ob % num_gpsref_tot = 0
2054 ob % trace_gpsref = 0
2055 if ( ob % num_gpsref > 0 ) then
2056 trace1 = 0.0
2057 count1 = 0
2058
2059 do n = 1, ob % num_gpsref
2060 do k = 1, ob % gpsref(n) % numlevs
2061 call da_calc_trace_single( ob % gpsref(n) % ref(k), count1, trace1 )
2062 end do
2063 end do
2064
2065 ob % jo_gpsref_ref = 0.5 * ( count1 - trace1 )
2066 ob % trace_gpsref = trace1
2067
2068 if ( ob % trace_gpsref < 0.0 ) &
2069 write(6,'(a,f15.5)') &
2070 ' Warning: ob % trace_gpsref < 0 = ', ob % trace_gpsref
2071
2072 ob % trace_total = ob % trace_total + ob % trace_gpsref
2073 ob % num_gpsref_tot = count1
2074 end if
2075
2076 ob % total_obs = ob % num_synop_tot + ob % num_metar_tot + &
2077 ob % num_ships_tot + ob % num_polaramv_tot + ob % num_geoamv_tot + &
2078 ob % num_gpspw_tot + ob % num_sound_tot + ob % num_airep_tot + &
2079 ob % num_pilot_tot + ob % num_ssmir_tot + ob % num_satem_tot + &
2080 ob % num_ssmt1_tot + ob % num_ssmt2_tot + ob % num_buoy_tot + &
2081 ob % num_sonde_sfc_tot + ob % num_qscat_tot + ob % num_profiler_tot + &
2082 ob% num_bogus_tot + ob% num_gpsref_tot
2083
2084 ! radiance part
2085 if ( rtminit_nsensor > 0 ) then
2086
2087 do n = 1, rtminit_nsensor
2088 do ichan = 1, ob % rad(n) % nchan
2089 if ( ob % rad(n) % num_rad_tot(ichan) > 0 ) then
2090 trace1 = 0.0
2091 count1 = 0
2092 do k=1,ob % rad(n) % num_rad_tot(ichan)
2093 call da_calc_trace_single( &
2094 ob % rad(n) % tb(ichan)%pixel(k), count1, trace1 )
2095 end do
2096 ob % rad(n) % jo_rad(ichan) = 0.5 * ( count1 - trace1 )
2097 ob % rad(n) % trace_rad(ichan) = trace1
2098 if ( ob % rad(n) % trace_rad(ichan) < 0.0 ) &
2099 write(6,'(a,i3,a,f15.5)') &
2100 ' Warning: '//trim(ob%rad(n)%rttovid_string), ichan, &
2101 ' Trace(HK) < 0 = ', ob%rad(n)%trace_rad(ichan)
2102 ob % trace_total = ob % trace_total + ob%rad(n)%trace_rad(ichan)
2103 ob % total_obs = ob % total_obs + ob%rad(n)%num_rad_tot(ichan)
2104 end if
2105 end do
2106 end do
2107 end if
2108
2109 end subroutine da_calc_jo_expected
2110
2111
2112
2113 subroutine da_calc_trace_single( field, count, trace )
2114
2115 implicit none
2116
2117 type (field_type), intent(in) :: field
2118 integer, intent(inout) :: count
2119 real, intent(inout) :: trace
2120
2121 if ( field % yp /= missing_r .and. field % y /= missing_r ) then
2122 count = count + 1
2123 trace = trace + ( field % yp - field % y ) * field % pert / field % error
2124 end if
2125
2126 end subroutine da_calc_trace_single
2127
2128
2129
2130 subroutine da_read_jo_actual( ob )
2131
2132 implicit none
2133
2134 type (iv_type), intent(inout) :: ob
2135
2136 character (len=46) :: str
2137 character (len=15) :: str1, str2, str3, str4, str5
2138 character (len=15) :: str6, str7, str8, str9, str10
2139 character (len=5) :: ob_name
2140 real :: dum1, dum2, dum3, dum4, dum5, jo
2141 integer :: n, num_obs
2142
2143 ob % joa_synop_u = 0.0
2144 ob % joa_synop_v = 0.0
2145 ob % joa_synop_t = 0.0
2146 ob % joa_synop_p = 0.0
2147 ob % joa_synop_q = 0.0
2148 ob % joa_metar_u = 0.0
2149 ob % joa_metar_v = 0.0
2150 ob % joa_metar_t = 0.0
2151 ob % joa_metar_p = 0.0
2152 ob % joa_metar_q = 0.0
2153 ob % joa_ships_u = 0.0
2154 ob % joa_ships_v = 0.0
2155 ob % joa_ships_t = 0.0
2156 ob % joa_ships_p = 0.0
2157 ob % joa_ships_q = 0.0
2158 ob % joa_polaramv_u = 0.0
2159 ob % joa_polaramv_v = 0.0
2160 ob % joa_geoamv_u = 0.0
2161 ob % joa_geoamv_v = 0.0
2162 ob % joa_gpspw_tpw = 0.0
2163 ob % joa_sound_u = 0.0
2164 ob % joa_sound_v = 0.0
2165 ob % joa_sound_t = 0.0
2166 ob % joa_sound_q = 0.0
2167 ob % joa_airep_u = 0.0
2168 ob % joa_airep_v = 0.0
2169 ob % joa_airep_t = 0.0
2170 ob % joa_pilot_u = 0.0
2171 ob % joa_pilot_v = 0.0
2172 ob % joa_ssmir_speed = 0.0
2173 ob % joa_ssmir_tpw = 0.0
2174 ob % joa_satem_thickness = 0.0
2175 ob % joa_ssmt1_t = 0.0
2176 ob % joa_ssmt2_rh = 0.0
2177 ob % joa_buoy_u = 0.0
2178 ob % joa_buoy_v = 0.0
2179 ob % joa_buoy_t = 0.0
2180 ob % joa_buoy_p = 0.0
2181 ob % joa_buoy_q = 0.0
2182 ob % joa_sonde_sfc_u = 0.0
2183 ob % joa_sonde_sfc_v = 0.0
2184 ob % joa_sonde_sfc_t = 0.0
2185 ob % joa_sonde_sfc_p = 0.0
2186 ob % joa_sonde_sfc_q = 0.0
2187 ob % joa_qscat_u = 0.0
2188 ob % joa_qscat_v = 0.0
2189 ob % joa_profiler_u = 0.0
2190 ob % joa_profiler_v = 0.0
2191 ob % joa_bogus_u = 0.0
2192 ob % joa_bogus_v = 0.0
2193 ob % joa_bogus_t = 0.0
2194 ob % joa_bogus_q = 0.0
2195 ob % joa_bogus_slp = 0.0
2196 ob % joa_airsr_t = 0.0
2197 ob % joa_airsr_q = 0.0
2198 ob % joa_gpsref_ref = 0.0
2199
2200 if ( rtminit_nsensor > 0 ) then
2201 do n = 1,rtminit_nsensor
2202 ob%rad(n)%num_rad_tot(:) = 0
2203 ob%rad(n)%joa_rad(:) = 0.0
2204 end do
2205 end if
2206
2207 rewind(jo_unit)
2208
2209 do
2210 read(jo_unit,'(a46,10a15)')str, str1, str2, str3, str4, str5, &
2211 str6, str7, str8, str9, str10
2212 ob_name = adjustl(str)
2213
2214 if ( ob_name == 'synop' .and. index(str,'actual') > 0 ) then
2215 call da_read_jo_actual1( str1, str2, str3, str4, str5, &
2216 str6, str7, str8, str9, str10, &
2217 ob % joa_synop_u, ob % joa_synop_v, &
2218 ob % joa_synop_t, ob % joa_synop_p, &
2219 ob % joa_synop_q )
2220
2221 else if ( ob_name == 'metar' .and. index(str,'actual') > 0 ) then
2222 call da_read_jo_actual1( str1, str2, str3, str4, str5, &
2223 str6, str7, str8, str9, str10, &
2224 ob % joa_metar_u, ob % joa_metar_v, &
2225 ob % joa_metar_t, ob % joa_metar_p, &
2226 ob % joa_metar_q )
2227
2228 else if ( ob_name == 'ships' .and. index(str,'actual') > 0 ) then
2229 call da_read_jo_actual1( str1, str2, str3, str4, str5, &
2230 str6, str7, str8, str9, str10, &
2231 ob % joa_ships_u, ob % joa_ships_v, &
2232 ob % joa_ships_t, ob % joa_ships_p, &
2233 ob % joa_ships_q )
2234
2235 else if ( ob_name == 'polar' .and. index(str,'actual') > 0 ) then
2236 call da_read_jo_actual1( str1, str2, str3, str4, str5, &
2237 str6, str7, str8, str9, str10, &
2238 ob % joa_polaramv_u, ob % joa_polaramv_v, &
2239 dum3, dum4, dum5 )
2240
2241 else if ( ob_name == 'geoam' .and. index(str,'actual') > 0 ) then
2242 call da_read_jo_actual1( str1, str2, str3, str4, str5, &
2243 str6, str7, str8, str9, str10, &
2244 ob % joa_geoamv_u, ob % joa_geoamv_v, &
2245 dum3, dum4, dum5 )
2246
2247 else if ( ob_name == 'gpspw' .and. index(str,'actual') > 0 ) then
2248 call da_read_jo_actual1( str1, str2, str3, str4, str5, &
2249 str6, str7, str8, str9, str10, &
2250 ob % joa_gpspw_tpw, dum2, dum3, dum4, dum5 )
2251
2252 else if ( ob_name == 'sound' .and. index(str,'actual') > 0 ) then
2253 call da_read_jo_actual1( str1, str2, str3, str4, str5, &
2254 str6, str7, str8, str9, str10, &
2255 ob % joa_sound_u, ob % joa_sound_v, ob % joa_sound_t, &
2256 ob % joa_sound_q, dum1 )
2257
2258 else if ( ob_name == 'airsr' .and. index(str,'actual') > 0 ) then
2259 call da_read_jo_actual1( str1, str2, str3, str4, str5, &
2260 str6, str7, str8, str9, str10, &
2261 ob % joa_airsr_t, ob % joa_airsr_q,&
2262 dum1, dum2, dum3 )
2263
2264 else if ( ob_name == 'airep' .and. index(str,'actual') > 0 ) then
2265 call da_read_jo_actual1( str1, str2, str3, str4, str5, &
2266 str6, str7, str8, str9, str10, &
2267 ob % joa_airep_u, ob % joa_airep_v, &
2268 ob % joa_airep_t, dum1, dum2 )
2269
2270 else if ( ob_name == 'pilot' .and. index(str,'actual') > 0 ) then
2271 call da_read_jo_actual1( str1, str2, str3, str4, str5, &
2272 str6, str7, str8, str9, str10, &
2273 ob % joa_pilot_u, ob % joa_pilot_v, &
2274 dum1, dum2, dum3 )
2275
2276 else if ( ob_name == 'ssmir' .and. index(str,'actual') > 0 ) then
2277 call da_read_jo_actual1( str1, str2, str3, str4, str5, &
2278 str6, str7, str8, str9, str10, &
2279 ob % joa_ssmir_speed, ob % joa_ssmir_tpw, &
2280 dum1, dum2, dum3 )
2281
2282 else if ( ob_name == 'satem' .and. index(str,'actual') > 0 ) then
2283 call da_read_jo_actual1( str1, str2, str3, str4, str5, &
2284 str6, str7, str8, str9, str10, &
2285 ob % joa_satem_thickness, dum1, dum2, &
2286 dum3, dum4 )
2287
2288 else if ( ob_name == 'ssmt1' .and. index(str,'actual') > 0 ) then
2289 call da_read_jo_actual1( str1, str2, str3, str4, str5, &
2290 str6, str7, str8, str9, str10, &
2291 ob % joa_ssmt1_t, dum1, dum2, dum3, dum4 )
2292
2293 else if ( ob_name == 'ssmt2' .and. index(str,'actual') > 0 ) then
2294 call da_read_jo_actual1( str1, str2, str3, str4, str5, &
2295 str6, str7, str8, str9, str10, &
2296 ob % joa_ssmt2_rh, dum1, dum2, dum3, dum4 )
2297
2298 else if ( ob_name == 'buoy ' .and. index(str,'actual') > 0 ) then
2299 call da_read_jo_actual1( str1, str2, str3, str4, str5, &
2300 str6, str7, str8, str9, str10, &
2301 ob % joa_buoy_u, ob % joa_buoy_v, &
2302 ob % joa_buoy_t, ob % joa_buoy_p, &
2303 ob % joa_buoy_q )
2304 else if ( ob_name == 'sonde' .and. index(str,'actual') > 0 ) then
2305 call da_read_jo_actual1( str1, str2, str3, str4, str5, &
2306 str6, str7, str8, str9, str10, &
2307 ob % joa_sonde_sfc_u, ob % joa_sonde_sfc_v, &
2308 ob % joa_sonde_sfc_t, ob % joa_sonde_sfc_p, &
2309 ob % joa_sonde_sfc_q )
2310
2311 else if ( ob_name == 'profi' .and. index(str,'actual') > 0 ) then
2312 call da_read_jo_actual1( str1, str2, str3, str4, str5, &
2313 str6, str7, str8, str9, str10, &
2314 ob % joa_profiler_u, ob % joa_profiler_v, &
2315 dum1, dum2, dum3 )
2316
2317 else if ( ob_name == 'qscat' .and. index(str,'actual') > 0 ) then
2318 call da_read_jo_actual1( str1, str2, str3, str4, str5, &
2319 str6, str7, str8, str9, str10, &
2320 ob % joa_qscat_u, ob % joa_qscat_v, &
2321 dum1, dum2, dum3 )
2322
2323 else if ( ob_name == 'bogus' .and. index(str,'actual') > 0 ) then
2324 call da_read_jo_actual1( str1, str2, str3, str4, str5, &
2325 str6, str7, str8, str9, str10, &
2326 ob % joa_bogus_u, ob % joa_bogus_v, &
2327 ob % joa_bogus_t, ob % joa_bogus_q, &
2328 ob % joa_bogus_slp)
2329 else if ( ob_name == 'gpsre' .and. index(str,'actual') > 0 ) then
2330 call da_read_jo_actual1( str1, str2, str3, str4, str5, &
2331 str6, str7, str8, str9, str10, &
2332 ob % joa_gpsref_ref, &
2333 dum1, dum2, dum3, dum4 )
2334 else if ( ob_name == 'radia' .and. index(str,'actual') > 0 ) then
2335 platform_id = 1
2336 read (str,'(30x,a4,1x,i2,1x,a5)')platform, satellite_id,sensor
2337 read (str1,'(i5,i10)')ichan, num_obs
2338 read (str2,'(f15.5)')jo
2339
2340 if ( sensor == 'amsua' ) then
2341 sensor_id = 3
2342 else if ( sensor == 'amsub' ) then
2343 sensor_id = 4
2344 else
2345 write(6,*) ' Unrecognized Sensor '
2346 end if
2347 do n = 1, rtminit_nsensor
2348 if ( platform_id == rtminit_platform(n) &
2349 .and. satellite_id == rtminit_satid(n) &
2350 .and. sensor_id == rtminit_sensor(n) ) then
2351 ob%rad(n)%joa_rad(ichan) = ob%rad(n)%joa_rad(ichan) + jo
2352 ob%rad(n)%num_rad_tot(ichan) = ob%rad(n)%num_rad_tot(ichan) + &
2353 num_obs
2354 exit
2355 end if
2356 end do
2357 else if ( str(1:5) == '*****' ) then
2358 exit
2359 end if
2360
2361 end do
2362
2363 end subroutine da_read_jo_actual
2364
2365
2366
2367 subroutine da_read_jo_actual1( str1, str2, str3, str4, str5, &
2368 str6, str7, str8, str9, str10, &
2369 j1, j2, j3, j4, j5 )
2370
2371 implicit none
2372
2373 character (len=15), intent(in) :: str1, str2, str3, str4, str5
2374 character (len=15), intent(in) :: str6, str7, str8, str9, str10
2375 real, intent(inout) :: j1, j2, j3, j4, j5
2376
2377 real :: j1n, j2n, j3n, j4n, j5n
2378 real :: f1n, f2n, f3n, f4n, f5n
2379
2380 read(str1,*)j1n
2381 read(str2,*)f1n
2382 read(str3,*)j2n
2383 read(str4,*)f2n
2384 read(str5,*)j3n
2385 read(str6,*)f3n
2386 read(str7,*)j4n
2387 read(str8,*)f4n
2388 read(str9,*)j5n
2389 read(str10,*)f5n
2390
2391 ! Jo(actual) is scaled by the square of the input error factor used.
2392
2393 j1 = j1 + f1n**2 * j1n
2394 j2 = j2 + f2n**2 * j2n
2395 j3 = j3 + f3n**2 * j3n
2396 j4 = j4 + f4n**2 * j4n
2397 j5 = j5 + f5n**2 * j5n
2398
2399 end subroutine da_read_jo_actual1
2400
2401
2402
2403 subroutine da_calc_new_factors( ob )
2404
2405 implicit none
2406
2407 type (iv_type), intent(inout) :: ob
2408 integer :: n, ichan
2409
2410 write(6,*)
2411
2412 if ( ob % num_synop > 0 ) then
2413 call da_calc_new_factors1( 'synop', ob % num_synop, ob % num_synop_tot, &
2414 ob % jo_synop_u, ob % jo_synop_v, ob % jo_synop_t, ob % jo_synop_p, &
2415 ob % jo_synop_q, &
2416 ob % joa_synop_u, ob % joa_synop_v, ob % joa_synop_t, &
2417 ob % joa_synop_p, ob % joa_synop_q )
2418 end if
2419
2420 if ( ob % num_metar > 0 ) then
2421 call da_calc_new_factors1( 'metar', ob % num_metar, ob % num_metar_tot, &
2422 ob % jo_metar_u, ob % jo_metar_v, ob % jo_metar_t, ob % jo_metar_p, &
2423 ob % jo_metar_q, &
2424 ob % joa_metar_u, ob % joa_metar_v, ob % joa_metar_t, &
2425 ob % joa_metar_p, ob % joa_metar_q )
2426 end if
2427
2428 if ( ob % num_ships > 0 ) then
2429 call da_calc_new_factors1( 'ships', ob % num_ships, ob % num_ships_tot, &
2430 ob % jo_ships_u, ob % jo_ships_v, ob % jo_ships_t, ob % jo_ships_p, &
2431 ob % jo_ships_q, &
2432 ob % joa_ships_u, ob % joa_ships_v, ob % joa_ships_t, &
2433 ob % joa_ships_p, ob % joa_ships_q )
2434 end if
2435
2436 if ( ob % num_geoamv > 0 ) then
2437 call da_calc_new_factors1( 'geoamv', ob % num_geoamv, &
2438 ob % num_geoamv_tot, &
2439 ob % jo_geoamv_u, ob % jo_geoamv_v, 0.0, 0.0, 0.0, &
2440 ob % joa_geoamv_u, ob % joa_geoamv_v, 0.0, 0.0, 0.0 )
2441 end if
2442
2443 if ( ob % num_polaramv > 0 ) then
2444 call da_calc_new_factors1( 'polaramv', ob % num_polaramv, &
2445 ob % num_polaramv_tot, &
2446 ob % jo_polaramv_u, ob % jo_polaramv_v, 0.0, 0.0, 0.0, &
2447 ob % joa_polaramv_u, ob % joa_polaramv_v, 0.0, 0.0, 0.0 )
2448 end if
2449
2450 if ( ob % num_gpspw > 0 ) then
2451 call da_calc_new_factors1( 'gpspw', ob % num_gpspw, ob % num_gpspw_tot, &
2452 ob % jo_gpspw_tpw, 0.0, 0.0, 0.0, 0.0, &
2453 ob % joa_gpspw_tpw, 0.0, 0.0, 0.0, 0.0 )
2454 end if
2455
2456 if ( ob % num_sound > 0 ) then
2457 call da_calc_new_factors1( 'sound', ob % num_sound, ob % num_sound_tot, &
2458 ob % jo_sound_u, ob % jo_sound_v, ob % jo_sound_t, ob % jo_sound_q, &
2459 0.0, &
2460 ob % joa_sound_u, ob % joa_sound_v, ob % joa_sound_t, &
2461 ob % joa_sound_q, 0.0)
2462 end if
2463
2464 if ( ob % num_airep > 0 ) then
2465 call da_calc_new_factors1( 'airep', ob % num_airep, ob % num_airep_tot, &
2466 ob % jo_airep_u, ob % jo_airep_v, ob % jo_airep_t, 0.0, 0.0, &
2467 ob % joa_airep_u, ob % joa_airep_v, ob % joa_airep_t, 0.0, 0.0 )
2468 end if
2469
2470 if ( ob % num_airsr > 0 ) then
2471 call da_calc_new_factors1( 'airsr', ob % num_airsr, ob % num_airsr_tot, &
2472 ob % jo_airsr_t, ob % jo_airsr_q, 0.0, 0.0, 0.0, &
2473 ob % joa_airsr_t, ob % joa_airsr_q, 0.0, 0.0, 0.0 )
2474 end if
2475
2476 if ( ob % num_pilot > 0 ) then
2477 call da_calc_new_factors1( 'pilot', ob % num_pilot, ob % num_pilot_tot, &
2478 ob % jo_pilot_u, ob % jo_pilot_v, 0.0, 0.0, 0.0, &
2479 ob % joa_pilot_u, ob % joa_pilot_v, 0.0, 0.0, 0.0 )
2480 end if
2481
2482 if ( ob % num_ssmir > 0 ) then
2483 call da_calc_new_factors1( 'ssmir', ob % num_ssmir, ob % num_ssmir_tot, &
2484 ob % jo_ssmir_speed, ob % jo_ssmir_tpw, 0.0, 0.0, 0.0, &
2485 ob % joa_ssmir_speed, ob % joa_ssmir_tpw, 0.0, 0.0, 0.0 )
2486 end if
2487
2488 if ( ob % num_satem > 0 ) then
2489 call da_calc_new_factors1( 'satem', ob % num_satem, ob % num_satem_tot, &
2490 ob % jo_satem_thickness, 0.0, 0.0, 0.0, 0.0, &
2491 ob % joa_satem_thickness, 0.0, 0.0, 0.0, 0.0 )
2492 end if
2493
2494 if ( ob % num_ssmt1 > 0 ) then
2495 call da_calc_new_factors1( 'ssmt1', ob % num_ssmt1, ob % num_ssmt1_tot, &
2496 ob % jo_ssmt1_t, 0.0, 0.0, 0.0, 0.0, &
2497 ob % joa_ssmt1_t, 0.0, 0.0, 0.0, 0.0 )
2498 end if
2499
2500 if ( ob % num_ssmt2 > 0 ) then
2501 call da_calc_new_factors1( 'ssmt2', ob % num_ssmt2, ob % num_ssmt2_tot, &
2502 ob % jo_ssmt2_rh, 0.0, 0.0, 0.0, 0.0, &
2503 ob % joa_ssmt2_rh, 0.0, 0.0, 0.0, 0.0 )
2504 end if
2505
2506 if ( ob % num_buoy > 0 ) then
2507 call da_calc_new_factors1( 'buoy ', ob % num_buoy, ob % num_buoy_tot, &
2508 ob % jo_buoy_u, ob % jo_buoy_v, ob % jo_buoy_t, ob % jo_buoy_p, &
2509 ob % jo_buoy_q, &
2510 ob % joa_buoy_u, ob % joa_buoy_v, ob % joa_buoy_t, ob % joa_buoy_p, &
2511 ob % joa_buoy_q )
2512 end if
2513
2514 if ( ob % num_sonde_sfc > 0 ) then
2515 call da_calc_new_factors1( 'sonde', ob % num_sonde_sfc, &
2516 ob % num_sonde_sfc_tot, &
2517 ob % jo_sonde_sfc_u, ob % jo_sonde_sfc_v, ob % jo_sonde_sfc_t, &
2518 ob % jo_sonde_sfc_p, ob % jo_sonde_sfc_q, &
2519 ob % joa_sonde_sfc_u, ob % joa_sonde_sfc_v, ob % joa_sonde_sfc_t, &
2520 ob % joa_sonde_sfc_p, ob % joa_sonde_sfc_q )
2521 end if
2522
2523 if ( ob % num_profiler > 0 ) then
2524 call da_calc_new_factors1( 'profi', ob % num_profiler, &
2525 ob % num_profiler_tot, &
2526 ob % jo_profiler_u, ob % jo_profiler_v, 0.0, 0.0, 0.0, &
2527 ob % joa_profiler_u, ob % joa_profiler_v, 0.0, 0.0, 0.0 )
2528 end if
2529
2530 if ( ob % num_qscat > 0 ) then
2531 call da_calc_new_factors1( 'qscat', ob % num_qscat, ob % num_qscat_tot, &
2532 ob % jo_qscat_u, ob % jo_qscat_v, 0.0, 0.0, 0.0, &
2533 ob % joa_qscat_u, ob % joa_qscat_v, 0.0, 0.0, 0.0 )
2534 end if
2535
2536 if ( ob % num_bogus > 0 ) then
2537 call da_calc_new_factors1( 'bogus', ob % num_bogus, ob % num_bogus_tot, &
2538 ob % jo_bogus_u, ob % jo_bogus_v, ob % jo_bogus_t, ob % jo_bogus_q, &
2539 ob % jo_bogus_slp, &
2540 ob % joa_bogus_u, ob % joa_bogus_v, ob % joa_bogus_t, &
2541 ob % joa_bogus_q, ob % joa_bogus_slp )
2542 end if
2543
2544 if ( ob % num_gpsref > 0 ) then
2545 call da_calc_new_factors1( 'gpsre', ob % num_gpsref, ob % num_gpsref_tot, &
2546 ob % jo_gpsref_ref, 0.0, 0.0, 0.0, 0.0, &
2547 ob % joa_gpsref_ref, 0.0, 0.0, 0.0, 0.0 )
2548 end if
2549
2550 if ( rtminit_nsensor > 0 ) then
2551 write(6,*) &
2552 ' sensor chan num Jo_mini Jo_exp trace(HK) factor'
2553 do n = 1, rtminit_nsensor
2554 do ichan = 1, ob % rad(n) % nchan
2555 if ( ob % rad(n) % num_rad_tot(ichan) > 0 ) then
2556 ob % rad(n) % factor_rad(ichan) = &
2557 sqrt( ob % rad(n) % joa_rad(ichan) / ob % rad(n) % jo_rad(ichan) )
2558 write(6,'(a15,i3,i8,3f15.5,f8.3)') &
2559 trim(ob%rad(n)%rttovid_string), &
2560 ichan, &
2561 ob%rad(n)%num_rad_tot(ichan), &
2562 ob%rad(n)%joa_rad(ichan), &
2563 ob%rad(n)%jo_rad(ichan), &
2564 ob%rad(n)%trace_rad(ichan), &
2565 ob%rad(n)%factor_rad(ichan)
2566 end if
2567 end do
2568 end do
2569 end if
2570
2571 end subroutine da_calc_new_factors
2572
2573
2574
2575 subroutine da_calc_new_factors1( ob_name, ob_num, ob_num_tot, &
2576 j1e, j2e, j3e, j4e, j5e, &
2577 j1a, j2a, j3a, j4a, j5a )
2578
2579 implicit none
2580
2581 character (len=5), intent(in) :: ob_name
2582 integer, intent(in) :: ob_num, ob_num_tot
2583 real, intent(in) :: j1e, j2e, j3e, j4e, j5e
2584 real, intent(in) :: j1a, j2a, j3a, j4a, j5a
2585
2586 real :: f1, f2, f3, f4, f5
2587
2588 f1 = 1.0; f2 = 1.0; f3 = 1.0; f4 = 1.0; f5 = 1.0
2589
2590 if ( j1e > 0.0 ) f1 = sqrt( j1a / j1e )
2591 if ( j2e > 0.0 ) f2 = sqrt( j2a / j2e )
2592 if ( j3e > 0.0 ) f3 = sqrt( j3a / j3e )
2593 if ( j4e > 0.0 ) f4 = sqrt( j4a / j4e )
2594 if ( j5e > 0.0 ) f5 = sqrt( j5a / j5e )
2595
2596 write(6,'(1x,a5,a21,2i8,6f15.5)')ob_name, ' obs, Jo (expected)= ', &
2597 ob_num, ob_num_tot, &
2598 j1e, j2e, j3e, j4e, j5e
2599
2600 write(6,'(1x,a5,a21,2i8,6f15.5)')ob_name, ' obs, Jo (actual) = ', &
2601 ob_num, ob_num_tot, &
2602 j1a, j2a, j3a, j4a, j5a
2603
2604 write(6,'(1x,a5,a21,2i8,6f15.5)')ob_name, ' obs, Error Factor = ', &
2605 ob_num, ob_num_tot, f1, f2, f3, f4, f5
2606 write(6,*)
2607
2608 end subroutine da_calc_new_factors1
2609
2610
2611
2612 subroutine da_get_j( ob )
2613
2614 implicit none
2615
2616 type (iv_type), intent(in) :: ob
2617
2618 character (len=80) :: str
2619 integer :: icount
2620 real :: j, jo, jb, jn, jon, jbn, j_e, jo_e, jb_e
2621 real :: jb_factor_old, jb_factor_oldn, jb_factor_new
2622
2623 rewind(in_unit)
2624
2625 j = 0.0; jo = 0.0; jb = 0.0; j_e = 0.0; jo_e = 0.0; jb_e = 0.0
2626 icount = 0
2627 jb_factor_old = 0
2628 jb_factor_new = 0
2629
2630 do
2631 read(in_unit,'(a80)')str
2632 if ( index(str,'Final value of J ') > 0 ) then
2633 read(str(index(str,'=')+1:80),*)jn
2634 j = j + jn
2635 else if ( index(str,'Final value of Jo') > 0 ) then
2636 read(str(index(str,'=')+1:80),*)jon
2637 jo = jo + jon
2638 else if ( index(str,'Final value of Jb') > 0 ) then
2639 read(str(index(str,'=')+1:80),*)jbn
2640 jb = jb + jbn
2641 else if ( index(str,'Jb factor used(1)') > 0 ) then
2642 read(str(index(str,'=')+1:80),*)jb_factor_oldn
2643 jb_factor_old = jb_factor_old + jb_factor_oldn
2644 icount = icount + 1
2645 else if ( str(1:5) == '*****' ) then
2646 exit
2647 end if
2648 end do
2649
2650 jb_factor_old = jb_factor_old / real(icount)
2651
2652 write(6,'(/a,i8)') ' Total number of obs. = ', ob % total_obs
2653 write(6,'(/a,i8)') ' Total number of cases = ', icount
2654 j_e = 0.5 * ob % total_obs
2655 jo_e = 0.5 * ( ob % total_obs - ob % trace_total )
2656 jb_e = 0.5 * ob % trace_total
2657
2658 write(6,'(a,f15.5)')' Total J (actual) = ', j
2659 write(6,'(a,f15.5)')' Total J (expected) = ', j_e
2660 write(6,*)
2661
2662 write(6,'(a,f15.5)')' Total Jo(actual) = ', jo
2663 write(6,'(a,f15.5)')' Total Jo(expected) = ', jo_e
2664 write(6,'(a,f15.5)')' Total Jo factor = ', sqrt(jo/jo_e)
2665
2666 write(6,*)
2667 write(6,'(a,f15.5)')' Total Jb(actual) = ', jb
2668 write(6,'(a,f15.5)')' Total Jb(expected) = ', jb_e
2669 write(6,'(a,f15.5)')' Total Jb factor (old) = ', jb_factor_old
2670
2671 if ( jb_e < 0.0 ) then
2672 write(6,'(a)')' Warning: Tr(HK) < 0.0 Too small a sample?'
2673 stop
2674 end if
2675 jb_factor_new = sqrt(jb/jb_e) * jb_factor_old
2676 write(6,'(a,f15.5)')' Total Jb factor (new) = ', jb_factor_new
2677 write(6,*)
2678
2679 end subroutine da_get_j
2680
2681
2682 subroutine read_namelist_radiance
2683
2684 implicit none
2685
2686 character(len=filename_len) :: namelist_file ! Input namelist filename.
2687 integer, parameter :: namelist_unit = 7 ! Input namelist unit.
2688 integer :: iost ! Error code.
2689
2690 namelist /rtminit/ rtminit_nsensor, rtminit_platform, &
2691 rtminit_satid, rtminit_sensor, rtminit_nchan
2692
2693 namelist_file = 'namelist.radiance'
2694 write (6, '(3x,a,a)' ) ' radiance namelist file : ', namelist_file
2695 iost = 0
2696
2697 open ( file = trim(namelist_file), unit = namelist_unit, &
2698 status = 'old' , access = 'sequential', &
2699 form = 'formatted', action = 'read', &
2700 iostat = iost )
2701 if ( iost /= 0 ) stop ' error in opening namelist file '
2702
2703 iost = 0
2704
2705 read ( unit = namelist_unit, nml = rtminit, iostat = iost)
2706 write(6,'(a,i4 )') ' rtminit_nsensor = ', rtminit_nsensor
2707 write(6,'(a,10i4)') ' rtminit_platform = ', rtminit_platform(1:rtminit_nsensor)
2708 write(6,'(a,10i4)') ' rtminit_satid = ', rtminit_satid (1:rtminit_nsensor)
2709 write(6,'(a,10i4)') ' rtminit_sensor = ', rtminit_sensor (1:rtminit_nsensor)
2710 write(6,'(a,10i4)') ' rtminit_nchan = ', rtminit_nchan (1:rtminit_nsensor)
2711
2712 if ( iost /= 0 ) stop ' error in reading namelist file '
2713
2714 close (namelist_unit)
2715
2716 if ( rtminit_nsensor > 0 ) then
2717 allocate ( ob % rad(rtminit_nsensor) )
2718 do n = 1,rtminit_nsensor
2719 ob % rad(n) % nchan = rtminit_nchan(n)
2720 allocate ( ob % rad(n) % num_rad_tot(ob % rad(n) % nchan) )
2721 allocate ( ob % rad(n) % jo_rad(ob % rad(n) % nchan) )
2722 allocate ( ob % rad(n) % joa_rad(ob % rad(n) % nchan) )
2723 allocate ( ob % rad(n) % factor_rad(ob % rad(n) % nchan) )
2724 allocate ( ob % rad(n) % trace_rad(ob % rad(n) % nchan) )
2725 allocate ( ob % rad(n) % tb(ob % rad(n) % nchan) )
2726 ob % rad(n) % num_rad_tot(:) = 0
2727 write(ob%rad(n)%rttovid_string, '(a,i2.2,a)') &
2728 trim( platform_name(rtminit_platform(n)) )//'-', &
2729 rtminit_satid(n), &
2730 '-'//trim( inst_name(rtminit_sensor(n)) )
2731 write(6,*) 'Tuning radiance error for ', trim(ob%rad(n)%rttovid_string)
2732 end do
2733 end if
2734
2735 end subroutine read_namelist_radiance
2736
2737 end program da_tune_obs_desroziers
2738
2739