Opened 3 years ago

Closed 3 years ago

#757 closed upgrade (Done)

Updated 4D-Var algorithms

Reported by: arango Owned by:
Priority: major Milestone: Release ROMS/TOMS 3.7
Component: Adjoint Version: 3.7
Keywords: Cc:

Description

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

Change History (1)

comment:1 Changed 3 years ago by arango

  • Resolution set to Done
  • Status changed from new to closed
Note: See TracTickets for help on using tickets.