Very IMPORTANT: Added new 4DVAR operator for HF radials

[arango]

This update is not public and the only users having access to our omlab svn repository may use it. This is an ongoing funded research. This is my private research svn branch and very few have access to it.

---

The 4DVar data assimilation data operators were updated to allow extra-observations type classes with multiple correspondences with the model state variables. Previously, the observation operator in obs_write.F (NLM, TLM, RPM), ad_misfit.F (ADM, primal formulation), and ad_htobs.F (ADM, dual formulation) have one-to-one correspondence with each state variable (zeta, u, v, T, S, and other passive tracers). Recall that the observation operator spatially interpolates the model at the observation locations at the appropiate time. Examples of extra-observation classes include HF radials, travel time, pressure, etc. This updade includes the observation operator to process HF radar surface velocity meassurements (radials) that affect both u- and v-momentum state variables.

• The input 4DVAR script s4dvar.in was modified 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 = 1
  
       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, etc. 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, however, only the radials operator is coded.

• 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. For example, it contains the velocity radials angle (radians) from HF Radar data. The values in obs_meta for any non-radial observation are ignored.

In curvilinear coordinates, the velocity radial is defined as:

    radial = u * COS(obs_meta - angler) + v * SIN(obs_meta - angler)

where angler is the curviliar rotation angle. If no curvlinear grid, angler=0.

• 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 routine initial_fourdvar of mod_fourdvar.F, these mapping indices are allocated and initialized as follows:

!
!  Allocate and initialize observation types names and indices.
!  Notice that a mapping from state-to-type (ObsState2Type) and its
!  inverse type-to-state (ObsType2State) indices are needed because
!  the User is allowed to add extra-observation operators with
!  nonsequential type enumeration.
!
!  Both mapping arrays ObsState2Type and ObsType2State have a zero
!  array element to allow applications with no extra-observations
!  to work with their zero associated state index (as initialized
!  in mod_ncparam.F).  For example, if the index "isRadial" is not
!  redefined below, the following assignment
!
!        ObsState2Type(isRadial)=ObsState2Type(0)=0
!        ObsType2State(isRadial)=ObsType2State(0)=0
!
!  is still legal with isRadial=0. It avoids a Fortran segmentation
!  violation (i.e., subscript #1 of the array ObsState2Type has
!  value 0 which is less than the lower bound of 1).  Sorry for
!  the awkward logic but we need a generic way to specify extra-
!  observation operators.
!
      allocate ( ObsName(MAXVAL(NobsVar)) )
      icount=MAXVAL(NstateVar)
      ObsState2Type(0)=0
      DO i=1,icount                               ! 5+NT
        ObsState2Type(i)=i
        ObsName(i)=TRIM(Vname(1,idSvar(i)))
      END DO
      IF (NextraObs.gt.0) THEN
        DO i=1,NextraObs
          icount=icount+1
          ObsName(icount)=TRIM(ExtraName(i))
          SELECT CASE (TRIM(uppercase(ExtraName(i))))
            CASE ('RADIAL')
              isRadial=icount
              ObsState2Type(icount)=ExtraIndex(i)
          END SELECT
        END DO
      END IF

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

          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)

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

• Modified routined obs_write.F, ad_misfit.F, and ad_htob.F to include the HF Radar data operator:

  !
  !  Compute NLM radial velocity.
  !
            DO iobs=Mstr,Mend
              IF (ObsType(iobs).eq.ObsState2Type(isRadial)) THEN
  #   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
  #   ifdef BGQC
                BgErr(iobs)=MAX(uBgErr(iobs), vBgErr(iobs))
  #   endif
              END IF
            END DO
  

  where uradial and vradial are the model velocities interpolated to the radial velocity locations in space and time, obs_angler is the curvilinear grid angle interpolated at the radial locations, and obs_meta is the HF radar radial angle.

The adjoint forcing is simply:

        IF (FOURDVAR(ng)%ObsCount(isRadial).gt.0) THEN
          DO iobs=1,Nobs(ng)
            ad_uradial(iobs)=IniVal
            ad_vradial(iobs)=IniVal
          END DO
          DO iobs=1,Nobs(ng)
            IF (ObsType(iobs).eq.ObsState2Type(isRadial)) THEN
#  ifdef CURVGRID
              angle=ObsMeta(iobs)-ObsAngler(iobs)
              ad_uradial(iobs)=ad_uradial(iobs)+                        &
     &                         ADmodVal(iobs)*COS(angle)
              ad_vradial(iobs)=ad_vradial(iobs)+                        &
     &                         ADmodVal(iobs)*SIN(angle)
#  else
              ad_uradial(iobs)=ad_uradial(iobs)+                        &
     &                         ADmodVal(iobs)*COS(ObsMeta(iobs))
              ad_vradial(iobs)=ad_vradial(iobs)+                        &
     &                         ADmodVal(iobs)*SIN(ObsMeta(iobs))
#  endif
            END IF
          END DO
          ...
        END IF

• The HF radial assimilation was tested with Brian Powell Hawaii application. Both IS4DVAR and W4DPSAS gives identical cost functions.

---

Many thanks to Brain Powell for his help in the coding and testing of the HF Radar operator to assimilate surface velocity radials.

---

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