ROMS
Loading...
Searching...
No Matches
Functions/Subroutines
stats_modobs_mod Module Reference

Functions/Subroutines

subroutine, public stats_modobs (ng, tile)
 
subroutine, private stats_modobs_nf90 (ng, tile, lbi, ubi, lbj, ubj)
 
subroutine, private stats_modobs_pio (ng, tile, lbi, ubi, lbj, ubj)
 

Function/Subroutine Documentation

◆ stats_modobs()

subroutine, public stats_modobs_mod::stats_modobs ( integer, intent(in) ng,
integer, intent(in) tile )

Definition at line 62 of file stats_modobs.F.

63!***********************************************************************
64!
65! Imported variable declarations.
66!
67 integer, intent(in) :: ng, tile
68!
69! Local variable declarations.
70!
71 integer :: LBi, UBi, LBj, UBj
72!
73 character (len=*), parameter :: MyFile = &
74 & __FILE__
75!
76!-----------------------------------------------------------------------
77! Compute and write out model-observation statistics.
78!-----------------------------------------------------------------------
79!
80! Skip if configuration error.
81!
82 IF ((exit_flag.eq.5).or.(exit_flag.eq.6)) RETURN
83!
84! Process model-observation statistics
85!
86 lbi=bounds(ng)%LBi(tile)
87 ubi=bounds(ng)%UBi(tile)
88 lbj=bounds(ng)%LBj(tile)
89 ubj=bounds(ng)%UBj(tile)
90!
91 SELECT CASE (dav(ng)%IOtype)
92 CASE (io_nf90)
93 CALL stats_modobs_nf90 (ng, tile, &
94 & lbi, ubi, lbj, ubj)
95
96# if defined PIO_LIB && defined DISTRIBUTE
97 CASE (io_pio)
98 CALL stats_modobs_pio (ng, tile, &
99 & lbi, ubi, lbj, ubj)
100# endif
101 CASE DEFAULT
102 IF (master) WRITE (stdout,10) dav(ng)%IOtype
103 exit_flag=3
104 END SELECT
105 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
106!
107 10 FORMAT (' STATS_MODOBS - Illegal output file type, io_type = ', &
108 & i0,/,16x,'Check KeyWord ''OUT_LIB'' in ''roms.in''.')
109!
110 RETURN

References mod_param::bounds, mod_iounits::dav, mod_scalars::exit_flag, strings_mod::founderror(), mod_ncparam::io_nf90, mod_ncparam::io_pio, mod_parallel::master, mod_scalars::noerror, stats_modobs_nf90(), stats_modobs_pio(), and mod_iounits::stdout.

Referenced by roms_kernel_mod::roms_run().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ stats_modobs_nf90()

subroutine, private stats_modobs_mod::stats_modobs_nf90 ( integer, intent(in) ng,
integer, intent(in) tile,
integer, intent(in) lbi,
integer, intent(in) ubi,
integer, intent(in) lbj,
integer, intent(in) ubj )
private

Definition at line 114 of file stats_modobs.F.

116!***********************************************************************
117!
118! Imported variable declarations.
119!
120 integer, intent(in) :: ng, tile
121 integer, intent(in) :: LBi, UBi, LBj, UBj
122!
123! Local variable declarations.
124!
125 integer :: i, ic, iobs, status, vindex
126 integer :: saved_exit_flag
127
128 integer, dimension(2) :: start, total
129 integer, dimension(NobsVar(ng)) :: Ncount, is
130
131 integer, dimension(Ndatum(ng)) :: obs_type
132 integer, dimension(Ndatum(ng)) :: provenance
133 integer, dimension(Nsurvey(ng)) :: survey_Nobs
134!
135 real(r8) :: cff1, cff2
136 real(r8), parameter :: IniVal = 0.0_r8
137
138 real(r8), dimension(NobsVar(ng)) :: CC, MB, MSE, SDE
139 real(r8), dimension(NobsVar(ng)) :: mod_min, mod_max
140 real(r8), dimension(NobsVar(ng)) :: mod_mean, mod_std
141 real(r8), dimension(NobsVar(ng)) :: obs_min, obs_max
142 real(r8), dimension(NobsVar(ng)) :: obs_mean, obs_std
143
144 real(r8), dimension(Ndatum(ng)) :: mod_value
145 real(r8), dimension(Ndatum(ng)) :: obs_depths
146 real(r8), dimension(Ndatum(ng)) :: obs_scale
147 real(dp), dimension(Ndatum(ng)) :: obs_time
148 real(r8), dimension(Ndatum(ng)) :: obs_value
149 real(r8), dimension(Ndatum(ng)) :: obs_Xgrid
150 real(r8), dimension(Ndatum(ng)) :: obs_Ygrid
151 real(r8), dimension(Ndatum(ng)) :: obs_work
152 real(dp), dimension(Nsurvey(ng)) :: survey
153!
154 character (len=11), dimension(NobsVar(ng)) :: text, svar_name
155
156 character (len=30) :: F1, F2, F3, F4
157
158 character (len=*), parameter :: MyFile = &
159 & __FILE__//", stats_modobs_nf90"
160!
161 sourcefile=myfile
162!
163!-----------------------------------------------------------------------
164! Read in model and observations data.
165!-----------------------------------------------------------------------
166!
167! If exit_flag > 0, save and reset to allow completion of the DAV
168! when ROMS blows-up.
169!
170 IF ((exit_flag.eq.1).or.(blowup.gt.0)) THEN
171 saved_exit_flag=exit_flag;
172 exit_flag=noerror
173 ELSE
174 saved_exit_flag=noerror;
175 END IF
176!
177! Inquire about input observations variables.
178!
179 CALL netcdf_inq_var (ng, inlm, obs(ng)%name, &
180 & ncid = obs(ng)%ncid)
181 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
182!
183! Read in and write out number of observations per survey.
184!
185 CALL netcdf_get_ivar (ng, inlm, obs(ng)%name, &
186 & vname(1,idnobs), survey_nobs, &
187 & ncid = obs(ng)%ncid, &
188 & start = (/1/), &
189 & total = (/nsurvey(ng)/))
190 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
191
192 CALL netcdf_put_ivar (ng, inlm, dav(ng)%name, &
193 & vname(1,idnobs), survey_nobs, &
194 & (/1/), (/nsurvey(ng)/), &
195 & ncid = dav(ng)%ncid)
196 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
197!
198! Read in and write out survey time.
199!
200 CALL netcdf_get_time (ng, inlm, obs(ng)%name, &
201 & vname(1,idoday), &
202 & rclock%DateNumber, survey, &
203 & ncid = obs(ng)%ncid, &
204 & start = (/1/), &
205 & total = (/nsurvey(ng)/))
206 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
207
208 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
209 & vname(1,idoday), survey, &
210 & (/1/), (/nsurvey(ng)/), &
211 & ncid = dav(ng)%ncid)
212 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
213!
214! Read in and write out observation type identifier.
215!
216 CALL netcdf_get_ivar (ng, inlm, obs(ng)%name, &
217 & vname(1,idotyp), obs_type, &
218 & ncid = obs(ng)%ncid, &
219 & start = (/1/), &
220 & total = (/ndatum(ng)/))
221 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
222
223 CALL netcdf_put_ivar (ng, inlm, dav(ng)%name, &
224 & vname(1,idotyp), obs_type, &
225 & (/1/), (/ndatum(ng)/), &
226 & ncid = dav(ng)%ncid)
227 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
228!
229! Read in and write out observation provenance.
230!
231 CALL netcdf_get_ivar (ng, inlm, obs(ng)%name, &
232 & vname(1,idopro), provenance, &
233 & ncid = obs(ng)%ncid, &
234 & start = (/1/), &
235 & total = (/ndatum(ng)/))
236 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
237
238 CALL netcdf_put_ivar (ng, inlm, dav(ng)%name, &
239 & vname(1,idopro), provenance, &
240 & (/1/), (/ndatum(ng)/), &
241 & ncid = dav(ng)%ncid)
242 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
243!
244! Read in and write out observation time.
245!
246 CALL netcdf_get_time (ng, inlm, obs(ng)%name, &
247 & vname(1,idobst), &
248 & rclock%DateNumber, obs_time, &
249 & ncid = obs(ng)%ncid, &
250 & start = (/1/), &
251 & total = (/ndatum(ng)/))
252 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
253
254 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
255 & vname(1,idobst), obs_time, &
256 & (/1/), (/ndatum(ng)/), &
257 & ncid = dav(ng)%ncid)
258 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
259!
260! Read in and write out observation longitude.
261!
262 IF (find_string(var_name,n_var,vname(1,idolon),vindex)) THEN
263 CALL netcdf_get_fvar (ng, inlm, obs(ng)%name, &
264 & vname(1,idolon), obs_work, &
265 & ncid = obs(ng)%ncid, &
266 & start = (/1/), &
267 & total = (/ndatum(ng)/))
268 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
269
270 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
271 & vname(1,idolon), obs_work, &
272 & (/1/), (/ndatum(ng)/), &
273 & ncid = dav(ng)%ncid)
274 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
275 END IF
276!
277! Read in and write out observation latitude.
278!
279 IF (find_string(var_name,n_var,vname(1,idolon),vindex)) THEN
280 CALL netcdf_get_fvar (ng, inlm, obs(ng)%name, &
281 & vname(1,idolat), obs_work, &
282 & ncid = obs(ng)%ncid, &
283 & start = (/1/), &
284 & total = (/ndatum(ng)/))
285 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
286
287 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
288 & vname(1,idolat), obs_work, &
289 & (/1/), (/ndatum(ng)/), &
290 & ncid = dav(ng)%ncid)
291 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
292 END IF
293!
294! Read in observation depth.
295!
296 CALL netcdf_get_fvar (ng, inlm, obs(ng)%name, &
297 & vname(1,idobsd), obs_depths, &
298 & ncid = obs(ng)%ncid, &
299 & start = (/1/), &
300 & total = (/ndatum(ng)/))
301 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
302!
303! Read in and write out X-fractional coordinate.
304!
305 CALL netcdf_get_fvar (ng, inlm, obs(ng)%name, &
306 & vname(1,idobsx), obs_xgrid, &
307 & ncid = obs(ng)%ncid, &
308 & start = (/1/), &
309 & total = (/ndatum(ng)/))
310 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
311
312 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
313 & vname(1,idobsx), obs_xgrid, &
314 & (/1/), (/ndatum(ng)/), &
315 & ncid = dav(ng)%ncid)
316 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
317!
318! Read in and write out Y-fractional coordinate.
319!
320 CALL netcdf_get_fvar (ng, inlm, obs(ng)%name, &
321 & vname(1,idobsy), obs_ygrid, &
322 & ncid = obs(ng)%ncid, &
323 & start = (/1/), &
324 & total = (/ndatum(ng)/))
325 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
326
327 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
328 & vname(1,idobsy), obs_ygrid, &
329 & (/1/), (/ndatum(ng)/), &
330 & ncid = dav(ng)%ncid)
331 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
332!
333! Read in and write out observations total error covariance
334! (instrument + sampling + representation).
335!
336 CALL netcdf_get_fvar (ng, inlm, obs(ng)%name, &
337 & vname(1,idoerr), obs_work, &
338 & ncid = obs(ng)%ncid, &
339 & start = (/1/), &
340 & total = (/ndatum(ng)/))
341 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
342
343 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
344 & vname(1,idoerr), obs_work, &
345 & (/1/), (/ndatum(ng)/), &
346 & ncid = dav(ng)%ncid)
347 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
348!
349! Read in and write out observation values.
350!
351 CALL netcdf_get_fvar (ng, inlm, obs(ng)%name, &
352 & vname(1,idoval), obs_value, &
353 & ncid = obs(ng)%ncid, &
354 & start = (/1/), &
355 & total = (/ndatum(ng)/))
356 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
357
358 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
359 & vname(1,idoval), obs_value, &
360 & (/1/), (/ndatum(ng)/), &
361 & ncid = dav(ng)%ncid)
362 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
363!
364! Read in and write out observation meta values.
365!
366 IF (haveobsmeta(ng)) THEN
367 CALL netcdf_get_fvar (ng, inlm, obs(ng)%name, &
368 & vname(1,idomet), obs_work, &
369 & ncid = obs(ng)%ncid, &
370 & start = (/1/), &
371 & total = (/ndatum(ng)/))
372 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
373
374 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
375 & vname(1,idomet), obs_work, &
376 & (/1/), (/ndatum(ng)/), &
377 & ncid = dav(ng)%ncid)
378 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
379 END IF
380!
381! Read in observation screening flag.
382!
383 CALL netcdf_get_fvar (ng, inlm, dav(ng)%name, &
384 & vname(1,idobss), obs_scale, &
385 & ncid = dav(ng)%ncid, &
386 & start = (/1/), &
387 & total = (/ndatum(ng)/))
388 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
389!
390! Read in model values at observation locations.
391!
392# if defined R4DVAR || defined R4DVAR_ANA_SENSITIVITY || \
393 defined tl_r4dvar
394 CALL netcdf_get_fvar (ng, inlm, dav(ng)%name, &
395 & vname(1,idtlmo), mod_value, &
396 & ncid = dav(ng)%ncid, &
397 & start = (/1/), &
398 & total = (/ndatum(ng)/))
399# else
400# ifdef VERIFICATION
401 CALL netcdf_get_fvar (ng, inlm, dav(ng)%name, &
402 & vname(1,idnlmo), mod_value, &
403 & ncid = dav(ng)%ncid, &
404 & start = (/1/), &
405 & total = (/ndatum(ng)/))
406# elif !defined I4DVAR_ANA_SENSITIVITY
407 CALL netcdf_get_fvar (ng, inlm, dav(ng)%name, &
408 & vname(1,idnlmo), mod_value, &
409 & ncid = dav(ng)%ncid, &
410 & start = (/1,max(1,outer)/), &
411 & total = (/ndatum(ng),1/))
412# endif
413# endif
414 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
415!
416!-----------------------------------------------------------------------
417! If appropriate, convert depth in vertical fractional coordinates to
418! actual depths in meters (negative: downwards).
419!-----------------------------------------------------------------------
420!
421 obs_work(1:ndatum(ng))=inival
422!
423 CALL obs_k2z (ng, 0, lm(ng)+1, 0, mm(ng)+1, &
424 & lbi, ubi, lbj, ubj, 0, n(ng), &
425 & rxmin(ng), rxmax(ng), &
426 & rymin(ng), rymax(ng), &
427 & ndatum(ng), &
428 & obs_xgrid, obs_ygrid, obs_depths, obs_scale, &
429 & grid(ng) % z0_w, &
430 & obs_work)
431
432# ifdef DISTRIBUTE
433!
434! Collect all the computed depths.
435!
436 CALL mp_collect (ng, inlm, ndatum(ng), inival, obs_work)
437# endif
438!
439! Write out to depths to NetCDF file.
440!
441 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
442 & vname(1,idobsd), obs_work, &
443 & (/1/), (/ndatum(ng)/), &
444 & ncid = dav(ng)%ncid)
445 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
446
447# if defined RBL4DVAR_ANA_SENSITIVITY || \
448 defined rbl4dvar_fct_sensitivity || \
449 defined r4dvar_ana_sensitivity
450!
451!-----------------------------------------------------------------------
452! Write out several variables from the 4D-Var estimation for
453! completeness.
454!-----------------------------------------------------------------------
455!
456! Inquire about input observations variables.
457!
458 CALL netcdf_inq_var (ng, inlm, lcz(ng)%name)
459 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
460!
461! 4D-Var innovation (observation minus background).
462!
463 IF (find_string(var_name,n_var,vname(1,idinno),vindex)) THEN
464 CALL netcdf_get_fvar (ng, inlm, lcz(ng)%name, &
465 & vname(1,idinno), obs_work, &
466 & start = (/1/), &
467 & total = (/ndatum(ng)/))
468 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
469
470 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
471 & vname(1,idinno), obs_work, &
472 & (/1/), (/ndatum(ng)/), &
473 & ncid = dav(ng)%ncid)
474 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
475 END IF
476!
477! 4D-Var increment (analysis minus background).
478!
479 IF (find_string(var_name,n_var,vname(1,idincr),vindex)) THEN
480 CALL netcdf_get_fvar (ng, inlm, lcz(ng)%name, &
481 & vname(1,idincr), obs_work, &
482 & start = (/1/), &
483 & total = (/ndatum(ng)/))
484 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
485
486 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
487 & vname(1,idincr), obs_work, &
488 & (/1/), (/ndatum(ng)/), &
489 & ncid = dav(ng)%ncid)
490 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
491 END IF
492!
493! 4D-Var residual (observation minus analysis).
494!
495 IF (find_string(var_name,n_var,vname(1,idresi),vindex)) THEN
496 CALL netcdf_get_fvar (ng, inlm, lcz(ng)%name, &
497 & vname(1,idresi), obs_work, &
498 & start = (/1/), &
499 & total = (/ndatum(ng)/))
500 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
501
502 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
503 & vname(1,idresi), obs_work, &
504 & (/1/), (/ndatum(ng)/), &
505 & ncid = dav(ng)%ncid)
506 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
507 END IF
508!
509! 4D-Var initial model-observation misfit.
510!
511 IF (find_string(var_name,n_var,vname(1,idmomi),vindex)) THEN
512 CALL netcdf_get_fvar (ng, inlm, lcz(ng)%name, &
513 & vname(1,idmomi), obs_work, &
514 & start = (/1/), &
515 & total = (/ndatum(ng)/))
516 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
517
518 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
519 & vname(1,idmomi), obs_work, &
520 & (/1/), (/ndatum(ng)/), &
521 & ncid = dav(ng)%ncid)
522 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
523 END IF
524!
525! 4D-Var final model-observation misfit.
526!
527 IF (find_string(var_name,n_var,vname(1,idmomf),vindex)) THEN
528 CALL netcdf_get_fvar (ng, inlm, lcz(ng)%name, &
529 & vname(1,idmomf), obs_work, &
530 & start = (/1/), &
531 & total = (/ndatum(ng)/))
532 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
533
534 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
535 & vname(1,idmomf), obs_work, &
536 & (/1/), (/ndatum(ng)/), &
537 & ncid = dav(ng)%ncid)
538 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
539 END IF
540!
541! Unvetted initial NLM at observation locations.
542!
543 IF (find_string(var_name,n_var,vname(1,idnlmp),vindex)) THEN
544 CALL netcdf_get_fvar (ng, inlm, lcz(ng)%name, &
545 & vname(1,idnlmp), obs_work, &
546 & start = (/1/), &
547 & total = (/ndatum(ng)/))
548 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
549
550 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
551 & vname(1,idnlmp), obs_work, &
552 & (/1/), (/ndatum(ng)/), &
553 & ncid = dav(ng)%ncid)
554 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
555 END IF
556!
557! Final NLM at observation locations.
558!
559 IF (find_string(var_name,n_var,vname(1,idnlmf),vindex)) THEN
560 CALL netcdf_get_fvar (ng, inlm, lcz(ng)%name, &
561 & vname(1,idnlmf), obs_work, &
562 & start = (/1/), &
563 & total = (/ndatum(ng)/))
564 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
565
566 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
567 & vname(1,idnlmf), obs_work, &
568 & (/1/), (/ndatum(ng)/), &
569 & ncid = dav(ng)%ncid)
570 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
571 END IF
572!
573! Unvetted NLM at observation locations per outer-loop.
574!
575 IF (find_string(var_name,n_var,vname(1,idnlmu),vindex)) THEN
576 DO i=1,nouter
577 CALL netcdf_get_fvar (ng, inlm, lcz(ng)%name, &
578 & vname(1,idnlmu), obs_work, &
579 & start = (/1,i/), &
580 total = (/ndatum(ng),1/))
581 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
582
583 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
584 & vname(1,idnlmu), obs_work, &
585 & (/1,i/), (/ndatum(ng),1/), &
586 & ncid = dav(ng)%ncid)
587 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
588 END DO
589 END IF
590!
591! NLM at observation locations per outer-loop
592!
593 IF (find_string(var_name,n_var,vname(1,idnlmo),vindex)) THEN
594 DO i=1,nouter
595 CALL netcdf_get_fvar (ng, inlm, lcz(ng)%name, &
596 & vname(1,idnlmo), obs_work, &
597 & start = (/1,i/), &
598 total = (/ndatum(ng),1/))
599 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
600
601 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
602 & vname(1,idnlmo), obs_work, &
603 & (/1,i/), (/ndatum(ng),1/), &
604 & ncid = dav(ng)%ncid)
605 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
606 END DO
607 END IF
608# endif
609
610# if defined I4DVAR || defined RBL4DVAR || \
611 defined r4dvar || defined sp4dvar || \
612 defined split_i4dvar || defined split_rbl4dvar || \
613 defined split_r4dvar || defined split_sp4dvar || \
614 defined verification
615!
616!-----------------------------------------------------------------------
617! Compute model and observations comparison statistics.
618!-----------------------------------------------------------------------
619!
620! Initialize.
621!
622 IF (saved_exit_flag.eq.noerror) THEN
623 DO i=1,nobsvar(ng)
624 cc(i)=0.0_r8
625 mb(i)=0.0_r8
626 mse(i)=0.0_r8
627 sde(i)=0.0_r8
628 mod_min(i)=large
629 mod_max(i)=-large
630 mod_mean(i)=0.0_r8
631 obs_min(i)=large
632 obs_max(i)=-large
633 obs_mean(i)=0.0_r8
634 mod_std(i)=0.0_r8
635 obs_std(i)=0.0_r8
636 ncount(i)=0
637 END DO
638!
639! Compute model and observations mean per each state variable.
640!
641 DO iobs=1,ndatum(ng)
642 IF (obs_scale(iobs).eq.1.0_r8) THEN
643 i=obstype2state(obs_type(iobs))
644 ncount(i)=ncount(i)+1
645 mod_min(i)=min(mod_min(i),mod_value(iobs))
646 obs_min(i)=min(obs_min(i),obs_value(iobs))
647 mod_max(i)=max(mod_max(i),mod_value(iobs))
648 obs_max(i)=max(obs_max(i),obs_value(iobs))
649 mod_mean(i)=mod_mean(i)+mod_value(iobs)
650 obs_mean(i)=obs_mean(i)+obs_value(iobs)
651 END IF
652 END DO
653 DO i=1,nobsvar(ng)
654 IF (ncount(i).gt.0) THEN
655 mod_mean(i)=mod_mean(i)/real(ncount(i),r8)
656 obs_mean(i)=obs_mean(i)/real(ncount(i),r8)
657 END IF
658 END DO
659!
660! Compute standard deviation and cross-correlation between model and
661! observations (CC).
662!
663 DO iobs=1,ndatum(ng)
664 IF (obs_scale(iobs).eq.1.0_r8) THEN
665 i=obstype2state(obs_type(iobs))
666 cff1=mod_value(iobs)-mod_mean(i)
667 cff2=obs_value(iobs)-obs_mean(i)
668 mod_std(i)=mod_std(i)+cff1*cff1
669 obs_std(i)=obs_std(i)+cff2*cff2
670 cc(i)=cc(i)+cff1*cff2
671 END IF
672 END DO
673 DO i=1,nobsvar(ng)
674 IF (ncount(i).gt.1) THEN
675 mod_std(i)=sqrt(mod_std(i)/real(ncount(i)-1,r8))
676 obs_std(i)=sqrt(obs_std(i)/real(ncount(i)-1,r8))
677 cc(i)=(cc(i)/real(ncount(i),r8))/(mod_std(i)*obs_std(i))
678 END IF
679 END DO
680!
681! Compute model bias (MB), standard deviation error (SDE), and mean
682! squared error (MSE).
683!
684 DO i=1,nobsvar(ng)
685 IF (ncount(i).gt.0) THEN
686 mb(i)=mod_mean(i)-obs_mean(i)
687 sde(i)=mod_std(i)-obs_std(i)
688 mse(i)=mb(i)*mb(i)+ &
689 & sde(i)*sde(i)+ &
690 & 2.0_r8*mod_std(i)*obs_std(i)*(1.0_r8-cc(i))
691 END IF
692 END DO
693!
694! Report model and observations comparison statistics.
695!
696
697 IF (master) THEN
698 ic=0
699 DO i=1,nobsvar(ng)
700 svar_name(i)=' '
701 text(i)=' '
702 IF (ncount(i).gt.0) THEN
703 ic=ic+1
704 is(ic)=i
705 svar_name(ic)=trim(obsname(i))
706 text(ic)='-----------'
707 END IF
708 END DO
709 WRITE(f1,10) max(1,ic)
710 WRITE(f2,20) max(1,ic)
711 WRITE(f3,30) max(1,ic)
712 WRITE(f4,40) max(1,ic)
713 WRITE(stdout,50)
714 WRITE(stdout,f1) (svar_name(i),i=1,ic)
715 WRITE(stdout,f2) (text(i),i=1,ic)
716 WRITE(stdout,f3) 'Observation Min ',(obs_min(is(i)),i=1,ic)
717 WRITE(stdout,f3) 'Observation Max ',(obs_max(is(i)),i=1,ic)
718 WRITE(stdout,f3) 'Observation Mean ',(obs_mean(is(i)),i=1,ic)
719 WRITE(stdout,f3) 'Observation STD ',(obs_std(is(i)),i=1,ic)
720 WRITE(stdout,f3) 'Model Min ',(mod_min(is(i)),i=1,ic)
721 WRITE(stdout,f3) 'Model Max ',(mod_max(is(i)),i=1,ic)
722 WRITE(stdout,f3) 'Model Mean ',(mod_mean(is(i)),i=1,ic)
723 WRITE(stdout,f3) 'Model STD ',(mod_std(is(i)),i=1,ic)
724 WRITE(stdout,f3) 'Model Bias ',(mb(is(i)),i=1,ic)
725 WRITE(stdout,f3) 'STD Error ',(sde(is(i)),i=1,ic)
726 WRITE(stdout,f3) 'Cross-Correlation ',(cc(is(i)),i=1,ic)
727 WRITE(stdout,f3) 'Mean Squared Error',(mse(is(i)),i=1,ic)
728 WRITE(stdout,f4) 'Observation Count ',(ncount(is(i)),i=1,ic)
729 END IF
730!
731! Write comparison statistics to NetCDF file.
732!
733 CALL netcdf_put_ivar (ng, inlm, dav(ng)%name, &
734 & 'Nused_obs', ncount, &
735 & (/1/), (/nobsvar(ng)/), &
736 & ncid = dav(ng)%ncid)
737 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
738
739 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
740 & 'obs_mean', obs_mean, &
741 & (/1/), (/nobsvar(ng)/), &
742 & ncid = dav(ng)%ncid)
743 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
744
745 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
746 & 'obs_std', obs_std, &
747 & (/1/), (/nobsvar(ng)/), &
748 & ncid = dav(ng)%ncid)
749 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
750
751 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
752 & 'model_mean', mod_mean, &
753 & (/1/), (/nobsvar(ng)/), &
754 & ncid = dav(ng)%ncid)
755 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
756
757 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
758 & 'model_std', mod_std, &
759 & (/1/), (/nobsvar(ng)/), &
760 & ncid = dav(ng)%ncid)
761 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
762
763 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
764 & 'model_bias', mb, &
765 & (/1/), (/nobsvar(ng)/), &
766 & ncid = dav(ng)%ncid)
767 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
768
769 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
770 & 'SDE', sde, &
771 & (/1/), (/nobsvar(ng)/), &
772 & ncid = dav(ng)%ncid)
773 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
774
775 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
776 & 'CC', cc, &
777 & (/1/), (/nobsvar(ng)/), &
778 & ncid = dav(ng)%ncid)
779 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
780
781 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
782 & 'MSE', mse, &
783 & (/1/), (/nobsvar(ng)/), &
784 & ncid = dav(ng)%ncid)
785 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
786 END IF
787# endif
788!
789! Synchronize NetCDF file to disk.
790!
791 CALL netcdf_sync (ng, inlm, dav(ng)%name, dav(ng)%ncid)
792 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
793!
794! Reset "exit_flag" to its saved value when entering this routine.
795!
796 exit_flag=saved_exit_flag
797!
798 10 FORMAT ('(t22,',i2,'(a11,1x))')
799 20 FORMAT ('(t22,',i2,'(a11,1x),/)')
800 30 FORMAT ('(a,3x,',i2,'(1p,e11.4,0p,1x))')
801 40 FORMAT ('(a,3x,',i2,'(i11,1x))')
802 50 FORMAT (/,' Model-Observations Comparison Statistics:',/)
803
804 RETURN

References mod_scalars::blowup, mod_iounits::dav, mod_scalars::exit_flag, strings_mod::find_string(), strings_mod::founderror(), mod_grid::grid, mod_fourdvar::haveobsmeta, mod_ncparam::idincr, mod_ncparam::idinno, mod_ncparam::idmomf, mod_ncparam::idmomi, mod_ncparam::idnlmf, mod_ncparam::idnlmo, mod_ncparam::idnlmp, mod_ncparam::idnlmu, mod_ncparam::idnobs, mod_ncparam::idobsd, mod_ncparam::idobss, mod_ncparam::idobst, mod_ncparam::idobsx, mod_ncparam::idobsy, mod_ncparam::idoday, mod_ncparam::idoerr, mod_ncparam::idolat, mod_ncparam::idolon, mod_ncparam::idomet, mod_ncparam::idopro, mod_ncparam::idotyp, mod_ncparam::idoval, mod_ncparam::idresi, mod_ncparam::idtlmo, mod_param::inlm, mod_scalars::large, mod_iounits::lcz, mod_param::lm, mod_parallel::master, mod_param::mm, mod_param::n, mod_netcdf::n_var, mod_fourdvar::ndatum, mod_netcdf::netcdf_inq_var(), mod_netcdf::netcdf_sync(), mod_fourdvar::nobsvar, mod_scalars::noerror, mod_scalars::nouter, mod_fourdvar::nsurvey, mod_iounits::obs, obs_k2z_mod::obs_k2z(), mod_fourdvar::obsname, mod_fourdvar::obstype2state, mod_scalars::outer, mod_scalars::rclock, mod_ncparam::rxmax, mod_ncparam::rxmin, mod_ncparam::rymax, mod_ncparam::rymin, mod_iounits::sourcefile, mod_iounits::stdout, mod_netcdf::var_name, and mod_ncparam::vname.

Referenced by stats_modobs().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ stats_modobs_pio()

subroutine, private stats_modobs_mod::stats_modobs_pio ( integer, intent(in) ng,
integer, intent(in) tile,
integer, intent(in) lbi,
integer, intent(in) ubi,
integer, intent(in) lbj,
integer, intent(in) ubj )
private

Definition at line 810 of file stats_modobs.F.

812!***********************************************************************
813!
814! Imported variable declarations.
815!
816 integer, intent(in) :: ng, tile
817 integer, intent(in) :: LBi, UBi, LBj, UBj
818!
819! Local variable declarations.
820!
821 integer :: i, ic, iobs, status, vindex
822 integer :: saved_exit_flag
823
824 integer, dimension(2) :: start, total
825 integer, dimension(NobsVar(ng)) :: Ncount, is
826
827 integer, dimension(Ndatum(ng)) :: obs_type
828 integer, dimension(Ndatum(ng)) :: provenance
829 integer, dimension(Nsurvey(ng)) :: survey_Nobs
830!
831 real(r8) :: cff1, cff2
832 real(r8), parameter :: IniVal = 0.0_r8
833
834 real(r8), dimension(NobsVar(ng)) :: CC, MB, MSE, SDE
835 real(r8), dimension(NobsVar(ng)) :: mod_min, mod_max
836 real(r8), dimension(NobsVar(ng)) :: mod_mean, mod_std
837 real(r8), dimension(NobsVar(ng)) :: obs_min, obs_max
838 real(r8), dimension(NobsVar(ng)) :: obs_mean, obs_std
839
840 real(r8), dimension(Ndatum(ng)) :: mod_value
841 real(r8), dimension(Ndatum(ng)) :: obs_depths
842 real(r8), dimension(Ndatum(ng)) :: obs_scale
843 real(dp), dimension(Ndatum(ng)) :: obs_time
844 real(r8), dimension(Ndatum(ng)) :: obs_value
845 real(r8), dimension(Ndatum(ng)) :: obs_Xgrid
846 real(r8), dimension(Ndatum(ng)) :: obs_Ygrid
847 real(r8), dimension(Ndatum(ng)) :: obs_work
848 real(dp), dimension(Nsurvey(ng)) :: survey
849!
850 character (len=11), dimension(NobsVar(ng)) :: text, svar_name
851
852 character (len=30) :: F1, F2, F3, F4
853
854 character (len=*), parameter :: MyFile = &
855 & __FILE__//", stats_modobs_pio"
856!
857 sourcefile=myfile
858!
859!-----------------------------------------------------------------------
860! Read in model and observations data.
861!-----------------------------------------------------------------------
862!
863! If exit_flag > 0, save and reset to allow completion of the DAV
864! when ROMS blows-up.
865!
866 IF ((exit_flag.eq.1).or.(blowup.gt.0)) THEN
867 saved_exit_flag=exit_flag;
868 exit_flag=noerror
869 ELSE
870 saved_exit_flag=noerror;
871 END IF
872!
873! Inquire about input observations variables.
874!
875 CALL pio_netcdf_inq_var (ng, inlm, obs(ng)%name, &
876 & piofile = obs(ng)%pioFile)
877 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
878!
879! Read in and write out number of observations per survey.
880!
881 CALL pio_netcdf_get_ivar (ng, inlm, obs(ng)%name, &
882 & vname(1,idnobs), survey_nobs, &
883 & piofile = obs(ng)%pioFile, &
884 & start = (/1/), &
885 & total = (/nsurvey(ng)/))
886 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
887
888 CALL pio_netcdf_put_ivar (ng, inlm, dav(ng)%name, &
889 & vname(1,idnobs), survey_nobs, &
890 & (/1/), (/nsurvey(ng)/), &
891 & piofile = dav(ng)%pioFile)
892 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
893!
894! Read in and write out survey time.
895!
896 CALL pio_netcdf_get_time (ng, inlm, obs(ng)%name, &
897 & vname(1,idoday), &
898 & rclock%DateNumber, survey, &
899 & piofile = obs(ng)%pioFile, &
900 & start = (/1/), &
901 & total = (/nsurvey(ng)/))
902 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
903
904 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
905 & vname(1,idoday), survey, &
906 & (/1/), (/nsurvey(ng)/), &
907 & piofile = dav(ng)%pioFile)
908 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
909!
910! Read in and write out observation type identifier.
911!
912 CALL pio_netcdf_get_ivar (ng, inlm, obs(ng)%name, &
913 & vname(1,idotyp), obs_type, &
914 & piofile = obs(ng)%pioFile, &
915 & start = (/1/), &
916 & total = (/ndatum(ng)/))
917 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
918
919 CALL pio_netcdf_put_ivar (ng, inlm, dav(ng)%name, &
920 & vname(1,idotyp), obs_type, &
921 & (/1/), (/ndatum(ng)/), &
922 & piofile = dav(ng)%pioFile)
923 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
924!
925! Read in and write out observation provenance.
926!
927 CALL pio_netcdf_get_ivar (ng, inlm, obs(ng)%name, &
928 & vname(1,idopro), provenance, &
929 & piofile = obs(ng)%pioFile, &
930 & start = (/1/), &
931 & total = (/ndatum(ng)/))
932 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
933
934 CALL pio_netcdf_put_ivar (ng, inlm, dav(ng)%name, &
935 & vname(1,idopro), provenance, &
936 & (/1/), (/ndatum(ng)/), &
937 & piofile = dav(ng)%pioFile)
938 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
939!
940! Read in and write out observation time.
941!
942 CALL pio_netcdf_get_time (ng, inlm, obs(ng)%name, &
943 & vname(1,idobst), &
944 & rclock%DateNumber, obs_time, &
945 & piofile = obs(ng)%pioFile, &
946 & start = (/1/), &
947 & total = (/ndatum(ng)/))
948 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
949
950 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
951 & vname(1,idobst), obs_time, &
952 & (/1/), (/ndatum(ng)/), &
953 & piofile = dav(ng)%pioFile)
954 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
955!
956! Read in and write out observation longitude.
957!
958 IF (find_string(var_name,n_var,vname(1,idolon),vindex)) THEN
959 CALL pio_netcdf_get_fvar (ng, inlm, obs(ng)%name, &
960 & vname(1,idolon), obs_work, &
961 & piofile = obs(ng)%pioFile, &
962 & start = (/1/), &
963 & total = (/ndatum(ng)/))
964 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
965
966 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
967 & vname(1,idolon), obs_work, &
968 & (/1/), (/ndatum(ng)/), &
969 & piofile = dav(ng)%pioFile)
970 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
971 END IF
972!
973! Read in and write out observation latitude.
974!
975 IF (find_string(var_name,n_var,vname(1,idolon),vindex)) THEN
976 CALL pio_netcdf_get_fvar (ng, inlm, obs(ng)%name, &
977 & vname(1,idolat), obs_work, &
978 & piofile = obs(ng)%pioFile, &
979 & start = (/1/), &
980 & total = (/ndatum(ng)/))
981 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
982
983 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
984 & vname(1,idolat), obs_work, &
985 & (/1/), (/ndatum(ng)/), &
986 & piofile = dav(ng)%pioFile)
987 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
988 END IF
989!
990! Read in observation depth.
991!
992 CALL pio_netcdf_get_fvar (ng, inlm, obs(ng)%name, &
993 & vname(1,idobsd), obs_depths, &
994 & piofile = obs(ng)%pioFile, &
995 & start = (/1/), &
996 & total = (/ndatum(ng)/))
997 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
998!
999! Read in and write out X-fractional coordinate.
1000!
1001 CALL pio_netcdf_get_fvar (ng, inlm, obs(ng)%name, &
1002 & vname(1,idobsx), obs_xgrid, &
1003 & piofile = obs(ng)%pioFile, &
1004 & start = (/1/), &
1005 & total = (/ndatum(ng)/))
1006 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1007
1008 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
1009 & vname(1,idobsx), obs_xgrid, &
1010 & (/1/), (/ndatum(ng)/), &
1011 & piofile = dav(ng)%pioFile)
1012 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1013!
1014! Read in and write out Y-fractional coordinate.
1015!
1016 CALL pio_netcdf_get_fvar (ng, inlm, obs(ng)%name, &
1017 & vname(1,idobsy), obs_ygrid, &
1018 & piofile = obs(ng)%pioFile, &
1019 & start = (/1/), &
1020 & total = (/ndatum(ng)/))
1021 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1022
1023 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
1024 & vname(1,idobsy), obs_ygrid, &
1025 & (/1/), (/ndatum(ng)/), &
1026 & piofile = dav(ng)%pioFile)
1027 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1028!
1029! Read in and write out observations total error covariance
1030! (instrument + sampling + representation).
1031!
1032 CALL pio_netcdf_get_fvar (ng, inlm, obs(ng)%name, &
1033 & vname(1,idoerr), obs_work, &
1034 & piofile = obs(ng)%pioFile, &
1035 & start = (/1/), &
1036 & total = (/ndatum(ng)/))
1037 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1038
1039 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
1040 & vname(1,idoerr), obs_work, &
1041 & (/1/), (/ndatum(ng)/), &
1042 & piofile = dav(ng)%pioFile)
1043 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1044!
1045! Read in and write out observation values.
1046!
1047 CALL pio_netcdf_get_fvar (ng, inlm, obs(ng)%name, &
1048 & vname(1,idoval), obs_value, &
1049 & piofile = obs(ng)%pioFile, &
1050 & start = (/1/), &
1051 & total = (/ndatum(ng)/))
1052 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1053
1054 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
1055 & vname(1,idoval), obs_value, &
1056 & (/1/), (/ndatum(ng)/), &
1057 & piofile = dav(ng)%pioFile)
1058 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1059!
1060! Read in and write out observation meta values.
1061!
1062 IF (haveobsmeta(ng)) THEN
1063 CALL pio_netcdf_get_fvar (ng, inlm, obs(ng)%name, &
1064 & vname(1,idomet), obs_work, &
1065 & piofile = obs(ng)%pioFile, &
1066 & start = (/1/), &
1067 & total = (/ndatum(ng)/))
1068 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1069
1070 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
1071 & vname(1,idomet), obs_work, &
1072 & (/1/), (/ndatum(ng)/), &
1073 & piofile = dav(ng)%pioFile)
1074 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1075 END IF
1076!
1077! Read in observation screening flag.
1078!
1079 CALL pio_netcdf_get_fvar (ng, inlm, dav(ng)%name, &
1080 & vname(1,idobss), obs_scale, &
1081 & piofile = dav(ng)%pioFile, &
1082 & start = (/1/), &
1083 & total = (/ndatum(ng)/))
1084 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1085!
1086! Read in model values at observation locations.
1087!
1088# if defined R4DVAR || defined R4DVAR_ANA_SENSITIVITY || \
1089 defined tl_r4dvar
1090 CALL pio_netcdf_get_fvar (ng, inlm, dav(ng)%name, &
1091 & vname(1,idtlmo), mod_value, &
1092 & piofile = dav(ng)%pioFile, &
1093 & start = (/1/), &
1094 & total = (/ndatum(ng)/))
1095# else
1096# ifdef VERIFICATION
1097 CALL pio_netcdf_get_fvar (ng, inlm, dav(ng)%name, &
1098 & vname(1,idnlmo), mod_value, &
1099 & piofile = dav(ng)%pioFile, &
1100 & start = (/1/), &
1101 & total = (/ndatum(ng)/))
1102# elif !defined I4DVAR_ANA_SENSITIVITY
1103 CALL pio_netcdf_get_fvar (ng, inlm, dav(ng)%name, &
1104 & vname(1,idnlmo), mod_value, &
1105 & piofile = dav(ng)%pioFile, &
1106 & start = (/1,max(1,outer)/), &
1107 & total = (/ndatum(ng),1/))
1108# endif
1109# endif
1110 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1111!
1112!-----------------------------------------------------------------------
1113! If appropriate, convert depth in vertical fractional coordinates to
1114! actual depths in meters (negative: downwards).
1115!-----------------------------------------------------------------------
1116!
1117 obs_work(1:ndatum(ng))=inival
1118!
1119 CALL obs_k2z (ng, 0, lm(ng)+1, 0, mm(ng)+1, &
1120 & lbi, ubi, lbj, ubj, 0, n(ng), &
1121 & rxmin(ng), rxmax(ng), &
1122 & rymin(ng), rymax(ng), &
1123 & ndatum(ng), &
1124 & obs_xgrid, obs_ygrid, obs_depths, obs_scale, &
1125 & grid(ng) % z0_w, &
1126 & obs_work)
1127
1128# ifdef DISTRIBUTE
1129!
1130! Collect all the computed depths.
1131!
1132 CALL mp_collect (ng, inlm, ndatum(ng), inival, obs_work)
1133# endif
1134!
1135! Write out to depths to NetCDF file.
1136!
1137 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, vname(1,idobsd), &
1138 & obs_work, (/1/), (/ndatum(ng)/), &
1139 & piofile = dav(ng)%pioFile)
1140 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1141
1142# if defined RBL4DVAR_ANA_SENSITIVITY || \
1143 defined rbl4dvar_fct_sensitivity || \
1144 defined r4dvar_ana_sensitivity
1145!
1146!-----------------------------------------------------------------------
1147! Write out several variables from the 4D-Var estimation for
1148! completeness.
1149!-----------------------------------------------------------------------
1150!
1151! Inquire about input observations variables.
1152!
1153 CALL pio_netcdf_inq_var (ng, inlm, lcz(ng)%name)
1154 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1155!
1156! 4D-Var innovation (observation minus background).
1157!
1158 IF (find_string(var_name,n_var,vname(1,idinno),vindex)) THEN
1159 CALL pio_netcdf_get_fvar (ng, inlm, lcz(ng)%name, &
1160 & vname(1,idinno), obs_work, &
1161 & start = (/1/), &
1162 & total = (/ndatum(ng)/))
1163 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1164
1165 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
1166 & vname(1,idinno), obs_work, &
1167 & (/1/), (/ndatum(ng)/), &
1168 & piofile = dav(ng)%pioFile)
1169 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1170 END IF
1171!
1172! 4D-Var increment (analysis minus background).
1173!
1174 IF (find_string(var_name,n_var,vname(1,idincr),vindex)) THEN
1175 CALL pio_netcdf_get_fvar (ng, inlm, lcz(ng)%name, &
1176 & vname(1,idincr), obs_work, &
1177 & start = (/1/), &
1178 & total = (/ndatum(ng)/))
1179 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1180
1181 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
1182 & vname(1,idincr), obs_work, &
1183 & (/1/), (/ndatum(ng)/), &
1184 & piofile = dav(ng)%pioFile)
1185 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1186 END IF
1187!
1188! 4D-Var residual (observation minus analysis).
1189!
1190 IF (find_string(var_name,n_var,vname(1,idresi),vindex)) THEN
1191 CALL pio_netcdf_get_fvar (ng, inlm, lcz(ng)%name, &
1192 & vname(1,idresi), obs_work, &
1193 & start = (/1/), &
1194 & total = (/ndatum(ng)/))
1195 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1196
1197 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
1198 & vname(1,idresi), obs_work, &
1199 & (/1/), (/ndatum(ng)/), &
1200 & piofile = dav(ng)%pioFile)
1201 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1202 END IF
1203!
1204! 4D-Var initial model-observation misfit.
1205!
1206 IF (find_string(var_name,n_var,vname(1,idmomi),vindex)) THEN
1207 CALL pio_netcdf_get_fvar (ng, inlm, lcz(ng)%name, &
1208 & vname(1,idmomi), obs_work, &
1209 & start = (/1/), &
1210 & total = (/ndatum(ng)/))
1211 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1212
1213 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
1214 & vname(1,idmomi), obs_work, &
1215 & (/1/), (/ndatum(ng)/), &
1216 & piofile = dav(ng)%pioFile)
1217 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1218 END IF
1219!
1220! 4D-Var final model-observation misfit.
1221!
1222 IF (find_string(var_name,n_var,vname(1,idmomf),vindex)) THEN
1223 CALL pio_netcdf_get_fvar (ng, inlm, lcz(ng)%name, &
1224 & vname(1,idmomf), obs_work, &
1225 & start = (/1/), &
1226 & total = (/ndatum(ng)/))
1227 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1228
1229 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
1230 & vname(1,idmomf), obs_work, &
1231 & (/1/), (/ndatum(ng)/), &
1232 & piofile = dav(ng)%pioFile)
1233 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1234 END IF
1235!
1236! Unvetted initial NLM at observation locations.
1237!
1238 IF (find_string(var_name,n_var,vname(1,idnlmp),vindex)) THEN
1239 CALL pio_netcdf_get_fvar (ng, inlm, lcz(ng)%name, &
1240 & vname(1,idnlmp), obs_work, &
1241 & start = (/1/), &
1242 & total = (/ndatum(ng)/))
1243 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1244
1245 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
1246 & vname(1,idnlmp), obs_work, &
1247 & (/1/), (/ndatum(ng)/), &
1248 & piofile = dav(ng)%pioFile)
1249 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1250 END IF
1251!
1252! Final NLM at observation locations.
1253!
1254 IF (find_string(var_name,n_var,vname(1,idnlmf),vindex)) THEN
1255 CALL pio_netcdf_get_fvar (ng, inlm, lcz(ng)%name, &
1256 & vname(1,idnlmf), obs_work, &
1257 & start = (/1/), &
1258 & total = (/ndatum(ng)/))
1259 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1260
1261 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
1262 & vname(1,idnlmf), obs_work, &
1263 & (/1/), (/ndatum(ng)/), &
1264 & piofile = dav(ng)%pioFile)
1265 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1266 END IF
1267!
1268! Unvetted NLM at observation locations per outer-loop.
1269!
1270 IF (find_string(var_name,n_var,vname(1,idnlmu),vindex)) THEN
1271 DO i=1,nouter
1272 CALL pio_netcdf_get_fvar (ng, inlm, lcz(ng)%name, &
1273 & vname(1,idnlmu), obs_work, &
1274 & start = (/1,i/), &
1275 total = (/ndatum(ng),1/))
1276 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1277
1278 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
1279 & vname(1,idnlmu), obs_work, &
1280 & (/1,i/), (/ndatum(ng),1/), &
1281 & piofile = dav(ng)%pioFile)
1282 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1283 END DO
1284 END IF
1285!
1286! NLM at observation locations per outer-loop
1287!
1288 IF (find_string(var_name,n_var,vname(1,idnlmo),vindex)) THEN
1289 DO i=1,nouter
1290 CALL pio_netcdf_get_fvar (ng, inlm, lcz(ng)%name, &
1291 & vname(1,idnlmo), obs_work, &
1292 & start = (/1,i/), &
1293 total = (/ndatum(ng),1/))
1294 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1295
1296 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
1297 & vname(1,idnlmo), obs_work, &
1298 & (/1,i/), (/ndatum(ng),1/), &
1299 & piofile = dav(ng)%pioFile)
1300 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1301 END DO
1302 END IF
1303# endif
1304
1305# if defined I4DVAR || defined RBL4DVAR || \
1306 defined r4dvar || defined sp4dvar || \
1307 defined split_i4dvar || defined split_rbl4dvar || \
1308 defined split_r4dvar || defined split_sp4dvar || \
1309 defined verification
1310!
1311!-----------------------------------------------------------------------
1312! Compute model and observations comparison statistics.
1313!-----------------------------------------------------------------------
1314!
1315! Initialize.
1316!
1317 IF (saved_exit_flag.eq.noerror) THEN
1318 DO i=1,nobsvar(ng)
1319 cc(i)=0.0_r8
1320 mb(i)=0.0_r8
1321 mse(i)=0.0_r8
1322 sde(i)=0.0_r8
1323 mod_min(i)=large
1324 mod_max(i)=-large
1325 mod_mean(i)=0.0_r8
1326 obs_min(i)=large
1327 obs_max(i)=-large
1328 obs_mean(i)=0.0_r8
1329 mod_std(i)=0.0_r8
1330 obs_std(i)=0.0_r8
1331 ncount(i)=0
1332 END DO
1333!
1334! Compute model and observations mean per each state variable.
1335!
1336 DO iobs=1,ndatum(ng)
1337 IF (obs_scale(iobs).eq.1.0_r8) THEN
1338 i=obstype2state(obs_type(iobs))
1339 ncount(i)=ncount(i)+1
1340 mod_min(i)=min(mod_min(i),mod_value(iobs))
1341 obs_min(i)=min(obs_min(i),obs_value(iobs))
1342 mod_max(i)=max(mod_max(i),mod_value(iobs))
1343 obs_max(i)=max(obs_max(i),obs_value(iobs))
1344 mod_mean(i)=mod_mean(i)+mod_value(iobs)
1345 obs_mean(i)=obs_mean(i)+obs_value(iobs)
1346 END IF
1347 END DO
1348 DO i=1,nobsvar(ng)
1349 IF (ncount(i).gt.0) THEN
1350 mod_mean(i)=mod_mean(i)/real(ncount(i),r8)
1351 obs_mean(i)=obs_mean(i)/real(ncount(i),r8)
1352 END IF
1353 END DO
1354!
1355! Compute standard deviation and cross-correlation between model and
1356! observations (CC).
1357!
1358 DO iobs=1,ndatum(ng)
1359 IF (obs_scale(iobs).eq.1.0_r8) THEN
1360 i=obstype2state(obs_type(iobs))
1361 cff1=mod_value(iobs)-mod_mean(i)
1362 cff2=obs_value(iobs)-obs_mean(i)
1363 mod_std(i)=mod_std(i)+cff1*cff1
1364 obs_std(i)=obs_std(i)+cff2*cff2
1365 cc(i)=cc(i)+cff1*cff2
1366 END IF
1367 END DO
1368 DO i=1,nobsvar(ng)
1369 IF (ncount(i).gt.1) THEN
1370 mod_std(i)=sqrt(mod_std(i)/real(ncount(i)-1,r8))
1371 obs_std(i)=sqrt(obs_std(i)/real(ncount(i)-1,r8))
1372 cc(i)=(cc(i)/real(ncount(i),r8))/(mod_std(i)*obs_std(i))
1373 END IF
1374 END DO
1375!
1376! Compute model bias (MB), standard deviation error (SDE), and mean
1377! squared error (MSE).
1378!
1379 DO i=1,nobsvar(ng)
1380 IF (ncount(i).gt.0) THEN
1381 mb(i)=mod_mean(i)-obs_mean(i)
1382 sde(i)=mod_std(i)-obs_std(i)
1383 mse(i)=mb(i)*mb(i)+ &
1384 & sde(i)*sde(i)+ &
1385 & 2.0_r8*mod_std(i)*obs_std(i)*(1.0_r8-cc(i))
1386 END IF
1387 END DO
1388!
1389! Report model and observations comparison statistics.
1390!
1391
1392 IF (master) THEN
1393 ic=0
1394 DO i=1,nobsvar(ng)
1395 svar_name(i)=' '
1396 text(i)=' '
1397 IF (ncount(i).gt.0) THEN
1398 ic=ic+1
1399 is(ic)=i
1400 svar_name(ic)=trim(obsname(i))
1401 text(ic)='-----------'
1402 END IF
1403 END DO
1404 WRITE(f1,10) max(1,ic)
1405 WRITE(f2,20) max(1,ic)
1406 WRITE(f3,30) max(1,ic)
1407 WRITE(f4,40) max(1,ic)
1408 WRITE(stdout,50)
1409 WRITE(stdout,f1) (svar_name(i),i=1,ic)
1410 WRITE(stdout,f2) (text(i),i=1,ic)
1411 WRITE(stdout,f3) 'Observation Min ',(obs_min(is(i)),i=1,ic)
1412 WRITE(stdout,f3) 'Observation Max ',(obs_max(is(i)),i=1,ic)
1413 WRITE(stdout,f3) 'Observation Mean ',(obs_mean(is(i)),i=1,ic)
1414 WRITE(stdout,f3) 'Observation STD ',(obs_std(is(i)),i=1,ic)
1415 WRITE(stdout,f3) 'Model Min ',(mod_min(is(i)),i=1,ic)
1416 WRITE(stdout,f3) 'Model Max ',(mod_max(is(i)),i=1,ic)
1417 WRITE(stdout,f3) 'Model Mean ',(mod_mean(is(i)),i=1,ic)
1418 WRITE(stdout,f3) 'Model STD ',(mod_std(is(i)),i=1,ic)
1419 WRITE(stdout,f3) 'Model Bias ',(mb(is(i)),i=1,ic)
1420 WRITE(stdout,f3) 'STD Error ',(sde(is(i)),i=1,ic)
1421 WRITE(stdout,f3) 'Cross-Correlation ',(cc(is(i)),i=1,ic)
1422 WRITE(stdout,f3) 'Mean Squared Error',(mse(is(i)),i=1,ic)
1423 WRITE(stdout,f4) 'Observation Count ',(ncount(is(i)),i=1,ic)
1424 END IF
1425!
1426! Write comparison statistics to NetCDF file.
1427!
1428 CALL pio_netcdf_put_ivar (ng, inlm, dav(ng)%name, &
1429 & 'Nused_obs', ncount, &
1430 & (/1/), (/nobsvar(ng)/), &
1431 & piofile = dav(ng)%pioFile)
1432 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1433
1434 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
1435 & 'obs_mean', obs_mean, &
1436 & (/1/), (/nobsvar(ng)/), &
1437 & piofile = dav(ng)%pioFile)
1438 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1439
1440 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
1441 & 'obs_std', obs_std, &
1442 & (/1/), (/nobsvar(ng)/), &
1443 & piofile = dav(ng)%pioFile)
1444 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1445
1446 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
1447 & 'model_mean', mod_mean, &
1448 & (/1/), (/nobsvar(ng)/), &
1449 & piofile = dav(ng)%pioFile)
1450 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1451
1452 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
1453 & 'model_std', mod_std, &
1454 & (/1/), (/nobsvar(ng)/), &
1455 & piofile = dav(ng)%pioFile)
1456 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1457
1458 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
1459 & 'model_bias', mb, &
1460 & (/1/), (/nobsvar(ng)/), &
1461 & piofile = dav(ng)%pioFile)
1462 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1463
1464 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
1465 & 'SDE', sde, &
1466 & (/1/), (/nobsvar(ng)/), &
1467 & piofile = dav(ng)%pioFile)
1468 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1469
1470 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
1471 & 'CC', cc, &
1472 & (/1/), (/nobsvar(ng)/), &
1473 & piofile = dav(ng)%pioFile)
1474 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1475
1476 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
1477 & 'MSE', mse, &
1478 & (/1/), (/nobsvar(ng)/), &
1479 & piofile = dav(ng)%pioFile)
1480 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1481 END IF
1482# endif
1483!
1484! Synchronize NetCDF file to disk.
1485!
1486 CALL pio_netcdf_sync (ng, inlm, dav(ng)%name, dav(ng)%pioFile)
1487 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1488!
1489! Reset "exit_flag" to its saved value when entering this routine.
1490!
1491 exit_flag=saved_exit_flag
1492!
1493 10 FORMAT ('(t22,',i2,'(a11,1x))')
1494 20 FORMAT ('(t22,',i2,'(a11,1x),/)')
1495 30 FORMAT ('(a,3x,',i2,'(1p,e11.4,0p,1x))')
1496 40 FORMAT ('(a,3x,',i2,'(i11,1x))')
1497 50 FORMAT (/,' Model-Observations Comparison Statistics:',/)
1498
1499 RETURN

References mod_scalars::blowup, mod_iounits::dav, mod_scalars::exit_flag, strings_mod::find_string(), strings_mod::founderror(), mod_grid::grid, mod_fourdvar::haveobsmeta, mod_ncparam::idincr, mod_ncparam::idinno, mod_ncparam::idmomf, mod_ncparam::idmomi, mod_ncparam::idnlmf, mod_ncparam::idnlmo, mod_ncparam::idnlmp, mod_ncparam::idnlmu, mod_ncparam::idnobs, mod_ncparam::idobsd, mod_ncparam::idobss, mod_ncparam::idobst, mod_ncparam::idobsx, mod_ncparam::idobsy, mod_ncparam::idoday, mod_ncparam::idoerr, mod_ncparam::idolat, mod_ncparam::idolon, mod_ncparam::idomet, mod_ncparam::idopro, mod_ncparam::idotyp, mod_ncparam::idoval, mod_ncparam::idresi, mod_ncparam::idtlmo, mod_param::inlm, mod_scalars::large, mod_iounits::lcz, mod_param::lm, mod_parallel::master, mod_param::mm, mod_param::n, mod_netcdf::n_var, mod_fourdvar::ndatum, mod_fourdvar::nobsvar, mod_scalars::noerror, mod_scalars::nouter, mod_fourdvar::nsurvey, mod_iounits::obs, obs_k2z_mod::obs_k2z(), mod_fourdvar::obsname, mod_fourdvar::obstype2state, mod_scalars::outer, mod_pio_netcdf::pio_netcdf_inq_var(), mod_pio_netcdf::pio_netcdf_sync(), mod_scalars::rclock, mod_ncparam::rxmax, mod_ncparam::rxmin, mod_ncparam::rymax, mod_ncparam::rymin, mod_iounits::sourcefile, mod_iounits::stdout, mod_netcdf::var_name, and mod_ncparam::vname.

Referenced by stats_modobs().

Here is the call graph for this function:
Here is the caller graph for this function: