IMPORTANT: Corrected Observation Screening and Quality Control in Primal Formulation

[arango]

• Corrected a bug in ad_misfit.F and obs_write.F in the managing of the screening and quality control variable obs_scale for the 4DVAR primal formulation. The issue is that in the primal formulation (IS4DVAR), the current time window (or survey) observations are loaded to the working arrays using local indices (array elements 1:Nobs) as opposed to global indices in the dual formulation (W4DPSAS and W4DVAR). 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.

  # ifndef WEAK_CONSTRAINT
  !
  !-----------------------------------------------------------------------
  !  Save observation reject/accept flag in GLOBAL screening variable.
  !-----------------------------------------------------------------------
  !
          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
  

  Then, in ad_misfit.F, we just need to load the correct ObsScale for the current time assimilation window:

  !
  !-----------------------------------------------------------------------
  !  Load observation reject and accept flag into screening variable.
  !-----------------------------------------------------------------------
  !
          ic=0
          DO iobs=NstrObs(ng),NendObs(ng)
            ic=ic+1
            ObsScale(ic)=ObsScaleGlobal(iobs)
          END DO
  

• Updated input 4DVAR script s4dvar.in to include the extra-observation parameters:

  ! In any application, the number of observation types is computed as:
  ! NobsVar = NstateVar + NextraObs.
  !
  ! NextraObs values are expected for keywords ExtraIndex and ExtraName.
  ! If NextraVar > 1, enter one observation type names per line and
  ! use a backslash as the continuation.
  
        NextraObs = 0
  
       ExtraIndex = 20
  
        ExtraName = radial
  
  

  Here, NextraObs is the number of extra-observation classes to consider in addition to those associated with the state variables (one-to-one correspondence). They are used in observation operators that require more that one state variable to evaluate a particular extra-observation type like HF radials, travel time, pressure, or any other complex observation operator to be implemented in the future. In any application, the number of observation types is now computed as:

       NobsVar(ng) = NstateVar(ng) + NextraObs
  

  If not processing extra-observation classes, set NextraObs=0.

The parameter ExtraIndex defines the extra-observation class identification indices, as specified in the input observation NetCDF file variable obs_type. The index has to be unique and a number greater than 7+2*NT, where NT is the total of active plus passive tracers. The 7 values include zeta, ubar, vbar, u, v, sustr, and svstr. The 2 factor for NT is to account for the surface tracer flux. Therefore, avoid using 1:7+2*NT numbers for indices because they are reserved for the state vector. Be aware that you may want to use biological data assimilation in the future. NextraObs values are expected for this keyword. This parameter is only processed when NextraObs > 0.

The keyword ExtraName defines the extra-observation class names. NextraObs values are expected. This parameter is only processed when NextraObs > 0. Enter one type name per line and use a backslash for continuation. For example:

     ExtraName = radial   \
                 pressure

Currently, we are exploring the implementation of complex observation operators. This update is a structural change to facilitate implementing such operators in the future.

• A new variable obs_meta is added to the input observations NetCDF file:

  	double obs_meta(datum) ;
  		obs_meta:long_name = "observation meta value" ;
  		obs_meta:units = "associated state variable units" ;
  

  The variable obs_meta contains additional data qualifiers for the extra-observation value specified in obs_value, which can be used in complex operators in the future.

• Introduced mapping indices from state variable to observation type (ObsState2Type) and its inverse from observation type to state variable (ObsType2State). It is done because the extra-observation index (ExtraIndex) may not be in sequential enumeration with respect to regular observation types.

In obs_write.F, the ObsState2Type is used as follows:

        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)

In obs_cost.F, the reverse mapping ObsType2State is used as follows:

       DO iobs=NstrObs(ng),NendObs(ng)
         ivar=ObsType2State(ObsType(iobs))                        !<==                                                       
         IF ((ivar.gt.0).and.(ObsScale(iobs).gt.0.0_r8).and.           &
    &        (ObsErr(iobs).ne.0.0_r8)) THEN
           cff=ObsScale(iobs)*(NLmodVal(iobs)-ObsVal(iobs))**2/        &
    &          ObsErr(iobs)
           my_ObsCost(0)=my_ObsCost(0)+cff
           my_ObsCost(ivar)=my_ObsCost(ivar)+cff
         END IF
       END DO