Updated 4D-Var algorithms

[arango]

This update includes several improvements to the 4D-Var data assimilation algorithms:

• Modified routine allocate_boundary in module mod_boundary.F to allow the allocation of the boundary arrays using both the nonlinear and adjoint ad_LBC and LBC switches. For example,

  #if defined ADJOINT || defined TANGENT || defined TL_IOMS
        IF (   LBC(iwest,isFsur,ng)%acquire.or.                           &
       &    ad_LBC(iwest,isFsur,ng)%acquire) THEN
  #else
        IF (LBC(iwest,isFsur,ng)%acquire) THEN
  #endif
          allocate ( BOUNDARY(ng) % zeta_west(LBj:UBj) )
          ...
        END IF
  

  It allows processing of the boundary data for the adjoint, tangent linear or representer models even if it is not needed by the nonlinear model when using different boundary conditions between models. Many thanks to Julia Levin for bringing this to my attention.

• Renumbered the model and message identifiers passed to routine get_state. Added new entries to differentiate the processing of various state tasks:

        character (len=48), dimension(17) :: StateMsg =                  &
      &    (/'state initial conditions,                       ',         & !01
      &      'previous state initial conditions,              ',         & !02
      &      'previous adjoint state solution,                ',         & !03
      &      'latest adjoint state solution,                  ',         & !04
      &      'surface forcing and or OBC increments,          ',         & !05
      &      'tangent linear model error forcing,             ',         & !06
      &      'impulse forcing,                                ',         & !07
      &      'v-space increments,                             ',         & !08
      &      'background state,                               ',         & !09
      &      'IC correlation standard deviation,              ',         & !10
      &      'model error correlation standard deviation,     ',         & !11
      &      'OBC correlation standard deviation,             ',         & !12
      &      'surface forcing correlation standard deviation, ',         & !13
      &      'IC normalization factors,                       ',         & !14
      &      'model error normalization factors,              ',         & !15
      &      'OBC normalization factors,                      ',         & !16
      &      'surface forcing normalization factors,          '/)          !17
  

  For Example, during the processing of the correlations standard deviations we now have:

  !  Initial conditions standard deviation. They are loaded in Tindex=1
  !  of the e_var(...,Tindex) state variables.
  !
        STDrec=1
        Tindex=1
        DO ng=1,Ngrids
          CALL get_state (ng, 10, 10, STD(1,ng)%name, STDrec, Tindex)
          IF (FoundError(exit_flag, NoError, __LINE__,                    &
       &                 __FILE__)) RETURN
        END DO
  !
  !  Model error standard deviation. They are loaded in Tindex=2
  !  of the e_var(...,Tindex) state variables.
  !
        STDrec=1
        Tindex=2
        IF (NSA.eq.2) THEN
          DO ng=1,Ngrids
            CALL get_state (ng, 11, 11, STD(2,ng)%name, STDrec, Tindex)
            IF (FoundError(exit_flag, NoError, __LINE__,                  &
       &                   __FILE__)) RETURN
          END DO
        END IF
  
  #ifdef ADJUST_BOUNDARY
  !
  !  Open boundary conditions standard deviation.
  !
        STDrec=1
        Tindex=1
        DO ng=1,Ngrids
          CALL get_state (ng, 12, 12, STD(3,ng)%name, STDrec, Tindex)
          IF (FoundError(exit_flag, NoError, __LINE__,                    &
       &                 __FILE__)) RETURN
        END DO
  #endif
  #if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
  !
  !  Surface forcing standard deviation.
  !
        STDrec=1
        Tindex=1
        DO ng=1,Ngrids
          CALL get_state (ng, 13, 13, STD(4,ng)%name, STDrec, Tindex)
          IF (FoundError(exit_flag, NoError, __LINE__,                    &
       &                 __FILE__)) RETURN
        END DO
  #endif
  

  Similarly, to read and process the error covariance normalization factors we have:

            NRMrec=1
            CALL get_state (ng, 14, 14, NRM(1,ng)%name, NRMrec, 1)
            IF (FoundError(exit_flag, NoError, __LINE__,                  &
       &                   __FILE__)) RETURN
  
            IF (NSA.eq.2) THEN
              CALL get_state (ng, 15, 15, NRM(2,ng)%name, NRMrec, 2)
              IF (FoundError(exit_flag, NoError, __LINE__,                &
       &                     __FILE__)) RETURN
            END IF
  
  #ifdef ADJUST_BOUNDARY
            CALL get_state (ng, 16, 16, NRM(3,ng)%name, NRMrec, 1)
            IF (FoundError(exit_flag, NoError, __LINE__,                  &
       &                   __FILE__)) RETURN
  #endif
  
  #if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
            CALL get_state (ng, 17, 17, NRM(4,ng)%name, NRMrec, 1)
            IF (FoundError(exit_flag, NoError, __LINE__,                  &
       &                   __FILE__)) RETURN
  #endif
          END IF
  

  WARNING: All the 4D-Var drivers and associated algorithms were modified because of the renumbering of the tasks. The new indices allow to turn off the time messages when processing both standard deviation and normalization factors. It was confusing for some users. In get_state.F, we just need to overwrite the string t_code as follows:

       IF (Master) THEN
         CALL time_string (INPtime, t_code)
         IF ((10.le.model).and.(model.le.17)) THEN
           t_code=' '              ! time is meaningless for these fields
         END IF
         WRITE (Tstring,'(f15.4)') tdays(ng)
         IF (ERend.gt.ERstr) THEN
           WRITE (stdout,40) string, 'Reading '//TRIM(StateMsg(msg)),    &
      &                      t_code, ng, Nrun, TRIM(ADJUSTL(Tstring)),   &
      &                      ncname(lstr:lend), InpRec, Tindex
         ELSE
           WRITE (stdout,50) string, 'Reading '//TRIM(StateMsg(msg)),    &
      &                      t_code, ng, TRIM(ADJUSTL(Tstring)),         &
      &                      ncname(lstr:lend), InpRec, Tindex
         END IF
       END IF
  

• Updated the observation impact and sensitivity algorithms to allow the saving of the analysis data for each inner loop. A new C-preprocessing option IMPACT_INNER is introduced to achieve the saving of analysis quantities into the output MOD NetCDF file. For Example, in def_mod.F we have:

  !
  !  Define observation impact due to initial condition increments.
  !
  #  ifdef IMPACT_INNER
          Vinfo( 1)='ObsImpact_IC'
          Vinfo( 2)='observation impact due to initial conditions'
          vardim(1)=datumDim
          vardim(2)=NinnerDim
          Vinfo(24)='_FillValue'
          Aval(6)=spval
          status=def_var(ng, iNLM, DAV(ng)%ncid, varid, nf90_double,      &
       &                 2, vardim, Aval, Vinfo, ncname)
          IF (FoundError(exit_flag, NoError, __LINE__,                    &
       &                 __FILE__)) RETURN
  #  else
          Vinfo( 1)='ObsImpact_IC'
          Vinfo( 2)='observation impact due to initial conditions'
          vardim(1)=datumDim
          Vinfo(24)='_FillValue'
          Aval(6)=spval
          status=def_var(ng, iNLM, DAV(ng)%ncid, varid, nf90_double,      &
       &                 1, vardim, Aval, Vinfo, ncname)
          IF (FoundError(exit_flag, NoError, __LINE__,                    &
       &                 __FILE__)) RETURN
  #  endif
  

  Notice that the output variable ObsImpact_IC can be either one- or two-dimension array.

• Added a new adjoint routine ad_wvelocity.F to allow observation impacts or sensitivities using vertical velocity functionals. Many thanks to Andy Moore for coding the adjoint of ROMS diagnostic routine wvelocity.F. Also, for modifying the algorithms to allow impacts and sensitivities functionals that use vertical velocity.

• Added new switch Lstate(isWvel) in standard input script ocean.in to allow the processing of functional using vertical velocity:

  ! Logical switches (TRUE/FALSE) to specify the adjoint state variables
  ! whose sensitivity is required.
  
  Lstate(isFsur) == F                        ! free-surface
  Lstate(isUbar) == F                        ! 2D U-momentum
  Lstate(isVbar) == F                        ! 2D V-momentum
  Lstate(isUvel) == T                        ! 3D U-momentum
  Lstate(isVvel) == T                        ! 3D V-momentum
  Lstate(isWvel) == T                        ! 3D W-momentum
  
  Lstate(isTvar) == F F                      ! NT tracers
  

  WARNING: All ROMS input script were modified to add Lstate(isWvel).