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

Functions/Subroutines

subroutine, public obs_write (ng, tile, model)
 
subroutine, private obs_operator (ng, tile, model, lbi, ubi, lbj, ubj, mstr, mend)
 
subroutine, private obs_write_nf90 (ng, tile, model, lbi, ubi, lbj, ubj, mstr, mend)
 
subroutine, private obs_write_pio (ng, tile, model, lbi, ubi, lbj, ubj, mstr, mend)
 

Function/Subroutine Documentation

◆ obs_operator()

subroutine, private obs_write_mod::obs_operator ( integer, intent(in) ng,
integer, intent(in) tile,
integer, intent(in) model,
integer, intent(in) lbi,
integer, intent(in) ubi,
integer, intent(in) lbj,
integer, intent(in) ubj,
integer, intent(out) mstr,
integer, intent(out) mend )
private

Definition at line 221 of file obs_write.F.

224!***********************************************************************
225!
226! Imported variable declarations.
227!
228 integer, intent(in) :: ng, tile, model
229 integer, intent(in) :: LBi, UBi, LBj, UBj
230
231 integer, intent(out) :: Mstr, Mend
232!
233! Local variable declarations.
234!
235 integer :: i, ic, iobs, itrc
236# ifdef SOLVE3D
237 integer :: j, k
238# endif
239# ifdef DISTRIBUTE
240 integer :: Ncollect
241# endif
242!
243 real(r8), parameter :: IniVal = 0.0_r8
244
245 real(r8) :: angle
246# ifdef BGQC
247 real(r8) :: df1, df2, Thresh
248# endif
249 real(r8) :: misfit(Mobs), unvetted(Mobs)
250 real(r8) :: uradial(Mobs), vradial(Mobs)
251# ifndef VERIFICATION
252 real(r8) :: uBgErr(Mobs), vBgErr(Mobs)
253# endif
254!
255 character (len=*), parameter :: MyFile = &
256 & __FILE__//", obs_operator"
257
258# include "set_bounds.h"
259!
260 sourcefile=myfile
261!
262!=======================================================================
263! Interpolate model state at observation locations.
264!=======================================================================
265!
266# ifdef WEAK_CONSTRAINT
267!
268! Set starting and ending indices of observations to process for the
269! current time survey. In the dual formulation, the entire observation
270! vector for the assimilation window is maintained and the data array
271! indexes are global.
272!
273 mstr=nstrobs(ng)
274 mend=nendobs(ng)
275# else
276!
277! Set starting and ending indices of observations to process for the
278! current time survey. In the primal formulation, only the observations
279! for the current assimilation window are maintained and the data array
280! indexes are local (1:Nobs; starting index is always one).
281!
282 mstr=1
283 mend=nobs(ng)
284# endif
285!
286! Initialize observation reject/accept processing flag.
287!
288 DO iobs=mstr,mend
289 unvetted(iobs)=inival
290 obsvetting(iobs)=inival
291 END DO
292
293# ifndef I4DVAR_ANA_SENSITIVITY
294!
295! Some entries are not computed in the extraction routine. Set values
296! to zero to avoid problems when writing non initialized values.
297# ifdef DISTRIBUTE
298! Notice that only the appropriate indices are zero-out to facilate
299! collecting all the extrated data as sum between all nodes.
300# endif
301!
302 IF (wrtnlmod(ng)) THEN
303 DO iobs=mstr,mend
304 nlmodval(iobs)=inival
305# ifndef VERIFICATION
306 bgerr(iobs)=inival
307 ubgerr(iobs)=inival
308 vbgerr(iobs)=inival
309# endif
310 END DO
311
312# ifdef BGQC
313!
314! Set background quality control check (BgThresh). The background
315! quality control is either in terms of the state variable index
316! (1:MstateVar) or observation provenance index.
317!
318 IF (bgqc_type(ng).eq.1) THEN ! state variable index QC
319 DO iobs=mstr,mend
320 DO i=1,mstatevar
321 IF (obstype(iobs).eq.obsstate2type(i)) THEN
322 bgthresh(iobs)=s_bgqc(i,ng)
323 EXIT
324 END IF
325 END DO
326 END DO
327 ELSE IF (bgqc_type(ng).eq.2) THEN ! provenance index QC
328 DO iobs=mstr,mend
329 bgthresh(iobs)=bgqc_large ! do not reject
330 DO i=1,nprovenance(ng)
331 IF (obsprov(iobs).eq.iprovenance(i,ng)) THEN
332 bgthresh(iobs)=p_bgqc(i,ng)
333 EXIT
334 END IF
335 END DO
336 END DO
337 END IF
338# endif
339 END IF
340
341# if defined TLM_OBS && !defined SP4DVAR
342 IF (wrttlmod(ng).or.wrtrpmod(ng)) THEN
343 DO iobs=mstr,mend
344 tlmodval(iobs)=inival
345 END DO
346 END IF
347# endif
348# endif
349!
350!-----------------------------------------------------------------------
351! Free-surface operator.
352!-----------------------------------------------------------------------
353!
354# ifndef I4DVAR_ANA_SENSITIVITY
355! Nonlinear free-surface at observation locations.
356!
357 IF (wrtnlmod(ng).and. &
358 & (fourdvar(ng)%ObsCount(isfsur).gt.0)) THEN
359 CALL extract_obs2d (ng, 0, lm(ng)+1, 0, mm(ng)+1, &
360 & lbi, ubi, lbj, ubj, &
361 & obsstate2type(isfsur), &
362 & mobs, mstr, mend, &
363 & rxmin(ng), rxmax(ng), &
364 & rymin(ng), rymax(ng), &
365 & time(ng), dt(ng), &
366 & obstype, obsvetting, &
367 & tobs, xobs, yobs, &
368 & ocean(ng)%zeta(:,:,kout), &
369# ifdef MASKING
370 & grid(ng)%rmask, &
371# endif
372 & nlmodval)
373# ifndef VERIFICATION
374!
375! Free-surface backgound error standard deviation at observation
376! locations.
377!
378 CALL extract_obs2d (ng, 0, lm(ng)+1, 0, mm(ng)+1, &
379 & lbi, ubi, lbj, ubj, &
380 & obsstate2type(isfsur), &
381 & mobs, mstr, mend, &
382 & rxmin(ng), rxmax(ng), &
383 & rymin(ng), rymax(ng), &
384 & time(ng), dt(ng), &
385 & obstype, obsvetting, &
386 & tobs, xobs, yobs, &
387 & ocean(ng)%e_zeta(:,:,1), &
388# ifdef MASKING
389 & grid(ng)%rmask, &
390# endif
391 & bgerr)
392# endif
393 END IF
394# endif
395
396# ifdef TLM_OBS
397!
398! Tangent linear free-surface at observation locations.
399!
400 IF ((wrttlmod(ng).or.(wrtrpmod(ng))).and. &
401 & (fourdvar(ng)%ObsCount(isfsur).gt.0)) THEN
402 CALL extract_obs2d (ng, 0, lm(ng)+1, 0, mm(ng)+1, &
403 & lbi, ubi, lbj, ubj, &
404 & obsstate2type(isfsur), &
405 & mobs, mstr, mend, &
406 & rxmin(ng), rxmax(ng), &
407 & rymin(ng), rymax(ng), &
408 & time(ng), dt(ng), &
409 & obstype, obsvetting, &
410 & tobs, xobs, yobs, &
411 & ocean(ng)%tl_zeta(:,:,kout), &
412# ifdef MASKING
413 & grid(ng)%rmask, &
414# endif
415 & tlmodval)
416 END IF
417# endif
418!
419!-----------------------------------------------------------------------
420! Vertically integrated u-velocity operator.
421!-----------------------------------------------------------------------
422!
423# ifndef I4DVAR_ANA_SENSITIVITY
424! Nonlinear 2D u-velocity at observation locations.
425!
426 IF (wrtnlmod(ng).and. &
427 & (fourdvar(ng)%ObsCount(isubar).gt.0)) THEN
428 CALL extract_obs2d (ng, 1, lm(ng)+1, 0, mm(ng)+1, &
429 & lbi, ubi, lbj, ubj, &
430 & obsstate2type(isubar), &
431 & mobs, mstr, mend, &
432 & uxmin(ng), uxmax(ng), &
433 & uymin(ng), uymax(ng), &
434 & time(ng), dt(ng), &
435 & obstype, obsvetting, &
436 & tobs, xobs, yobs, &
437 & ocean(ng)%ubar(:,:,kout), &
438# ifdef MASKING
439 & grid(ng)%umask, &
440# endif
441 & nlmodval)
442
443# ifndef VERIFICATION
444!
445! 2D u-velocity backgound error standard deviation at observation
446! locations.
447!
448 CALL extract_obs2d (ng, 1, lm(ng)+1, 0, mm(ng)+1, &
449 & lbi, ubi, lbj, ubj, &
450 & obsstate2type(isubar), &
451 & mobs, mstr, mend, &
452 & uxmin(ng), uxmax(ng), &
453 & uymin(ng), uymax(ng), &
454 & time(ng), dt(ng), &
455 & obstype, obsvetting, &
456 & tobs, xobs, yobs, &
457 & ocean(ng)%e_ubar(:,:,1), &
458# ifdef MASKING
459 & grid(ng)%umask, &
460# endif
461 & bgerr)
462# endif
463 END IF
464# endif
465
466# ifdef TLM_OBS
467!
468! Tangent linear 2D u-velocity at observation locations.
469!
470 IF ((wrttlmod(ng).or.(wrtrpmod(ng))).and. &
471 & (fourdvar(ng)%ObsCount(isubar).gt.0)) THEN
472 CALL extract_obs2d (ng, 1, lm(ng)+1, 0, mm(ng)+1, &
473 & lbi, ubi, lbj, ubj, &
474 & obsstate2type(isubar), &
475 & mobs, mstr, mend, &
476 & uxmin(ng), uxmax(ng), &
477 & uymin(ng), uymax(ng), &
478 & time(ng), dt(ng), &
479 & obstype, obsvetting, &
480 & tobs, xobs, yobs, &
481 & ocean(ng)%tl_ubar(:,:,kout), &
482# ifdef MASKING
483 & grid(ng)%umask, &
484# endif
485 & tlmodval)
486 END IF
487# endif
488!
489!-----------------------------------------------------------------------
490! Vertically integrated v-velocity observations.
491!-----------------------------------------------------------------------
492!
493# ifndef I4DVAR_ANA_SENSITIVITY
494! Nonlinear 2D v-velocity at observation locations.
495!
496 IF (wrtnlmod(ng).and. &
497 & (fourdvar(ng)%ObsCount(isvbar).gt.0)) THEN
498 CALL extract_obs2d (ng, 0, lm(ng)+1, 1, mm(ng)+1, &
499 & lbi, ubi, lbj, ubj, &
500 & obsstate2type(isvbar), &
501 & mobs, mstr, mend, &
502 & vxmin(ng), vxmax(ng), &
503 & vymin(ng), vymax(ng), &
504 & time(ng), dt(ng), &
505 & obstype, obsvetting, &
506 & tobs, xobs, yobs, &
507 & ocean(ng)%vbar(:,:,kout), &
508# ifdef MASKING
509 & grid(ng)%vmask, &
510# endif
511 & nlmodval)
512
513# ifndef VERIFICATION
514!
515! 2D v-velocity backgound error standard deviation at observation
516! locations.
517!
518 CALL extract_obs2d (ng, 0, lm(ng)+1, 1, mm(ng)+1, &
519 & lbi, ubi, lbj, ubj, &
520 & obsstate2type(isvbar), &
521 & mobs, mstr, mend, &
522 & vxmin(ng), vxmax(ng), &
523 & vymin(ng), vymax(ng), &
524 & time(ng), dt(ng), &
525 & obstype, obsvetting, &
526 & tobs, xobs, yobs, &
527 & ocean(ng)%e_vbar(:,:,1), &
528# ifdef MASKING
529 & grid(ng)%vmask, &
530# endif
531 & bgerr)
532# endif
533 END IF
534# endif
535
536# ifdef TLM_OBS
537!
538! Tangent linear 2D v-velocity at observation locations.
539!
540 IF ((wrttlmod(ng).or.(wrtrpmod(ng))).and. &
541 & (fourdvar(ng)%ObsCount(isvbar).gt.0)) THEN
542 CALL extract_obs2d (ng, 0, lm(ng)+1, 1, mm(ng)+1, &
543 & lbi, ubi, lbj, ubj, &
544 & obsstate2type(isvbar), &
545 & mobs, mstr, mend, &
546 & vxmin(ng), vxmax(ng), &
547 & vymin(ng), vymax(ng), &
548 & time(ng), dt(ng), &
549 & obstype, obsvetting, &
550 & tobs, xobs, yobs, &
551 & ocean(ng)%tl_vbar(:,:,kout), &
552# ifdef MASKING
553 & grid(ng)%vmask, &
554# endif
555 & tlmodval)
556 END IF
557# endif
558
559# ifdef SOLVE3D
560!
561!-----------------------------------------------------------------------
562! 3D u-velocity observations.
563!-----------------------------------------------------------------------
564!
565# ifndef I4DVAR_ANA_SENSITIVITY
566! Nonlinear 3D u-velocity at observation locations.
567!
568 IF (wrtnlmod(ng).and. &
569 & (fourdvar(ng)%ObsCount(isuvel).gt.0)) THEN
570 DO k=1,n(ng)
571 DO j=jstr-1,jend+1
572 DO i=istru-1,iend+1
573 grid(ng)%z_v(i,j,k)=0.5_r8*(grid(ng)%z_r(i-1,j,k)+ &
574 & grid(ng)%z_r(i ,j,k))
575 END DO
576 END DO
577 END DO
578 CALL extract_obs3d (ng, 1, lm(ng)+1, 0, mm(ng)+1, &
579 & lbi, ubi, lbj, ubj, 1, n(ng), &
580 & obsstate2type(isuvel), &
581 & mobs, mstr, mend, &
582 & uxmin(ng), uxmax(ng), &
583 & uymin(ng), uymax(ng), &
584 & time(ng), dt(ng), &
585 & obstype, obsvetting, &
586 & tobs, xobs, yobs, zobs, &
587 & ocean(ng)%u(:,:,:,nout), &
588 & grid(ng)%z_v, &
589# ifdef MASKING
590 & grid(ng)%umask, &
591# endif
592 & nlmodval)
593
594# ifndef VERIFICATION
595!
596! 3D u-velocity backgound error standard deviation at observation
597! locations.
598!
599 CALL extract_obs3d (ng, 1, lm(ng)+1, 0, mm(ng)+1, &
600 & lbi, ubi, lbj, ubj, 1, n(ng), &
601 & obsstate2type(isuvel), &
602 & mobs, mstr, mend, &
603 & uxmin(ng), uxmax(ng), &
604 & uymin(ng), uymax(ng), &
605 & time(ng), dt(ng), &
606 & obstype, obsvetting, &
607 & tobs, xobs, yobs, zobs, &
608 & ocean(ng)%e_u(:,:,:,1), &
609 & grid(ng)%z_v, &
610# ifdef MASKING
611 & grid(ng)%umask, &
612# endif
613 & bgerr)
614# endif
615 END IF
616# endif
617
618# ifdef TLM_OBS
619!
620! Tangent linear 3D u-velocity at observation locations.
621!
622 IF ((wrttlmod(ng).or.(wrtrpmod(ng))).and. &
623 & (fourdvar(ng)%ObsCount(isuvel).gt.0)) THEN
624 DO k=1,n(ng)
625 DO j=jstr-1,jend+1
626 DO i=istru-1,iend+1
627 grid(ng)%z_v(i,j,k)=0.5_r8*(grid(ng)%z_r(i-1,j,k)+ &
628 & grid(ng)%z_r(i ,j,k))
629 END DO
630 END DO
631 END DO
632 CALL extract_obs3d (ng, 1, lm(ng)+1, 0, mm(ng)+1, &
633 & lbi, ubi, lbj, ubj, 1, n(ng), &
634 & obsstate2type(isuvel), &
635 & mobs, mstr, mend, &
636 & uxmin(ng), uxmax(ng), &
637 & uymin(ng), uymax(ng), &
638 & time(ng), dt(ng), &
639 & obstype, obsvetting, &
640 & tobs, xobs, yobs, zobs, &
641 & ocean(ng)%tl_u(:,:,:,nout), &
642 & grid(ng)%z_v, &
643# ifdef MASKING
644 & grid(ng)%umask, &
645# endif
646 & tlmodval)
647 END IF
648# endif
649!
650!-----------------------------------------------------------------------
651! 3D v-velocity observations.
652!-----------------------------------------------------------------------
653!
654# ifndef I4DVAR_ANA_SENSITIVITY
655! Nonlinear 3D v-velocity at observation locations.
656!
657 IF (wrtnlmod(ng).and. &
658 & (fourdvar(ng)%ObsCount(isvvel).gt.0)) THEN
659 DO k=1,n(ng)
660 DO j=jstrv-1,jend+1
661 DO i=istr-1,iend+1
662 grid(ng)%z_v(i,j,k)=0.5_r8*(grid(ng)%z_r(i,j-1,k)+ &
663 & grid(ng)%z_r(i,j ,k))
664 END DO
665 END DO
666 END DO
667 CALL extract_obs3d (ng, 0, lm(ng)+1, 1, mm(ng)+1, &
668 & lbi, ubi, lbj, ubj, 1, n(ng), &
669 & obsstate2type(isvvel), &
670 & mobs, mstr, mend, &
671 & vxmin(ng), vxmax(ng), &
672 & vymin(ng), vymax(ng), &
673 & time(ng), dt(ng), &
674 & obstype, obsvetting, &
675 & tobs, xobs, yobs, zobs, &
676 & ocean(ng)%v(:,:,:,nout), &
677 & grid(ng)%z_v, &
678# ifdef MASKING
679 & grid(ng)%vmask, &
680# endif
681 & nlmodval)
682
683# ifndef VERIFICATION
684!
685! 3D v-velocity backgound error standard deviation at observation
686! locations.
687!
688 CALL extract_obs3d (ng, 0, lm(ng)+1, 1, mm(ng)+1, &
689 & lbi, ubi, lbj, ubj, 1, n(ng), &
690 & obsstate2type(isvvel), &
691 & mobs, mstr, mend, &
692 & vxmin(ng), vxmax(ng), &
693 & vymin(ng), vymax(ng), &
694 & time(ng), dt(ng), &
695 & obstype, obsvetting, &
696 & tobs, xobs, yobs, zobs, &
697 & ocean(ng)%e_v(:,:,:,1), &
698 & grid(ng)%z_v, &
699# ifdef MASKING
700 & grid(ng)%vmask, &
701# endif
702 & bgerr)
703# endif
704 END IF
705# endif
706
707# ifdef TLM_OBS
708!
709! Tangent linear 3D v-velocity at observation locations.
710!
711 IF ((wrttlmod(ng).or.(wrtrpmod(ng))).and. &
712 & (fourdvar(ng)%ObsCount(isvvel).gt.0)) THEN
713 DO k=1,n(ng)
714 DO j=jstrv-1,jend+1
715 DO i=istr-1,iend+1
716 grid(ng)%z_v(i,j,k)=0.5_r8*(grid(ng)%z_r(i,j-1,k)+ &
717 & grid(ng)%z_r(i,j ,k))
718 END DO
719 END DO
720 END DO
721 CALL extract_obs3d (ng, 0, lm(ng)+1, 1, mm(ng)+1, &
722 & lbi, ubi, lbj, ubj, 1, n(ng), &
723 & obsstate2type(isvvel), &
724 & mobs, mstr, mend, &
725 & vxmin(ng), vxmax(ng), &
726 & vymin(ng), vymax(ng), &
727 & time(ng), dt(ng), &
728 & obstype, obsvetting, &
729 & tobs, xobs, yobs, zobs, &
730 & ocean(ng)%tl_v(:,:,:,nout), &
731 & grid(ng)%z_v, &
732# ifdef MASKING
733 & grid(ng)%vmask, &
734# endif
735 & tlmodval)
736 END IF
737# endif
738!
739!-----------------------------------------------------------------------
740! Radial Velocity. The observations are in terms of radial speed and
741! angle (stored in obs_meta). The observation angle converts the
742! velocity components to geographical EAST and North components.
743# ifdef RADIAL_ANGLE_CCW_EAST
744! The radial velocity observations are processed as magnitude and
745! heading angle (obs_meta; radians) in the math convention: an
746! azimuth that is counterclockwise from TRUE East.
747!
748! In curvilinear coordinates, the radial forward problem is:
749!
750! radial = u * COS(obs_meta - angler) + v * SIN(obs_meta - angler)
751# else
752! By default, the radial velocity observations are processed as
753! magnitude and heading angle (obs_meta; radians) in the navigation
754! convention: an azimuth that is clockwise from TRUE North.
755!
756! In curvilinear coordinates, the radial forward problem is:
757!
758! radial = u * SIN(obs_meta + angler) + v * COS(obs_meta + angler)
759# endif
760!-----------------------------------------------------------------------
761!
762# ifndef I4DVAR_ANA_SENSITIVITY
763 IF (wrtnlmod(ng).and. &
764 & (fourdvar(ng)%ObsCount(isradial).gt.0)) THEN
765 DO iobs=mstr,mend
766 uradial(iobs)=inival
767 vradial(iobs)=inival
768 END DO
769
770# ifdef CURVGRID
771!
772! If curvilinear coordinates, interpolate grid rotation angle (radians)
773! for each observation location.
774!
775 CALL extract_obs2d (ng, 0, lm(ng)+1, 0, mm(ng)+1, &
776 & lbi, ubi, lbj, ubj, &
777 & obsstate2type(isradial), &
778 & mobs, mstr, mend, &
779 & rxmin(ng), rxmax(ng), &
780 & rymin(ng), rymax(ng), &
781 & time(ng), dt(ng), &
782 & obstype, obsvetting, &
783 & tobs, xobs, yobs, &
784 & grid(ng)%angler, &
785# ifdef MASKING
786 & grid(ng)%rmask, &
787# endif
788 & obsangler)
789# endif
790!
791! Radial nonlinear u-component at observation locations.
792!
793 DO k=1,n(ng)
794 DO j=jstr-1,jend+1
795 DO i=istru-1,iend+1
796 grid(ng)%z_v(i,j,k)=0.5_r8*(grid(ng)%z_r(i-1,j,k)+ &
797 & grid(ng)%z_r(i ,j,k))
798 END DO
799 END DO
800 END DO
801 CALL extract_obs3d (ng, 1, lm(ng)+1, 0, mm(ng)+1, &
802 & lbi, ubi, lbj, ubj, 1, n(ng), &
803 & obsstate2type(isradial), &
804 & mobs, mstr, mend, &
805 & uxmin(ng), uxmax(ng), &
806 & uymin(ng), uymax(ng), &
807 & time(ng), dt(ng), &
808 & obstype, obsvetting, &
809 & tobs, xobs, yobs, zobs, &
810 & ocean(ng)%u(:,:,:,nout), &
811 & grid(ng)%z_v, &
812# ifdef MASKING
813 & grid(ng)%umask, &
814# endif
815 & uradial)
816
817# ifndef VERIFICATION
818!
819! Radial u-component backgound error standard deviation at observation
820! locations.
821!
822 CALL extract_obs3d (ng, 1, lm(ng)+1, 0, mm(ng)+1, &
823 & lbi, ubi, lbj, ubj, 1, n(ng), &
824 & obsstate2type(isradial), &
825 & mobs, mstr, mend, &
826 & uxmin(ng), uxmax(ng), &
827 & uymin(ng), uymax(ng), &
828 & time(ng), dt(ng), &
829 & obstype, obsvetting, &
830 & tobs, xobs, yobs, zobs, &
831 & ocean(ng)%e_u(:,:,:,1), &
832 & grid(ng)%z_v, &
833# ifdef MASKING
834 & grid(ng)%umask, &
835# endif
836 & ubgerr)
837# endif
838!
839! Radial nonlinear v-component at observation locations.
840!
841 DO k=1,n(ng)
842 DO j=jstrv-1,jend+1
843 DO i=istr-1,iend+1
844 grid(ng)%z_v(i,j,k)=0.5_r8*(grid(ng)%z_r(i,j-1,k)+ &
845 & grid(ng)%z_r(i,j ,k))
846 END DO
847 END DO
848 END DO
849 CALL extract_obs3d (ng, 0, lm(ng)+1, 1, mm(ng)+1, &
850 & lbi, ubi, lbj, ubj, 1, n(ng), &
851 & obsstate2type(isradial), &
852 & mobs, mstr, mend, &
853 & vxmin(ng), vxmax(ng), &
854 & vymin(ng), vymax(ng), &
855 & time(ng), dt(ng), &
856 & obstype, obsvetting, &
857 & tobs, xobs, yobs, zobs, &
858 & ocean(ng)%v(:,:,:,nout), &
859 & grid(ng)%z_v, &
860# ifdef MASKING
861 & grid(ng)%vmask, &
862# endif
863 & vradial)
864
865# ifndef VERIFICATION
866!
867! Radial v-component backgound error standard deviation at observation
868! locations.
869!
870 CALL extract_obs3d (ng, 0, lm(ng)+1, 1, mm(ng)+1, &
871 & lbi, ubi, lbj, ubj, 1, n(ng), &
872 & obsstate2type(isradial), &
873 & mobs, mstr, mend, &
874 & vxmin(ng), vxmax(ng), &
875 & vymin(ng), vymax(ng), &
876 & time(ng), dt(ng), &
877 & obstype, obsvetting, &
878 & tobs, xobs, yobs, zobs, &
879 & ocean(ng)%e_v(:,:,:,1), &
880 & grid(ng)%z_v, &
881# ifdef MASKING
882 & grid(ng)%vmask, &
883# endif
884 & vbgerr)
885# endif
886!
887! Compute nonlinear radial velocity.
888!
889 DO iobs=mstr,mend
890 IF (obstype(iobs).eq.obsstate2type(isradial)) THEN
891# ifdef RADIAL_ANGLE_CCW_EAST
892# ifdef CURVGRID
893 angle=obsmeta(iobs)-obsangler(iobs)
894 nlmodval(iobs)=uradial(iobs)*cos(angle)+ &
895 & vradial(iobs)*sin(angle)
896# else
897 nlmodval(iobs)=uradial(iobs)*cos(obsmeta(iobs))+ &
898 & vradial(iobs)*sin(obsmeta(iobs))
899# endif
900# else
901# ifdef CURVGRID
902 angle=obsmeta(iobs)+obsangler(iobs)
903 nlmodval(iobs)=uradial(iobs)*sin(angle)+ &
904 & vradial(iobs)*cos(angle)
905# else
906 nlmodval(iobs)=uradial(iobs)*sin(obsmeta(iobs))+ &
907 & vradial(iobs)*cos(obsmeta(iobs))
908# endif
909# endif
910# ifndef VERIFICATION
911 bgerr(iobs)=max(ubgerr(iobs), vbgerr(iobs))
912# endif
913 END IF
914 END DO
915 END IF
916# endif
917
918# ifdef TLM_OBS
919!
920! Tangent linear radial velocities.
921!
922 IF ((wrttlmod(ng).or.(wrtrpmod(ng))).and. &
923 & (fourdvar(ng)%ObsCount(isradial).gt.0)) THEN
924 DO iobs=mstr,mend
925 uradial(iobs)=inival
926 vradial(iobs)=inival
927 END DO
928
929# ifdef CURVGRID
930!
931! If curvilinear coordinates, interpolate grid rotation angle (radians)
932! for each observation location.
933!
934 CALL extract_obs2d (ng, 0, lm(ng)+1, 0, mm(ng)+1, &
935 & lbi, ubi, lbj, ubj, &
936 & obsstate2type(isradial), &
937 & mobs, mstr, mend, &
938 & rxmin(ng), rxmax(ng), &
939 & rymin(ng), rymax(ng), &
940 & time(ng), dt(ng), &
941 & obstype, obsvetting, &
942 & tobs, xobs, yobs, &
943 & grid(ng)%angler, &
944# ifdef MASKING
945 & grid(ng)%rmask, &
946# endif
947 & obsangler)
948# endif
949!
950! Tangent linear radial u-component at the observation locations.
951!
952 DO k=1,n(ng)
953 DO j=jstr-1,jend+1
954 DO i=istru-1,iend+1
955 grid(ng)%z_v(i,j,k)=0.5_r8*(grid(ng)%z_r(i-1,j,k)+ &
956 & grid(ng)%z_r(i ,j,k))
957 END DO
958 END DO
959 END DO
960 CALL extract_obs3d (ng, 1, lm(ng)+1, 0, mm(ng)+1, &
961 & lbi, ubi, lbj, ubj, 1, n(ng), &
962 & obsstate2type(isradial), &
963 & mobs, mstr, mend, &
964 & uxmin(ng), uxmax(ng), &
965 & uymin(ng), uymax(ng), &
966 & time(ng), dt(ng), &
967 & obstype, obsvetting, &
968 & tobs, xobs, yobs, zobs, &
969 & ocean(ng)%tl_u(:,:,:,nout), &
970 & grid(ng)%z_v, &
971# ifdef MASKING
972 & grid(ng)%umask, &
973# endif
974 & uradial)
975!
976! Tangent linear radial v-component at the observation locations.
977!
978 DO k=1,n(ng)
979 DO j=jstrv-1,jend+1
980 DO i=istr-1,iend+1
981 grid(ng)%z_v(i,j,k)=0.5_r8*(grid(ng)%z_r(i,j-1,k)+ &
982 & grid(ng)%z_r(i,j ,k))
983 END DO
984 END DO
985 END DO
986 CALL extract_obs3d (ng, 0, lm(ng)+1, 1, mm(ng)+1, &
987 & lbi, ubi, lbj, ubj, 1, n(ng), &
988 & obsstate2type(isradial), &
989 & mobs, mstr, mend, &
990 & vxmin(ng), vxmax(ng), &
991 & vymin(ng), vymax(ng), &
992 & time(ng), dt(ng), &
993 & obstype, obsvetting, &
994 & tobs, xobs, yobs, zobs, &
995 & ocean(ng)%tl_v(:,:,:,nout), &
996 & grid(ng)%z_v, &
997# ifdef MASKING
998 & grid(ng)%vmask, &
999# endif
1000 & vradial)
1001!
1002! Compute tangent linear radial velocity.
1003!
1004 DO iobs=mstr,mend
1005 IF (obstype(iobs).eq.obsstate2type(isradial)) THEN
1006# ifdef RADIAL_ANGLE_CCW_EAST
1007# ifdef CURVGRID
1008 angle=obsmeta(iobs)-obsangler(iobs)
1009 tlmodval(iobs)=uradial(iobs)*cos(angle)+ &
1010 & vradial(iobs)*sin(angle)
1011# else
1012 tlmodval(iobs)=uradial(iobs)*cos(obsmeta(iobs))+ &
1013 & vradial(iobs)*sin(obsmeta(iobs))
1014# endif
1015# else
1016# ifdef CURVGRID
1017 angle=obsmeta(iobs)+obsangler(iobs)
1018 tlmodval(iobs)=uradial(iobs)*sin(angle)+ &
1019 & vradial(iobs)*cos(angle)
1020# else
1021 tlmodval(iobs)=uradial(iobs)*sin(obsmeta(iobs))+ &
1022 & vradial(iobs)*cos(obsmeta(iobs))
1023# endif
1024# endif
1025 END IF
1026 END DO
1027 END IF
1028# endif
1029!
1030!-----------------------------------------------------------------------
1031! Tracers variable observations.
1032!-----------------------------------------------------------------------
1033!
1034 DO itrc=1,nt(ng)
1035
1036# ifndef I4DVAR_ANA_SENSITIVITY
1037!
1038! Nonlinear tracer at the observation locations.
1039
1040 IF (wrtnlmod(ng).and. &
1041 & (fourdvar(ng)%ObsCount(istvar(itrc)).gt.0)) THEN
1042 CALL extract_obs3d (ng, 0, lm(ng)+1, 0, mm(ng)+1, &
1043 & lbi, ubi, lbj, ubj, 1, n(ng), &
1044 & obsstate2type(istvar(itrc)), &
1045 & mobs, mstr, mend, &
1046 & rxmin(ng), rxmax(ng), &
1047 & rymin(ng), rymax(ng), &
1048 & time(ng), dt(ng), &
1049 & obstype, obsvetting, &
1050 & tobs, xobs, yobs, zobs, &
1051 & ocean(ng)%t(:,:,:,nout,itrc), &
1052 & grid(ng)%z_r, &
1053# ifdef MASKING
1054 & grid(ng)%rmask, &
1055# endif
1056 & nlmodval)
1057
1058# ifndef VERIFICATION
1059!
1060! Tracer backgound error standard deviation at observation locations.
1061!
1062 CALL extract_obs3d (ng, 0, lm(ng)+1, 0, mm(ng)+1, &
1063 & lbi, ubi, lbj, ubj, 1, n(ng), &
1064 & obsstate2type(istvar(itrc)), &
1065 & mobs, mstr, mend, &
1066 & rxmin(ng), rxmax(ng), &
1067 & rymin(ng), rymax(ng), &
1068 & time(ng), dt(ng), &
1069 & obstype, obsvetting, &
1070 & tobs, xobs, yobs, zobs, &
1071 & ocean(ng)%e_t(:,:,:,1,itrc), &
1072 & grid(ng)%z_r, &
1073# ifdef MASKING
1074 & grid(ng)%rmask, &
1075# endif
1076 & bgerr)
1077# endif
1078 END IF
1079# endif
1080
1081# ifdef TLM_OBS
1082!
1083! Tangent linear tracer at the observation locations.
1084!
1085 IF ((wrttlmod(ng).or.(wrtrpmod(ng))).and. &
1086 & (fourdvar(ng)%ObsCount(istvar(itrc)).gt.0)) THEN
1087 CALL extract_obs3d (ng, 0, lm(ng)+1, 0, mm(ng)+1, &
1088 & lbi, ubi, lbj, ubj, 1, n(ng), &
1089 & obsstate2type(istvar(itrc)), &
1090 & mobs, mstr, mend, &
1091 & rxmin(ng), rxmax(ng), &
1092 & rymin(ng), rymax(ng), &
1093 & time(ng), dt(ng), &
1094 & obstype, obsvetting, &
1095 & tobs, xobs, yobs, zobs, &
1096 & ocean(ng)%tl_t(:,:,:,nout,itrc), &
1097 & grid(ng)%z_r, &
1098# ifdef MASKING
1099 & grid(ng)%rmask, &
1100# endif
1101 & tlmodval)
1102 END IF
1103# endif
1104 END DO
1105# endif
1106!
1107!-----------------------------------------------------------------------
1108! If processing nonlinear trajectory, load observation reject/accept
1109! flag into screening variable. This needs to be done once and it is
1110! controlled by the main driver. Notice that in routine "extract_obs"
1111! the observations are only accepted if they are at the appropriate
1112! time window and grid bounded at water points.
1113!-----------------------------------------------------------------------
1114!
1115 IF (wrtobsscale(ng).and.wrtnlmod(ng)) THEN
1116 DO iobs=mstr,mend
1117 obsscale(iobs)=obsvetting(iobs)
1118 END DO
1119 END IF
1120
1121# if defined BGQC && !defined I4DVAR_ANA_SENSITIVITY
1122!
1123!-----------------------------------------------------------------------
1124! Perform the background quality control: check and reject observations
1125! that fail. Only observations that are bounded in time and space
1126! (ObsScale .ne. 0) are considered.
1127!-----------------------------------------------------------------------
1128!
1129 IF (wrtobsscale(ng).and.wrtnlmod(ng)) THEN
1130 DO iobs=mstr,mend
1131 IF (obsscale(iobs).ne.inival) THEN
1132# if defined I4DVAR
1133 df1=(obsval(iobs)-nlmodval(iobs))/bgerr(iobs)
1134 df2=(1.0_r8/obserr(iobs))/(bgerr(iobs)*bgerr(iobs))
1135# elif defined WEAK_CONSTRAINT
1136# ifdef R4DVAR
1137 df1=(obsval(iobs)-tlmodval(iobs))/bgerr(iobs)
1138 df2=obserr(iobs)/(bgerr(iobs)*bgerr(iobs))
1139# else
1140 df1=(obsval(iobs)-nlmodval(iobs))/bgerr(iobs)
1141 df2=obserr(iobs)/(bgerr(iobs)*bgerr(iobs))
1142# endif
1143# endif
1144 thresh=bgthresh(iobs)*(1.0_r8+df2)
1145 IF (df1*df1.gt.thresh) THEN
1146 obsscale(iobs)=0.0_r8
1147 END IF
1148 END IF
1149 END DO
1150 END IF
1151# endif
1152
1153# ifdef DISTRIBUTE
1154!
1155!-----------------------------------------------------------------------
1156! Collect screening variable for all extracted data.
1157!-----------------------------------------------------------------------
1158!
1159# ifdef WEAK_CONSTRAINT
1160 ncollect=mend-mstr+1
1161# else
1162 ncollect=mobs
1163# endif
1164 IF (wrtobsscale(ng).and.wrtnlmod(ng)) THEN
1165 CALL mp_collect (ng, model, ncollect, inival, &
1166# ifdef WEAK_CONSTRAINT
1167 & obsscale(mstr:))
1168# else
1169 & obsscale)
1170# endif
1171 END IF
1172# endif
1173
1174# ifndef WEAK_CONSTRAINT
1175!
1176!-----------------------------------------------------------------------
1177! Save observation reject/accept flag into GLOBAL screening variable.
1178!-----------------------------------------------------------------------
1179!
1180! In the primal formulation, the current time window (or survey)
1181! observations are loaded to the working arrays using local indexes
1182! (array elements 1:Nobs) as opposed to global indexes in the dual
1183! formulation. Recall that the screening variable ObsScale is computed
1184! only once to facilitate Background Quality Control on the first pass
1185! of the NLM (WrtObsScale=T and wrtNLmod=T). Therefore, we need to
1186! load and save values into global array ObsScaleGlobal so it can be
1187! used correctly by the TLM, RPM, and ADM.
1188!
1189 IF (wrtobsscale(ng).and.wrtnlmod(ng)) THEN
1190 ic=0
1191 DO iobs=nstrobs(ng),nendobs(ng)
1192 ic=ic+1
1193 obsscaleglobal(iobs)=obsscale(ic)
1194 END DO
1195 ELSE
1196 ic=0
1197 DO iobs=nstrobs(ng),nendobs(ng)
1198 ic=ic+1
1199 obsscale(ic)=obsscaleglobal(iobs)
1200 END DO
1201 END IF
1202# endif
1203!
1204!-----------------------------------------------------------------------
1205! Scale extracted data at observations location. Save NLM unvetted
1206! values before they are rejected by quality control. The unvetted
1207! vector is used only for output purposes to analyze why a datum was
1208! rejected.
1209!-----------------------------------------------------------------------
1210!
1211! Nonlinear state vector.
1212!
1213 IF (wrtnlmod(ng)) THEN
1214 DO iobs=mstr,mend
1215 unvetted(iobs)=nlmodval(iobs)
1216 nlmodval(iobs)=obsscale(iobs)*nlmodval(iobs)
1217 END DO
1218 END IF
1219
1220# ifdef TLM_OBS
1221!
1222! Tangent linear state vector.
1223!
1224 IF (wrttlmod(ng).or.wrtrpmod(ng)) THEN
1225 DO iobs=mstr,mend
1226# ifdef I4DVAR_ANA_SENSITIVITY
1227 tlmodval(iobs)=obsscale(iobs)*tlmodval(iobs)*obserr(iobs)
1228# else
1229 tlmodval(iobs)=obsscale(iobs)*tlmodval(iobs)
1230# endif
1231 END DO
1232 END IF
1233# endif
1234
1235# ifdef DISTRIBUTE
1236!
1237!-----------------------------------------------------------------------
1238! Collect extracted data.
1239!-----------------------------------------------------------------------
1240!
1241# ifndef I4DVAR_ANA_SENSITIVITY
1242!
1243! Nonlinear state vector.
1244!
1245 IF (wrtnlmod(ng)) THEN
1246 CALL mp_collect (ng, model, ncollect, inival, &
1247# if defined WEAK_CONSTRAINT
1248 & nlmodval(mstr:))
1249# else
1250 & nlmodval)
1251# endif
1252!
1253 CALL mp_collect (ng, model, ncollect, inival, &
1254# if defined WEAK_CONSTRAINT
1255 & unvetted(mstr:))
1256# else
1257 & unvetted)
1258# endif
1259 END IF
1260
1261# ifndef VERIFICATION
1262!
1263! Background error standard deviation state vector.
1264!
1265 IF (wrtobsscale(ng).and.wrtnlmod(ng)) THEN
1266 CALL mp_collect (ng, model, ncollect, inival, &
1267# if defined WEAK_CONSTRAINT
1268 & bgerr(mstr:))
1269# else
1270 & bgerr)
1271# endif
1272 END IF
1273# endif
1274# endif
1275
1276# ifdef TLM_OBS
1277!
1278! Tangent linear state vector.
1279!
1280 IF (wrttlmod(ng).or.wrtrpmod(ng)) THEN
1281 CALL mp_collect (ng, model, ncollect, inival, &
1282# if defined WEAK_CONSTRAINT
1283 & tlmodval(mstr:))
1284# else
1285 & tlmodval)
1286# endif
1287 END IF
1288# endif
1289
1290# ifdef SOLVE3D
1291!
1292! State vector depths.
1293!
1294 IF (load_zobs(ng)) THEN
1295 CALL mp_collect (ng, model, ncollect, inival, &
1296# ifdef WEAK_CONSTRAINT
1297 & zobs(mstr:))
1298# else
1299 & zobs)
1300# endif
1301 END IF
1302# endif
1303# endif
1304
1305 RETURN

References mod_fourdvar::bgerr, mod_scalars::bgqc_large, mod_scalars::bgqc_type, mod_fourdvar::bgthresh, mod_scalars::dt, extract_obs_mod::extract_obs2d(), extract_obs_mod::extract_obs3d(), mod_fourdvar::fourdvar, mod_grid::grid, mod_scalars::iprovenance, mod_ncparam::isfsur, mod_ncparam::isradial, mod_ncparam::istvar, mod_ncparam::isubar, mod_ncparam::isuvel, mod_ncparam::isvbar, mod_ncparam::isvvel, mod_param::lm, mod_fourdvar::load_zobs, mod_param::mm, mod_fourdvar::mobs, mod_scalars::mstatevar, mod_param::n, mod_fourdvar::nendobs, mod_fourdvar::nlmodval, mod_fourdvar::nobs, mod_scalars::nprovenance, mod_fourdvar::nstrobs, mod_param::nt, mod_fourdvar::obsangler, mod_fourdvar::obserr, mod_fourdvar::obsmeta, mod_fourdvar::obsprov, mod_fourdvar::obsscale, mod_fourdvar::obsstate2type, mod_fourdvar::obstype, mod_fourdvar::obsval, mod_fourdvar::obsvetting, mod_ocean::ocean, mod_scalars::p_bgqc, mod_ncparam::rxmax, mod_ncparam::rxmin, mod_ncparam::rymax, mod_ncparam::rymin, mod_scalars::s_bgqc, mod_iounits::sourcefile, mod_scalars::time, mod_fourdvar::tobs, mod_fourdvar::unvetted, mod_fourdvar::uradial, mod_ncparam::uxmax, mod_ncparam::uxmin, mod_ncparam::uymax, mod_ncparam::uymin, mod_fourdvar::vradial, mod_ncparam::vxmax, mod_ncparam::vxmin, mod_ncparam::vymax, mod_ncparam::vymin, mod_fourdvar::wrtnlmod, mod_fourdvar::wrtobsscale, mod_fourdvar::wrtrpmod, mod_fourdvar::wrttlmod, mod_fourdvar::xobs, mod_fourdvar::yobs, and mod_fourdvar::zobs.

Referenced by obs_write().

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

◆ obs_write()

subroutine, public obs_write_mod::obs_write ( integer, intent(in) ng,
integer, intent(in) tile,
integer, intent(in) model )

Definition at line 55 of file obs_write.F.

56!***********************************************************************
57!
58! Imported variable declarations.
59!
60 integer, intent(in) :: ng, tile, model
61!
62! Local variable declarations.
63!
64 integer :: Mstr, Mend
65 integer :: LBi, UBi, LBj, UBj
66 integer :: ObsSum, ObsVoid
67 integer :: i, ie, is, iobs, itrc
68!
69 character (len=50) :: string
70 character (len=*), parameter :: MyFile = &
71 & __FILE__
72!
73!-----------------------------------------------------------------------
74! Extract state vector at the observation locations in space and
75! time.
76!-----------------------------------------------------------------------
77!
78 h_operator : IF (processobs(ng)) THEN
79 lbi=bounds(ng)%LBi(tile)
80 ubi=bounds(ng)%UBi(tile)
81 lbj=bounds(ng)%LBj(tile)
82 ubj=bounds(ng)%UBj(tile)
83!
84 CALL obs_operator (ng, tile, model, &
85 & lbi, ubi, lbj, ubj, &
86 & mstr, mend)
87!
88!-----------------------------------------------------------------------
89! Write out state vector to output Data Assimilation Variables (DAV)
90! NetCDF file.
91!-----------------------------------------------------------------------
92!
93 SELECT CASE (dav(ng)%IOtype)
94 CASE (io_nf90)
95 CALL obs_write_nf90 (ng, tile, model, &
96 & lbi, ubi, lbj, ubj, &
97 & mstr, mend)
98
99# if defined PIO_LIB && defined DISTRIBUTE
100 CASE (io_pio)
101 CALL obs_write_pio (ng, tile, model, &
102 & lbi, ubi, lbj, ubj, &
103 & mstr, mend)
104# endif
105 CASE DEFAULT
106 IF (master) WRITE (stdout,10) dav(ng)%IOtype
107 exit_flag=3
108 END SELECT
109 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
110
111!
112!-----------------------------------------------------------------------
113! Set counters for the number of rejected observations for each state
114! variable.
115!-----------------------------------------------------------------------
116!
117 DO iobs=mstr,mend
118 IF (obsscale(iobs).lt.1.0) THEN
119 IF (obstype(iobs).eq.obsstate2type(isfsur)) THEN
120 fourdvar(ng)%ObsReject(isfsur)= &
121 & fourdvar(ng)%ObsReject(isfsur)+1
122 ELSE IF (obstype(iobs).eq.obsstate2type(isubar)) THEN
123 fourdvar(ng)%ObsReject(isubar)= &
124 & fourdvar(ng)%ObsReject(isubar)+1
125 ELSE IF (obstype(iobs).eq.obsstate2type(isvbar)) THEN
126 fourdvar(ng)%ObsReject(isvbar)= &
127 & fourdvar(ng)%ObsReject(isvbar)+1
128# ifdef SOLVE3D
129 ELSE IF (obstype(iobs).eq.obsstate2type(isuvel)) THEN
130 fourdvar(ng)%ObsReject(isuvel)= &
131 & fourdvar(ng)%ObsReject(isuvel)+1
132 ELSE IF (obstype(iobs).eq.obsstate2type(isvvel)) THEN
133 fourdvar(ng)%ObsReject(isvvel)= &
134 & fourdvar(ng)%ObsReject(isvvel)+1
135 ELSE IF (obstype(iobs).eq.obsstate2type(isradial)) THEN
136 fourdvar(ng)%ObsReject(isradial)= &
137 & fourdvar(ng)%ObsReject(isradial)+1
138 ELSE
139 DO itrc=1,nt(ng)
140 IF (obstype(iobs).eq.obsstate2type(istvar(itrc))) THEN
141 i=istvar(itrc)
142 fourdvar(ng)%ObsReject(i)=fourdvar(ng)%ObsReject(i)+1
143 END IF
144 END DO
145# endif
146 END IF
147 END IF
148 END DO
149!
150! Load total available and rejected observations into structure
151! array.
152!
153 DO i=1,nobsvar(ng)
154 fourdvar(ng)%ObsCount(0)=fourdvar(ng)%ObsCount(0)+ &
155 & fourdvar(ng)%ObsCount(i)
156 fourdvar(ng)%ObsReject(0)=fourdvar(ng)%ObsReject(0)+ &
157 & fourdvar(ng)%ObsReject(i)
158 END DO
159!
160!-----------------------------------------------------------------------
161! Report observation processing information.
162!-----------------------------------------------------------------------
163!
164 IF (master) THEN
165 obssum=0
166 obsvoid=0
167 is=nstrobs(ng)
168 DO i=1,nobsvar(ng)
169 IF (fourdvar(ng)%ObsCount(i).gt.0) THEN
170 ie=is+fourdvar(ng)%ObsCount(i)-1
171 WRITE (stdout,20) trim(obsname(i)), is, ie, &
172 & ie-is+1, fourdvar(ng)%ObsReject(i)
173 is=ie+1
174 obssum=obssum+fourdvar(ng)%ObsCount(i)
175 obsvoid=obsvoid+fourdvar(ng)%ObsReject(i)
176 END IF
177 END DO
178 WRITE (stdout,30) obssum, obsvoid, &
179 & fourdvar(ng)%ObsCount(0), &
180 & fourdvar(ng)%ObsReject(0)
181 END IF
182!
183 IF (wrtnlmod(ng)) THEN
184 string='Wrote NLM state at observation locations, '
185# ifdef TLM_OBS
186# ifdef WEAK_CONSTRAINT
187 ELSE IF (wrttlmod(ng)) THEN
188 string='Wrote TLM state at observation locations, '
189 ELSE IF (wrtrpmod(ng)) THEN
190 string='Wrote RPM state at observation locations, '
191# else
192# ifdef I4DVAR_ANA_SENSITIVITY
193 ELSE IF (wrttlmod(ng)) THEN
194 string='Wrote 4DVAR observation sensitivity, '
195# else
196 ELSE IF (wrttlmod(ng)) THEN
197 string='Wrote TLM increments at observation locations, '
198# endif
199# endif
200# endif
201 END IF
202 IF (master) THEN
203 IF (wrtnlmod(ng).or.wrttlmod(ng).or.wrtrpmod(ng)) THEN
204 WRITE (stdout,40) trim(string), nstrobs(ng), nendobs(ng)
205 END IF
206 END IF
207
208 END IF h_operator
209!
210 10 FORMAT (' OBS_WRITE - Illegal output file type, io_type = ',i0, &
211 & /,13x,'Check KeyWord ''OUT_LIB'' in ''roms.in''.')
212 20 FORMAT (10x,a,t25,4(1x,i10))
213 30 FORMAT (/,10x,'Total',t47,2(1x,i10), &
214 & /,10x,'Obs Tally',t47,2(1x,i10),/)
215 40 FORMAT (1x,a,' datum = ',i0,' - ',i0,/)
216!
217 RETURN

References mod_param::bounds, mod_iounits::dav, mod_scalars::exit_flag, strings_mod::founderror(), mod_fourdvar::fourdvar, mod_ncparam::io_nf90, mod_ncparam::io_pio, mod_ncparam::isfsur, mod_ncparam::isradial, mod_ncparam::istvar, mod_ncparam::isubar, mod_ncparam::isuvel, mod_ncparam::isvbar, mod_ncparam::isvvel, mod_parallel::master, mod_fourdvar::nendobs, mod_fourdvar::nobsvar, mod_scalars::noerror, mod_fourdvar::nstrobs, mod_param::nt, obs_operator(), obs_write_nf90(), obs_write_pio(), mod_fourdvar::obsname, mod_fourdvar::obsscale, mod_fourdvar::obsstate2type, mod_fourdvar::obstype, mod_fourdvar::processobs, mod_iounits::stdout, mod_fourdvar::wrtnlmod, mod_fourdvar::wrtrpmod, and mod_fourdvar::wrttlmod.

Referenced by output(), rp_output(), and tl_output().

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

◆ obs_write_nf90()

subroutine, private obs_write_mod::obs_write_nf90 ( integer, intent(in) ng,
integer, intent(in) tile,
integer, intent(in) model,
integer, intent(in) lbi,
integer, intent(in) ubi,
integer, intent(in) lbj,
integer, intent(in) ubj,
integer, intent(in) mstr,
integer, intent(in) mend )
private

Definition at line 1309 of file obs_write.F.

1312!***********************************************************************
1313!
1314 USE mod_netcdf
1315!
1316! Imported variable declarations.
1317!
1318 integer, intent(in) :: ng, tile, model
1319 integer, intent(in) :: LBi, UBi, LBj, UBj
1320 integer, intent(in) :: Mstr, Mend
1321!
1322! Local variable declarations.
1323!
1324 logical :: Lwrote
1325!
1326 integer :: Tindex, iobs, status
1327!
1328 real(dp) :: scale
1329!
1330 character (len=*), parameter :: MyFile = &
1331 & __FILE__//", obs_write_nf90"
1332!
1333 sourcefile=myfile
1334
1335# if defined FOUR_DVAR && !defined I4DVAR_ANA_SENSITIVITY
1336!
1337!-----------------------------------------------------------------------
1338! Compute and write initial and final model-observation misfit
1339! (innovation) vector for output purposes only. Write also initial
1340! nonlinear model at observation locations.
1341!-----------------------------------------------------------------------
1342!
1343 IF (wrtmisfit(ng)) THEN
1344 DO iobs=mstr,mend
1345# if defined I4DVAR
1346 misfit(iobs)=obsscale(iobs)*sqrt(obserr(iobs))* &
1347 & (nlmodval(iobs)-obsval(iobs))
1348# elif defined WEAK_CONSTRAINT
1349# ifdef R4DVAR
1350 misfit(iobs)=obsscale(iobs)/sqrt(obserr(iobs))* &
1351 & (tlmodval(iobs)-obsval(iobs))
1352# else
1353 misfit(iobs)=obsscale(iobs)/sqrt(obserr(iobs))* &
1354 & (nlmodval(iobs)-obsval(iobs))
1355# endif
1356# endif
1357 END DO
1358!
1359! Write out initial or final values into DAV NetCDF file.
1360!
1361 IF (nrun.eq.1) THEN
1362 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1363 & vname(1,idnlmp), &
1364 & unvetted(mstr:), &
1365 & (/nstrobs(ng)/), (/nobs(ng)/), &
1366 & ncid = dav(ng)%ncid, &
1367 & varid = dav(ng)%Vid(idnlmp))
1368 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1369!
1370 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1371 & vname(1,idnlmi), &
1372 & nlmodval(mstr:), &
1373 & (/nstrobs(ng)/), (/nobs(ng)/), &
1374 & ncid = dav(ng)%ncid, &
1375 & varid = dav(ng)%Vid(idnlmi))
1376 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1377!
1378 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1379 & vname(1,idmomi), &
1380 & misfit(mstr:), &
1381 & (/nstrobs(ng)/), (/nobs(ng)/), &
1382 & ncid = dav(ng)%ncid, &
1383 & varid = dav(ng)%Vid(idmomi))
1384 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1385 ELSE
1386 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1387 & vname(1,idmomf), &
1388 & misfit(mstr:), &
1389 & (/nstrobs(ng)/), (/nobs(ng)/), &
1390 & ncid = dav(ng)%ncid, &
1391 & varid = dav(ng)%Vid(idmomf))
1392 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1393 END IF
1394 END IF
1395
1396# ifdef BGQC
1397!
1398! Write out bacground error standard deviation.
1399!
1400 IF (wrtobsscale(ng).and.wrtnlmod(ng)) THEN
1401 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1402 & vname(1,idbgth), &
1403 & bgthresh(mstr:), &
1404 & (/nstrobs(ng)/), (/nobs(ng)/), &
1405 & ncid = dav(ng)%ncid, &
1406 & varid = dav(ng)%Vid(idbgth))
1407 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1408 END IF
1409# endif
1410# endif
1411
1412# if defined FOUR_DVAR && !defined I4DVAR_ANA_SENSITIVITY
1413!
1414!-----------------------------------------------------------------------
1415! Write out variables needed to compute Desroziers et al. (2005)
1416! statistics.
1417!-----------------------------------------------------------------------
1418!
1419 IF (wrtnlmod(ng)) THEN
1420!
1421! 4D-Var background error standard deviation at observation locations.
1422!
1423 IF (wrtobsscale(ng)) THEN
1424 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1425 & vname(1,idbger), &
1426 & bgerr(mstr:), &
1427 & (/nstrobs(ng)/), (/nobs(ng)/), &
1428 & ncid = dav(ng)%ncid, &
1429 & varid = dav(ng)%Vid(idbger))
1430 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1431 END IF
1432!
1433! 4D-Var innovation vector: observation minus background. Save initial
1434! background in the NLincrement vector.
1435!
1436 IF (nrun.eq.1) THEN
1437 DO iobs=mstr,mend
1438 innovation(iobs)=obsscale(iobs)* &
1439 & (obsval(iobs)-nlmodval(iobs))
1440 nlincrement(iobs)=nlmodval(iobs)
1441 END DO
1442!
1443 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1444 & vname(1,idinno), &
1445 & innovation(mstr:), &
1446 & (/nstrobs(ng)/), (/nobs(ng)/), &
1447 & ncid = dav(ng)%ncid, &
1448 & varid = dav(ng)%Vid(idinno))
1449 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1450 END IF
1451!
1452! 4D-Var increment (analysis minus background) and residual
1453! (observation minus analysis).
1454!
1455 IF (outer.eq.nouter) THEN
1456 DO iobs=mstr,mend
1457 nlincrement(iobs)=obsscale(iobs)* &
1458 & (nlmodval(iobs)-nlincrement(iobs))
1459 residual(iobs)=obsscale(iobs)* &
1460 & (obsval(iobs)-nlmodval(iobs))
1461 END DO
1462!
1463 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1464 & vname(1,idincr), &
1465 & nlincrement(mstr:), &
1466 & (/nstrobs(ng)/), (/nobs(ng)/), &
1467 & ncid = dav(ng)%ncid, &
1468 & varid = dav(ng)%Vid(idincr))
1469 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1470!
1471 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1472 & vname(1,idresi), &
1473 & residual(mstr:), &
1474 & (/nstrobs(ng)/), (/nobs(ng)/), &
1475 & ncid = dav(ng)%ncid, &
1476 & varid = dav(ng)%Vid(idresi))
1477 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1478 END IF
1479 END IF
1480# endif
1481!
1482!-----------------------------------------------------------------------
1483! Write out other variable into output DAC NetCDF file.
1484!-----------------------------------------------------------------------
1485
1486# if defined I4DVAR || defined WEAK_CONSTRAINT && \
1487 !defined I4DVAR_ANA_SENSITIVITY
1488!
1489! Current outer and inner loop.
1490!
1491 IF (wrtnlmod(ng).or.wrttlmod(ng).or.wrtrpmod(ng)) THEN
1492 CALL netcdf_put_ivar (ng, model, dav(ng)%name, 'outer', &
1493 & outer, (/0/), (/0/), &
1494 & ncid = dav(ng)%ncid)
1495 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1496!
1497 CALL netcdf_put_ivar (ng, model, dav(ng)%name, 'inner', &
1498 & inner, (/0/), (/0/), &
1499 & ncid = dav(ng)%ncid)
1500 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1501 END IF
1502#endif
1503!
1504! Observation screening flag: accept or reject.
1505!
1506 IF (wrtobsscale(ng).and.wrtnlmod(ng)) THEN
1507 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1508 & vname(1,idobss), &
1509 & obsscale(mstr:), &
1510 & (/nstrobs(ng)/), (/nobs(ng)/), &
1511 & ncid = dav(ng)%ncid, &
1512 & varid = dav(ng)%Vid(idobss))
1513 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1514 END IF
1515
1516# ifndef I4DVAR_ANA_SENSITIVITY
1517!
1518! Nonlinear model or first guess (background) state at observation
1519! locations. If 4D-Var, the values are written for every outer loop.
1520! Also, the final value is written as individual variable to
1521! facilitate ERDDAP processing. Notice that it is always overwritten
1522! to avoid complicated logic with all the 4D-Var drivers. The last
1523! one written correspont to the final values.
1524!
1525# if defined VERIFICATION || defined TLM_CHECK
1526 IF (wrtnlmod(ng)) THEN
1527 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1528 & vname(1,idnlmo), &
1529 & nlmodval(mstr:), &
1530 & (/nstrobs(ng)/), (/nobs(ng)/), &
1531 & ncid = dav(ng)%ncid, &
1532 & varid = dav(ng)%Vid(idnlmo))
1533 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1534 havenlmod(ng)=.true.
1535 END IF
1536# elif defined I4DVAR || defined WEAK_CONSTRAINT
1537 IF (wrtnlmod(ng).and.outer.gt.0) THEN
1538 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1539 & vname(1,idnlmo), &
1540 & nlmodval(mstr:), &
1541 & (/nstrobs(ng),outer/), (/nobs(ng),1/), &
1542 & ncid = dav(ng)%ncid, &
1543 & varid = dav(ng)%Vid(idnlmo))
1544 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1545 END IF
1546# endif
1547# ifdef FOUR_DVAR
1548!
1549 IF (wrtnlmod(ng).and.outer.gt.0) THEN
1550 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1551 & vname(1,idnlmf), &
1552 & nlmodval(mstr:), &
1553 & (/nstrobs(ng)/), (/nobs(ng)/), &
1554 & ncid = dav(ng)%ncid, &
1555 & varid = dav(ng)%Vid(idnlmf))
1556 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1557 havenlmod(ng)=.true.
1558 END IF
1559# endif
1560# endif
1561# if !defined I4DVAR_ANA_SENSITIVITY && \
1562 (defined i4dvar || defined weak_constraint)
1563!
1564! Write out unvetted nonlinear model at the observation locations per
1565! outer loop.
1566!
1567 IF (wrtnlmod(ng).and.outer.gt.0) THEN
1568 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1569 & vname(1,idnlmu), &
1570 & unvetted(mstr:), &
1571 & (/nstrobs(ng),outer/), (/nobs(ng),1/), &
1572 & ncid = dav(ng)%ncid, &
1573 & varid = dav(ng)%Vid(idnlmu))
1574 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1575 END IF
1576# endif
1577# if defined I4DVAR || defined I4DVAR_ANA_SENSITIVITY || \
1578 defined weak_constraint
1579!
1580! Tangent linear model state increments at observation locations.
1581!
1582 IF (wrttlmod(ng).or.wrtrpmod(ng)) THEN
1583 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1584 & vname(1,idtlmo), &
1585 & tlmodval(mstr:), &
1586 & (/nstrobs(ng)/), (/nobs(ng)/), &
1587 & ncid = dav(ng)%ncid, &
1588 & varid = dav(ng)%Vid(idtlmo))
1589 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1590 havetlmod(ng)=.true.
1591 END IF
1592# endif
1593# ifdef TLM_OBS
1594# ifdef I4DVAR_ANA_SENSITIVITY
1595# ifdef OBS_IMPACT
1596!
1597! Write the total observation impact.
1598!
1599 IF (wrtimpact_tot(ng)) THEN
1600 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1601 & 'ObsImpact_total', &
1602 & tlmodval(mstr:), &
1603 & (/nstrobs(ng)/), (/nobs(ng)/), &
1604 & ncid = dav(ng)%ncid)
1605 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1606 END IF
1607# endif
1608
1609# ifdef OBS_IMPACT_SPLIT
1610!
1611! Write the observation impact of initial conditions.
1612!
1613 IF (wrtimpact_ic(ng)) THEN
1614 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1615 & 'ObsImpact_IC', &
1616 & tlmodval(mstr:), &
1617 & (/nstrobs(ng)/), (/nobs(ng)/), &
1618 & ncid = dav(ng)%ncid)
1619 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1620 END IF
1621
1622# if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
1623!
1624! Write the observation impact of surface forcing.
1625!
1626 IF (wrtimpact_fc(ng)) THEN
1627 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1628 & 'ObsImpact_FC', &
1629 & tlmodval(mstr:), &
1630 & (/nstrobs(ng)/), (/nobs(ng)/), &
1631 & ncid = dav(ng)%ncid)
1632 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1633 END IF
1634# endif
1635
1636# if defined ADJUST_BOUNDARY
1637!
1638! Write the observation impact of boundary conditions.
1639!
1640 IF (wrtimpact_bc(ng)) THEN
1641 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1642 & 'ObsImpact_BC', &
1643 & tlmodval(mstr:), &
1644 & (/nstrobs(ng)/), (/nobs(ng)/), &
1645 & ncid = dav(ng)%ncid)
1646 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1647 END IF
1648# endif
1649# endif
1650# endif
1651# endif
1652
1653# if defined R4DVAR || defined R4DVAR_ANA_SENSITIVITY || \
1654 defined tl_r4dvar
1655!
1656! Write initial representer model increments at observation locations.
1657!
1658 IF (nrun.eq.1.and.wrtrpmod(ng)) THEN
1659 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1660 & 'RPmodel_initial', &
1661 & tlmodval(mstr:), &
1662 & (/nstrobs(ng)/), (/nobs(ng)/) , &
1663 & ncid = dav(ng)%ncid)
1664 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1665 END IF
1666# endif
1667
1668# ifdef SOLVE3D
1669!
1670! Now, that Zobs has all the Z-locations in grid coordinates, write
1671! "obs_Zgrid" into data assimilation output file.
1672!
1673 IF ((nrun.eq.1).and. &
1674 & (wrtnlmod(ng).or.wrttlmod(ng).or.wrtrpmod(ng))) THEN
1675 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1676 & vname(1,idobsz), &
1677 & zobs(mstr:), &
1678 & (/nstrobs(ng)/), (/nobs(ng)/), &
1679 & ncid = dav(ng)%ncid, &
1680 & varid = dav(ng)%Vid(idobsz))
1681 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1682 END IF
1683!
1684! Write Z-location of observation in grid coordinates, if applicable.
1685! This values are needed elsewhere when using the interpolation
1686! weight matrix. Recall that the depth of observations can be in
1687! meters or grid coordinates. Notice that since the model levels
1688! evolve in time, the fractional level coordinate is unknow during
1689! the processing of the observations.
1690!
1691 IF ((nrun.eq.1).and. &
1692 & (wrtnlmod(ng).or.wrttlmod(ng).or.wrtrpmod(ng))) THEN
1693 CALL netcdf_put_fvar (ng, model, obs(ng)%name, &
1694 & vname(1,idobsz), &
1695 & zobs(mstr:), &
1696 & (/nstrobs(ng)/), (/nobs(ng)/), &
1697 & ncid = obs(ng)%ncid, &
1698 & varid = obs(ng)%Vid(idobsz))
1699 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1700
1701 IF (model.eq.iadm) THEN
1702 lwrote=obssurvey(ng).eq.1
1703 ELSE
1704 lwrote=obssurvey(ng).eq.nsurvey(ng)
1705 END IF
1706 IF (lwrote) wrote_zobs(ng)=.true.
1707 END IF
1708# endif
1709
1710# if defined BALANCE_OPERATOR && defined ZETA_ELLIPTIC
1711!
1712! Write out reference free-surface used in the balance operator.
1713!
1714 IF (wrtzetaref(ng).and.balance(isfsur)) THEN
1715 tindex=0
1716 scale=1.0_dp
1717 status=nf_fwrite2d(ng, model, dav(ng)%ncid, idfsur, &
1718 & dav(ng)%Vid(idfsur), tindex, r2dvar, &
1719 & lbi, ubi, lbj, ubj, scale, &
1720# ifdef MASKING
1721 & grid(ng) % rmask, &
1722# endif
1723 & fourdvar(ng) % zeta_ref, &
1724 & setfillval = .false.)
1725 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1726 wrtzetaref(ng)=.false.
1727 END IF
1728# endif
1729!
1730!-----------------------------------------------------------------------
1731! Synchronize observations NetCDF file to disk.
1732!-----------------------------------------------------------------------
1733!
1734 IF (wrtnlmod(ng).or.wrttlmod(ng).or.wrtrpmod(ng)) THEN
1735 CALL netcdf_sync (ng, model, dav(ng)%name, dav(ng)%ncid)
1736 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1737
1738# ifdef SOLVE3D
1739 IF (nrun.eq.1) THEN
1740 CALL netcdf_sync (ng, model, obs(ng)%name, obs(ng)%ncid)
1741 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1742 END IF
1743# endif
1744 END IF
1745!
1746 RETURN
mod_netcdf::netcdf_put_fvar
Definition mod_netcdf.F:82
mod_netcdf::netcdf_put_ivar
Definition mod_netcdf.F:96
mod_netcdf
Definition mod_netcdf.F:2
mod_netcdf::netcdf_sync
subroutine, public netcdf_sync(ng, model, ncname, ncid)
Definition mod_netcdf.F:9505

References mod_scalars::balance, mod_fourdvar::bgerr, mod_fourdvar::bgthresh, mod_iounits::dav, mod_scalars::exit_flag, strings_mod::founderror(), mod_fourdvar::fourdvar, mod_grid::grid, mod_fourdvar::havenlmod, mod_fourdvar::havetlmod, mod_param::iadm, mod_ncparam::idbger, mod_ncparam::idbgth, mod_ncparam::idfsur, mod_ncparam::idincr, mod_ncparam::idinno, mod_ncparam::idmomf, mod_ncparam::idmomi, mod_ncparam::idnlmf, mod_ncparam::idnlmi, mod_ncparam::idnlmo, mod_ncparam::idnlmp, mod_ncparam::idnlmu, mod_ncparam::idobss, mod_ncparam::idobsz, mod_ncparam::idresi, mod_ncparam::idtlmo, mod_scalars::inner, mod_fourdvar::innovation, mod_ncparam::isfsur, mod_fourdvar::misfit, mod_netcdf::netcdf_sync(), mod_fourdvar::nlincrement, mod_fourdvar::nlmodval, mod_fourdvar::nobs, mod_scalars::noerror, mod_scalars::nouter, mod_scalars::nrun, mod_fourdvar::nstrobs, mod_fourdvar::nsurvey, mod_iounits::obs, mod_fourdvar::obserr, mod_fourdvar::obsscale, mod_fourdvar::obssurvey, mod_fourdvar::obsval, mod_scalars::outer, mod_param::r2dvar, mod_fourdvar::residual, mod_iounits::sourcefile, mod_fourdvar::unvetted, mod_ncparam::vname, mod_fourdvar::wrote_zobs, mod_fourdvar::wrtimpact_bc, mod_fourdvar::wrtimpact_fc, mod_fourdvar::wrtimpact_ic, mod_fourdvar::wrtimpact_tot, mod_fourdvar::wrtmisfit, mod_fourdvar::wrtnlmod, mod_fourdvar::wrtobsscale, mod_fourdvar::wrtrpmod, mod_fourdvar::wrttlmod, mod_fourdvar::wrtzetaref, and mod_fourdvar::zobs.

Referenced by obs_write().

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

◆ obs_write_pio()

subroutine, private obs_write_mod::obs_write_pio ( integer, intent(in) ng,
integer, intent(in) tile,
integer, intent(in) model,
integer, intent(in) lbi,
integer, intent(in) ubi,
integer, intent(in) lbj,
integer, intent(in) ubj,
integer, intent(in) mstr,
integer, intent(in) mend )
private

Definition at line 1752 of file obs_write.F.

1755!***********************************************************************
1756!
1757 USE mod_pio_netcdf
1758!
1759! Imported variable declarations.
1760!
1761 integer, intent(in) :: ng, tile, model
1762 integer, intent(in) :: LBi, UBi, LBj, UBj
1763 integer, intent(in) :: Mstr, Mend
1764!
1765! Local variable declarations.
1766!
1767 logical :: Lwrote
1768!
1769 integer :: Tindex, iobs, status
1770!
1771 real(dp) :: scale
1772!
1773 character (len=*), parameter :: MyFile = &
1774 & __FILE__//", obs_write_pio"
1775
1776# if defined BALANCE_OPERATOR && defined ZETA_ELLIPTIC
1777!
1778 TYPE (IO_desc_t), pointer :: ioDesc
1779# endif
1780!
1781 sourcefile=myfile
1782
1783# if defined FOUR_DVAR && !defined I4DVAR_ANA_SENSITIVITY
1784!
1785!-----------------------------------------------------------------------
1786! Compute and write initial and final model-observation misfit
1787! (innovation) vector for output purposes only. Write also initial
1788! nonlinear model at observation locations.
1789!-----------------------------------------------------------------------
1790!
1791 IF (wrtmisfit(ng)) THEN
1792 DO iobs=mstr,mend
1793# if defined I4DVAR
1794 misfit(iobs)=obsscale(iobs)*sqrt(obserr(iobs))* &
1795 & (nlmodval(iobs)-obsval(iobs))
1796# elif defined WEAK_CONSTRAINT
1797# ifdef R4DVAR
1798 misfit(iobs)=obsscale(iobs)/sqrt(obserr(iobs))* &
1799 & (tlmodval(iobs)-obsval(iobs))
1800# else
1801 misfit(iobs)=obsscale(iobs)/sqrt(obserr(iobs))* &
1802 & (nlmodval(iobs)-obsval(iobs))
1803# endif
1804# endif
1805 END DO
1806!
1807! Write out initial or final values into DAV NetCDF file.
1808!
1809 IF (nrun.eq.1) THEN
1810 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1811 & vname(1,idnlmp), &
1812 & unvetted(mstr:mend), &
1813 & (/nstrobs(ng)/), (/nobs(ng)/), &
1814 & piofile = dav(ng)%pioFile, &
1815 & piovar = dav(ng)%pioVar(idnlmp)%vd)
1816 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1817!
1818 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1819 & vname(1,idnlmi), &
1820 & nlmodval(mstr:mend), &
1821 & (/nstrobs(ng)/), (/nobs(ng)/), &
1822 & piofile = dav(ng)%pioFile, &
1823 & piovar = dav(ng)%pioVar(idnlmi)%vd)
1824 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1825!
1826 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1827 & vname(1,idmomi), &
1828 & misfit(mstr:mend), &
1829 & (/nstrobs(ng)/), (/nobs(ng)/), &
1830 & piofile = dav(ng)%pioFile, &
1831 & piovar = dav(ng)%pioVar(idmomi)%vd)
1832 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1833 ELSE
1834 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1835 & vname(1,idmomf), &
1836 & misfit(mstr:mend), &
1837 & (/nstrobs(ng)/), (/nobs(ng)/), &
1838 & piofile = dav(ng)%pioFile, &
1839 & piovar = dav(ng)%pioVar(idmomf)%vd)
1840 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1841 END IF
1842 END IF
1843
1844# ifdef BGQC
1845!
1846! Write out bacground error standard deviation threshold.
1847!
1848 IF (wrtobsscale(ng).and.wrtnlmod(ng)) THEN
1849 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1850 & vname(1,idbgth), &
1851 & bgthresh(mstr:mend), &
1852 & (/nstrobs(ng)/), (/nobs(ng)/), &
1853 & piofile = dav(ng)%pioFile, &
1854 & piovar = dav(ng)%pioVar(idbgth)%vd)
1855 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1856 END IF
1857# endif
1858# endif
1859
1860# if defined FOUR_DVAR && !defined I4DVAR_ANA_SENSITIVITY
1861!
1862!-----------------------------------------------------------------------
1863! Write out variables needed to compute Desroziers et al. (2005)
1864! statistics.
1865!-----------------------------------------------------------------------
1866!
1867 IF (wrtnlmod(ng)) THEN
1868!
1869! 4D-Var background error standard deviation at observation locations.
1870!
1871 IF (wrtobsscale(ng)) THEN
1872 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1873 & vname(1,idbger), &
1874 & bgerr(mstr:mend), &
1875 & (/nstrobs(ng)/), (/nobs(ng)/), &
1876 & piofile = dav(ng)%pioFile, &
1877 & piovar = dav(ng)%pioVar(idbger)%vd)
1878 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1879 END IF
1880!
1881! 4D-Var innovation vector: observation minus background. Save initial
1882! background in the NLincrement vector.
1883!
1884 IF (nrun.eq.1) THEN
1885 DO iobs=mstr,mend
1886 innovation(iobs)=obsscale(iobs)* &
1887 & (obsval(iobs)-nlmodval(iobs))
1888 nlincrement(iobs)=nlmodval(iobs)
1889 END DO
1890!
1891 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1892 & vname(1,idinno), &
1893 & innovation(mstr:mend), &
1894 & (/nstrobs(ng)/), (/nobs(ng)/), &
1895 & piofile = dav(ng)%pioFile, &
1896 & piovar = dav(ng)%pioVar(idinno)%vd)
1897 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1898 END IF
1899!
1900! 4D-Var increment (analysis minus background) and residual
1901! (observation minus analysis).
1902!
1903 IF (outer.eq.nouter) THEN
1904 DO iobs=mstr,mend
1905 nlincrement(iobs)=obsscale(iobs)* &
1906 & (nlmodval(iobs)-nlincrement(iobs))
1907 residual(iobs)=obsscale(iobs)* &
1908 & (obsval(iobs)-nlmodval(iobs))
1909 END DO
1910!
1911 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1912 & vname(1,idincr), &
1913 & nlincrement(mstr:mend), &
1914 & (/nstrobs(ng)/), (/nobs(ng)/), &
1915 & piofile = dav(ng)%pioFile, &
1916 & piovar = dav(ng)%pioVar(idincr)%vd)
1917 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1918!
1919 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1920 & vname(1,idresi), &
1921 & residual(mstr:mend), &
1922 & (/nstrobs(ng)/), (/nobs(ng)/), &
1923 & piofile = dav(ng)%pioFile, &
1924 & piovar = dav(ng)%pioVar(idresi)%vd)
1925 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1926 END IF
1927 END IF
1928# endif
1929!
1930!-----------------------------------------------------------------------
1931! Write out other variable into output DAC NetCDF file.
1932!-----------------------------------------------------------------------
1933
1934# if defined I4DVAR || defined WEAK_CONSTRAINT && \
1935 !defined I4DVAR_ANA_SENSITIVITY
1936!
1937! Current outer and inner loop.
1938!
1939 IF (wrtnlmod(ng).or.wrttlmod(ng).or.wrtrpmod(ng)) THEN
1940 CALL pio_netcdf_put_ivar (ng, model, dav(ng)%name, 'outer', &
1941 & outer, (/0/), (/0/), &
1942 & piofile = dav(ng)%pioFile)
1943 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1944!
1945 CALL pio_netcdf_put_ivar (ng, model, dav(ng)%name, 'inner', &
1946 & inner, (/0/), (/0/), &
1947 & piofile = dav(ng)%pioFile)
1948 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1949 END IF
1950# endif
1951!
1952! Observation screening flag: accept or reject.
1953!
1954 IF (wrtobsscale(ng).and.wrtnlmod(ng)) THEN
1955 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1956 & vname(1,idobss), &
1957 & obsscale(mstr:mend), &
1958 & (/nstrobs(ng)/), (/nobs(ng)/), &
1959 & piofile = dav(ng)%pioFile, &
1960 & piovar = dav(ng)%pioVar(idobss)%vd)
1961 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1962 END IF
1963
1964# ifndef I4DVAR_ANA_SENSITIVITY
1965!
1966! Nonlinear model or first guess (background) state at observation
1967! locations. If 4D-Var, the values are written for every outer loop.
1968! Also, the final value is written as individual variable to
1969! facilitate ERDDAP processing. Notice that it is always overwritten
1970! to avoid complicated logic with all the 4D-Var drivers. The last
1971! one written correspont to the final values.
1972!
1973# if defined VERIFICATION || defined TLM_CHECK
1974 IF (wrtnlmod(ng)) THEN
1975 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1976 & vname(1,idnlmo), &
1977 & nlmodval(mstr:mend), &
1978 & (/nstrobs(ng)/), (/nobs(ng)/), &
1979 & piofile = dav(ng)%pioFile, &
1980 & piovar = dav(ng)%pioVar(idnlmo)%vd)
1981 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1982 havenlmod(ng)=.true.
1983 END IF
1984# elif defined I4DVAR || defined WEAK_CONSTRAINT
1985 IF (wrtnlmod(ng).and.outer.gt.0) THEN
1986 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1987 & vname(1,idnlmo), &
1988 & nlmodval(mstr:mend), &
1989 & (/nstrobs(ng),outer/), (/nobs(ng),1/),&
1990 & piofile = dav(ng)%pioFile, &
1991 & piovar = dav(ng)%pioVar(idnlmo)%vd)
1992 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1993 END IF
1994# endif
1995# ifdef FOUR_DVAR
1996!
1997 IF (wrtnlmod(ng).and.outer.gt.0) THEN
1998 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1999 & vname(1,idnlmf), &
2000 & nlmodval(mstr:mend), &
2001 & (/nstrobs(ng)/), (/nobs(ng)/), &
2002 & piofile = dav(ng)%pioFile, &
2003 & piovar = dav(ng)%pioVar(idnlmf)%vd)
2004 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2005 havenlmod(ng)=.true.
2006 END IF
2007# endif
2008# endif
2009# if !defined I4DVAR_ANA_SENSITIVITY && \
2010 (defined i4dvar || defined weak_constraint)
2011!
2012! Write out unvetted nonlinear model at the observation locations per
2013! outer loop.
2014!
2015 IF (wrtnlmod(ng).and.outer.gt.0) THEN
2016 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2017 & vname(1,idnlmu), &
2018 & unvetted(mstr:mend), &
2019 & (/nstrobs(ng),outer/), (/nobs(ng),1/),&
2020 & piofile = dav(ng)%pioFile, &
2021 & piovar = dav(ng)%pioVar(idnlmu)%vd)
2022 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2023 END IF
2024# endif
2025# if defined I4DVAR || defined I4DVAR_ANA_SENSITIVITY || \
2026 defined weak_constraint
2027!
2028! Tangent linear model state increments at observation locations.
2029!
2030 IF (wrttlmod(ng).or.wrtrpmod(ng)) THEN
2031 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2032 & vname(1,idtlmo), &
2033 & tlmodval(mstr:mend), &
2034 & (/nstrobs(ng)/), (/nobs(ng)/), &
2035 & piofile = dav(ng)%pioFile, &
2036 & piovar = dav(ng)%pioVar(idtlmo)%vd)
2037 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2038 havetlmod(ng)=.true.
2039 END IF
2040# endif
2041# ifdef TLM_OBS
2042# ifdef I4DVAR_ANA_SENSITIVITY
2043# ifdef OBS_IMPACT
2044!
2045! Write the total observation impact.
2046!
2047 IF (wrtimpact_tot(ng)) THEN
2048 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2049 & 'ObsImpact_total', &
2050 & tlmodval(mstr:mend), &
2051 & (/nstrobs(ng)/), (/nobs(ng)/), &
2052 & piofile = dav(ng)%pioFile)
2053 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2054 END IF
2055# endif
2056
2057# ifdef OBS_IMPACT_SPLIT
2058!
2059! Write the observation impact of initial conditions.
2060!
2061 IF (wrtimpact_ic(ng)) THEN
2062 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2063 & 'ObsImpact_IC', &
2064 & tlmodval(mstr:mend), &
2065 & (/nstrobs(ng)/), (/nobs(ng)/), &
2066 & piofile = dav(ng)%pioFile)
2067 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2068 END IF
2069
2070# if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
2071!
2072! Write the observation impact of surface forcing.
2073!
2074 IF (wrtimpact_fc(ng)) THEN
2075 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2076 & 'ObsImpact_FC', &
2077 & tlmodval(mstr:mend), &
2078 & (/nstrobs(ng)/), (/nobs(ng)/), &
2079 & piofile = dav(ng)%pioFile)
2080 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2081 END IF
2082# endif
2083
2084# if defined ADJUST_BOUNDARY
2085!
2086! Write the observation impact of boundary conditions.
2087!
2088 IF (wrtimpact_bc(ng)) THEN
2089 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2090 & 'ObsImpact_BC', &
2091 & tlmodval(mstr:mend), &
2092 & (/nstrobs(ng)/), (/nobs(ng)/), &
2093 & piofile = dav(ng)%pioFile)
2094 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2095 END IF
2096# endif
2097# endif
2098# endif
2099# endif
2100
2101# if defined R4DVAR || defined R4DVAR_ANA_SENSITIVITY || \
2102 defined tl_r4dvar
2103!
2104! Write initial representer model increments at observation locations.
2105!
2106 IF (nrun.eq.1.and.wrtrpmod(ng)) THEN
2107 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2108 & 'RPmodel_initial', &
2109 & tlmodval(mstr:mend), &
2110 & (/nstrobs(ng)/), (/nobs(ng)/) , &
2111 & piofile = dav(ng)%pioFile)
2112 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2113 END IF
2114# endif
2115
2116# ifdef SOLVE3D
2117!
2118! Now, that Zobs has all the Z-locations in grid coordinates, write
2119! "obs_Zgrid" into data assimilation output file.
2120!
2121 IF ((nrun.eq.1).and. &
2122 & (wrtnlmod(ng).or.wrttlmod(ng).or.wrtrpmod(ng))) THEN
2123 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2124 & vname(1,idobsz), &
2125 & zobs(mstr:mend), &
2126 & (/nstrobs(ng)/), (/nobs(ng)/), &
2127 & piofile = dav(ng)%pioFile, &
2128 & piovar = dav(ng)%pioVar(idobsz)%vd)
2129 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2130 END IF
2131!
2132! Write Z-location of observation in grid coordinates, if applicable.
2133! This values are needed elsewhere when using the interpolation
2134! weight matrix. Recall that the depth of observations can be in
2135! meters or grid coordinates. Notice that since the model levels
2136! evolve in time, the fractional level coordinate is unknow during
2137! the processing of the observations.
2138!
2139 IF ((nrun.eq.1).and. &
2140 & (wrtnlmod(ng).or.wrttlmod(ng).or.wrtrpmod(ng))) THEN
2141 CALL pio_netcdf_put_fvar (ng, model, obs(ng)%name, &
2142 & vname(1,idobsz), &
2143 & zobs(mstr:mend), &
2144 & (/nstrobs(ng)/), (/nobs(ng)/), &
2145 & piofile = obs(ng)%pioFile, &
2146 & piovar = obs(ng)%pioVar(idobsz)%vd)
2147 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2148
2149 IF (model.eq.iadm) THEN
2150 lwrote=obssurvey(ng).eq.1
2151 ELSE
2152 lwrote=obssurvey(ng).eq.nsurvey(ng)
2153 END IF
2154 IF (lwrote) wrote_zobs(ng)=.true.
2155 END IF
2156# endif
2157
2158# if defined BALANCE_OPERATOR && defined ZETA_ELLIPTIC
2159!
2160! Write out reference free-surface used in the balance operator.
2161!
2162 IF (wrtzetaref(ng).and.balance(isfsur)) THEN
2163 IF (dav(ng)%pioVar(idfsur)%dkind.eq.pio_double) THEN
2164 iodesc => iodesc_dp_r2dvar(ng)
2165 ELSE
2166 iodesc => iodesc_sp_r2dvar(ng)
2167 END IF
2168 tindex=0
2169 scale=1.0_dp
2170 status=nf_fwrite2d(ng, model, dav(ng)%pioFile, idfsur, &
2171 & dav(ng)%pioVar(idfsur), tindex, &
2172 & iodesc, &
2173 & lbi, ubi, lbj, ubj, scale, &
2174# ifdef MASKING
2175 & grid(ng) % rmask, &
2176# endif
2177 & fourdvar(ng) % zeta_ref, &
2178 & setfillval = .false.)
2179 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2180 wrtzetaref(ng)=.false.
2181 END IF
2182# endif
2183!
2184!-----------------------------------------------------------------------
2185! Synchronize observations NetCDF file to disk.
2186!-----------------------------------------------------------------------
2187!
2188 IF (wrtnlmod(ng).or.wrttlmod(ng).or.wrtrpmod(ng)) THEN
2189 CALL pio_netcdf_sync (ng, model, dav(ng)%name, &
2190 & dav(ng)%pioFile)
2191 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2192
2193# ifdef SOLVE3D
2194 IF (nrun.eq.1) THEN
2195 CALL pio_netcdf_sync (ng, model, obs(ng)%name, &
2196 & obs(ng)%pioFile)
2197 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2198 END IF
2199# endif
2200 END IF
2201!
2202 RETURN
mod_pio_netcdf::pio_netcdf_put_fvar
Definition mod_pio_netcdf.F:94
mod_pio_netcdf::pio_netcdf_put_ivar
Definition mod_pio_netcdf.F:108
mod_pio_netcdf
Definition mod_pio_netcdf.F:2
mod_pio_netcdf::pio_netcdf_sync
subroutine, public pio_netcdf_sync(ng, model, ncname, piofile)
Definition mod_pio_netcdf.F:9275
mod_pio_netcdf::iodesc_sp_r2dvar
type(io_desc_t), dimension(:), pointer iodesc_sp_r2dvar
Definition mod_pio_netcdf.F:194
mod_pio_netcdf::iodesc_dp_r2dvar
type(io_desc_t), dimension(:), pointer iodesc_dp_r2dvar
Definition mod_pio_netcdf.F:233

References mod_scalars::balance, mod_fourdvar::bgerr, mod_fourdvar::bgthresh, mod_iounits::dav, mod_scalars::exit_flag, strings_mod::founderror(), mod_fourdvar::fourdvar, mod_grid::grid, mod_fourdvar::havenlmod, mod_fourdvar::havetlmod, mod_param::iadm, mod_ncparam::idbger, mod_ncparam::idbgth, mod_ncparam::idfsur, mod_ncparam::idincr, mod_ncparam::idinno, mod_ncparam::idmomf, mod_ncparam::idmomi, mod_ncparam::idnlmf, mod_ncparam::idnlmi, mod_ncparam::idnlmo, mod_ncparam::idnlmp, mod_ncparam::idnlmu, mod_ncparam::idobss, mod_ncparam::idobsz, mod_ncparam::idresi, mod_ncparam::idtlmo, mod_scalars::inner, mod_fourdvar::innovation, mod_pio_netcdf::iodesc_dp_r2dvar, mod_pio_netcdf::iodesc_sp_r2dvar, mod_ncparam::isfsur, mod_fourdvar::misfit, mod_fourdvar::nlincrement, mod_fourdvar::nlmodval, mod_fourdvar::nobs, mod_scalars::noerror, mod_scalars::nouter, mod_scalars::nrun, mod_fourdvar::nstrobs, mod_fourdvar::nsurvey, mod_iounits::obs, mod_fourdvar::obserr, mod_fourdvar::obsscale, mod_fourdvar::obssurvey, mod_fourdvar::obsval, mod_scalars::outer, mod_pio_netcdf::pio_netcdf_sync(), mod_fourdvar::residual, mod_iounits::sourcefile, mod_fourdvar::unvetted, mod_ncparam::vname, mod_fourdvar::wrote_zobs, mod_fourdvar::wrtimpact_bc, mod_fourdvar::wrtimpact_fc, mod_fourdvar::wrtimpact_ic, mod_fourdvar::wrtimpact_tot, mod_fourdvar::wrtmisfit, mod_fourdvar::wrtnlmod, mod_fourdvar::wrtobsscale, mod_fourdvar::wrtrpmod, mod_fourdvar::wrttlmod, mod_fourdvar::wrtzetaref, and mod_fourdvar::zobs.

Referenced by obs_write().

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