da_get_var_diagnostics.inc

References to this file elsewhere.
1 subroutine da_get_var_diagnostics( iv, j)
2 
3    !-----------------------------------------------------------------------
4    ! Purpose: TBD
5    !-----------------------------------------------------------------------
6 
7    implicit none
8 
9    type(iv_type), intent(inout):: iv      ! innovation vector.
10    type(j_type), intent(inout) :: j       ! Cost function.
11 
12    integer                      :: num_stats_tot
13    integer                      :: i,k
14    real                         :: jo_radiance
15    real                         :: temp(63)
16 
17    if (trace_use) call da_trace_entry("da_get_var_diagnostics")
18 
19    !--------------------------------------------------------------------------
20    ! [1.0] Sum up Jo across processors:
21    !--------------------------------------------------------------------------
22 
23    num_stats_tot = sum(iv%nstats(:))
24 
25    temp(1)  = j % jo % synop_u
26    temp(2)  = j % jo % synop_v
27    temp(3)  = j % jo % synop_t
28    temp(4)  = j % jo % synop_p
29    temp(5)  = j % jo % synop_q
30    temp(6)  = j % jo % metar_u
31    temp(7)  = j % jo % metar_v
32    temp(8)  = j % jo % metar_t
33    temp(9)  = j % jo % metar_p
34    temp(10) = j % jo % metar_q
35    temp(11) = j % jo % ships_u
36    temp(12) = j % jo % ships_v
37    temp(13) = j % jo % ships_t
38    temp(14) = j % jo % ships_p
39    temp(15) = j % jo % ships_q
40    temp(16) = j % jo % geoamv_u
41    temp(17) = j % jo % geoamv_v
42    temp(18) = j % jo % polaramv_u
43    temp(19) = j % jo % polaramv_v      
44    temp(20) = j % jo % gpspw_tpw       
45    temp(21) = j % jo % gpsref_ref      
46    temp(22) = j % jo % sound_u         
47    temp(23) = j % jo % sound_v         
48    temp(24) = j % jo % sound_t         
49    temp(25) = j % jo % sound_q         
50    temp(26) = j % jo % sonde_sfc_u     
51    temp(27) = j % jo % sonde_sfc_v     
52    temp(28) = j % jo % sonde_sfc_t     
53    temp(29) = j % jo % sonde_sfc_p     
54    temp(30) = j % jo % sonde_sfc_q     
55    temp(31) = j % jo % airep_u         
56    temp(32) = j % jo % airep_v         
57    temp(33) = j % jo % airep_t         
58    temp(34) = j % jo % pilot_u         
59    temp(35) = j % jo % pilot_v         
60    temp(36) = j % jo % bogus_u         
61    temp(37) = j % jo % bogus_v         
62    temp(38) = j % jo % bogus_t         
63    temp(39) = j % jo % bogus_q         
64    temp(40) = j % jo % bogus_slp       
65    temp(41) = j % jo % ssmir_speed     
66    temp(42) = j % jo % ssmir_tpw       
67    temp(43) = j % jo % ssmi_tb19v      
68    temp(44) = j % jo % ssmi_tb19h      
69    temp(45) = j % jo % ssmi_tb22v      
70    temp(46) = j % jo % ssmi_tb37v      
71    temp(47) = j % jo % ssmi_tb37h      
72    temp(48) = j % jo % ssmi_tb85v      
73    temp(49) = j % jo % ssmi_tb85h      
74    temp(50) = j % jo % satem_thickness 
75    temp(51) = j % jo % ssmt1_t         
76    temp(52) = j % jo % ssmt2_rh        
77    temp(53) = j % jo % qscat_u         
78    temp(54) = j % jo % qscat_v         
79    temp(55) = j % jo % profiler_u      
80    temp(56) = j % jo % profiler_v      
81    temp(57) = j % jo % buoy_u          
82    temp(58) = j % jo % buoy_v          
83    temp(59) = j % jo % buoy_t          
84    temp(60) = j % jo % buoy_p          
85    temp(61) = j % jo % buoy_q          
86    temp(62) = j % jo % airsr_t         
87    temp(63) = j % jo % airsr_q         
88 
89    call da_proc_sum_real(temp(:))
90 
91    j % jo % synop_u         = temp(1)  
92    j % jo % synop_v         = temp(2)  
93    j % jo % synop_t         = temp(3)  
94    j % jo % synop_p         = temp(4)  
95    j % jo % synop_q         = temp(5)  
96    j % jo % metar_u         = temp(6)  
97    j % jo % metar_v         = temp(7)  
98    j % jo % metar_t         = temp(8)  
99    j % jo % metar_p         = temp(9)  
100    j % jo % metar_q         = temp(10) 
101    j % jo % ships_u         = temp(11) 
102    j % jo % ships_v         = temp(12) 
103    j % jo % ships_t         = temp(13) 
104    j % jo % ships_p         = temp(14) 
105    j % jo % ships_q         = temp(15) 
106    j % jo % geoamv_u        = temp(16) 
107    j % jo % geoamv_v        = temp(17) 
108    j % jo % polaramv_u      = temp(18) 
109    j % jo % polaramv_v      = temp(19) 
110    j % jo % gpspw_tpw       = temp(20) 
111    j % jo % gpsref_ref      = temp(21) 
112    j % jo % sound_u         = temp(22) 
113    j % jo % sound_v         = temp(23) 
114    j % jo % sound_t         = temp(24) 
115    j % jo % sound_q         = temp(25) 
116    j % jo % sonde_sfc_u     = temp(26) 
117    j % jo % sonde_sfc_v     = temp(27) 
118    j % jo % sonde_sfc_t     = temp(28) 
119    j % jo % sonde_sfc_p     = temp(29) 
120    j % jo % sonde_sfc_q     = temp(30) 
121    j % jo % airep_u         = temp(31) 
122    j % jo % airep_v         = temp(32) 
123    j % jo % airep_t         = temp(33) 
124    j % jo % pilot_u         = temp(34) 
125    j % jo % pilot_v         = temp(35) 
126    j % jo % bogus_u         = temp(36) 
127    j % jo % bogus_v         = temp(37) 
128    j % jo % bogus_t         = temp(38) 
129    j % jo % bogus_q         = temp(39) 
130    j % jo % bogus_slp       = temp(40) 
131    j % jo % ssmir_speed     = temp(41) 
132    j % jo % ssmir_tpw       = temp(42) 
133    j % jo % ssmi_tb19v      = temp(43) 
134    j % jo % ssmi_tb19h      = temp(44) 
135    j % jo % ssmi_tb22v      = temp(45) 
136    j % jo % ssmi_tb37v      = temp(46) 
137    j % jo % ssmi_tb37h      = temp(47) 
138    j % jo % ssmi_tb85v      = temp(48) 
139    j % jo % ssmi_tb85h      = temp(49) 
140    j % jo % satem_thickness = temp(50) 
141    j % jo % ssmt1_t         = temp(51) 
142    j % jo % ssmt2_rh        = temp(52) 
143    j % jo % qscat_u         = temp(53) 
144    j % jo % qscat_v         = temp(54) 
145    j % jo % profiler_u      = temp(55) 
146    j % jo % profiler_v      = temp(56) 
147    j % jo % buoy_u          = temp(57) 
148    j % jo % buoy_v          = temp(58) 
149    j % jo % buoy_t          = temp(59) 
150    j % jo % buoy_p          = temp(60) 
151    j % jo % buoy_q          = temp(61) 
152    j % jo % airsr_t         = temp(62) 
153    j % jo % airsr_q         = temp(63) 
154 
155    if (use_rad) then
156       jo_radiance = 0.0
157       do i = 1, iv%num_inst                 ! loop for sensor
158          call da_proc_sum_ints(j % jo % rad(i)% num_ichan(:))
159          call da_proc_sum_real(j % jo % rad(i) % jo_ichan(:))
160          jo_radiance = jo_radiance + sum(j % jo % rad(i) % jo_ichan(:))
161       end do
162    end if
163 
164    !-----------------------------------------------------------------------------
165    ! [2.0] Print out VAR diagnostics:
166    !-----------------------------------------------------------------------------
167 
168    if (rootproc) then
169 
170       write(unit=stdout,fmt=*) ' '
171       write(unit=stdout,fmt='(A)') 'Diagnostics'
172       write(unit=stdout,fmt='(A,F12.2)')   '   Final cost function J       = ', j % total
173       write(unit=stdout,fmt=*) ' '
174 
175       write(unit=stdout,fmt='(a,i8)')    '   Total number of obs.        = ', num_stats_tot
176       write(unit=stdout,fmt='(a,f15.5)') '   Final value of J            = ', j % total
177       write(unit=stdout,fmt='(a,f15.5)') '   Final value of Jo           = ', j % jo % total
178       write(unit=stdout,fmt='(a,f15.5)') '   Final value of Jb           = ', j % jb
179       write(unit=stdout,fmt='(a,f15.5)') '   Final value of Jc           = ', j % jc
180       write(unit=stdout,fmt='(a,f15.5)') '   Final value of Je           = ', j % je
181       if (num_stats_tot > 0) &
182          write(unit=stdout,fmt='(a,f15.5)') '   Final J / total num_obs     = ', j % total / &
183                                                           real(num_stats_tot)
184       write(unit=stdout,fmt='(a,f15.5)') '   Jb factor used(1)           = ', var_scaling1
185       write(unit=stdout,fmt='(a,f15.5)') '   Jb factor used(2)           = ', var_scaling2
186       write(unit=stdout,fmt='(a,f15.5)') '   Jb factor used(3)           = ', var_scaling3
187       write(unit=stdout,fmt='(a,f15.5)') '   Jb factor used(4)           = ', var_scaling4
188       write(unit=stdout,fmt='(a,f15.5)') '   Jb factor used(5)           = ', var_scaling5
189       write(unit=stdout,fmt='(a,f15.5)') '   Jb factor used              = ', jb_factor
190       write(unit=stdout,fmt='(a,f15.5)') '   Je factor used              = ', je_factor
191       write(unit=stdout,fmt=*) ' '
192 
193       if (use_rad) then
194          write(unit=stdout,fmt='(a,i8)')    '   Total number of radiances    = ', iv%nstats(radiance)
195          write(unit=stdout,fmt='(a,f15.5)') '   Cost function for radiances  = ', jo_radiance
196          write(unit=stdout,fmt=*) ' '
197       end if
198 
199       ! [4.2] Output components of Jo:
200 
201       if (iv%info(synop)%ntotal > 0) then
202          write(unit=jo_unit,fmt='(a30,2i8,10f15.5)')'    synop obs, Jo(actual)  = ', &
203                                   iv%info(synop)%ntotal, iv%nstats(synop), &
204                                   j % jo % synop_u, iv % synop_ef_u, &
205                                   j % jo % synop_v, iv % synop_ef_v, &
206                                   j % jo % synop_t, iv % synop_ef_t, &
207                                   j % jo % synop_p, iv % synop_ef_p, &
208                                   j % jo % synop_q, iv % synop_ef_q
209 
210       end if
211 
212       if (trace_use) call da_trace("da_get_var_diagnostics", &
213          message="Memory increase from internal write")
214 
215       if (iv%info(metar)%ntotal > 0) then
216          write(unit=jo_unit,fmt='(a30,2i8,10f15.5)')'    metar obs, Jo(actual)  = ', &
217                                iv%info(metar)%ntotal, iv%nstats(metar), &
218                                j % jo % metar_u, iv % metar_ef_u, &
219                                j % jo % metar_v, iv % metar_ef_v, &
220                                j % jo % metar_t, iv % metar_ef_t, &
221                                j % jo % metar_p, iv % metar_ef_p, &
222                                j % jo % metar_q, iv % metar_ef_q    
223       end if
224 
225       if (iv%info(ships)%ntotal > 0) then
226          write(unit=jo_unit,fmt='(a30,2i8,10f15.5)')'    ships obs, Jo(actual)  = ', &
227                                iv%info(ships)%ntotal, iv%nstats(ships), &
228                                j % jo % ships_u, iv % ships_ef_u, &
229                                j % jo % ships_v, iv % ships_ef_v, &
230                                j % jo % ships_t, iv % ships_ef_t, &
231                                j % jo % ships_p, iv % ships_ef_p, &
232                                j % jo % ships_q, iv % ships_ef_q                                
233       end if
234 
235 
236       if (iv%info(geoamv)%ntotal > 0) then
237          write(unit=jo_unit,fmt='(a30,2i8,10f15.5)')'    geoamv ob, Jo(actual)  = ', &
238                                iv%info(geoamv)%ntotal, iv%nstats(geoamv), &
239                                j % jo % geoamv_u, iv % geoamv_ef_u, &
240                                j % jo % geoamv_v, iv % geoamv_ef_v, &
241                                0.0, 1.0, 0.0, 1.0, 0.0, 1.0
242       end if
243 
244       if (iv%info(polaramv)%ntotal > 0) then
245          write(unit=jo_unit,fmt='(a30,2i8,10f15.5)')'    polaramv,  Jo(actual)  = ', &
246                                iv%info(polaramv)%ntotal, iv%nstats(polaramv), &
247                                j % jo % polaramv_u, iv % polaramv_ef_u, &
248                                j % jo % polaramv_v, iv % polaramv_ef_v, &
249                                0.0, 1.0, 0.0, 1.0, 0.0, 1.0
250       end if
251 
252 
253       if (iv%info(gpspw)%ntotal > 0) then
254          write(unit=jo_unit,fmt='(a30,2i8,10f15.5)')'    gpspw obs, Jo(actual)  = ', &
255                                iv%info(gpspw)%ntotal, iv%nstats(gpspw), &
256                                j % jo % gpspw_tpw, iv % gpspw_ef_tpw, &
257                                0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0
258       end if
259 
260       if (iv%info(gpsref)%ntotal > 0) then
261          write(unit=jo_unit,fmt='(a30,2i8,10f15.5)')'    gpsref obs, Jo(actual)  = ', &
262                                iv%info(gpsref)%ntotal, iv%nstats(gpsref), &
263                                j % jo % gpsref_ref, iv % gpsref_ef_ref, &
264                                0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0
265       end if
266 
267       if (iv%info(sound)%ntotal > 0) then
268          write(unit=jo_unit,fmt='(a30,2i8,10f15.5)')'    sound obs, Jo(actual)  = ', &
269                                iv%info(sound)%ntotal, iv%nstats(sound), &
270                                j % jo % sound_u, iv % sound_ef_u, &
271                                j % jo % sound_v, iv % sound_ef_v, &
272                                j % jo % sound_t, iv % sound_ef_t, &
273                                j % jo % sound_q, iv % sound_ef_q, 0.0, 1.0 
274          write(unit=jo_unit,fmt='(a30,2i8,10f15.5)')'    sonde obs, Jo(actual)  = ', &
275                                iv%info(sound)%ntotal, iv%nstats(sound), &
276                                j % jo % sonde_sfc_u, iv % synop_ef_u, &
277                                j % jo % sonde_sfc_v, iv % synop_ef_v, &
278                                j % jo % sonde_sfc_t, iv % synop_ef_t, &
279                                j % jo % sonde_sfc_p, iv % synop_ef_p, &
280                                j % jo % sonde_sfc_q, iv % synop_ef_q
281       end if
282 
283       if (iv%info(airep)%ntotal > 0) then
284          write(unit=jo_unit,fmt='(a30,2i8,10f15.5)')'    airep obs, Jo(actual)  = ', &
285                                iv%info(airep)%ntotal, iv%nstats(airep), &
286                                j % jo % airep_u, iv % airep_ef_u, &
287                                j % jo % airep_v, iv % airep_ef_v, &
288                                j % jo % airep_t, iv % airep_ef_t, &
289                                0.0, 1.0, 0.0, 1.0
290       end if
291 
292       if (iv%info(bogus)%ntotal > 0) then
293          write(unit=jo_unit,fmt='(a30,2i8,10f15.5)')'    bogus obs, Jo(actual)  = ', &
294                                iv%info(bogus)%ntotal, iv%nstats(bogus), &
295                                j % jo % bogus_u, iv % bogus_ef_u, &
296                                j % jo % bogus_v, iv % bogus_ef_v, &
297                                j % jo % bogus_t, iv % bogus_ef_t, &
298                                j % jo % bogus_q, iv % bogus_ef_q, &
299                                j % jo % bogus_slp, iv % bogus_ef_slp
300       end if
301 
302       if (iv%info(pilot)%ntotal > 0) then
303          write(unit=jo_unit,fmt='(a30,2i8,10f15.5)')'    pilot obs, Jo(actual)  = ', &
304                                iv%info(pilot)%ntotal, iv%nstats(pilot), &
305                                j % jo % pilot_u, iv % pilot_ef_u, &
306                                j % jo % pilot_v, iv % pilot_ef_v, &
307                                0.0, 1.0, 0.0, 1.0, 0.0, 1.0
308       end if
309 
310       if (iv%info(ssmi_rv)%ntotal > 0) then
311          write(unit=jo_unit,fmt='(a30,2i8,10f15.5)')'    ssmir obs, Jo(actual) = ', &
312                                   iv%info(ssmi_rv)%ntotal, iv%nstats(ssmi_rv), &
313                                   j % jo % ssmir_speed, iv % ssmir_ef_speed, &
314                                   j % jo % ssmir_tpw, iv % ssmir_ef_tpw, &
315                                   0.0, 1.0, 0.0, 1.0, 0.0, 1.0
316       end if
317 
318       if (iv%info(satem)%ntotal > 0) then
319          write(unit=jo_unit,fmt='(a30,2i8,10f15.5)')'    satem obs, Jo(actual)  = ', &
320                                iv%info(satem)%ntotal, iv%nstats(satem), &
321                                j % jo % satem_thickness, iv % satem_ef_thickness, &
322                                0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0
323       end if
324 
325       if (iv%info(ssmt1)%ntotal > 0) then
326          write(unit=jo_unit,fmt='(a30,2i8,10f15.5)')'    ssmt1 obs, Jo(actual)  = ', &
327                                iv%info(ssmt1)%ntotal, iv%nstats(ssmt1), &
328                                j % jo % ssmt1_t, iv % ssmt1_ef_t, &
329                                0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0
330       end if
331 
332       if (iv%info(ssmt2)%ntotal > 0) then
333          write(unit=jo_unit,fmt='(a30,2i8,10f15.5)')'    ssmt2 obs, Jo(actual)  = ', &
334                                iv%info(ssmt2)%ntotal, iv%nstats(ssmt2), &
335                                j % jo % ssmt2_rh, iv % ssmt2_ef_rh, &
336                                0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0
337       end if
338 
339       if (iv%info(qscat)%ntotal > 0) then
340          write(unit=jo_unit,fmt='(a30,2i8,10f15.5)')'    qscat obs, Jo(actual)  = ', &
341                                iv%info(qscat)%ntotal, iv%nstats(qscat), &
342                                j % jo % qscat_u, iv % qscat_ef_u, &
343                                j % jo % qscat_v, iv % qscat_ef_v, &
344                                0.0, 1.0, 0.0, 1.0, 0.0, 1.0
345       end if
346 
347       if (iv%info(buoy)%ntotal > 0) then
348          write(unit=jo_unit,fmt='(a30,2i8,10f15.5)')'    buoy  obs, Jo(actual)  = ', &
349                                iv%info(buoy)%ntotal, iv%nstats(buoy), &
350                                j % jo % buoy_u, iv % buoy_ef_u, &
351                                j % jo % buoy_v, iv % buoy_ef_v, &
352                                j % jo % buoy_t, iv % buoy_ef_t, &
353                                j % jo % buoy_p, iv % buoy_ef_p, &
354                                j % jo % buoy_q, iv % buoy_ef_q                                
355       end if
356 
357       if (iv%info(profiler)%ntotal > 0) then
358          write(unit=jo_unit,fmt='(a30,2i8,10f15.5)')'    profiler,  Jo(actual)  = ', &
359                                iv%info(profiler)%ntotal, iv%nstats(profiler), &
360                                j % jo % profiler_u, iv % profiler_ef_u, &
361                                j % jo % profiler_v, iv % profiler_ef_v, &
362                                0.0, 1.0, 0.0, 1.0, 0.0, 1.0
363       end if
364       if (iv%info(airsr)%ntotal > 0) then
365          write(unit=jo_unit,fmt='(a30,2i8,10f15.5)')'    airsr obs, Jo(actual)  = ', &
366                                iv%info(airsr)%ntotal, iv%nstats(airsr), &
367                                j % jo % airsr_t, iv % airsr_ef_t, &
368                                j % jo % airsr_q, iv % airsr_ef_q, &
369                                0.0, 1.0, 0.0, 1.0, 0.0, 1.0
370       end if
371       do i = 1, iv%num_inst                 ! loop for sensor
372          do k = 1, iv%instid(i)%nchan
373             if (j % jo % rad(i) % num_ichan(k) > 0) then
374                write(unit=jo_unit,fmt='(a30,a16,i5,i10,8f15.5)')'    radiance,  Jo(actual)  = ', &
375                   iv%instid(i)%rttovid_string, k , &
376                   j % jo % rad(i) % num_ichan(k), &
377                   j % jo % rad(i) % jo_ichan(k)
378             end if
379          end do
380       end do
381    end if
382 
383    if (trace_use) call da_trace_exit("da_get_var_diagnostics")
384       
385 end subroutine da_get_var_diagnostics
386 
387