Updated 4D-Var observation impact and sensitivity

[arango]

In the past, the observation impact/sensitivity in the 4D-Var dual formulation algorithms (W4D-PSAS, W4D-VAR) was restricted for a single outer loop (Nouter = 1) because the minimization solvers are very tricky for multiple outer loops. In ​src:ticket:702 4D-Var major update, we added the new Restricted B-preconditioned Lanczos (RBLanczos) solver by Gurol et al. (2014). The RBLanczos solver (rpcg_lanczos.F) included the augmented terms so we can run over multiple outer loops.

Still, the observations impact/sensitivity requires particular care in multiple outer loops applications (say, Nouter = 2). We require running separately each outer loop to compute the desired observation impact or sensitivity.

For example, we need to run first the regular W4D-PSAS assimilation cycle and save all the output files. Then, run the observation impact/sensitivity by activating W4DPSAS_SENSITIVITY and SKIP_NLM. Notice that the nonlinear model (NLM) forward trajectory is not computed (since SKIP_NLM is activatted) but read from the regular data assimilation run. That is, we read files *_fwd_0001.nc, *_fwd_0001.nc, and so on. We need to specify wich outer loop is beeing processed for the observation impact/sensitivity by setting Nimpact in s4dvar.in to its sppropriate outer loop value.

This update simplifies the running of observation impact/sensitivities in multiple outer loop applications:

• Added new parameter Nimpact in 4D-Var input script s4dvar.in:

  ! If observations impact or observations sensitivity, set the 4D-Var
  ! outer loop to process. It must be less or equal that Nouter.
  
           Nimpact = 1
  

It specifies the outer loop to consider in the computation of the observations impact or observation sensitivity, Nimpact =< Nouter. The observation analysis needs to be computed separately for each outer loop. The full analysis for all outer loops is combined offline.

• The computation of the observation screening and quality control flag (obs_scale NetCDF output variable and ObsScale ROMS variable) was redesigned so it is computed and written just once to facilitate the background quality control (BGQC) option. A new working variable ObsVetting is introduced and passed to the extraction routines (extract_obs*d and ad_extract_obs*d) to avoid the continuous rewriting of the important ObsScale variable that is used elsewhere. In obs_write.F, we now have:

  !
  !-----------------------------------------------------------------------
  !  If processing nonlinear trajectory, load observation reject/accept
  !  flag into screening variable.  This needs to be done once and
  !  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
  
  

  A new logical switch wrtObsScale(:) is introduced to control when to process the ObsScale variable in ROMS. This is done during the first pass to the nonlinear model when computing the background trajectory and interpolating the solution to the observation locations.

• If the background quality control (BGQC) is activated, it is done after loading the screening values in ObsScale and also computed once in obs_write.F:

  # if defined BGQC && !defined IS4DVAR_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 IS4DVAR
                df1=(ObsVal(iobs)-NLmodVal(iobs))/BgErr(iobs)
                df2=(1.0_r8/ObsErr(iobs))/(BgErr(iobs)*BgErr(iobs))
  #  elif defined WEAK_CONSTRAINT
  #    ifdef W4DVAR
                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
  

  This fixes the problem of running observations impact/sensitivity when the BGQC is activated over multiple outer loops. Before, we were getting zero impacts because the continuous internal manipulation of the ObsScale variable when running the adjoint and tangent linear models.

• Added reading the obs_scale variable from input NetCDF in the observation impact/sensitivity drivers:

  #ifdef SKIP_NLM
  !
  !-----------------------------------------------------------------------
  !  If skiping runing nonlinear model, read in observation screening and
  !  quality control flag.
  !-----------------------------------------------------------------------
  !
        wrtObsScale(1:Ngrids)=.FALSE.
        DO ng=1,Ngrids
          CALL netcdf_get_fvar (ng, iTLM, LCZ(ng)%name, Vname(1,idObsS),  &
       &                        ObsScale)
          IF (exit_flag.ne.NoError) RETURN
        END DO
  #endif
  

  This allows computing observation impact/sensitivity for each outer loop separately since we cannot recompute it by running the nonlinear model. Recall that we need SKIP_NLM activated for the multiple outer loop applications.

---

WARNING:

The computation of the 4D-Var observations impact or sensitivity are complex and requires expertise in setting the impact/sensitivity functional, running the application, and post-processing the results.

[arango]

Fixed a missing CPP option in obs_write.F:

        IF (wrtNLmod(ng)) THEN
          DO iobs=Mstr,Mend
            NLmodVal(iobs)=ObsScale(iobs)*NLmodVal(iobs)
          END DO
        END IF

# ifdef TLM_OBS
        IF (wrtTLmod(ng).or.wrtRPmod(ng)) THEN
          DO iobs=Mstr,Mend
#  ifdef IS4DVAR_SENSITIVITY
            TLmodVal(iobs)=ObsScale(iobs)*TLmodVal(iobs)*ObsErr(iobs)
#  else
            TLmodVal(iobs)=ObsScale(iobs)*TLmodVal(iobs)
#  endif
          END DO
        END IF
# endif