Custom Query (964 matches)

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.

This update contains several minor corrections:

• check_multifile.F: Added the source/sink files (SSF structure) for multifile management. Splitting the SSF input data into multiple files is now possible. Also, corrected the forcing data time range checking in the adjoint model. We need to have instead:

          IF (Lcheck) THEN
            IF (model.eq.iADM) THEN
              Tmin=Tscale*FRC(i,ng)%time_min(1)
              ...
            END IF
          END IF
  

  Many thanks to Julia Levin for bringing this problem to my attention.

• qet_cycle.F: Replaced the multifile LowerBound logic with an ELSE statement to also allow computing the LOWER time-snapshot in the time cycling case:

  !
  !  If not processing a multi-file field, find LOWER time-snapshot
  !  and its associated time record (Tindex).
  !
        ELSE
          IF (job.gt.0) THEN     ! forward:   Tval(i) =< mday =< Tval(i+1)
            IF (Linside) THEN
              tstart=Tmin
              DO i=2,ntime
                IF ((tstart.le.mday).and.(mday.le.Tval(i))) THEN
                  Tstr=tstart
                  Tindex=i-1     ! one is added when processing
                  EXIT
                END IF
                tstart=Tval(i)
              END DO
            ELSE
              Tstr=Tmax          ! LowerBound for next multifile or
              Tindex=ntime       ! time cycling
            END IF
          ELSE
            ...
          ENDIF
        ENDIF
  

  Many thanks to Dezhou Yang for ​reporting this issue.

• read_asspar.F: Corrected loading of 4D-Var parameters HdecayB and VdecayB. The order of this input values in input script s4dvar.in is (4,1:NT,Ngrids):

             CASE ('HdecayB(isTvar)')
               Npts=load_r(Nval, Rval, 4*MT*Ngrids, Rboundary)
               DO ng=1,Ngrids
                 DO itrc=1,NT(ng)
                   DO ib=1,4
                     HdecayB(isTvar(itrc),ib,ng)=Rboundary(ib,itrc,ng)
                   END DO
                 END DO
               END DO
               ...
             CASE ('VdecayB(isTvar)')
               Npts=load_r(Nval, Rval, 4*MT*Ngrids, Rboundary)
               DO ng=1,Ngrids
                 DO itrc=1,NT(ng)
                   DO ib=1,4
                     VdecayB(isTvar(itrc),ib,ng)=Rboundary(ib,itrc,ng)
                   END DO
                 END DO
               END DO
  

  The documentation in s4dvar.in indicates that the tracer boundary decorrelation scales need to be specified in the (1:4,1:NT,1:Ngrids) to facilitate having different values for each tracer. For example for Ngrids=1 and NT=2, we can have:

  ! Boundary conditions error covariance: horizontal, isotropic decorrelation
  ! scales (m). A value is expected for each boundary edge in the following
  ! order:
  !                  1: west  2: south  3: east  4: north
  
  HdecayB(isFsur) == 40.0d+3  40.0d+3   40.0d+3  40.0d+3  ! free-surface
  HdecayB(isUbar) == 40.0d+3  40.0d+3   40.0d+3  40.0d+3  ! 2D U-momentum
  HdecayB(isVbar) == 40.0d+3  40.0d+3   40.0d+3  40.0d+3  ! 2D V-momentum
  HdecayB(isUvel) == 40.0d+3  40.0d+3   40.0d+3  40.0d+3  ! 3D U-momentum
  HdecayB(isVvel) == 40.0d+3  40.0d+3   40.0d+3  40.0d+3  ! 3D V-momentum
  HdecayB(isTvar) == 4*40.0d+3  4*15.0d+3                 ! (4,1:NT) tracers
  
  ! Boundary conditions error covariance: vertical, isotropic decorrelation
  ! scales (m). A value is expected for each boundary edge in the following
  ! order:
  !                 1: west  2: south  3: east  4: north
  
  VdecayB(isUvel) == 10.0d0  10.0d0  10.0d0  10.0d0       ! 3D U-momentum
  VdecayB(isVvel) == 10.0d0  10.0d0  10.0d0  10.0d0       ! 3D V-momentum
  VdecayB(isTvar) == 4*10.d0  4*20.d0                     ! (4,1:NT) tracers
  

  It allows to use the compact notation to specify different values for temperature and salinity.

• Added new makefile configurations files Linux-ftn-cray.mk and Linux-ftn-intel.mk. Many thnaks to Mark Hadfield for providing us these two new compiling files.