obs_write.F

Go to the documentation of this file.

#include "cppdefs.h"
      MODULE obs_write_mod
#if (defined FOUR_DVAR || defined VERIFICATION) && defined OBSERVATIONS
!
!git $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2025 The ROMS Group                              !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  This routine interpolates nonlinear (background) and/or tangent     !
!  linear model (increments)  state at observations location  when     !
!  appropriate.                                                        !
!                                                                      !
!  It writes the model values at observation locations using the       !
!  standard NetCDF library or the Parallel-IO (PIO) library.           !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_fourdvar
      USE mod_grid
      USE mod_iounits
      USE mod_ncparam
      USE mod_ocean
      USE mod_scalars
      USE mod_stepping
!
# ifdef DISTRIBUTE
      USE distribute_mod,  ONLY : mp_collect
# endif
      USE extract_obs_mod, ONLY : extract_obs2d
# ifdef SOLVE3D
      USE extract_obs_mod, ONLY : extract_obs3d
# endif
# if defined BALANCE_OPERATOR && defined ZETA_ELLIPTIC
      USE nf_fwrite2d_mod, ONLY : nf_fwrite2d
# endif
      USE strings_mod,     ONLY : founderror
!
      implicit none
!
      PUBLIC  :: obs_write
      PRIVATE :: obs_operator
      PRIVATE :: obs_write_nf90
# if defined PIO_LIB && defined DISTRIBUTE
      PRIVATE :: obs_write_pio
# endif
!
      CONTAINS
!
!***********************************************************************
      SUBROUTINE obs_write (ng, tile, model)
!***********************************************************************
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile, model
!
!  Local variable declarations.
!
      integer :: mstr, mend
      integer :: lbi, ubi, lbj, ubj
      integer :: obssum, obsvoid
      integer :: i, ie, is, iobs, itrc
!
      character (len=50) :: string
      character (len=*), parameter :: myfile =                          &
     &  __FILE__
!
!-----------------------------------------------------------------------
!  Extract state vector at the observation locations in space and
!  time.
!-----------------------------------------------------------------------
!
      h_operator : IF (processobs(ng)) THEN
        lbi=bounds(ng)%LBi(tile)
        ubi=bounds(ng)%UBi(tile)
        lbj=bounds(ng)%LBj(tile)
        ubj=bounds(ng)%UBj(tile)
!
        CALL obs_operator (ng, tile, model,                             &
     &                     lbi, ubi, lbj, ubj,                          &
     &                     mstr, mend)
!
!-----------------------------------------------------------------------
!  Write out state vector to output Data Assimilation Variables (DAV)
!  NetCDF file.
!-----------------------------------------------------------------------
!
        SELECT CASE (dav(ng)%IOtype)
          CASE (io_nf90)
            CALL obs_write_nf90 (ng, tile, model,                       &
     &                           lbi, ubi, lbj, ubj,                    &
     &                           mstr, mend)
 
# if defined PIO_LIB && defined DISTRIBUTE
          CASE (io_pio)
            CALL obs_write_pio (ng, tile, model,                        &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          mstr, mend)
# endif
          CASE DEFAULT
            IF (master) WRITE (stdout,10) dav(ng)%IOtype
            exit_flag=3
        END SELECT
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
!
!-----------------------------------------------------------------------
!  Set counters for the number of rejected observations for each state
!  variable.
!-----------------------------------------------------------------------
!
        DO iobs=mstr,mend
          IF (obsscale(iobs).lt.1.0) THEN
            IF  (obstype(iobs).eq.obsstate2type(isfsur)) THEN
              fourdvar(ng)%ObsReject(isfsur)=                           &
     &                              fourdvar(ng)%ObsReject(isfsur)+1
            ELSE IF (obstype(iobs).eq.obsstate2type(isubar)) THEN
              fourdvar(ng)%ObsReject(isubar)=                           &
     &                              fourdvar(ng)%ObsReject(isubar)+1
            ELSE IF (obstype(iobs).eq.obsstate2type(isvbar)) THEN
              fourdvar(ng)%ObsReject(isvbar)=                           &
     &                              fourdvar(ng)%ObsReject(isvbar)+1
# ifdef SOLVE3D
            ELSE IF (obstype(iobs).eq.obsstate2type(isuvel)) THEN
              fourdvar(ng)%ObsReject(isuvel)=                           &
     &                              fourdvar(ng)%ObsReject(isuvel)+1
            ELSE IF (obstype(iobs).eq.obsstate2type(isvvel)) THEN
              fourdvar(ng)%ObsReject(isvvel)=                           &
     &                              fourdvar(ng)%ObsReject(isvvel)+1
            ELSE IF (obstype(iobs).eq.obsstate2type(isradial)) THEN
              fourdvar(ng)%ObsReject(isradial)=                         &
     &                              fourdvar(ng)%ObsReject(isradial)+1
            ELSE
              DO itrc=1,nt(ng)
                IF (obstype(iobs).eq.obsstate2type(istvar(itrc))) THEN
                  i=istvar(itrc)
                  fourdvar(ng)%ObsReject(i)=fourdvar(ng)%ObsReject(i)+1
                END IF
              END DO
# endif
            END IF
          END IF
        END DO
!
!  Load total available and rejected observations into structure
!  array.
!
        DO i=1,nobsvar(ng)
          fourdvar(ng)%ObsCount(0)=fourdvar(ng)%ObsCount(0)+            &
     &                             fourdvar(ng)%ObsCount(i)
          fourdvar(ng)%ObsReject(0)=fourdvar(ng)%ObsReject(0)+          &
     &                              fourdvar(ng)%ObsReject(i)
        END DO
!
!-----------------------------------------------------------------------
!  Report observation processing information.
!-----------------------------------------------------------------------
!
        IF (master) THEN
          obssum=0
          obsvoid=0
          is=nstrobs(ng)
          DO i=1,nobsvar(ng)
            IF (fourdvar(ng)%ObsCount(i).gt.0) THEN
              ie=is+fourdvar(ng)%ObsCount(i)-1
              WRITE (stdout,20) trim(obsname(i)), is, ie,               &
     &                          ie-is+1, fourdvar(ng)%ObsReject(i)
              is=ie+1
              obssum=obssum+fourdvar(ng)%ObsCount(i)
              obsvoid=obsvoid+fourdvar(ng)%ObsReject(i)
            END IF
          END DO
          WRITE (stdout,30) obssum, obsvoid,                            &
     &                      fourdvar(ng)%ObsCount(0),                   &
     &                      fourdvar(ng)%ObsReject(0)
        END IF
!
        IF (wrtnlmod(ng)) THEN
          string='Wrote NLM state at observation locations,         '
# ifdef TLM_OBS
#  ifdef WEAK_CONSTRAINT
        ELSE IF (wrttlmod(ng)) THEN
          string='Wrote TLM state at observation locations,         '
        ELSE IF (wrtrpmod(ng)) THEN
          string='Wrote RPM state at observation locations,         '
#  else
#   ifdef I4DVAR_ANA_SENSITIVITY
        ELSE IF (wrttlmod(ng)) THEN
          string='Wrote 4DVAR observation sensitivity,              '
#   else
        ELSE IF (wrttlmod(ng)) THEN
          string='Wrote TLM increments at observation locations,    '
#   endif
#  endif
# endif
        END IF
        IF (master) THEN
          IF (wrtnlmod(ng).or.wrttlmod(ng).or.wrtrpmod(ng)) THEN
             WRITE (stdout,40) trim(string), nstrobs(ng), nendobs(ng)
          END IF
        END IF
 
      END IF h_operator
!
  10  FORMAT (' OBS_WRITE - Illegal output file type, io_type = ',i0,   &
     &        /,13x,'Check KeyWord ''OUT_LIB'' in ''roms.in''.')
  20  FORMAT (10x,a,t25,4(1x,i10))
  30  FORMAT (/,10x,'Total',t47,2(1x,i10),                              &
     &        /,10x,'Obs Tally',t47,2(1x,i10),/)
  40  FORMAT (1x,a,' datum = ',i0,' - ',i0,/)
!
      RETURN
      END SUBROUTINE obs_write
!
!***********************************************************************
      SUBROUTINE obs_operator (ng, tile, model,                         &
     &                         LBi, UBi, LBj, UBj,                      &
     &                         Mstr, Mend)
!***********************************************************************
!
!  Imported variable declarations.
!
      integer, intent(in)  :: ng, tile, model
      integer, intent(in)  :: lbi, ubi, lbj, ubj
 
      integer, intent(out) :: mstr, mend
!
!  Local variable declarations.
!
      integer :: i, ic, iobs, itrc
# ifdef SOLVE3D
      integer :: j, k
# endif
# ifdef DISTRIBUTE
      integer :: ncollect
# endif
!
      real(r8), parameter :: inival = 0.0_r8
 
      real(r8) :: angle
# ifdef BGQC
      real(r8) :: df1, df2, thresh
# endif
      real(r8) :: misfit(mobs), unvetted(mobs)
      real(r8) :: uradial(mobs), vradial(mobs)
# ifndef VERIFICATION
      real(r8) :: ubgerr(mobs), vbgerr(mobs)
# endif
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", obs_operator"
 
# include "set_bounds.h"
!
      sourcefile=myfile
!
!=======================================================================
!  Interpolate model state at observation locations.
!=======================================================================
!
# ifdef WEAK_CONSTRAINT
!
!  Set starting and ending indices of observations to process for the
!  current time survey. In the dual formulation, the entire observation
!  vector for the assimilation window is maintained and the data array
!  indexes are global.
!
      mstr=nstrobs(ng)
      mend=nendobs(ng)
# else
!
!  Set starting and ending indices of observations to process for the
!  current time survey. In the primal formulation, only the observations
!  for the current assimilation window are maintained and the data array
!  indexes are local (1:Nobs; starting index is always one).
!
      mstr=1
      mend=nobs(ng)
# endif
!
!  Initialize observation reject/accept processing flag.
!
      DO iobs=mstr,mend
        unvetted(iobs)=inival
        obsvetting(iobs)=inival
      END DO
 
# ifndef I4DVAR_ANA_SENSITIVITY
!
!  Some entries are not computed in the extraction routine.  Set values
!  to zero to avoid problems when writing non initialized values.
#  ifdef DISTRIBUTE
!  Notice that only the appropriate indices are zero-out to facilate
!  collecting all the extrated data as sum between all nodes.
#  endif
!
      IF (wrtnlmod(ng)) THEN
        DO iobs=mstr,mend
          nlmodval(iobs)=inival
#  ifndef VERIFICATION
          bgerr(iobs)=inival
          ubgerr(iobs)=inival
          vbgerr(iobs)=inival
#  endif
        END DO
 
#  ifdef BGQC
!
!  Set background quality control check (BgThresh). The background
!  quality control is either in terms of the state variable index
!  (1:MstateVar) or observation provenance index.
!
        IF (bgqc_type(ng).eq.1) THEN          ! state variable index QC
          DO iobs=mstr,mend
            DO i=1,mstatevar
              IF (obstype(iobs).eq.obsstate2type(i)) THEN
                bgthresh(iobs)=s_bgqc(i,ng)
                EXIT
              END IF
            END DO
          END DO
        ELSE IF (bgqc_type(ng).eq.2) THEN     ! provenance index QC
          DO iobs=mstr,mend
            bgthresh(iobs)=bgqc_large         ! do not reject
            DO i=1,nprovenance(ng)
              IF (obsprov(iobs).eq.iprovenance(i,ng)) THEN
                bgthresh(iobs)=p_bgqc(i,ng)
                EXIT
              END IF
            END DO
          END DO
        END IF
#  endif
      END IF
 
#  if defined TLM_OBS && !defined SP4DVAR
      IF (wrttlmod(ng).or.wrtrpmod(ng)) THEN
        DO iobs=mstr,mend
          tlmodval(iobs)=inival
        END DO
      END IF
#  endif
# endif
!
!-----------------------------------------------------------------------
!  Free-surface operator.
!-----------------------------------------------------------------------
!
# ifndef I4DVAR_ANA_SENSITIVITY
!  Nonlinear free-surface at observation locations.
!
      IF (wrtnlmod(ng).and.                                             &
     &    (fourdvar(ng)%ObsCount(isfsur).gt.0)) THEN
        CALL extract_obs2d (ng, 0, lm(ng)+1, 0, mm(ng)+1,               &
     &                      lbi, ubi, lbj, ubj,                         &
     &                      obsstate2type(isfsur),                      &
     &                      mobs, mstr, mend,                           &
     &                      rxmin(ng), rxmax(ng),                       &
     &                      rymin(ng), rymax(ng),                       &
     &                      time(ng), dt(ng),                           &
     &                      obstype, obsvetting,                        &
     &                      tobs, xobs, yobs,                           &
     &                      ocean(ng)%zeta(:,:,kout),                   &
#  ifdef MASKING
     &                      grid(ng)%rmask,                             &
#  endif
     &                      nlmodval)
#  ifndef VERIFICATION
!
!  Free-surface backgound error standard deviation at observation
!  locations.
!
        CALL extract_obs2d (ng, 0, lm(ng)+1, 0, mm(ng)+1,               &
     &                      lbi, ubi, lbj, ubj,                         &
     &                      obsstate2type(isfsur),                      &
     &                      mobs, mstr, mend,                           &
     &                      rxmin(ng), rxmax(ng),                       &
     &                      rymin(ng), rymax(ng),                       &
     &                      time(ng), dt(ng),                           &
     &                      obstype, obsvetting,                        &
     &                      tobs, xobs, yobs,                           &
     &                      ocean(ng)%e_zeta(:,:,1),                    &
#   ifdef MASKING
     &                      grid(ng)%rmask,                             &
#   endif
     &                      bgerr)
#  endif
      END IF
# endif
 
# ifdef TLM_OBS
!
!  Tangent linear free-surface at observation locations.
!
      IF ((wrttlmod(ng).or.(wrtrpmod(ng))).and.                         &
     &    (fourdvar(ng)%ObsCount(isfsur).gt.0)) THEN
        CALL extract_obs2d (ng, 0, lm(ng)+1, 0, mm(ng)+1,               &
     &                      lbi, ubi, lbj, ubj,                         &
     &                      obsstate2type(isfsur),                      &
     &                      mobs, mstr, mend,                           &
     &                      rxmin(ng), rxmax(ng),                       &
     &                      rymin(ng), rymax(ng),                       &
     &                      time(ng), dt(ng),                           &
     &                      obstype, obsvetting,                        &
     &                      tobs, xobs, yobs,                           &
     &                      ocean(ng)%tl_zeta(:,:,kout),                &
#  ifdef MASKING
     &                      grid(ng)%rmask,                             &
#  endif
     &                      tlmodval)
      END IF
# endif
!
!-----------------------------------------------------------------------
!  Vertically integrated u-velocity operator.
!-----------------------------------------------------------------------
!
# ifndef I4DVAR_ANA_SENSITIVITY
!  Nonlinear 2D u-velocity at observation locations.
!
      IF (wrtnlmod(ng).and.                                             &
     &    (fourdvar(ng)%ObsCount(isubar).gt.0)) THEN
        CALL extract_obs2d (ng, 1, lm(ng)+1, 0, mm(ng)+1,               &
     &                      lbi, ubi, lbj, ubj,                         &
     &                      obsstate2type(isubar),                      &
     &                      mobs, mstr, mend,                           &
     &                      uxmin(ng), uxmax(ng),                       &
     &                      uymin(ng), uymax(ng),                       &
     &                      time(ng), dt(ng),                           &
     &                      obstype, obsvetting,                        &
     &                      tobs, xobs, yobs,                           &
     &                      ocean(ng)%ubar(:,:,kout),                   &
#  ifdef MASKING
     &                      grid(ng)%umask,                             &
#  endif
     &                      nlmodval)
 
#  ifndef VERIFICATION
!
!  2D u-velocity backgound error standard deviation at observation
!  locations.
!
        CALL extract_obs2d (ng, 1, lm(ng)+1, 0, mm(ng)+1,               &
     &                      lbi, ubi, lbj, ubj,                         &
     &                      obsstate2type(isubar),                      &
     &                      mobs, mstr, mend,                           &
     &                      uxmin(ng), uxmax(ng),                       &
     &                      uymin(ng), uymax(ng),                       &
     &                      time(ng), dt(ng),                           &
     &                      obstype, obsvetting,                        &
     &                      tobs, xobs, yobs,                           &
     &                      ocean(ng)%e_ubar(:,:,1),                    &
#   ifdef MASKING
     &                      grid(ng)%umask,                             &
#   endif
     &                      bgerr)
#  endif
      END IF
# endif
 
# ifdef TLM_OBS
!
!  Tangent linear 2D u-velocity at observation locations.
!
      IF ((wrttlmod(ng).or.(wrtrpmod(ng))).and.                         &
     &    (fourdvar(ng)%ObsCount(isubar).gt.0)) THEN
        CALL extract_obs2d (ng, 1, lm(ng)+1, 0, mm(ng)+1,               &
     &                      lbi, ubi, lbj, ubj,                         &
     &                      obsstate2type(isubar),                      &
     &                      mobs, mstr, mend,                           &
     &                      uxmin(ng), uxmax(ng),                       &
     &                      uymin(ng), uymax(ng),                       &
     &                      time(ng), dt(ng),                           &
     &                      obstype, obsvetting,                        &
     &                      tobs, xobs, yobs,                           &
     &                      ocean(ng)%tl_ubar(:,:,kout),                &
#  ifdef MASKING
     &                      grid(ng)%umask,                             &
#  endif
     &                      tlmodval)
        END IF
# endif
!
!-----------------------------------------------------------------------
!  Vertically integrated v-velocity observations.
!-----------------------------------------------------------------------
!
# ifndef I4DVAR_ANA_SENSITIVITY
!  Nonlinear 2D v-velocity at observation locations.
!
      IF (wrtnlmod(ng).and.                                             &
     &    (fourdvar(ng)%ObsCount(isvbar).gt.0)) THEN
        CALL extract_obs2d (ng, 0, lm(ng)+1, 1, mm(ng)+1,               &
     &                      lbi, ubi, lbj, ubj,                         &
     &                      obsstate2type(isvbar),                      &
     &                      mobs, mstr, mend,                           &
     &                      vxmin(ng), vxmax(ng),                       &
     &                      vymin(ng), vymax(ng),                       &
     &                      time(ng), dt(ng),                           &
     &                      obstype, obsvetting,                        &
     &                      tobs, xobs, yobs,                           &
     &                      ocean(ng)%vbar(:,:,kout),                   &
#  ifdef MASKING
     &                      grid(ng)%vmask,                             &
#  endif
     &                      nlmodval)
 
#  ifndef VERIFICATION
!
!  2D v-velocity backgound error standard deviation at observation
!  locations.
!
        CALL extract_obs2d (ng, 0, lm(ng)+1, 1, mm(ng)+1,               &
     &                      lbi, ubi, lbj, ubj,                         &
     &                      obsstate2type(isvbar),                      &
     &                      mobs, mstr, mend,                           &
     &                      vxmin(ng), vxmax(ng),                       &
     &                      vymin(ng), vymax(ng),                       &
     &                      time(ng), dt(ng),                           &
     &                      obstype, obsvetting,                        &
     &                      tobs, xobs, yobs,                           &
     &                      ocean(ng)%e_vbar(:,:,1),                    &
#   ifdef MASKING
     &                      grid(ng)%vmask,                             &
#   endif
     &                      bgerr)
#  endif
      END IF
# endif
 
# ifdef TLM_OBS
!
!  Tangent linear 2D v-velocity at observation locations.
!
      IF ((wrttlmod(ng).or.(wrtrpmod(ng))).and.                         &
     &    (fourdvar(ng)%ObsCount(isvbar).gt.0)) THEN
        CALL extract_obs2d (ng, 0, lm(ng)+1, 1, mm(ng)+1,               &
     &                      lbi, ubi, lbj, ubj,                         &
     &                      obsstate2type(isvbar),                      &
     &                      mobs, mstr, mend,                           &
     &                      vxmin(ng), vxmax(ng),                       &
     &                      vymin(ng), vymax(ng),                       &
     &                      time(ng), dt(ng),                           &
     &                      obstype, obsvetting,                        &
     &                      tobs, xobs, yobs,                           &
     &                      ocean(ng)%tl_vbar(:,:,kout),                &
#  ifdef MASKING
     &                      grid(ng)%vmask,                             &
#  endif
     &                      tlmodval)
      END IF
# endif
 
# ifdef SOLVE3D
!
!-----------------------------------------------------------------------
!  3D u-velocity observations.
!-----------------------------------------------------------------------
!
#  ifndef I4DVAR_ANA_SENSITIVITY
!  Nonlinear 3D u-velocity at observation locations.
!
      IF (wrtnlmod(ng).and.                                             &
     &    (fourdvar(ng)%ObsCount(isuvel).gt.0)) THEN
        DO k=1,n(ng)
          DO j=jstr-1,jend+1
            DO i=istru-1,iend+1
              grid(ng)%z_v(i,j,k)=0.5_r8*(grid(ng)%z_r(i-1,j,k)+        &
     &                                    grid(ng)%z_r(i  ,j,k))
            END DO
          END DO
        END DO
        CALL extract_obs3d (ng, 1, lm(ng)+1, 0, mm(ng)+1,               &
     &                      lbi, ubi, lbj, ubj, 1, n(ng),               &
     &                      obsstate2type(isuvel),                      &
     &                      mobs, mstr, mend,                           &
     &                      uxmin(ng), uxmax(ng),                       &
     &                      uymin(ng), uymax(ng),                       &
     &                      time(ng), dt(ng),                           &
     &                      obstype, obsvetting,                        &
     &                      tobs, xobs, yobs, zobs,                     &
     &                      ocean(ng)%u(:,:,:,nout),                    &
     &                      grid(ng)%z_v,                               &
#   ifdef MASKING
     &                      grid(ng)%umask,                             &
#   endif
     &                      nlmodval)
 
#   ifndef VERIFICATION
!
!  3D u-velocity backgound error standard deviation at observation
!  locations.
!
        CALL extract_obs3d (ng, 1, lm(ng)+1, 0, mm(ng)+1,               &
     &                      lbi, ubi, lbj, ubj, 1, n(ng),               &
     &                      obsstate2type(isuvel),                      &
     &                      mobs, mstr, mend,                           &
     &                      uxmin(ng), uxmax(ng),                       &
     &                      uymin(ng), uymax(ng),                       &
     &                      time(ng), dt(ng),                           &
     &                      obstype,  obsvetting,                       &
     &                      tobs, xobs, yobs, zobs,                     &
     &                      ocean(ng)%e_u(:,:,:,1),                     &
     &                      grid(ng)%z_v,                               &
#    ifdef MASKING
     &                      grid(ng)%umask,                             &
#    endif
     &                      bgerr)
#   endif
      END IF
#  endif
 
#  ifdef TLM_OBS
!
!  Tangent linear 3D u-velocity at observation locations.
!
      IF ((wrttlmod(ng).or.(wrtrpmod(ng))).and.                         &
     &    (fourdvar(ng)%ObsCount(isuvel).gt.0)) THEN
        DO k=1,n(ng)
          DO j=jstr-1,jend+1
            DO i=istru-1,iend+1
              grid(ng)%z_v(i,j,k)=0.5_r8*(grid(ng)%z_r(i-1,j,k)+        &
     &                                    grid(ng)%z_r(i  ,j,k))
            END DO
          END DO
        END DO
        CALL extract_obs3d (ng, 1, lm(ng)+1, 0, mm(ng)+1,               &
     &                      lbi, ubi, lbj, ubj, 1, n(ng),               &
     &                      obsstate2type(isuvel),                      &
     &                      mobs, mstr, mend,                           &
     &                      uxmin(ng), uxmax(ng),                       &
     &                      uymin(ng), uymax(ng),                       &
     &                      time(ng), dt(ng),                           &
     &                      obstype, obsvetting,                        &
     &                      tobs, xobs, yobs, zobs,                     &
     &                      ocean(ng)%tl_u(:,:,:,nout),                 &
     &                      grid(ng)%z_v,                               &
#   ifdef MASKING
     &                      grid(ng)%umask,                             &
#   endif
     &                      tlmodval)
      END IF
#  endif
!
!-----------------------------------------------------------------------
!  3D v-velocity observations.
!-----------------------------------------------------------------------
!
#  ifndef I4DVAR_ANA_SENSITIVITY
!  Nonlinear 3D v-velocity at observation locations.
!
      IF (wrtnlmod(ng).and.                                             &
     &    (fourdvar(ng)%ObsCount(isvvel).gt.0)) THEN
        DO k=1,n(ng)
          DO j=jstrv-1,jend+1
            DO i=istr-1,iend+1
              grid(ng)%z_v(i,j,k)=0.5_r8*(grid(ng)%z_r(i,j-1,k)+        &
     &                                    grid(ng)%z_r(i,j  ,k))
            END DO
          END DO
        END DO
        CALL extract_obs3d (ng, 0, lm(ng)+1, 1, mm(ng)+1,               &
     &                      lbi, ubi, lbj, ubj, 1, n(ng),               &
     &                      obsstate2type(isvvel),                      &
     &                      mobs, mstr, mend,                           &
     &                      vxmin(ng), vxmax(ng),                       &
     &                      vymin(ng), vymax(ng),                       &
     &                      time(ng), dt(ng),                           &
     &                      obstype, obsvetting,                        &
     &                      tobs, xobs, yobs, zobs,                     &
     &                      ocean(ng)%v(:,:,:,nout),                    &
     &                      grid(ng)%z_v,                               &
#   ifdef MASKING
     &                      grid(ng)%vmask,                             &
#   endif
     &                      nlmodval)
 
#   ifndef VERIFICATION
!
!  3D v-velocity backgound error standard deviation at observation
!  locations.
!
        CALL extract_obs3d (ng, 0, lm(ng)+1, 1, mm(ng)+1,               &
     &                      lbi, ubi, lbj, ubj, 1, n(ng),               &
     &                      obsstate2type(isvvel),                      &
     &                      mobs, mstr, mend,                           &
     &                      vxmin(ng), vxmax(ng),                       &
     &                      vymin(ng), vymax(ng),                       &
     &                      time(ng), dt(ng),                           &
     &                      obstype, obsvetting,                        &
     &                      tobs, xobs, yobs, zobs,                     &
     &                      ocean(ng)%e_v(:,:,:,1),                     &
     &                      grid(ng)%z_v,                               &
#    ifdef MASKING
     &                      grid(ng)%vmask,                             &
#    endif
     &                      bgerr)
#   endif
      END IF
#  endif
 
#  ifdef TLM_OBS
!
!  Tangent linear 3D v-velocity at observation locations.
!
      IF ((wrttlmod(ng).or.(wrtrpmod(ng))).and.                         &
     &    (fourdvar(ng)%ObsCount(isvvel).gt.0)) THEN
        DO k=1,n(ng)
          DO j=jstrv-1,jend+1
            DO i=istr-1,iend+1
              grid(ng)%z_v(i,j,k)=0.5_r8*(grid(ng)%z_r(i,j-1,k)+        &
     &                                    grid(ng)%z_r(i,j  ,k))
            END DO
          END DO
        END DO
        CALL extract_obs3d (ng, 0, lm(ng)+1, 1, mm(ng)+1,               &
     &                      lbi, ubi, lbj, ubj, 1, n(ng),               &
     &                      obsstate2type(isvvel),                      &
     &                      mobs, mstr, mend,                           &
     &                      vxmin(ng), vxmax(ng),                       &
     &                      vymin(ng), vymax(ng),                       &
     &                      time(ng), dt(ng),                           &
     &                      obstype, obsvetting,                        &
     &                      tobs, xobs, yobs, zobs,                     &
     &                      ocean(ng)%tl_v(:,:,:,nout),                 &
     &                      grid(ng)%z_v,                               &
#   ifdef MASKING
     &                      grid(ng)%vmask,                             &
#   endif
     &                      tlmodval)
      END IF
#  endif
!
!-----------------------------------------------------------------------
!  Radial Velocity. The observations are in terms of radial speed and
!  angle (stored in obs_meta).  The observation angle converts the
!  velocity components to geographical EAST and North components.
#  ifdef RADIAL_ANGLE_CCW_EAST
!  The radial velocity observations are processed as magnitude and
!  heading angle (obs_meta; radians) in the math convention: an
!  azimuth that is counterclockwise from TRUE East.
!
!  In curvilinear coordinates, the radial forward problem is:
!
!     radial = u * COS(obs_meta - angler) + v * SIN(obs_meta - angler)
#  else
!  By default, the radial velocity observations are processed as
!  magnitude and heading angle (obs_meta; radians) in the navigation
!  convention: an azimuth that is clockwise from TRUE North.
!
!  In curvilinear coordinates, the radial forward problem is:
!
!     radial = u * SIN(obs_meta + angler) + v * COS(obs_meta + angler)
#  endif
!-----------------------------------------------------------------------
!
#  ifndef I4DVAR_ANA_SENSITIVITY
      IF (wrtnlmod(ng).and.                                             &
     &    (fourdvar(ng)%ObsCount(isradial).gt.0)) THEN
        DO iobs=mstr,mend
          uradial(iobs)=inival
          vradial(iobs)=inival
        END DO
 
#   ifdef CURVGRID
!
!  If curvilinear coordinates, interpolate grid rotation angle (radians)
!  for each observation location.
!
        CALL extract_obs2d (ng, 0, lm(ng)+1, 0, mm(ng)+1,               &
     &                      lbi, ubi, lbj, ubj,                         &
     &                      obsstate2type(isradial),                    &
     &                      mobs, mstr, mend,                           &
     &                      rxmin(ng), rxmax(ng),                       &
     &                      rymin(ng), rymax(ng),                       &
     &                      time(ng), dt(ng),                           &
     &                      obstype, obsvetting,                        &
     &                      tobs, xobs, yobs,                           &
     &                      grid(ng)%angler,                            &
#    ifdef MASKING
     &                      grid(ng)%rmask,                             &
#    endif
     &                      obsangler)
#   endif
!
!  Radial nonlinear u-component at observation locations.
!
        DO k=1,n(ng)
          DO j=jstr-1,jend+1
            DO i=istru-1,iend+1
              grid(ng)%z_v(i,j,k)=0.5_r8*(grid(ng)%z_r(i-1,j,k)+        &
     &                                    grid(ng)%z_r(i  ,j,k))
            END DO
          END DO
        END DO
        CALL extract_obs3d (ng, 1, lm(ng)+1, 0, mm(ng)+1,               &
     &                      lbi, ubi, lbj, ubj, 1, n(ng),               &
     &                      obsstate2type(isradial),                    &
     &                      mobs, mstr, mend,                           &
     &                      uxmin(ng), uxmax(ng),                       &
     &                      uymin(ng), uymax(ng),                       &
     &                      time(ng), dt(ng),                           &
     &                      obstype,  obsvetting,                       &
     &                      tobs, xobs, yobs, zobs,                     &
     &                      ocean(ng)%u(:,:,:,nout),                    &
     &                      grid(ng)%z_v,                               &
#   ifdef MASKING
     &                      grid(ng)%umask,                             &
#   endif
     &                      uradial)
 
#   ifndef VERIFICATION
!
!  Radial u-component backgound error standard deviation at observation
!  locations.
!
        CALL extract_obs3d (ng, 1, lm(ng)+1, 0, mm(ng)+1,               &
     &                      lbi, ubi, lbj, ubj, 1, n(ng),               &
     &                      obsstate2type(isradial),                    &
     &                      mobs, mstr, mend,                           &
     &                      uxmin(ng), uxmax(ng),                       &
     &                      uymin(ng), uymax(ng),                       &
     &                      time(ng), dt(ng),                           &
     &                      obstype, obsvetting,                        &
     &                      tobs, xobs, yobs, zobs,                     &
     &                      ocean(ng)%e_u(:,:,:,1),                     &
     &                      grid(ng)%z_v,                               &
#    ifdef MASKING
     &                      grid(ng)%umask,                             &
#    endif
     &                      ubgerr)
#   endif
!
!  Radial nonlinear v-component at observation locations.
!
        DO k=1,n(ng)
          DO j=jstrv-1,jend+1
            DO i=istr-1,iend+1
              grid(ng)%z_v(i,j,k)=0.5_r8*(grid(ng)%z_r(i,j-1,k)+        &
     &                                    grid(ng)%z_r(i,j  ,k))
            END DO
          END DO
        END DO
        CALL extract_obs3d (ng, 0, lm(ng)+1, 1, mm(ng)+1,               &
     &                      lbi, ubi, lbj, ubj, 1, n(ng),               &
     &                      obsstate2type(isradial),                    &
     &                      mobs, mstr, mend,                           &
     &                      vxmin(ng), vxmax(ng),                       &
     &                      vymin(ng), vymax(ng),                       &
     &                      time(ng), dt(ng),                           &
     &                      obstype, obsvetting,                        &
     &                      tobs, xobs, yobs, zobs,                     &
     &                      ocean(ng)%v(:,:,:,nout),                    &
     &                      grid(ng)%z_v,                               &
#   ifdef MASKING
     &                      grid(ng)%vmask,                             &
#   endif
     &                      vradial)
 
#   ifndef VERIFICATION
!
!  Radial v-component backgound error standard deviation at observation
!  locations.
!
        CALL extract_obs3d (ng, 0, lm(ng)+1, 1, mm(ng)+1,               &
     &                      lbi, ubi, lbj, ubj, 1, n(ng),               &
     &                      obsstate2type(isradial),                    &
     &                      mobs, mstr, mend,                           &
     &                      vxmin(ng), vxmax(ng),                       &
     &                      vymin(ng), vymax(ng),                       &
     &                      time(ng), dt(ng),                           &
     &                      obstype, obsvetting,                        &
     &                      tobs, xobs, yobs, zobs,                     &
     &                      ocean(ng)%e_v(:,:,:,1),                     &
     &                      grid(ng)%z_v,                               &
#    ifdef MASKING
     &                      grid(ng)%vmask,                             &
#    endif
     &                      vbgerr)
#   endif
!
!  Compute nonlinear radial velocity.
!
        DO iobs=mstr,mend
          IF (obstype(iobs).eq.obsstate2type(isradial)) THEN
#   ifdef RADIAL_ANGLE_CCW_EAST
#    ifdef CURVGRID
            angle=obsmeta(iobs)-obsangler(iobs)
            nlmodval(iobs)=uradial(iobs)*cos(angle)+                    &
     &                     vradial(iobs)*sin(angle)
#    else
            nlmodval(iobs)=uradial(iobs)*cos(obsmeta(iobs))+            &
     &                     vradial(iobs)*sin(obsmeta(iobs))
#    endif
#   else
#    ifdef CURVGRID
            angle=obsmeta(iobs)+obsangler(iobs)
            nlmodval(iobs)=uradial(iobs)*sin(angle)+                    &
     &                     vradial(iobs)*cos(angle)
#    else
            nlmodval(iobs)=uradial(iobs)*sin(obsmeta(iobs))+            &
     &                     vradial(iobs)*cos(obsmeta(iobs))
#    endif
#   endif
#   ifndef VERIFICATION
            bgerr(iobs)=max(ubgerr(iobs), vbgerr(iobs))
#   endif
          END IF
        END DO
      END IF
#  endif
 
#  ifdef TLM_OBS
!
!  Tangent linear radial velocities.
!
      IF ((wrttlmod(ng).or.(wrtrpmod(ng))).and.                         &
     &    (fourdvar(ng)%ObsCount(isradial).gt.0)) THEN
        DO iobs=mstr,mend
          uradial(iobs)=inival
          vradial(iobs)=inival
        END DO
 
#   ifdef CURVGRID
!
!  If curvilinear coordinates, interpolate grid rotation angle (radians)
!  for each observation location.
!
        CALL extract_obs2d (ng, 0, lm(ng)+1, 0, mm(ng)+1,               &
     &                      lbi, ubi, lbj, ubj,                         &
     &                      obsstate2type(isradial),                    &
     &                      mobs, mstr, mend,                           &
     &                      rxmin(ng), rxmax(ng),                       &
     &                      rymin(ng), rymax(ng),                       &
     &                      time(ng), dt(ng),                           &
     &                      obstype, obsvetting,                        &
     &                      tobs, xobs, yobs,                           &
     &                      grid(ng)%angler,                            &
#    ifdef MASKING
     &                      grid(ng)%rmask,                             &
#    endif
     &                      obsangler)
#   endif
!
!  Tangent linear radial u-component at the observation locations.
!
        DO k=1,n(ng)
          DO j=jstr-1,jend+1
            DO i=istru-1,iend+1
              grid(ng)%z_v(i,j,k)=0.5_r8*(grid(ng)%z_r(i-1,j,k)+        &
     &                                    grid(ng)%z_r(i  ,j,k))
            END DO
          END DO
        END DO
        CALL extract_obs3d (ng, 1, lm(ng)+1, 0, mm(ng)+1,               &
     &                      lbi, ubi, lbj, ubj, 1, n(ng),               &
     &                      obsstate2type(isradial),                    &
     &                      mobs, mstr, mend,                           &
     &                      uxmin(ng), uxmax(ng),                       &
     &                      uymin(ng), uymax(ng),                       &
     &                      time(ng), dt(ng),                           &
     &                      obstype,  obsvetting,                       &
     &                      tobs, xobs, yobs, zobs,                     &
     &                      ocean(ng)%tl_u(:,:,:,nout),                 &
     &                      grid(ng)%z_v,                               &
#   ifdef MASKING
     &                      grid(ng)%umask,                             &
#   endif
     &                      uradial)
!
!  Tangent linear radial v-component at the observation locations.
!
        DO k=1,n(ng)
          DO j=jstrv-1,jend+1
            DO i=istr-1,iend+1
              grid(ng)%z_v(i,j,k)=0.5_r8*(grid(ng)%z_r(i,j-1,k)+        &
     &                                    grid(ng)%z_r(i,j  ,k))
            END DO
          END DO
        END DO
        CALL extract_obs3d (ng, 0, lm(ng)+1, 1, mm(ng)+1,               &
     &                      lbi, ubi, lbj, ubj, 1, n(ng),               &
     &                      obsstate2type(isradial),                    &
     &                      mobs, mstr, mend,                           &
     &                      vxmin(ng), vxmax(ng),                       &
     &                      vymin(ng), vymax(ng),                       &
     &                      time(ng), dt(ng),                           &
     &                      obstype, obsvetting,                        &
     &                      tobs, xobs, yobs, zobs,                     &
     &                      ocean(ng)%tl_v(:,:,:,nout),                 &
     &                      grid(ng)%z_v,                               &
#   ifdef MASKING
     &                      grid(ng)%vmask,                             &
#   endif
     &                      vradial)
!
!  Compute tangent linear radial velocity.
!
        DO iobs=mstr,mend
          IF (obstype(iobs).eq.obsstate2type(isradial)) THEN
#   ifdef RADIAL_ANGLE_CCW_EAST
#    ifdef CURVGRID
            angle=obsmeta(iobs)-obsangler(iobs)
            tlmodval(iobs)=uradial(iobs)*cos(angle)+                    &
     &                     vradial(iobs)*sin(angle)
#    else
            tlmodval(iobs)=uradial(iobs)*cos(obsmeta(iobs))+            &
     &                     vradial(iobs)*sin(obsmeta(iobs))
#    endif
#   else
#    ifdef CURVGRID
            angle=obsmeta(iobs)+obsangler(iobs)
            tlmodval(iobs)=uradial(iobs)*sin(angle)+                    &
     &                     vradial(iobs)*cos(angle)
#    else
            tlmodval(iobs)=uradial(iobs)*sin(obsmeta(iobs))+            &
     &                     vradial(iobs)*cos(obsmeta(iobs))
#    endif
#   endif
          END IF
        END DO
      END IF
#  endif
!
!-----------------------------------------------------------------------
!  Tracers variable observations.
!-----------------------------------------------------------------------
!
      DO itrc=1,nt(ng)
 
#  ifndef I4DVAR_ANA_SENSITIVITY
!
!  Nonlinear tracer at the observation locations.
 
        IF (wrtnlmod(ng).and.                                           &
     &      (fourdvar(ng)%ObsCount(istvar(itrc)).gt.0)) THEN
          CALL extract_obs3d (ng, 0, lm(ng)+1, 0, mm(ng)+1,             &
     &                        lbi, ubi, lbj, ubj, 1, n(ng),             &
     &                        obsstate2type(istvar(itrc)),              &
     &                        mobs, mstr, mend,                         &
     &                        rxmin(ng), rxmax(ng),                     &
     &                        rymin(ng), rymax(ng),                     &
     &                        time(ng), dt(ng),                         &
     &                        obstype, obsvetting,                      &
     &                        tobs, xobs, yobs, zobs,                   &
     &                        ocean(ng)%t(:,:,:,nout,itrc),             &
     &                        grid(ng)%z_r,                             &
#   ifdef MASKING
     &                        grid(ng)%rmask,                           &
#   endif
     &                        nlmodval)
 
#   ifndef VERIFICATION
!
!  Tracer backgound error standard deviation at observation locations.
!
          CALL extract_obs3d (ng, 0, lm(ng)+1, 0, mm(ng)+1,             &
     &                        lbi, ubi, lbj, ubj, 1, n(ng),             &
     &                        obsstate2type(istvar(itrc)),              &
     &                        mobs, mstr, mend,                         &
     &                        rxmin(ng), rxmax(ng),                     &
     &                        rymin(ng), rymax(ng),                     &
     &                        time(ng), dt(ng),                         &
     &                        obstype, obsvetting,                      &
     &                        tobs, xobs, yobs, zobs,                   &
     &                        ocean(ng)%e_t(:,:,:,1,itrc),              &
     &                        grid(ng)%z_r,                             &
#    ifdef MASKING
     &                        grid(ng)%rmask,                           &
#    endif
     &                        bgerr)
#   endif
        END IF
#  endif
 
#  ifdef TLM_OBS
!
!  Tangent linear tracer at the observation locations.
!
        IF ((wrttlmod(ng).or.(wrtrpmod(ng))).and.                       &
     &      (fourdvar(ng)%ObsCount(istvar(itrc)).gt.0)) THEN
          CALL extract_obs3d (ng, 0, lm(ng)+1, 0, mm(ng)+1,             &
     &                        lbi, ubi, lbj, ubj, 1, n(ng),             &
     &                        obsstate2type(istvar(itrc)),              &
     &                        mobs, mstr, mend,                         &
     &                        rxmin(ng), rxmax(ng),                     &
     &                        rymin(ng), rymax(ng),                     &
     &                        time(ng), dt(ng),                         &
     &                        obstype, obsvetting,                      &
     &                        tobs, xobs, yobs, zobs,                   &
     &                        ocean(ng)%tl_t(:,:,:,nout,itrc),          &
     &                        grid(ng)%z_r,                             &
#   ifdef MASKING
     &                        grid(ng)%rmask,                           &
#   endif
     &                        tlmodval)
        END IF
#  endif
      END DO
# endif
!
!-----------------------------------------------------------------------
!  If processing nonlinear trajectory, load observation reject/accept
!  flag into screening variable.  This needs to be done once and it is
!  controlled by the main driver. Notice that in routine "extract_obs"
!  the observations are only accepted if they are at the appropriate
!  time window and grid bounded at water points.
!-----------------------------------------------------------------------
!
      IF (wrtobsscale(ng).and.wrtnlmod(ng)) THEN
        DO iobs=mstr,mend
          obsscale(iobs)=obsvetting(iobs)
        END DO
      END IF
 
# if defined BGQC && !defined I4DVAR_ANA_SENSITIVITY
!
!-----------------------------------------------------------------------
!  Perform the background quality control: check and reject observations
!  that fail. Only observations that are bounded in time and space
!  (ObsScale .ne. 0) are considered.
!-----------------------------------------------------------------------
!
      IF (wrtobsscale(ng).and.wrtnlmod(ng)) THEN
        DO iobs=mstr,mend
          IF (obsscale(iobs).ne.inival) THEN
#  if defined I4DVAR
            df1=(obsval(iobs)-nlmodval(iobs))/bgerr(iobs)
            df2=(1.0_r8/obserr(iobs))/(bgerr(iobs)*bgerr(iobs))
#  elif defined WEAK_CONSTRAINT
#    ifdef R4DVAR
            df1=(obsval(iobs)-tlmodval(iobs))/bgerr(iobs)
            df2=obserr(iobs)/(bgerr(iobs)*bgerr(iobs))
#    else
            df1=(obsval(iobs)-nlmodval(iobs))/bgerr(iobs)
            df2=obserr(iobs)/(bgerr(iobs)*bgerr(iobs))
#    endif
#  endif
            thresh=bgthresh(iobs)*(1.0_r8+df2)
            IF (df1*df1.gt.thresh) THEN
              obsscale(iobs)=0.0_r8
            END IF
          END IF
        END DO
      END IF
# endif
 
# ifdef DISTRIBUTE
!
!-----------------------------------------------------------------------
!  Collect screening variable for all extracted data.
!-----------------------------------------------------------------------
!
#  ifdef WEAK_CONSTRAINT
      ncollect=mend-mstr+1
#  else
      ncollect=mobs
#  endif
      IF (wrtobsscale(ng).and.wrtnlmod(ng)) THEN
        CALL mp_collect (ng, model, ncollect, inival,                   &
#  ifdef WEAK_CONSTRAINT
     &                   obsscale(mstr:))
#  else
     &                   obsscale)
#  endif
      END IF
# endif
 
# ifndef WEAK_CONSTRAINT
!
!-----------------------------------------------------------------------
!  Save observation reject/accept flag into GLOBAL screening variable.
!-----------------------------------------------------------------------
!
!  In the primal formulation, the current time window (or survey)
!  observations are loaded to the working arrays using local indexes
!  (array elements 1:Nobs) as opposed to global indexes in the dual
!  formulation. Recall that the screening variable ObsScale is computed
!  only once to facilitate Background Quality Control on the first pass
!  of the NLM (WrtObsScale=T and wrtNLmod=T). Therefore, we need to
!  load and save values into global array ObsScaleGlobal so it can be
!  used correctly by the TLM, RPM, and ADM.
!
      IF (wrtobsscale(ng).and.wrtnlmod(ng)) THEN
        ic=0
        DO iobs=nstrobs(ng),nendobs(ng)
          ic=ic+1
          obsscaleglobal(iobs)=obsscale(ic)
        END DO
      ELSE
        ic=0
        DO iobs=nstrobs(ng),nendobs(ng)
          ic=ic+1
          obsscale(ic)=obsscaleglobal(iobs)
        END DO
      END IF
# endif
!
!-----------------------------------------------------------------------
!  Scale extracted data at observations location.  Save NLM unvetted
!  values before they are rejected by quality control.  The unvetted
!  vector is used only for output purposes to analyze why a datum was
!  rejected.
!-----------------------------------------------------------------------
!
!  Nonlinear state vector.
!
      IF (wrtnlmod(ng)) THEN
        DO iobs=mstr,mend
          unvetted(iobs)=nlmodval(iobs)
          nlmodval(iobs)=obsscale(iobs)*nlmodval(iobs)
        END DO
      END IF
 
# ifdef TLM_OBS
!
!  Tangent linear state vector.
!
      IF (wrttlmod(ng).or.wrtrpmod(ng)) THEN
        DO iobs=mstr,mend
#  ifdef I4DVAR_ANA_SENSITIVITY
          tlmodval(iobs)=obsscale(iobs)*tlmodval(iobs)*obserr(iobs)
#  else
          tlmodval(iobs)=obsscale(iobs)*tlmodval(iobs)
#  endif
        END DO
      END IF
# endif
 
# ifdef DISTRIBUTE
!
!-----------------------------------------------------------------------
!  Collect extracted data.
!-----------------------------------------------------------------------
!
#  ifndef I4DVAR_ANA_SENSITIVITY
!
!  Nonlinear state vector.
!
      IF (wrtnlmod(ng)) THEN
        CALL mp_collect (ng, model, ncollect, inival,                   &
#   if defined WEAK_CONSTRAINT
     &                   nlmodval(mstr:))
#   else
     &                   nlmodval)
#   endif
!
        CALL mp_collect (ng, model, ncollect, inival,                   &
#   if defined WEAK_CONSTRAINT
     &                   unvetted(mstr:))
#   else
     &                   unvetted)
#   endif
      END IF
 
#   ifndef VERIFICATION
!
!  Background error standard deviation state vector.
!
      IF (wrtobsscale(ng).and.wrtnlmod(ng)) THEN
        CALL mp_collect (ng, model, ncollect, inival,                   &
#    if defined WEAK_CONSTRAINT
     &                   bgerr(mstr:))
#    else
     &                   bgerr)
#    endif
      END IF
#   endif
#  endif
 
#  ifdef TLM_OBS
!
!  Tangent linear state vector.
!
      IF (wrttlmod(ng).or.wrtrpmod(ng)) THEN
        CALL mp_collect (ng, model, ncollect, inival,                   &
#   if defined WEAK_CONSTRAINT
     &                   tlmodval(mstr:))
#   else
     &                   tlmodval)
#   endif
      END IF
#  endif
 
#  ifdef SOLVE3D
!
!  State vector depths.
!
      IF (load_zobs(ng)) THEN
        CALL mp_collect (ng, model, ncollect, inival,                   &
#   ifdef WEAK_CONSTRAINT
     &                   zobs(mstr:))
#   else
     &                   zobs)
#   endif
      END IF
#  endif
# endif
 
      RETURN
      END SUBROUTINE obs_operator
!
!***********************************************************************
      SUBROUTINE obs_write_nf90 (ng, tile, model,                       &
     &                           LBi, UBi, LBj, UBj,                    &
     &                           Mstr, Mend)
!***********************************************************************
!
      USE mod_netcdf
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile, model
      integer, intent(in) :: lbi, ubi, lbj, ubj
      integer, intent(in) :: mstr, mend
!
!  Local variable declarations.
!
      logical :: lwrote
!
      integer :: tindex, iobs, status
!
      real(dp) :: scale
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", obs_write_nf90"
!
      sourcefile=myfile
 
# if defined FOUR_DVAR && !defined I4DVAR_ANA_SENSITIVITY
!
!-----------------------------------------------------------------------
!  Compute and write initial and final model-observation misfit
!  (innovation) vector for output purposes only. Write also initial
!  nonlinear model at observation locations.
!-----------------------------------------------------------------------
!
      IF (wrtmisfit(ng)) THEN
        DO iobs=mstr,mend
#  if defined I4DVAR
          misfit(iobs)=obsscale(iobs)*sqrt(obserr(iobs))*               &
     &                 (nlmodval(iobs)-obsval(iobs))
#  elif defined WEAK_CONSTRAINT
#   ifdef R4DVAR
          misfit(iobs)=obsscale(iobs)/sqrt(obserr(iobs))*               &
     &                 (tlmodval(iobs)-obsval(iobs))
#   else
          misfit(iobs)=obsscale(iobs)/sqrt(obserr(iobs))*               &
     &                   (nlmodval(iobs)-obsval(iobs))
#   endif
#  endif
        END DO
!
!  Write out initial or final values into DAV NetCDF file.
!
        IF (nrun.eq.1) THEN
          CALL netcdf_put_fvar (ng, model, dav(ng)%name,                &
     &                          vname(1,idnlmp),                        &
     &                          unvetted(mstr:),                        &
     &                          (/nstrobs(ng)/), (/nobs(ng)/),          &
     &                          ncid = dav(ng)%ncid,                    &
     &                          varid = dav(ng)%Vid(idnlmp))
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
          CALL netcdf_put_fvar (ng, model, dav(ng)%name,                &
     &                          vname(1,idnlmi),                        &
     &                          nlmodval(mstr:),                        &
     &                          (/nstrobs(ng)/), (/nobs(ng)/),          &
     &                          ncid = dav(ng)%ncid,                    &
     &                          varid = dav(ng)%Vid(idnlmi))
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
          CALL netcdf_put_fvar (ng, model, dav(ng)%name,                &
     &                          vname(1,idmomi),                        &
     &                          misfit(mstr:),                          &
     &                          (/nstrobs(ng)/), (/nobs(ng)/),          &
     &                          ncid = dav(ng)%ncid,                    &
     &                          varid = dav(ng)%Vid(idmomi))
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        ELSE
          CALL netcdf_put_fvar (ng, model, dav(ng)%name,                &
     &                          vname(1,idmomf),                        &
     &                          misfit(mstr:),                          &
     &                          (/nstrobs(ng)/), (/nobs(ng)/),          &
     &                          ncid = dav(ng)%ncid,                    &
     &                          varid = dav(ng)%Vid(idmomf))
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
      END IF
 
#  ifdef BGQC
!
!  Write out bacground error standard deviation.
!
      IF (wrtobsscale(ng).and.wrtnlmod(ng)) THEN
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        vname(1,idbgth),                          &
     &                        bgthresh(mstr:),                          &
     &                        (/nstrobs(ng)/), (/nobs(ng)/),            &
     &                        ncid = dav(ng)%ncid,                      &
     &                        varid = dav(ng)%Vid(idbgth))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
#  endif
# endif
 
# if defined FOUR_DVAR && !defined I4DVAR_ANA_SENSITIVITY
!
!-----------------------------------------------------------------------
!  Write out variables needed to compute Desroziers et al. (2005)
!  statistics.
!-----------------------------------------------------------------------
!
      IF (wrtnlmod(ng)) THEN
!
! 4D-Var background error standard deviation at observation locations.
!
        IF (wrtobsscale(ng)) THEN
          CALL netcdf_put_fvar (ng, model, dav(ng)%name,                &
     &                          vname(1,idbger),                        &
     &                          bgerr(mstr:),                           &
     &                          (/nstrobs(ng)/), (/nobs(ng)/),          &
     &                          ncid = dav(ng)%ncid,                    &
     &                          varid = dav(ng)%Vid(idbger))
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
!
!  4D-Var innovation vector: observation minus background. Save initial
!  background in the NLincrement vector.
!
        IF (nrun.eq.1) THEN
          DO iobs=mstr,mend
            innovation(iobs)=obsscale(iobs)*                            &
     &                       (obsval(iobs)-nlmodval(iobs))
            nlincrement(iobs)=nlmodval(iobs)
          END DO
!
          CALL netcdf_put_fvar (ng, model, dav(ng)%name,                &
     &                          vname(1,idinno),                        &
     &                          innovation(mstr:),                      &
     &                          (/nstrobs(ng)/), (/nobs(ng)/),          &
     &                          ncid = dav(ng)%ncid,                    &
     &                          varid = dav(ng)%Vid(idinno))
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
!
!  4D-Var increment (analysis minus background) and residual
!  (observation minus analysis).
!
        IF (outer.eq.nouter) THEN
          DO iobs=mstr,mend
            nlincrement(iobs)=obsscale(iobs)*                           &
     &                        (nlmodval(iobs)-nlincrement(iobs))
            residual(iobs)=obsscale(iobs)*                              &
     &                     (obsval(iobs)-nlmodval(iobs))
          END DO
!
          CALL netcdf_put_fvar (ng, model, dav(ng)%name,                &
     &                          vname(1,idincr),                        &
     &                          nlincrement(mstr:),                     &
     &                          (/nstrobs(ng)/), (/nobs(ng)/),          &
     &                          ncid = dav(ng)%ncid,                    &
     &                          varid = dav(ng)%Vid(idincr))
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
          CALL netcdf_put_fvar (ng, model, dav(ng)%name,                &
     &                          vname(1,idresi),                        &
     &                          residual(mstr:),                        &
     &                          (/nstrobs(ng)/), (/nobs(ng)/),          &
     &                          ncid = dav(ng)%ncid,                    &
     &                          varid = dav(ng)%Vid(idresi))
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
      END IF
# endif
!
!-----------------------------------------------------------------------
!  Write out other variable into output DAC NetCDF file.
!-----------------------------------------------------------------------
 
# if defined I4DVAR                  || defined WEAK_CONSTRAINT && \
    !defined I4DVAR_ANA_SENSITIVITY
!
!  Current outer and inner loop.
!
      IF (wrtnlmod(ng).or.wrttlmod(ng).or.wrtrpmod(ng)) THEN
        CALL netcdf_put_ivar (ng, model, dav(ng)%name, 'outer',         &
     &                        outer, (/0/), (/0/),                      &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
        CALL netcdf_put_ivar (ng, model, dav(ng)%name, 'inner',         &
     &                        inner, (/0/), (/0/),                      &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
#endif
!
!  Observation screening flag: accept or reject.
!
      IF (wrtobsscale(ng).and.wrtnlmod(ng)) THEN
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        vname(1,idobss),                          &
     &                        obsscale(mstr:),                          &
     &                        (/nstrobs(ng)/), (/nobs(ng)/),            &
     &                        ncid = dav(ng)%ncid,                      &
     &                        varid = dav(ng)%Vid(idobss))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
 
# ifndef I4DVAR_ANA_SENSITIVITY
!
!  Nonlinear model or first guess (background) state at observation
!  locations. If 4D-Var, the values are written for every outer loop.
!  Also, the final value is written as individual variable to
!  facilitate ERDDAP processing. Notice that it is always overwritten
!  to avoid complicated logic with all the 4D-Var drivers. The last
!  one written correspont to the final values.
!
#  if defined VERIFICATION || defined TLM_CHECK
      IF (wrtnlmod(ng)) THEN
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        vname(1,idnlmo),                          &
     &                        nlmodval(mstr:),                          &
     &                        (/nstrobs(ng)/), (/nobs(ng)/),            &
     &                        ncid = dav(ng)%ncid,                      &
     &                        varid = dav(ng)%Vid(idnlmo))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        havenlmod(ng)=.true.
      END IF
#  elif defined I4DVAR || defined WEAK_CONSTRAINT
      IF (wrtnlmod(ng).and.outer.gt.0) THEN
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        vname(1,idnlmo),                          &
     &                        nlmodval(mstr:),                          &
     &                        (/nstrobs(ng),outer/), (/nobs(ng),1/),    &
     &                        ncid = dav(ng)%ncid,                      &
     &                        varid = dav(ng)%Vid(idnlmo))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
#  endif
#  ifdef FOUR_DVAR
!
      IF (wrtnlmod(ng).and.outer.gt.0) THEN
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        vname(1,idnlmf),                          &
     &                        nlmodval(mstr:),                          &
     &                        (/nstrobs(ng)/), (/nobs(ng)/),            &
     &                        ncid = dav(ng)%ncid,                      &
     &                        varid = dav(ng)%Vid(idnlmf))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        havenlmod(ng)=.true.
      END IF
#  endif
# endif
# if !defined I4DVAR_ANA_SENSITIVITY && \
     (defined i4dvar                 || defined weak_constraint)
!
!  Write out unvetted nonlinear model at the observation locations per
!  outer loop.
!
      IF (wrtnlmod(ng).and.outer.gt.0) THEN
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        vname(1,idnlmu),                          &
     &                        unvetted(mstr:),                          &
     &                        (/nstrobs(ng),outer/), (/nobs(ng),1/),    &
     &                        ncid = dav(ng)%ncid,                      &
     &                        varid = dav(ng)%Vid(idnlmu))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
# endif
# if defined I4DVAR          || defined I4DVAR_ANA_SENSITIVITY || \
     defined weak_constraint
!
!  Tangent linear model state increments at observation locations.
!
      IF (wrttlmod(ng).or.wrtrpmod(ng)) THEN
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        vname(1,idtlmo),                          &
     &                        tlmodval(mstr:),                          &
     &                        (/nstrobs(ng)/), (/nobs(ng)/),            &
     &                        ncid = dav(ng)%ncid,                      &
     &                        varid = dav(ng)%Vid(idtlmo))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        havetlmod(ng)=.true.
      END IF
# endif
# ifdef TLM_OBS
#  ifdef I4DVAR_ANA_SENSITIVITY
#   ifdef OBS_IMPACT
!
!  Write the total observation impact.
!
      IF (wrtimpact_tot(ng)) THEN
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        'ObsImpact_total',                        &
     &                        tlmodval(mstr:),                          &
     &                        (/nstrobs(ng)/), (/nobs(ng)/),            &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
#   endif
 
#   ifdef OBS_IMPACT_SPLIT
!
!  Write the observation impact of initial conditions.
!
      IF (wrtimpact_ic(ng)) THEN
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        'ObsImpact_IC',                           &
     &                        tlmodval(mstr:),                          &
     &                        (/nstrobs(ng)/), (/nobs(ng)/),            &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
 
#    if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
!
!  Write the observation impact of surface forcing.
!
      IF (wrtimpact_fc(ng)) THEN
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        'ObsImpact_FC',                           &
     &                        tlmodval(mstr:),                          &
     &                        (/nstrobs(ng)/), (/nobs(ng)/),            &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
#    endif
 
#    if defined ADJUST_BOUNDARY
!
!  Write the observation impact of boundary conditions.
!
      IF (wrtimpact_bc(ng)) THEN
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        'ObsImpact_BC',                           &
     &                        tlmodval(mstr:),                          &
     &                        (/nstrobs(ng)/), (/nobs(ng)/),            &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
#    endif
#   endif
#  endif
# endif
 
# if defined R4DVAR    || defined R4DVAR_ANA_SENSITIVITY || \
     defined tl_r4dvar
!
!  Write initial representer model increments at observation locations.
!
      IF (nrun.eq.1.and.wrtrpmod(ng)) THEN
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        'RPmodel_initial',                        &
     &                        tlmodval(mstr:),                          &
     &                        (/nstrobs(ng)/), (/nobs(ng)/) ,           &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
# endif
 
# ifdef SOLVE3D
!
!  Now, that Zobs has all the Z-locations in grid coordinates, write
!  "obs_Zgrid" into data assimilation output file.
!
      IF ((nrun.eq.1).and.                                              &
     &    (wrtnlmod(ng).or.wrttlmod(ng).or.wrtrpmod(ng))) THEN
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        vname(1,idobsz),                          &
     &                        zobs(mstr:),                              &
     &                        (/nstrobs(ng)/), (/nobs(ng)/),            &
     &                        ncid = dav(ng)%ncid,                      &
     &                        varid = dav(ng)%Vid(idobsz))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Write Z-location of observation in grid coordinates, if applicable.
!  This values are needed elsewhere when using the interpolation
!  weight matrix.  Recall that the depth of observations can be in
!  meters or grid coordinates.  Notice that since the model levels
!  evolve in time, the fractional level coordinate is unknow during
!  the processing of the observations.
!
      IF ((nrun.eq.1).and.                                              &
     &    (wrtnlmod(ng).or.wrttlmod(ng).or.wrtrpmod(ng))) THEN
        CALL netcdf_put_fvar (ng, model, obs(ng)%name,                  &
     &                        vname(1,idobsz),                          &
     &                        zobs(mstr:),                              &
     &                        (/nstrobs(ng)/), (/nobs(ng)/),            &
     &                        ncid = obs(ng)%ncid,                      &
     &                        varid = obs(ng)%Vid(idobsz))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        IF (model.eq.iadm) THEN
          lwrote=obssurvey(ng).eq.1
        ELSE
          lwrote=obssurvey(ng).eq.nsurvey(ng)
        END IF
        IF (lwrote) wrote_zobs(ng)=.true.
      END IF
# endif
 
# if defined BALANCE_OPERATOR && defined ZETA_ELLIPTIC
!
!  Write out reference free-surface used in the balance operator.
!
      IF (wrtzetaref(ng).and.balance(isfsur)) THEN
        tindex=0
        scale=1.0_dp
        status=nf_fwrite2d(ng, model, dav(ng)%ncid, idfsur,             &
     &                     dav(ng)%Vid(idfsur), tindex, r2dvar,         &
     &                     lbi, ubi, lbj, ubj, scale,                   &
#  ifdef MASKING
     &                     grid(ng) % rmask,                            &
#  endif
     &                     fourdvar(ng) % zeta_ref,                     &
     &                     setfillval = .false.)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        wrtzetaref(ng)=.false.
      END IF
# endif
!
!-----------------------------------------------------------------------
!  Synchronize observations NetCDF file to disk.
!-----------------------------------------------------------------------
!
      IF (wrtnlmod(ng).or.wrttlmod(ng).or.wrtrpmod(ng)) THEN
        CALL netcdf_sync (ng, model, dav(ng)%name, dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
# ifdef SOLVE3D
        IF (nrun.eq.1) THEN
          CALL netcdf_sync (ng, model, obs(ng)%name, obs(ng)%ncid)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
# endif
      END IF
!
      RETURN
      END SUBROUTINE obs_write_nf90
 
# if defined PIO_LIB && defined DISTRIBUTE
!
!***********************************************************************
      SUBROUTINE obs_write_pio (ng, tile, model,                        &
     &                          LBi, UBi, LBj, UBj,                     &
     &                          Mstr, Mend)
!***********************************************************************
!
      USE mod_pio_netcdf
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile, model
      integer, intent(in) :: lbi, ubi, lbj, ubj
      integer, intent(in) :: mstr, mend
!
!  Local variable declarations.
!
      logical :: lwrote
!
      integer :: tindex, iobs, status
!
      real(dp) :: scale
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", obs_write_pio"
 
#  if defined BALANCE_OPERATOR && defined ZETA_ELLIPTIC
!
      TYPE (io_desc_t), pointer :: iodesc
#  endif
!
      sourcefile=myfile
 
#  if defined FOUR_DVAR && !defined I4DVAR_ANA_SENSITIVITY
!
!-----------------------------------------------------------------------
!  Compute and write initial and final model-observation misfit
!  (innovation) vector for output purposes only. Write also initial
!  nonlinear model at observation locations.
!-----------------------------------------------------------------------
!
      IF (wrtmisfit(ng)) THEN
        DO iobs=mstr,mend
#   if defined I4DVAR
          misfit(iobs)=obsscale(iobs)*sqrt(obserr(iobs))*               &
     &                 (nlmodval(iobs)-obsval(iobs))
#   elif defined WEAK_CONSTRAINT
#    ifdef R4DVAR
          misfit(iobs)=obsscale(iobs)/sqrt(obserr(iobs))*               &
     &                 (tlmodval(iobs)-obsval(iobs))
#    else
          misfit(iobs)=obsscale(iobs)/sqrt(obserr(iobs))*               &
     &                   (nlmodval(iobs)-obsval(iobs))
#    endif
#   endif
        END DO
!
!  Write out initial or final values into DAV NetCDF file.
!
        IF (nrun.eq.1) THEN
          CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,            &
     &                              vname(1,idnlmp),                    &
     &                              unvetted(mstr:mend),                &
     &                              (/nstrobs(ng)/), (/nobs(ng)/),   &
     &                              piofile = dav(ng)%pioFile,          &
     &                              piovar = dav(ng)%pioVar(idnlmp)%vd)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
          CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,            &
     &                              vname(1,idnlmi),                    &
     &                              nlmodval(mstr:mend),                &
     &                              (/nstrobs(ng)/), (/nobs(ng)/),      &
     &                              piofile = dav(ng)%pioFile,          &
     &                              piovar = dav(ng)%pioVar(idnlmi)%vd)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
          CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,            &
     &                              vname(1,idmomi),                    &
     &                              misfit(mstr:mend),                  &
     &                              (/nstrobs(ng)/), (/nobs(ng)/),      &
     &                              piofile = dav(ng)%pioFile,          &
     &                              piovar = dav(ng)%pioVar(idmomi)%vd)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        ELSE
          CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,            &
     &                              vname(1,idmomf),                    &
     &                              misfit(mstr:mend),                  &
     &                              (/nstrobs(ng)/), (/nobs(ng)/),      &
     &                              piofile = dav(ng)%pioFile,          &
     &                              piovar = dav(ng)%pioVar(idmomf)%vd)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
      END IF
 
#   ifdef BGQC
!
!  Write out bacground error standard deviation threshold.
!
      IF (wrtobsscale(ng).and.wrtnlmod(ng)) THEN
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            vname(1,idbgth),                      &
     &                            bgthresh(mstr:mend),                  &
     &                            (/nstrobs(ng)/), (/nobs(ng)/),        &
     &                            piofile = dav(ng)%pioFile,            &
     &                            piovar = dav(ng)%pioVar(idbgth)%vd)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
#   endif
#  endif
 
#  if defined FOUR_DVAR && !defined I4DVAR_ANA_SENSITIVITY
!
!-----------------------------------------------------------------------
!  Write out variables needed to compute Desroziers et al. (2005)
!  statistics.
!-----------------------------------------------------------------------
!
      IF (wrtnlmod(ng)) THEN
!
! 4D-Var background error standard deviation at observation locations.
!
        IF (wrtobsscale(ng)) THEN
          CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,            &
     &                              vname(1,idbger),                    &
     &                              bgerr(mstr:mend),                   &
     &                              (/nstrobs(ng)/), (/nobs(ng)/),      &
     &                              piofile = dav(ng)%pioFile,          &
     &                              piovar = dav(ng)%pioVar(idbger)%vd)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
!
!  4D-Var innovation vector: observation minus background. Save initial
!  background in the NLincrement vector.
!
        IF (nrun.eq.1) THEN
          DO iobs=mstr,mend
            innovation(iobs)=obsscale(iobs)*                            &
     &                       (obsval(iobs)-nlmodval(iobs))
            nlincrement(iobs)=nlmodval(iobs)
          END DO
!
          CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,            &
     &                              vname(1,idinno),                    &
     &                              innovation(mstr:mend),              &
     &                              (/nstrobs(ng)/), (/nobs(ng)/),      &
     &                              piofile = dav(ng)%pioFile,          &
     &                              piovar = dav(ng)%pioVar(idinno)%vd)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
!
!  4D-Var increment (analysis minus background) and residual
!  (observation minus analysis).
!
        IF (outer.eq.nouter) THEN
          DO iobs=mstr,mend
            nlincrement(iobs)=obsscale(iobs)*                           &
     &                        (nlmodval(iobs)-nlincrement(iobs))
            residual(iobs)=obsscale(iobs)*                              &
     &                     (obsval(iobs)-nlmodval(iobs))
          END DO
!
          CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,            &
     &                              vname(1,idincr),                    &
     &                              nlincrement(mstr:mend),             &
     &                              (/nstrobs(ng)/), (/nobs(ng)/),      &
     &                              piofile = dav(ng)%pioFile,          &
     &                              piovar = dav(ng)%pioVar(idincr)%vd)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
          CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,            &
     &                              vname(1,idresi),                    &
     &                              residual(mstr:mend),                &
     &                              (/nstrobs(ng)/), (/nobs(ng)/),      &
     &                              piofile = dav(ng)%pioFile,          &
     &                              piovar = dav(ng)%pioVar(idresi)%vd)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
      END IF
#  endif
!
!-----------------------------------------------------------------------
!  Write out other variable into output DAC NetCDF file.
!-----------------------------------------------------------------------
 
# if defined I4DVAR                  || defined WEAK_CONSTRAINT && \
    !defined I4DVAR_ANA_SENSITIVITY
!
!  Current outer and inner loop.
!
      IF (wrtnlmod(ng).or.wrttlmod(ng).or.wrtrpmod(ng)) THEN
        CALL pio_netcdf_put_ivar (ng, model, dav(ng)%name, 'outer',     &
     &                            outer, (/0/), (/0/),                  &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
        CALL pio_netcdf_put_ivar (ng, model, dav(ng)%name, 'inner',     &
     &                            inner, (/0/), (/0/),                  &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
# endif
!
!  Observation screening flag: accept or reject.
!
      IF (wrtobsscale(ng).and.wrtnlmod(ng)) THEN
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            vname(1,idobss),                      &
     &                            obsscale(mstr:mend),                  &
     &                            (/nstrobs(ng)/), (/nobs(ng)/),        &
     &                            piofile = dav(ng)%pioFile,            &
     &                            piovar = dav(ng)%pioVar(idobss)%vd)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
 
#  ifndef I4DVAR_ANA_SENSITIVITY
!
!  Nonlinear model or first guess (background) state at observation
!  locations. If 4D-Var, the values are written for every outer loop.
!  Also, the final value is written as individual variable to
!  facilitate ERDDAP processing. Notice that it is always overwritten
!  to avoid complicated logic with all the 4D-Var drivers. The last
!  one written correspont to the final values.
!
#   if defined VERIFICATION || defined TLM_CHECK
      IF (wrtnlmod(ng)) THEN
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            vname(1,idnlmo),                      &
     &                            nlmodval(mstr:mend),                  &
     &                            (/nstrobs(ng)/), (/nobs(ng)/),        &
     &                            piofile = dav(ng)%pioFile,            &
     &                            piovar = dav(ng)%pioVar(idnlmo)%vd)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        havenlmod(ng)=.true.
      END IF
#   elif defined I4DVAR || defined WEAK_CONSTRAINT
      IF (wrtnlmod(ng).and.outer.gt.0) THEN
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            vname(1,idnlmo),                      &
     &                            nlmodval(mstr:mend),                  &
     &                            (/nstrobs(ng),outer/), (/nobs(ng),1/),&
     &                            piofile = dav(ng)%pioFile,            &
     &                            piovar = dav(ng)%pioVar(idnlmo)%vd)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
#   endif
#   ifdef FOUR_DVAR
!
      IF (wrtnlmod(ng).and.outer.gt.0) THEN
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            vname(1,idnlmf),                      &
     &                            nlmodval(mstr:mend),                  &
     &                            (/nstrobs(ng)/), (/nobs(ng)/),        &
     &                            piofile = dav(ng)%pioFile,            &
     &                            piovar = dav(ng)%pioVar(idnlmf)%vd)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        havenlmod(ng)=.true.
      END IF
#   endif
#  endif
#  if !defined I4DVAR_ANA_SENSITIVITY && \
     (defined i4dvar                 || defined weak_constraint)
!
!  Write out unvetted nonlinear model at the observation locations per
!  outer loop.
!
      IF (wrtnlmod(ng).and.outer.gt.0) THEN
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            vname(1,idnlmu),                      &
     &                            unvetted(mstr:mend),                  &
     &                            (/nstrobs(ng),outer/), (/nobs(ng),1/),&
     &                            piofile = dav(ng)%pioFile,            &
     &                            piovar = dav(ng)%pioVar(idnlmu)%vd)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
#  endif
#  if defined I4DVAR          || defined I4DVAR_ANA_SENSITIVITY || \
      defined weak_constraint
!
!  Tangent linear model state increments at observation locations.
!
      IF (wrttlmod(ng).or.wrtrpmod(ng)) THEN
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            vname(1,idtlmo),                      &
     &                            tlmodval(mstr:mend),                  &
     &                            (/nstrobs(ng)/), (/nobs(ng)/),        &
     &                            piofile = dav(ng)%pioFile,            &
     &                            piovar = dav(ng)%pioVar(idtlmo)%vd)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        havetlmod(ng)=.true.
      END IF
#  endif
#  ifdef TLM_OBS
#   ifdef I4DVAR_ANA_SENSITIVITY
#    ifdef OBS_IMPACT
!
!  Write the total observation impact.
!
      IF (wrtimpact_tot(ng)) THEN
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            'ObsImpact_total',                    &
     &                            tlmodval(mstr:mend),                  &
     &                            (/nstrobs(ng)/), (/nobs(ng)/),        &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
#    endif
 
#    ifdef OBS_IMPACT_SPLIT
!
!  Write the observation impact of initial conditions.
!
      IF (wrtimpact_ic(ng)) THEN
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            'ObsImpact_IC',                       &
     &                            tlmodval(mstr:mend),                  &
     &                            (/nstrobs(ng)/), (/nobs(ng)/),        &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
 
#     if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
!
!  Write the observation impact of surface forcing.
!
      IF (wrtimpact_fc(ng)) THEN
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            'ObsImpact_FC',                       &
     &                            tlmodval(mstr:mend),                  &
     &                            (/nstrobs(ng)/), (/nobs(ng)/),        &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
#     endif
 
#     if defined ADJUST_BOUNDARY
!
!  Write the observation impact of boundary conditions.
!
      IF (wrtimpact_bc(ng)) THEN
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            'ObsImpact_BC',                       &
     &                            tlmodval(mstr:mend),                  &
     &                            (/nstrobs(ng)/), (/nobs(ng)/),        &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
#     endif
#    endif
#   endif
#  endif
 
#  if defined R4DVAR    || defined R4DVAR_ANA_SENSITIVITY || \
     defined tl_r4dvar
!
!  Write initial representer model increments at observation locations.
!
      IF (nrun.eq.1.and.wrtrpmod(ng)) THEN
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            'RPmodel_initial',                    &
     &                            tlmodval(mstr:mend),                  &
     &                            (/nstrobs(ng)/), (/nobs(ng)/) ,       &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
#  endif
 
#  ifdef SOLVE3D
!
!  Now, that Zobs has all the Z-locations in grid coordinates, write
!  "obs_Zgrid" into data assimilation output file.
!
      IF ((nrun.eq.1).and.                                              &
     &    (wrtnlmod(ng).or.wrttlmod(ng).or.wrtrpmod(ng))) THEN
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            vname(1,idobsz),                      &
     &                            zobs(mstr:mend),                      &
     &                            (/nstrobs(ng)/), (/nobs(ng)/),        &
     &                            piofile = dav(ng)%pioFile,            &
     &                            piovar = dav(ng)%pioVar(idobsz)%vd)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Write Z-location of observation in grid coordinates, if applicable.
!  This values are needed elsewhere when using the interpolation
!  weight matrix.  Recall that the depth of observations can be in
!  meters or grid coordinates.  Notice that since the model levels
!  evolve in time, the fractional level coordinate is unknow during
!  the processing of the observations.
!
      IF ((nrun.eq.1).and.                                              &
     &    (wrtnlmod(ng).or.wrttlmod(ng).or.wrtrpmod(ng))) THEN
        CALL pio_netcdf_put_fvar (ng, model, obs(ng)%name,              &
     &                            vname(1,idobsz),                      &
     &                            zobs(mstr:mend),                      &
     &                            (/nstrobs(ng)/), (/nobs(ng)/),        &
     &                            piofile = obs(ng)%pioFile,            &
     &                            piovar = obs(ng)%pioVar(idobsz)%vd)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        IF (model.eq.iadm) THEN
          lwrote=obssurvey(ng).eq.1
        ELSE
          lwrote=obssurvey(ng).eq.nsurvey(ng)
        END IF
        IF (lwrote) wrote_zobs(ng)=.true.
      END IF
#  endif
 
#  if defined BALANCE_OPERATOR && defined ZETA_ELLIPTIC
!
!  Write out reference free-surface used in the balance operator.
!
      IF (wrtzetaref(ng).and.balance(isfsur)) THEN
        IF (dav(ng)%pioVar(idfsur)%dkind.eq.pio_double) THEN
          iodesc => iodesc_dp_r2dvar(ng)
        ELSE
          iodesc => iodesc_sp_r2dvar(ng)
        END IF
        tindex=0
        scale=1.0_dp
        status=nf_fwrite2d(ng, model, dav(ng)%pioFile, idfsur,          &
     &                     dav(ng)%pioVar(idfsur), tindex,              &
     &                     iodesc,                                      &
     &                     lbi, ubi, lbj, ubj, scale,                   &
#   ifdef MASKING
     &                     grid(ng) % rmask,                            &
#   endif
     &                     fourdvar(ng) % zeta_ref,                     &
     &                     setfillval = .false.)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        wrtzetaref(ng)=.false.
      END IF
#  endif
!
!-----------------------------------------------------------------------
!  Synchronize observations NetCDF file to disk.
!-----------------------------------------------------------------------
!
      IF (wrtnlmod(ng).or.wrttlmod(ng).or.wrtrpmod(ng)) THEN
        CALL pio_netcdf_sync (ng, model, dav(ng)%name,                  &
     &                        dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
#  ifdef SOLVE3D
        IF (nrun.eq.1) THEN
          CALL pio_netcdf_sync (ng, model, obs(ng)%name,                &
     &                          obs(ng)%pioFile)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
#  endif
      END IF
!
      RETURN
      END SUBROUTINE obs_write_pio
# endif
#endif
      END MODULE obs_write_mod