Custom Query (969 matches)

I tracked a very nasty bug when using the multifile options in iterative algorithms (4D-Var and friends). The multifile NetCDF option is used for open boundaries, climatology, and forcing data. For example, we can have in ocean.in:

    BRYNAME == ../om/doppio_bry_MercatorV3_2014_atmpress_NAVD88.nc |
               ../om/doppio_bry_MercatorV3_2015_atmpress_NAVD88.nc |
               ../om/doppio_bry_MercatorV3_2016_atmpress_NAVD88.nc

    CLMNAME == ../om/doppio_clm_MercatorV3_2014_atmpress_NAVD88.nc |
               ../om/doppio_clm_MercatorV3_2015_atmpress_NAVD88.nc |
               ../om/doppio_clm_MercatorV3_2016_atmpress_NAVD88.nc

    ...

    NFFILES == 8                          ! number of unique forcing files

    FRCNAME == ../om/lwrad_down_nam_3hourly_MAB_and_GoM_2014.nc |
               ../om/lwrad_down_nam_3hourly_MAB_and_GoM_2015.nc |
               ../om/lwrad_down_nam_3hourly_MAB_and_GoM_2016.nc \

               ../om/Pair_nam_3hourly_MAB_and_GoM_2014.nc |
               ../om/Pair_nam_3hourly_MAB_and_GoM_2015.nc |
               ../om/Pair_nam_3hourly_MAB_and_GoM_2016.nc \

               ../om/Qair_nam_3hourly_MAB_and_GoM_2014.nc |
               ../om/Qair_nam_3hourly_MAB_and_GoM_2015.nc |
               ../om/Qair_nam_3hourly_MAB_and_GoM_2016.nc \

               ../om/rain_nam_3hourly_MAB_and_GoM_2014.nc |
               ../om/rain_nam_3hourly_MAB_and_GoM_2015.nc |
               ../om/rain_nam_3hourly_MAB_and_GoM_2016.nc \

               ../om/swrad_daily_nam_3hourly_MAB_and_GoM_2014.nc |
               ../om/swrad_daily_nam_3hourly_MAB_and_GoM_2015.nc |
               ../om/swrad_daily_nam_3hourly_MAB_and_GoM_2016.nc \

               ../om/Tair_nam_3hourly_MAB_and_GoM_2014.nc |
               ../om/Tair_nam_3hourly_MAB_and_GoM_2015.nc |
               ../om/Tair_nam_3hourly_MAB_and_GoM_2016.nc \

               ../om/Uwind_nam_3hourly_MAB_and_GoM_2014.nc |
               ../om/Uwind_nam_3hourly_MAB_and_GoM_2015.nc |
               ../om/Uwind_nam_3hourly_MAB_and_GoM_2016.nc \

               ../om/Vwind_nam_3hourly_MAB_and_GoM_2014.nc |
               ../om/Vwind_nam_3hourly_MAB_and_GoM_2015.nc |
               ../om/Vwind_nam_3hourly_MAB_and_GoM_2016.nc

Here, the input data is split into annual NetCDF files.

The problem was the location of close_inp in initial.F, tl_initial.F, rp_initial.F, and ad_initial.F. It needs to be before the call to check_multifile and not after because it erased all the initialization in BRY(:), CLM(:), and FRC(:,:) structure that is essential for processing that appropriate multifile. For example, in tl_initia.F we need to have instead:

!
!-----------------------------------------------------------------------
!  If applicable, close all input boundary, climatology, and forcing
!  NetCDF files and set associated parameters to the closed state. This
!  step is essential in iterative algorithms that run the full TLM
!  repetitively. Then, Initialize several parameters in their file
!  structure, so the appropriate input single or multi-file is selected
!  during initialization/restart.
!-----------------------------------------------------------------------
!
!$OMP MASTER
      CALL close_inp (ng, iTLM)
      CALL check_multifile (ng, iTLM)
!$OMP END MASTER
# ifdef DISTRIBUTE
      CALL mp_bcasti (ng, iTLM, exit_flag)
# endif
!$OMP BARRIER
      IF (FoundError(exit_flag, NoError, __LINE__,                      &
     &               __FILE__)) RETURN

The call to close_inp is necessary for iterative algorithms that call the nonlinear, tangent, representer, and adjoint models repetitively during a simulation because the input data needs to be processed from start (end) records during forward (backward) time-stepping.

Also, I corrected an out bounds error in routine close_inp of file close_io.F when closing the boundary, climatology, and forcing NetCDF files. For example, we needed to have:

         FRC(i,ng)%Fcount=Fcount
         FRC(i,ng)%name=TRIM(FRC(i,ng)%files(Fcount))
         lstr=LEN_TRIM(FRC(i,ng)%name)
         FRC(i,ng)%base=FRC(i,ng)%name(1:lstr-3)
         CALL netcdf_close (ng, model, FRC(i,ng)%ncid,                 &
    &                       FRC(i,ng)%name, .FALSE.)
         IF (FoundError(exit_flag, NoError, __LINE__,                  &
    &                   __FILE__)) RETURN

         ...

         BRY(ng)%Fcount=Fcount
         BRY(ng)%name=TRIM(BRY(ng)%files(Fcount))
         lstr=LEN_TRIM(BRY(ng)%name)
         BRY(ng)%base=BRY(ng)%name(1:lstr-3)
         CALL netcdf_close (ng, model, BRY(ng)%ncid,                   &
    &                       BRY(ng)%name, .FALSE.)
         IF (FoundError(exit_flag, NoError, __LINE__,                  &
    &                   __FILE__)) RETURN

         ...

         CLM(ng)%Fcount=Fcount
         CLM(ng)%name=TRIM(CLM(ng)%files(Fcount))
         lstr=LEN_TRIM(CLM(ng)%name)
         CLM(ng)%base=CLM(ng)%name(1:lstr-3)
         CALL netcdf_close (ng, model, CLM(ng)%ncid,                   &
    &                       CLM(ng)%name, .FALSE.)
         IF (FoundError(exit_flag, NoError, __LINE__,                  &
    &                   __FILE__)) RETURN

We need to pass FRC(i,ng)%name, BRY(ng)%name, and CLM(ng)%name arguments to routine netcdf_close.

Many thanks to Andy Moore for bringing this to my attention and helping me to set-up a case that I can follow in the debugger.

In ROMS/Utility/checkvars.F, the argument ncid is not declared in the body of the function. Line 53 of this file needs to be changed to

      integer, intent(in) :: ng, model, Nvar, ncid

Oddly, the Cray ftn compiler is the only one I have that picked up this one. Something to do with compiler options, I guess.

I have been looking at the processing of input multi-files in the adjoint model (backward time-stepping) for a while. The logic breaks down at the transition between files in the list. Usually, we use annual NetCDF files for model boundary conditions, climatology, and surface forcing in our applications. The adjoint logic for backward time-stepping is very tricky. I rewrote this logic six times over the last two-weeks and spent lots of hours in the debugger looking for a robust and generic strategy.

An example of the multifile option in ocean.in is shown below:

   NFFILES == 8                          ! number of unique forcing files

   FRCNAME == ../om/lwrad_down_nam_3hourly_MAB_and_GoM_2014.nc |
              ../om/lwrad_down_nam_3hourly_MAB_and_GoM_2015.nc |
              ../om/lwrad_down_nam_3hourly_MAB_and_GoM_2016.nc \

              ../om/Pair_nam_3hourly_MAB_and_GoM_2014.nc |
              ../om/Pair_nam_3hourly_MAB_and_GoM_2015.nc |
              ../om/Pair_nam_3hourly_MAB_and_GoM_2016.nc \

              ../om/Qair_nam_3hourly_MAB_and_GoM_2014.nc |
              ../om/Qair_nam_3hourly_MAB_and_GoM_2015.nc |
              ../om/Qair_nam_3hourly_MAB_and_GoM_2016.nc \

              ../om/rain_nam_3hourly_MAB_and_GoM_2014.nc |
              ../om/rain_nam_3hourly_MAB_and_GoM_2015.nc |
              ../om/rain_nam_3hourly_MAB_and_GoM_2016.nc \

              ../om/swrad_daily_nam_3hourly_MAB_and_GoM_2014.nc |
              ../om/swrad_daily_nam_3hourly_MAB_and_GoM_2015.nc |
              ../om/swrad_daily_nam_3hourly_MAB_and_GoM_2016.nc \

              ../om/Tair_nam_3hourly_MAB_and_GoM_2014.nc |
              ../om/Tair_nam_3hourly_MAB_and_GoM_2015.nc |
              ../om/Tair_nam_3hourly_MAB_and_GoM_2016.nc \

              ../om/Uwind_nam_3hourly_MAB_and_GoM_2014.nc |
              ../om/Uwind_nam_3hourly_MAB_and_GoM_2015.nc |
              ../om/Uwind_nam_3hourly_MAB_and_GoM_2016.nc \

              ../om/Vwind_nam_3hourly_MAB_and_GoM_2014.nc |
              ../om/Vwind_nam_3hourly_MAB_and_GoM_2015.nc |
              ../om/Vwind_nam_3hourly_MAB_and_GoM_2016.nc

Each record of every multifile is processed for the data time interpolation during ROMS time stepping.

It is highly recommended that the time variable in the NetCDF multifile to be monotonic increasing without repeating dates. However, the new logic also works for repetitive dates cases as shown below.

As an example, the time variable associated with the data shows non-repetitive and repetitive dates between the multi-files in the list:

>> check_records(Pair, 'time');              % non-repetitive dates
 
File : Pair_nam_3hourly_MAB_and_GoM_2014.nc
       Record(1)     = 01-Jan-2014 00:00:00  (2922)
       Record(2)     = 01-Jan-2014 03:00:00  (2922.125)
       ...
       Record(end-1) = 31-Dec-2014 18:00:00  (3286.75)
       Record(end  ) = 31-Dec-2014 21:00:00  (3286.875)
 
File : Pair_nam_3hourly_MAB_and_GoM_2015.nc
       Record(1)     = 01-Jan-2015 00:00:00  (3287)
       Record(2)     = 01-Jan-2015 03:00:00  (3287.125)
       ...
       Record(end-1) = 31-Dec-2015 18:00:00  (3651.75)
       Record(end  ) = 31-Dec-2015 21:00:00  (3651.875)
 
File : Pair_nam_3hourly_MAB_and_GoM_2016.nc
       Record(1)     = 01-Jan-2016 00:00:00  (3652)
       Record(2)     = 01-Jan-2016 03:00:00  (3652.125)
       ...
       Record(end-1) = 31-Dec-2016 18:00:00  (4017.75)
       Record(end  ) = 31-Dec-2016 21:00:00  (4017.875)

>> check_records(Pair, 'time');              % repetitive dates
 
File : Pair_nam_3hourly_MAB_and_GoM_2014.nc
       Record(1)     = 01-Jan-2014 00:00:00  (2922)
       Record(2)     = 01-Jan-2014 03:00:00  (2922.125)
       ...
       Record(end-1) = 31-Dec-2014 21:00:00  (3286.875)
       Record(end  ) = 01-Jan-2015 00:00:00  (3287)
 
File : Pair_nam_3hourly_MAB_and_GoM_2015.nc
       Record(1)     = 01-Jan-2015 00:00:00  (3287)
       Record(2)     = 01-Jan-2015 03:00:00  (3287.125)
       ...
       Record(end-1) = 31-Dec-2015 21:00:00  (3651.875)
       Record(end  ) = 01-Jan-2016 00:00:00  (3652)
 
File : Pair_nam_3hourly_MAB_and_GoM_2016.nc
       Record(1)     = 01-Jan-2016 00:00:00  (3652)
       Record(2)     = 01-Jan-2016 03:00:00  (3652.125)
       ...
       Record(end-1) = 31-Dec-2016 21:00:00  (4017.875)
       Record(end  ) = 01-Jan-2017 00:00:00  (4018)

Notice that in the files repetitive dates, the data for 01-Jan-2015 00:00:00 appears in both the 2014 and 2015 file. We recommend you to have it only in the 2015 file.

Designing a sound logic that supports both non-repetitive and repetitive dates in the files, as shown above, was difficult. The changes were documented well in the code for future reference. During initialization (initial, tl_initial, rp_initial, or ad_initial), the routines check_multifile.F, inquiry.F, and get_cycle.F process the information for the first snapshot of the time interpolation. In the forward time stepping (NLM, TLM, and RPM), the first data snapshot corresponds to the LOWER interpolant. Contrarily, in backward time stepping (ADM), the first data snapshot corresponds to the UPPER interpolant. At initialization, there is no information yet if a particular field is from a multifile or not. So the strategy is to treat such field differently and resolve its multifile status latter when processing the second snapshot (UPPER interpolant for the forward and LOWER interpolant for the backward time stepping).

Two new logical switches were added to information variables, Linfo(5,ifield,ng) and Linfo(6,ifield,ng):

• Linfo(5,ifield,ng) is activated in get_ngfld.F, get_2dfld.F, and get_3dfld.F during forward time-stepping to indicate the processing of the LAST record of the file for the LOWER time snapshot. We need to get the UPPER time snapshot from NEXT multifile.

• Linfo(6,ifield,ng) is activated in get_ngfldr.F, get_2dfldr.F, and get_3dfldr.F during backward time-stepping to indicate the processing of the FIRST record of the file for the UPPER time snapshot. We need to get the LOWER time snapshot from PREVIOUS multifile.

This update to the code is essential because in the future we may want to split the 4D-Var forward trajectory into multi files in applications with large grids to avoid the basic state NetCDF files getting too big.

Many thanks to Andy Moore and Julia Levin for reporting the issue of processing multi-files across the date interface in their 4D-Var applications.