IMPORTANT: Split 4D-Var algorithms, Phase III

[arango]

The following is a very critical update of ROMS 4D-Var algorithms. Several permanent changes are done on how the surface forcing fields for all the adjoint-based algorithms are specified. Please read the following information carefully.

• WARNING: Due to the complexities of the split 4D-Var, coupled 4dvar, and regular 4D-Var algorithms and their associated analysis tools, the processing logic of how the surface atmospheric forcing, in general, is changed in routines get_data.F and set_data.F and respective TLM, RPM, and ADM versions.

We have incredible convoluted CPP options with BULK_FLUXES on and off, direct atmosphere coupling (FRC_COUPLING) on and off, and forward trajectory forcing of the TLM, RPM, and ADM kernels. Also, we may desire to impose or not the NLM prior surface forcing in subsequent outer-loops (if Nouter>1) and in the analysis 4D-Var phase.

The state variable for the surface (stflx) and bottom (btlfx) tracer fluxes are no longer used for I/O processing. New variables (stflux) and (btflux) are introduced for such tasks. Therefore, we now have:

!  Surface tracer fluxes.                                              !
!                                                                      !
!  stflux       Forcing surface flux of tracer type variables from     !
!                 data, coupling, bulk flux parameterization, or       !
!                 analytical formulas.                                 !
!                                                                      !
!                 stflux(:,:,itemp)  surface net heat flux             !
!                 stflux(:,:,isalt)  surface net freshwater flux (E-P) !
!                                                                      !
!  stfluxG      Latest two-time snapshots of input "stflux" grided     !
!                 data used for interpolation.                         !
!                                                                      !
!  stflx        ROMS state surface flux of tracer type variables       !
!                 (TracerUnits m/s) at horizontal RHO-points, as used  !
!                 in the governing equations.                          !                                                                                                !
!                                                                      !
!  Bottom tracer fluxes.                                               !
!                                                                      !
!  btflux       Forcing bottom flux of tracer type variables from      !
!                 data or analytical formulas. Usually, the bottom     !
!                 flux of tracer is zero.                              !
!                                                                      !
!                 btflux(:,:,itemp)  bottom heat flux                  !
!                 btflux(:,:,isalt)  bottom freshwater flux            !
!                                                                      !
!  btfluxG      Latest two-time snapshots of input "vtflux" grided     !
!                 data used for interpolation.                         !
!                                                                      !
!  btflx        ROMS state bottom flux of tracer type variables        !
!                 (TracerUnits m/s) at horizontal RHO-points, as used  !
!                 in the governing equations.                          !

Recall that stflx and btflx are the surface and bottom boundary conditions for the vertical tracer diffusion term. The new strategy removes the issues we were having when manipulating such fluxes in set_vbc.F, like surface heat and freshwater corrections. In particular, the scaling of freshwater flux (E-P)/ρ, which needs to be multiplied by the surface salinity, was problematic. This issue was addressed earlier in ​src:ticket:806, but it is further reworked here.

This change simplified the logic in all the get_data and set_data routines. It is a generic solution to all combinations of options and algorithms. For example, in tl_set_data we have:

#  ifdef SALINITY
#   ifdef ANA_SSFLUX
!
!  Surface freshwater (E-P) flux (m/s) from analytical function.
!
      CALL ana_stflux (ng, tile, iNLM, isalt)
#   else

#    if !(defined BULK_FLUXES  || defined EMINUSP      || \
          defined FRC_COUPLING || defined SRELAXATION)
!
!  Surface freshwater (E-P) flux (m/s) from NetCDF variable "swflux".
!
      CALL set_2dfld_tile (ng, tile, iNLM, idsfwf,                      &
     &                     LBi, UBi, LBj, UBj,                          &
     &                     FORCES(ng)%stfluxG(:,:,:,isalt),             &
     &                     FORCES(ng)%stflux (:,:,isalt),               &
     &                     update)
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN

#    elif defined BULK_FLUXES  && !defined EMINUSP
!
!  Surface freshwater (E-P) flux (m/s) from NetCDF variable "EminusP".
!
      IF (Lprocess) THEN
        CALL set_2dfld_tile (ng, tile, iNLM, idEmPf,                    &
     &                       LBi, UBi, LBj, UBj,                        &
     &                       FORCES(ng)%stfluxG(:,:,:,isalt),           &
     &                       FORCES(ng)%stflux (:,:,isalt),             &
     &                       update)
        IF (FoundError(exit_flag, NoError, __LINE__,                    &
     &                 __FILE__)) RETURN
      END IF
#    endif     
#   endif
      ...
#  endif

Again, the important difference is that data processing is done with variables stflux and stfluxG and not with state variable stflx.

• WARNING: The CPP option NL_BULK_FLUXES is deprecated and replaced with FORWARD_FLUXES and PRIOR_BULK_FLUXES.

• The option FORWARD_FLUXES indicates that the surface forcing fields for the TLM, RPM, and ADM kernels is from the background nonlinear model trajectory (BLK file structure) instead of input atmosphere forcing files (FRC file structure). Such fields are now stored in the quicksave (QCK) output NetCDF instead of history (FWD) files. It allows saving the surface forcing frequently (say every 1, 2, or 3 hours) without increasing the history file size. Currently, we have the following statement in mosty all the adjoint-based drivers:

  # ifdef FORWARD_FLUXES
  !
  !  Set the BLK structure to contain the nonlinear model surface fluxes
  !  needed by the tangent linear and adjoint models. Also, set switches
  !  to process that structure in routine "check_multifile". Notice that
  !  it is possible to split the solution into multiple NetCDF files to
  !  reduce their size.
  !
  !  The switch LreadFRC is deactivated because all the atmospheric
  !  forcing, including shortwave radiation, is read from the NLM
  !  surface fluxes or is assigned during ESM coupling.  Such fluxes
  !  are available from the QCK structure. There is no need for reading
  !  and processing from the FRC structure input forcing-files.
  !
        CALL edit_multifile ('QCK2BLK')
        IF (FoundError(exit_flag, NoError, __LINE__,                      &
       &               __FILE__)) RETURN
        DO ng=1,Ngrids
          LreadBLK(ng)=.TRUE.
          LreadFRC(ng)=.FALSE.
          LreadQCK(ng)=.FALSE.
        END DO
  # endif
  

  The User needs to provide a reasonable value for NQCK in roms.in to correctly sample the daily cycle. Please avoid activating the Qout switches for lots of fields written in the quicksave file, so it is relatively small. Avoid activating switches for 3D variables. Optimally, we need to activate the following switches:

  Qout(idsurT) == T T     ! temp_sur, salt_sur surface temperature and salinity
  
  Qout(idUsms) == T       ! sustr              surface U-stress
  Qout(idVsms) == T       ! svstr              surface V-stress
  
  Qout(idPair) == T       ! Pair               surface air pressure
  Qout(idTair) == T       ! Tair               surface air temperature
  Qout(idUair) == T       ! Uair               surface U-wind component
  Qout(idVair) == T       ! Vair               surface V-wind component
  
  Qout(idTsur) == T T     ! shflux, ssflux     surface net heat and salt flux
  Qout(idLhea) == T       ! latent             latent heat flux
  Qout(idShea) == T       ! sensible           sensible heat flux
  Qout(idLrad) == T       ! lwrad              longwave radiation flux
  Qout(idSrad) == T       ! swrad              shortwave radiation flux
  Qout(idEmPf) == T       ! EminusP            E-P flux
  Qout(idevap) == T       ! evaporation        evaporation rate
  Qout(idrain) == T       ! rain               precipitation rate
  

  Some of these fields depend on the activated CPP options. The appropriate fields are written to the quicksave NetCDF file.

Previously, the time snapshots of the background trajectory used to linearize the TLM and ADM were incompatible with the direct forcing from atmospheric input files if SOLAR_SOURCE and DIURNAL_SRFLUX are activated. The shortwave contribution to the net surface heat flux is not the same. The new strategy removes the inconsistency.

• The option PRIOR_BULK_FLUXES is used to impose the initial background trajectory (prior) surface forcing fields stored in QCK in subsequent outer-loops (Nouter>1) and the analysis phase. In main3d.F, the call to bulk_flux is avoided in such cases:

              DO ig=1,GridsInLayer(nl)
                ng=GridNumber(ig,nl)
                DO tile=first_tile(ng),last_tile(ng),+1
  # ifdef BULK_FLUXES
  #  if defined FOUR_DVAR && defined PRIOR_BULK_FLUXES
                  IF (Nrun.eq.1) CALL bulk_flux (ng, tile)
  #  else
                  CALL bulk_flux (ng, tile)
  #  endif
  # endif
                  ...
                END DO
              END DO
  

  Notice that get_data.F will read the required surface forcing fields from the QCK NetCDF file.

• The observation operator was updated to include new variables in the output MOD NetCDF file:

         double NLmodel_unvetted(datum) ;
                 NLmodel_unvetted:long_name = "initial nonlinear model at observation locations unvetted by quality control" ;
                 NLmodel_unvetted:units = "state variable units" ;
  
         double NLmodel_value_unvetted(Nouter, datum) ;
                 NLmodel_value_unvetted:long_name = "nonlinear model at observation locations unvetted by quality control per outer-loop" ;
                 NLmodel_value_unvetted:units = "state variable units" ;
  

  The values rejected during background quality control (BGQC) are saved before multiplying by obs_scale screening. It facilitates analyzing and comparing the observation vector. If BGQC is activated, NLmodel_unvetted will have values for all the observations bounded in space and time for the 4D-Var cycle, whereas NLmodel_initial will have only the values after the quality control screening.

All the _FillValue attributes for output variables in the MOD NetCDF file were removed. Many thanks to John Wilkin for his suggestions to improve the output MOD NetCDF file.

• Added option REGRID_SHPIRO to apply a Shapiro filter in routine regrid.F called by nf_fread2d to remove spikes and noise due to the bilinear interpolation of very coarse 2D input data.

---

The primary drivers to the split 4D-Var algorithms are released:

• split_i4dvar_ocean.h: Split I4D-Var algorithm. It uses submit_split_i4dvar.sh to run the driver.

• split_r4dvar_ocean.h: Split R4D-Var algorithm. It uses submit_split_r4dvar.sh to run the driver.

• split_rbl4dvar_ocean.h: Split RBL4D-Var algorithm. It uses submit_split_rbl4dvar.sh to run the driver.

The bash scripts in ROMS/Bin are complex and well documented. Several test cases were added to the test repository for the WC13 application. It took me several months to code and to test these new drivers. We will add documentation to wikiROMS about the new 4D-Var drivers soon.

More complexity is being added and tested for these drivers like coarser grid resolution in the increment phase during the inner-loops. The increment phase may be run at a lower precision.

I am currently testing the coupled 4D-Var that allows the background phase to be part of a coupling system. The coupling between ROMS and the atmosphere model is done with the ESMF/NUOPC library. Such capabilities will be released in the future.

---

WARNING: The metadata file varinfo.dat was modified to include additional variables.