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