#include "cppdefs.h"
#if (defined FOUR_DVAR || defined VERIFICATION) && defined OBSERVATIONS
      SUBROUTINE stats_modobs (ng)
!
!svn $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2020 The ROMS/TOMS Group                         !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.txt                                              !
!=======================================================================
!                                                                      !
!  This subroutine computes several statistical quantities between     !
!  model and observations:                                             !
!                                                                      !
!     CC         Cross-Correlation                                     !
!     MB         Model Bias                                            !
!     MSE        Mean Squared Error                                    !
!     SDE        Standard Deviation Error                              !
!                                                                      !
!  Reference:                                                          !
!                                                                      !
!  Oke, P.R., J.S. Allen, R.N. Miller, G.D. Egbert, J.A. Austin,       !
!    J.A. Barth, T.J. Boyd, P.M. Kosro, and M.D. Levine, 2002: A       !
!    Modeling Study of the Three-Dimensional Continental Shelf         !
!    Circulation off Oregon. Part I: Model-Data Comparison, J.         !
!    Phys. Oceanogr., 32, 1360-1382.                                   !
!                                                                      !
!  Notice that this routine is never called inside of a parallel       !
!  region. Therefore, parallel reductions are not required.            !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_fourdvar
      USE mod_grid
      USE mod_iounits
      USE mod_ncparam
      USE mod_netcdf
      USE mod_scalars
!
# ifdef DISTRIBUTE
      USE distribute_mod, ONLY :  mp_collect
# endif
      USE obs_k2z_mod, ONLY : obs_k2z
      USE strings_mod, ONLY : FoundError, find_string
!
      implicit none
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng
!
!  Local variable declarations.
!
      integer :: i, ic, iobs, status, varid, vindex
      integer :: LBi, UBi, LBj, UBj
      integer, dimension(2) :: start, total
      integer, dimension(NobsVar(ng)) :: Ncount, is

      integer, dimension(Ndatum(ng)) :: obs_type
      integer, dimension(Ndatum(ng)) :: provenance
      integer, dimension(Nsurvey(ng)) :: survey_Nobs

      real(r8) :: cff1, cff2
      real(r8), parameter :: IniVal = 0.0_r8

      real(r8), dimension(NobsVar(ng)) :: CC, MB, MSE, SDE
      real(r8), dimension(NobsVar(ng)) :: mod_min, mod_max
      real(r8), dimension(NobsVar(ng)) :: mod_mean, mod_std
      real(r8), dimension(NobsVar(ng)) :: obs_min, obs_max
      real(r8), dimension(NobsVar(ng)) :: obs_mean, obs_std

      real(r8), dimension(Ndatum(ng)) :: mod_value
      real(r8), dimension(Ndatum(ng)) :: obs_depths
      real(r8), dimension(Ndatum(ng)) :: obs_scale
      real(dp), dimension(Ndatum(ng)) :: obs_time
      real(r8), dimension(Ndatum(ng)) :: obs_value
      real(r8), dimension(Ndatum(ng)) :: obs_Xgrid
      real(r8), dimension(Ndatum(ng)) :: obs_Ygrid
      real(r8), dimension(Ndatum(ng)) :: obs_work
      real(dp), dimension(Nsurvey(ng)) :: survey

      character (len=11), dimension(NobsVar(ng)) :: text, svar_name

      character (len=30) :: F1, F2, F3, F4
!
      SourceFile=__FILE__
!
      LBi=LBOUND(GRID(ng)%h,DIM=1)
      UBi=UBOUND(GRID(ng)%h,DIM=1)
      LBj=LBOUND(GRID(ng)%h,DIM=2)
      UBj=UBOUND(GRID(ng)%h,DIM=2)
!
!-----------------------------------------------------------------------
!  Read in model and observations data.
!-----------------------------------------------------------------------
!
!  Inquire about input observations variables.
!
      CALL netcdf_inq_var (ng, iNLM, OBS(ng)%name,                      &
     &                     ncid = OBS(ng)%ncid)
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN
!
!  Read in and write out number of observations per survey.
!
      CALL netcdf_get_ivar (ng, iNLM, OBS(ng)%name, Vname(1,idNobs),    &
     &                      survey_Nobs,                                &
     &                      ncid = OBS(ng)%ncid,                        &
     &                      start = (/1/), total = (/Nsurvey(ng)/))
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN

      CALL netcdf_put_ivar (ng, iNLM, DAV(ng)%name, Vname(1,idNobs),    &
     &                      survey_Nobs, (/1/), (/Nsurvey(ng)/),        &
     &                      ncid = DAV(ng)%ncid)
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN
!
!  Read in and write out survey time.
!
      CALL netcdf_get_time (ng, iNLM, OBS(ng)%name, Vname(1,idOday),    &
     &                      Rclock%DateNumber, survey,                  &
     &                      ncid = OBS(ng)%ncid,                        &
     &                      start = (/1/), total = (/Nsurvey(ng)/))
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN

      CALL netcdf_put_fvar (ng, iNLM, DAV(ng)%name, Vname(1,idOday),    &
     &                      survey, (/1/), (/Nsurvey(ng)/),             &
     &                      ncid = DAV(ng)%ncid)
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN
!
!  Read in and write out observation type identifier.
!
      CALL netcdf_get_ivar (ng, iNLM, OBS(ng)%name, Vname(1,idOtyp),    &
     &                      obs_type,                                   &
     &                      ncid = OBS(ng)%ncid,                        &
     &                      start = (/1/), total = (/Ndatum(ng)/))
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN

      CALL netcdf_put_ivar (ng, iNLM, DAV(ng)%name, Vname(1,idOtyp),    &
     &                      obs_type, (/1/), (/Ndatum(ng)/),            &
     &                      ncid = DAV(ng)%ncid)
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN
!
!  Read in and write out observation provenance.
!
      CALL netcdf_get_ivar (ng, iNLM, OBS(ng)%name, Vname(1,idOpro),    &
     &                      provenance,                                 &
     &                      ncid = OBS(ng)%ncid,                        &
     &                      start = (/1/), total = (/Ndatum(ng)/))
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN

      CALL netcdf_put_ivar (ng, iNLM, DAV(ng)%name, Vname(1,idOpro),    &
     &                      provenance, (/1/), (/Ndatum(ng)/),          &
     &                      ncid = DAV(ng)%ncid)
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN
!
!  Read in and write out observation time.
!
      CALL netcdf_get_time (ng, iNLM, OBS(ng)%name, Vname(1,idObsT),    &
     &                      Rclock%DateNumber, obs_time,                &
     &                      ncid = OBS(ng)%ncid,                        &
     &                      start = (/1/), total = (/Ndatum(ng)/))
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN

      CALL netcdf_put_fvar (ng, iNLM, DAV(ng)%name, Vname(1,idObsT),    &
     &                      obs_time, (/1/), (/Ndatum(ng)/),            &
     &                      ncid = DAV(ng)%ncid)
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN
!
!  Read in and write out observation longitude.
!
      IF (find_string(var_name,n_var,Vname(1,idOlon),vindex)) THEN
        CALL netcdf_get_fvar (ng, iNLM, OBS(ng)%name, Vname(1,idOlon),  &
     &                        obs_work,                                 &
     &                        ncid = OBS(ng)%ncid,                      &
     &                        start = (/1/), total = (/Ndatum(ng)/))
        IF (FoundError(exit_flag, NoError, __LINE__,                    &
     &                 __FILE__)) RETURN

        CALL netcdf_put_fvar (ng, iNLM, DAV(ng)%name, Vname(1,idOlon),  &
     &                        obs_work, (/1/), (/Ndatum(ng)/),          &
     &                        ncid = DAV(ng)%ncid)
        IF (FoundError(exit_flag, NoError, __LINE__,                    &
     &                 __FILE__)) RETURN
      END IF
!
!  Read in and write out observation latitude.
!
      IF (find_string(var_name,n_var,Vname(1,idOlon),vindex)) THEN
        CALL netcdf_get_fvar (ng, iNLM, OBS(ng)%name, Vname(1,idOlat),  &
     &                        obs_work,                                 &
     &                        ncid = OBS(ng)%ncid,                      &
     &                        start = (/1/), total = (/Ndatum(ng)/))
        IF (FoundError(exit_flag, NoError, __LINE__,                    &
     &                 __FILE__)) RETURN

        CALL netcdf_put_fvar (ng, iNLM, DAV(ng)%name, Vname(1,idOlat),  &
     &                        obs_work, (/1/), (/Ndatum(ng)/),          &
     &                        ncid = DAV(ng)%ncid)
        IF (FoundError(exit_flag, NoError, __LINE__,                    &
     &                 __FILE__)) RETURN
      END IF
!
!  Read in observation depth.
!
      CALL netcdf_get_fvar (ng, iNLM, OBS(ng)%name, Vname(1,idObsD),    &
     &                      obs_depths,                                 &
     &                      ncid = OBS(ng)%ncid,                        &
     &                      start = (/1/), total = (/Ndatum(ng)/))
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN
!
!  Read in and write out X-fractional coordinate.
!
      CALL netcdf_get_fvar (ng, iNLM, OBS(ng)%name, Vname(1,idObsX),    &
     &                      obs_Xgrid,                                  &
     &                      ncid = OBS(ng)%ncid,                        &
     &                      start = (/1/), total = (/Ndatum(ng)/))
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN

      CALL netcdf_put_fvar (ng, iNLM, DAV(ng)%name, Vname(1,idObsX),    &
     &                      obs_Xgrid, (/1/), (/Ndatum(ng)/),           &
     &                      ncid = DAV(ng)%ncid)
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN
!
!  Read in and write out Y-fractional coordinate.
!
      CALL netcdf_get_fvar (ng, iNLM, OBS(ng)%name, Vname(1,idObsY),    &
     &                      obs_Ygrid,                                  &
     &                      ncid = OBS(ng)%ncid,                        &
     &                      start = (/1/), total = (/Ndatum(ng)/))
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN

      CALL netcdf_put_fvar (ng, iNLM, DAV(ng)%name, Vname(1,idObsY),    &
     &                      obs_Ygrid, (/1/), (/Ndatum(ng)/),           &
     &                      ncid = DAV(ng)%ncid)
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN
!
!  Read in and write out observations total error covariance
!  (instrument + sampling + representation).
!
      CALL netcdf_get_fvar (ng, iNLM, OBS(ng)%name, Vname(1,idOerr),    &
     &                      obs_work,                                   &
     &                      ncid = OBS(ng)%ncid,                        &
     &                      start = (/1/), total = (/Ndatum(ng)/))
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN

      CALL netcdf_put_fvar (ng, iNLM, DAV(ng)%name, Vname(1,idOerr),    &
     &                      obs_work, (/1/), (/Ndatum(ng)/),            &
     &                      ncid = DAV(ng)%ncid)
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN
!
!  Read in and write out observation values.
!
      CALL netcdf_get_fvar (ng, iNLM, OBS(ng)%name, Vname(1,idOval),    &
     &                      obs_value,                                  &
     &                      ncid = OBS(ng)%ncid,                        &
     &                      start = (/1/), total = (/Ndatum(ng)/))
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN

      CALL netcdf_put_fvar (ng, iNLM, DAV(ng)%name, Vname(1,idOval),    &
     &                      obs_value, (/1/), (/Ndatum(ng)/),           &
     &                      ncid = DAV(ng)%ncid)
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN
!
!  Read in and write out observation meta values.
!
      IF (haveObsMeta(ng)) THEN
        CALL netcdf_get_fvar (ng, iNLM, OBS(ng)%name, Vname(1,idOmet),  &
     &                        obs_work,                                 &
     &                        ncid = OBS(ng)%ncid,                      &
     &                        start = (/1/), total = (/Ndatum(ng)/))
        IF (FoundError(exit_flag, NoError, __LINE__,                    &
     &                 __FILE__)) RETURN

        CALL netcdf_put_fvar (ng, iNLM, DAV(ng)%name, Vname(1,idOmet),  &
     &                        obs_work, (/1/), (/Ndatum(ng)/),          &
     &                        ncid = DAV(ng)%ncid)
        IF (FoundError(exit_flag, NoError, __LINE__,                    &
     &                 __FILE__)) RETURN
      END IF
!
!  Read in observation screening flag.
!
      CALL netcdf_get_fvar (ng, iNLM, DAV(ng)%name, Vname(1,idObsS),    &
     &                      obs_scale,                                  &
     &                      ncid = DAV(ng)%ncid,                        &
     &                      start = (/1/), total = (/Ndatum(ng)/))
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN
!
!  Read in model values at observation locations.
!
# if defined TL_W4DVAR          || defined W4DVAR || \
     defined W4DVAR_SENSITIVITY
      CALL netcdf_get_fvar (ng, iNLM, DAV(ng)%name, Vname(1,idTLmo),    &
     &                      mod_value,                                  &
     &                      ncid = DAV(ng)%ncid,                        &
     &                      start = (/1/), total = (/Ndatum(ng)/))
# else
#  ifdef VERIFICATION
      CALL netcdf_get_fvar (ng, iNLM, DAV(ng)%name, Vname(1,idNLmo),    &
     &                      mod_value,                                  &
     &                      ncid = DAV(ng)%ncid,                        &
     &                      start = (/1/), total = (/Ndatum(ng)/))
#  else
      CALL netcdf_get_fvar (ng, iNLM, DAV(ng)%name, Vname(1,idNLmo),    &
     &                      mod_value,                                  &
     &                      ncid = DAV(ng)%ncid,                        &
     &                      start = (/1,outer/),                        &
     &                      total = (/Ndatum(ng),1/))
#  endif
# endif
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN
!
!-----------------------------------------------------------------------
!  If appropriate, convert depth in vertical fractional coordinates to
!  actual depths in meters (negative: downwards).
!-----------------------------------------------------------------------
!
      obs_work(1:Ndatum(ng))=IniVal
!
      CALL obs_k2z (ng, 0, Lm(ng)+1, 0, Mm(ng)+1,                       &
     &              LBi, UBi, LBj, UBj, 0, N(ng),                       &
     &              rXmin(ng), rXmax(ng),                               &
     &              rYmin(ng), rYmax(ng),                               &
     &              Ndatum(ng),                                         &
     &              obs_Xgrid, obs_Ygrid, obs_depths, obs_scale,        &
     &              GRID(ng) % z0_w,                                    &
     &              obs_work)

# ifdef DISTRIBUTE
!
!  Collect all the computed depths.
!
      CALL mp_collect (ng, iNLM, Ndatum(ng), IniVal, obs_work)
# endif
!
!  Write out to depths to NetCDF file.
!
      CALL netcdf_put_fvar (ng, iNLM, DAV(ng)%name, Vname(1,idObsD),    &
     &                      obs_work, (/1/), (/Ndatum(ng)/),            &
     &                      ncid = DAV(ng)%ncid)
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN
!
!-----------------------------------------------------------------------
!  Compute model and observations comparison statistics.
!-----------------------------------------------------------------------
!
!  Initialize.
!
      DO i=1,NobsVar(ng)
        CC(i)=0.0_r8
        MB(i)=0.0_r8
        MSE(i)=0.0_r8
        SDE(i)=0.0_r8
        mod_min(i)=Large
        mod_max(i)=-Large
        mod_mean(i)=0.0_r8
        obs_min(i)=Large
        obs_max(i)=-Large
        obs_mean(i)=0.0_r8
        mod_std(i)=0.0_r8
        obs_std(i)=0.0_r8
        Ncount(i)=0
      END DO
!
!  Compute model and observations mean per each state variable.
!
      DO iobs=1,Ndatum(ng)
        IF (obs_scale(iobs).eq.1.0_r8) THEN
          i=ObsType2State(obs_type(iobs))
          Ncount(i)=Ncount(i)+1
          mod_min(i)=MIN(mod_min(i),mod_value(iobs))
          obs_min(i)=MIN(obs_min(i),obs_value(iobs))
          mod_max(i)=MAX(mod_max(i),mod_value(iobs))
          obs_max(i)=MAX(obs_max(i),obs_value(iobs))
          mod_mean(i)=mod_mean(i)+mod_value(iobs)
          obs_mean(i)=obs_mean(i)+obs_value(iobs)
        END IF
      END DO
      DO i=1,NobsVar(ng)
        IF (Ncount(i).gt.0) THEN
          mod_mean(i)=mod_mean(i)/REAL(Ncount(i),r8)
          obs_mean(i)=obs_mean(i)/REAL(Ncount(i),r8)
        END IF
      END DO
!
!  Compute standard deviation and cross-correlation between model and
!  observations (CC).
!
      DO iobs=1,Ndatum(ng)
        IF (obs_scale(iobs).eq.1.0_r8) THEN
          i=ObsType2State(obs_type(iobs))
          cff1=mod_value(iobs)-mod_mean(i)
          cff2=obs_value(iobs)-obs_mean(i)
          mod_std(i)=mod_std(i)+cff1*cff1
          obs_std(i)=obs_std(i)+cff2*cff2
          CC(i)=CC(i)+cff1*cff2
        END IF
      END DO
      DO i=1,NobsVar(ng)
        IF (Ncount(i).gt.1) THEN
          mod_std(i)=SQRT(mod_std(i)/REAL(Ncount(i)-1,r8))
          obs_std(i)=SQRT(obs_std(i)/REAL(Ncount(i)-1,r8))
          CC(i)=(CC(i)/REAL(Ncount(i),r8))/(mod_std(i)*obs_std(i))
        END IF
      END DO
!
!  Compute model bias (MB), standard deviation error (SDE), and mean
!  squared error (MSE).
!
      DO i=1,NobsVar(ng)
        IF (Ncount(i).gt.0) THEN
          MB(i)=mod_mean(i)-obs_mean(i)
          SDE(i)=mod_std(i)-obs_std(i)
          MSE(i)=MB(i)*MB(i)+                                           &
     &           SDE(i)*SDE(i)+                                         &
     &           2.0_r8*mod_std(i)*obs_std(i)*(1.0_r8-CC(i))
        END IF
      END DO
!
!  Report model and observations comparison statistics.
!

      IF (Master) THEN
        ic=0
        DO i=1,NobsVar(ng)
          svar_name(i)='           '
          text(i)='           '
          IF (Ncount(i).gt.0) THEN
            ic=ic+1
            is(ic)=i
            svar_name(ic)=TRIM(ObsName(i))
            text(ic)='-----------'
          END IF
        END DO
        WRITE (F1,10) MAX(1,ic)
        WRITE (F2,20) MAX(1,ic)
        WRITE (F3,30) MAX(1,ic)
        WRITE (F4,40) MAX(1,ic)
        WRITE (stdout,50)
        WRITE (stdout,F1) (svar_name(i),i=1,ic)
        WRITE (stdout,F2) (text(i),i=1,ic)
        WRITE (stdout,F3) 'Observation Min   ', (obs_min (is(i)),i=1,ic)
        WRITE (stdout,F3) 'Observation Max   ', (obs_max (is(i)),i=1,ic)
        WRITE (stdout,F3) 'Observation Mean  ', (obs_mean(is(i)),i=1,ic)
        WRITE (stdout,F3) 'Observation STD   ', (obs_std (is(i)),i=1,ic)
        WRITE (stdout,F3) 'Model Min         ', (mod_min (is(i)),i=1,ic)
        WRITE (stdout,F3) 'Model Max         ', (mod_max (is(i)),i=1,ic)
        WRITE (stdout,F3) 'Model Mean        ', (mod_mean(is(i)),i=1,ic)
        WRITE (stdout,F3) 'Model STD         ', (mod_std (is(i)),i=1,ic)
        WRITE (stdout,F3) 'Model Bias        ', (MB(is(i)),i=1,ic)
        WRITE (stdout,F3) 'STD Error         ', (SDE(is(i)),i=1,ic)
        WRITE (stdout,F3) 'Cross-Correlation ', (CC(is(i)),i=1,ic)
        WRITE (stdout,F3) 'Mean Squared Error', (MSE(is(i)),i=1,ic)
        WRITE (stdout,F4) 'Observation Count ', (Ncount(is(i)),i=1,ic)
      END IF
!
!  Write comparison statistics to NetCDF file.
!
      CALL netcdf_put_ivar (ng, iNLM, DAV(ng)%name, 'Nused_obs',        &
     &                      Ncount, (/1/), (/NobsVar(ng)/),             &
     &                      ncid = DAV(ng)%ncid)
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN

      CALL netcdf_put_fvar (ng, iNLM, DAV(ng)%name, 'obs_mean',         &
     &                      obs_mean, (/1/), (/NobsVar(ng)/),           &
     &                      ncid = DAV(ng)%ncid)
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN

      CALL netcdf_put_fvar (ng, iNLM, DAV(ng)%name, 'obs_std',          &
     &                      obs_std, (/1/), (/NobsVar(ng)/),            &
     &                      ncid = DAV(ng)%ncid)
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN

      CALL netcdf_put_fvar (ng, iNLM, DAV(ng)%name, 'model_mean',       &
     &                      mod_mean, (/1/), (/NobsVar(ng)/),           &
     &                      ncid = DAV(ng)%ncid)
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN

      CALL netcdf_put_fvar (ng, iNLM, DAV(ng)%name, 'model_std',        &
     &                      mod_std, (/1/), (/NobsVar(ng)/),            &
     &                      ncid = DAV(ng)%ncid)
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN

      CALL netcdf_put_fvar (ng, iNLM, DAV(ng)%name, 'model_bias',       &
     &                      MB, (/1/), (/NobsVar(ng)/),                 &
     &                      ncid = DAV(ng)%ncid)
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN

      CALL netcdf_put_fvar (ng, iNLM, DAV(ng)%name, 'SDE',              &
     &                      SDE, (/1/), (/NobsVar(ng)/),                &
     &                      ncid = DAV(ng)%ncid)
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN

      CALL netcdf_put_fvar (ng, iNLM, DAV(ng)%name, 'CC',               &
     &                      CC, (/1/), (/NobsVar(ng)/),                 &
     &                      ncid = DAV(ng)%ncid)
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN

      CALL netcdf_put_fvar (ng, iNLM, DAV(ng)%name, 'MSE',              &
     &                      MSE, (/1/), (/NobsVar(ng)/),                &
     &                      ncid = DAV(ng)%ncid)
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN
!
!  Synchronize NetCDF file to disk.
!
      CALL netcdf_sync (ng, iNLM, DAV(ng)%name, DAV(ng)%ncid)
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN
!
 10   FORMAT ('(t22,',i2,'(a11,1x))')
 20   FORMAT ('(t22,',i2,'(a11,1x),/)')
 30   FORMAT ('(a,3x,',i2,'(1p,e11.4,0p,1x))')
 40   FORMAT ('(a,3x,',i2,'(i11,1x))')
 50   FORMAT (/,' Model-Observations Comparison Statistics:',/)

      RETURN
      END SUBROUTINE stats_modobs
#else
      SUBROUTINE stats_modobs
      RETURN
      END SUBROUTINE stats_modobs
#endif

