VERY IMPORTANT: Added multi-file capability to adjoint-based drivers

[arango]

There is a lot of information in this update, please read carefully and understand if you use any of the algorithms mentioned below.

---

This update is crucial because it includes the capability allows the forward nonlinear solution to be time split into multiple files to avoid creating huge NetCDF files in applications with large grids. Recall that the nonlinear trajectory is used to linearize the tangent linear and adjoint models used in 4D-Var data assimilation, observation impacts, observation sensitivities, adjoint sensitivities, generalized stability analysis, and other adjoint-based algorithms.

• Two subroutines (edit_multifile and edit_file_struct) are included in edit_multifile.F to manipulate the basic state forward trajectory used in various adjoint-based drivers. The TYPE(T_IO) structure can be deallocated and allocate with multiple file capabilities if so desired. For example in edit_multifile, we have:

        DO ng=1,Ngrids
  !
          SELECT CASE (TRIM(ADJUSTL(task)))
  !
  !  Load HIS information into the FWD structure so it can be used to
  !  process the NLM background trajectory by the ADM and TLM kernels.
  !
            CASE ('HIS2FWD')
              IF (ndefHIS(ng).gt.0) THEN
                IF (HIS(ng)%ncid.ne.-1) THEN
                  CALL netcdf_close (ng, iNLM, HIS(ng)%ncid)
                END IF
                Nfiles=ntimes(ng)/ndefHIS(ng)
                IF (nHIS(ng).eq.ndefHIS(ng)) Nfiles=Nfiles+1
                CALL edit_file_struct (ng, Nfiles, FWD)
                DO ifile=1,Nfiles
                  FWD(ng)%files(ifile)=TRIM(HIS(ng)%files(ifile))
                END DO
                FWD(ng)%name=TRIM(FWD(ng)%files(1))
              ELSE
                FWD(ng)%ncid=HIS(ng)%ncid
                FWD(ng)%name=TRIM(HIS(ng)%name)
                FWD(ng)%files(1)=TRIM(HIS(ng)%name)
              END IF
  
            CASE ('HIS2BLK')
              ...
  
            CASE ('FWD2BLK')
              ...
  
            CASE ('FWD2HIS')
              ...
  
            CASE ('FCTA2FWD')
              ...
  
            CASE ('FCTA2BLK')
              ...
  
            CASE ('FCTB2FWD')
              ...
  
            CASE ('FCTB2BLK')
              ...
  
            CASE ('TLM2FWD')
              ...            
  
          END SELECT
        END DO
  

As you can see, there are several possibilities in the manipulations between file structures. In the calling driver, we need to add the specific instructions for such handling. For example in w4dpsas_ocean.h, we have:

!
!  Set structure for the nonlinear forward trajectory to be processed
!  by the tangent linear and adjoint models. Also, set switches to
!  process the FWD structure in routine "check_multifile". Notice that
!  it is possible to split solution into multiple NetCDF files to reduce
!  their size.
!
      CALL edit_multifile ('HIS2FWD')
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN
      DO ng=1,Ngrids
        LreadFWD(ng)=.TRUE.
      END DO

#if defined BULK_FLUXES && defined NL_BULK_FLUXES
!
!  Set structure for the nonlinear surface fluxes to be processed by
!  by the tangent linear and adjoint models. Also, set switches to
!  process the BLK structure in routine "check_multifile".  Notice that
!  it is possible to split solution into multiple NetCDF files to reduce
!  their size.
!
      CALL edit_multifile ('HIS2BLK')
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN
      DO ng=1,Ngrids
        LreadBLK(ng)=.TRUE.
      END DO
#endif

• The routine check_multifile was updated to include more TYPE(T_IO) structures to process for information about multiple files or not. Two new generic subroutines multifile_info_s1d and multifile_info_s2d were written to inquire about 1D and 2D TYPE(T_IO) structures, repectively. They are called in a compact way in check_multifile:

  #ifdef FORWARD_READ
  !
  !-----------------------------------------------------------------------
  !  Nonlinear model forward trajectory data.
  !-----------------------------------------------------------------------
  !
        IF (LreadFWD(ng)) THEN
          CALL multifile_info_s1d (ng, model, 'Forward', FWD)
          IF (FoundError(exit_flag, NoError, __LINE__,                    &
       &                 __FILE__)) RETURN
        END IF
  #endif
  
        ...
  
  #ifdef FRC_FILE
  !
  !-----------------------------------------------------------------------
  !  Input forcing data.
  !-----------------------------------------------------------------------
  !
        file_type='Forcing'
        max_files=MAXVAL(nFfiles)
        CALL multifile_info_s2d (ng, model, file_type, nFfiles,           &
       &                         max_files, FRC)
        IF (FoundError(exit_flag, NoError, __LINE__,                      &
       &               __FILE__)) RETURN
  #endif
  

• Several switches (LreadADM, LreadBLK, LreadFWD, and LreadTLM) are introduced in mod_scalar.F to activate the processing of information about file structures at the right time in the algorithms.

• The TYPE(T_IO) file structure definition in mod_param.F now includes the head field, which remains unchanged after initialization and has the I/O heading prefix filename. It is used as follows:

         WRITE (HIS(ng)%name,10) TRIM(FWD(ng)%head), outer
         lstr=LEN_TRIM(HIS(ng)%name)
         HIS(ng)%base=HIS(ng)%name(1:lstr-3)
  

  The HIS(ng)%base is used in output.F to create the multile files. For example, in 4D-Var we can have:

    wc13_fwd_outer0_0001.nc
    wc13_fwd_outer0_0002.nc
    wc13_fwd_outer0_0003.nc
    wc13_fwd_outer0_0004.nc
  
    wc13_fwd_outer1_0001.nc
    wc13_fwd_outer1_0002.nc
    wc13_fwd_outer1_0003.nc
    wc13_fwd_outer1_0004.nc
  

  Here, the basic state forward trajectory is split in daily files having data records every two-hours in a 4D-Var 4-days data assimilation cycle for each outer loop.

WARNING: If not using output history multi-files, NDEFHIS=0 in roms.in, we will have

  wc13_fwd_outer0.nc   instead of   wc13_fwd_0000.nc
  wc13_fwd_outer1.nc   instead of   wc13_fwd_0001.nc

as in previous versions of the code. So you need to update your post-processing scripts.

• Several files were updated to prevent having Too many open files error. In the UNIX environment, there is a limit to the number of open files during program execution. Use the commands:

   % ulimit -a
   % ulimit -S -n
  

  to check the limit in your system. Usually, 256 files can be opened by default. A new C-preprocessing option (CHECK_OPEN_FILES) is introduced to report the number of files created, opened, and closed for an application. The report is written to Fortran file fort.1000. One can easily check that report by typing:

   % grep CREATE fort.1000 | wc -l
   % grep OPEN   fort.1000 | wc -l
   % grep CLOSE  fort.1000 | wc -l
  

  or using provided scripts check_openfiles.bash or check_openfiles.sh located in ROMS/Bin.

• Use the multiple filename entries in roms.in and end with the vertical bar (|) symbol. For example:

      FWDNAME == wc13_fwd_outer0_0001.nc |
                 wc13_fwd_outer0_0002.nc |
                 wc13_fwd_outer0_0003.nc |
                 wc13_fwd_outer0_0004.nc
  
      FWDNAME == wc13_his_0001.nc |
                 wc13_his_0002.nc |
                 wc13_his_0003.nc |
                 wc13_his_0004.nc
  
     FCTnameA == wc13_fct_A_0001.nc |
                 wc13_fct_A_0002.nc |
                 wc13_fct_A_0003.nc |
                 wc13_fct_A_0004.nc |
                 wc13_fct_A_0005.nc |
                 wc13_fct_A_0006.nc |
                 wc13_fct_A_0007.nc
  
     FCTnameB == wc13_fct_B_0001.nc |
                 wc13_fct_B_0002.nc |
                 wc13_fct_B_0003.nc |
                 wc13_fct_B_0004.nc |
                 wc13_fct_B_0005.nc |
                 wc13_fct_B_0006.nc |
                 wc13_fct_B_0007.nc
  

  if appropriate. The files need to be ordered in ascending time coordinate.

• The multipe output files are activated in roms.in by setting any of their parameters (NDEFHIS, NDEFQCK, NDEFAVG, NDEFDIA) to be greater than zero. ROMS creates new files every specified number of timesteps. The parameters NDEFADJ and NDEFTLM are special and some algorithms do not allow its values to be greater than zero.

• Started to solve the issue with DSTART in the adjoint-based algorithms. The management of DSTART has evolved over the years. It is set in roms.in standard input file. Nowadays, we need to set its value to the actual initialization time (in days) every time that we run any of the 4D-Var drivers or their associated analysis tools. The reason for it is that tl_initial, rp_initial, and ad_initial rarely reads an initial conditions NetCDF file for the TLM, RPM, and ADM, respectively. Usually, those initial conditions are initialized from rest (zeros) or computed by the 4D-Var driver. Therefore, we need to use DSTART to compute the correct and crucial time-stepping parameters: time(ng), tdays(ng), ntstart(ng), ntend(ng), and ntfirst(ng) since we are not reading an initial NetCDF file to set time(ng) for those kernels.

My current strategy is to use the initialization time from the nonlinear model or forward trajectory, which is now saved in INItime(ng). It is already coded but not activated because it requires rigorous testing. I am not fully satisfied yet. This capability is activated with the new GENERIC_DSTART C-preprocessing option. Please do not activate this option! You may get different and wrong results. There is a tricky bug somewhere that I haven't been able to find. In tl_initial.F, we have:

# ifdef GENERIC_DSTART
!
!  Rarely, the tangent linear model is initialized from a NetCDF file,
!  so we do not know its actual initialization time. Usually, it is
!  computed from DSTART, implying that its value is correct in the ROMS
!  input script. Therefore, the user needs to check and update its value
!  to every time that ROMS is executed. Alternatively, if available, we
!  can use the initialization time from the nonlinear model, INItime.
!  This variable is assigned when computing or processing the basic
!  state trajectory needed to linearize the adjoint model.
!
      IF (INItime(ng).lt.0.0_dp) THEN
        my_dstart=dstart                    ! ROMS input script
      ELSE
        my_dstart=INItime(ng)/86400.0_dp    ! NLM IC time is known
      END IF
      tdays(ng)=my_dstart
# else
      tdays(ng)=dstart
# endif
      time(ng)=tdays(ng)*day2sec

The reason why I am looking at it is that we want an unrestrained value of DSTART because of the multiple-file option that now is allowed in the adjoint-based algorithms. Recall that DSTART controls the split file enumeration and the starting value for iic(ng). If DSTART is set to the initial conditions time value, the count always starts at zero.

• Corrected few tiny bugs in the 2D applications set-up. Many thanks to Brian Powell for bringing these issues to my attention.