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