VERY IMPORTANT: Corrected HF radials 4D-Var Operator

[arango]

The following update corrects a parallel bug in the HF radial observation operator, corrects CPP options typo in distribute.F, and adds the variables needed to compute Desroziers et al. (2005) statistics for a 4D-Var analysis.

• The routine obs_write.F is corrected to remove a parallel bug in the processing of HF radial observations. The fractional coordinate tile range argument passed to extract_obs3d was increased by value of 0.5. It is completely incorrect and caused a parallel bug that yielded an observation screening value of obs_scale=2.0 because the datum is accounted in two different tiles:

            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)+0.5_r8, uXmax(ng)+0.5_r8,       &
       &                        uYmin(ng), uYmax(ng),                     &
       &                        time(ng), dt(ng),                         &
       &                        ObsType,  ObsVetting,                     &
       &                        Tobs, Xobs+0.5_r8, Yobs, Zobs,            &
       &                        OCEAN(ng)%u(:,:,:,NOUT),                  &
       &                        GRID(ng)%z_v,                             &
  #   ifdef MASKING
       &                        GRID(ng)%umask,                           &
  #   endif
       &                        uradial)
  
            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)+0.5_r8, vYmax(ng)+0.5_r8,       &
       &                        time(ng), dt(ng),                         &
       &                        ObsType, ObsVetting,                      &
       &                        Tobs, Xobs, Yobs+0.5_r8, Zobs,            &
       &                        OCEAN(ng)%v(:,:,:,NOUT),                  &
       &                        GRID(ng)%z_v,                             &
  #   ifdef MASKING
       &                        GRID(ng)%vmask,                           &
  #   endif
       &                        vradial)
  

  The values for uXmin, uXmax, and Xobs where wrong when computing uradial. Similarly, the values of vYmin, vYmax, and Yobs were wrong when computing vradial. There is a confusion of staggered grid location here. We need to have instead:

            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)
  
            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)
  

  The bug caused that such observation lying on the tile edge have the observation error covariance increased by a factor of two, which affected the observation impact and sensitivity computations.

Similar corrections are done to adjoint observation operator routines ad_htobs.F and ad_misfit.F.

Many thanks to Julia Levin for reporting issues about the assimilation of HF radials and Andy Moore for his help in tracking the problem.

• The module distribute.F was corrected because CPP option BOUNDARY_ALLGAHTER was miss-spelled. We need to use BOUNDARY_ALLGATHER instead. The spelling did not affect the solution because the alternated methodology was used for the collective MPI communications of input lateral boundary conditions.

Many thanks to Frank Colberg for ​reporting this bug.

---

The output ROMS 4D-Var NetCDF file (MODname) was enhanced to include the variables necessary to compute the Desroziers et al. (2005) statistics. It now includes innovation (observation minus background), increment (analysis minus background), and residual (observation minus analysis), which are used to evaluate the error covariance hypothesis of the 4D-Var analysis. The NetCDF file has the following new variables:

       double innovation(datum) ;
               innovation:long_name = "4D-Var innovation: observation minus background, d_b = y - H(x_b)" ;
               innovation:units = "state variable units" ;
               innovation:_FillValue = 1.e+37 ;
       double increment(datum) ;
               increment:long_name = "4D-Var increment: analysis minus background, dx_a = H(x_a) - H(x_b)" ;
               increment:units = "state variable units" ;
               increment:_FillValue = 1.e+37 ;
       double residual(datum) ;
               residual:long_name = "4D-Var residual: observation minus analysis, d_a = y - H(x_b + dx_a)" ;
               residual:units = "state variable units" ;
               residual:_FillValue = 1.e+37 ;

The module mod_fourdvar.F and mod_ncparam.F were modified to include the new variables.

Many thanks to John Wilkin for suggesting the new variables. We can use third-party analysis software like NOAA's ERDDAP to browse and plot these critical variables and statistics.

Desroziers, G, L. Berre, B. Chapnik, and P. Poli, 2005: Diagnosis of observation, background and analysis-error statistics in observation space, Q. J. R. Meteorol. Soc., 131, 3385-3395, doi: 10.1256/qj.05.108.

Also, the output variable BgError_value is processed and written into MODname even if BGQC is not activated:

       double BgError_value(datum) ;
               BgError_value:long_name = "Background error standard deviation at observation locations" ;
               BgError_value:units = "state variable units" ;
               BgError_value:_FillValue = 1.e+37 ;

Warning: the metadata for some of the MODname output NetCDF variables were updated for clarity in ROMS's varinfo.dat.